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