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