Geant4  10.02.p01
G4LorentzConvertor.cc
Go to the documentation of this file.
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 //
26 // $Id: G4LorentzConvertor.cc 66241 2012-12-13 18:34:42Z gunter $
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 // 20100617 M. Kelsey -- Add more diagnostic messages with multiple levels
35 // 20100713 M. Kelsey -- All diagnostic messages should be verbose > 1
36 // 20110525 M. Kelsey -- Add some diagnostic messages in both ::rotate()
37 // 20110602 M. Kelsey -- Simplify some kinematics, dropping intermediate calcs
38 
39 #include "G4LorentzConvertor.hh"
40 #include "G4ThreeVector.hh"
41 #include "G4HadronicException.hh"
42 #include "G4InuclParticle.hh"
43 
44 
45 const G4double G4LorentzConvertor::small = 1.0e-10;
46 
48  : verboseLevel(0), v2(0.), ecm_tot(0.), valong(0.), degenerated(false) {}
49 
52  const G4LorentzVector& tmom, G4double tmass)
53  : verboseLevel(0), v2(0.), ecm_tot(0.), valong(0.), degenerated(false) {
54  setBullet(bmom, bmass);
55  setTarget(tmom, tmass);
56 }
57 
60  const G4InuclParticle* target)
61  : verboseLevel(0), v2(0.), ecm_tot(0.), valong(0.), degenerated(false) {
62  setBullet(bullet);
63  setTarget(target);
64 }
65 
66 // Extract four-vectors from input particles
68  setBullet(bullet->getMomentum());
69 }
70 
72  setTarget(target->getMomentum());
73 }
74 
75 
76 // Boost bullet and target four-vectors into desired frame
77 
79  if (verboseLevel > 2)
80  G4cout << " >>> G4LorentzConvertor::toTheCenterOfMass" << G4endl;
81 
82  velocity = (target_mom+bullet_mom).boostVector();
83  if (verboseLevel > 3) G4cout << " boost " << velocity << G4endl;
84 
85  // "SCM" is reverse target momentum in the CM frame
87  scm_momentum.boost(-velocity);
88  scm_momentum.setVect(-scm_momentum.vect());
89 
90  if (verboseLevel > 3) G4cout << " pscm " << scm_momentum.vect() << G4endl;
91 
93 }
94 
96  if (verboseLevel > 2)
97  G4cout << " >>> G4LorentzConvertor::toTheTargetRestFrame" << G4endl;
98 
99  velocity = target_mom.boostVector();
100  if (verboseLevel > 3) G4cout << " boost " << velocity << G4endl;
101 
102  // "SCM" is bullet momentum in the target's frame
104  scm_momentum.boost(-velocity);
105 
106  if (verboseLevel > 3) G4cout << " pseudo-pscm " << scm_momentum.vect() << G4endl;
107 
108  fillKinematics();
109 }
110 
111 // Compute kinematic quantities for rotate() functions
112 
115 
116  scm_direction = scm_momentum.vect().unit();
118 
119  v2 = velocity.mag2();
120 
121  G4double pvsq = v2 - valong*valong; // velocity perp to scm_momentum
122  if (verboseLevel > 3) G4cout << " pvsq " << pvsq << G4endl;
123 
124  degenerated = (pvsq < small);
125  if (degenerated && verboseLevel > 2)
126  G4cout << " degenerated case (already along Z) " << G4endl;
127 
128  if (verboseLevel > 3) {
129  G4cout << " v2 " << v2 << " valong " << valong
130  << " valong*valong " << valong*valong << G4endl;
131  }
132 }
133 
136  if (verboseLevel > 2)
137  G4cout << " >>> G4LorentzConvertor::backToTheLab" << G4endl;
138 
139  if (verboseLevel > 3)
140  G4cout << " at rest: px " << mom.x() << " py " << mom.y() << " pz "
141  << mom.z() << " e " << mom.e() << G4endl
142  << " v2 " << v2 << G4endl;
143 
144  G4LorentzVector mom1 = mom;
145  if (v2 > small) mom1.boost(velocity);
146 
147  if (verboseLevel > 3)
148  G4cout << " at lab: px " << mom1.x() << " py " << mom1.y() << " pz "
149  << mom1.z() << G4endl;
150 
151  return mom1;
152 }
153 
154 
155 // Bullet kinematics in target rest frame (LAB frame, usually)
156 
158  if (verboseLevel > 2)
159  G4cout << " >>> G4LorentzConvertor::getKinEnergyInTheTRS" << G4endl;
160 
162  bmom.boost(-target_mom.boostVector());
163  return bmom.e()-bmom.m();
164 }
165 
167  if (verboseLevel > 2)
168  G4cout << " >>> G4LorentzConvertor::getTRSMomentum" << G4endl;
169 
171  bmom.boost(-target_mom.boostVector());
172  return bmom.rho();
173 }
174 
176  if (verboseLevel > 2)
177  G4cout << " >>> G4LorentzConvertor::rotate(G4LorentzVector)" << G4endl;
178 
179  if (verboseLevel > 3) {
180  G4cout << " valong " << valong << " degenerated " << degenerated << G4endl
181  << " before rotation: px " << mom.x() << " py " << mom.y()
182  << " pz " << mom.z() << G4endl;
183  }
184 
185  G4LorentzVector mom_rot = mom;
186  if (!degenerated) {
187  if (verboseLevel > 2)
188  G4cout << " rotating to align with reference z axis " << G4endl;
189 
191  G4ThreeVector vxcm = scm_direction.cross(velocity);
192 
193  if (vscm.mag() > small && vxcm.mag() > small) { // Double check
194  if (verboseLevel > 3) {
195  G4cout << " reference z axis " << scm_direction
196  << " vscm " << vscm << " vxcm " << vxcm << G4endl;
197  }
198 
199  mom_rot.setVect(mom.x()*vscm.unit() + mom.y()*vxcm.unit() +
200  mom.z()*scm_direction);
201  } else {
202  if (verboseLevel)
203  G4cerr << ">>> G4LorentzVector::rotate zero with !degenerated" << G4endl;
204  }
205  }
206 
207  if (verboseLevel > 3) {
208  G4cout << " after rotation: px " << mom_rot.x() << " py " << mom_rot.y()
209  << " pz " << mom_rot.z() << G4endl;
210  }
211 
212  return mom_rot;
213 }
214 
216  const G4LorentzVector& mom) const {
217  if (verboseLevel > 2)
218  G4cout << " >>> G4LorentzConvertor::rotate(G4LorentzVector,G4LorentzVector)"
219  << G4endl;
220 
221  if (verboseLevel > 3) {
222  G4cout << " before rotation: px " << mom.x() << " py " << mom.y()
223  << " pz " << mom.z() << G4endl;
224  }
225 
226  G4ThreeVector mom1_dir = mom1.vect().unit();
227  G4double pv = velocity.dot(mom1_dir);
228 
229  G4double vperp = v2 - pv*pv; // velocity perpendicular to mom1
230  if (verboseLevel > 3) {
231  G4cout << " vperp " << vperp << " small? " << (vperp <= small) << G4endl;
232  }
233 
234  G4LorentzVector mom_rot = mom;
235 
236  if (vperp > small) {
237  if (verboseLevel > 2)
238  G4cout << " rotating to align with first z axis " << G4endl;
239 
240  G4ThreeVector vmom1 = velocity - pv*mom1_dir;
241  G4ThreeVector vxm1 = mom1_dir.cross(velocity);
242 
243  if (vmom1.mag() > small && vxm1.mag() > small) { // Double check
244  if (verboseLevel > 3) {
245  G4cout << " first z axis " << mom1_dir << G4endl
246  << " vmom1 " << vmom1 << " vxm1 " << vxm1 << G4endl;
247  }
248 
249  mom_rot.setVect(mom.x()*vmom1.unit() + mom.y()*vxm1.unit() +
250  mom.z()*mom1_dir );
251  } else {
252  if (verboseLevel)
253  G4cerr << ">>> G4LorentzVector::rotate zero with !degenerated" << G4endl;
254  }
255  }
256 
257  if (verboseLevel > 3) {
258  G4cout << " after rotation: px " << mom_rot.x() << " py " << mom_rot.y()
259  << " pz " << mom_rot.z() << G4endl;
260  }
261 
262  return mom_rot;
263 }
264 
266  if (verboseLevel > 2)
267  G4cout << " >>> G4LorentzConvertor::reflectionNeeded (query)" << G4endl;
268 
269  if (verboseLevel > 3) {
270  G4cout << " v2 = " << v2 << " SCM z = " << scm_momentum.z()
271  << " degenerated? " << degenerated << G4endl;
272  }
273 
274  if (v2 < small && !degenerated)
275  throw G4HadronicException(__FILE__, __LINE__, "G4LorentzConvertor::reflectionNeeded - return value undefined");
276 
277  if (verboseLevel > 2) {
278  G4cout << " reflection across XY is"
279  << ((v2>=small && (!degenerated || scm_momentum.z()<0.0))?"":" NOT")
280  << " needed" << G4endl;
281  }
282 
283  return (v2>=small && (!degenerated || scm_momentum.z()<0.0));
284 }
285 
286 
287 // Reporting functions for diagnostics
288 
290  G4cout << " G4LC bullet: px " << bullet_mom.px() << " py " << bullet_mom.py()
291  << " pz " << bullet_mom.pz() << " e " << bullet_mom.e()
292  << " mass " << bullet_mom.m() << G4endl;
293  }
294 
296  G4cout << " G4LC target: px " << target_mom.px() << " py " << target_mom.py()
297  << " pz " << target_mom.pz() << " e " << target_mom.e()
298  << " mass " << target_mom.m() << G4endl;
299 }
300 
G4LorentzVector bullet_mom
G4double getTRSMomentum() const
CLHEP::Hep3Vector G4ThreeVector
G4LorentzVector getMomentum() const
G4ThreeVector scm_direction
G4LorentzVector rotate(const G4LorentzVector &mom) const
void setBullet(const G4InuclParticle *bullet)
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4LorentzVector scm_momentum
G4double getKinEnergyInTheTRS() const
G4bool reflectionNeeded() const
#define G4endl
Definition: G4ios.hh:61
static const double m
Definition: G4SIunits.hh:128
G4LorentzVector target_mom
double G4double
Definition: G4Types.hh:76
static const G4double small
G4GLOB_DLL std::ostream G4cerr
void setTarget(const G4InuclParticle *target)
CLHEP::HepLorentzVector G4LorentzVector