Geant4_10
G4MesonAbsorption.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 
27 #include "G4MesonAbsorption.hh"
28 #include "G4PhysicalConstants.hh"
29 #include "G4SystemOfUnits.hh"
30 #include "G4LorentzRotation.hh"
31 #include "G4LorentzVector.hh"
32 #include "Randomize.hh"
33 #include "G4KineticTrackVector.hh"
35 #include <cmath>
36 #include "G4PionPlus.hh"
37 #include "G4PionMinus.hh"
38 #include "G4ParticleDefinition.hh"
39 #include "G4HadTmpUtil.hh"
40 
41 // first prototype
42 
43 const std::vector<G4CollisionInitialState *> & G4MesonAbsorption::
45  std::vector<G4KineticTrack *> & someCandidates,
46  G4double aCurrentTime)
47 {
48  theCollisions.clear();
49  if(someCandidates.size() >1)
50  {
51  std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
52  for(; j != someCandidates.end(); ++j)
53  {
54  G4double collisionTime = GetTimeToAbsorption(*aProjectile, **j);
55  if(collisionTime == DBL_MAX)
56  {
57  continue;
58  }
59  G4KineticTrackVector aTarget;
60  aTarget.push_back(*j);
61  FindAndFillCluster(aTarget, aProjectile, someCandidates);
62  if(aTarget.size()>=2)
63  {
64  theCollisions.push_back(
65  new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) );
66  }
67  }
68  }
69  return theCollisions;
70 }
71 
72 
73 void G4MesonAbsorption::
74 FindAndFillCluster(G4KineticTrackVector & result,
75  G4KineticTrack * aProjectile, std::vector<G4KineticTrack *> & someCandidates)
76 {
77  std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
78  G4KineticTrack * aTarget = result[0];
79  G4int chargeSum = G4lrint(aTarget->GetDefinition()->GetPDGCharge());
80  chargeSum+=G4lrint(aProjectile->GetDefinition()->GetPDGCharge());
81  G4ThreeVector firstBase = aTarget->GetPosition();
83  G4KineticTrack * partner = 0;
84  for(; j != someCandidates.end(); ++j)
85  {
86  if(*j == aTarget) continue;
87  G4int cCharge = G4lrint((*j)->GetDefinition()->GetPDGCharge());
88  if (chargeSum+cCharge > 2) continue;
89  if (chargeSum+cCharge < 0) continue;
90  // get the one with the smallest distance.
91  G4ThreeVector secodeBase = (*j)->GetPosition();
92  if((firstBase+secodeBase).mag()<min)
93  {
94  min=(firstBase+secodeBase).mag();
95  partner = *j;
96  }
97  }
98  if(partner) result.push_back(partner);
99  else result.clear();
100 }
101 
102 
105  std::vector<G4KineticTrack *> & targets)
106 {
107  // G4cout << "We have an absorption !!!!!!!!!!!!!!!!!!!!!!"<<G4endl;
108  // Only 2-body absorption for the time being.
109  // If insufficient, use 2-body absorption and clusterization to emulate 3-,4-,etc body absorption.
110  G4LorentzVector thePro = projectile->Get4Momentum();
111  G4LorentzVector theT1 = targets[0]->Get4Momentum();
112  G4LorentzVector theT2 = targets[1]->Get4Momentum();
113  G4LorentzVector incoming = thePro + theT1 + theT2;
114  G4double energyBalance = incoming.t();
115  G4int chargeBalance = G4lrint(projectile->GetDefinition()->GetPDGCharge()
116  + targets[0]->GetDefinition()->GetPDGCharge()
117  + targets[1]->GetDefinition()->GetPDGCharge());
118 
119  G4int baryonBalance = projectile->GetDefinition()->GetBaryonNumber()
120  + targets[0]->GetDefinition()->GetBaryonNumber()
121  + targets[1]->GetDefinition()->GetBaryonNumber();
122 
123 
124  // boost all to MMS
125  G4LorentzRotation toSPS((-1)*(thePro + theT1 + theT2).boostVector());
126  theT1 = toSPS * theT1;
127  theT2 = toSPS * theT2;
128  thePro = toSPS * thePro;
129  G4LorentzRotation fromSPS(toSPS.inverse());
130 
131  // rotate to z
132  G4LorentzRotation toZ;
133  G4LorentzVector Ptmp=projectile->Get4Momentum();
134  toZ.rotateZ(-1*Ptmp.phi());
135  toZ.rotateY(-1*Ptmp.theta());
136  theT1 = toZ * theT1;
137  theT2 = toZ * theT2;
138  thePro = toZ * thePro;
139  G4LorentzRotation toLab(toZ.inverse());
140 
141  // Get definitions
142  G4ParticleDefinition * d1 = targets[0]->GetDefinition();
143  G4ParticleDefinition * d2 = targets[1]->GetDefinition();
144  if(0.5>G4UniformRand())
145  {
146  G4ParticleDefinition * temp;
147  temp=d1;d1=d2;d2=temp;
148  }
149  G4ParticleDefinition * dp = projectile->GetDefinition();
150  if(dp->GetPDGCharge()<-.5)
151  {
152  if(d1->GetPDGCharge()>.5)
153  {
154  if(d2->GetPDGCharge()>.5 && 0.5>G4UniformRand())
155  {
157  }
158  else
159  {
161  }
162  }
163  else if(d2->GetPDGCharge()>.5)
164  {
166  }
167  }
168  else if(dp->GetPDGCharge()>.5)
169  {
170  if(d1->GetPDGCharge()<.5)
171  {
172  if(d2->GetPDGCharge()<.5 && 0.5>G4UniformRand())
173  {
175  }
176  else
177  {
179  }
180  }
181  else if(d2->GetPDGCharge()<.5)
182  {
184  }
185  }
186 
187  // calculate the momenta.
188  G4double M_sq = (thePro+theT1+theT2).mag2();
189  G4double m1_sq = sqr(d1->GetPDGMass());
190  G4double m2_sq = sqr(d2->GetPDGMass());
191  G4double m_sq = M_sq-m1_sq-m2_sq;
192  G4double p = std::sqrt((m_sq*m_sq - 4.*m1_sq * m2_sq)/(4.*M_sq));
193  G4double costh = 2.*G4UniformRand()-1.;
194  G4double phi = 2.*pi*G4UniformRand();
195  G4ThreeVector pFinal(p*std::sin(std::acos(costh))*std::cos(phi), p*std::sin(std::acos(costh))*std::sin(phi), p*costh);
196 
197  // G4cout << "testing p "<<p-pFinal.mag()<<G4endl;
198  // construct the final state particles lorentz momentum.
199  G4double eFinal1 = std::sqrt(m1_sq+pFinal.mag2());
200  G4LorentzVector final1(pFinal, eFinal1);
201  G4double eFinal2 = std::sqrt(m2_sq+pFinal.mag2());
202  G4LorentzVector final2(-1.*pFinal, eFinal2);
203 
204  // rotate back.
205  final1 = toLab * final1;
206  final2 = toLab * final2;
207 
208  // boost back to LAB
209  final1 = fromSPS * final1;
210  final2 = fromSPS * final2;
211 
212  // make the final state
213  G4KineticTrack * f1 =
214  new G4KineticTrack(d1, 0., targets[0]->GetPosition(), final1);
215  G4KineticTrack * f2 =
216  new G4KineticTrack(d2, 0., targets[1]->GetPosition(), final2);
218  result->push_back(f1);
219  result->push_back(f2);
220 
221  for(size_t hpw=0; hpw<result->size(); hpw++)
222  {
223  energyBalance-=result->operator[](hpw)->Get4Momentum().t();
224  chargeBalance-=G4lrint(result->operator[](hpw)->GetDefinition()->GetPDGCharge());
225  baryonBalance-=result->operator[](hpw)->GetDefinition()->GetBaryonNumber();
226  }
227  if(getenv("AbsorptionEnergyBalanceCheck"))
228  std::cout << "DEBUGGING energy balance B: "
229  <<energyBalance<<" "
230  <<chargeBalance<<" "
231  <<baryonBalance<<" "
232  <<G4endl;
233  return result;
234 }
235 
236 
237 G4double G4MesonAbsorption::
238 GetTimeToAbsorption(const G4KineticTrack& trk1, const G4KineticTrack& trk2)
239 {
244  {
245  return DBL_MAX;
246  }
247  G4double time = DBL_MAX;
248  G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
249 
250  // Check whether there is enough energy for elastic scattering
251  // (to put the particles on to mass shell
252 
253  if (trk1.GetActualMass() + trk2.GetActualMass() < sqrtS)
254  {
255  G4LorentzVector mom1 = trk1.GetTrackingMomentum();
256  G4ThreeVector position = trk1.GetPosition() - trk2.GetPosition();
257  if ( mom1.mag2() < -1.*eV )
258  {
259  G4cout << "G4MesonAbsorption::GetTimeToInteraction(): negative m2:" << mom1.mag2() << G4endl;
260  }
261  G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light;
262  G4double collisionTime = - (position * velocity) / (velocity * velocity);
263 
264  if (collisionTime > 0)
265  {
266  G4LorentzVector mom2(0,0,0,trk2.Get4Momentum().mag());
267  G4LorentzRotation toCMSFrame((-1)*(mom1 + mom2).boostVector());
268  mom1 = toCMSFrame * mom1;
269  mom2 = toCMSFrame * mom2;
270 
271  G4LorentzVector coordinate1(trk1.GetPosition(), 100.);
272  G4LorentzVector coordinate2(trk2.GetPosition(), 100.);
273  G4ThreeVector pos = ((toCMSFrame * coordinate1).vect() -
274  (toCMSFrame * coordinate2).vect());
275  G4ThreeVector mom = mom1.vect() - mom2.vect();
276 
277  G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom*mom);
278 
279  // global optimization
280  static const G4double maxCrossSection = 500*millibarn;
281  if(pi*distance>maxCrossSection) return time;
282  // charged particles special optimization
283  static const G4double maxChargedCrossSection = 200*millibarn;
284  if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 &&
285  std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 &&
286  pi*distance>maxChargedCrossSection) return time;
287  // neutrons special optimization
288  if(( trk1.GetDefinition() == G4Neutron::Neutron() ||
289  trk2.GetDefinition() == G4Neutron::Neutron() ) &&
290  sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time;
291 
292  G4double totalCrossSection = AbsorptionCrossSection(trk1,trk2);
293  if ( totalCrossSection > 0 )
294  {
295  if (distance <= totalCrossSection / pi)
296  {
297  time = collisionTime;
298  }
299  }
300  }
301  }
302  return time;
303 }
304 
305 G4double G4MesonAbsorption::
306 AbsorptionCrossSection(const G4KineticTrack & aT, const G4KineticTrack & bT)
307 {
308  G4double t = 0;
311  {
312  t = aT.Get4Momentum().t()-aT.Get4Momentum().mag()/MeV;
313  }
316  {
317  t = bT.Get4Momentum().t()-bT.Get4Momentum().mag()/MeV;
318  }
319 
320  static const G4double it [26] =
321  {0,4,50,5.5,75,8,95,10,120,11.5,140,12,160,11.5,180,10,190,8,210,6,235,4,260,3,300,2};
322 
323  G4double aCross(0);
324  if(t<=it[24])
325  {
326  G4int count = 0;
327  while(t>it[count])count+=2;
328  G4double x1 = it[count-2];
329  G4double x2 = it[count];
330  G4double y1 = it[count-1];
331  G4double y2 = it[count+1];
332  aCross = y1+(y2-y1)/(x2-x1)*(t-x1);
333  // G4cout << "Printing the absorption crosssection "
334  // <<x1<< " "<<x2<<" "<<t<<" "<<y1<<" "<<y2<<" "<<0.5*aCross<<G4endl;
335  }
336  return .5*aCross*millibarn;
337 }
Double_t y2[nxs]
Definition: Style.C:21
Double_t y1[nxs]
Definition: Style.C:20
Double_t x2[nxs]
Definition: Style.C:19
const G4ThreeVector & GetPosition() const
const char * p
Definition: xmltok.h:285
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
G4double G4NeutronHPJENDLHEData::G4double result
G4double GetActualMass() const
HepLorentzVector & rotateZ(double)
int G4int
Definition: G4Types.hh:78
int millibarn
Definition: hepunit.py:40
double phi() const
Float_t f1
G4ParticleDefinition * GetDefinition() const
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
Float_t f2
double theta() const
double mag() const
static G4PionMinus * PionMinusDefinition()
Definition: G4PionMinus.cc:93
Double_t x1[nxs]
Definition: Style.C:18
static G4PionPlus * PionPlusDefinition()
Definition: G4PionPlus.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetPDGMass() const
int G4lrint(double ad)
Definition: templates.hh:163
const G4LorentzVector & GetTrackingMomentum() const
double mag2() const
HepLorentzVector & rotateY(double)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double mag2() const
virtual G4KineticTrackVector * GetFinalState(G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &theTargets)
#define G4endl
Definition: G4ios.hh:61
HepLorentzRotation inverse() const
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
const G4LorentzVector & Get4Momentum() const
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
virtual const std::vector< G4CollisionInitialState * > & GetCollisions(G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &someCandidates, G4double aCurrentTime)
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:99
float c_light
Definition: hepunit.py:257