Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm6/src/RunAction.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 /// \file electromagnetic/TestEm6/src/RunAction.cc
 27 /// \brief Implementation of the RunAction class
 28 //
 29 //
 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32 
 33 #include "RunAction.hh"
 34 
 35 #include "G4AnnihiToMuPair.hh"
 36 #include "G4EmCalculator.hh"
 37 #include "G4MuonMinus.hh"
 38 #include "G4ParticleDefinition.hh"
 39 #include "G4ParticleTable.hh"
 40 #include "G4PhysicalConstants.hh"
 41 #include "G4Positron.hh"
 42 #include "G4Run.hh"
 43 #include "G4RunManager.hh"
 44 #include "G4SystemOfUnits.hh"
 45 #include "G4eBremsstrahlung.hh"
 46 #include "G4eeToHadrons.hh"
 47 #include "G4eeToHadronsModel.hh"
 48 #include "Randomize.hh"
 49 
 50 #include <sstream>
 51 
 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 53 
 54 RunAction::RunAction(DetectorConstruction* det)
 55   : G4UserRunAction(), fDetector(det), fProcCounter(0), fAnalysis(0), fMat(0)
 56 {
 57   fMinE = 40 * GeV;
 58   fMaxE = 10000 * GeV;
 59   fnBin = 10000;
 60 }
 61 
 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 63 
 64 RunAction::~RunAction() {}
 65 
 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 67 
 68 void RunAction::BeginOfRunAction(const G4Run* aRun)
 69 {
 70   G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl;
 71 
 72   // get material
 73   //
 74   fMat = fDetector->GetMaterial();
 75   G4cout << "###RunAction::BeginOfRunAction  Material:" << fMat->GetName() << G4endl;
 76 
 77   fProcCounter = new ProcessesCount;
 78 
 79   fAnalysis = G4AnalysisManager::Instance();
 80   fAnalysis->SetDefaultFileType("root");
 81 
 82   // Open an output file
 83   //
 84   std::stringstream tmp;
 85   tmp << "testem6_" << aRun->GetRunID();
 86   G4String fileName = tmp.str();
 87   fAnalysis->OpenFile(fileName);
 88   fAnalysis->SetVerboseLevel(2);
 89   fAnalysis->SetActivation(true);
 90 
 91   // Creating histograms 1,2,3,4,5,6
 92   //
 93   fAnalysis->SetFirstHistoId(1);
 94   fAnalysis->CreateH1("h1", "1/(1+(theta+[g]+)**2)", 100, 0, 1.);
 95   fAnalysis->CreateH1("h2", "log10(theta+ [g]+)", 100, -3., 1.);
 96   fAnalysis->CreateH1("h3", "log10(theta- [g]-)", 100, -3., 1.);
 97   fAnalysis->CreateH1("h4", "log10(theta+ [g]+ -theta- [g]-)", 100, -3., 1.);
 98   fAnalysis->CreateH1("h5", "xPlus", 100, 0., 1.);
 99   fAnalysis->CreateH1("h6", "xMinus", 100, 0., 1.);
100 
101   // creating histogram 7,8,9,10,11 (CrossSectionPerAtom)
102   //
103   G4double minBin = std::log10(fMinE / GeV);
104   G4double maxBin = std::log10(fMaxE / GeV);
105   fAnalysis->CreateH1("h7", "CrossSectionPerAtom of AnnihiToMuMu (microbarn)", fnBin, minBin,
106                       maxBin);
107   fAnalysis->CreateH1("h8", "CrossSectionPerAtom of AnnihiToTwoGamma (microbarn)", fnBin, minBin,
108                       maxBin);
109   fAnalysis->CreateH1("h9", "CrossSectionPerAtom of AnnihiToHadrons (microbarn)", fnBin, minBin,
110                       maxBin);
111   fAnalysis->CreateH1("h10", "Theoretical CrossSectionPerAtom of AnnihiToTwoGamma (microbarn)",
112                       fnBin, minBin, maxBin);
113   fAnalysis->CreateH1("h11", "Theoretical CrossSectionPerAtom of AnnihiToMuMu (microbarn)", fnBin,
114                       minBin, maxBin);
115 
116   // creating histogram 12,13,14,15,16(CrossSectionPerVolume)
117   //
118   fAnalysis->CreateH1("h12", "CrossSectionPerVol of Bremsstraulung (1/mm) ", fnBin, minBin, maxBin);
119   fAnalysis->CreateH1("h13", "CrossSectionPerVol of Ionization (1/mm)", fnBin, minBin, maxBin);
120   fAnalysis->CreateH1("h14", "CrossSectionPerVol of AnnihiToMuMu (1/mm)", fnBin, minBin, maxBin);
121   fAnalysis->CreateH1("h15", "CrossSectionPerVol of AnnihiToTwoGamma (1/mm)", fnBin, minBin,
122                       maxBin);
123   fAnalysis->CreateH1("h16", "CrossSectionPerVol of AnnihiToHadrons (1/mm)", fnBin, minBin, maxBin);
124 
125   // creating histogram 17 (R ratio)
126   fAnalysis->CreateH1("h17", "R : eeToHadr/eeToMu", fnBin, minBin, maxBin);
127 
128   G4cout << "\n----> Histogram file is opened in " << fileName << G4endl;
129 }
130 
131 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
132 
133 void RunAction::CountProcesses(G4String procName)
134 {
135   // does the process  already encounted ?
136   //
137   size_t nbProc = fProcCounter->size();
138   size_t i = 0;
139   while ((i < nbProc) && ((*fProcCounter)[i]->GetName() != procName))
140     i++;
141   if (i == nbProc) fProcCounter->push_back(new OneProcessCount(procName));
142 
143   (*fProcCounter)[i]->Count();
144 }
145 
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147 
148 void RunAction::EndOfRunAction(const G4Run*)
149 {
150   G4cout << "### RunAction::EndOfRunAction" << G4endl;
151   // total number of process calls
152   //
153   G4cout << "\n Number of process calls --->";
154   for (size_t i = 0; i < fProcCounter->size(); ++i) {
155     G4String procName = (*fProcCounter)[i]->GetName();
156     if (procName != "Transportation") {
157       G4int count = (*fProcCounter)[i]->GetCounter();
158       G4cout << "\t" << procName << " : " << count;
159     }
160   }
161   G4cout << G4endl;
162 
163   // instance EmCalculator
164   //
165   G4EmCalculator emCal;
166   emCal.SetVerbose(0);
167 
168   // get positron
169   //
170   G4String positronName = "e+";
171   G4ParticleDefinition* positron = G4ParticleTable::GetParticleTable()->FindParticle(positronName);
172 
173   // process name
174   //
175   G4String annihilName = "annihil";
176   G4String annihiToMuName = "AnnihiToMuPair";
177   G4String annihiToHadrName = "ee2hadr";
178   G4String BremName = "eBrem";
179   G4String IoniName = "eIoni";
180 
181   // get AnnihiToMuPair
182   //
183   G4AnnihiToMuPair* annihiToMu =
184     reinterpret_cast<G4AnnihiToMuPair*>(emCal.FindProcess(positron, annihiToMuName));
185 
186   // parameters for ComputeCrossSection
187   //
188   G4double atomicZ = 1.;
189   G4double atomicA = 2.;
190 
191   // set parameters for theory
192   //
193   const G4ParticleDefinition* muon = G4MuonMinus::MuonMinus();
194   G4double Mu = muon->GetPDGMass();
195   G4double Me = electron_mass_c2;
196   G4double Re = classic_electr_radius;
197   G4double Ru = Re * Me / Mu;
198   G4double Eth = 2 * Mu * Mu / Me - Me;
199   G4PhysicsLogVector v(fMinE, fMaxE, fnBin, false);
200 
201   // Compute CrossSections and Fill histgrams
202   //
203   for (G4int i = 0; i <= fnBin; ++i) {
204     G4double energy = v.Energy(i);
205     G4double x = std::log10(energy / GeV);
206 
207     // CrossSectionPerAtom
208     //
209     G4double crs_annihiToMu = annihiToMu->ComputeCrossSectionPerAtom(energy, atomicZ);
210     // G4cout << "crs_annihiToMu(mkb)=" << crs_annihiToMu/microbarn << G4endl;
211     fAnalysis->FillH1(7, x, crs_annihiToMu / microbarn);
212 
213     G4double crs_annihil =
214       emCal.ComputeCrossSectionPerAtom(energy, positron, annihilName, atomicZ, atomicA);
215     fAnalysis->FillH1(8, x, crs_annihil / microbarn);
216 
217     G4double crs_annihiToHadr =
218       emCal.ComputeCrossSectionPerAtom(energy, positron, annihiToHadrName, atomicZ, atomicA);
219     fAnalysis->FillH1(9, x, crs_annihiToHadr / microbarn);
220 
221     // CrossSectionPerVolume
222     //
223     G4double crsVol_Brem =
224       emCal.ComputeCrossSectionPerVolume(energy, positron, BremName, fMat, 100 * keV);
225     fAnalysis->FillH1(12, x, crsVol_Brem * mm);
226 
227     G4double crsVol_Ioni =
228       emCal.ComputeCrossSectionPerVolume(energy, positron, IoniName, fMat, 100 * keV);
229     fAnalysis->FillH1(13, x, crsVol_Ioni * mm);
230 
231     G4double crsVol_annihiToMu = annihiToMu->CrossSectionPerVolume(energy, fMat);
232     fAnalysis->FillH1(14, x, crsVol_annihiToMu * mm);
233 
234     G4double crsVol_annihil =
235       emCal.ComputeCrossSectionPerVolume(energy, positron, annihilName, fMat);
236     fAnalysis->FillH1(15, x, crsVol_annihil * mm);
237 
238     G4double crsVol_annihiToHadr =
239       emCal.ComputeCrossSectionPerVolume(energy, positron, annihiToHadrName, fMat);
240     fAnalysis->FillH1(16, x, crsVol_annihiToHadr * mm);
241 
242     // R ratio
243     //
244     G4double RR = 0.0;
245     if (crsVol_annihiToMu > 0.) RR = crsVol_annihiToHadr / crsVol_annihiToMu;
246     fAnalysis->FillH1(17, x, RR);
247 
248     // Theoretical calculation
249     //
250     G4double X1 = energy / Me;
251     if (X1 > 1 && i % 1000 == 0) {
252       G4double crs_annihil_theory =
253         atomicZ * pi * Re * Re
254         * ((X1 * X1 + 4 * X1 + 1) * G4Log(X1 + std::sqrt(X1 * X1 - 1)) / (X1 * X1 - 1)
255            - (X1 + 3) / std::sqrt(X1 * X1 - 1))
256         / (X1 + 1);
257       fAnalysis->FillH1(10, x, crs_annihil_theory / microbarn);
258     }
259 
260     G4double X2 = Eth / energy;
261     if (X2 < 1. && i % 1000 == 0) {
262       G4double crs_annihiToMu_theory =
263         atomicZ * pi * Ru * Ru / 3 * X2 * (1 + X2 / 2) * std::sqrt(1 - X2);
264       fAnalysis->FillH1(11, x, crs_annihiToMu_theory / microbarn);
265     }
266 
267     // if(i%1000==0)G4cout <<"###energy:" << energy << "/crs_ToMuMu:"
268     //         << crs_annihiToMu << "/crs_ToTwoGamma:"<< crs_annihil
269     //         <<"/crs_ToToHadr:"<<crs_annihiToHadr<< G4endl;
270   }
271 
272   fAnalysis->Write();
273   fAnalysis->CloseFile();
274   fAnalysis->Clear();
275 
276   G4cout << G4endl;
277 }
278 
279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
280