Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Public Member Functions | List of all members
G4Scatterer Class Reference

#include <G4Scatterer.hh>

Inheritance diagram for G4Scatterer:
Inheritance graph
[legend]
Collaboration diagram for G4Scatterer:
Collaboration graph
[legend]

Public Member Functions

 G4Scatterer ()
 
virtual ~G4Scatterer ()
 
virtual G4double GetTimeToInteraction (const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
 
G4double GetCrossSection (const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
 
virtual G4KineticTrackVector * Scatter (const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
 
virtual const std::vector
< G4CollisionInitialState * > & 		GetCollisions (G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &someCandidates, G4double aCurrentTime)
 
virtual G4KineticTrackVector * GetFinalState (G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &theTargets)
 
- Public Member Functions inherited from G4VScatterer
 G4VScatterer ()
 
virtual ~G4VScatterer ()
 
- Public Member Functions inherited from G4BCAction
 G4BCAction ()
 
virtual ~G4BCAction ()
 

Detailed Description

Definition at line 44 of file G4Scatterer.hh.

Constructor & Destructor Documentation

G4Scatterer::G4Scatterer ( )

Definition at line 63 of file G4Scatterer.cc.

64 {
65  G4AutoLock l(&collisions_mutex);
66  if ( ! setupDone )
67  {
68  Register aR;
69  G4ForEach<theChannels>::Apply(&aR, &collisions);
70  setupDone = true;
71  }
72 }
G4ForEach::Apply
static void Apply()
Definition: G4Pair.hh:179
G4AutoDelete::Register
void Register(T *inst)
Definition: G4AutoDelete.hh:65

Here is the call graph for this function:

G4Scatterer::~G4Scatterer ( )
virtual

Definition at line 76 of file G4Scatterer.cc.

77 {
78  G4AutoLock l(&collisions_mutex);
79  std::for_each(collisions.begin(), collisions.end(), G4Delete());
80  collisions.clear();
81 }

Member Function Documentation

const std::vector< G4CollisionInitialState * > & G4Scatterer::GetCollisions ( G4KineticTrack *  aProjectile,
std::vector< G4KineticTrack * > &  someCandidates,
G4double  aCurrentTime 
)
virtual

Implements G4BCAction.

Definition at line 439 of file G4Scatterer.cc.

442 {
443  theCollisions.clear();
444  std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
445  for(; j != someCandidates.end(); ++j)
446  {
447  G4double collisionTime = GetTimeToInteraction(*aProjectile, **j);
448  if(collisionTime == DBL_MAX) // no collision
449  {
450  continue;
451  }
452  G4KineticTrackVector aTarget;
453  aTarget.push_back(*j);
454  theCollisions.push_back(
455  new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) );
456 // G4cerr <<" !!!!!! debug collisions "<<collisionTime<<" "<<pkt->GetDefinition()->GetParticleName()<<G4endl;
457  }
458  return theCollisions;
459 }
G4Scatterer::GetTimeToInteraction
virtual G4double GetTimeToInteraction(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
Definition: G4Scatterer.cc:85
G4double
double G4double
Definition: G4Types.hh:76
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

G4double G4Scatterer::GetCrossSection ( const G4KineticTrack &  trk1,
const G4KineticTrack &  trk2 
) const

Definition at line 424 of file G4Scatterer.cc.

426 {
427  const G4VCollision* collision = FindCollision(trk1,trk2);
428  G4double aCrossSection = 0;
429  if (collision != 0)
430  {
431  aCrossSection = collision->CrossSection(trk1,trk2);
432  }
433  return aCrossSection;
434 }
G4VCollision::CrossSection
virtual G4double CrossSection(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
Definition: G4VCollision.cc:55
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4KineticTrackVector * G4Scatterer::GetFinalState ( G4KineticTrack *  aProjectile,
std::vector< G4KineticTrack * > &  theTargets 
)
virtual

Implements G4BCAction.

Definition at line 463 of file G4Scatterer.cc.

465 {
466  G4KineticTrack target_reloc(*(theTargets[0]));
467  return Scatter(*aProjectile, target_reloc);
468 }
G4Scatterer::Scatter
virtual G4KineticTrackVector * Scatter(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
Definition: G4Scatterer.cc:284

Here is the call graph for this function:

G4double G4Scatterer::GetTimeToInteraction ( const G4KineticTrack &  trk1,
const G4KineticTrack &  trk2 
) const
virtual

Implements G4VScatterer.

Definition at line 85 of file G4Scatterer.cc.

87 {
88  G4double time = DBL_MAX;
89  G4double distance_fast;
90  G4LorentzVector mom1 = trk1.GetTrackingMomentum();
91 // G4cout << "zcomp=" << std::abs(mom1.vect().unit().z() -1 ) << G4endl;
92  G4double collisionTime;
93 
94  if ( std::abs(mom1.vect().unit().z() -1 ) < 1e-6 )
95  {
96  G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
97  G4double deltaz=position.z();
98  G4double velocity = mom1.z()/mom1.e() * c_light;
99 
100  collisionTime=deltaz/velocity;
101  distance_fast=position.x()*position.x() + position.y()*position.y();
102  } else {
103 
104  // The nucleons of the nucleus are FROZEN, ie. do not move..
105 
106  G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
107 
108  G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light; // mom1.boostVector() will exit on slightly negative mass
109  collisionTime = (position * velocity) / velocity.mag2(); // can't divide by /c_light;
110  position -= velocity * collisionTime;
111  distance_fast=position.mag2();
112 
113 // if ( collisionTime>0 ) G4cout << " dis1/2 square" << dis1 <<" "<< dis2 << G4endl;
114 // collisionTime = GetTimeToClosestApproach(trk1,trk2);
115  }
116  if (collisionTime > 0)
117  {
118  static const G4double maxCrossSection = 500*millibarn;
119  if(0.7*pi*distance_fast>maxCrossSection) return time;
120 
121 
122  G4LorentzVector mom2(0,0,0,trk2.Get4Momentum().mag());
123 
124 // G4ThreeVector momLab = mom1.vect();// frozen Nucleus - mom2.vect();
125 // G4ThreeVector posLab = trk1.GetPosition() - trk2.GetPosition();
126 // G4double disLab=posLab * posLab - (posLab*momLab) * (posLab*momLab) /(momLab.mag2());
127 
128  G4LorentzRotation toCMSFrame((-1)*(mom1 + mom2).boostVector());
129  mom1 = toCMSFrame * mom1;
130  mom2 = toCMSFrame * mom2;
131 
132  G4LorentzVector coordinate1(trk1.GetPosition(), 100.);
133  G4LorentzVector coordinate2(trk2.GetPosition(), 100.);
134  G4ThreeVector pos = ((toCMSFrame * coordinate1).vect() -
135  (toCMSFrame * coordinate2).vect());
136 
137  G4ThreeVector mom = mom1.vect() - mom2.vect();
138 
139  // Calculate the impact parameter
140 
141  G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom.mag2());
142 
143 // G4cout << " disDiff " << distance-disLab << " " << disLab
144 // << " " << std::abs(distance-disLab)/distance << G4endl
145 // << " mom/Lab " << mom << " " << momLab << G4endl
146 // << " pos/Lab " << pos << " " << posLab
147 // << G4endl;
148 
149  if(pi*distance>maxCrossSection) return time;
150 
151  // charged particles special
152  static const G4double maxChargedCrossSection = 200*millibarn;
153  if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 &&
154  std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 &&
155  pi*distance>maxChargedCrossSection) return time;
156 
157  G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
158  // neutrons special pn is largest cross-section, but above 1.91 GeV is less than 200 mb
159  if(( trk1.GetDefinition() == G4Neutron::Neutron() ||
160  trk2.GetDefinition() == G4Neutron::Neutron() ) &&
161  sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time;
162 
163 /*
164  * if(distance <= sqr(1.14*fermi))
165  * {
166  * time = collisionTime;
167  *
168  * *
169  * * G4cout << "Scatter distance/time: " << std::sqrt(distance)/fermi <<
170  * * " / "<< time/ns << G4endl;
171  * * G4ThreeVector pos1=trk1.GetPosition();
172  * * G4ThreeVector pos2=trk2.GetPosition();
173  * * G4LorentzVector xmom1 = trk1.Get4Momentum();
174  * * G4LorentzVector xmom2 = trk2.Get4Momentum();
175  * * G4cout << "position1: " << pos1.x() << " " << pos1.y() << " "
176  * * << pos1.z();
177  * * pos1+=(collisionTime*c_light/xmom1.e())*xmom1.vect();
178  * * G4cout << " straight line trprt: "
179  * * << pos1.x() << " " << pos1.y() << " "
180  * * << pos1.z() << G4endl;
181  * * G4cout << "position2: " << pos2.x() << " " << pos2.y() << " "
182  * * << pos2.z() << G4endl;
183  * * G4cout << "straight line distance 2 fixed:" << (pos1-pos2).mag()/fermi << G4endl;
184  * * pos2+= (collisionTime*c_light/xmom2.e())*xmom2.vect();
185  * * G4cout<< " straight line trprt: "
186  * * << pos2.x() << " " << pos2.y() << " "
187  * * << pos2.z() << G4endl;
188  * * G4cout << "straight line distance :" << (pos1-pos2).mag()/fermi << G4endl;
189  * *
190  * }
191  *
192  * if(1)
193  * return time;
194  */
195 
196  if ((trk1.GetActualMass()+trk2.GetActualMass()) > sqrtS) return time;
197 
198 
199 
200  const G4VCollision* collision = FindCollision(trk1,trk2);
201  G4double totalCrossSection;
202  // The cross section is interpreted geometrically as an area
203  // Two particles are assumed to collide if their distance is < (totalCrossSection/pi)
204 
205  if (collision != 0)
206  {
207  totalCrossSection = collision->CrossSection(trk1,trk2);
208  if ( totalCrossSection > 0 )
209  {
210 /* G4cout << " totalCrossection = "<< totalCrossSection << ", trk1/2, s, e-m: "
211  * << trk1.GetDefinition()->GetParticleName()
212  * << " / "
213  * << trk2.GetDefinition()->GetParticleName()
214  * << ", "
215  * << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
216  * << ", "
217  * << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
218  * trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
219  * << G4endl;
220  */
221  if (distance <= totalCrossSection / pi)
222  {
223  time = collisionTime;
224  }
225  } else
226  {
227 
228  // For debugging...
229  // G4cout << " totalCrossection = 0, trk1/2, s, e-m: "
230  // << trk1.GetDefinition()->GetParticleName()
231  // << " / "
232  // << trk2.GetDefinition()->GetParticleName()
233  // << ", "
234  // << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
235  // << ", "
236  // << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
237  // trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
238  // << G4endl;
239 
240  }
241 /*
242  * if(distance <= sqr(5.*fermi))
243  * {
244  * G4cout << " distance,xsect, std::sqrt(xsect/pi) : " << std::sqrt(distance)/fermi
245  * << " " << totalCrossSection/sqr(fermi)
246  * << " " << std::sqrt(totalCrossSection / pi)/fermi << G4endl;
247  * }
248  */
249 
250  }
251  else
252  {
253  time = DBL_MAX;
254 // /*
255  // For debugging
256 //hpw G4cout << "G4Scatterer - collision not found: "
257 //hpw << trk1.GetDefinition()->GetParticleName()
258 //hpw << " - "
259 //hpw << trk2.GetDefinition()->GetParticleName()
260 //hpw << G4endl;
261  // End of debugging
262 // */
263  }
264  }
265 
266  else
267  {
268  /*
269  // For debugging
270  G4cout << "G4Scatterer - negative collisionTime"
271  << ": collisionTime = " << collisionTime
272  << ", position = " << position
273  << ", velocity = " << velocity
274  << G4endl;
275  // End of debugging
276  */
277  }
278 
279  return time;
280 }
CLHEP::Hep3Vector::x
double x() const
G4KineticTrack::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4KineticTrack.hh:230
G4KineticTrack::GetActualMass
G4double GetActualMass() const
Definition: G4KineticTrack.hh:319
CLHEP::Hep3Vector::z
double z() const
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
CLHEP::HepLorentzVector::mag
double mag() const
G4VCollision::CrossSection
virtual G4double CrossSection(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
Definition: G4VCollision.cc:55
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4KineticTrack::GetTrackingMomentum
const G4LorentzVector & GetTrackingMomentum() const
Definition: G4KineticTrack.hh:245
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::y
double y() const
GeV
static constexpr double GeV
Definition: G4SIunits.hh:217
CLHEP::Hep3Vector::mag2
double mag2() const
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:124
G4KineticTrack::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4KineticTrack.hh:240
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4KineticTrack::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4KineticTrack.hh:210
python.hepunit.c_light
float c_light
Definition: hepunit.py:257
pos
static const G4double pos
Definition: G4ElectroNuclearCrossSection.cc:66
millibarn
static constexpr double millibarn
Definition: G4SIunits.hh:106

Here is the call graph for this function:

Here is the caller graph for this function:

G4KineticTrackVector * G4Scatterer::Scatter ( const G4KineticTrack &  trk1,
const G4KineticTrack &  trk2 
) const
virtual

Implements G4VScatterer.

Definition at line 284 of file G4Scatterer.cc.

286 {
287 // G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
288  G4LorentzVector pInitial=trk1.Get4Momentum() + trk2.Get4Momentum();
289  G4double energyBalance = pInitial.t();
290  G4double pxBalance = pInitial.vect().x();
291  G4double pyBalance = pInitial.vect().y();
292  G4double pzBalance = pInitial.vect().z();
293  G4int chargeBalance = G4lrint(trk1.GetDefinition()->GetPDGCharge()
294  + trk2.GetDefinition()->GetPDGCharge());
295  G4int baryonBalance = trk1.GetDefinition()->GetBaryonNumber()
296  + trk2.GetDefinition()->GetBaryonNumber();
297 
298  const G4VCollision* collision = FindCollision(trk1,trk2);
299  if (collision != 0)
300  {
301  G4double aCrossSection = collision->CrossSection(trk1,trk2);
302  if (aCrossSection > 0.0)
303  {
304 
305 
306  #ifdef debug_G4Scatterer
307  G4cout << "be4 FinalState 1(p,e,m): "
308  << trk1.Get4Momentum() << " "
309  << trk1.Get4Momentum().mag()
310  << ", 2: "
311  << trk2.Get4Momentum()<< " "
312  << trk2.Get4Momentum().mag() << " "
313  << G4endl;
314  #endif
315 
316 
317  G4KineticTrackVector* products = collision->FinalState(trk1,trk2);
318  if(!products || products->size() == 0) return products;
319 
320  #ifdef debug_G4Scatterer
321  G4cout << "size of FS: "<<products->size()<<G4endl;
322  #endif
323 
324  G4KineticTrack *final= products->operator[](0);
325 
326 
327  #ifdef debug_G4Scatterer
328  G4cout << " FinalState 1: "
329  << final->Get4Momentum()<< " "
330  << final->Get4Momentum().mag() ;
331  #endif
332 
333  if(products->size() == 1) return products;
334  final=products->operator[](1);
335  #ifdef debug_G4Scatterer
336  G4cout << ", 2: "
337  << final->Get4Momentum() << " "
338  << final->Get4Momentum().mag() << " " << G4endl;
339  #endif
340 
341  final= products->operator[](0);
342  G4LorentzVector pFinal=final->Get4Momentum();
343  if(products->size()==2)
344  {
345  final=products->operator[](1);
346  pFinal +=final->Get4Momentum();
347  }
348 
349  #ifdef debug_G4Scatterer
350  if ( (pInitial-pFinal).mag() > 0.1*MeV )
351  {
352  G4cout << "G4Scatterer: momentum imbalance, pInitial= " <<pInitial << " pFinal= " <<pFinal<< G4endl;
353  }
354  G4cout << "Scatterer costh= " << trk1.Get4Momentum().vect().unit() *(products->operator[](0))->Get4Momentum().vect().unit()<< G4endl;
355  #endif
356 
357  for(size_t hpw=0; hpw<products->size(); hpw++)
358  {
359  energyBalance-=products->operator[](hpw)->Get4Momentum().t();
360  pxBalance-=products->operator[](hpw)->Get4Momentum().vect().x();
361  pyBalance-=products->operator[](hpw)->Get4Momentum().vect().y();
362  pzBalance-=products->operator[](hpw)->Get4Momentum().vect().z();
363  chargeBalance-=G4lrint(products->operator[](hpw)->GetDefinition()->GetPDGCharge());
364  baryonBalance-=products->operator[](hpw)->GetDefinition()->GetBaryonNumber();
365  }
366  if(getenv("ScattererEnergyBalanceCheck"))
367  std::cout << "DEBUGGING energy balance A: "
368  <<energyBalance<<" "
369  <<pxBalance<<" "
370  <<pyBalance<<" "
371  <<pzBalance<<" "
372  <<chargeBalance<<" "
373  <<baryonBalance<<" "
374  <<G4endl;
375  if(chargeBalance !=0 )
376  {
377  G4cout << "track 1"<<trk1.GetDefinition()->GetParticleName()<<G4endl;
378  G4cout << "track 2"<<trk2.GetDefinition()->GetParticleName()<<G4endl;
379  for(size_t hpw=0; hpw<products->size(); hpw++)
380  {
381  G4cout << products->operator[](hpw)->GetDefinition()->GetParticleName()<<G4endl;
382  }
383  G4Exception("G4Scatterer", "im_r_matrix001", FatalException,
384  "Problem in ChargeBalance");
385  }
386  return products;
387  }
388  }
389 
390  return NULL;
391 }
CLHEP::Hep3Vector::x
double x() const
G4int
int G4int
Definition: G4Types.hh:78
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:120
CLHEP::Hep3Vector::z
double z() const
CLHEP::HepLorentzVector::vect
Hep3Vector vect() const
G4cout
G4GLOB_DLL std::ostream G4cout
CLHEP::HepLorentzVector::mag
double mag() const
G4VCollision::CrossSection
virtual G4double CrossSection(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const
Definition: G4VCollision.cc:55
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4lrint
int G4lrint(double ad)
Definition: templates.hh:163
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::y
double y() const
FatalException
Definition: G4ExceptionSeverity.hh:60
G4VCollision::FinalState
virtual G4KineticTrackVector * FinalState(const G4KineticTrack &trk1, const G4KineticTrack &trk2) const =0
G4endl
#define G4endl
Definition: G4ios.hh:61
MeV
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:124
G4KineticTrack::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4KineticTrack.hh:240
G4KineticTrack::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4KineticTrack.hh:210

Here is the call graph for this function:

Here is the caller graph for this function:

The documentation for this class was generated from the following files: