Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/exoticphysics/monopole/src/Run.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 exoticphysics/monopole/src/Run.cc
 27 /// \brief Implementation of the Run class
 28 //
 29 //
 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32 
 33 #include "Run.hh"
 34 
 35 #include "DetectorConstruction.hh"
 36 #include "PrimaryGeneratorAction.hh"
 37 
 38 #include "G4EmCalculator.hh"
 39 #include "G4Proton.hh"
 40 #include "G4SystemOfUnits.hh"
 41 #include "G4UnitsTable.hh"
 42 
 43 #include <iomanip>
 44 
 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 46 
 47 Run::Run(DetectorConstruction* det, PrimaryGeneratorAction* prim) : fDetector(det), fPrimary(prim)
 48 {
 49   fAnalysisManager = G4AnalysisManager::Instance();
 50 
 51   G4double length = fDetector->GetAbsorSizeX();
 52   fOffsetX = -0.5 * length;
 53 
 54   fVerboseLevel = 1;
 55   fNevt = 0;
 56   fProjRange = fProjRange2 = 0.;
 57 }
 58 
 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 60 
 61 Run::~Run() {}
 62 
 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 64 
 65 void Run::Merge(const G4Run* run)
 66 {
 67   const Run* localRun = static_cast<const Run*>(run);
 68 
 69   fNevt += localRun->GetNumberOfEvent();
 70   fProjRange += localRun->fProjRange;
 71   fProjRange2 += localRun->fProjRange2;
 72 
 73   G4Run::Merge(run);
 74 }
 75 
 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 77 
 78 void Run::EndOfRun(G4double binLength)
 79 {
 80   if (!G4Threading::IsMultithreadedApplication()) {
 81     fNevt += this->GetNumberOfEvent();
 82   }
 83 
 84   G4int nEvents = fNevt;
 85   if (nEvents == 0) {
 86     return;
 87   }
 88 
 89   // run conditions
 90   //
 91   const G4Material* material = fDetector->GetAbsorMaterial();
 92   G4double density = material->GetDensity();
 93   G4String matName = material->GetName();
 94 
 95   const G4ParticleDefinition* part = fPrimary->GetParticleGun()->GetParticleDefinition();
 96   G4String particle = part->GetParticleName();
 97   const G4ParticleDefinition* proton = G4Proton::Proton();
 98 
 99   G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
100 
101   if (GetVerbose() > 0) {
102     G4cout << "\n The run consists of " << nEvents << " " << particle << " of "
103            << G4BestUnit(energy, "Energy") << "\n through "
104            << G4BestUnit(fDetector->GetAbsorSizeX(), "Length") << " of " << matName
105            << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
106     // G4cout<<"Proj "<<fProjRange<<" "<<fProjRange2<<G4endl;
107   };
108 
109   // compute projected range and straggling
110   fProjRange /= nEvents;
111   fProjRange2 /= nEvents;
112   G4double rms = fProjRange2 - fProjRange * fProjRange;
113   if (rms > 0.) {
114     rms = std::sqrt(rms);
115   }
116   else {
117     rms = 0.;
118   }
119 
120   if (GetVerbose() > 0) {
121     G4cout.precision(5);
122     G4cout << " Projected Range= " << G4BestUnit(fProjRange, "Length")
123            << "   rms= " << G4BestUnit(rms, "Length") << "\n"
124            << G4endl;
125   };
126 
127   G4double ekin[100], dedxp[100], dedxmp[100], tdedxp[100], tdedxmp[100], xsp[100], xsmp[100];
128   G4EmCalculator calc;
129   // calc.SetVerbose(2);
130   G4int i;
131   for (i = 0; i < 100; ++i) {
132     ekin[i] = std::pow(10., 0.1 * G4double(i)) * keV;
133     dedxp[i] = calc.GetDEDX(ekin[i], proton, material);
134     xsp[i] = calc.GetCrossSectionPerVolume(ekin[i], proton, "hIoni", material);
135     tdedxp[i] = calc.ComputeElectronicDEDX(ekin[i], proton, material);
136     dedxmp[i] = calc.GetDEDX(ekin[i], part, material);
137     xsmp[i] = calc.GetCrossSectionPerVolume(ekin[i], part, "mplIoni", material);
138     tdedxmp[i] = calc.ComputeElectronicDEDX(ekin[i], part, material);
139   }
140 
141   if (GetVerbose() > 0) {
142     G4int prec = G4cout.precision(3);
143     G4cout << "##################################################################" << G4endl;
144     G4cout << "### Stopping Powers and Cross Sections" << G4endl;
145     G4cout << "##################################################################" << G4endl;
146 
147     G4cout << "# N   E(MeV)  p_dEdx(MeV/mm) mpl_dEdx(MeV/mm)           xs(1/mm)" << G4endl;
148     G4cout << "               restr    tot      restr    tot          p      mpl" << G4endl;
149     G4cout << "##################################################################" << G4endl;
150     for (i = 0; i < 100; ++i) {
151       G4cout << std::setw(2) << i << "." << std::setw(9) << ekin[i] << std::setw(8) << dedxp[i]
152              << std::setw(8) << tdedxp[i] << std::setw(9) << dedxmp[i] << std::setw(9) << tdedxmp[i]
153              << std::setw(10) << xsp[i] << std::setw(10) << xsmp[i] << G4endl;
154     }
155     G4cout.precision(prec);
156     G4cout << "##################################################################" << G4endl;
157   }
158 
159   // normalize histogram
160   G4double fac = (mm / MeV) / (nEvents * binLength);
161   fAnalysisManager->ScaleH1(1, fac);
162 
163   for (i = 0; i < 100; ++i) {
164     G4double e = std::log10(ekin[i] / MeV) + 0.05;
165     fAnalysisManager->FillH1(2, e, tdedxp[i]);
166     fAnalysisManager->FillH1(3, e, tdedxmp[i]);
167     fAnalysisManager->FillH1(4, e, std::log10(calc.GetRange(ekin[i], "proton", matName) / mm));
168     fAnalysisManager->FillH1(5, e, std::log10(calc.GetRange(ekin[i], "monopole", matName) / mm));
169     fAnalysisManager->FillH1(6, e, dedxp[i]);
170     fAnalysisManager->FillH1(7, e, dedxmp[i]);
171     fAnalysisManager->FillH1(8, e, xsp[i]);
172     fAnalysisManager->FillH1(9, e, xsmp[i]);
173   }
174 }
175 
176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
177 
178 void Run::FillHisto(G4int histoId, G4double v1, G4double v2)
179 {
180   fAnalysisManager->FillH1(histoId, v1, v2);
181 }
182 
183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
184