TestEm12/src/Run.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/TestEm12/src/Run.cc
0027 /// \brief Implementation of the Run class
0028 //
0029 //
0030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032
0033 #include "Run.hh"
0034
0035 #include "DetectorConstruction.hh"
0036 #include "HistoManager.hh"
0037 #include "PrimaryGeneratorAction.hh"
0038
0039 #include "G4Material.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4UnitsTable.hh"
0042
0043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0044
0045 Run::Run(DetectorConstruction* detector) : fDetector(detector) {}
0046
0047 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0048
0049 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0050 {
0051   fParticle = particle;
0052   fEkin = energy;
0053 }
0054
0055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0056
0057 void Run::AddEdep(G4double e)
0058 {
0059   fEdeposit += e;
0060   fEdeposit2 += e * e;
0061 }
0062
0063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0064
0065 void Run::AddTrackLength(G4double t)
0066 {
0067   fTrackLen += t;
0068   fTrackLen2 += t * t;
0069 }
0070
0071 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0072
0073 void Run::AddProjRange(G4double x)
0074 {
0075   fProjRange += x;
0076   fProjRange2 += x * x;
0077 }
0078
0079 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0080
0081 void Run::AddStepSize(G4int nb, G4double st)
0082 {
0083   fNbOfSteps += nb;
0084   fNbOfSteps2 += nb * nb;
0085   fStepSize += st;
0086   fStepSize2 += st * st;
0087 }
0088
0089 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0090
0091 void Run::SetCsdaRange(G4double value)
0092 {
0093   fCsdaRange = value;
0094 }
0095
0096 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0097
0098 G4double Run::GetCsdaRange()
0099 {
0100   return fCsdaRange;
0101 }
0102
0103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0104
0105 void Run::Merge(const G4Run* run)
0106 {
0107   const Run* localRun = static_cast<const Run*>(run);
0108
0109   // pass information about primary particle
0110   fParticle = localRun->fParticle;
0111   fEkin = localRun->fEkin;
0112
0113   // accumulate sums
0114   fEdeposit += localRun->fEdeposit;
0115   fEdeposit2 += localRun->fEdeposit2;
0116   fTrackLen += localRun->fTrackLen;
0117   fTrackLen2 += localRun->fTrackLen2;
0118   fProjRange += localRun->fProjRange;
0119   fProjRange2 += localRun->fProjRange2;
0120   fNbOfSteps += localRun->fNbOfSteps;
0121   fNbOfSteps2 += localRun->fNbOfSteps2;
0122   fStepSize += localRun->fStepSize;
0123   fStepSize2 += localRun->fStepSize2;
0124
0125   fCsdaRange = localRun->fCsdaRange;
0126
0127   G4Run::Merge(run);
0128 }
0129
0130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0131
0132 void Run::EndOfRun()
0133 {
0134   std::ios::fmtflags mode = G4cout.flags();
0135   G4cout.setf(std::ios::fixed, std::ios::floatfield);
0136   G4int prec = G4cout.precision(2);
0137
0138   // run conditions
0139   //
0140   G4Material* material = fDetector->GetAbsorMaterial();
0141   G4double density = material->GetDensity();
0142   G4String partName = fParticle->GetParticleName();
0143
0144   G4cout << "\n ======================== run summary =====================\n";
0145   G4cout << "\n The run is " << numberOfEvent << " " << partName << " of "
0146          << G4BestUnit(fEkin, "Energy") << " through "
0147          << G4BestUnit(fDetector->GetAbsorRadius(), "Length") << " of " << material->GetName()
0148          << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
0149
0150   if (numberOfEvent == 0) {
0151     G4cout.setf(mode, std::ios::floatfield);
0152     G4cout.precision(prec);
0153     return;
0154   }
0155
0156   fEdeposit /= numberOfEvent;
0157   fEdeposit2 /= numberOfEvent;
0158   G4double rms = fEdeposit2 - fEdeposit * fEdeposit;
0159   if (rms > 0.)
0160     rms = std::sqrt(rms);
0161   else
0162     rms = 0.;
0163
0164   G4cout.precision(3);
0165   G4cout << "\n Total Energy deposited        = " << G4BestUnit(fEdeposit, "Energy") << " +- "
0166          << G4BestUnit(rms, "Energy") << G4endl;
0167
0168   // compute track length of primary track
0169   //
0170   fTrackLen /= numberOfEvent;
0171   fTrackLen2 /= numberOfEvent;
0172   rms = fTrackLen2 - fTrackLen * fTrackLen;
0173   if (rms > 0.)
0174     rms = std::sqrt(rms);
0175   else
0176     rms = 0.;
0177
0178   G4cout.precision(3);
0179   G4cout << "\n Track length of primary track = " << G4BestUnit(fTrackLen, "Length") << " +- "
0180          << G4BestUnit(rms, "Length");
0181
0182   // compare with csda range
0183   //
0184   G4cout << "\n Range from EmCalculator = " << G4BestUnit(fCsdaRange, "Length")
0185          << " (from full dE/dx)" << G4endl;
0186
0187   // compute projected range of primary track
0188   //
0189   fProjRange /= numberOfEvent;
0190   fProjRange2 /= numberOfEvent;
0191   rms = fProjRange2 - fProjRange * fProjRange;
0192   if (rms > 0.)
0193     rms = std::sqrt(rms);
0194   else
0195     rms = 0.;
0196
0197   G4cout << "\n Projected range               = " << G4BestUnit(fProjRange, "Length") << " +- "
0198          << G4BestUnit(rms, "Length") << G4endl;
0199
0200   // nb of steps and step size of primary track
0201   //
0202   G4double dNofEvents = double(numberOfEvent);
0203   G4double fNbSteps = fNbOfSteps / dNofEvents, fNbSteps2 = fNbOfSteps2 / dNofEvents;
0204   rms = fNbSteps2 - fNbSteps * fNbSteps;
0205   if (rms > 0.)
0206     rms = std::sqrt(rms);
0207   else
0208     rms = 0.;
0209
0210   G4cout.precision(2);
0211   G4cout << "\n Nb of steps of primary track  = " << fNbSteps << " +- " << rms;
0212
0213   fStepSize /= numberOfEvent;
0214   fStepSize2 /= numberOfEvent;
0215   rms = fStepSize2 - fStepSize * fStepSize;
0216   if (rms > 0.)
0217     rms = std::sqrt(rms);
0218   else
0219     rms = 0.;
0220
0221   G4cout.precision(3);
0222   G4cout << "\t Step size= " << G4BestUnit(fStepSize, "Length") << " +- "
0223          << G4BestUnit(rms, "Length") << G4endl;
0224
0225   // normalize histograms of longitudinal energy profile
0226   //
0227   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0228   G4int ih = 1;
0229   G4double binWidth = analysisManager->GetH1Width(ih) * analysisManager->GetH1Unit(ih);
0230   G4double fac = (1. / (numberOfEvent * binWidth)) * (mm / MeV);
0231   analysisManager->ScaleH1(ih, fac);
0232
0233   // normalize histogram d(E/E0)/d(r/r0)
0234   //
0235   ih = 8;
0236   binWidth = analysisManager->GetH1Width(ih);
0237   fac = 1. / (numberOfEvent * binWidth * fEkin);
0238   analysisManager->ScaleH1(ih, fac);
0239
0240   // reset default formats
0241   G4cout.setf(mode, std::ios::floatfield);
0242   G4cout.precision(prec);
0243 }
0244
0245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......