#include <G4Scatterer.hh>
Inheritance diagram for G4Scatterer:
Public Member Functions | |
G4Scatterer () | |
virtual | ~G4Scatterer () |
virtual G4double | GetTimeToInteraction (const G4KineticTrack &trk1, const G4KineticTrack &trk2) |
G4double | GetCrossSection (const G4KineticTrack &trk1, const G4KineticTrack &trk2) |
virtual G4KineticTrackVector * | Scatter (const G4KineticTrack &trk1, const G4KineticTrack &trk2) |
virtual const std::vector< G4CollisionInitialState * > & | GetCollisions (G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &someCandidates, G4double aCurrentTime) |
virtual G4KineticTrackVector * | GetFinalState (G4KineticTrack *aProjectile, std::vector< G4KineticTrack * > &theTargets) |
Data Structures | |
struct | Register |
Definition at line 44 of file G4Scatterer.hh.
G4Scatterer::G4Scatterer | ( | ) |
Definition at line 54 of file G4Scatterer.cc.
References G4ForEach< group >::Apply().
00055 { 00056 Register aR; 00057 G4ForEach<theChannels>::Apply(&aR, &collisions); 00058 }
G4Scatterer::~G4Scatterer | ( | ) | [virtual] |
Definition at line 62 of file G4Scatterer.cc.
00063 { 00064 std::for_each(collisions.begin(), collisions.end(), G4Delete()); 00065 collisions.clear(); 00066 }
const std::vector< G4CollisionInitialState * > & G4Scatterer::GetCollisions | ( | G4KineticTrack * | aProjectile, | |
std::vector< G4KineticTrack * > & | someCandidates, | |||
G4double | aCurrentTime | |||
) | [virtual] |
Implements G4BCAction.
Definition at line 424 of file G4Scatterer.cc.
References DBL_MAX, and GetTimeToInteraction().
00427 { 00428 theCollisions.clear(); 00429 std::vector<G4KineticTrack *>::iterator j=someCandidates.begin(); 00430 for(; j != someCandidates.end(); ++j) 00431 { 00432 G4double collisionTime = GetTimeToInteraction(*aProjectile, **j); 00433 if(collisionTime == DBL_MAX) // no collision 00434 { 00435 continue; 00436 } 00437 G4KineticTrackVector aTarget; 00438 aTarget.push_back(*j); 00439 theCollisions.push_back( 00440 new G4CollisionInitialState(collisionTime+aCurrentTime, aProjectile, aTarget, this) ); 00441 // G4cerr <<" !!!!!! debug collisions "<<collisionTime<<" "<<pkt->GetDefinition()->GetParticleName()<<G4endl; 00442 } 00443 return theCollisions; 00444 }
G4double G4Scatterer::GetCrossSection | ( | const G4KineticTrack & | trk1, | |
const G4KineticTrack & | trk2 | |||
) |
Definition at line 409 of file G4Scatterer.cc.
References G4VCollision::CrossSection().
00411 { 00412 G4VCollision* collision = FindCollision(trk1,trk2); 00413 G4double aCrossSection = 0; 00414 if (collision != 0) 00415 { 00416 aCrossSection = collision->CrossSection(trk1,trk2); 00417 } 00418 return aCrossSection; 00419 }
G4KineticTrackVector * G4Scatterer::GetFinalState | ( | G4KineticTrack * | aProjectile, | |
std::vector< G4KineticTrack * > & | theTargets | |||
) | [virtual] |
Implements G4BCAction.
Definition at line 448 of file G4Scatterer.cc.
References Scatter().
00450 { 00451 G4KineticTrack target_reloc(*(theTargets[0])); 00452 return Scatter(*aProjectile, target_reloc); 00453 }
G4double G4Scatterer::GetTimeToInteraction | ( | const G4KineticTrack & | trk1, | |
const G4KineticTrack & | trk2 | |||
) | [virtual] |
Implements G4VScatterer.
Definition at line 70 of file G4Scatterer.cc.
References G4VCollision::CrossSection(), DBL_MAX, G4KineticTrack::Get4Momentum(), G4KineticTrack::GetActualMass(), G4KineticTrack::GetDefinition(), G4ParticleDefinition::GetPDGCharge(), G4KineticTrack::GetPosition(), G4KineticTrack::GetTrackingMomentum(), G4Neutron::Neutron(), G4INCL::Math::pi, and position.
Referenced by GetCollisions().
00072 { 00073 G4double time = DBL_MAX; 00074 G4double distance_fast; 00075 G4LorentzVector mom1 = trk1.GetTrackingMomentum(); 00076 // G4cout << "zcomp=" << std::abs(mom1.vect().unit().z() -1 ) << G4endl; 00077 G4double collisionTime; 00078 00079 if ( std::abs(mom1.vect().unit().z() -1 ) < 1e-6 ) 00080 { 00081 G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition(); 00082 G4double deltaz=position.z(); 00083 G4double velocity = mom1.z()/mom1.e() * c_light; 00084 00085 collisionTime=deltaz/velocity; 00086 distance_fast=position.x()*position.x() + position.y()*position.y(); 00087 } else { 00088 00089 // The nucleons of the nucleus are FROZEN, ie. do not move.. 00090 00091 G4ThreeVector position = trk2.GetPosition() - trk1.GetPosition(); 00092 00093 G4ThreeVector velocity = mom1.vect()/mom1.e() * c_light; // mom1.boostVector() will exit on slightly negative mass 00094 collisionTime = (position * velocity) / velocity.mag2(); // can't divide by /c_light; 00095 position -= velocity * collisionTime; 00096 distance_fast=position.mag2(); 00097 00098 // if ( collisionTime>0 ) G4cout << " dis1/2 square" << dis1 <<" "<< dis2 << G4endl; 00099 // collisionTime = GetTimeToClosestApproach(trk1,trk2); 00100 } 00101 if (collisionTime > 0) 00102 { 00103 static const G4double maxCrossSection = 500*millibarn; 00104 if(0.7*pi*distance_fast>maxCrossSection) return time; 00105 00106 00107 G4LorentzVector mom2(0,0,0,trk2.Get4Momentum().mag()); 00108 00109 // G4ThreeVector momLab = mom1.vect();// frozen Nucleus - mom2.vect(); 00110 // G4ThreeVector posLab = trk1.GetPosition() - trk2.GetPosition(); 00111 // G4double disLab=posLab * posLab - (posLab*momLab) * (posLab*momLab) /(momLab.mag2()); 00112 00113 G4LorentzRotation toCMSFrame((-1)*(mom1 + mom2).boostVector()); 00114 mom1 = toCMSFrame * mom1; 00115 mom2 = toCMSFrame * mom2; 00116 00117 G4LorentzVector coordinate1(trk1.GetPosition(), 100.); 00118 G4LorentzVector coordinate2(trk2.GetPosition(), 100.); 00119 G4ThreeVector pos = ((toCMSFrame * coordinate1).vect() - 00120 (toCMSFrame * coordinate2).vect()); 00121 00122 G4ThreeVector mom = mom1.vect() - mom2.vect(); 00123 00124 // Calculate the impact parameter 00125 00126 G4double distance = pos * pos - (pos*mom) * (pos*mom) / (mom.mag2()); 00127 00128 // G4cout << " disDiff " << distance-disLab << " " << disLab 00129 // << " " << std::abs(distance-disLab)/distance << G4endl 00130 // << " mom/Lab " << mom << " " << momLab << G4endl 00131 // << " pos/Lab " << pos << " " << posLab 00132 // << G4endl; 00133 00134 if(pi*distance>maxCrossSection) return time; 00135 00136 // charged particles special 00137 static const G4double maxChargedCrossSection = 200*millibarn; 00138 if(std::abs(trk1.GetDefinition()->GetPDGCharge())>0.1 && 00139 std::abs(trk2.GetDefinition()->GetPDGCharge())>0.1 && 00140 pi*distance>maxChargedCrossSection) return time; 00141 00142 G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag(); 00143 // neutrons special pn is largest cross-section, but above 1.91 GeV is less than 200 mb 00144 if(( trk1.GetDefinition() == G4Neutron::Neutron() || 00145 trk2.GetDefinition() == G4Neutron::Neutron() ) && 00146 sqrtS>1.91*GeV && pi*distance>maxChargedCrossSection) return time; 00147 00148 /* 00149 * if(distance <= sqr(1.14*fermi)) 00150 * { 00151 * time = collisionTime; 00152 * 00153 * * 00154 * * G4cout << "Scatter distance/time: " << std::sqrt(distance)/fermi << 00155 * * " / "<< time/ns << G4endl; 00156 * * G4ThreeVector pos1=trk1.GetPosition(); 00157 * * G4ThreeVector pos2=trk2.GetPosition(); 00158 * * G4LorentzVector xmom1 = trk1.Get4Momentum(); 00159 * * G4LorentzVector xmom2 = trk2.Get4Momentum(); 00160 * * G4cout << "position1: " << pos1.x() << " " << pos1.y() << " " 00161 * * << pos1.z(); 00162 * * pos1+=(collisionTime*c_light/xmom1.e())*xmom1.vect(); 00163 * * G4cout << " straight line trprt: " 00164 * * << pos1.x() << " " << pos1.y() << " " 00165 * * << pos1.z() << G4endl; 00166 * * G4cout << "position2: " << pos2.x() << " " << pos2.y() << " " 00167 * * << pos2.z() << G4endl; 00168 * * G4cout << "straight line distance 2 fixed:" << (pos1-pos2).mag()/fermi << G4endl; 00169 * * pos2+= (collisionTime*c_light/xmom2.e())*xmom2.vect(); 00170 * * G4cout<< " straight line trprt: " 00171 * * << pos2.x() << " " << pos2.y() << " " 00172 * * << pos2.z() << G4endl; 00173 * * G4cout << "straight line distance :" << (pos1-pos2).mag()/fermi << G4endl; 00174 * * 00175 * } 00176 * 00177 * if(1) 00178 * return time; 00179 */ 00180 00181 if ((trk1.GetActualMass()+trk2.GetActualMass()) > sqrtS) return time; 00182 00183 00184 00185 G4VCollision* collision = FindCollision(trk1,trk2); 00186 G4double totalCrossSection; 00187 // The cross section is interpreted geometrically as an area 00188 // Two particles are assumed to collide if their distance is < (totalCrossSection/pi) 00189 00190 if (collision != 0) 00191 { 00192 totalCrossSection = collision->CrossSection(trk1,trk2); 00193 if ( totalCrossSection > 0 ) 00194 { 00195 /* G4cout << " totalCrossection = "<< totalCrossSection << ", trk1/2, s, e-m: " 00196 * << trk1.GetDefinition()->GetParticleName() 00197 * << " / " 00198 * << trk2.GetDefinition()->GetParticleName() 00199 * << ", " 00200 * << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag() 00201 * << ", " 00202 * << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()- 00203 * trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag() 00204 * << G4endl; 00205 */ 00206 if (distance <= totalCrossSection / pi) 00207 { 00208 time = collisionTime; 00209 } 00210 } else 00211 { 00212 00213 // For debugging... 00214 // G4cout << " totalCrossection = 0, trk1/2, s, e-m: " 00215 // << trk1.GetDefinition()->GetParticleName() 00216 // << " / " 00217 // << trk2.GetDefinition()->GetParticleName() 00218 // << ", " 00219 // << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag() 00220 // << ", " 00221 // << (trk1.Get4Momentum()+trk2.Get4Momentum()).mag()- 00222 // trk1.Get4Momentum().mag() - trk2.Get4Momentum().mag() 00223 // << G4endl; 00224 00225 } 00226 /* 00227 * if(distance <= sqr(5.*fermi)) 00228 * { 00229 * G4cout << " distance,xsect, std::sqrt(xsect/pi) : " << std::sqrt(distance)/fermi 00230 * << " " << totalCrossSection/sqr(fermi) 00231 * << " " << std::sqrt(totalCrossSection / pi)/fermi << G4endl; 00232 * } 00233 */ 00234 00235 } 00236 else 00237 { 00238 time = DBL_MAX; 00239 // /* 00240 // For debugging 00241 //hpw G4cout << "G4Scatterer - collision not found: " 00242 //hpw << trk1.GetDefinition()->GetParticleName() 00243 //hpw << " - " 00244 //hpw << trk2.GetDefinition()->GetParticleName() 00245 //hpw << G4endl; 00246 // End of debugging 00247 // */ 00248 } 00249 } 00250 00251 else 00252 { 00253 /* 00254 // For debugging 00255 G4cout << "G4Scatterer - negative collisionTime" 00256 << ": collisionTime = " << collisionTime 00257 << ", position = " << position 00258 << ", velocity = " << velocity 00259 << G4endl; 00260 // End of debugging 00261 */ 00262 } 00263 00264 return time; 00265 }
G4KineticTrackVector * G4Scatterer::Scatter | ( | const G4KineticTrack & | trk1, | |
const G4KineticTrack & | trk2 | |||
) | [virtual] |
Implements G4VScatterer.
Definition at line 269 of file G4Scatterer.cc.
References G4VCollision::CrossSection(), FatalException, G4VCollision::FinalState(), G4cout, G4endl, G4Exception(), G4lrint(), G4KineticTrack::Get4Momentum(), G4ParticleDefinition::GetBaryonNumber(), G4KineticTrack::GetDefinition(), G4ParticleDefinition::GetParticleName(), and G4ParticleDefinition::GetPDGCharge().
Referenced by G4QMDCollision::CalFinalStateOfTheBinaryCollision(), and GetFinalState().
00271 { 00272 // G4double sqrtS = (trk1.Get4Momentum() + trk2.Get4Momentum()).mag(); 00273 G4LorentzVector pInitial=trk1.Get4Momentum() + trk2.Get4Momentum(); 00274 G4double energyBalance = pInitial.t(); 00275 G4double pxBalance = pInitial.vect().x(); 00276 G4double pyBalance = pInitial.vect().y(); 00277 G4double pzBalance = pInitial.vect().z(); 00278 G4int chargeBalance = G4lrint(trk1.GetDefinition()->GetPDGCharge() 00279 + trk2.GetDefinition()->GetPDGCharge()); 00280 G4int baryonBalance = trk1.GetDefinition()->GetBaryonNumber() 00281 + trk2.GetDefinition()->GetBaryonNumber(); 00282 00283 G4VCollision* collision = FindCollision(trk1,trk2); 00284 if (collision != 0) 00285 { 00286 G4double aCrossSection = collision->CrossSection(trk1,trk2); 00287 if (aCrossSection > 0.0) 00288 { 00289 00290 00291 #ifdef debug_G4Scatterer 00292 G4cout << "be4 FinalState 1(p,e,m): " 00293 << trk1.Get4Momentum() << " " 00294 << trk1.Get4Momentum().mag() 00295 << ", 2: " 00296 << trk2.Get4Momentum()<< " " 00297 << trk2.Get4Momentum().mag() << " " 00298 << G4endl; 00299 #endif 00300 00301 00302 G4KineticTrackVector* products = collision->FinalState(trk1,trk2); 00303 if(!products || products->size() == 0) return products; 00304 00305 #ifdef debug_G4Scatterer 00306 G4cout << "size of FS: "<<products->size()<<G4endl; 00307 #endif 00308 00309 G4KineticTrack *final= products->operator[](0); 00310 00311 00312 #ifdef debug_G4Scatterer 00313 G4cout << " FinalState 1: " 00314 << final->Get4Momentum()<< " " 00315 << final->Get4Momentum().mag() ; 00316 #endif 00317 00318 if(products->size() == 1) return products; 00319 final=products->operator[](1); 00320 #ifdef debug_G4Scatterer 00321 G4cout << ", 2: " 00322 << final->Get4Momentum() << " " 00323 << final->Get4Momentum().mag() << " " << G4endl; 00324 #endif 00325 00326 final= products->operator[](0); 00327 G4LorentzVector pFinal=final->Get4Momentum(); 00328 if(products->size()==2) 00329 { 00330 final=products->operator[](1); 00331 pFinal +=final->Get4Momentum(); 00332 } 00333 00334 #ifdef debug_G4Scatterer 00335 if ( (pInitial-pFinal).mag() > 0.1*MeV ) 00336 { 00337 G4cout << "G4Scatterer: momentum imbalance, pInitial= " <<pInitial << " pFinal= " <<pFinal<< G4endl; 00338 } 00339 G4cout << "Scatterer costh= " << trk1.Get4Momentum().vect().unit() *(products->operator[](0))->Get4Momentum().vect().unit()<< G4endl; 00340 #endif 00341 00342 for(size_t hpw=0; hpw<products->size(); hpw++) 00343 { 00344 energyBalance-=products->operator[](hpw)->Get4Momentum().t(); 00345 pxBalance-=products->operator[](hpw)->Get4Momentum().vect().x(); 00346 pyBalance-=products->operator[](hpw)->Get4Momentum().vect().y(); 00347 pzBalance-=products->operator[](hpw)->Get4Momentum().vect().z(); 00348 chargeBalance-=G4lrint(products->operator[](hpw)->GetDefinition()->GetPDGCharge()); 00349 baryonBalance-=products->operator[](hpw)->GetDefinition()->GetBaryonNumber(); 00350 } 00351 if(getenv("ScattererEnergyBalanceCheck")) 00352 std::cout << "DEBUGGING energy balance A: " 00353 <<energyBalance<<" " 00354 <<pxBalance<<" " 00355 <<pyBalance<<" " 00356 <<pzBalance<<" " 00357 <<chargeBalance<<" " 00358 <<baryonBalance<<" " 00359 <<G4endl; 00360 if(chargeBalance !=0 ) 00361 { 00362 G4cout << "track 1"<<trk1.GetDefinition()->GetParticleName()<<G4endl; 00363 G4cout << "track 2"<<trk2.GetDefinition()->GetParticleName()<<G4endl; 00364 for(size_t hpw=0; hpw<products->size(); hpw++) 00365 { 00366 G4cout << products->operator[](hpw)->GetDefinition()->GetParticleName()<<G4endl; 00367 } 00368 G4Exception("G4Scatterer", "im_r_matrix001", FatalException, 00369 "Problem in ChargeBalance"); 00370 } 00371 return products; 00372 } 00373 } 00374 00375 return NULL; 00376 }