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