Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/hadronic/Hadr09/Hadr09.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/hadronic/Hadr09/Hadr09.cc (Version 11.3.0) and /examples/extended/hadronic/Hadr09/Hadr09.cc (Version 6.1)


  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 /// \file Hadr09.cc                               
 28 /// \brief Main program of the hadronic/Hadr09    
 29 //                                                
 30 //--------------------------------------------    
 31 // This program shows how to use the class Had    
 32 // The class HadronicGenerator is a kind of "h    
 33 // provides Geant4 final states (i.e. secondar    
 34 // hadron-nuclear inelastic collisions.           
 35 // Please see the class itself for more inform    
 36 //                                                
 37 // The use of the class Hadronic Generator is     
 38 // the constructor needs to be invoked only on    
 39 // of the Geant4 "physics case" to consider ("    
 40 // considered as default is the name is not sp    
 41 // method needs to be called at each collision    
 42 // collision (hadron, energy, direction, mater    
 43 // The class HadronicGenerator is expected to     
 44 // multi-threaded environment with "external"     
 45 // that are not necessarily managed by Geant4     
 46 // each thread should have its own instance of    
 47 //                                                
 48 // See the string "***LOOKHERE***" below for t    
 49 // of this example: the "physics case", the se    
 50 // which to sample the projectile (i.e. whethe    
 51 // hadron or an ion - in the case of hadron pr    
 52 // is possible from which to sample at each co    
 53 // ion projectile, only one type of ion needs     
 54 // the kinetic energy of the projectile (which    
 55 // an interval), whether the direction of the     
 56 // sampled at each collision, the target mater    
 57 // is possible, from which the target material    
 58 // collision, and then from this target materi    
 59 // will be chosen randomly by Geant4 itself),     
 60 // some information or not and how frequently.    
 61 // Once a well-defined type of hadron-nucleus     
 62 // inelastic collision has been chosen, the me    
 63 //   HadronicGenerator::GenerateInteraction       
 64 // returns the secondaries produced by that in    
 65 // of a G4VParticleChange object).                
 66 // Some information about this final-state is     
 67 //                                                
 68 // Usage:  Hadr09                                 
 69 //--------------------------------------------    
 70                                                   
 71 //....oooOO0OOooo........oooOO0OOooo........oo    
 72 //....oooOO0OOooo........oooOO0OOooo........oo    
 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........oo    
 95                                                   
 96 int main(int, char**)                             
 97 {                                                 
 98   G4cout << "=== Test of the HadronicGenerator    
 99                                                   
100   // Enable light hypernuclei and anti-hypernu    
101   G4HadronicParameters::Instance()->SetEnableH    
102                                                   
103   // See the HadronicGenerator class for the p    
104   // ( In short, it is the name of the Geant4     
105   //   the collision, with the possibility of     
106   //   a given energy interval, as in physics     
107   const G4String namePhysics = "FTFP_BERT";  /    
108   // const G4String namePhysics = "FTFP_BERT_A    
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    
120   // in the interval [minEnergy, maxEnergy].      
121   G4double minEnergy = 1.0 * CLHEP::GeV;  //**    
122   G4double maxEnergy = 30.0 * CLHEP::GeV;  //*    
123                                                   
124   const G4int numCollisions = 1000;  //***LOOK    
125                                                   
126   // Enable or disable the print out of this p    
127   // produced in each collisions is printed ou    
128   // collisions, the list of secondaries is pr    
129   const G4bool isPrintingEnabled = true;  //**    
130   const G4int printingGap = 100;  //***LOOKHER    
131                                                   
132   // Vector of Geant4 names of hadron projecti    
133   // (with uniform probability) for each colli    
134   // (note that the 6 light hypernuclei and an    
135   // other hadrons, not as generic ions).         
136   // Note: comment out the corresponding line     
137   std::vector<G4String> vecProjectiles;  //***    
138   vecProjectiles.push_back("pi-");                
139   // Note: vecProjectiles.push_back( "pi0" );     
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" );     
146   // Note: vecProjectiles.push_back( "eta_prim    
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"     
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_sig    
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 int    
175   if (namePhysics == "FTFP_BERT" || namePhysic    
176       || namePhysics == "QGSP_BIC" || namePhys    
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"     
186     // Note: vecProjectiles.push_back( "J/psi"    
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( "Upsilo    
196     vecProjectiles.push_back("lambda_c+");        
197     vecProjectiles.push_back("anti_lambda_c+")    
198     // Note: vecProjectiles.push_back( "sigma_    
199     // Note: vecProjectiles.push_back( "anti_s    
200     // Note: vecProjectiles.push_back( "sigma_    
201     // Note: vecProjectiles.push_back( "anti_s    
202     // Note: vecProjectiles.push_back( "sigma_    
203     // Note: vecProjectiles.push_back( "anti_s    
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_    
213     // Note: vecProjectiles.push_back( "anti_s    
214     // Note: vecProjectiles.push_back( "sigma_    
215     // Note: vecProjectiles.push_back( "sigma_    
216     // Note: vecProjectiles.push_back( "sigma_    
217     // Note: vecProjectiles.push_back( "anti_s    
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 hyp    
227   // are swtiched on, then only FTFP and INCL     
228   // of light hypernuclei, and only FTFP is ca    
229   // interactions of light anti-hypernuclei.      
230   if (G4HadronicParameters::Instance()->Enable    
231     if (namePhysics == "FTFP_BERT" || namePhys    
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("doublehyperdou    
240       vecProjectiles.push_back("hyperHe5");       
241     }                                             
242     if (namePhysics == "FTFP_BERT" || namePhys    
243       // Light anti-hypernuclei                   
244       vecProjectiles.push_back("anti_hypertrit    
245       vecProjectiles.push_back("anti_hyperalph    
246       vecProjectiles.push_back("anti_hyperH4")    
247       vecProjectiles.push_back("anti_doublehyp    
248       vecProjectiles.push_back("anti_doublehyp    
249       vecProjectiles.push_back("anti_hyperHe5"    
250     }                                             
251   }                                               
252                                                   
253   G4ParticleDefinition* projectileNucleus = nu    
254   G4GenericIon* gion = G4GenericIon::GenericIo    
255   gion->SetProcessManager(new G4ProcessManager    
256   G4ParticleTable* partTable = G4ParticleTable    
257   G4IonTable* ions = partTable->GetIonTable();    
258   partTable->SetReadiness();                      
259   ions->CreateAllIon();                           
260   ions->CreateAllIsomer();                        
261                                                   
262   const G4bool isProjectileIon = false;  //***    
263   if (isProjectileIon) {                          
264     minEnergy = 40.0 * 13.0 * CLHEP::GeV;  //*    
265     maxEnergy = 40.0 * 13.0 * CLHEP::GeV;  //*    
266     G4int ionZ = 18, ionA = 40;  //***LOOKHERE    
267     projectileNucleus = partTable->GetIonTable    
268   }                                               
269                                                   
270   // Vector of Geant4 NIST names of materials:    
271   // (with uniform probability) for each colli    
272   // Note: comment out the corresponding line     
273   //       or, vice versa, add a new line to e    
274   std::vector<G4String> vecMaterials;  //***LO    
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.    
289   const G4int numMaterials = vecMaterials.size    
290                                                   
291   G4cout << G4endl << "=================  Conf    
292          << "Model: " << namePhysics << G4endl    
293          << maxEnergy / CLHEP::GeV << " ] GeV"    
294          << "Number of collisions:  " << numCo    
295          << "Number of hadron projectiles: " <    
296          << "Number of materials:   " << numMa    
297          << "IsIonProjectile: " << (projectile    
298          << (projectileNucleus != nullptr ? pr    
299          << G4endl << "=======================    
300                                                   
301   CLHEP::Ranlux64Engine defaultEngine(1234567,    
302   CLHEP::HepRandom::setTheEngine(&defaultEngin    
303   G4int seed = time(NULL);                        
304   CLHEP::HepRandom::setTheSeed(seed);             
305   G4cout << G4endl << " Initial seed = " << se    
306                                                   
307   // Instanciate the HadronicGenerator providi    
308   HadronicGenerator* theHadronicGenerator = ne    
309   //******************************************    
310                                                   
311   if (theHadronicGenerator == nullptr) {          
312     G4cerr << "ERROR: theHadronicGenerator is     
313     return 1;                                     
314   }                                               
315   else if (!theHadronicGenerator->IsPhysicsCas    
316     G4cerr << "ERROR: this physics case is NOT    
317     return 2;                                     
318   }                                               
319                                                   
320   // Loop over the collisions                     
321   G4double rnd1, rnd2, rnd3, rnd4, rnd5, rnd6,    
322   G4VParticleChange* aChange = nullptr;           
323   for (G4int i = 0; i < numCollisions; ++i) {     
324     // Draw some random numbers to select the     
325     // projectile hadron, projectile kinetic e    
326     rnd1 = CLHEP::HepRandom::getTheEngine()->f    
327     rnd2 = CLHEP::HepRandom::getTheEngine()->f    
328     rnd3 = CLHEP::HepRandom::getTheEngine()->f    
329     rnd4 = CLHEP::HepRandom::getTheEngine()->f    
330     rnd5 = CLHEP::HepRandom::getTheEngine()->f    
331     rnd6 = CLHEP::HepRandom::getTheEngine()->f    
332     // Sample the projectile kinetic energy       
333     projectileEnergy = minEnergy + rnd1 * (max    
334     if (projectileEnergy <= 0.0) projectileEne    
335     // Sample the projectile direction            
336     normalization = 1.0 / std::sqrt(rnd2 * rnd    
337     const G4bool isOnSmearingDirection = true;    
338     G4ThreeVector aDirection = G4ThreeVector(0    
339     if (isOnSmearingDirection) {                  
340       aDirection = G4ThreeVector(normalization    
341     }                                             
342     // Sample the projectile hadron from the v    
343     G4int index_projectile = std::trunc(rnd5 *    
344     G4String nameProjectile = vecProjectiles[i    
345     G4ParticleDefinition* projectile = partTab    
346     if (projectileNucleus) {                      
347       nameProjectile = projectileNucleus->GetP    
348       projectile = projectileNucleus;             
349     }                                             
350     // Sample the target material from the vec    
351     // (Note: the target nucleus will be sampl    
352     G4int index_material = std::trunc(rnd6 * n    
353     G4String nameMaterial = vecMaterials[index    
354     G4Material* material = G4NistManager::Inst    
355     if (material == nullptr) {                    
356       G4cerr << "ERROR: Material " << nameMate    
357       return 3;                                   
358     }                                             
359     if (isPrintingEnabled) {                      
360       G4cout << "\t Collision " << i << " ; pr    
361       if (projectileNucleus) {                    
362         G4cout << " ; Ekin[MeV]/nucleon="         
363                << projectileEnergy                
364                     / static_cast<G4double>(st    
365       }                                           
366       else {                                      
367         G4cout << " ; Ekin[MeV]=" << projectil    
368       }                                           
369       G4cout << " ; direction=" << aDirection     
370     }                                             
371                                                   
372     // Call here the "hadronic generator" to g    
373     aChange = theHadronicGenerator->GenerateIn    
374       projectile, projectileEnergy,               
375       /* *************************************    
376                                                   
377     G4int nsec = aChange ? aChange->GetNumberO    
378     G4bool isPrintingOfSecondariesEnabled = fa    
379     if (isPrintingEnabled) {                      
380       G4cout << G4endl << "\t --> #secondaries    
381              << " ; impactParameter[fm]=" << t    
382              << " ; #projectileSpectatorNucleo    
383              << theHadronicGenerator->GetNumbe    
384              << " ; #targetSpectatorNucleons="    
385              << theHadronicGenerator->GetNumbe    
386              << " ; #NNcollisions=" << theHadr    
387       if (i % printingGap == 0) {                 
388         isPrintingOfSecondariesEnabled = true;    
389         G4cout << "\t \t List of produced seco    
390       }                                           
391     }                                             
392     // Loop over produced secondaries and even    
393     for (G4int j = 0; j < nsec; ++j) {            
394       const G4DynamicParticle* sec = aChange->    
395       if (isPrintingOfSecondariesEnabled) {       
396         G4cout << "\t \t \t j=" << j << "\t" <    
397                << "\t p=" << sec->Get4Momentum    
398       }                                           
399       delete aChange->GetSecondary(j);            
400     }                                             
401     if (aChange) aChange->Clear();                
402   }                                               
403                                                   
404   G4cout << G4endl << " Final random number =     
405          << G4endl << "=== End of test ===" <<    
406 }                                                 
407                                                   
408 //....oooOO0OOooo........oooOO0OOooo........oo    
409