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 // 27 /// \file Hadr09.cc 28 /// \brief Main program of the hadronic/Hadr09 example 29 // 30 //------------------------------------------------------------------------ 31 // This program shows how to use the class Hadronic Generator. 32 // The class HadronicGenerator is a kind of "hadronic generator", i.e. 33 // provides Geant4 final states (i.e. secondary particles) produced by 34 // hadron-nuclear inelastic collisions. 35 // Please see the class itself for more information. 36 // 37 // The use of the class Hadronic Generator is very simple: 38 // the constructor needs to be invoked only once - specifying the name 39 // of the Geant4 "physics case" to consider ("FTFP_BERT" will be 40 // considered as default is the name is not specified) - and then one 41 // method needs to be called at each collision, specifying the type of 42 // collision (hadron, energy, direction, material) to be simulated. 43 // The class HadronicGenerator is expected to work also in a 44 // multi-threaded environment with "external" threads (i.e. threads 45 // that are not necessarily managed by Geant4 run-manager): 46 // each thread should have its own instance of the class. 47 // 48 // See the string "***LOOKHERE***" below for the setting of parameters 49 // of this example: the "physics case", the set of possibilities from 50 // which to sample the projectile (i.e. whether the projectile is a 51 // hadron or an ion - in the case of hadron projectile, a list of hadrons 52 // is possible from which to sample at each collision; in the case of 53 // ion projectile, only one type of ion needs to be specified), 54 // the kinetic energy of the projectile (which can be sampled within 55 // an interval), whether the direction of the projectile is fixed or 56 // sampled at each collision, the target material (a list of materials 57 // is possible, from which the target material can be sampled at each 58 // collision, and then from this target material, the target nucleus 59 // will be chosen randomly by Geant4 itself), and whether to print out 60 // some information or not and how frequently. 61 // Once a well-defined type of hadron-nucleus or nucleus-nucleus 62 // inelastic collision has been chosen, the method 63 // HadronicGenerator::GenerateInteraction 64 // returns the secondaries produced by that interaction (in the form 65 // of a G4VParticleChange object). 66 // Some information about this final-state is printed out as an example. 67 // 68 // Usage: Hadr09 69 //------------------------------------------------------------------------ 70 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 73 74 #include "CLHEP/Random/Randomize.h" 75 #include "CLHEP/Random/Ranlux64Engine.h" 76 #include "HadronicGenerator.hh" 77 78 #include "G4GenericIon.hh" 79 #include "G4HadronicParameters.hh" 80 #include "G4IonTable.hh" 81 #include "G4Material.hh" 82 #include "G4NistManager.hh" 83 #include "G4ParticleTable.hh" 84 #include "G4PhysicalConstants.hh" 85 #include "G4ProcessManager.hh" 86 #include "G4SystemOfUnits.hh" 87 #include "G4UnitsTable.hh" 88 #include "G4VParticleChange.hh" 89 #include "G4ios.hh" 90 #include "globals.hh" 91 92 #include <iomanip> 93 94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 95 96 int main(int, char**) 97 { 98 G4cout << "=== Test of the HadronicGenerator ===" << G4endl; 99 100 // Enable light hypernuclei and anti-hypernuclei 101 G4HadronicParameters::Instance()->SetEnableHyperNuclei(true); 102 103 // See the HadronicGenerator class for the possibilities and meaning of the "physics cases". 104 // ( In short, it is the name of the Geant4 hadronic model used for the simulation of 105 // the collision, with the possibility of having a transition between two models in 106 // a given energy interval, as in physics lists. ) 107 const G4String namePhysics = "FTFP_BERT"; //***LOOKHERE*** PHYSICS CASE 108 // const G4String namePhysics = "FTFP_BERT_ATL"; 109 // const G4String namePhysics = "QGSP_BERT"; 110 // const G4String namePhysics = "QGSP_BIC"; 111 // const G4String namePhysics = "FTFP_INCLXX"; 112 // const G4String namePhysics = "FTFP"; 113 // const G4String namePhysics = "QGSP"; 114 // const G4String namePhysics = "BERT"; 115 // const G4String namePhysics = "BIC"; 116 // const G4String namePhysics = "IonBIC"; 117 // const G4String namePhysics = "INCL"; 118 119 // The kinetic energy of the projectile will be sampled randomly, with flat probability 120 // in the interval [minEnergy, maxEnergy]. 121 G4double minEnergy = 1.0 * CLHEP::GeV; //***LOOKHERE*** HADRON PROJECTILE MIN Ekin 122 G4double maxEnergy = 30.0 * CLHEP::GeV; //***LOOKHERE*** HADRON PROJECTILE MAX Ekin 123 124 const G4int numCollisions = 1000; //***LOOKHERE*** NUMBER OF COLLISIONS 125 126 // Enable or disable the print out of this program: if enabled, the number of secondaries 127 // produced in each collisions is printed out; moreover, once every "printingGap" 128 // collisions, the list of secondaries is printed out. 129 const G4bool isPrintingEnabled = true; //***LOOKHERE*** PRINT OUT ON/OFF 130 const G4int printingGap = 100; //***LOOKHERE*** GAP IN PRINTING 131 132 // Vector of Geant4 names of hadron projectiles: one of this will be sampled randomly 133 // (with uniform probability) for each collision, when the projectile is not a generic ion 134 // (note that the 6 light hypernuclei and anti-hypernuclei are treated here as for the 135 // other hadrons, not as generic ions). 136 // Note: comment out the corresponding line in order to exclude a particle. 137 std::vector<G4String> vecProjectiles; //***LOOKHERE*** POSSIBLE HADRON PROJECTILES 138 vecProjectiles.push_back("pi-"); 139 // Note: vecProjectiles.push_back( "pi0" ); // Excluded because too short-lived 140 vecProjectiles.push_back("pi+"); 141 vecProjectiles.push_back("kaon-"); 142 vecProjectiles.push_back("kaon+"); 143 vecProjectiles.push_back("kaon0L"); 144 vecProjectiles.push_back("kaon0S"); 145 // Note: vecProjectiles.push_back( "eta" ); // Excluded because too short-lived 146 // Note: vecProjectiles.push_back( "eta_prime" ); // Excluded because too short-lived 147 vecProjectiles.push_back("proton"); 148 vecProjectiles.push_back("neutron"); 149 vecProjectiles.push_back("deuteron"); 150 vecProjectiles.push_back("triton"); 151 vecProjectiles.push_back("He3"); 152 vecProjectiles.push_back("alpha"); 153 vecProjectiles.push_back("lambda"); 154 vecProjectiles.push_back("sigma-"); 155 // Note: vecProjectiles.push_back( "sigma0" ); // Excluded because too short-lived 156 vecProjectiles.push_back("sigma+"); 157 vecProjectiles.push_back("xi-"); 158 vecProjectiles.push_back("xi0"); 159 vecProjectiles.push_back("omega-"); 160 vecProjectiles.push_back("anti_proton"); 161 vecProjectiles.push_back("anti_neutron"); 162 vecProjectiles.push_back("anti_lambda"); 163 vecProjectiles.push_back("anti_sigma-"); 164 // Note: vecProjectiles.push_back( "anti_sigma0" ); // Excluded because too short-lived 165 vecProjectiles.push_back("anti_sigma+"); 166 vecProjectiles.push_back("anti_xi-"); 167 vecProjectiles.push_back("anti_xi0"); 168 vecProjectiles.push_back("anti_omega-"); 169 vecProjectiles.push_back("anti_deuteron"); 170 vecProjectiles.push_back("anti_triton"); 171 vecProjectiles.push_back("anti_He3"); 172 vecProjectiles.push_back("anti_alpha"); 173 174 // Only FTFP and QGSP can handle nuclear interaction of charm and bottom hadrons 175 if (namePhysics == "FTFP_BERT" || namePhysics == "FTFP_BERT_ATL" || namePhysics == "QGSP_BERT" 176 || namePhysics == "QGSP_BIC" || namePhysics == "FTFP" || namePhysics == "QGSP") 177 { 178 // Charm and bottom hadrons 179 vecProjectiles.push_back("D+"); 180 vecProjectiles.push_back("D-"); 181 vecProjectiles.push_back("D0"); 182 vecProjectiles.push_back("anti_D0"); 183 vecProjectiles.push_back("Ds+"); 184 vecProjectiles.push_back("Ds-"); 185 // Note: vecProjectiles.push_back( "etac" ); // Excluded because too short-lived 186 // Note: vecProjectiles.push_back( "J/psi" ); // Excluded because too short-lived 187 vecProjectiles.push_back("B+"); 188 vecProjectiles.push_back("B-"); 189 vecProjectiles.push_back("B0"); 190 vecProjectiles.push_back("anti_B0"); 191 vecProjectiles.push_back("Bs0"); 192 vecProjectiles.push_back("anti_Bs0"); 193 vecProjectiles.push_back("Bc+"); 194 vecProjectiles.push_back("Bc-"); 195 // Note: vecProjectiles.push_back( "Upsilon" ); // Excluded because too short-lived 196 vecProjectiles.push_back("lambda_c+"); 197 vecProjectiles.push_back("anti_lambda_c+"); 198 // Note: vecProjectiles.push_back( "sigma_c+" ); // Excluded because too short-lived 199 // Note: vecProjectiles.push_back( "anti_sigma_c+" ); // Excluded because too short-lived 200 // Note: vecProjectiles.push_back( "sigma_c0" ); // Excluded because too short-lived 201 // Note: vecProjectiles.push_back( "anti_sigma_c0" ); // Excluded because too short-lived 202 // Note: vecProjectiles.push_back( "sigma_c++" ); // Excluded because too short-lived 203 // Note: vecProjectiles.push_back( "anti_sigma_c++" ); // Excluded because too short-lived 204 vecProjectiles.push_back("xi_c+"); 205 vecProjectiles.push_back("anti_xi_c+"); 206 vecProjectiles.push_back("xi_c0"); 207 vecProjectiles.push_back("anti_xi_c0"); 208 vecProjectiles.push_back("omega_c0"); 209 vecProjectiles.push_back("anti_omega_c0"); 210 vecProjectiles.push_back("lambda_b"); 211 vecProjectiles.push_back("anti_lambda_b"); 212 // Note: vecProjectiles.push_back( "sigma_b+" ); // Excluded because too short-lived 213 // Note: vecProjectiles.push_back( "anti_sigma_b+" ); // Excluded because too short-lived 214 // Note: vecProjectiles.push_back( "sigma_b0" ); // Excluded because too short-lived 215 // Note: vecProjectiles.push_back( "sigma_b0" ); // Excluded because too short-lived 216 // Note: vecProjectiles.push_back( "sigma_b-" ); // Excluded because too short-lived 217 // Note: vecProjectiles.push_back( "anti_sigma_b-" ); // Excluded because too short-lived 218 vecProjectiles.push_back("xi_b0"); 219 vecProjectiles.push_back("anti_xi_b0"); 220 vecProjectiles.push_back("xi_b-"); 221 vecProjectiles.push_back("anti_xi_b-"); 222 vecProjectiles.push_back("omega_b-"); 223 vecProjectiles.push_back("anti_omega_b-"); 224 } 225 226 // If the hadronic interactions of light hypernuclei and anti-hypernuclei 227 // are swtiched on, then only FTFP and INCL can handle the nuclear interactions 228 // of light hypernuclei, and only FTFP is capable of handling the nuclear 229 // interactions of light anti-hypernuclei. 230 if (G4HadronicParameters::Instance()->EnableHyperNuclei()) { 231 if (namePhysics == "FTFP_BERT" || namePhysics == "FTFP_INCLXX" || namePhysics == "FTFP" 232 || namePhysics == "INCL") 233 { 234 // Light hypernuclei 235 vecProjectiles.push_back("hypertriton"); 236 vecProjectiles.push_back("hyperalpha"); 237 vecProjectiles.push_back("hyperH4"); 238 vecProjectiles.push_back("doublehyperH4"); 239 vecProjectiles.push_back("doublehyperdoubleneutron"); 240 vecProjectiles.push_back("hyperHe5"); 241 } 242 if (namePhysics == "FTFP_BERT" || namePhysics == "FTFP_INCLXX" || namePhysics == "FTFP") { 243 // Light anti-hypernuclei 244 vecProjectiles.push_back("anti_hypertriton"); 245 vecProjectiles.push_back("anti_hyperalpha"); 246 vecProjectiles.push_back("anti_hyperH4"); 247 vecProjectiles.push_back("anti_doublehyperH4"); 248 vecProjectiles.push_back("anti_doublehyperdoubleneutron"); 249 vecProjectiles.push_back("anti_hyperHe5"); 250 } 251 } 252 253 G4ParticleDefinition* projectileNucleus = nullptr; 254 G4GenericIon* gion = G4GenericIon::GenericIon(); 255 gion->SetProcessManager(new G4ProcessManager(gion)); 256 G4ParticleTable* partTable = G4ParticleTable::GetParticleTable(); 257 G4IonTable* ions = partTable->GetIonTable(); 258 partTable->SetReadiness(); 259 ions->CreateAllIon(); 260 ions->CreateAllIsomer(); 261 262 const G4bool isProjectileIon = false; //***LOOKHERE*** HADRON (false) OR ION (true) PROJECTILE? 263 if (isProjectileIon) { 264 minEnergy = 40.0 * 13.0 * CLHEP::GeV; //***LOOKHERE*** ION PROJECTILE MIN Ekin 265 maxEnergy = 40.0 * 13.0 * CLHEP::GeV; //***LOOKHERE*** ION PROJECTILE MAX Ekin 266 G4int ionZ = 18, ionA = 40; //***LOOKHERE*** ION PROJECTILE (Z, A) 267 projectileNucleus = partTable->GetIonTable()->GetIon(ionZ, ionA, 0.0); 268 } 269 270 // Vector of Geant4 NIST names of materials: one of this will be sampled randomly 271 // (with uniform probability) for each collision and used as target material. 272 // Note: comment out the corresponding line in order to exclude a material; 273 // or, vice versa, add a new line to extend the list with another material. 274 std::vector<G4String> vecMaterials; //***LOOKHERE*** : NIST TARGET MATERIALS 275 vecMaterials.push_back("G4_H"); 276 vecMaterials.push_back("G4_He"); 277 vecMaterials.push_back("G4_Be"); 278 vecMaterials.push_back("G4_C"); 279 vecMaterials.push_back("G4_Al"); 280 vecMaterials.push_back("G4_Si"); 281 // vecMaterials.push_back( "G4_Sc" ); 282 vecMaterials.push_back("G4_Ar"); 283 vecMaterials.push_back("G4_Fe"); 284 vecMaterials.push_back("G4_Cu"); 285 vecMaterials.push_back("G4_W"); 286 vecMaterials.push_back("G4_Pb"); 287 288 const G4int numProjectiles = vecProjectiles.size(); 289 const G4int numMaterials = vecMaterials.size(); 290 291 G4cout << G4endl << "================= Configuration ==================" << G4endl 292 << "Model: " << namePhysics << G4endl << "Ekin: [ " << minEnergy / CLHEP::GeV << " , " 293 << maxEnergy / CLHEP::GeV << " ] GeV" << G4endl 294 << "Number of collisions: " << numCollisions << G4endl 295 << "Number of hadron projectiles: " << numProjectiles << G4endl 296 << "Number of materials: " << numMaterials << G4endl 297 << "IsIonProjectile: " << (projectileNucleus != nullptr ? "true \t" : "false") 298 << (projectileNucleus != nullptr ? projectileNucleus->GetParticleName() : G4String("")) 299 << G4endl << "===================================================" << G4endl << G4endl; 300 301 CLHEP::Ranlux64Engine defaultEngine(1234567, 4); 302 CLHEP::HepRandom::setTheEngine(&defaultEngine); 303 G4int seed = time(NULL); 304 CLHEP::HepRandom::setTheSeed(seed); 305 G4cout << G4endl << " Initial seed = " << seed << G4endl << G4endl; 306 307 // Instanciate the HadronicGenerator providing the name of the "physics case" 308 HadronicGenerator* theHadronicGenerator = new HadronicGenerator(namePhysics); 309 //**************************************************************************** 310 311 if (theHadronicGenerator == nullptr) { 312 G4cerr << "ERROR: theHadronicGenerator is NULL !" << G4endl; 313 return 1; 314 } 315 else if (!theHadronicGenerator->IsPhysicsCaseSupported()) { 316 G4cerr << "ERROR: this physics case is NOT supported !" << G4endl; 317 return 2; 318 } 319 320 // Loop over the collisions 321 G4double rnd1, rnd2, rnd3, rnd4, rnd5, rnd6, normalization, projectileEnergy; 322 G4VParticleChange* aChange = nullptr; 323 for (G4int i = 0; i < numCollisions; ++i) { 324 // Draw some random numbers to select the hadron-nucleus interaction: 325 // projectile hadron, projectile kinetic energy, projectile direction, and target material. 326 rnd1 = CLHEP::HepRandom::getTheEngine()->flat(); 327 rnd2 = CLHEP::HepRandom::getTheEngine()->flat(); 328 rnd3 = CLHEP::HepRandom::getTheEngine()->flat(); 329 rnd4 = CLHEP::HepRandom::getTheEngine()->flat(); 330 rnd5 = CLHEP::HepRandom::getTheEngine()->flat(); 331 rnd6 = CLHEP::HepRandom::getTheEngine()->flat(); 332 // Sample the projectile kinetic energy 333 projectileEnergy = minEnergy + rnd1 * (maxEnergy - minEnergy); 334 if (projectileEnergy <= 0.0) projectileEnergy = minEnergy; 335 // Sample the projectile direction 336 normalization = 1.0 / std::sqrt(rnd2 * rnd2 + rnd3 * rnd3 + rnd4 * rnd4); 337 const G4bool isOnSmearingDirection = true; //***LOOKHERE*** 338 G4ThreeVector aDirection = G4ThreeVector(0.0, 0.0, 1.0); //***LOOKHERE*** 339 if (isOnSmearingDirection) { 340 aDirection = G4ThreeVector(normalization * rnd2, normalization * rnd3, normalization * rnd4); 341 } 342 // Sample the projectile hadron from the vector vecProjectiles 343 G4int index_projectile = std::trunc(rnd5 * numProjectiles); 344 G4String nameProjectile = vecProjectiles[index_projectile]; 345 G4ParticleDefinition* projectile = partTable->FindParticle(nameProjectile); 346 if (projectileNucleus) { 347 nameProjectile = projectileNucleus->GetParticleName(); 348 projectile = projectileNucleus; 349 } 350 // Sample the target material from the vector vecMaterials 351 // (Note: the target nucleus will be sampled by Geant4) 352 G4int index_material = std::trunc(rnd6 * numMaterials); 353 G4String nameMaterial = vecMaterials[index_material]; 354 G4Material* material = G4NistManager::Instance()->FindOrBuildMaterial(nameMaterial); 355 if (material == nullptr) { 356 G4cerr << "ERROR: Material " << nameMaterial << " is not found !" << G4endl; 357 return 3; 358 } 359 if (isPrintingEnabled) { 360 G4cout << "\t Collision " << i << " ; projectile=" << nameProjectile; 361 if (projectileNucleus) { 362 G4cout << " ; Ekin[MeV]/nucleon=" 363 << projectileEnergy 364 / static_cast<G4double>(std::abs(projectileNucleus->GetBaryonNumber())); 365 } 366 else { 367 G4cout << " ; Ekin[MeV]=" << projectileEnergy; 368 } 369 G4cout << " ; direction=" << aDirection << " ; material=" << nameMaterial; 370 } 371 372 // Call here the "hadronic generator" to get the secondaries produced by the hadronic collision 373 aChange = theHadronicGenerator->GenerateInteraction( 374 projectile, projectileEnergy, 375 /* ********************************************** */ aDirection, material); 376 377 G4int nsec = aChange ? aChange->GetNumberOfSecondaries() : 0; 378 G4bool isPrintingOfSecondariesEnabled = false; 379 if (isPrintingEnabled) { 380 G4cout << G4endl << "\t --> #secondaries=" << nsec 381 << " ; impactParameter[fm]=" << theHadronicGenerator->GetImpactParameter() / fermi 382 << " ; #projectileSpectatorNucleons=" 383 << theHadronicGenerator->GetNumberOfProjectileSpectatorNucleons() 384 << " ; #targetSpectatorNucleons=" 385 << theHadronicGenerator->GetNumberOfTargetSpectatorNucleons() 386 << " ; #NNcollisions=" << theHadronicGenerator->GetNumberOfNNcollisions() << G4endl; 387 if (i % printingGap == 0) { 388 isPrintingOfSecondariesEnabled = true; 389 G4cout << "\t \t List of produced secondaries: " << G4endl; 390 } 391 } 392 // Loop over produced secondaries and eventually print out some information. 393 for (G4int j = 0; j < nsec; ++j) { 394 const G4DynamicParticle* sec = aChange->GetSecondary(j)->GetDynamicParticle(); 395 if (isPrintingOfSecondariesEnabled) { 396 G4cout << "\t \t \t j=" << j << "\t" << sec->GetDefinition()->GetParticleName() 397 << "\t p=" << sec->Get4Momentum() << " MeV" << G4endl; 398 } 399 delete aChange->GetSecondary(j); 400 } 401 if (aChange) aChange->Clear(); 402 } 403 404 G4cout << G4endl << " Final random number = " << CLHEP::HepRandom::getTheEngine()->flat() 405 << G4endl << "=== End of test ===" << G4endl; 406 } 407 408 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 409