source: trunk/source/processes/hadronic/models/im_r_matrix/src/G4Scatterer.cc@ 1330

Last change on this file since 1330 was 1315, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 14.7 KB
Line 
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// $Id: G4Scatterer.cc,v 1.16 2010/03/12 15:45:18 gunter Exp $ //
27//
28
29#include "globals.hh"
30#include <vector>
31#include "G4ios.hh"
32#include "G4Scatterer.hh"
33#include "G4KineticTrack.hh"
34#include "G4ThreeVector.hh"
35#include "G4LorentzRotation.hh"
36#include "G4LorentzVector.hh"
37
38#include "G4CollisionNN.hh"
39#include "G4CollisionPN.hh"
40#include "G4CollisionMesonBaryon.hh"
41
42#include "G4CollisionInitialState.hh"
43#include "G4HadTmpUtil.hh"
44#include "G4Pair.hh"
45
46
47// Declare the categories of collisions the Scatterer can handle
48typedef GROUP2(G4CollisionNN, G4CollisionMesonBaryon) theChannels;
49
50
51G4Scatterer::G4Scatterer()
52{
53 Register aR;
54 G4ForEach<theChannels>::Apply(&aR, &collisions);
55}
56
57
58G4Scatterer::~G4Scatterer()
59{
60 std::for_each(collisions.begin(), collisions.end(), G4Delete());
61 collisions.clear();
62}
63
64G4double G4Scatterer::GetTimeToInteraction(const G4KineticTrack& trk1,
65 const G4KineticTrack& trk2)
66{
67 G4double time = DBL_MAX;
68 G4double distance_fast;
69 G4LorentzVector mom1 = trk1.GetTrackingMomentum();
70// G4cout << "zcomp=" << std::abs(mom1.vect().unit().z() -1 ) << G4endl;
71 G4double collisionTime;
72
73 if ( std::abs(mom1.vect().unit().z() -1 ) < 1e-6 )
74 {
75 G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
76 G4double deltaz=position.z();
77 G4double velocity = mom1.z()/mom1.e() * c_light;
78
79 collisionTime=deltaz/velocity;
80 distance_fast=position.x()*position.x() + position.y()*position.y();
81 } else {
82
83 // The nucleons of the nucleus are FROZEN, ie. do not move..
84
85 G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition();
86
87 G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light; // mom1.boostVector() will exit on slightly negative mass
88 collisionTime = (position * velocity) / velocity.mag2(); // can't divide by /c_light;
89 position -= velocity * collisionTime;
90 distance_fast=position.mag2();
91
92// if ( collisionTime>0 ) G4cout << " dis1/2 square" << dis1 <<" "<< dis2 << G4endl;
93// collisionTime = GetTimeToClosestApproach(trk1,trk2);
94 }
95 if (collisionTime > 0)
96 {
97 static const G4double maxCrossSection = 500*millibarn;
98 if(0.7*pi*distance_fast>maxCrossSection) return time;
99
100
101 G4LorentzVector mom2(0,0,0,trk2.Get4Momentum().mag());
102
103// G4ThreeVector momLab = mom1.vect();// frozen Nucleus - mom2.vect();
104// G4ThreeVector posLab = trk1.GetPosition() - trk2.GetPosition();
105// G4double disLab=posLab * posLab - (posLab*momLab) * (posLab*momLab) /(momLab.mag2());
106
107 G4LorentzRotation toCMSFrame((-1)*(mom1 + mom2).boostVector());
108 mom1 = toCMSFrame * mom1;
109 mom2 = toCMSFrame * mom2;
110
111 G4LorentzVector coordinate1(trk1.GetPosition(), 100.);
112 G4LorentzVector coordinate2(trk2.GetPosition(), 100.);
113 G4ThreeVector pos = ((toCMSFrame * coordinate1).vect() -
114 (toCMSFrame * coordinate2).vect());
115
116 G4ThreeVector mom = mom1.vect() - mom2.vect();
117
118 // Calculate the impact parameter
119
120 G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom.mag2());
121
122// G4cout << " disDiff " << distance-disLab << " " << disLab
123// << " " << std::abs(distance-disLab)/distance << G4endl
124// << " mom/Lab " << mom << " " << momLab << G4endl
125// << " pos/Lab " << pos << " " << posLab
126// << G4endl;
127
128 if(pi*distance>maxCrossSection) return time;
129
130 // charged particles special
131 static const G4double maxChargedCrossSection = 200*millibarn;
132 if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 &&
133 std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 &&
134 pi*distance>maxChargedCrossSection) return time;
135
136 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
137 // neutrons special
138 if(( trk1.GetDefinition() == G4Neutron::Neutron() ||
139 trk1.GetDefinition() == G4Neutron::Neutron() ) &&
140 sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time;
141
142/*
143 * if(distance <= sqr(1.14*fermi))
144 * {
145 * time = collisionTime;
146 *
147 * *
148 * * G4cout << "Scatter distance/time: " << std::sqrt(distance)/fermi <<
149 * * " / "<< time/ns << G4endl;
150 * * G4ThreeVector pos1=trk1.GetPosition();
151 * * G4ThreeVector pos2=trk2.GetPosition();
152 * * G4LorentzVector xmom1 = trk1.Get4Momentum();
153 * * G4LorentzVector xmom2 = trk2.Get4Momentum();
154 * * G4cout << "position1: " << pos1.x() << " " << pos1.y() << " "
155 * * << pos1.z();
156 * * pos1+=(collisionTime*c_light/xmom1.e())*xmom1.vect();
157 * * G4cout << " straight line trprt: "
158 * * << pos1.x() << " " << pos1.y() << " "
159 * * << pos1.z() << G4endl;
160 * * G4cout << "position2: " << pos2.x() << " " << pos2.y() << " "
161 * * << pos2.z() << G4endl;
162 * * G4cout << "straight line distance 2 fixed:" << (pos1-pos2).mag()/fermi << G4endl;
163 * * pos2+= (collisionTime*c_light/xmom2.e())*xmom2.vect();
164 * * G4cout<< " straight line trprt: "
165 * * << pos2.x() << " " << pos2.y() << " "
166 * * << pos2.z() << G4endl;
167 * * G4cout << "straight line distance :" << (pos1-pos2).mag()/fermi << G4endl;
168 * *
169 * }
170 *
171 * if(1)
172 * return time;
173 */
174
175 if ((trk1.GetActualMass()+trk2.GetActualMass()) > sqrtS) return time;
176
177
178
179 G4VCollision* collision = FindCollision(trk1,trk2);
180 G4double totalCrossSection;
181 // The cross section is interpreted geometrically as an area
182 // Two particles are assumed to collide if their distance is < (totalCrossSection/pi)
183
184 if (collision != 0)
185 {
186 totalCrossSection = collision->CrossSection(trk1,trk2);
187 if ( totalCrossSection > 0 )
188 {
189/* G4cout << " totalCrossection = "<< totalCrossSection << ", trk1/2, s, e-m: "
190 * << trk1.GetDefinition()->GetParticleName()
191 * << " / "
192 * << trk2.GetDefinition()->GetParticleName()
193 * << ", "
194 * << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
195 * << ", "
196 * << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
197 * trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
198 * << G4endl;
199 */
200 if (distance <= totalCrossSection / pi)
201 {
202 time = collisionTime;
203 }
204 } else
205 {
206
207 // For debugging...
208 // G4cout << " totalCrossection = 0, trk1/2, s, e-m: "
209 // << trk1.GetDefinition()->GetParticleName()
210 // << " / "
211 // << trk2.GetDefinition()->GetParticleName()
212 // << ", "
213 // << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
214 // << ", "
215 // << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
216 // trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
217 // << G4endl;
218
219 }
220/*
221 * if(distance <= sqr(5.*fermi))
222 * {
223 * G4cout << " distance,xsect, std::sqrt(xsect/pi) : " << std::sqrt(distance)/fermi
224 * << " " << totalCrossSection/sqr(fermi)
225 * << " " << std::sqrt(totalCrossSection / pi)/fermi << G4endl;
226 * }
227 */
228
229 }
230 else
231 {
232 time = DBL_MAX;
233// /*
234 // For debugging
235//hpw G4cout << "G4Scatterer - collision not found: "
236//hpw << trk1.GetDefinition()->GetParticleName()
237//hpw << " - "
238//hpw << trk2.GetDefinition()->GetParticleName()
239//hpw << G4endl;
240 // End of debugging
241// */
242 }
243 }
244
245 else
246 {
247 /*
248 // For debugging
249 G4cout << "G4Scatterer - negative collisionTime"
250 << ": collisionTime = " << collisionTime
251 << ", position = " << position
252 << ", velocity = " << velocity
253 << G4endl;
254 // End of debugging
255 */
256 }
257
258 return time;
259}
260
261G4KineticTrackVector* G4Scatterer::Scatter(const G4KineticTrack& trk1,
262 const G4KineticTrack& trk2)
263{
264// G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
265 G4LorentzVector pInitial=trk1.Get4Momentum() + trk2.Get4Momentum();
266 G4double energyBalance = pInitial.t();
267 G4double pxBalance = pInitial.vect().x();
268 G4double pyBalance = pInitial.vect().y();
269 G4double pzBalance = pInitial.vect().z();
270 G4int chargeBalance = G4lrint(trk1.GetDefinition()->GetPDGCharge()
271 + trk2.GetDefinition()->GetPDGCharge());
272 G4int baryonBalance = trk1.GetDefinition()->GetBaryonNumber()
273 + trk2.GetDefinition()->GetBaryonNumber();
274
275 G4VCollision* collision = FindCollision(trk1,trk2);
276 if (collision != 0)
277 {
278 G4double aCrossSection = collision->CrossSection(trk1,trk2);
279 if (aCrossSection > 0.0)
280 {
281
282
283 #ifdef debug_G4Scatterer
284 G4cout << "be4 FinalState 1(p,e,m): "
285 << trk1.Get4Momentum() << " "
286 << trk1.Get4Momentum().mag()
287 << ", 2: "
288 << trk2.Get4Momentum()<< " "
289 << trk2.Get4Momentum().mag() << " "
290 << G4endl;
291 #endif
292
293
294 G4KineticTrackVector* products = collision->FinalState(trk1,trk2);
295 if(!products || products->size() == 0) return products;
296
297 #ifdef debug_G4Scatterer
298 G4cout << "size of FS: "<<products->size()<<G4endl;
299 #endif
300
301 G4KineticTrack *final= products->operator[](0);
302
303
304 #ifdef debug_G4Scatterer
305 G4cout << " FinalState 1: "
306 << final->Get4Momentum()<< " "
307 << final->Get4Momentum().mag() ;
308 #endif
309
310 if(products->size() == 1) return products;
311 final=products->operator[](1);
312 #ifdef debug_G4Scatterer
313 G4cout << ", 2: "
314 << final->Get4Momentum() << " "
315 << final->Get4Momentum().mag() << " " << G4endl;
316 #endif
317
318 final= products->operator[](0);
319 G4LorentzVector pFinal=final->Get4Momentum();
320 if(products->size()==2)
321 {
322 final=products->operator[](1);
323 pFinal +=final->Get4Momentum();
324 }
325
326 #ifdef debug_G4Scatterer
327 if ( (pInitial-pFinal).mag() > 0.1*MeV )
328 {
329 G4cout << "G4Scatterer: momentum imbalance, pInitial= " <<pInitial << " pFinal= " <<pFinal<< G4endl;
330 }
331 G4cout << "Scatterer costh= " << trk1.Get4Momentum().vect().unit() *(products->operator[](0))->Get4Momentum().vect().unit()<< G4endl;
332 #endif
333
334 for(size_t hpw=0; hpw<products->size(); hpw++)
335 {
336 energyBalance-=products->operator[](hpw)->Get4Momentum().t();
337 pxBalance-=products->operator[](hpw)->Get4Momentum().vect().x();
338 pyBalance-=products->operator[](hpw)->Get4Momentum().vect().y();
339 pzBalance-=products->operator[](hpw)->Get4Momentum().vect().z();
340 chargeBalance-=G4lrint(products->operator[](hpw)->GetDefinition()->GetPDGCharge());
341 baryonBalance-=products->operator[](hpw)->GetDefinition()->GetBaryonNumber();
342 }
343 if(getenv("ScattererEnergyBalanceCheck"))
344 std::cout << "DEBUGGING energy balance A: "
345 <<energyBalance<<" "
346 <<pxBalance<<" "
347 <<pyBalance<<" "
348 <<pzBalance<<" "
349 <<chargeBalance<<" "
350 <<baryonBalance<<" "
351 <<G4endl;
352 if(chargeBalance !=0 )
353 {
354 G4cout << "track 1"<<trk1.GetDefinition()->GetParticleName()<<G4endl;
355 G4cout << "track 2"<<trk2.GetDefinition()->GetParticleName()<<G4endl;
356 for(size_t hpw=0; hpw<products->size(); hpw++)
357 {
358 G4cout << products->operator[](hpw)->GetDefinition()->GetParticleName()<<G4endl;
359 }
360 G4Exception("We have the problem");
361 }
362 return products;
363 }
364 }
365
366 return NULL;
367}
368
369
370G4VCollision* G4Scatterer::FindCollision(const G4KineticTrack& trk1,
371 const G4KineticTrack& trk2)
372{
373 G4VCollision* collisionInCharge = 0;
374
375 size_t i;
376 for (i=0; i<collisions.size(); i++)
377 {
378 G4VCollision* component = collisions[i];
379 if (component->IsInCharge(trk1,trk2))
380 {
381 collisionInCharge = component;
382 break;
383 }
384 }
385// if(collisionInCharge)
386// {
387// G4cout << "found collision : "
388// << collisionInCharge->GetName()<< " "
389// << "for "
390// << trk1.GetDefinition()->GetParticleName()<<" + "
391// << trk2.GetDefinition()->GetParticleName()<<" "
392// << G4endl;;
393// }
394 return collisionInCharge;
395}
396
397G4double G4Scatterer::GetCrossSection(const G4KineticTrack& trk1,
398 const G4KineticTrack& trk2)
399{
400 G4VCollision* collision = FindCollision(trk1,trk2);
401 G4double aCrossSection = 0;
402 if (collision != 0)
403 {
404 aCrossSection = collision->CrossSection(trk1,trk2);
405 }
406 return aCrossSection;
407}
408
409
410const std::vector<G4CollisionInitialState *> & G4Scatterer::
411GetCollisions(G4KineticTrack * aProjectile,
412 std::vector<G4KineticTrack *> & someCandidates,
413 G4double aCurrentTime)
414{
415 theCollisions.clear();
416 std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
417 for(; j != someCandidates.end(); ++j)
418 {
419 G4double collisionTime = GetTimeToInteraction(*aProjectile, **j);
420 if(collisionTime == DBL_MAX) // no collision
421 {
422 continue;
423 }
424 G4KineticTrackVector aTarget;
425 aTarget.push_back(*j);
426 theCollisions.push_back(
427 new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) );
428// G4cerr <<" !!!!!! debug collisions "<<collisionTime<<" "<<pkt->GetDefinition()->GetParticleName()<<G4endl;
429 }
430 return theCollisions;
431}
432
433
434G4KineticTrackVector * G4Scatterer::
435GetFinalState(G4KineticTrack * aProjectile,
436 std::vector<G4KineticTrack *> & theTargets)
437{
438 G4KineticTrack target_reloc(*(theTargets[0]));
439 return Scatter(*aProjectile, target_reloc);
440}
Note: See TracBrowser for help on using the repository browser.