ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Scatterer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Scatterer.cc
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 //
27 
28 #include <vector>
29 
30 #include "globals.hh"
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4ios.hh"
34 #include "G4Scatterer.hh"
35 #include "G4KineticTrack.hh"
36 #include "G4ThreeVector.hh"
37 #include "G4LorentzRotation.hh"
38 #include "G4LorentzVector.hh"
39 
40 #include "G4CollisionNN.hh"
41 #include "G4CollisionPN.hh"
42 #include "G4CollisionMesonBaryon.hh"
43 
44 #include "G4CollisionInitialState.hh"
45 #include "G4HadTmpUtil.hh"
46 #include "G4Pair.hh"
47 #include "G4AutoLock.hh"
48 
49 //Mutex for control of shared resource
50 namespace {
51  G4Mutex collisions_mutex = G4MUTEX_INITIALIZER;
52  G4bool setupDone = false;
53 }
54 
55 // Declare the categories of collisions the Scatterer can handle
56 typedef GROUP2(G4CollisionNN, G4CollisionMesonBaryon) theChannels;
57 
58 G4CollisionVector G4Scatterer::collisions;
59 
60 //----------------------------------------------------------------------------
61 
62 G4Scatterer::G4Scatterer()
63 {
64  G4AutoLock l(&collisions_mutex);
65  if ( ! setupDone )
66  {
67  Register aR;
68  G4ForEach<theChannels>::Apply(&aR, &collisions);
69  setupDone = true;
70  }
71 }
72 
73 //----------------------------------------------------------------------------
74 
75 G4Scatterer::~G4Scatterer()
76 {
77  G4AutoLock l(&collisions_mutex);
78  std::for_each(collisions.begin(), collisions.end(), G4Delete());
79  collisions.clear();
80 }
81 
82 //----------------------------------------------------------------------------
83 
84 G4double G4Scatterer::GetTimeToInteraction(const G4KineticTrack& trk1,
85  const G4KineticTrack& trk2) const
86 {
87  G4double time = DBL_MAX;
88  G4double distance_fast;
89  G4LorentzVector mom1 = trk1.GetTrackingMomentum();
90 // G4cout << "zcomp=" << std::abs(mom1.vect().unit().z() -1 ) << G4endl;
91  G4double collisionTime;
92 
93  if ( std::abs(mom1.vect().unit().z() -1 ) < 1e-6 )
94  {
95  G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
96  G4double deltaz=position.z();
97  G4double velocity = mom1.z()/mom1.e() * c_light;
98 
99  collisionTime=deltaz/velocity;
100  distance_fast=position.x()*position.x() + position.y()*position.y();
101  } else {
102 
103  // The nucleons of the nucleus are FROZEN, ie. do not move..
104 
105  G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
106 
107  G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light; // mom1.boostVector() will exit on slightly negative mass
108  collisionTime = (position * velocity) / velocity.mag2(); // can't divide by /c_light;
109  position -= velocity * collisionTime;
110  distance_fast=position.mag2();
111 
112 // if ( collisionTime>0 ) G4cout << " dis1/2 square" << dis1 <<" "<< dis2 << G4endl;
113 // collisionTime = GetTimeToClosestApproach(trk1,trk2);
114  }
115  if (collisionTime > 0)
116  {
117  static const G4double maxCrossSection = 500*millibarn;
118  if(0.7*pi*distance_fast>maxCrossSection) return time;
119 
120 
121  G4LorentzVector mom2(0,0,0,trk2.Get4Momentum().mag());
122 
123 // G4ThreeVector momLab = mom1.vect();// frozen Nucleus - mom2.vect();
124 // G4ThreeVector posLab = trk1.GetPosition() - trk2.GetPosition();
125 // G4double disLab=posLab * posLab - (posLab*momLab) * (posLab*momLab) /(momLab.mag2());
126 
127  G4LorentzRotation toCMSFrame((-1)*(mom1 + mom2).boostVector());
128  mom1 = toCMSFrame * mom1;
129  mom2 = toCMSFrame * mom2;
130 
131  G4LorentzVector coordinate1(trk1.GetPosition(), 100.);
132  G4LorentzVector coordinate2(trk2.GetPosition(), 100.);
133  G4ThreeVector pos = ((toCMSFrame * coordinate1).vect() -
134  (toCMSFrame * coordinate2).vect());
135 
136  G4ThreeVector mom = mom1.vect() - mom2.vect();
137 
138  // Calculate the impact parameter
139 
140  G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom.mag2());
141 
142 // G4cout << " disDiff " << distance-disLab << " " << disLab
143 // << " " << std::abs(distance-disLab)/distance << G4endl
144 // << " mom/Lab " << mom << " " << momLab << G4endl
145 // << " pos/Lab " << pos << " " << posLab
146 // << G4endl;
147 
148  if(pi*distance>maxCrossSection) return time;
149 
150  // charged particles special
151  static const G4double maxChargedCrossSection = 200*millibarn;
152  if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 &&
153  std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 &&
154  pi*distance>maxChargedCrossSection) return time;
155 
156  G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
157  // neutrons special pn is largest cross-section, but above 1.91 GeV is less than 200 mb
158  if(( trk1.GetDefinition() == G4Neutron::Neutron() ||
159  trk2.GetDefinition() == G4Neutron::Neutron() ) &&
160  sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time;
161 
162 /*
163  * if(distance <= sqr(1.14*fermi))
164  * {
165  * time = collisionTime;
166  *
167  * *
168  * * G4cout << "Scatter distance/time: " << std::sqrt(distance)/fermi <<
169  * * " / "<< time/ns << G4endl;
170  * * G4ThreeVector pos1=trk1.GetPosition();
171  * * G4ThreeVector pos2=trk2.GetPosition();
172  * * G4LorentzVector xmom1 = trk1.Get4Momentum();
173  * * G4LorentzVector xmom2 = trk2.Get4Momentum();
174  * * G4cout << "position1: " << pos1.x() << " " << pos1.y() << " "
175  * * << pos1.z();
176  * * pos1+=(collisionTime*c_light/xmom1.e())*xmom1.vect();
177  * * G4cout << " straight line trprt: "
178  * * << pos1.x() << " " << pos1.y() << " "
179  * * << pos1.z() << G4endl;
180  * * G4cout << "position2: " << pos2.x() << " " << pos2.y() << " "
181  * * << pos2.z() << G4endl;
182  * * G4cout << "straight line distance 2 fixed:" << (pos1-pos2).mag()/fermi << G4endl;
183  * * pos2+= (collisionTime*c_light/xmom2.e())*xmom2.vect();
184  * * G4cout<< " straight line trprt: "
185  * * << pos2.x() << " " << pos2.y() << " "
186  * * << pos2.z() << G4endl;
187  * * G4cout << "straight line distance :" << (pos1-pos2).mag()/fermi << G4endl;
188  * *
189  * }
190  *
191  * if(1)
192  * return time;
193  */
194 
195  if ((trk1.GetActualMass()+trk2.GetActualMass()) > sqrtS) return time;
196 
197 
198 
199  const G4VCollision* collision = FindCollision(trk1,trk2);
200  G4double totalCrossSection;
201  // The cross section is interpreted geometrically as an area
202  // Two particles are assumed to collide if their distance is < (totalCrossSection/pi)
203 
204  if (collision != 0)
205  {
206  totalCrossSection = collision->CrossSection(trk1,trk2);
207  if ( totalCrossSection > 0 )
208  {
209 /* G4cout << " totalCrossection = "<< totalCrossSection << ", trk1/2, s, e-m: "
210  * << trk1.GetDefinition()->GetParticleName()
211  * << " / "
212  * << trk2.GetDefinition()->GetParticleName()
213  * << ", "
214  * << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
215  * << ", "
216  * << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
217  * trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
218  * << G4endl;
219  */
220  if (distance <= totalCrossSection / pi)
221  {
222  time = collisionTime;
223  }
224  } else
225  {
226 
227  // For debugging...
228  // G4cout << " totalCrossection = 0, trk1/2, s, e-m: "
229  // << trk1.GetDefinition()->GetParticleName()
230  // << " / "
231  // << trk2.GetDefinition()->GetParticleName()
232  // << ", "
233  // << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
234  // << ", "
235  // << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
236  // trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
237  // << G4endl;
238 
239  }
240 /*
241  * if(distance <= sqr(5.*fermi))
242  * {
243  * G4cout << " distance,xsect, std::sqrt(xsect/pi) : " << std::sqrt(distance)/fermi
244  * << " " << totalCrossSection/sqr(fermi)
245  * << " " << std::sqrt(totalCrossSection / pi)/fermi << G4endl;
246  * }
247  */
248 
249  }
250  else
251  {
252  time = DBL_MAX;
253 // /*
254  // For debugging
255 //hpw G4cout << "G4Scatterer - collision not found: "
256 //hpw << trk1.GetDefinition()->GetParticleName()
257 //hpw << " - "
258 //hpw << trk2.GetDefinition()->GetParticleName()
259 //hpw << G4endl;
260  // End of debugging
261 // */
262  }
263  }
264 
265  else
266  {
267  /*
268  // For debugging
269  G4cout << "G4Scatterer - negative collisionTime"
270  << ": collisionTime = " << collisionTime
271  << ", position = " << position
272  << ", velocity = " << velocity
273  << G4endl;
274  // End of debugging
275  */
276  }
277 
278  return time;
279 }
280 
281 //----------------------------------------------------------------------------
282 
283 G4KineticTrackVector* G4Scatterer::Scatter(const G4KineticTrack& trk1,
284  const G4KineticTrack& trk2) const
285 {
286 // G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
287  G4LorentzVector pInitial=trk1.Get4Momentum() + trk2.Get4Momentum();
288  G4double energyBalance = pInitial.t();
289  G4double pxBalance = pInitial.vect().x();
290  G4double pyBalance = pInitial.vect().y();
291  G4double pzBalance = pInitial.vect().z();
292  G4int chargeBalance = G4lrint(trk1.GetDefinition()->GetPDGCharge()
293  + trk2.GetDefinition()->GetPDGCharge());
294  G4int baryonBalance = trk1.GetDefinition()->GetBaryonNumber()
295  + trk2.GetDefinition()->GetBaryonNumber();
296 
297  const G4VCollision* collision = FindCollision(trk1,trk2);
298  if (collision != 0)
299  {
300  G4double aCrossSection = collision->CrossSection(trk1,trk2);
301  if (aCrossSection > 0.0)
302  {
303 
304 
305  #ifdef debug_G4Scatterer
306  G4cout << "be4 FinalState 1(p,e,m): "
307  << trk1.Get4Momentum() << " "
308  << trk1.Get4Momentum().mag()
309  << ", 2: "
310  << trk2.Get4Momentum()<< " "
311  << trk2.Get4Momentum().mag() << " "
312  << G4endl;
313  #endif
314 
315 
316  G4KineticTrackVector* products = collision->FinalState(trk1,trk2);
317  if(!products || products->size() == 0) return products;
318 
319  #ifdef debug_G4Scatterer
320  G4cout << "size of FS: "<<products->size()<<G4endl;
321  #endif
322 
323  G4KineticTrack *final= products->operator[](0);
324 
325 
326  #ifdef debug_G4Scatterer
327  G4cout << " FinalState 1: "
328  << final->Get4Momentum()<< " "
329  << final->Get4Momentum().mag() ;
330  #endif
331 
332  if(products->size() == 1) return products;
333  final=products->operator[](1);
334  #ifdef debug_G4Scatterer
335  G4cout << ", 2: "
336  << final->Get4Momentum() << " "
337  << final->Get4Momentum().mag() << " " << G4endl;
338  #endif
339 
340  final= products->operator[](0);
341  G4LorentzVector pFinal=final->Get4Momentum();
342  if(products->size()==2)
343  {
344  final=products->operator[](1);
345  pFinal +=final->Get4Momentum();
346  }
347 
348  #ifdef debug_G4Scatterer
349  if ( (pInitial-pFinal).mag() > 0.1*MeV )
350  {
351  G4cout << "G4Scatterer: momentum imbalance, pInitial= " <<pInitial << " pFinal= " <<pFinal<< G4endl;
352  }
353  G4cout << "Scatterer costh= " << trk1.Get4Momentum().vect().unit() *(products->operator[](0))->Get4Momentum().vect().unit()<< G4endl;
354  #endif
355 
356  for(size_t hpw=0; hpw<products->size(); hpw++)
357  {
358  energyBalance-=products->operator[](hpw)->Get4Momentum().t();
359  pxBalance-=products->operator[](hpw)->Get4Momentum().vect().x();
360  pyBalance-=products->operator[](hpw)->Get4Momentum().vect().y();
361  pzBalance-=products->operator[](hpw)->Get4Momentum().vect().z();
362  chargeBalance-=G4lrint(products->operator[](hpw)->GetDefinition()->GetPDGCharge());
363  baryonBalance-=products->operator[](hpw)->GetDefinition()->GetBaryonNumber();
364  }
365  if(std::getenv("ScattererEnergyBalanceCheck"))
366  std::cout << "DEBUGGING energy balance A: "
367  <<energyBalance<<" "
368  <<pxBalance<<" "
369  <<pyBalance<<" "
370  <<pzBalance<<" "
371  <<chargeBalance<<" "
372  <<baryonBalance<<" "
373  <<G4endl;
374  if(chargeBalance !=0 )
375  {
376  G4cout << "track 1"<<trk1.GetDefinition()->GetParticleName()<<G4endl;
377  G4cout << "track 2"<<trk2.GetDefinition()->GetParticleName()<<G4endl;
378  for(size_t hpw=0; hpw<products->size(); hpw++)
379  {
380  G4cout << products->operator[](hpw)->GetDefinition()->GetParticleName()<<G4endl;
381  }
382  G4Exception("G4Scatterer", "im_r_matrix001", FatalException,
383  "Problem in ChargeBalance");
384  }
385  return products;
386  }
387  }
388 
389  return NULL;
390 }
391 
392 //----------------------------------------------------------------------------
393 
394 const G4VCollision* G4Scatterer::FindCollision(const G4KineticTrack& trk1,
395  const G4KineticTrack& trk2) const
396 {
397  G4VCollision* collisionInCharge = 0;
398 
399  size_t i;
400  for (i=0; i<collisions.size(); i++)
401  {
402  G4VCollision* component = collisions[i];
403  if (component->IsInCharge(trk1,trk2))
404  {
405  collisionInCharge = component;
406  break;
407  }
408  }
409 // if(collisionInCharge)
410 // {
411 // G4cout << "found collision : "
412 // << collisionInCharge->GetName()<< " "
413 // << "for "
414 // << trk1.GetDefinition()->GetParticleName()<<" + "
415 // << trk2.GetDefinition()->GetParticleName()<<" "
416 // << G4endl;;
417 // }
418  return collisionInCharge;
419 }
420 
421 //----------------------------------------------------------------------------
422 
423 G4double G4Scatterer::GetCrossSection(const G4KineticTrack& trk1,
424  const G4KineticTrack& trk2) const
425 {
426  const G4VCollision* collision = FindCollision(trk1,trk2);
427  G4double aCrossSection = 0;
428  if (collision != 0)
429  {
430  aCrossSection = collision->CrossSection(trk1,trk2);
431  }
432  return aCrossSection;
433 }
434 
435 //----------------------------------------------------------------------------
436 
437 const std::vector<G4CollisionInitialState *> & G4Scatterer::
438 GetCollisions(G4KineticTrack * aProjectile,
439  std::vector<G4KineticTrack *> & someCandidates,
440  G4double aCurrentTime)
441 {
442  theCollisions.clear();
443  std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
444  for(; j != someCandidates.end(); ++j)
445  {
446  G4double collisionTime = GetTimeToInteraction(*aProjectile, **j);
447  if(collisionTime == DBL_MAX) // no collision
448  {
449  continue;
450  }
451  G4KineticTrackVector aTarget;
452  aTarget.push_back(*j);
453  theCollisions.push_back(
454  new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) );
455 // G4cerr <<" !!!!!! debug collisions "<<collisionTime<<" "<<pkt->GetDefinition()->GetParticleName()<<G4endl;
456  }
457  return theCollisions;
458 }
459 
460 
461 G4KineticTrackVector * G4Scatterer::
462 GetFinalState(G4KineticTrack * aProjectile,
463  std::vector<G4KineticTrack *> & theTargets)
464 {
465  G4KineticTrack target_reloc(*(theTargets[0]));
466  return Scatter(*aProjectile, target_reloc);
467 }
468 //----------------------------------------------------------------------------