Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-04-01 07:50:18 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file RunAction.cc
0027 /// \brief Implementation of the RunAction class
0028 
0029 #include "RunAction.hh"
0030 
0031 #include "DetectorConstruction.hh"
0032 #include "PrimaryGeneratorAction.hh"
0033 
0034 #include "G4Electron.hh"
0035 #include "G4EmCalculator.hh"
0036 #include "G4LossTableManager.hh"
0037 #include "G4PhysicalConstants.hh"
0038 #include "G4Positron.hh"
0039 #include "G4ProcessManager.hh"
0040 #include "G4Run.hh"
0041 #include "G4SystemOfUnits.hh"
0042 #include "G4UnitsTable.hh"
0043 
0044 #include <vector>
0045 
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0047 
0048 RunAction::RunAction(DetectorConstruction* det, PrimaryGeneratorAction* kin)
0049   : fDetector(det), fPrimary(kin)
0050 {}
0051 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0052 
0053 void RunAction::BeginOfRunAction(const G4Run*)
0054 {
0055   // set precision for printing
0056   G4int prec = G4cout.precision(6);
0057 
0058   // instanciate EmCalculator
0059   G4EmCalculator emCal;
0060   //  emCal.SetVerbose(2);
0061 
0062   // get particle
0063   G4ParticleDefinition* particle = fPrimary->GetParticleGun()->GetParticleDefinition();
0064   G4String partName = particle->GetParticleName();
0065   G4double charge = particle->GetPDGCharge();
0066   G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
0067 
0068   // get material
0069   const G4Material* material = fDetector->GetMaterial();
0070   G4String matName = material->GetName();
0071   G4double density = material->GetDensity();
0072   G4double radl = material->GetRadlen();
0073 
0074   G4cout << "\n " << partName << " (" << G4BestUnit(energy, "Energy") << ") in "
0075          << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass")
0076          << ";   radiation length: " << G4BestUnit(radl, "Length") << ")" << G4endl;
0077 
0078   // get cuts
0079   GetCuts();
0080   if (charge != 0.) {
0081     G4cout << "\n  Range cuts: \t gamma " << std::setw(12) << G4BestUnit(fRangeCut[0], "Length")
0082            << "\t e- " << std::setw(12) << G4BestUnit(fRangeCut[1], "Length");
0083     G4cout << "\n Energy cuts: \t gamma " << std::setw(12) << G4BestUnit(fEnergyCut[0], "Energy")
0084            << "\t e- " << std::setw(12) << G4BestUnit(fEnergyCut[1], "Energy") << G4endl;
0085   }
0086 
0087   // max energy transfert
0088   if (charge != 0.) {
0089     G4double Mass_c2 = particle->GetPDGMass();
0090     G4double moverM = electron_mass_c2 / Mass_c2;
0091     G4double gamM1 = energy / Mass_c2, gam = gamM1 + 1., gamP1 = gam + 1.;
0092     G4double Tmax = energy;
0093     if (particle == G4Electron::Electron()) {
0094       Tmax *= 0.5;
0095     }
0096     else if (particle != G4Positron::Positron()) {
0097       Tmax = (2 * electron_mass_c2 * gamM1 * gamP1) / (1. + 2 * gam * moverM + moverM * moverM);
0098     }
0099     G4double range = emCal.GetCSDARange(Tmax, G4Electron::Electron(), material);
0100 
0101     G4cout << "\n  Max_energy _transferable  : " << G4BestUnit(Tmax, "Energy") << " ("
0102            << G4BestUnit(range, "Length") << ")" << G4endl;
0103   }
0104 
0105   // get processList and extract EM processes (but not MultipleScattering)
0106   G4ProcessVector* plist = particle->GetProcessManager()->GetProcessList();
0107   G4String procName;
0108   G4double cut;
0109   std::vector<G4String> emName;
0110   std::vector<G4double> enerCut;
0111   size_t length = plist->size();
0112   for (size_t j = 0; j < length; j++) {
0113     procName = (*plist)[j]->GetProcessName();
0114     cut = fEnergyCut[1];
0115     if ((procName == "eBrem") || (procName == "muBrems")) cut = fEnergyCut[0];
0116     if (((*plist)[j]->GetProcessType() == fElectromagnetic) && (procName != "msc")) {
0117       emName.push_back(procName);
0118       enerCut.push_back(cut);
0119     }
0120   }
0121 
0122   // write html documentation, if requested
0123   char* htmlDocName = std::getenv("G4PhysListName");  // file name
0124   char* htmlDocDir = std::getenv("G4PhysListDocDir");  // directory
0125   if (htmlDocName && htmlDocDir) {
0126     G4LossTableManager::Instance()->DumpHtml();
0127   }
0128 
0129   // print list of processes
0130   G4cout << "\n  processes :                ";
0131   for (size_t j = 0; j < emName.size(); ++j) {
0132     G4cout << "\t" << std::setw(14) << emName[j] << "\t";
0133   }
0134   G4cout << "\t" << std::setw(14) << "total";
0135 
0136   // compute cross section per atom (only for single material)
0137   if (material->GetNumberOfElements() == 1) {
0138     G4double Z = material->GetZ();
0139     G4double A = material->GetA();
0140 
0141     std::vector<G4double> sigma0;
0142     G4double sig, sigtot = 0.;
0143 
0144     for (size_t j = 0; j < emName.size(); j++) {
0145       sig = emCal.ComputeCrossSectionPerAtom(energy, particle, emName[j], Z, A, enerCut[j]);
0146       sigtot += sig;
0147       sigma0.push_back(sig);
0148     }
0149     sigma0.push_back(sigtot);
0150 
0151     G4cout << "\n \n  cross section per atom   : ";
0152     for (size_t j = 0; j < sigma0.size(); ++j) {
0153       G4cout << "\t" << std::setw(9) << G4BestUnit(sigma0[j], "Surface");
0154     }
0155     G4cout << G4endl;
0156   }
0157 
0158   // get cross section per volume
0159   std::vector<G4double> sigma0;
0160   std::vector<G4double> sigma1;
0161   std::vector<G4double> sigma2;
0162   G4double Sig, SigtotComp = 0., Sigtot = 0.;
0163 
0164   for (size_t j = 0; j < emName.size(); ++j) {
0165     Sig = emCal.ComputeCrossSectionPerVolume(energy, particle, emName[j], material, enerCut[j]);
0166     SigtotComp += Sig;
0167     sigma0.push_back(Sig);
0168     Sig = emCal.GetCrossSectionPerVolume(energy, particle, emName[j], material);
0169     Sigtot += Sig;
0170     sigma1.push_back(Sig);
0171     sigma2.push_back(Sig / density);
0172   }
0173   sigma0.push_back(SigtotComp);
0174   sigma1.push_back(Sigtot);
0175   sigma2.push_back(Sigtot / density);
0176 
0177   // print cross sections
0178   G4cout << "\n  compCrossSectionPerVolume: ";
0179   for (size_t j = 0; j < sigma0.size(); ++j) {
0180     G4cout << "\t" << std::setw(9) << sigma0[j] * cm << " cm^-1\t";
0181   }
0182   G4cout << "\n  cross section per volume : ";
0183   for (size_t j = 0; j < sigma1.size(); ++j) {
0184     G4cout << "\t" << std::setw(9) << sigma1[j] * cm << " cm^-1\t";
0185   }
0186 
0187   G4cout << "\n  cross section per mass   : ";
0188   for (size_t j = 0; j < sigma2.size(); ++j) {
0189     G4cout << "\t" << std::setw(9) << G4BestUnit(sigma2[j], "Surface/Mass");
0190   }
0191 
0192   // print mean free path
0193 
0194   G4double lambda;
0195 
0196   G4cout << "\n \n  mean free path           : ";
0197   for (size_t j = 0; j < sigma1.size(); ++j) {
0198     lambda = DBL_MAX;
0199     if (sigma1[j] > 0.) lambda = 1 / sigma1[j];
0200     G4cout << "\t" << std::setw(9) << G4BestUnit(lambda, "Length") << "   ";
0201   }
0202 
0203   // mean free path (g/cm2)
0204   G4cout << "\n        (g/cm2)            : ";
0205   for (size_t j = 0; j < sigma2.size(); ++j) {
0206     lambda = DBL_MAX;
0207     if (sigma2[j] > 0.) lambda = 1 / sigma2[j];
0208     G4cout << "\t" << std::setw(9) << G4BestUnit(lambda, "Mass/Surface");
0209   }
0210   G4cout << G4endl;
0211 
0212   if (charge == 0.) {
0213     G4cout.precision(prec);
0214     G4cout << "\n-----------------------------------------------------------\n" << G4endl;
0215     return;
0216   }
0217 
0218   // get stopping power
0219   std::vector<G4double> dedx1;
0220   std::vector<G4double> dedx2;
0221   G4double dedx, dedxtot = 0.;
0222   size_t nproc = emName.size();
0223 
0224   for (size_t j = 0; j < nproc; ++j) {
0225     dedx = emCal.ComputeDEDX(energy, particle, emName[j], material, enerCut[j]);
0226     dedxtot += dedx;
0227     dedx1.push_back(dedx);
0228     dedx2.push_back(dedx / density);
0229   }
0230   dedx1.push_back(dedxtot);
0231   dedx2.push_back(dedxtot / density);
0232 
0233   // print stopping power
0234   G4cout << "\n  restricted dE/dx         : ";
0235   for (size_t j = 0; j <= nproc; ++j) {
0236     G4cout << "\t" << std::setw(9) << G4BestUnit(dedx1[j], "Energy/Length");
0237   }
0238   G4cout << G4endl;
0239 
0240   G4cout << "\n      (MeV/g/cm2)          : ";
0241   for (size_t j = 0; j <= nproc; ++j) {
0242     G4cout << "\t" << std::setw(9) << G4BestUnit(dedx2[j], "Energy*Surface/Mass");
0243   }
0244   G4cout << G4endl;
0245 
0246   dedxtot = 0.;
0247 
0248   for (size_t j = 0; j < nproc; ++j) {
0249     dedx = emCal.ComputeDEDX(energy, particle, emName[j], material, energy);
0250     dedxtot += dedx;
0251     dedx1[j] = dedx;
0252     dedx2[j] = dedx / density;
0253   }
0254   dedx1[nproc] = dedxtot;
0255   dedx2[nproc] = dedxtot / density;
0256 
0257   // print stopping power
0258   G4cout << "\n  unrestricted dE/dx       : ";
0259   for (size_t j = 0; j <= nproc; ++j) {
0260     G4cout << "\t" << std::setw(9) << G4BestUnit(dedx1[j], "Energy/Length");
0261   }
0262   G4cout << G4endl;
0263 
0264   G4cout << "      (MeV/g/cm2)          : ";
0265   for (size_t j = 0; j <= nproc; ++j) {
0266     G4cout << "\t" << std::setw(9) << G4BestUnit(dedx2[j], "Energy*Surface/Mass");
0267   }
0268   G4cout << G4endl;
0269 
0270   // get range from restricted dedx
0271   G4double range1 = emCal.GetRangeFromRestricteDEDX(energy, particle, material);
0272   G4double range2 = range1 * density;
0273 
0274   // print range
0275   G4cout << "\n  range from restrict dE/dx: "
0276          << "\t" << std::setw(9) << G4BestUnit(range1, "Length") << " (" << std::setw(9)
0277          << G4BestUnit(range2, "Mass/Surface") << ")";
0278 
0279   // get range from full dedx
0280   G4double EmaxTable = G4EmParameters::Instance()->MaxEnergyForCSDARange();
0281   if (energy < EmaxTable) {
0282     G4double Range1 = emCal.GetCSDARange(energy, particle, material);
0283     G4double Range2 = Range1 * density;
0284 
0285     G4cout << "\n  range from full dE/dx    : "
0286            << "\t" << std::setw(9) << G4BestUnit(Range1, "Length") << " (" << std::setw(9)
0287            << G4BestUnit(Range2, "Mass/Surface") << ")";
0288   }
0289 
0290   // get transport mean free path (for multiple scattering)
0291   G4double MSmfp1 = emCal.GetMeanFreePath(energy, particle, "msc", material);
0292   G4double MSmfp2 = MSmfp1 * density;
0293 
0294   // print transport mean free path
0295   G4cout << "\n \n  transport mean free path : "
0296          << "\t" << std::setw(9) << G4BestUnit(MSmfp1, "Length") << " (" << std::setw(9)
0297          << G4BestUnit(MSmfp2, "Mass/Surface") << ")";
0298 
0299   if (particle == G4Electron::Electron()) CriticalEnergy();
0300 
0301   G4cout << "\n-------------------------------------------------------------\n";
0302   G4cout << G4endl;
0303 
0304   // reset default precision
0305   G4cout.precision(prec);
0306 }
0307 
0308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0309 
0310 void RunAction::EndOfRunAction(const G4Run*) {}
0311 
0312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0313 
0314 #include "G4ProductionCutsTable.hh"
0315 
0316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0317 
0318 void RunAction::GetCuts()
0319 {
0320   G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
0321 
0322   size_t numOfCouples = theCoupleTable->GetTableSize();
0323   const G4MaterialCutsCouple* couple = 0;
0324   G4int index = 0;
0325   for (size_t i = 0; i < numOfCouples; i++) {
0326     couple = theCoupleTable->GetMaterialCutsCouple(i);
0327     if (couple->GetMaterial() == fDetector->GetMaterial()) {
0328       index = i;
0329       break;
0330     }
0331   }
0332 
0333   fRangeCut[0] = (*(theCoupleTable->GetRangeCutsVector(idxG4GammaCut)))[index];
0334   fRangeCut[1] = (*(theCoupleTable->GetRangeCutsVector(idxG4ElectronCut)))[index];
0335   fRangeCut[2] = (*(theCoupleTable->GetRangeCutsVector(idxG4PositronCut)))[index];
0336 
0337   fEnergyCut[0] = (*(theCoupleTable->GetEnergyCutsVector(idxG4GammaCut)))[index];
0338   fEnergyCut[1] = (*(theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut)))[index];
0339   fEnergyCut[2] = (*(theCoupleTable->GetEnergyCutsVector(idxG4PositronCut)))[index];
0340 }
0341 
0342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0343 
0344 void RunAction::CriticalEnergy()
0345 {
0346   // compute e- critical energy (Rossi definition) and Moliere radius.
0347   // Review of Particle Physics - Eur. Phys. J. C3 (1998) page 147
0348   //
0349   G4EmCalculator emCal;
0350 
0351   const G4Material* material = fDetector->GetMaterial();
0352   const G4double radl = material->GetRadlen();
0353   G4double ekin = 5 * MeV;
0354   G4double deioni;
0355   G4double err = 1., errmax = 0.001;
0356   G4int iter = 0, itermax = 10;
0357   while (err > errmax && iter < itermax) {
0358     iter++;
0359     deioni = radl * emCal.ComputeDEDX(ekin, G4Electron::Electron(), "eIoni", material);
0360     err = std::abs(deioni - ekin) / ekin;
0361     ekin = deioni;
0362   }
0363   G4cout << "\n \n  critical energy (Rossi)  : "
0364          << "\t" << std::setw(8) << G4BestUnit(ekin, "Energy");
0365 
0366   // Pdg formula (only for single material)
0367   G4double pdga[2] = {610 * MeV, 710 * MeV};
0368   G4double pdgb[2] = {1.24, 0.92};
0369   G4double EcPdg = 0.;
0370 
0371   if (material->GetNumberOfElements() == 1) {
0372     G4int istat = 0;
0373     if (material->GetState() == kStateGas) istat = 1;
0374     G4double Zeff = material->GetZ() + pdgb[istat];
0375     EcPdg = pdga[istat] / Zeff;
0376     G4cout << "\t\t\t (from Pdg formula : " << std::setw(8) << G4BestUnit(EcPdg, "Energy") << ")";
0377   }
0378 
0379   const G4double Es = 21.2052 * MeV;
0380   G4double rMolier1 = Es / ekin, rMolier2 = rMolier1 * radl;
0381   G4cout << "\n  Moliere radius           : "
0382          << "\t" << std::setw(8) << rMolier1 << " X0 "
0383          << "= " << std::setw(8) << G4BestUnit(rMolier2, "Length");
0384 
0385   if (material->GetNumberOfElements() == 1) {
0386     G4double rMPdg = radl * Es / EcPdg;
0387     G4cout << "\t (from Pdg formula : " << std::setw(8) << G4BestUnit(rMPdg, "Length") << ")";
0388   }
0389 }
0390 
0391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!