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