Geant4 Cross Reference |
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 // 080602 Fix memory leaks by T. Koi 27 // 081120 Add deltaT in signature of CalKinematicsOfBinaryCollisions 28 // Add several required updating of Mean Filed 29 // Modified handling of absorption case by T. Koi 30 // 090126 Fix in absorption case by T. Koi 31 // 090331 Fix for gamma participant by T. Koi 32 // 33 // 230308 Tentative modified in a short-lived particle production by Y-H. Sato and A. Haga 34 // 230308 Energy difference evaluated by "GetTotalEnergy" calculated in Mean Field (Y-H. Sato and A. Haga) 35 // 36 #include "G4LightIonQMDCollision.hh" 37 #include "G4Scatterer.hh" 38 #include "G4Pow.hh" 39 #include "G4Exp.hh" 40 #include "G4Log.hh" 41 #include "G4PhysicalConstants.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "Randomize.hh" 44 45 G4LightIonQMDCollision::G4LightIonQMDCollision() 46 : fdeltar ( 4.0 ) 47 , fbcmax0 ( 1.323142 ) // NN maximum impact parameter 48 , fbcmax1 ( 2.523 ) // others maximum impact parameter 49 // , sig0 ( 55 ) // NN cross section 50 //110617 fix for gcc 4.6 compilation warnings 51 //, sig1 ( 200 ) // others cross section 52 , fepse ( 0.0001 ) 53 { 54 //These two pointers will be set through SetMeanField method 55 theSystem=NULL; 56 theMeanField=NULL; 57 theScatterer = new G4Scatterer(); 58 } 59 60 /* 61 G4LightIonQMDCollision::G4LightIonQMDCollision( const G4LightIonQMDCollision& obj ) 62 : fdeltar ( obj.fdeltar ) 63 , fbcmax0 ( obj.fbcmax0 ) // NN maximum impact parameter 64 , fbcmax1 ( obj.fbcmax1 ) // others maximum impact parameter 65 , fepse ( obj.fepse ) 66 { 67 68 if ( obj.theSystem != NULL ) { 69 theSystem = new G4QMDSystem; 70 *theSystem = *obj.theSystem; 71 } else { 72 theSystem = NULL; 73 } 74 if ( obj.theMeanField != NULL ) { 75 theMeanField = new G4LightIonQMDMeanField; 76 *theMeanField = *obj.theMeanField; 77 } else { 78 theMeanField = NULL; 79 } 80 theScatterer = new G4Scatterer(); 81 *theScatterer = *obj.theScatterer; 82 } 83 84 G4LightIonQMDCollision & G4LightIonQMDCollision::operator= ( const G4LightIonQMDCollision& obj) 85 { 86 fdeltar = obj.fdeltar; 87 fbcmax0 = obj.fbcmax1; 88 fepse = obj.fepse; 89 90 if ( obj.theSystem != NULL ) { 91 delete theSystem; 92 theSystem = new G4QMDSystem; 93 *theSystem = *obj.theSystem; 94 } else { 95 theSystem = NULL; 96 } 97 if ( obj.theMeanField != NULL ) { 98 delete theMeanField; 99 theMeanField = new G4LightIonQMDMeanField; 100 *theMeanField = *obj.theMeanField; 101 } else { 102 theMeanField = NULL; 103 } 104 delete theScatterer; 105 theScatterer = new G4Scatterer(); 106 *theScatterer = *obj.theScatterer; 107 108 return *this; 109 } 110 */ 111 112 113 G4LightIonQMDCollision::~G4LightIonQMDCollision() 114 { 115 //if ( theSystem != NULL ) delete theSystem; 116 //if ( theMeanField != NULL ) delete theMeanField; 117 delete theScatterer; 118 } 119 120 121 void G4LightIonQMDCollision::CalKinematicsOfBinaryCollisions( G4double dt ) 122 { 123 G4double deltaT = dt; 124 125 G4int n = theSystem->GetTotalNumberOfParticipant(); 126 //081118 127 //G4int nb = 0; 128 for ( G4int i = 0 ; i < n ; i++ ) 129 { 130 theSystem->GetParticipant( i )->UnsetHitMark(); 131 theSystem->GetParticipant( i )->UnsetHitMark(); 132 //nb += theSystem->GetParticipant( i )->GetBaryonNumber(); 133 } 134 //G4cout << "nb = " << nb << " n = " << n << G4endl; 135 136 137 //071101 138 for ( G4int i = 0 ; i < n ; i++ ) 139 { 140 141 //std::cout << i << " " << theSystem->GetParticipant( i )->GetDefinition()->GetParticleName() << " " << theSystem->GetParticipant( i )->GetPosition() << std::endl; 142 143 if ( theSystem->GetParticipant( i )->GetDefinition()->IsShortLived() ) 144 { 145 146 G4bool decayed = false; 147 148 const G4ParticleDefinition* pd0 = theSystem->GetParticipant( i )->GetDefinition(); 149 G4ThreeVector p0 = theSystem->GetParticipant( i )->GetMomentum(); 150 G4ThreeVector r0 = theSystem->GetParticipant( i )->GetPosition(); 151 152 G4LorentzVector p40 = theSystem->GetParticipant( i )->Get4Momentum(); 153 154 G4double eini = theMeanField->GetTotalEnergy(); // R-JQMD, Skyrme, RMF 155 156 G4int n0 = theSystem->GetTotalNumberOfParticipant(); 157 G4int i0 = 0; 158 159 G4bool isThisEnergyOK = false; 160 161 G4int maximumNumberOfTrial=4; 162 for ( G4int ii = 0 ; ii < maximumNumberOfTrial ; ii++ ) 163 { 164 165 //G4LorentzVector p4 = theSystem->GetParticipant( i )->Get4Momentum(); 166 G4LorentzVector p400 = p40; 167 168 p400 *= GeV; 169 //G4KineticTrack kt( theSystem->GetParticipant( i )->GetDefinition() , 0.0 , (theSystem->GetParticipant( i )->GetPosition())*fermi , p4 ); 170 G4KineticTrack kt( pd0 , 0.0 , r0*fermi , p400 ); 171 //std::cout << "G4KineticTrack " << i << " " << kt.GetDefinition()->GetParticleName() << kt.GetPosition() << std::endl; 172 G4KineticTrackVector* secs = NULL; 173 secs = kt.Decay(); 174 G4int id = 0; 175 //G4double et = 0; 176 if ( secs ) 177 { 178 for ( G4KineticTrackVector::iterator it 179 = secs->begin() ; it != secs->end() ; it++ ) 180 { 181 /* 182 G4cout << "G4KineticTrack" 183 << " " << (*it)->GetDefinition()->GetParticleName() 184 << " " << (*it)->Get4Momentum() 185 << " " << (*it)->GetPosition()/fermi 186 << G4endl; 187 */ 188 if ( id == 0 ) 189 { 190 theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() ); 191 theSystem->GetParticipant( i )->SetMomentum( (*it)->Get4Momentum().v()/GeV ); 192 theSystem->GetParticipant( i )->SetPosition( (*it)->GetPosition()/fermi ); 193 //theMeanField->Cal2BodyQuantities( i ); 194 //et += (*it)->Get4Momentum().e()/GeV; 195 } 196 if ( id > 0 ) 197 { 198 // Append end; 199 theSystem->SetParticipant ( new G4QMDParticipant ( (*it)->GetDefinition() , (*it)->Get4Momentum().v()/GeV , (*it)->GetPosition()/fermi ) ); 200 //et += (*it)->Get4Momentum().e()/GeV; 201 if ( id > 1 ) 202 { 203 //081118 204 //G4cout << "G4LightIonQMDCollision id >2; id= " << id << G4endl; 205 } 206 } 207 id++; // number of daughter particles 208 209 delete *it; 210 } 211 212 theMeanField->Update(); 213 i0 = id-1; // 0 enter to i 214 215 delete secs; 216 } 217 218 // EnergyCheck 219 220 G4double efin = theMeanField->GetTotalEnergy(); // R-JQMD, Skyrme, RMF 221 //std::cout << std::abs ( eini - efin ) - fepse << std::endl; 222 // std::cout << std::abs ( eini - efin ) - fepse*10 << std::endl; 223 // *10 TK 224 if ( std::abs ( eini - efin ) < fepse*10 ) 225 { 226 // Energy OK 227 isThisEnergyOK = true; 228 break; 229 } 230 else 231 { 232 233 theSystem->GetParticipant( i )->SetDefinition( pd0 ); 234 theSystem->GetParticipant( i )->SetPosition( r0 ); 235 theSystem->GetParticipant( i )->SetMomentum( p0 ); 236 237 //for ( G4int i0i = 0 ; i0i < id-1 ; i0i++ ) 238 //160210 deletion must be done in descending order 239 for ( G4int i0i = id-2 ; 0 <= i0i ; i0i-- ) { 240 //081118 241 //std::cout << "Decay Energitically Blocked deleteing " << i0i+n0 << std::endl; 242 theSystem->DeleteParticipant( i0i+n0 ); 243 } 244 //081103 245 theMeanField->Update(); 246 } 247 248 } 249 250 251 // Pauli Check 252 if ( isThisEnergyOK == true ) 253 { 254 if ( theMeanField->IsPauliBlocked ( i ) != true ) 255 { 256 257 G4bool allOK = true; 258 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ ) 259 { 260 if ( theMeanField->IsPauliBlocked ( i0i+n0 ) == true ) 261 { 262 allOK = false; 263 break; 264 } 265 } 266 267 if ( allOK ) 268 { 269 decayed = true; //Decay Succeeded 270 } 271 } 272 273 } 274 // 275 276 if ( decayed ) 277 { 278 //081119 279 //G4cout << "Decay Suceeded! " << std::endl; 280 theSystem->GetParticipant( i )->SetHitMark(); 281 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ ) 282 { 283 theSystem->GetParticipant( i0i+n0 )->SetHitMark(); 284 } 285 286 } 287 else 288 { 289 290 // Decay Blocked and re-enter orginal participant; 291 292 if ( isThisEnergyOK == true ) // for false case already done 293 { 294 295 theSystem->GetParticipant( i )->SetDefinition( pd0 ); 296 theSystem->GetParticipant( i )->SetPosition( r0 ); 297 theSystem->GetParticipant( i )->SetMomentum( p0 ); 298 299 for ( G4int i0i = 0 ; i0i < i0 ; i0i++ ) 300 { 301 //081118 302 //std::cout << "Decay Blocked deleteing " << i0i+n0 << std::endl; 303 //160210 adding commnet: deletion must be done in descending order 304 theSystem->DeleteParticipant( i0+n0-i0i-1 ); 305 } 306 //081103 307 theMeanField->Update(); 308 } 309 310 } 311 312 } //shortlive 313 } // go next participant 314 //071101 315 316 317 n = theSystem->GetTotalNumberOfParticipant(); 318 319 //081118 320 //for ( G4int i = 1 ; i < n ; i++ ) 321 for ( G4int i = 1 ; i < theSystem->GetTotalNumberOfParticipant() ; i++ ) 322 { 323 324 //std::cout << "Collision i " << i << std::endl; 325 if ( theSystem->GetParticipant( i )->IsThisHit() ) continue; 326 327 328 G4ThreeVector ri = theSystem->GetParticipant( i )->GetPosition(); 329 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum(); 330 G4double rmi = theSystem->GetParticipant( i )->GetMass(); 331 const G4ParticleDefinition* pdi = theSystem->GetParticipant( i )->GetDefinition(); 332 //090331 gamma 333 if ( pdi->GetPDGMass() == 0.0 ) continue; 334 335 //std::cout << " p4i00 " << p4i << std::endl; 336 for ( G4int j = 0 ; j < i ; j++ ) 337 { 338 339 340 /* 341 G4cout << "Collision " << i << " " << theSystem->GetParticipant( i )->IsThisProjectile() << G4endl; 342 G4cout << "Collision " << j << " " << theSystem->GetParticipant( j )->IsThisProjectile() << G4endl; 343 G4cout << "Collision " << i << " " << theSystem->GetParticipant( i )->IsThisTarget() << G4endl; 344 G4cout << "Collision " << j << " " << theSystem->GetParticipant( j )->IsThisTarget() << G4endl; 345 */ 346 347 // Only 1 Collision allowed for each particle in a time step. 348 //081119 349 if ( theSystem->GetParticipant( i )->IsThisHit() ) continue; 350 if ( theSystem->GetParticipant( j )->IsThisHit() ) continue; 351 352 //std::cout << "Collision " << i << " " << j << std::endl; 353 354 // Do not allow collision between nucleons in target/projectile til its first collision. 355 if ( theSystem->GetParticipant( i )->IsThisProjectile() ) 356 { 357 if ( theSystem->GetParticipant( j )->IsThisProjectile() ) continue; 358 } 359 else if ( theSystem->GetParticipant( i )->IsThisTarget() ) 360 { 361 if ( theSystem->GetParticipant( j )->IsThisTarget() ) continue; 362 } 363 364 365 G4ThreeVector rj = theSystem->GetParticipant( j )->GetPosition(); 366 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum(); 367 G4double rmj = theSystem->GetParticipant( j )->GetMass(); 368 const G4ParticleDefinition* pdj = theSystem->GetParticipant( j )->GetDefinition(); 369 //090331 gamma 370 if ( pdj->GetPDGMass() == 0.0 ) continue; 371 372 G4double rr2 = theMeanField->GetRR2( i , j ); 373 374 // Here we assume elab (beam momentum less than 5 GeV/n ) 375 if ( rr2 > fdeltar*fdeltar ) continue; 376 377 //G4double s = (p4i+p4j)*(p4i+p4j); 378 //G4double srt = std::sqrt ( s ); 379 380 G4double srt = std::sqrt( (p4i+p4j)*(p4i+p4j) ); 381 382 G4double cutoff = 0.0; 383 G4double fbcmax = 0.0; 384 //110617 fix for gcc 4.6 compilation warnings 385 //G4double sig = 0.0; 386 387 if ( rmi < 0.94 && rmj < 0.94 ) 388 { 389 // nucleon or pion case 390 cutoff = rmi + rmj + 0.02; 391 fbcmax = fbcmax0; 392 //110617 fix for gcc 4.6 compilation warnings 393 //sig = sig0; 394 } 395 else 396 { 397 cutoff = rmi + rmj; 398 fbcmax = fbcmax1; 399 //110617 fix for gcc compilation warnings 400 //sig = sig1; 401 } 402 403 //std::cout << "Collision cutoff " << i << " " << j << " " << cutoff << std::endl; 404 if ( srt < cutoff ) continue; 405 406 G4ThreeVector dr = ri - rj; 407 G4double rsq = dr*dr; 408 409 G4double pij = p4i*p4j; 410 G4double pidr = p4i.vect()*dr; 411 G4double pjdr = p4j.vect()*dr; 412 413 G4double aij = 1.0 - ( rmi*rmj /pij ) * ( rmi*rmj /pij ); 414 G4double bij = pidr / rmi - pjdr*rmi/pij; 415 G4double cij = rsq + ( pidr / rmi ) * ( pidr / rmi ); 416 G4double brel = std::sqrt ( std::abs ( cij - bij*bij/aij ) ); 417 418 if ( brel > fbcmax ) continue; 419 //std::cout << "collisions3 " << std::endl; 420 421 G4double bji = -pjdr/rmj + pidr * rmj /pij; 422 423 G4double ti = ( pidr/rmi - bij / aij ) * p4i.e() / rmi; 424 G4double tj = (-pjdr/rmj - bji / aij ) * p4j.e() / rmj; 425 426 427 /* 428 G4cout << "collisions4 p4i " << p4i << G4endl; 429 G4cout << "collisions4 ri " << ri << G4endl; 430 G4cout << "collisions4 p4j " << p4j << G4endl; 431 G4cout << "collisions4 rj " << rj << G4endl; 432 G4cout << "collisions4 dr " << dr << G4endl; 433 G4cout << "collisions4 pij " << pij << G4endl; 434 G4cout << "collisions4 aij " << aij << G4endl; 435 G4cout << "collisions4 bij bji " << bij << " " << bji << G4endl; 436 G4cout << "collisions4 pidr pjdr " << pidr << " " << pjdr << G4endl; 437 G4cout << "collisions4 p4i.e() p4j.e() " << p4i.e() << " " << p4j.e() << G4endl; 438 G4cout << "collisions4 rmi rmj " << rmi << " " << rmj << G4endl; 439 G4cout << "collisions4 " << ti << " " << tj << G4endl; 440 */ 441 if ( std::abs ( ti + tj ) > deltaT ) continue; 442 //std::cout << "collisions4 " << std::endl; 443 444 G4ThreeVector beta = ( p4i + p4j ).boostVector(); 445 446 G4LorentzVector p = p4i; 447 G4LorentzVector p4icm = p.boost( p.findBoostToCM ( p4j ) ); 448 G4ThreeVector pcm = p4icm.vect(); 449 450 G4double prcm = pcm.mag(); 451 452 if ( prcm <= 0.00001 ) continue; 453 //std::cout << "collisions5 " << std::endl; 454 455 G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollision ( i , j ) ); // Use Geant4 Collision Library 456 //G4bool energetically_forbidden = !( CalFinalStateOfTheBinaryCollisionJQMD ( sig , cutoff , pcm , prcm , srt, beta , gamma , i , j ) ); // JQMD Elastic 457 458 /* 459 G4bool pauli_blocked = false; 460 if ( energetically_forbidden == false ) // result true 461 { 462 if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) 463 { 464 pauli_blocked = true; 465 //std::cout << "G4QMDRESULT Collsion Pauli Blocked " << std::endl; 466 } 467 } 468 else 469 { 470 if ( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) 471 pauli_blocked = false; 472 //std::cout << "G4QMDRESULT Collsion Blocked " << std::endl; 473 } 474 */ 475 476 /* 477 G4cout << "G4QMDRESULT Collsion initial p4 i and j " 478 << p4i << " " << p4j 479 << G4endl; 480 */ 481 // 081118 482 //if ( energetically_forbidden == true || pauli_blocked == true ) 483 if ( energetically_forbidden == true ) 484 { 485 486 //G4cout << " energetically_forbidden " << G4endl; 487 // Collsion not allowed then re enter orginal participants 488 // Now only momentum, becasuse we only consider elastic scattering of nucleons 489 490 theSystem->GetParticipant( i )->SetMomentum( p4i.vect() ); 491 theSystem->GetParticipant( i )->SetDefinition( pdi ); 492 theSystem->GetParticipant( i )->SetPosition( ri ); 493 494 theSystem->GetParticipant( j )->SetMomentum( p4j.vect() ); 495 theSystem->GetParticipant( j )->SetDefinition( pdj ); 496 theSystem->GetParticipant( j )->SetPosition( rj ); 497 498 theMeanField->Cal2BodyQuantities( i ); 499 theMeanField->Cal2BodyQuantities( j ); 500 501 } 502 else 503 { 504 505 506 G4bool absorption = false; 507 if ( n == theSystem->GetTotalNumberOfParticipant()+1 ) absorption = true; 508 if ( absorption ) 509 { 510 //G4cout << "Absorption happend " << G4endl; 511 i = i-1; 512 n = n-1; 513 } 514 515 // Collsion allowed (really happened) 516 517 // Unset Projectile/Target flag 518 theSystem->GetParticipant( i )->UnsetInitialMark(); 519 if ( !absorption ) theSystem->GetParticipant( j )->UnsetInitialMark(); 520 521 theSystem->GetParticipant( i )->SetHitMark(); 522 if ( !absorption ) theSystem->GetParticipant( j )->SetHitMark(); 523 524 theSystem->IncrementCollisionCounter(); 525 526 /* 527 G4cout << "G4QMDRESULT Collsion Really Happened between " 528 << i << " and " << j 529 << G4endl; 530 G4cout << "G4QMDRESULT Collsion initial p4 i and j " 531 << p4i << " " << p4j 532 << G4endl; 533 G4cout << "G4QMDRESULT Collsion after p4 i and j " 534 << theSystem->GetParticipant( i )->Get4Momentum() 535 << " " 536 << theSystem->GetParticipant( j )->Get4Momentum() 537 << G4endl; 538 G4cout << "G4QMDRESULT Collsion Diff " 539 << p4i + p4j - theSystem->GetParticipant( i )->Get4Momentum() - theSystem->GetParticipant( j )->Get4Momentum() 540 << G4endl; 541 G4cout << "G4QMDRESULT Collsion initial r i and j " 542 << ri << " " << rj 543 << G4endl; 544 G4cout << "G4QMDRESULT Collsion after r i and j " 545 << theSystem->GetParticipant( i )->GetPosition() 546 << " " 547 << theSystem->GetParticipant( j )->GetPosition() 548 << G4endl; 549 */ 550 551 552 } 553 554 } 555 556 } 557 558 559 } 560 561 562 563 G4bool G4LightIonQMDCollision::CalFinalStateOfTheBinaryCollision( G4int i , G4int j ) 564 { 565 566 //081103 567 //G4cout << "CalFinalStateOfTheBinaryCollision " << i << " " << j << " " << theSystem->GetTotalNumberOfParticipant() << G4endl; 568 569 G4int k = theSystem->GetTotalNumberOfParticipant(); // added by Y-H. Sato and A.H. 570 571 G4bool result = false; 572 G4bool energyOK = false; 573 G4bool pauliOK = false; 574 G4bool abs = false; 575 G4bool pion_prod = false; // added by Y-H. Sato and A.H. 576 //G4bool pion_abs = false; // added by Y-H. Sato and A.H. 577 G4QMDParticipant* absorbed = NULL; 578 579 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum(); 580 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum(); 581 582 //071031 583 584 //G4double epot = theMeanField->GetTotalPotential(); 585 //G4double eini = epot + p4i.e() + p4j.e(); 586 G4double eini = theMeanField->GetTotalEnergy(); // R-JQMD, Skyrme, RMF 587 588 //071031 589 // will use KineticTrack 590 const G4ParticleDefinition* pdi0 =theSystem->GetParticipant( i )->GetDefinition(); 591 const G4ParticleDefinition* pdj0 =theSystem->GetParticipant( j )->GetDefinition(); 592 G4LorentzVector p4i0 = p4i*GeV; 593 G4LorentzVector p4j0 = p4j*GeV; 594 G4ThreeVector ri0 = ( theSystem->GetParticipant( i )->GetPosition() )*fermi; 595 G4ThreeVector rj0 = ( theSystem->GetParticipant( j )->GetPosition() )*fermi; 596 597 for ( G4int iitry = 0 ; iitry < 4 ; iitry++ ) 598 { 599 600 abs = false; 601 602 G4KineticTrack kt1( pdi0 , 0.0 , ri0 , p4i0 ); 603 G4KineticTrack kt2( pdj0 , 0.0 , rj0 , p4j0 ); 604 605 G4LorentzVector p4ix_new; 606 G4LorentzVector p4jx_new; 607 G4LorentzVector p4kx_new(G4ThreeVector(0,0,0) , 0 ); // added by Y-H. S. and A.H. 608 G4KineticTrackVector* secs = NULL; 609 secs = theScatterer->Scatter( kt1 , kt2 ); 610 611 //std::cout << "G4QMDSCATTERER BEFORE " << kt1.GetDefinition()->GetParticleName() << " " << kt1.Get4Momentum()/GeV << " " << kt1.GetPosition()/fermi << std::endl; 612 //std::cout << "G4QMDSCATTERER BEFORE " << kt2.GetDefinition()->GetParticleName() << " " << kt2.Get4Momentum()/GeV << " " << kt2.GetPosition()/fermi << std::endl; 613 //std::cout << "THESCATTERER " << theScatterer->GetCrossSection ( kt1 , kt2 )/millibarn << " " << elastic << " " << sig << std::endl; 614 615 616 if ( secs ) 617 { 618 G4int iti = 0; 619 if ( secs->size() == 2 ) 620 { 621 for ( G4KineticTrackVector::iterator it 622 = secs->begin() ; it != secs->end() ; it++ ) 623 { 624 if ( iti == 0 ) 625 { 626 theSystem->GetParticipant( i )->SetDefinition( (*it)->GetDefinition() ); 627 p4ix_new = (*it)->Get4Momentum()/GeV; 628 //std::cout << "THESCATTERER " << (*it)->GetDefinition()->GetParticleName() << std::endl; 629 theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() ); 630 } 631 if ( iti == 1 ) 632 { 633 //theSystem->GetParticipant( j )->SetDefinition( (*it)->GetDefinition() ); 634 //p4jx_new = (*it)->Get4Momentum()/GeV; 635 //std::cout << "THESCATTERER " << p4jx_new.e()-p4jx_new.m() << std::endl; 636 //theSystem->GetParticipant( j )->SetMomentum( p4jx_new.v() ); 637 638 // added by Y-H. S and A.H. 639 if((*it)->GetDefinition()->IsShortLived()) 640 { 641 G4KineticTrackVector * dec = (*it)->Decay(); 642 643 G4int ita = 0; 644 for(G4KineticTrackVector::iterator jter=dec->begin(); jter != dec->end(); jter++) 645 { 646 //G4cout << "decay " << (*jter)->GetDefinition()->GetParticleName() << " " << (*jter)->Get4Momentum()/GeV << " " << (*jter)->GetDefinition()->GetBaryonNumber() << G4endl; 647 648 if(ita == 0) 649 { 650 theSystem->SetParticipant( new G4QMDParticipant( (*jter)->GetDefinition() , (*jter)->Get4Momentum().v()/GeV , (*jter)->GetPosition()/fermi ) ); 651 //G4cout << "decay " << (*jter)->GetDefinition()->GetParticleName() << " " << (*jter)->Get4Momentum()/GeV << G4endl; 652 //theMeanField->Update(); 653 //G4cout << "decay " << (*jter)->GetDefinition()->GetParticleName() << " " << theMeanField->GetTotalEnergy() << " " << eini << G4endl; 654 theSystem->GetParticipant( k )->SetDefinition( (*jter)->GetDefinition() ); 655 //theMeanField->Update(); 656 //G4cout << "decay " << (*jter)->GetDefinition()->GetParticleName() << " " << theSystem->GetParticipant( k )->GetNuc() << " " << theMeanField->GetTotalEnergy() << G4endl; 657 p4kx_new = (*jter)->Get4Momentum()/GeV; 658 theSystem->GetParticipant( k )->SetMomentum( p4kx_new.v() ); 659 //theSystem->ShowParticipants(); 660 661 } 662 else if(ita == 1) 663 { 664 theSystem->GetParticipant( j )->SetDefinition( (*jter)->GetDefinition() ); 665 p4jx_new = (*jter)->Get4Momentum()/GeV; 666 theSystem->GetParticipant( j )->SetMomentum( p4jx_new.v() ); 667 pion_prod = true; 668 //theMeanField->Update(); 669 //G4cout << "decay " << (*jter)->GetDefinition()->GetParticleName() << " " << theMeanField->GetTotalEnergy() << " " << eini << G4endl; 670 //theSystem->ShowParticipants(); 671 //exit(0); 672 } 673 else 674 { 675 std::cout << "************ Multi-particle decay ************" << std::endl; 676 } 677 ita ++; 678 679 } 680 delete dec; 681 682 //std::cout << "THESCATTERER " << "dec "<< dec << std::endl; 683 //std::cout << "THESCATTERER " << p4j.m() << " " << p4j << " " << p4i << std::endl; 684 //std::cout << "THESCATTERER " << pdi0->GetParticleName() << " " << pdj0->GetParticleName() << " " << p4i.e() + p4j.e() << std::endl; 685 //std::cout << "THESCATTERER " << theSystem->GetParticipant( k )->GetDefinition()->GetParticleName() << " " << p4kx_new.m() << " " << (*it)->Get4Momentum() << " " << theSystem->GetParticipant( k )->Get4Momentum() << " " << theSystem->GetParticipant( k )->Get4Momentum().m() << std::endl; 686 //std::cout << "THESCATTERER " << theSystem->GetParticipant( i )->GetDefinition()->GetParticleName() << " " << theSystem->GetParticipant( j )->GetDefinition()->GetParticleName() << " " << p4ix_new.e() + p4jx_new.e() << std::endl; 687 } 688 else 689 { 690 theSystem->GetParticipant( j )->SetDefinition( (*it)->GetDefinition() ); 691 p4jx_new = (*it)->Get4Momentum()/GeV; 692 //std::cout << "THESCATTERER " << p4jx_new.e()-p4jx_new.m() << std::endl; 693 theSystem->GetParticipant( j )->SetMomentum( p4jx_new.v() ); 694 } 695 696 } 697 //std::cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << std::endl; 698 iti++; 699 } 700 } 701 else if ( secs->size() == 1 ) 702 { 703 //081118 704 abs = true; 705 //G4cout << "G4LightIonQMDCollision pion absrorption " << secs->front()->GetDefinition()->GetParticleName() << G4endl; 706 707 // added by Y-H. S and A.H. 708 if(secs->front()->GetDefinition()->IsShortLived()) 709 { 710 G4KineticTrackVector * dec = secs->front()->Decay(); 711 G4int ita = 0; 712 for(G4KineticTrackVector::iterator jter=dec->begin(); jter != dec->end(); jter++) 713 { 714 //G4cout << "decay "<<(*jter)->GetDefinition()->GetParticleName()<< " " << (*jter)->Get4Momentum()/GeV << G4endl; 715 if(ita == 0) 716 { 717 theSystem->GetParticipant( i )->SetDefinition( (*jter)->GetDefinition() ); 718 p4ix_new = (*jter)->Get4Momentum()/GeV; 719 theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() ); 720 } 721 else if(ita == 1) 722 { 723 theSystem->GetParticipant( j )->SetDefinition( (*jter)->GetDefinition() ); 724 p4jx_new = (*jter)->Get4Momentum()/GeV; 725 theSystem->GetParticipant( j )->SetMomentum( p4jx_new.v() ); 726 abs = false; 727 //pion_abs = true; 728 } 729 else 730 { 731 std::cout << "************ Multi-particle decay ************" << std::endl; 732 } 733 ita ++; 734 } 735 delete dec; 736 } 737 else 738 { 739 theSystem->GetParticipant( i )->SetDefinition( secs->front()->GetDefinition() ); 740 p4ix_new = secs->front()->Get4Momentum()/GeV; 741 theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() ); 742 } 743 // added by Y-H. S and A.H. -- end 744 745 //secs->front()->Decay(); 746 /* 747 theSystem->GetParticipant( i )->SetDefinition( secs->front()->GetDefinition() ); 748 p4ix_new = secs->front()->Get4Momentum()/GeV; 749 theSystem->GetParticipant( i )->SetMomentum( p4ix_new.v() ); 750 */ 751 //std::cout << "THESCATTERER " << p4i.e()+p4j.e() << " " << p4j << " " << p4i << std::endl; 752 //std::cout << "THESCATTERER " << p4ix_new.e()+p4jx_new.e() << " " << p4ix_new << " " << p4jx_new << std::endl; 753 //exit(0); 754 } 755 756 //081118 757 if ( secs->size() > 2 ) 758 { 759 760 G4cout << "G4LightIonQMDCollision secs size > 2; " << secs->size() << G4endl; 761 762 for ( G4KineticTrackVector::iterator it 763 = secs->begin() ; it != secs->end() ; it++ ) 764 { 765 G4cout << "G4QMDSCATTERER AFTER " << (*it)->GetDefinition()->GetParticleName() << " " << (*it)->Get4Momentum()/GeV << G4endl; 766 } 767 768 } 769 770 // deleteing KineticTrack 771 for ( G4KineticTrackVector::iterator it 772 = secs->begin() ; it != secs->end() ; it++ ) 773 { 774 delete *it; 775 } 776 777 delete secs; 778 } 779 //071031 780 781 if ( !abs ) 782 { 783 //theMeanField->Cal2BodyQuantities( i ); 784 //theMeanField->Cal2BodyQuantities( j ); 785 theMeanField->Update(); 786 787 } 788 else 789 { 790 absorbed = theSystem->EraseParticipant( j ); 791 theMeanField->Update(); 792 } 793 794 //epot = theMeanField->GetTotalPotential(); 795 //G4double efin = epot + p4ix_new.e() + p4jx_new.e(); 796 G4double efin = theMeanField->GetTotalEnergy(); // R-JQMD, Skyrme, RMF 797 //std::cout << "Collision NEW epot " << i << " " << j << " " << epot << " " << std::abs ( eini - efin ) - fepse << std::endl; 798 799 /* 800 G4cout << "Collision efin " << i << " " << j << " " << efin << G4endl; 801 G4cout << "Collision " << i << " " << j << " " << std::abs ( eini - efin ) << " " << fepse << G4endl; 802 G4cout << "Collision " << std::abs ( eini - efin ) << " " << fepse << G4endl; 803 */ 804 805 //071031 806 //G4double fepse_change = fepse; // Added by Y-H. S and A.H. 807 //if(pion_prod || pion_abs) fepse_change = fepse*10; // Added by Y-H. S and A.H. 808 809 if ( std::abs ( eini - efin ) < fepse ) 810 { 811 // Collison OK 812 //std::cout << "collisions6" << std::endl; 813 //std::cout << "collisions before " << p4i << " " << p4j << std::endl; 814 //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl; 815 //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl; 816 //std::cout << "collisions before " << ri0/fermi << " " << rj0/fermi << std::endl; 817 //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl; 818 819 energyOK = true; 820 break; 821 } 822 else 823 { 824 //if(pion_prod || pion_abs) G4cout << "EnergyNotOK " << std::abs ( eini - efin )*1000 << " orig " << std::abs ( eini0 - efin0 )*1000 << G4endl; 825 //G4cout << "EnergyNotOK " << std::abs ( eini - efin )*1000 << G4endl; 826 //if(iitry == 3) G4cout << "Energy non-conserved and go through" << G4endl; 827 /* 828 if(std::abs ( eini - efin )*1000 > 20) 829 { 830 //G4cout << "EnergyNotOK " << std::abs ( eini - efin )*1000 << " orig " << std::abs ( eini0 - efin0 )*1000 << G4endl; 831 G4cout << p4ix_new + p4jx_new << " " << theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() << G4endl; 832 G4cout << p4ix_new.m() << " " << p4jx_new.m() << " " << theSystem->GetParticipant( i )->Get4Momentum().m() << " " << theSystem->GetParticipant( j )->Get4Momentum().m() << G4endl; 833 G4cout << "G4QMDRESULT Collsion Really Happened between " 834 << i << " and " << j 835 << G4endl; 836 G4cout << "G4QMDRESULT Collsion initial p4 i and j " 837 << p4i << " " << p4j 838 << G4endl; 839 G4cout << "G4QMDRESULT Collsion after p4 i and j " 840 << theSystem->GetParticipant( i )->Get4Momentum() 841 << " " 842 << theSystem->GetParticipant( j )->Get4Momentum() 843 << G4endl; 844 G4cout << "G4QMDRESULT Collsion Diff " 845 << p4i + p4j - theSystem->GetParticipant( i )->Get4Momentum() - theSystem->GetParticipant( j )->Get4Momentum() 846 << G4endl; 847 G4cout << "G4QMDRESULT Particle Name " 848 << theSystem->GetParticipant( i )->GetDefinition()->GetParticleName() << " " << theSystem->GetParticipant( j )->GetDefinition()->GetParticleName() << " mj " <<theSystem->GetParticipant( j )->Get4Momentum().m() 849 << G4endl; 850 //G4cout << "G4QMDRESULT Collsion initial r i and j " 851 //<< ri << " " << rj 852 //<< G4endl; 853 //G4cout << "G4QMDRESULT Collsion after r i and j " 854 //<< theSystem->GetParticipant( i )->GetPosition() 855 //<< " " 856 //<< theSystem->GetParticipant( j )->GetPosition() 857 //<< G4endl; 858 } 859 */ 860 861 if ( abs ) 862 { 863 //G4cout << "TKDB reinsert j " << G4endl; 864 theSystem->InsertParticipant( absorbed , j ); 865 theMeanField->Update(); 866 } 867 else if ( pion_prod ) // added by Y-H. S and A.H. 868 { 869 //G4cout << "TKDB reinsert j " << G4endl; 870 theSystem->EraseParticipant( k ); 871 theMeanField->Update(); 872 theSystem->GetParticipant( i )->SetDefinition( pdi0 ); 873 theSystem->GetParticipant( i )->SetMomentum( p4i.v() ); 874 theSystem->GetParticipant( j )->SetDefinition( pdj0 ); 875 theSystem->GetParticipant( j )->SetMomentum( p4j.v() ); 876 theMeanField->Update(); 877 pion_prod = false; 878 } 879 else 880 { 881 theSystem->GetParticipant( i )->SetDefinition( pdi0 ); 882 theSystem->GetParticipant( i )->SetMomentum( p4i.v() ); 883 884 theSystem->GetParticipant( j )->SetDefinition( pdj0 ); 885 theSystem->GetParticipant( j )->SetMomentum( p4j.v() ); 886 theMeanField->Update(); 887 } // added by Y-H. S and A.H. -- end 888 // do not need reinsert in no absroption case 889 } 890 //071031 891 } 892 893 // Energetically forbidden collision 894 895 if ( energyOK ) 896 { 897 // Pauli Check 898 //G4cout << "Pauli Checking " << theSystem->GetTotalNumberOfParticipant() << G4endl; 899 if ( !abs ) 900 { 901 if ( !( theMeanField->IsPauliBlocked ( i ) == true || theMeanField->IsPauliBlocked ( j ) == true ) ) 902 { 903 //G4cout << "Binary Collision Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl; 904 pauliOK = true; 905 } 906 } 907 else 908 { 909 //if ( theMeanField->IsPauliBlocked ( i ) == false ) 910 //090126 i-1 cause jth is erased 911 if ( theMeanField->IsPauliBlocked ( i-1 ) == false ) 912 { 913 //G4cout << "Absorption Happen " << theSystem->GetTotalNumberOfParticipant() << G4endl; 914 delete absorbed; 915 pauliOK = true; 916 } 917 } 918 919 920 if ( pauliOK ) 921 { 922 result = true; 923 } 924 else 925 { 926 //G4cout << "Pauli Blocked" << G4endl; 927 if ( abs ) 928 { 929 //G4cout << "TKDB reinsert j pauli block" << G4endl; 930 theSystem->InsertParticipant( absorbed , j ); 931 theMeanField->Update(); 932 } 933 else if ( pion_prod ) // added by Y-H. S and A.H. 934 { 935 //G4cout << "TKDB reinsert j " << G4endl; 936 theSystem->EraseParticipant( k ); 937 theMeanField->Update(); 938 theSystem->GetParticipant( i )->SetDefinition( pdi0 ); 939 theSystem->GetParticipant( i )->SetMomentum( p4i.v() ); 940 theSystem->GetParticipant( j )->SetDefinition( pdj0 ); 941 theSystem->GetParticipant( j )->SetMomentum( p4j.v() ); 942 theMeanField->Update(); 943 pion_prod = false; 944 } 945 else 946 { 947 theSystem->GetParticipant( i )->SetDefinition( pdi0 ); 948 theSystem->GetParticipant( i )->SetMomentum( p4i.v() ); 949 950 theSystem->GetParticipant( j )->SetDefinition( pdj0 ); 951 theSystem->GetParticipant( j )->SetMomentum( p4j.v() ); 952 theMeanField->Update(); 953 } // added by Y-H. S and A.H. -- end 954 } 955 } 956 957 return result; 958 959 } 960 961 962 963 G4bool G4LightIonQMDCollision::CalFinalStateOfTheBinaryCollisionJQMD( G4double sig , G4double cutoff , G4ThreeVector pcm , G4double prcm , G4double srt , G4ThreeVector beta , G4double gamma , G4int i , G4int j ) 964 { 965 966 //G4cout << "CalFinalStateOfTheBinaryCollisionJQMD" << G4endl; 967 968 G4bool result = true; 969 970 G4LorentzVector p4i = theSystem->GetParticipant( i )->Get4Momentum(); 971 G4double rmi = theSystem->GetParticipant( i )->GetMass(); 972 G4int zi = theSystem->GetParticipant( i )->GetChargeInUnitOfEplus(); 973 974 G4LorentzVector p4j = theSystem->GetParticipant( j )->Get4Momentum(); 975 G4double rmj = theSystem->GetParticipant( j )->GetMass(); 976 G4int zj = theSystem->GetParticipant( j )->GetChargeInUnitOfEplus(); 977 978 G4double pr = prcm; 979 980 G4double c2 = pcm.z()/pr; 981 982 G4double csrt = srt - cutoff; 983 984 //G4double pri = prcm; 985 //G4double prf = sqrt ( 0.25 * srt*srt -rm2 ); 986 987 G4double asrt = srt - rmi - rmj; 988 G4double pra = prcm; 989 990 991 992 G4double elastic = 0.0; 993 994 if ( zi == zj ) 995 { 996 if ( csrt < 0.4286 ) 997 { 998 elastic = 35.0 / ( 1. + csrt * 100.0 ) + 20.0; 999 } 1000 else 1001 { 1002 elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 ) 1003 * 2. / pi + 1.0 ) * 9.65 + 7.0; 1004 } 1005 } 1006 else 1007 { 1008 if ( csrt < 0.4286 ) 1009 { 1010 elastic = 28.0 / ( 1. + csrt * 100.0 ) + 27.0; 1011 } 1012 else 1013 { 1014 elastic = ( - std::atan( ( csrt - 0.4286 ) * 1.5 - 0.8 ) 1015 * 2. / pi + 1.0 ) * 12.34 + 10.0; 1016 } 1017 } 1018 1019 // std::cout << "Collision csrt " << i << " " << j << " " << csrt << std::endl; 1020 // std::cout << "Collision elstic " << i << " " << j << " " << elastic << std::endl; 1021 1022 1023 // std::cout << "Collision sig " << i << " " << j << " " << sig << std::endl; 1024 if ( G4UniformRand() > elastic / sig ) 1025 { 1026 //std::cout << "Inelastic " << std::endl; 1027 //std::cout << "elastic/sig " << elastic/sig << std::endl; 1028 return result; 1029 } 1030 else 1031 { 1032 //std::cout << "Elastic " << std::endl; 1033 } 1034 // std::cout << "Collision ELSTIC " << i << " " << j << std::endl; 1035 1036 1037 G4double as = G4Pow::GetInstance()->powN ( 3.65 * asrt , 6 ); 1038 G4double a = 6.0 * as / (1.0 + as); 1039 G4double ta = -2.0 * pra*pra; 1040 G4double x = G4UniformRand(); 1041 G4double t1 = G4Log( (1-x) * G4Exp(2.*a*ta) + x ) / a; 1042 G4double c1 = 1.0 - t1/ta; 1043 1044 if( std::abs(c1) > 1.0 ) c1 = 2.0 * x - 1.0; 1045 1046 /* 1047 G4cout << "Collision as " << i << " " << j << " " << as << G4endl; 1048 G4cout << "Collision a " << i << " " << j << " " << a << G4endl; 1049 G4cout << "Collision ta " << i << " " << j << " " << ta << G4endl; 1050 G4cout << "Collision x " << i << " " << j << " " << x << G4endl; 1051 G4cout << "Collision t1 " << i << " " << j << " " << t1 << G4endl; 1052 G4cout << "Collision c1 " << i << " " << j << " " << c1 << G4endl; 1053 */ 1054 t1 = 2.0*pi*G4UniformRand(); 1055 // std::cout << "Collision t1 " << i << " " << j << " " << t1 << std::endl; 1056 G4double t2 = 0.0; 1057 if ( pcm.x() == 0.0 && pcm.y() == 0 ) 1058 { 1059 t2 = 0.0; 1060 } 1061 else 1062 { 1063 t2 = std::atan2( pcm.y() , pcm.x() ); 1064 } 1065 // std::cout << "Collision t2 " << i << " " << j << " " << t2 << std::endl; 1066 1067 G4double s1 = std::sqrt ( 1.0 - c1*c1 ); 1068 G4double s2 = std::sqrt ( 1.0 - c2*c2 ); 1069 1070 G4double ct1 = std::cos(t1); 1071 G4double st1 = std::sin(t1); 1072 1073 G4double ct2 = std::cos(t2); 1074 G4double st2 = std::sin(t2); 1075 1076 G4double ss = c2*s1*ct1 + s2*c1; 1077 1078 pcm.setX( pr * ( ss*ct2 - s1*st1*st2) ); 1079 pcm.setY( pr * ( ss*st2 + s1*st1*ct2) ); 1080 pcm.setZ( pr * ( c1*c2 - s1*s2*ct1) ); 1081 1082 // std::cout << "Collision pcm " << i << " " << j << " " << pcm << std::endl; 1083 1084 G4double epot = theMeanField->GetTotalPotential(); 1085 1086 G4double eini = epot + p4i.e() + p4j.e(); 1087 G4double etwo = p4i.e() + p4j.e(); 1088 1089 /* 1090 G4cout << "Collision epot " << i << " " << j << " " << epot << G4endl; 1091 G4cout << "Collision eini " << i << " " << j << " " << eini << G4endl; 1092 G4cout << "Collision etwo " << i << " " << j << " " << etwo << G4endl; 1093 */ 1094 1095 1096 for ( G4int itry = 0 ; itry < 4 ; itry++ ) 1097 { 1098 1099 G4double eicm = std::sqrt ( rmi*rmi + pcm*pcm ); 1100 G4double pibeta = pcm*beta; 1101 1102 G4double trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + eicm ); 1103 1104 G4ThreeVector pi_new = beta*trans + pcm; 1105 1106 G4double ejcm = std::sqrt ( rmj*rmj + pcm*pcm ); 1107 trans = gamma * ( gamma * pibeta / ( gamma + 1 ) + ejcm ); 1108 1109 G4ThreeVector pj_new = beta*trans - pcm; 1110 1111 // 1112 // Delete old 1113 // Add new Particitipants 1114 // 1115 // Now only change momentum ( Beacuse we only have elastic sctter of nucleon 1116 // In future Definition also will be change 1117 // 1118 1119 theSystem->GetParticipant( i )->SetMomentum( pi_new ); 1120 theSystem->GetParticipant( j )->SetMomentum( pj_new ); 1121 1122 //G4double pi_new_e = (theSystem->GetParticipant( i )->Get4Momentum()).e(); 1123 //G4double pj_new_e = (theSystem->GetParticipant( j )->Get4Momentum()).e(); 1124 1125 theMeanField->Cal2BodyQuantities( i ); 1126 theMeanField->Cal2BodyQuantities( j ); 1127 1128 //epot = theMeanField->GetTotalPotential(); 1129 //G4double efin = epot + pi_new_e + pj_new_e ; 1130 G4double efin = theMeanField->GetTotalEnergy(); 1131 //std::cout << "Collision NEW epot " << i << " " << j << " " << epot << " " << std::abs ( eini - efin ) - fepse << std::endl; 1132 /* 1133 G4cout << "Collision efin " << i << " " << j << " " << efin << G4endl; 1134 G4cout << "Collision " << i << " " << j << " " << std::abs ( eini - efin ) << " " << fepse << G4endl; 1135 G4cout << "Collision " << std::abs ( eini - efin ) << " " << fepse << G4endl; 1136 */ 1137 1138 //071031 1139 if ( std::abs ( eini - efin ) < fepse ) 1140 { 1141 // Collison OK 1142 //std::cout << "collisions6" << std::endl; 1143 //std::cout << "collisions before " << p4i << " " << p4j << std::endl; 1144 //std::cout << "collisions after " << theSystem->GetParticipant( i )->Get4Momentum() << " " << theSystem->GetParticipant( j )->Get4Momentum() << std::endl; 1145 //std::cout << "collisions dif " << ( p4i + p4j ) - ( theSystem->GetParticipant( i )->Get4Momentum() + theSystem->GetParticipant( j )->Get4Momentum() ) << std::endl; 1146 //std::cout << "collisions before " << rix/fermi << " " << rjx/fermi << std::endl; 1147 //std::cout << "collisions after " << theSystem->GetParticipant( i )->GetPosition() << " " << theSystem->GetParticipant( j )->GetPosition() << std::endl; 1148 } 1149 //071031 1150 1151 if ( std::abs ( eini - efin ) < fepse ) return result; // Collison OK 1152 1153 G4double cona = ( eini - efin + etwo ) / gamma; 1154 G4double fac2 = 1.0 / ( 4.0 * cona*cona * pr*pr ) * 1155 ( ( cona*cona - ( rmi*rmi + rmj*rmj ) )*( cona*cona - ( rmi*rmi + rmj*rmj ) ) 1156 - 4.0 * rmi*rmi * rmj*rmj ); 1157 1158 if ( fac2 > 0 ) 1159 { 1160 G4double fact = std::sqrt ( fac2 ); 1161 pcm = fact*pcm; 1162 } 1163 1164 1165 } 1166 1167 // Energetically forbidden collision 1168 result = false; 1169 1170 return result; 1171 1172 } 1173