| Geant4 Cross Reference |
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"
34 #include "G4CollisionInitialState.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::
44 GetCollisions(G4KineticTrack * aProjectile,
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();
82 G4double min = DBL_MAX;
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
103 G4KineticTrackVector * G4MesonAbsorption::
104 GetFinalState(G4KineticTrack * projectile,
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 const G4ParticleDefinition * d1 = targets[0]->GetDefinition();
143 const G4ParticleDefinition * d2 = targets[1]->GetDefinition();
144 if(0.5>G4UniformRand())
145 {
146 const G4ParticleDefinition * temp;
147 temp=d1;d1=d2;d2=temp;
148 }
149 const 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 {
156 d2 = G4Neutron::NeutronDefinition();
157 }
158 else
159 {
160 d1 = G4Neutron::NeutronDefinition();
161 }
162 }
163 else if(d2->GetPDGCharge()>.5)
164 {
165 d2 = G4Neutron::NeutronDefinition();
166 }
167 }
168 else if(dp->GetPDGCharge()>.5)
169 {
170 if(d1->GetPDGCharge()<.5)
171 {
172 if(d2->GetPDGCharge()<.5 && 0.5>G4UniformRand())
173 {
174 d2 = G4Proton::ProtonDefinition();
175 }
176 else
177 {
178 d1 = G4Proton::ProtonDefinition();
179 }
180 }
181 else if(d2->GetPDGCharge()<.5)
182 {
183 d2 = G4Proton::ProtonDefinition();
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);
217 G4KineticTrackVector * result = new G4KineticTrackVector;
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(std::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 {
240 if(trk1.GetDefinition()!=G4PionPlus::PionPlusDefinition() &&
241 trk1.GetDefinition()!=G4PionMinus::PionMinusDefinition() &&
242 trk2.GetDefinition()!=G4PionPlus::PionPlusDefinition() &&
243 trk2.GetDefinition()!=G4PionMinus::PionMinusDefinition())
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;
309 if(aT.GetDefinition()==G4PionPlus::PionPlusDefinition() ||
310 aT.GetDefinition()==G4PionMinus::PionMinusDefinition() )
311 {
312 t = aT.Get4Momentum().t()-aT.Get4Momentum().mag()/MeV;
313 }
314 else if(bT.GetDefinition()==G4PionPlus::PionPlusDefinition() ||
315 bT.GetDefinition()!=G4PionMinus::PionMinusDefinition())
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; /* Loop checking, 30-Oct-2015, G.Folger */
328
329 G4double x1 = it[count-2];
330 G4double x2 = it[count];
331 G4double y1 = it[count-1];
332 G4double y2 = it[count+1];
333 aCross = y1+(y2-y1)/(x2-x1)*(t-x1);
334 // G4cout << "Printing the absorption crosssection "
335 // <<x1<< " "<<x2<<" "<<t<<" "<<y1<<" "<<y2<<" "<<0.5*aCross<<G4endl;
336 }
337 return .5*aCross*millibarn;
338 }
339