| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27
28 #include <vector>
29
30 #include "globals.hh"
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4ios.hh"
34 #include "G4Scatterer.hh"
35 #include "G4KineticTrack.hh"
36 #include "G4ThreeVector.hh"
37 #include "G4LorentzRotation.hh"
38 #include "G4LorentzVector.hh"
39
40 #include "G4CollisionNN.hh"
41 #include "G4CollisionPN.hh"
42 #include "G4CollisionMesonBaryon.hh"
43
44 #include "G4CollisionInitialState.hh"
45 #include "G4HadTmpUtil.hh"
46 #include "G4Pair.hh"
47 #include "G4AutoLock.hh"
48
49 //Mutex for control of shared resource
50 namespace {
51 G4Mutex collisions_mutex = G4MUTEX_INITIAL
52 G4bool setupDone = false;
53 }
54
55 // Declare the categories of collisions the Sc
56 typedef GROUP2(G4CollisionNN, G4CollisionMeson
57
58 G4CollisionVector G4Scatterer::collisions;
59
60 //--------------------------------------------
61
62 G4Scatterer::G4Scatterer()
63 {
64 G4AutoLock l(&collisions_mutex);
65 if ( ! setupDone )
66 {
67 Register aR;
68 G4ForEach<theChannels>::Apply(&aR, &coll
69 setupDone = true;
70 }
71 }
72
73 //--------------------------------------------
74
75 G4Scatterer::~G4Scatterer()
76 {
77 G4AutoLock l(&collisions_mutex);
78 std::for_each(collisions.begin(), collisions
79 collisions.clear();
80 }
81
82 //--------------------------------------------
83
84 G4double G4Scatterer::GetTimeToInteraction(con
85 const G4KineticTrack& trk2) const
86 {
87 G4double time = DBL_MAX;
88 G4double distance_fast;
89 G4LorentzVector mom1 = trk1.GetTrackingMomen
90 // G4cout << "zcomp=" << std::abs(mom1.vect()
91 G4double collisionTime;
92
93 if ( std::abs(mom1.vect().unit().z() -1 ) <
94 {
95 G4ThreeVector position = trk2.GetPosition
96 G4double deltaz=position.z();
97 G4double velocity = mom1.z()/mom1.e() * c
98
99 collisionTime=deltaz/velocity;
100 distance_fast=position.x()*position.x() +
101 } else {
102
103 // The nucleons of the nucleus are FROZEN
104
105 G4ThreeVector position = trk2.GetPosition(
106
107 G4ThreeVector velocity = mom1.vect()/mom1.
108 collisionTime = (position * velocity) / ve
109 position -= velocity * collisionTime;
110 distance_fast=position.mag2();
111
112 // if ( collisionTime>0 ) G4cout << " dis1/
113 // collisionTime = GetTimeToClosestApproac
114 }
115 if (collisionTime > 0)
116 {
117 static const G4double maxCrossSection = 5
118 if(0.7*pi*distance_fast>maxCrossSection)
119
120
121 G4LorentzVector mom2(0,0,0,trk2.Get
122
123 // G4ThreeVector momLab = mom1.vect();// f
124 // G4ThreeVector posLab = trk1.GetPosition
125 // G4double disLab=posLab * posLab - (posL
126
127 G4LorentzRotation toCMSFrame((-1)*(mom1 +
128 mom1 = toCMSFrame * mom1;
129 mom2 = toCMSFrame * mom2;
130
131 G4LorentzVector coordinate1(trk1.GetPosit
132 G4LorentzVector coordinate2(trk2.GetPosit
133 G4ThreeVector pos = ((toCMSFrame * coordi
134 (toCMSFrame * coordinate2).vect());
135
136 G4ThreeVector mom = mom1.vect() - mom2.ve
137
138 // Calculate the impact parameter
139
140 G4double distance = pos * pos - (pos*mom)
141
142 // G4cout << " disDiff " << distance-disLa
143 // << " " << std::abs(distance-disL
144 // << " mom/Lab " << mom << " " << momLab
145 // << " pos/Lab " << pos << " " << posLab
146 // << G4endl;
147
148 if(pi*distance>maxCrossSection) return ti
149
150 // charged particles special
151 static const G4double maxChargedCrossSect
152 if(std::abs(trk1.GetDefinition()->GetPDGC
153 std::abs(trk2.GetDefinition()->GetPDGC
154 pi*distance>maxChargedCrossSection) re
155
156 G4double sqrtS = (trk1.Get4Momentum
157 // neutrons special pn is largest cross-
158 if(( trk1.GetDefinition() == G4Neutron::N
159 trk2.GetDefinition() == G4Neutron::N
160 sqrtS>1.91*GeV && pi*distance>maxChar
161
162 /*
163 * if(distance <= sqr(1.14*fermi))
164 * {
165 * time = collisionTime;
166 *
167 * *
168 * * G4cout << "Scatter distance/time:
169 * * " / "<< time/ns << G4endl;
170 * * G4ThreeVector pos1=trk1.GetPosit
171 * * G4ThreeVector pos2=trk2.GetPosit
172 * * G4LorentzVector xmom1 = trk1.Get
173 * * G4LorentzVector xmom2 = trk2.Get
174 * * G4cout << "position1: " << pos1
175 * * << pos1.z();
176 * * pos1+=(collisionTime*c_light/xmo
177 * * G4cout << " straight line trprt:
178 * * << pos1.x() << " " << pos1.
179 * * << pos1.z() << G4endl;
180 * * G4cout << "position2: " << pos2
181 * * << pos2.z() << G4endl;
182 * * G4cout << "straight line distanc
183 * * pos2+= (collisionTime*c_light/xm
184 * * G4cout<< " straight line trprt:
185 * * << pos2.x() << " " << pos2.
186 * * << pos2.z() << G4endl;
187 * * G4cout << "straight line distanc
188 * *
189 * }
190 *
191 * if(1)
192 * return time;
193 */
194
195 if ((trk1.GetActualMass()+trk2.GetActualM
196
197
198
199 const G4VCollision* collision = FindCollis
200 G4double totalCrossSection;
201 // The cross section is interpreted geomet
202 // Two particles are assumed to collide if
203
204 if (collision != 0)
205 {
206 totalCrossSection = collision->CrossSe
207 if ( totalCrossSection > 0 )
208 {
209 /* G4cout << " totalCrossection = "<< t
210 * << trk1.GetDefinition()->GetP
211 * << " / "
212 * << trk2.GetDefinition()->GetP
213 * << ", "
214 * << (trk1.Get4Momentum()+trk2.Get4Mo
215 * << ", "
216 * << (trk1.Get4Momentum()+trk2.Get4Mo
217 * trk1.Get4Momentum().mag() - trk
218 * << G4endl;
219 */
220 if (distance <= totalCrossSection / pi)
221 {
222 time = collisionTime;
223 }
224 } else
225 {
226
227 // For debugging...
228 // G4cout << " totalCrossection = 0, tr
229 // << trk1.GetDefinition()->GetP
230 // << " / "
231 // << trk2.GetDefinition()->GetP
232 // << ", "
233 // << (trk1.Get4Momentum()+trk2.Get4Mo
234 // << ", "
235 // << (trk1.Get4Momentum()+trk2.Get4Mo
236 // trk1.Get4Momentum().mag() - trk
237 // << G4endl;
238
239 }
240 /*
241 * if(distance <= sqr(5.*fermi))
242 * {
243 * G4cout << " distance,xsect, std::sqrt(
244 * << " " << totalCrossSection/sqr(fermi
245 * << " " << std::sqrt(totalCrossSection
246 * }
247 */
248
249 }
250 else
251 {
252 time = DBL_MAX;
253 // /*
254 // For debugging
255 //hpw G4cout << "G4Scatterer - c
256 //hpw << trk1.GetDefinition()->GetParti
257 //hpw << " - "
258 //hpw << trk2.GetDefinition()->GetParti
259 //hpw << G4endl;
260 // End of debugging
261 // */
262 }
263 }
264
265 else
266 {
267 /*
268 // For debugging
269 G4cout << "G4Scatterer - negative collisio
270 << ": collisionTime = " << collisionTime
271 << ", position = " << position
272 << ", velocity = " << velocity
273 << G4endl;
274 // End of debugging
275 */
276 }
277
278 return time;
279 }
280
281 //--------------------------------------------
282
283 G4KineticTrackVector* G4Scatterer::Scatter(con
284 const G4KineticTrack& trk2) co
285 {
286 // G4double sqrtS = (trk1.Get4Momentum() + t
287 G4LorentzVector pInitial=trk1.Get4Momentum(
288 G4double energyBalance = pInitial.t();
289 G4double pxBalance = pInitial.vect().x();
290 G4double pyBalance = pInitial.vect().y();
291 G4double pzBalance = pInitial.vect().z();
292 G4int chargeBalance = G4lrint(trk1.GetDefin
293 + trk2.GetDefinition()-
294 G4int baryonBalance = trk1.GetDefinition()-
295 + trk2.GetDefinition()-
296
297 const G4VCollision* collision = FindCollisi
298 if (collision != 0)
299 {
300 G4double aCrossSection = collision->Cross
301 if (aCrossSection > 0.0)
302 {
303
304
305 #ifdef debug_G4Scatterer
306 G4cout << "be4 FinalState 1(p,e,m): "
307 << trk1.Get4Momentum() << " "
308 << trk1.Get4Momentum().mag()
309 << ", 2: "
310 << trk2.Get4Momentum()<< " "
311 << trk2.Get4Momentum().mag() << " "
312 << G4endl;
313 #endif
314
315
316 G4KineticTrackVector* products = collis
317 if(!products || products->size() == 0)
318
319 #ifdef debug_G4Scatterer
320 G4cout << "size of FS: "<<products->siz
321 #endif
322
323 G4KineticTrack *final= products->operat
324
325
326 #ifdef debug_G4Scatterer
327 G4cout << " FinalState 1: "
328 << final->Get4Momentum()<< " "
329 << final->Get4Momentum().mag() ;
330 #endif
331
332 if(products->size() == 1) return produ
333 final=products->operator[](1);
334 #ifdef debug_G4Scatterer
335 G4cout << ", 2: "
336 << final->Get4Momentum() << " "
337 << final->Get4Momentum().mag() << "
338 #endif
339
340 final= products->operator[](0);
341 G4LorentzVector pFinal=final->Get4Momen
342 if(products->size()==2)
343 {
344 final=products->operator[](1);
345 pFinal +=final->Get4Momentum();
346 }
347
348 #ifdef debug_G4Scatterer
349 if ( (pInitial-pFinal).mag() > 0.1*MeV
350 {
351 G4cout << "G4Scatterer: momentum imb
352 }
353 G4cout << "Scatterer costh= " << trk1.G
354 #endif
355
356 for(size_t hpw=0; hpw<products->size();
357 {
358 energyBalance-=products->operator[](h
359 pxBalance-=products->operator[](hpw)-
360 pyBalance-=products->operator[](hpw)-
361 pzBalance-=products->operator[](hpw)-
362 chargeBalance-=G4lrint(products->operator[]
363 baryonBalance-=products->operator[](h
364 }
365 if(std::getenv("ScattererEnergyBalanceC
366 std::cout << "DEBUGGING energy balanc
367 <<energyBalance<<" "
368 <<pxBalance<<" "
369 <<pyBalance<<" "
370 <<pzBalance<<" "
371 <<chargeBalance<<" "
372 <<baryonBalance<<" "
373 <<G4endl;
374 if(chargeBalance !=0 )
375 {
376 G4cout << "track 1"<<trk1.GetDefiniti
377 G4cout << "track 2"<<trk2.GetDefiniti
378 for(size_t hpw=0; hpw<products->size(
379 {
380 G4cout << products->operator[](hpw
381 }
382 G4Exception("G4Scatterer", "im_r_matr
383 "Problem in ChargeBalance");
384 }
385 return products;
386 }
387 }
388
389 return NULL;
390 }
391
392 //--------------------------------------------
393
394 const G4VCollision* G4Scatterer::FindCollision
395 const G4KineticTrack& trk2) const
396 {
397 G4VCollision* collisionInCharge = 0;
398
399 size_t i;
400 for (i=0; i<collisions.size(); i++)
401 {
402 G4VCollision* component = collisions[i];
403 if (component->IsInCharge(trk1,trk2))
404 {
405 collisionInCharge = component;
406 break;
407 }
408 }
409 // if(collisionInCharge)
410 // {
411 // G4cout << "found collision : "
412 // << collisionInCharge->GetName()<< "
413 // << "for "
414 // << trk1.GetDefinition()->GetParticleName()
415 // << trk2.GetDefinition()->GetParticleName()
416 // << G4endl;;
417 // }
418 return collisionInCharge;
419 }
420
421 //--------------------------------------------
422
423 G4double G4Scatterer::GetCrossSection(const G4
424 const G4KineticTrack& trk2) cons
425 {
426 const G4VCollision* collision = FindCollisi
427 G4double aCrossSection = 0;
428 if (collision != 0)
429 {
430 aCrossSection = collision->CrossSection(t
431 }
432 return aCrossSection;
433 }
434
435 //--------------------------------------------
436
437 const std::vector<G4CollisionInitialState *> &
438 GetCollisions(G4KineticTrack * aProjectile,
439 std::vector<G4KineticTrack *> &
440 G4double aCurrentTime)
441 {
442 theCollisions.clear();
443 std::vector<G4KineticTrack *>::iterator j=so
444 for(; j != someCandidates.end(); ++j)
445 {
446 G4double collisionTime = GetTimeToInteract
447 if(collisionTime == DBL_MAX) // no collis
448 {
449 continue;
450 }
451 G4KineticTrackVector aTarget;
452 aTarget.push_back(*j);
453 theCollisions.push_back(
454 new G4CollisionInitialState(collisionTim
455 // G4cerr <<" !!!!!! debug collisions "<<
456 }
457 return theCollisions;
458 }
459
460
461 G4KineticTrackVector * G4Scatterer::
462 GetFinalState(G4KineticTrack * aProjectile,
463 std::vector<G4KineticTrack *> & theTar
464 {
465 G4KineticTrack target_reloc(*(theTargets[0
466 return Scatter(*aProjectile, target_reloc)
467 }
468 //--------------------------------------------
469