source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4LorentzConvertor.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 8.6 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25// $Id: G4LorentzConvertor.cc,v 1.23 2010/05/21 17:56:34 mkelsey Exp $
26// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
27//
28// 20100108  Michael Kelsey -- Use G4LorentzVector internally
29// 20100112  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
30// 20100308  M. Kelsey -- Bug fix in toTheTargetRestFrame: scm_momentum should
31//              be data member, not local.
32// 20100409  M. Kelsey -- Protect std::sqrt(ga) against round-off negatives
33// 20100519  M. Kelsey -- Add interfaces to pass G4InuclParticles directly
34
35#include "G4LorentzConvertor.hh"
36#include "G4ThreeVector.hh"
37#include "G4HadronicException.hh"
38#include "G4InuclParticle.hh"
39
40
41const G4double G4LorentzConvertor::small = 1.0e-10;
42
43G4LorentzConvertor::G4LorentzConvertor() 
44  : verboseLevel(0), degenerated(false) {
45
46  if (verboseLevel > 3) {
47    G4cout << " >>> G4LorentzConvertor::G4LorentzConvertor" << G4endl;
48  }
49}
50
51// Extract four-vectors from input particles
52void G4LorentzConvertor::setBullet(const G4InuclParticle* bullet) {
53  setBullet(bullet->getMomentum());
54}
55
56void G4LorentzConvertor::setTarget(const G4InuclParticle* target) {
57  setTarget(target->getMomentum());
58}
59
60// Boost bullet and target four-vectors into destired frame
61
62void G4LorentzConvertor::toTheCenterOfMass() {
63  if (verboseLevel > 3) {
64    G4cout << " >>> G4LorentzConvertor::toTheCenterOfMass" << G4endl;
65  }
66
67  G4LorentzVector cm4v = target_mom + bullet_mom;
68  velocity = cm4v.boostVector();
69
70  // "SCM" is reverse target momentum in the CM frame
71  scm_momentum = target_mom;
72  scm_momentum.boost(-velocity);
73  scm_momentum.setVect(-scm_momentum.vect());
74
75  if (verboseLevel > 3)
76    G4cout << " i 1 pscm(i) " << scm_momentum.x() << G4endl
77           << " i 2 pscm(i) " << scm_momentum.y() << G4endl
78           << " i 3 pscm(i) " << scm_momentum.z() << G4endl;
79
80  // Compute kinematic quantities for rotate() functions
81  v2 = velocity.mag2();
82  gamma = cm4v.e()/cm4v.m();
83
84  ecm_tot = cm4v.m();
85
86  G4double pscm = scm_momentum.rho();
87  G4double pa   = scm_momentum.vect().mag2();
88  G4double pb   = scm_momentum.vect().dot(velocity);
89
90  G4double gasq = v2-pb*pb/pa;
91  ga = (gasq > small*small) ? std::sqrt(gasq) : 0.;
92
93  gb = pb / pscm;
94  gbpp = gb / pscm;
95  gapp = ga * pscm;
96
97  degenerated = (ga < small);
98  if (degenerated && verboseLevel > 3) 
99    G4cout << " degenerated case " << G4endl; 
100
101  if (verboseLevel > 3) {
102    G4double pv = target_mom.vect().dot(velocity);
103    G4cout << " ga " << ga << " v2 " << v2 << " pb " << pb
104           << " pb * pb / pa " << pb * pb / pa << " pv " << pv << G4endl;
105  }
106}
107
108void G4LorentzConvertor::toTheTargetRestFrame() {
109  if (verboseLevel > 3) {
110    G4cout << " >>> G4LorentzConvertor::toTheTargetRestFrame" << G4endl;
111  }
112
113  velocity = target_mom.boostVector();
114
115  // "SCM" is bullet momentum in the target's frame
116  scm_momentum = bullet_mom;
117  scm_momentum.boost(-velocity);
118
119  if (verboseLevel > 3)
120    G4cout << " rf: i 1 pscm(i) " << scm_momentum.x() << G4endl
121           << " rf: i 2 pscm(i) " << scm_momentum.y() << G4endl
122           << " rf: i 3 pscm(i) " << scm_momentum.z() << G4endl;
123
124  // Compute kinematic quantities for rotate() functions
125  v2 = velocity.mag2();
126  gamma = target_mom.e() / target_mom.m();
127
128  G4double pscm = scm_momentum.rho();
129  G4double pa   = scm_momentum.vect().mag2();
130  G4double pb   = velocity.dot(scm_momentum.vect());
131
132  G4double gasq = v2-pb*pb/pa;
133  ga = (gasq > small*small) ? std::sqrt(gasq) : 0.;
134
135  gb = pb/pscm;
136  gbpp = gb/pscm;
137  gapp = ga*pscm;
138
139  degenerated = (ga < small);
140  if (degenerated && verboseLevel > 3) 
141    G4cout << " degenerated case " << G4endl; 
142}
143
144G4LorentzVector
145G4LorentzConvertor::backToTheLab(const G4LorentzVector& mom) const {
146  if (verboseLevel > 3) {
147    G4cout << " >>> G4LorentzConvertor::backToTheLab" << G4endl
148           << " at rest: px " << mom.x() << " py " << mom.y() << " pz "
149           << mom.z() << " e " << mom.e() << G4endl
150           << " v2 " << v2 << G4endl;
151  }
152
153  G4LorentzVector mom1 = mom;
154  if (v2 > small) mom1.boost(velocity);
155
156  if (verboseLevel > 3)
157    G4cout << " at lab: px " << mom1.x() << " py " << mom1.y() << " pz "
158           << mom1.z() << G4endl;
159
160  return mom1;
161}
162
163
164// Bullet kinematics in target rest frame (LAB frame, usually)
165
166G4double G4LorentzConvertor::getKinEnergyInTheTRS() const {
167  G4double pv = bullet_mom.vect().dot(target_mom.vect());
168 
169   G4double ekin_trf =
170    (target_mom.e() * bullet_mom.e() - pv) / target_mom.m()
171    - bullet_mom.m();
172 
173  return ekin_trf; 
174}
175
176G4double G4LorentzConvertor::getTRSMomentum() const {
177  G4LorentzVector bmom = bullet_mom;
178  bmom.boost(-target_mom.boostVector());
179  return bmom.rho();
180}
181
182G4LorentzVector G4LorentzConvertor::rotate(const G4LorentzVector& mom) const {
183  if (verboseLevel > 3) {
184    G4cout << " >>> G4LorentzConvertor::rotate(G4LorentzVector)" << G4endl
185           << " ga " << ga << " gbpp " << gbpp << " gapp " << gapp << G4endl
186           << " degenerated " << degenerated << G4endl
187           << " before rotation: px " << mom.x() << " py " << mom.y()
188           << " pz " << mom.z() << G4endl;
189  }
190
191  G4LorentzVector mom_rot = mom;
192  if (!degenerated) {
193    G4ThreeVector vscm = velocity - gbpp*scm_momentum.vect();
194    G4ThreeVector vxcm = scm_momentum.vect().cross(velocity);
195
196    mom_rot.setVect(mom.x()*vscm/ga + mom.y()*vxcm/gapp +
197                    mom.z()*scm_momentum.vect().unit() );
198  };
199
200  if (verboseLevel > 3) {
201    G4cout << " after rotation: px " << mom_rot.x() << " py " << mom_rot.y()
202           << " pz " << mom_rot.z() << G4endl;
203  }
204
205  return mom_rot;
206}
207
208G4LorentzVector G4LorentzConvertor::rotate(const G4LorentzVector& mom1, 
209                                           const G4LorentzVector& mom) const {
210  if (verboseLevel > 3) {
211    G4cout << " >>> G4LorentzConvertor::rotate(G4LorentzVector,G4LorentzVector)"
212           << G4endl
213           << " before rotation: px " << mom.x() << " py " << mom.y()
214           << " pz " << mom.z() << G4endl;
215  }
216
217  G4double pp = mom1.vect().mag2();
218  G4double pv = mom1.vect().dot(velocity);
219
220  G4double ga1 = v2 - pv * pv / pp;
221  if (verboseLevel > 3) {
222    G4cout << " ga1 " << ga1 << " small? " << (ga1 <= small) << G4endl;
223  }
224
225  G4LorentzVector mom_rot = mom;
226
227  if (ga1 > small) {
228    ga1 = std::sqrt(ga1);
229
230    G4double gb1 = pv / pp;
231
232    pp = std::sqrt(pp);
233
234    G4double ga1pp = ga1 * pp;
235
236    if (verboseLevel > 3) {
237      G4cout << " gb1 " << gb1 << " ga1pp " << ga1pp << G4endl;
238    }
239
240    G4ThreeVector vmom1 = velocity - gb1*mom1;
241    G4ThreeVector vxm1  = mom1.vect().cross(velocity);
242
243    mom_rot.setVect(mom.x()*vmom1/ga1 + mom.y()*vxm1/ga1pp +
244                    mom.z()*mom1.vect().unit() );
245  };
246
247  if (verboseLevel > 3) {
248    G4cout << " after rotation: px " << mom_rot.x() << " py " << mom_rot.y()
249           << " pz " << mom_rot.z() << G4endl;
250  }
251
252  return mom_rot;
253}
254
255G4bool G4LorentzConvertor::reflectionNeeded() const {
256  if (verboseLevel > 3) {
257    G4cout << " >>> G4LorentzConvertor::reflectionNeeded: v2 = " << v2
258           << " SCM z = " << scm_momentum.z() << " degenerated? "
259           << degenerated << G4endl;
260  }
261
262  if (v2 < small && !degenerated) 
263    throw G4HadronicException(__FILE__, __LINE__, "G4LorentzConvertor::reflectionNeeded - return value undefined");
264
265  return (v2>=small && (!degenerated || scm_momentum.z() < 0.0));
266}
267
268
269
270
271
272
273
Note: See TracBrowser for help on using the repository browser.