TestEm6/src/RunAction.cc

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/TestEm6/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 "G4AnnihiToMuPair.hh"
0036 #include "G4EmCalculator.hh"
0037 #include "G4MuonMinus.hh"
0038 #include "G4ParticleDefinition.hh"
0039 #include "G4ParticleTable.hh"
0040 #include "G4PhysicalConstants.hh"
0041 #include "G4Positron.hh"
0042 #include "G4Run.hh"
0043 #include "G4RunManager.hh"
0044 #include "G4SystemOfUnits.hh"
0045 #include "G4eBremsstrahlung.hh"
0046 #include "G4eeToHadrons.hh"
0047 #include "G4eeToHadronsModel.hh"
0048 #include "Randomize.hh"
0049
0050 #include <sstream>
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053
0054 RunAction::RunAction(DetectorConstruction* det)
0055   : G4UserRunAction(), fDetector(det), fProcCounter(0), fAnalysis(0), fMat(0)
0056 {
0057   fMinE = 40 * GeV;
0058   fMaxE = 10000 * GeV;
0059   fnBin = 10000;
0060 }
0061
0062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0063
0064 RunAction::~RunAction() {}
0065
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0067
0068 void RunAction::BeginOfRunAction(const G4Run* aRun)
0069 {
0070   G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
0071
0072   // get material
0073   //
0074   fMat = fDetector->GetMaterial();
0075   G4cout << "###RunAction::BeginOfRunAction  Material:" << fMat->GetName() << G4endl;
0076
0077   fProcCounter = new ProcessesCount;
0078
0079   fAnalysis = G4AnalysisManager::Instance();
0080   fAnalysis->SetDefaultFileType("root");
0081
0082   // Open an output file
0083   //
0084   std::stringstream tmp;
0085   tmp << "testem6_" << aRun->GetRunID();
0086   G4String fileName = tmp.str();
0087   fAnalysis->OpenFile(fileName);
0088   fAnalysis->SetVerboseLevel(2);
0089   fAnalysis->SetActivation(true);
0090
0091   // Creating histograms 1,2,3,4,5,6
0092   //
0093   fAnalysis->SetFirstHistoId(1);
0094   fAnalysis->CreateH1("h1", "1/(1+(theta+[g]+)**2)", 100, 0, 1.);
0095   fAnalysis->CreateH1("h2", "log10(theta+ [g]+)", 100, -3., 1.);
0096   fAnalysis->CreateH1("h3", "log10(theta- [g]-)", 100, -3., 1.);
0097   fAnalysis->CreateH1("h4", "log10(theta+ [g]+ -theta- [g]-)", 100, -3., 1.);
0098   fAnalysis->CreateH1("h5", "xPlus", 100, 0., 1.);
0099   fAnalysis->CreateH1("h6", "xMinus", 100, 0., 1.);
0100
0101   // creating histogram 7,8,9,10,11 (CrossSectionPerAtom)
0102   //
0103   G4double minBin = std::log10(fMinE / GeV);
0104   G4double maxBin = std::log10(fMaxE / GeV);
0105   fAnalysis->CreateH1("h7", "CrossSectionPerAtom of AnnihiToMuMu (microbarn)", fnBin, minBin,
0106                       maxBin);
0107   fAnalysis->CreateH1("h8", "CrossSectionPerAtom of AnnihiToTwoGamma (microbarn)", fnBin, minBin,
0108                       maxBin);
0109   fAnalysis->CreateH1("h9", "CrossSectionPerAtom of AnnihiToHadrons (microbarn)", fnBin, minBin,
0110                       maxBin);
0111   fAnalysis->CreateH1("h10", "Theoretical CrossSectionPerAtom of AnnihiToTwoGamma (microbarn)",
0112                       fnBin, minBin, maxBin);
0113   fAnalysis->CreateH1("h11", "Theoretical CrossSectionPerAtom of AnnihiToMuMu (microbarn)", fnBin,
0114                       minBin, maxBin);
0115
0116   // creating histogram 12,13,14,15,16(CrossSectionPerVolume)
0117   //
0118   fAnalysis->CreateH1("h12", "CrossSectionPerVol of Bremsstraulung (1/mm) ", fnBin, minBin, maxBin);
0119   fAnalysis->CreateH1("h13", "CrossSectionPerVol of Ionization (1/mm)", fnBin, minBin, maxBin);
0120   fAnalysis->CreateH1("h14", "CrossSectionPerVol of AnnihiToMuMu (1/mm)", fnBin, minBin, maxBin);
0121   fAnalysis->CreateH1("h15", "CrossSectionPerVol of AnnihiToTwoGamma (1/mm)", fnBin, minBin,
0122                       maxBin);
0123   fAnalysis->CreateH1("h16", "CrossSectionPerVol of AnnihiToHadrons (1/mm)", fnBin, minBin, maxBin);
0124
0125   // creating histogram 17 (R ratio)
0126   fAnalysis->CreateH1("h17", "R : eeToHadr/eeToMu", fnBin, minBin, maxBin);
0127
0128   G4cout << "\n----> Histogram file is opened in " << fileName << G4endl;
0129 }
0130
0131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0132
0133 void RunAction::CountProcesses(G4String procName)
0134 {
0135   // does the process  already encounted ?
0136   //
0137   size_t nbProc = fProcCounter->size();
0138   size_t i = 0;
0139   while ((i < nbProc) && ((*fProcCounter)[i]->GetName() != procName))
0140     i++;
0141   if (i == nbProc) fProcCounter->push_back(new OneProcessCount(procName));
0142
0143   (*fProcCounter)[i]->Count();
0144 }
0145
0146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0147
0148 void RunAction::EndOfRunAction(const G4Run*)
0149 {
0150   G4cout << "### RunAction::EndOfRunAction" << G4endl;
0151   // total number of process calls
0152   //
0153   G4cout << "\n Number of process calls --->";
0154   for (size_t i = 0; i < fProcCounter->size(); ++i) {
0155     G4String procName = (*fProcCounter)[i]->GetName();
0156     if (procName != "Transportation") {
0157       G4int count = (*fProcCounter)[i]->GetCounter();
0158       G4cout << "\t" << procName << " : " << count;
0159     }
0160   }
0161   G4cout << G4endl;
0162
0163   // instance EmCalculator
0164   //
0165   G4EmCalculator emCal;
0166   emCal.SetVerbose(0);
0167
0168   // get positron
0169   //
0170   G4String positronName = "e+";
0171   G4ParticleDefinition* positron = G4ParticleTable::GetParticleTable()->FindParticle(positronName);
0172
0173   // process name
0174   //
0175   G4String annihilName = "annihil";
0176   G4String annihiToMuName = "AnnihiToMuPair";
0177   G4String annihiToHadrName = "ee2hadr";
0178   G4String BremName = "eBrem";
0179   G4String IoniName = "eIoni";
0180
0181   // get AnnihiToMuPair
0182   //
0183   G4AnnihiToMuPair* annihiToMu =
0184     reinterpret_cast<G4AnnihiToMuPair*>(emCal.FindProcess(positron, annihiToMuName));
0185
0186   // parameters for ComputeCrossSection
0187   //
0188   G4double atomicZ = 1.;
0189   G4double atomicA = 2.;
0190
0191   // set parameters for theory
0192   //
0193   const G4ParticleDefinition* muon = G4MuonMinus::MuonMinus();
0194   G4double Mu = muon->GetPDGMass();
0195   G4double Me = electron_mass_c2;
0196   G4double Re = classic_electr_radius;
0197   G4double Ru = Re * Me / Mu;
0198   G4double Eth = 2 * Mu * Mu / Me - Me;
0199   G4PhysicsLogVector v(fMinE, fMaxE, fnBin, false);
0200
0201   // Compute CrossSections and Fill histgrams
0202   //
0203   for (G4int i = 0; i <= fnBin; ++i) {
0204     G4double energy = v.Energy(i);
0205     G4double x = std::log10(energy / GeV);
0206
0207     // CrossSectionPerAtom
0208     //
0209     G4double crs_annihiToMu = annihiToMu->ComputeCrossSectionPerAtom(energy, atomicZ);
0210     // G4cout << "crs_annihiToMu(mkb)=" << crs_annihiToMu/microbarn << G4endl;
0211     fAnalysis->FillH1(7, x, crs_annihiToMu / microbarn);
0212
0213     G4double crs_annihil =
0214       emCal.ComputeCrossSectionPerAtom(energy, positron, annihilName, atomicZ, atomicA);
0215     fAnalysis->FillH1(8, x, crs_annihil / microbarn);
0216
0217     G4double crs_annihiToHadr =
0218       emCal.ComputeCrossSectionPerAtom(energy, positron, annihiToHadrName, atomicZ, atomicA);
0219     fAnalysis->FillH1(9, x, crs_annihiToHadr / microbarn);
0220
0221     // CrossSectionPerVolume
0222     //
0223     G4double crsVol_Brem =
0224       emCal.ComputeCrossSectionPerVolume(energy, positron, BremName, fMat, 100 * keV);
0225     fAnalysis->FillH1(12, x, crsVol_Brem * mm);
0226
0227     G4double crsVol_Ioni =
0228       emCal.ComputeCrossSectionPerVolume(energy, positron, IoniName, fMat, 100 * keV);
0229     fAnalysis->FillH1(13, x, crsVol_Ioni * mm);
0230
0231     G4double crsVol_annihiToMu = annihiToMu->CrossSectionPerVolume(energy, fMat);
0232     fAnalysis->FillH1(14, x, crsVol_annihiToMu * mm);
0233
0234     G4double crsVol_annihil =
0235       emCal.ComputeCrossSectionPerVolume(energy, positron, annihilName, fMat);
0236     fAnalysis->FillH1(15, x, crsVol_annihil * mm);
0237
0238     G4double crsVol_annihiToHadr =
0239       emCal.ComputeCrossSectionPerVolume(energy, positron, annihiToHadrName, fMat);
0240     fAnalysis->FillH1(16, x, crsVol_annihiToHadr * mm);
0241
0242     // R ratio
0243     //
0244     G4double RR = 0.0;
0245     if (crsVol_annihiToMu > 0.) RR = crsVol_annihiToHadr / crsVol_annihiToMu;
0246     fAnalysis->FillH1(17, x, RR);
0247
0248     // Theoretical calculation
0249     //
0250     G4double X1 = energy / Me;
0251     if (X1 > 1 && i % 1000 == 0) {
0252       G4double crs_annihil_theory =
0253         atomicZ * pi * Re * Re
0254         * ((X1 * X1 + 4 * X1 + 1) * G4Log(X1 + std::sqrt(X1 * X1 - 1)) / (X1 * X1 - 1)
0255            - (X1 + 3) / std::sqrt(X1 * X1 - 1))
0256         / (X1 + 1);
0257       fAnalysis->FillH1(10, x, crs_annihil_theory / microbarn);
0258     }
0259
0260     G4double X2 = Eth / energy;
0261     if (X2 < 1. && i % 1000 == 0) {
0262       G4double crs_annihiToMu_theory =
0263         atomicZ * pi * Ru * Ru / 3 * X2 * (1 + X2 / 2) * std::sqrt(1 - X2);
0264       fAnalysis->FillH1(11, x, crs_annihiToMu_theory / microbarn);
0265     }
0266
0267     // if(i%1000==0)G4cout <<"###energy:" << energy << "/crs_ToMuMu:"
0268     //         << crs_annihiToMu << "/crs_ToTwoGamma:"<< crs_annihil
0269     //         <<"/crs_ToToHadr:"<<crs_annihiToHadr<< G4endl;
0270   }
0271
0272   fAnalysis->Write();
0273   fAnalysis->CloseFile();
0274   fAnalysis->Clear();
0275
0276   G4cout << G4endl;
0277 }
0278
0279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......