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

Last change on this file since 1007 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 14.8 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.13.2.4 2006/06/29 20:40:31 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 // global optimization
129// static const G4double maxCrossSection = 500*millibarn;
130 if(pi*distance>maxCrossSection) return time;
131
132 // charged particles special
133 static const G4double maxChargedCrossSection = 200*millibarn;
134 if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 &&
135 std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 &&
136 pi*distance>maxChargedCrossSection) return time;
137
138 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
139 // neutrons special
140 if(( trk1.GetDefinition() == G4Neutron::Neutron() ||
141 trk1.GetDefinition() == G4Neutron::Neutron() ) &&
142 sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time;
143
144/*
145 * if(distance <= sqr(1.14*fermi))
146 * {
147 * time = collisionTime;
148 *
149 * *
150 * * G4cout << "Scatter distance/time: " << std::sqrt(distance)/fermi <<
151 * * " / "<< time/ns << G4endl;
152 * * G4ThreeVector pos1=trk1.GetPosition();
153 * * G4ThreeVector pos2=trk2.GetPosition();
154 * * G4LorentzVector xmom1 = trk1.Get4Momentum();
155 * * G4LorentzVector xmom2 = trk2.Get4Momentum();
156 * * G4cout << "position1: " << pos1.x() << " " << pos1.y() << " "
157 * * << pos1.z();
158 * * pos1+=(collisionTime*c_light/xmom1.e())*xmom1.vect();
159 * * G4cout << " straight line trprt: "
160 * * << pos1.x() << " " << pos1.y() << " "
161 * * << pos1.z() << G4endl;
162 * * G4cout << "position2: " << pos2.x() << " " << pos2.y() << " "
163 * * << pos2.z() << G4endl;
164 * * G4cout << "straight line distance 2 fixed:" << (pos1-pos2).mag()/fermi << G4endl;
165 * * pos2+= (collisionTime*c_light/xmom2.e())*xmom2.vect();
166 * * G4cout<< " straight line trprt: "
167 * * << pos2.x() << " " << pos2.y() << " "
168 * * << pos2.z() << G4endl;
169 * * G4cout << "straight line distance :" << (pos1-pos2).mag()/fermi << G4endl;
170 * *
171 * }
172 *
173 * if(1)
174 * return time;
175 */
176
177 if ((trk1.GetActualMass()+trk2.GetActualMass()) > sqrtS) return time;
178
179
180
181 G4VCollision* collision = FindCollision(trk1,trk2);
182 G4double totalCrossSection;
183 // The cross section is interpreted geometrically as an area
184 // Two particles are assumed to collide if their distance is < (totalCrossSection/pi)
185
186 if (collision != 0)
187 {
188 totalCrossSection = collision->CrossSection(trk1,trk2);
189 if ( totalCrossSection > 0 )
190 {
191/* G4cout << " totalCrossection = "<< totalCrossSection << ", trk1/2, s, e-m: "
192 * << trk1.GetDefinition()->GetParticleName()
193 * << " / "
194 * << trk2.GetDefinition()->GetParticleName()
195 * << ", "
196 * << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
197 * << ", "
198 * << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
199 * trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
200 * << G4endl;
201 */
202 if (distance <= totalCrossSection / pi)
203 {
204 time = collisionTime;
205 }
206 } else
207 {
208
209 // For debugging...
210 // G4cout << " totalCrossection = 0, trk1/2, s, e-m: "
211 // << trk1.GetDefinition()->GetParticleName()
212 // << " / "
213 // << trk2.GetDefinition()->GetParticleName()
214 // << ", "
215 // << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()
216 // << ", "
217 // << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()-
218 // trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag()
219 // << G4endl;
220
221 }
222/*
223 * if(distance <= sqr(5.*fermi))
224 * {
225 * G4cout << " distance,xsect, std::sqrt(xsect/pi) : " << std::sqrt(distance)/fermi
226 * << " " << totalCrossSection/sqr(fermi)
227 * << " " << std::sqrt(totalCrossSection / pi)/fermi << G4endl;
228 * }
229 */
230
231 }
232 else
233 {
234 time = DBL_MAX;
235// /*
236 // For debugging
237//hpw G4cout << "G4Scatterer - collision not found: "
238//hpw << trk1.GetDefinition()->GetParticleName()
239//hpw << " - "
240//hpw << trk2.GetDefinition()->GetParticleName()
241//hpw << G4endl;
242 // End of debugging
243// */
244 }
245 }
246
247 else
248 {
249 /*
250 // For debugging
251 G4cout << "G4Scatterer - negative collisionTime"
252 << ": collisionTime = " << collisionTime
253 << ", position = " << position
254 << ", velocity = " << velocity
255 << G4endl;
256 // End of debugging
257 */
258 }
259
260 return time;
261}
262
263G4KineticTrackVector* G4Scatterer::Scatter(const G4KineticTrack& trk1,
264 const G4KineticTrack& trk2)
265{
266// G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag();
267 G4LorentzVector pInitial=trk1.Get4Momentum() + trk2.Get4Momentum();
268 G4double energyBalance = pInitial.t();
269 G4double pxBalance = pInitial.vect().x();
270 G4double pyBalance = pInitial.vect().y();
271 G4double pzBalance = pInitial.vect().z();
272 G4int chargeBalance = G4lrint(trk1.GetDefinition()->GetPDGCharge()
273 + trk2.GetDefinition()->GetPDGCharge());
274 G4int baryonBalance = trk1.GetDefinition()->GetBaryonNumber()
275 + trk2.GetDefinition()->GetBaryonNumber();
276
277 G4VCollision* collision = FindCollision(trk1,trk2);
278 if (collision != 0)
279 {
280 G4double aCrossSection = collision->CrossSection(trk1,trk2);
281 if (aCrossSection > 0.0)
282 {
283
284
285 #ifdef debug_G4Scatterer
286 G4cout << "be4 FinalState 1(p,e,m): "
287 << trk1.Get4Momentum() << " "
288 << trk1.Get4Momentum().mag()
289 << ", 2: "
290 << trk2.Get4Momentum()<< " "
291 << trk2.Get4Momentum().mag() << " "
292 << G4endl;
293 #endif
294
295
296 G4KineticTrackVector* products = collision->FinalState(trk1,trk2);
297 if(!products || products->size() == 0) return products;
298
299 #ifdef debug_G4Scatterer
300 G4cout << "size of FS: "<<products->size()<<G4endl;
301 #endif
302
303 G4KineticTrack *final= products->operator[](0);
304
305
306 #ifdef debug_G4Scatterer
307 G4cout << " FinalState 1: "
308 << final->Get4Momentum()<< " "
309 << final->Get4Momentum().mag() ;
310 #endif
311
312 if(products->size() == 1) return products;
313 final=products->operator[](1);
314 #ifdef debug_G4Scatterer
315 G4cout << ", 2: "
316 << final->Get4Momentum() << " "
317 << final->Get4Momentum().mag() << " " << G4endl;
318 #endif
319
320 final= products->operator[](0);
321 G4LorentzVector pFinal=final->Get4Momentum();
322 if(products->size()==2)
323 {
324 final=products->operator[](1);
325 pFinal +=final->Get4Momentum();
326 }
327
328 #ifdef debug_G4Scatterer
329 if ( (pInitial-pFinal).mag() > 0.1*MeV )
330 {
331 G4cout << "G4Scatterer: momentum imbalance, pInitial= " <<pInitial << " pFinal= " <<pFinal<< G4endl;
332 }
333 G4cout << "Scatterer costh= " << trk1.Get4Momentum().vect().unit() *(products->operator[](0))->Get4Momentum().vect().unit()<< G4endl;
334 #endif
335
336 for(size_t hpw=0; hpw<products->size(); hpw++)
337 {
338 energyBalance-=products->operator[](hpw)->Get4Momentum().t();
339 pxBalance-=products->operator[](hpw)->Get4Momentum().vect().x();
340 pyBalance-=products->operator[](hpw)->Get4Momentum().vect().y();
341 pzBalance-=products->operator[](hpw)->Get4Momentum().vect().z();
342 chargeBalance-=G4lrint(products->operator[](hpw)->GetDefinition()->GetPDGCharge());
343 baryonBalance-=products->operator[](hpw)->GetDefinition()->GetBaryonNumber();
344 }
345 if(getenv("ScattererEnergyBalanceCheck"))
346 std::cout << "DEBUGGING energy balance A: "
347 <<energyBalance<<" "
348 <<pxBalance<<" "
349 <<pyBalance<<" "
350 <<pzBalance<<" "
351 <<chargeBalance<<" "
352 <<baryonBalance<<" "
353 <<G4endl;
354 if(chargeBalance !=0 )
355 {
356 G4cout << "track 1"<<trk1.GetDefinition()->GetParticleName()<<G4endl;
357 G4cout << "track 2"<<trk2.GetDefinition()->GetParticleName()<<G4endl;
358 for(size_t hpw=0; hpw<products->size(); hpw++)
359 {
360 G4cout << products->operator[](hpw)->GetDefinition()->GetParticleName()<<G4endl;
361 }
362 G4Exception("We have the problem");
363 }
364 return products;
365 }
366 }
367
368 return NULL;
369}
370
371
372G4VCollision* G4Scatterer::FindCollision(const G4KineticTrack& trk1,
373 const G4KineticTrack& trk2)
374{
375 G4VCollision* collisionInCharge = 0;
376
377 size_t i;
378 for (i=0; i<collisions.size(); i++)
379 {
380 G4VCollision* component = collisions[i];
381 if (component->IsInCharge(trk1,trk2))
382 {
383 collisionInCharge = component;
384 break;
385 }
386 }
387// if(collisionInCharge)
388// {
389// G4cout << "found collision : "
390// << collisionInCharge->GetName()<< " "
391// << "for "
392// << trk1.GetDefinition()->GetParticleName()<<" + "
393// << trk2.GetDefinition()->GetParticleName()<<" "
394// << G4endl;;
395// }
396 return collisionInCharge;
397}
398
399G4double G4Scatterer::GetCrossSection(const G4KineticTrack& trk1,
400 const G4KineticTrack& trk2)
401{
402 G4VCollision* collision = FindCollision(trk1,trk2);
403 G4double aCrossSection = 0;
404 if (collision != 0)
405 {
406 aCrossSection = collision->CrossSection(trk1,trk2);
407 }
408 return aCrossSection;
409}
410
411
412const std::vector<G4CollisionInitialState *> & G4Scatterer::
413GetCollisions(G4KineticTrack * aProjectile,
414 std::vector<G4KineticTrack *> & someCandidates,
415 G4double aCurrentTime)
416{
417 theCollisions.clear();
418 std::vector<G4KineticTrack *>::iterator j=someCandidates.begin();
419 for(; j != someCandidates.end(); ++j)
420 {
421 G4double collisionTime = GetTimeToInteraction(*aProjectile, **j);
422 if(collisionTime == DBL_MAX) // no collision
423 {
424 continue;
425 }
426 G4KineticTrackVector aTarget;
427 aTarget.push_back(*j);
428 theCollisions.push_back(
429 new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) );
430// G4cerr <<" !!!!!! debug collisions "<<collisionTime<<" "<<pkt->GetDefinition()->GetParticleName()<<G4endl;
431 }
432 return theCollisions;
433}
434
435
436G4KineticTrackVector * G4Scatterer::
437GetFinalState(G4KineticTrack * aProjectile,
438 std::vector<G4KineticTrack *> & theTargets)
439{
440 G4KineticTrack target_reloc(*(theTargets[0]));
441 return Scatter(*aProjectile, target_reloc);
442}
Note: See TracBrowser for help on using the repository browser.