| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 /// \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