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 // 080505 Fixed and changed sampling method of impact parameter by T. Koi 27 // 080602 Fix memory leaks by T. Koi 28 // 080612 Delete unnecessary dependency and unused functions 29 // Change criterion of reaction by T. Koi 30 // 081107 Add UnUseGEM (then use the default channel of G4Evaporation) 31 // UseFrag (chage criterion of a inelastic reaction) 32 // Fix bug in nucleon projectiles by T. Koi 33 // 090122 Be8 -> Alpha + Alpha 34 // 090331 Change member shenXS and genspaXS object to pointer 35 // 091119 Fix for incidence of neutral particles 36 // 37 // 230306 Fix in the judgement of elasticLike_system for nucleon-nucleon, pion-nucleon collistion 38 // in line 450 by Y-H. Sato and A. Haga. 39 // 230306 Fix for nucleon deplication 40 // added system->Clear() in line 522 by Y-H. Sato and A. Haga. 41 // 230306 Allowing to simlate nucleon-nucleon, pion-nucleon scattering 42 // pion is accepted in the Ratherford parameter setting by Y-H. Sato and A. Haga. 43 // 44 #include "G4LightIonQMDReaction.hh" 45 #include "G4LightIonQMDNucleus.hh" 46 #include "G4LightIonQMDGroundStateNucleus.hh" 47 #include "G4Pow.hh" 48 #include "G4PhysicalConstants.hh" 49 #include "G4SystemOfUnits.hh" 50 #include "G4NistManager.hh" 51 52 #include "G4CrossSectionDataSetRegistry.hh" 53 #include "G4BGGPionElasticXS.hh" 54 #include "G4BGGPionInelasticXS.hh" 55 #include "G4VCrossSectionDataSet.hh" 56 #include "G4CrossSectionInelastic.hh" 57 #include "G4ComponentGGNuclNuclXsc.hh" 58 #include "G4PhysicsModelCatalog.hh" 59 60 // Fpr inelastic cross section check 61 #include "G4NuclearRadii.hh" 62 #include "G4HadronNucleonXsc.hh" 63 // test.csv (writting reaction data (particle, position, momentum)) 64 #include <iostream> 65 #include <fstream> 66 using std::endl; // *** 67 using std::ofstream; // *** 68 // -- test.csv 69 70 G4LightIonQMDReaction::G4LightIonQMDReaction() 71 : G4HadronicInteraction("LightIonQMDModel") 72 , system ( NULL ) 73 , deltaT ( 1 ) // in fsec (c=1) 74 , maxTime ( 100 ) // will have maxTime-th time step 75 , envelopF ( 1.05 ) // 10% for Peripheral reactions 76 , gem ( true ) 77 , frag ( false ) 78 , secID( -1 ) 79 { 80 G4cout << "G4LightIonQMDReaction::G4LightIonQMDReaction" << G4endl; 81 G4cout << "Recommended Energy of LightIonQMD: 30MeV/u - 500MeV/u" << G4endl; 82 83 theXS = new G4CrossSectionInelastic( new G4ComponentGGNuclNuclXsc ); 84 pipElNucXS = new G4BGGPionElasticXS(G4PionPlus::PionPlus() ); 85 pipElNucXS->BuildPhysicsTable(*(G4PionPlus::PionPlus() ) ); 86 87 pimElNucXS = new G4BGGPionElasticXS(G4PionMinus::PionMinus() ); 88 pimElNucXS->BuildPhysicsTable(*(G4PionMinus::PionMinus() ) ); 89 90 pipInelNucXS = new G4BGGPionInelasticXS(G4PionPlus::PionPlus() ); 91 pipInelNucXS->BuildPhysicsTable(*(G4PionPlus::PionPlus() ) ); 92 93 pimInelNucXS = new G4BGGPionInelasticXS(G4PionMinus::PionMinus() ); 94 pimInelNucXS->BuildPhysicsTable(*(G4PionMinus::PionMinus() ) ); 95 96 meanField = new G4LightIonQMDMeanField(); 97 collision = new G4LightIonQMDCollision(); 98 99 excitationHandler = new G4ExcitationHandler(); 100 setEvaporationCh(); 101 102 coulomb_collision_gamma_proj = 0.0; 103 coulomb_collision_rx_proj = 0.0; 104 coulomb_collision_rz_proj = 0.0; 105 coulomb_collision_px_proj = 0.0; 106 coulomb_collision_pz_proj = 0.0; 107 108 coulomb_collision_gamma_targ = 0.0; 109 coulomb_collision_rx_targ = 0.0; 110 coulomb_collision_rz_targ = 0.0; 111 coulomb_collision_px_targ = 0.0; 112 coulomb_collision_pz_targ = 0.0; 113 114 secID = G4PhysicsModelCatalog::GetModelID( "model_LightIonQMDModel" ); 115 } 116 117 118 G4LightIonQMDReaction::~G4LightIonQMDReaction() 119 { 120 delete excitationHandler; 121 delete collision; 122 delete meanField; 123 } 124 125 126 G4HadFinalState* G4LightIonQMDReaction::ApplyYourself( const G4HadProjectile & projectile , G4Nucleus & target ) 127 { 128 //G4cout << "G4LightIonQMDReaction::ApplyYourself" << G4endl; 129 130 theParticleChange.Clear(); 131 132 system = new G4QMDSystem; 133 134 G4int proj_Z = 0; 135 G4int proj_A = 0; 136 const G4ParticleDefinition* proj_pd = ( const G4ParticleDefinition* ) projectile.GetDefinition(); 137 if ( proj_pd->GetParticleType() == "nucleus" ) 138 { 139 proj_Z = proj_pd->GetAtomicNumber(); 140 proj_A = proj_pd->GetAtomicMass(); 141 } 142 else 143 { 144 proj_Z = (int)( proj_pd->GetPDGCharge()/eplus ); 145 proj_A = 1; 146 } 147 //G4int targ_Z = int ( target.GetZ() + 0.5 ); 148 //G4int targ_A = int ( target.GetN() + 0.5 ); 149 //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 150 G4int targ_Z = target.GetZ_asInt(); 151 G4int targ_A = target.GetA_asInt(); 152 const G4ParticleDefinition* targ_pd = G4IonTable::GetIonTable()->GetIon( targ_Z , targ_A , 0.0 ); 153 154 155 //G4NistManager* nistMan = G4NistManager::Instance(); 156 // G4Element* G4NistManager::FindOrBuildElement( targ_Z ); 157 158 const G4DynamicParticle* proj_dp = new G4DynamicParticle ( proj_pd , projectile.Get4Momentum() ); 159 //const G4Element* targ_ele = nistMan->FindOrBuildElement( targ_Z ); 160 //G4double aTemp = projectile.GetMaterial()->GetTemperature(); 161 162 // Glauber-Gribov nucleus-nucleus cross section does not have GetIsoCrossSection, 163 // therefore call GetElementCrossSection instead. 164 //G4double xs_0 = theXS->GetIsoCrossSection ( proj_dp , targ_Z , targ_A ); 165 G4double xs_0 = theXS->GetElementCrossSection( proj_dp , targ_Z , projectile.GetMaterial() ); 166 167 // When the projectile is a pion 168 if (proj_pd == G4PionPlus::PionPlus() ) { 169 xs_0 = pipElNucXS->GetElementCrossSection(proj_dp, targ_Z, projectile.GetMaterial() ) + 170 pipInelNucXS->GetElementCrossSection(proj_dp, targ_Z, projectile.GetMaterial() ); 171 } else if (proj_pd == G4PionMinus::PionMinus() ) { 172 xs_0 = pimElNucXS->GetElementCrossSection(proj_dp, targ_Z, projectile.GetMaterial() ) + 173 pimInelNucXS->GetElementCrossSection(proj_dp, targ_Z, projectile.GetMaterial() ); 174 } 175 176 //G4double xs_0 = genspaXS->GetCrossSection ( proj_dp , targ_ele , aTemp ); 177 //G4double xs_0 = theXS->GetCrossSection ( proj_dp , targ_ele , aTemp ); 178 //110822 179 180 G4double bmax_0 = std::sqrt( xs_0 / pi ); 181 //std::cout << "bmax_0 in fm (fermi) " << bmax_0/fermi << std::endl; 182 183 //delete proj_dp; 184 185 G4bool elastic = true; 186 187 std::vector< G4LightIonQMDNucleus* > nucleuses; // Secondary nuceluses 188 G4ThreeVector boostToReac; // ReactionSystem (CM or NN); 189 G4ThreeVector boostBackToLAB; // Reaction System to LAB; 190 191 G4LorentzVector targ4p( G4ThreeVector( 0.0 ) , targ_pd->GetPDGMass()/GeV ); 192 G4ThreeVector boostLABtoCM = targ4p.findBoostToCM( proj_dp->Get4Momentum()/GeV ); // CM of target and proj; 193 194 G4double p1 = proj_dp->GetMomentum().mag()/GeV/proj_A; 195 G4double m1 = proj_dp->GetDefinition()->GetPDGMass()/GeV/proj_A; 196 G4double e1 = std::sqrt( p1*p1 + m1*m1 ); 197 G4double e2 = targ_pd->GetPDGMass()/GeV/targ_A; 198 G4double beta_nn = -p1 / ( e1+e2 ); 199 200 G4ThreeVector boostLABtoNN ( 0. , 0. , beta_nn ); // CM of NN; 201 202 G4double beta_nncm = ( - boostLABtoCM.beta() + boostLABtoNN.beta() ) / ( 1 - boostLABtoCM.beta() * boostLABtoNN.beta() ) ; 203 204 //std::cout << targ4p << std::endl; 205 //std::cout << proj_dp->Get4Momentum()<< std::endl; 206 //std::cout << beta_nncm << std::endl; 207 G4ThreeVector boostNNtoCM( 0. , 0. , beta_nncm ); // 208 G4ThreeVector boostCMtoNN( 0. , 0. , -beta_nncm ); // 209 210 boostToReac = boostLABtoNN; 211 boostBackToLAB = -boostLABtoNN; 212 213 delete proj_dp; 214 G4int icounter = 0; 215 G4int icounter_max = 1024; 216 while ( elastic ) // Loop checking, 11.03.2015, T. Koi 217 { 218 icounter++; 219 if ( icounter > icounter_max ) { 220 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl; 221 break; 222 } 223 224 // impact parameter 225 //G4double bmax = 1.05*(bmax_0/fermi); // 10% for Peripheral reactions 226 G4double bmax = envelopF*(bmax_0/fermi); 227 G4double b = bmax * std::sqrt ( G4UniformRand() ); 228 //071112 229 //G4double b = 0; 230 //G4double b = bmax; 231 //G4double b = bmax/1.05 * 0.7 * G4UniformRand(); 232 233 //G4cout << "G4QMDRESULT bmax_0 = " << bmax_0/fermi << " fm, bmax = " << bmax << " fm , b = " << b << " fm " << G4endl; 234 235 G4double plab = projectile.GetTotalMomentum()/GeV; 236 G4double elab = ( projectile.GetKineticEnergy() + proj_pd->GetPDGMass() + targ_pd->GetPDGMass() )/GeV; 237 238 calcOffSetOfCollision( b , proj_pd , targ_pd , plab , elab , bmax , boostCMtoNN ); 239 240 // Projectile 241 G4LorentzVector proj4pLAB = projectile.Get4Momentum()/GeV; 242 243 G4LightIonQMDGroundStateNucleus* proj(NULL); 244 if ( projectile.GetDefinition()->GetParticleType() == "nucleus" 245 || projectile.GetDefinition()->GetParticleName() == "proton" 246 || projectile.GetDefinition()->GetParticleName() == "neutron" ) 247 { 248 249 proj_Z = proj_pd->GetAtomicNumber(); 250 proj_A = proj_pd->GetAtomicMass(); 251 proj = new G4LightIonQMDGroundStateNucleus( proj_Z , proj_A ); 252 //proj->ShowParticipants(); 253 254 255 meanField->SetSystem ( proj ); 256 if ( proj_A != 1 ) 257 { 258 proj->SetTotalPotential( meanField->GetTotalPotential() ); 259 proj->CalEnergyAndAngularMomentumInCM(); 260 } 261 } 262 263 // Target 264 //G4int iz = int ( target.GetZ() ); 265 //G4int ia = int ( target.GetN() ); 266 //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 267 G4int iz = int ( target.GetZ_asInt() ); 268 G4int ia = int ( target.GetA_asInt() ); 269 G4LightIonQMDGroundStateNucleus* targ = new G4LightIonQMDGroundStateNucleus( iz , ia ); 270 271 meanField->SetSystem (targ ); 272 if ( ia != 1 ) 273 { 274 targ->SetTotalPotential( meanField->GetTotalPotential() ); 275 targ->CalEnergyAndAngularMomentumInCM(); 276 } 277 278 //G4LorentzVector targ4p( G4ThreeVector( 0.0 ) , targ->GetNuclearMass()/GeV ); 279 // Boost Vector to CM 280 //boostToCM = targ4p.findBoostToCM( proj4pLAB ); 281 282 // Target 283 for ( G4int i = 0 ; i < targ->GetTotalNumberOfParticipant() ; i++ ) 284 { 285 286 G4ThreeVector p0 = targ->GetParticipant( i )->GetMomentum(); 287 G4ThreeVector r0 = targ->GetParticipant( i )->GetPosition(); 288 289 G4ThreeVector p ( p0.x() + coulomb_collision_px_targ 290 , p0.y() 291 , p0.z() * coulomb_collision_gamma_targ + coulomb_collision_pz_targ ); 292 293 G4ThreeVector r ( r0.x() + coulomb_collision_rx_targ 294 , r0.y() 295 , r0.z() / coulomb_collision_gamma_targ + coulomb_collision_rz_targ ); 296 297 system->SetParticipant( new G4QMDParticipant( targ->GetParticipant( i )->GetDefinition() , p , r ) ); 298 system->GetParticipant( i )->SetTarget(); 299 300 } 301 302 G4LorentzVector proj4pCM = CLHEP::boostOf ( proj4pLAB , boostToReac ); 303 G4LorentzVector targ4pCM = CLHEP::boostOf ( targ4p , boostToReac ); 304 305 // Projectile 306 //G4cout << "proj : " << proj << G4endl; 307 //if ( proj != NULL ) 308 if ( proj_A != 1 ) 309 { 310 311 // projectile is nucleus 312 313 for ( G4int i = 0 ; i < proj->GetTotalNumberOfParticipant() ; i++ ) 314 { 315 316 G4ThreeVector p0 = proj->GetParticipant( i )->GetMomentum(); 317 G4ThreeVector r0 = proj->GetParticipant( i )->GetPosition(); 318 319 G4ThreeVector p ( p0.x() + coulomb_collision_px_proj 320 , p0.y() 321 , p0.z() * coulomb_collision_gamma_proj + coulomb_collision_pz_proj ); 322 323 G4ThreeVector r ( r0.x() + coulomb_collision_rx_proj 324 , r0.y() 325 , r0.z() / coulomb_collision_gamma_proj + coulomb_collision_rz_proj ); 326 327 system->SetParticipant( new G4QMDParticipant( proj->GetParticipant( i )->GetDefinition() , p , r ) ); 328 system->GetParticipant ( i + targ->GetTotalNumberOfParticipant() )->SetProjectile(); 329 } 330 331 } 332 else 333 { 334 335 // projectile is particle 336 337 // avoid multiple set in "elastic" loop 338 //G4cout << "system Total Participants : " << system->GetTotalNumberOfParticipant() << ", target : " << targ->GetTotalNumberOfParticipant() << G4endl; 339 if ( system->GetTotalNumberOfParticipant() == targ->GetTotalNumberOfParticipant() ) 340 { 341 342 G4int i = targ->GetTotalNumberOfParticipant(); 343 344 G4ThreeVector p0( 0 ); 345 G4ThreeVector r0( 0 ); 346 347 G4ThreeVector p ( p0.x() + coulomb_collision_px_proj 348 , p0.y() 349 , p0.z() * coulomb_collision_gamma_proj + coulomb_collision_pz_proj ); 350 351 G4ThreeVector r ( r0.x() + coulomb_collision_rx_proj 352 , r0.y() 353 , r0.z() / coulomb_collision_gamma_proj + coulomb_collision_rz_proj ); 354 355 system->SetParticipant( new G4QMDParticipant( (G4ParticleDefinition*)projectile.GetDefinition() , p , r ) ); 356 // This is not important becase only 1 projectile particle. 357 system->GetParticipant ( i )->SetProjectile(); 358 } 359 360 } 361 //system->ShowParticipants(); 362 363 delete targ; 364 delete proj; 365 366 meanField->SetSystem ( system ); 367 collision->SetMeanField ( meanField ); 368 369 // Time Evolution 370 //std::cout << "Start time evolution " << std::endl; 371 //system->ShowParticipants(); 372 for ( G4int i = 0 ; i < maxTime ; i++ ) 373 { 374 //G4cout << " do Paropagate " << i << " th time step. " << G4endl; 375 meanField->DoPropagation( deltaT ); 376 //system->ShowParticipants(); 377 collision->CalKinematicsOfBinaryCollisions( deltaT ); 378 379 //if ( i / 10 * 10 == i ) 380 //{ 381 //G4cout << i << " th time step. " << G4endl; 382 //system->ShowParticipants(); 383 //} 384 //system->ShowParticipants(); 385 } 386 //system->ShowParticipants(); 387 388 389 //std::cout << "Doing Cluster Judgment " << std::endl; 390 391 nucleuses = meanField->DoClusterJudgment(); 392 393 // Elastic Judgment 394 395 G4int numberOfSecondary = int ( nucleuses.size() ) + system->GetTotalNumberOfParticipant(); 396 397 G4int sec_a_Z = 0; 398 G4int sec_a_A = 0; 399 const G4ParticleDefinition* sec_a_pd = NULL; 400 G4int sec_b_Z = 0; 401 G4int sec_b_A = 0; 402 const G4ParticleDefinition* sec_b_pd = NULL; 403 404 if ( numberOfSecondary == 2 ) 405 { 406 407 G4bool elasticLike_system = false; 408 if ( nucleuses.size() == 2 ) 409 { 410 411 sec_a_Z = nucleuses[0]->GetAtomicNumber(); 412 sec_a_A = nucleuses[0]->GetMassNumber(); 413 sec_b_Z = nucleuses[1]->GetAtomicNumber(); 414 sec_b_A = nucleuses[1]->GetMassNumber(); 415 416 if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_Z == targ_Z && sec_b_A == targ_A ) 417 || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_Z == proj_Z && sec_b_A == proj_A ) ) 418 { 419 elasticLike_system = true; 420 } 421 422 } 423 else if ( nucleuses.size() == 1 ) 424 { 425 426 sec_a_Z = nucleuses[0]->GetAtomicNumber(); 427 sec_a_A = nucleuses[0]->GetMassNumber(); 428 sec_b_pd = system->GetParticipant( 0 )->GetDefinition(); 429 430 if ( ( sec_a_Z == proj_Z && sec_a_A == proj_A && sec_b_pd == targ_pd ) 431 || ( sec_a_Z == targ_Z && sec_a_A == targ_A && sec_b_pd == proj_pd ) ) 432 { 433 elasticLike_system = true; 434 } 435 436 } 437 else 438 { 439 440 sec_a_pd = system->GetParticipant( 0 )->GetDefinition(); 441 sec_b_pd = system->GetParticipant( 1 )->GetDefinition(); 442 443 if ( ( sec_a_pd == proj_pd && sec_b_pd == targ_pd ) 444 || ( sec_a_pd == targ_pd && sec_b_pd == proj_pd ) ) 445 { 446 elasticLike_system = true; 447 } 448 // QMD should be inelastic collision, so that nucleon-nucleon collision should also be inelastic in this phase. by Y-H. S and A. H, Mar. 6, 2023. 449 if ( (proj_pd->GetParticleName() == "proton" && targ_pd->GetParticleName() == "proton") 450 || (proj_pd->GetParticleName() == "neutron" && targ_pd->GetParticleName() == "proton") 451 || (proj_pd->GetParticleName() == "pi+" && targ_pd->GetParticleName() == "proton") 452 || (proj_pd->GetParticleName() == "pi-" && targ_pd->GetParticleName() == "proton")) 453 { 454 elasticLike_system = false; 455 //G4cout << "elasticLike_system = false proton NOCollision " << system->GetNOCollision() << G4endl; 456 if ( system->GetNOCollision() == 1 || icounter+900 > icounter_max) elastic = false; 457 } 458 // Addition -- end 459 } 460 461 if ( elasticLike_system == true ) 462 { 463 464 G4bool elasticLike_energy = true; 465 // Cal ExcitationEnergy 466 for ( G4int i = 0 ; i < int ( nucleuses.size() ) ; i++ ) 467 { 468 469 //meanField->SetSystem( nucleuses[i] ); 470 meanField->SetNucleus( nucleuses[i] ); 471 //nucleuses[i]->SetTotalPotential( meanField->GetTotalPotential() ); 472 //nucleuses[i]->CalEnergyAndAngularMomentumInCM(); 473 474 if ( nucleuses[i]->GetExcitationEnergy()*GeV > 1.0*MeV ) elasticLike_energy = false; 475 476 } 477 478 // Check Collision 479 G4bool withCollision = true; 480 if ( system->GetNOCollision() == 0 ) withCollision = false; 481 482 // Final judegement for Inelasitc or Elastic; 483 // 484 // ElasticLike without Collision 485 //if ( elasticLike_energy == true && withCollision == false ) elastic = true; // ielst = 0 486 // ElasticLike with Collision 487 //if ( elasticLike_energy == true && withCollision == true ) elastic = true; // ielst = 1 488 // InelasticLike without Collision 489 //if ( elasticLike_energy == false ) elastic = false; // ielst = 2 490 if ( frag == true ) 491 if ( elasticLike_energy == false ) elastic = false; 492 // InelasticLike with Collision 493 if ( elasticLike_energy == false && withCollision == true ) elastic = false; // ielst = 3 494 495 } 496 497 } 498 else 499 { 500 501 // numberOfSecondary != 2 502 elastic = false; 503 504 } 505 506 //071115 507 //G4cout << elastic << G4endl; 508 // if elastic is true try again from sampling of impact parameter 509 510 if ( elastic == true ) 511 { 512 // delete this nucleues 513 for ( std::vector< G4LightIonQMDNucleus* >::iterator 514 it = nucleuses.begin() ; it != nucleuses.end() ; it++ ) 515 { 516 delete *it; 517 } 518 nucleuses.clear(); 519 // system->Clear() should be included here. Otherwise, the nucleon is repeatedly regstered if the nucleon is the projectile. by Y-H. S. and A. H, Mar. 6, 2023. 520 system->Clear(); 521 } 522 523 } 524 525 526 // Statical Decay Phase 527 528 for ( std::vector< G4LightIonQMDNucleus* >::iterator it 529 = nucleuses.begin() ; it != nucleuses.end() ; it++ ) 530 { 531 532 /* 533 G4cout << "G4QMDRESULT " 534 << (*it)->GetAtomicNumber() 535 << " " 536 << (*it)->GetMassNumber() 537 << " " 538 << (*it)->Get4Momentum() 539 << " " 540 << (*it)->Get4Momentum().vect() 541 << " " 542 << (*it)->Get4Momentum().restMass() 543 << " " 544 << (*it)->GetNuclearMass()/GeV 545 << G4endl; 546 */ 547 548 meanField->SetNucleus ( *it ); 549 550 if ( (*it)->GetAtomicNumber() == 0 // neutron cluster 551 || (*it)->GetAtomicNumber() == (*it)->GetMassNumber() ) // proton cluster 552 { 553 // push back system 554 for ( G4int i = 0 ; i < (*it)->GetTotalNumberOfParticipant() ; i++ ) 555 { 556 G4QMDParticipant* aP = new G4QMDParticipant( ( (*it)->GetParticipant( i ) )->GetDefinition() , ( (*it)->GetParticipant( i ) )->GetMomentum() , ( (*it)->GetParticipant( i ) )->GetPosition() ); 557 system->SetParticipant ( aP ); 558 } 559 continue; 560 } 561 562 G4double nucleus_e = std::sqrt ( G4Pow::GetInstance()->powN ( (*it)->GetNuclearMass()/GeV , 2 ) + G4Pow::GetInstance()->powN ( (*it)->Get4Momentum().vect().mag() , 2 ) ); 563 G4LorentzVector nucleus_p4CM ( (*it)->Get4Momentum().vect() , nucleus_e ); 564 565 // std::cout << "G4QMDRESULT nucleus deltaQ " << deltaQ << std::endl; 566 567 G4int ia = (*it)->GetMassNumber(); 568 G4int iz = (*it)->GetAtomicNumber(); 569 570 G4LorentzVector lv ( G4ThreeVector( 0.0 ) , (*it)->GetExcitationEnergy()*GeV + G4IonTable::GetIonTable()->GetIonMass( iz , ia ) ); 571 572 G4Fragment* aFragment = new G4Fragment( ia , iz , lv ); 573 574 G4ReactionProductVector* rv; 575 rv = excitationHandler->BreakItUp( *aFragment ); 576 G4bool notBreak = true; 577 for ( G4ReactionProductVector::iterator itt 578 = rv->begin() ; itt != rv->end() ; itt++ ) 579 { 580 581 notBreak = false; 582 // Secondary from this nucleus (*it) 583 const G4ParticleDefinition* pd = (*itt)->GetDefinition(); 584 585 G4LorentzVector p4 ( (*itt)->GetMomentum()/GeV , (*itt)->GetTotalEnergy()/GeV ); //in nucleus(*it) rest system 586 G4LorentzVector p4_CM = CLHEP::boostOf( p4 , -nucleus_p4CM.findBoostToCM() ); // Back to CM 587 G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); // Back to LAB 588 589 590 //090122 591 //theParticleChange.AddSecondary( dp ); 592 if ( !( pd->GetAtomicNumber() == 4 && pd->GetAtomicMass() == 8 ) ) 593 { 594 //G4cout << "pd out of notBreak loop : " << pd->GetParticleName() << G4endl; 595 G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV ); 596 theParticleChange.AddSecondary( dp ); 597 } 598 else 599 { 600 //Be8 -> Alpha + Alpha + Q 601 G4ThreeVector randomized_direction( G4UniformRand() , G4UniformRand() , G4UniformRand() ); 602 randomized_direction = randomized_direction.unit(); 603 G4double q_decay = (*itt)->GetMass() - 2*G4Alpha::Alpha()->GetPDGMass(); 604 G4double p_decay = std::sqrt ( G4Pow::GetInstance()->powN(G4Alpha::Alpha()->GetPDGMass()+q_decay/2,2) - G4Pow::GetInstance()->powN(G4Alpha::Alpha()->GetPDGMass() , 2 ) ); 605 G4LorentzVector p4_a1 ( p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 ); //in Be8 rest system 606 607 G4LorentzVector p4_a1_Be8 = CLHEP::boostOf ( p4_a1/GeV , -p4.findBoostToCM() ); 608 G4LorentzVector p4_a1_CM = CLHEP::boostOf ( p4_a1_Be8 , -nucleus_p4CM.findBoostToCM() ); 609 G4LorentzVector p4_a1_LAB = CLHEP::boostOf ( p4_a1_CM , boostBackToLAB ); 610 611 G4LorentzVector p4_a2 ( -p_decay*randomized_direction , G4Alpha::Alpha()->GetPDGMass()+q_decay/2 ); //in Be8 rest system 612 613 G4LorentzVector p4_a2_Be8 = CLHEP::boostOf ( p4_a2/GeV , -p4.findBoostToCM() ); 614 G4LorentzVector p4_a2_CM = CLHEP::boostOf ( p4_a2_Be8 , -nucleus_p4CM.findBoostToCM() ); 615 G4LorentzVector p4_a2_LAB = CLHEP::boostOf ( p4_a2_CM , boostBackToLAB ); 616 617 G4DynamicParticle* dp1 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a1_LAB*GeV ); 618 G4DynamicParticle* dp2 = new G4DynamicParticle( G4Alpha::Alpha() , p4_a2_LAB*GeV ); 619 theParticleChange.AddSecondary( dp1 ); 620 theParticleChange.AddSecondary( dp2 ); 621 } 622 //090122 623 624 /* 625 G4cout 626 << "Regist Secondary " 627 << (*itt)->GetDefinition()->GetParticleName() 628 << " " 629 << (*itt)->GetMomentum()/GeV 630 << " " 631 << (*itt)->GetKineticEnergy()/GeV 632 << " " 633 << (*itt)->GetMass()/GeV 634 << " " 635 << (*itt)->GetTotalEnergy()/GeV 636 << " " 637 << (*itt)->GetTotalEnergy()/GeV * (*itt)->GetTotalEnergy()/GeV 638 - (*itt)->GetMomentum()/GeV * (*itt)->GetMomentum()/GeV 639 << " " 640 << nucleus_p4CM.findBoostToCM() 641 << " " 642 << p4 643 << " " 644 << p4_CM 645 << " " 646 << p4_LAB 647 << G4endl; 648 */ 649 650 } 651 if ( notBreak == true ) 652 { 653 654 const G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon( (*it)->GetAtomicNumber() , (*it)->GetMassNumber(), (*it)->GetExcitationEnergy()*GeV ); 655 //G4cout << "pd in notBreak loop : " << pd->GetParticleName() << G4endl; 656 G4LorentzVector p4_CM = nucleus_p4CM; 657 G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); // Back to LAB 658 G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV ); 659 theParticleChange.AddSecondary( dp ); 660 661 } 662 663 for ( G4ReactionProductVector::iterator itt 664 = rv->begin() ; itt != rv->end() ; itt++ ) 665 { 666 delete *itt; 667 } 668 delete rv; 669 670 delete aFragment; 671 672 } 673 674 675 676 for ( G4int i = 0 ; i < system->GetTotalNumberOfParticipant() ; i++ ) 677 { 678 // Secondary particles 679 680 const G4ParticleDefinition* pd = system->GetParticipant( i )->GetDefinition(); 681 G4LorentzVector p4_CM = system->GetParticipant( i )->Get4Momentum(); 682 G4LorentzVector p4_LAB = CLHEP::boostOf( p4_CM , boostBackToLAB ); 683 G4DynamicParticle* dp = new G4DynamicParticle( pd , p4_LAB*GeV ); 684 theParticleChange.AddSecondary( dp ); 685 //G4cout << "In the last theParticleChange loop : " << pd->GetParticleName() << G4endl; 686 687 /* 688 G4cout << "G4QMDRESULT " 689 << "r" << i << " " << system->GetParticipant ( i ) -> GetPosition() << " " 690 << "p" << i << " " << system->GetParticipant ( i ) -> Get4Momentum() 691 << G4endl; 692 */ 693 694 } 695 696 for ( std::vector< G4LightIonQMDNucleus* >::iterator it 697 = nucleuses.begin() ; it != nucleuses.end() ; it++ ) 698 { 699 delete *it; // delete nulceuse 700 } 701 nucleuses.clear(); 702 703 system->Clear(); 704 delete system; 705 706 theParticleChange.SetStatusChange( stopAndKill ); 707 708 for (G4int i = 0; i < G4int(theParticleChange.GetNumberOfSecondaries() ); i++) 709 { 710 //G4cout << "Particle : " << theParticleChange.GetSecondary(i)->GetParticle()->GetParticleDefinition()->GetParticleName() << G4endl; 711 //G4cout << "KEnergy : " << theParticleChange.GetSecondary(i)->GetParticle()->GetKineticEnergy() << G4endl; 712 //G4cout << "modelID : " << theParticleChange.GetSecondary(i)->GetCreatorModelID() << G4endl; 713 theParticleChange.GetSecondary(i)->SetCreatorModelID(secID); 714 } 715 716 return &theParticleChange; 717 718 } 719 720 721 722 void G4LightIonQMDReaction::calcOffSetOfCollision( G4double b , 723 const G4ParticleDefinition* pd_proj , 724 const G4ParticleDefinition* pd_targ , 725 G4double ptot , G4double etot , G4double bmax , G4ThreeVector boostToCM ) 726 { 727 728 G4double mass_proj = pd_proj->GetPDGMass()/GeV; 729 G4double mass_targ = pd_targ->GetPDGMass()/GeV; 730 731 G4double stot = std::sqrt ( etot*etot - ptot*ptot ); 732 733 G4double pstt = std::sqrt ( ( stot*stot - ( mass_proj + mass_targ ) * ( mass_proj + mass_targ ) 734 ) * ( stot*stot - ( mass_proj - mass_targ ) * ( mass_proj - mass_targ ) ) ) 735 / ( 2.0 * stot ); 736 737 G4double pzcc = pstt; 738 G4double eccm = stot - ( mass_proj + mass_targ ); 739 740 G4int zp = 1; 741 G4int ap = 1; 742 if ( pd_proj->GetParticleType() == "nucleus" ) 743 { 744 zp = pd_proj->GetAtomicNumber(); 745 ap = pd_proj->GetAtomicMass(); 746 } 747 else 748 { 749 // proton, neutron, mesons 750 zp = int ( pd_proj->GetPDGCharge()/eplus + 0.5 ); 751 // ap = 1; 752 } 753 754 755 G4int zt = pd_targ->GetAtomicNumber(); 756 G4int at = pd_targ->GetAtomicMass(); 757 758 759 // Check the ramx0 value 760 //G4double rmax0 = 8.0; // T.K dicide parameter value // for low energy 761 G4double rmax0 = bmax + 4.0; 762 G4double rmax = std::sqrt( rmax0*rmax0 + b*b ); 763 764 G4double ccoul = 0.001439767; 765 G4double pcca = 1.0 - double ( zp * zt ) * ccoul / eccm / rmax - ( b / rmax )*( b / rmax ); 766 767 G4double pccf = std::sqrt( pcca ); 768 769 //Fix for neutral particles 770 G4double aas1 = 0.0; 771 G4double bbs = 0.0; 772 773 if ( zp != 0 ) 774 { 775 G4double aas = 2.0 * eccm * b / double ( zp * zt ) / ccoul; 776 bbs = 1.0 / std::sqrt ( 1.0 + aas*aas ); 777 aas1 = ( 1.0 + aas * b / rmax ) * bbs; 778 } 779 780 G4double cost = 0.0; 781 G4double sint = 0.0; 782 G4double thet1 = 0.0; 783 G4double thet2 = 0.0; 784 if ( 1.0 - aas1*aas1 <= 0 || 1.0 - bbs*bbs <= 0.0 ) 785 { 786 cost = 1.0; 787 sint = 0.0; 788 } 789 else 790 { 791 G4double aat1 = aas1 / std::sqrt ( 1.0 - aas1*aas1 ); 792 G4double aat2 = bbs / std::sqrt ( 1.0 - bbs*bbs ); 793 794 thet1 = std::atan ( aat1 ); 795 thet2 = std::atan ( aat2 ); 796 797 // TK enter to else block 798 G4double theta = thet1 - thet2; 799 cost = std::cos( theta ); 800 sint = std::sin( theta ); 801 } 802 803 G4double rzpr = -rmax * cost * ( mass_targ ) / ( mass_proj + mass_targ ); 804 G4double rzta = rmax * cost * ( mass_proj ) / ( mass_proj + mass_targ ); 805 806 G4double rxpr = rmax / 2.0 * sint; 807 808 G4double rxta = -rxpr; 809 810 811 G4double pzpc = pzcc * ( cost * pccf + sint * b / rmax ); 812 G4double pxpr = pzcc * ( -sint * pccf + cost * b / rmax ); 813 814 G4double pztc = - pzpc; 815 G4double pxta = - pxpr; 816 817 G4double epc = std::sqrt ( pzpc*pzpc + pxpr*pxpr + mass_proj*mass_proj ); 818 G4double etc = std::sqrt ( pztc*pztc + pxta*pxta + mass_targ*mass_targ ); 819 820 G4double pzpr = pzpc; 821 G4double pzta = pztc; 822 G4double epr = epc; 823 G4double eta = etc; 824 825 // CM -> NN 826 G4double gammacm = boostToCM.gamma(); 827 //G4double betacm = -boostToCM.beta(); 828 G4double betacm = boostToCM.z(); 829 pzpr = pzpc + betacm * gammacm * ( gammacm / ( 1. + gammacm ) * pzpc * betacm + epc ); 830 pzta = pztc + betacm * gammacm * ( gammacm / ( 1. + gammacm ) * pztc * betacm + etc ); 831 epr = gammacm * ( epc + betacm * pzpc ); 832 eta = gammacm * ( etc + betacm * pztc ); 833 834 //G4double betpr = pzpr / epr; 835 //G4double betta = pzta / eta; 836 837 G4double gammpr = epr / ( mass_proj ); 838 G4double gammta = eta / ( mass_targ ); 839 840 pzta = pzta / double ( at ); 841 pxta = pxta / double ( at ); 842 843 pzpr = pzpr / double ( ap ); 844 pxpr = pxpr / double ( ap ); 845 846 G4double zeroz = 0.0; 847 848 rzpr = rzpr -zeroz; 849 rzta = rzta -zeroz; 850 851 // Set results 852 coulomb_collision_gamma_proj = gammpr; 853 coulomb_collision_rx_proj = rxpr; 854 coulomb_collision_rz_proj = rzpr; 855 coulomb_collision_px_proj = pxpr; 856 coulomb_collision_pz_proj = pzpr; 857 858 coulomb_collision_gamma_targ = gammta; 859 coulomb_collision_rx_targ = rxta; 860 coulomb_collision_rz_targ = rzta; 861 coulomb_collision_px_targ = pxta; 862 coulomb_collision_pz_targ = pzta; 863 864 } 865 866 void G4LightIonQMDReaction::setEvaporationCh() 867 { 868 //fEvaporation - 8 default channels 869 //fCombined - 8 default + 60 GEM 870 //fGEM - 2 default + 66 GEM 871 G4DeexChannelType ctype = gem ? fGEM : fCombined; 872 excitationHandler->SetDeexChannelsType(ctype); 873 } 874 875 void G4LightIonQMDReaction::ModelDescription(std::ostream& outFile) const 876 { 877 outFile << "Lorentz covarianted Quantum Molecular Dynamics model for nucleus (particle) vs nucleus reactions\n"; 878 } 879