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

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

update ti head
File size: 10.8 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.29 2010/09/15 20:16:16 mkelsey Exp $
26// Geant4 tag: $Name: hadr-casc-V09-03-85 $
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
37#include "G4LorentzConvertor.hh"
38#include "G4ThreeVector.hh"
39#include "G4HadronicException.hh"
40#include "G4InuclParticle.hh"
41
42
43const G4double G4LorentzConvertor::small = 1.0e-10;
44
45G4LorentzConvertor::G4LorentzConvertor() 
46  : verboseLevel(0), velocity(0.), gamma(0.), v2(0.), ecm_tot(0.),
47    ga(0.), gb(0.), gbpp(0.), gapp(0.), degenerated(false) {
48  if (verboseLevel)
49    G4cout << " >>> G4LorentzConvertor::G4LorentzConvertor" << G4endl;
50}
51
52G4LorentzConvertor::
53G4LorentzConvertor(const G4LorentzVector& bmom, G4double bmass, 
54                   const G4LorentzVector& tmom, G4double tmass) 
55  : verboseLevel(0), velocity(0.), gamma(0.), v2(0.), ecm_tot(0.),
56    ga(0.), gb(0.), gbpp(0.), gapp(0.), degenerated(false) {
57  setBullet(bmom, bmass);
58  setTarget(tmom, tmass);
59}
60
61G4LorentzConvertor::
62G4LorentzConvertor(const G4InuclParticle* bullet, 
63                   const G4InuclParticle* target) 
64  : verboseLevel(0), velocity(0.), gamma(0.), v2(0.), ecm_tot(0.),
65    ga(0.), gb(0.), gbpp(0.), gapp(0.), degenerated(false) {
66  setBullet(bullet);
67  setTarget(target);
68}
69
70// Extract four-vectors from input particles
71void G4LorentzConvertor::setBullet(const G4InuclParticle* bullet) {
72  setBullet(bullet->getMomentum());
73}
74
75void G4LorentzConvertor::setTarget(const G4InuclParticle* target) {
76  setTarget(target->getMomentum());
77}
78
79
80// Boost bullet and target four-vectors into destired frame
81
82void G4LorentzConvertor::toTheCenterOfMass() {
83  if (verboseLevel > 2)
84    G4cout << " >>> G4LorentzConvertor::toTheCenterOfMass" << G4endl;
85
86  G4LorentzVector cm4v = target_mom + bullet_mom;
87  velocity = cm4v.boostVector();
88  if (verboseLevel > 3)
89    G4cout << " boost " << velocity.x() << " " << velocity.y() << " "
90           << velocity.z() << G4endl;
91
92  // "SCM" is reverse target momentum in the CM frame
93  scm_momentum = target_mom;
94  scm_momentum.boost(-velocity);
95  scm_momentum.setVect(-scm_momentum.vect());
96
97  if (verboseLevel > 3)
98    G4cout << " pscm " << scm_momentum.x() << " " << scm_momentum.y()
99           << " " << scm_momentum.z() << G4endl;
100
101  // Compute kinematic quantities for rotate() functions
102  v2 = velocity.mag2();
103  ecm_tot = cm4v.m();
104  gamma = cm4v.e()/ecm_tot;
105
106  G4double pscm = scm_momentum.rho();
107  G4double pa   = pscm*pscm;
108  G4double pb   = scm_momentum.vect().dot(velocity);
109
110  G4double gasq = v2-pb*pb/pa;
111  ga = (gasq > small*small) ? std::sqrt(gasq) : 0.;
112
113  gb = pb / pscm;
114  gbpp = gb / pscm;
115  gapp = ga * pscm;
116
117  degenerated = (ga < small);
118  if (degenerated && verboseLevel > 2) 
119    G4cout << " degenerated case (already in CM frame) " << G4endl; 
120
121  if (verboseLevel > 3) {
122    G4double pv = target_mom.vect().dot(velocity);
123    G4cout << " ga " << ga << " v2 " << v2 << " pb " << pb
124           << " pb * pb / pa " << pb * pb / pa << " pv " << pv << G4endl;
125  }
126}
127
128void G4LorentzConvertor::toTheTargetRestFrame() {
129  if (verboseLevel > 2)
130    G4cout << " >>> G4LorentzConvertor::toTheTargetRestFrame" << G4endl;
131
132  velocity = target_mom.boostVector();
133  if (verboseLevel > 3)
134    G4cout << " boost " << velocity.x() << " " << velocity.y() << " "
135           << velocity.z() << G4endl;
136
137  // "SCM" is bullet momentum in the target's frame
138  scm_momentum = bullet_mom;
139  scm_momentum.boost(-velocity);
140
141  if (verboseLevel > 3)
142    G4cout << " pseudo-pscm " << scm_momentum.x() << " " << scm_momentum.y()
143           << " " << scm_momentum.z() << G4endl;
144
145  // Compute kinematic quantities for rotate() functions
146  v2 = velocity.mag2();
147  gamma = target_mom.e() / target_mom.m();
148
149  G4double pscm = scm_momentum.rho();
150  G4double pa   = pscm*pscm;
151  G4double pb   = velocity.dot(scm_momentum.vect());
152
153  G4double gasq = v2-pb*pb/pa;
154  ga = (gasq > small*small) ? std::sqrt(gasq) : 0.;
155
156  gb = pb/pscm;
157  gbpp = gb/pscm;
158  gapp = ga*pscm;
159
160  degenerated = (ga < small);
161  if (degenerated && verboseLevel > 2) 
162    G4cout << " degenerated case (already in target frame) " << G4endl; 
163
164  if (verboseLevel > 3) {
165    G4cout << " ga " << ga << " v2 " << v2 << " pb " << pb
166           << " pb * pb / pa " << pb * pb / pa << G4endl;
167  }
168}
169
170G4LorentzVector
171G4LorentzConvertor::backToTheLab(const G4LorentzVector& mom) const {
172  if (verboseLevel > 2)
173    G4cout << " >>> G4LorentzConvertor::backToTheLab" << G4endl;
174
175  if (verboseLevel > 3)
176    G4cout << " at rest: px " << mom.x() << " py " << mom.y() << " pz "
177           << mom.z() << " e " << mom.e() << G4endl
178           << " v2 " << v2 << G4endl;
179
180  G4LorentzVector mom1 = mom;
181  if (v2 > small) mom1.boost(velocity);
182
183  if (verboseLevel > 3)
184    G4cout << " at lab: px " << mom1.x() << " py " << mom1.y() << " pz "
185           << mom1.z() << G4endl;
186
187  return mom1;
188}
189
190
191// Bullet kinematics in target rest frame (LAB frame, usually)
192
193G4double G4LorentzConvertor::getKinEnergyInTheTRS() const {
194  if (verboseLevel > 2)
195    G4cout << " >>> G4LorentzConvertor::getKinEnergyInTheTRS" << G4endl;
196
197  G4double pv = bullet_mom.vect().dot(target_mom.vect());
198 
199   G4double ekin_trf =
200    (target_mom.e() * bullet_mom.e() - pv) / target_mom.m()
201    - bullet_mom.m();
202 
203  return ekin_trf; 
204}
205
206G4double G4LorentzConvertor::getTRSMomentum() const {
207  if (verboseLevel > 2)
208    G4cout << " >>> G4LorentzConvertor::getTRSMomentum" << G4endl;
209
210  G4LorentzVector bmom = bullet_mom;
211  bmom.boost(-target_mom.boostVector());
212  return bmom.rho();
213}
214
215G4LorentzVector G4LorentzConvertor::rotate(const G4LorentzVector& mom) const {
216  if (verboseLevel > 2)
217    G4cout << " >>> G4LorentzConvertor::rotate(G4LorentzVector)" << G4endl;
218
219  if (verboseLevel > 3) {
220    G4cout << " ga " << ga << " gbpp " << gbpp << " gapp " << gapp << G4endl
221           << " degenerated " << degenerated << G4endl
222           << " before rotation: px " << mom.x() << " py " << mom.y()
223           << " pz " << mom.z() << G4endl;
224  }
225
226  G4LorentzVector mom_rot = mom;
227  if (!degenerated) {
228    if (verboseLevel > 2)
229      G4cout << " rotating to align with reference z axis " << G4endl;
230
231    G4ThreeVector vscm = velocity - gbpp*scm_momentum.vect();
232    G4ThreeVector vxcm = scm_momentum.vect().cross(velocity);
233
234    mom_rot.setVect(mom.x()*vscm/ga + mom.y()*vxcm/gapp +
235                    mom.z()*scm_momentum.vect().unit() );
236  };
237
238  if (verboseLevel > 3) {
239    G4cout << " after rotation: px " << mom_rot.x() << " py " << mom_rot.y()
240           << " pz " << mom_rot.z() << G4endl;
241  }
242
243  return mom_rot;
244}
245
246G4LorentzVector G4LorentzConvertor::rotate(const G4LorentzVector& mom1, 
247                                           const G4LorentzVector& mom) const {
248  if (verboseLevel > 2)
249    G4cout << " >>> G4LorentzConvertor::rotate(G4LorentzVector,G4LorentzVector)"
250           << G4endl;
251
252  if (verboseLevel > 3) {
253    G4cout << " before rotation: px " << mom.x() << " py " << mom.y()
254           << " pz " << mom.z() << G4endl;
255  }
256
257  G4double pp = mom1.vect().mag2();
258  G4double pv = mom1.vect().dot(velocity);
259
260  G4double ga1 = v2 - pv * pv / pp;
261  if (verboseLevel > 3) {
262    G4cout << " ga1 " << ga1 << " small? " << (ga1 <= small) << G4endl;
263  }
264
265  G4LorentzVector mom_rot = mom;
266
267  if (ga1 > small) {
268    ga1 = std::sqrt(ga1);
269    G4double gb1 = pv / pp;
270
271    pp = std::sqrt(pp);
272    G4double ga1pp = ga1 * pp;
273
274    if (verboseLevel > 3) {
275      G4cout << " gb1 " << gb1 << " ga1pp " << ga1pp << G4endl;
276    }
277
278    if (verboseLevel > 2)
279      G4cout << " rotating to align with first z axis " << G4endl;
280
281    G4ThreeVector vmom1 = velocity - gb1*mom1;
282    G4ThreeVector vxm1  = mom1.vect().cross(velocity);
283
284    mom_rot.setVect(mom.x()*vmom1/ga1 + mom.y()*vxm1/ga1pp +
285                    mom.z()*mom1.vect().unit() );
286  };
287
288  if (verboseLevel > 3) {
289    G4cout << " after rotation: px " << mom_rot.x() << " py " << mom_rot.y()
290           << " pz " << mom_rot.z() << G4endl;
291  }
292
293  return mom_rot;
294}
295
296G4bool G4LorentzConvertor::reflectionNeeded() const {
297  if (verboseLevel > 2)
298    G4cout << " >>> G4LorentzConvertor::reflectionNeeded (query)" << G4endl;
299
300  if (verboseLevel > 3) {
301    G4cout << " v2 = " << v2 << " SCM z = " << scm_momentum.z()
302           << " degenerated? " << degenerated << G4endl;
303  }
304
305  if (v2 < small && !degenerated) 
306    throw G4HadronicException(__FILE__, __LINE__, "G4LorentzConvertor::reflectionNeeded - return value undefined");
307
308  if (verboseLevel > 2) {
309    G4cout << " reflection across XY is"
310           << ((v2>=small && (!degenerated || scm_momentum.z()<0.0))?"":" NOT")
311           << " needed" << G4endl;
312  }
313
314  return (v2>=small && (!degenerated || scm_momentum.z()<0.0));
315}
316
317
318// Reporting functions for diagnostics
319
320void G4LorentzConvertor::printBullet() const {
321  G4cout << " G4LC bullet: px " << bullet_mom.px() << " py " << bullet_mom.py()
322         << " pz " << bullet_mom.pz() << " e " << bullet_mom.e()
323         << " mass " << bullet_mom.m() << G4endl;
324  }
325
326void G4LorentzConvertor::printTarget() const {
327  G4cout << " G4LC target: px " << target_mom.px() << " py " << target_mom.py()
328         << " pz " << target_mom.pz() << " e " << target_mom.e()
329         << " mass " << target_mom.m() << G4endl;
330}
331
Note: See TracBrowser for help on using the repository browser.