Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/hadronic/FlukaCern/ProcessLevel/FinalState/HadNucIneEvents.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 ]

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