Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-02-23 09:20:48 

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