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