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