dmparticle/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 exoticphysics/dmparticle/src/Run.cc
0027 /// \brief Implementation of the Run class
0028 //
0029
0030 #include "Run.hh"
0031
0032 #include "TestParameters.hh"
0033
0034 #include "G4ElectronIonPair.hh"
0035 #include "G4LossTableManager.hh"
0036 #include "G4PhysicalConstants.hh"
0037 #include "G4Run.hh"
0038 #include "G4Step.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "Randomize.hh"
0041
0042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0043
0044 Run::Run() : G4Run(), fParam(TestParameters::GetPointer())
0045 {
0046   fTotStepGas = fOverflow = fTotEdep = fStepGas = fMaxEnergy = 0.0;
0047   fEvt = fNbins = 0;
0048 }
0049
0050 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0051
0052 Run::~Run() {}
0053
0054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0055
0056 void Run::BeginOfRun()
0057 {
0058   // initilise scoring
0059   fTotStepGas = fOverflow = fTotEdep = fStepGas = fMaxEnergy = 0.0;
0060   fEvt = 0;
0061
0062   SetVerbose(1);
0063
0064   fNbins = fParam->GetNumberBins();
0065   fMaxEnergy = fParam->GetMaxEnergy();
0066
0067   fEgas.resize(fNbins, 0.0);
0068   fEdep.reset();
0069
0070   if (fVerbose > 0) {
0071     G4cout << "    BinsE= " << fNbins << "   Emax(keV)= " << fMaxEnergy / keV << G4endl;
0072   }
0073 }
0074
0075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0076
0077 void Run::EndOfRun()
0078 {
0079   G4int nEvt = GetNumberOfEvent();
0080
0081   G4double norm = (nEvt > 0) ? 1.0 / (G4double)nEvt : 0.0;
0082
0083   fTotStepGas *= norm;
0084   fOverflow *= norm;
0085
0086   G4double y1 = fEdep.mean();
0087   G4double y2 = fEdep.rms();
0088
0089   // G4double de = fMaxEnergy/G4double(fNbins);
0090   // G4double x1 = -de*0.5;
0091
0092   G4cout << " ====================================================" << G4endl;
0093   G4cout << "   Beam Particle: " << fParam->GetBeamParticle()->GetParticleName() << G4endl
0094          << "   Ekin(GeV)    = " << fParam->GetBeamEnergy() / GeV << G4endl
0095          << "   Z(mm)        = " << fParam->GetPositionZ() / mm << G4endl;
0096   G4cout << " ================== run summary =====================" << G4endl;
0097   G4int prec = G4cout.precision(5);
0098   G4cout << "   End of Run TotNbofEvents    = " << nEvt << G4endl;
0099
0100   G4cout << G4endl;
0101   G4cout << "   Mean energy deposit in absorber = " << y1 / GeV << " +- "
0102          << y2 * std::sqrt(norm) / GeV << " GeV; ";
0103   if (y1 > 0.0) {
0104     G4cout << "   RMS/Emean = " << y2 / y1;
0105   }
0106   G4cout << G4endl;
0107   G4cout << "   Mean number of steps in absorber= " << fTotStepGas << G4endl;
0108   G4cout << G4endl;
0109   /*
0110   G4cout << " ====== Energy deposit distribution   Noverflows= " << fOverflow
0111          << " ====== " << G4endl ;
0112   G4cout << " bin nb      Elow      entries     normalized " << G4endl;
0113
0114   x1 = 0.0;
0115
0116   for(G4int j=0; j<fNbins; ++j)
0117   {
0118     G4cout << std::setw(5) << j << std::setw(10) << x1/keV
0119            << std::setw(12) << fEgas[j] << std::setw(12) << fEgas[j]*norm
0120            << G4endl ;
0121     x1 += de;
0122   }
0123   */
0124   G4cout.precision(prec);
0125
0126   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0127
0128   // normalize histograms
0129
0130   analysisManager->ScaleH1(1, norm);
0131
0132   G4cout << " ================== run end ==========================" << G4endl;
0133 }
0134
0135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0136
0137 void Run::BeginOfEvent()
0138 {
0139   fTotEdep = 0.0;
0140   fStepGas = 0;
0141   ++fEvt;
0142 }
0143
0144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0145
0146 void Run::EndOfEvent()
0147 {
0148   /*
0149   G4cout << "fStepGas = " << fStepGas << " fTotStepGas= " << fTotStepGas
0150          << " fTotEdep= " << fTotEdep << "  fNbins= " << fNbins
0151          << " fMaxEnergy= " << fMaxEnergy << G4endl;
0152   */
0153   fTotStepGas += fStepGas;
0154
0155   G4int idx = G4lrint(fTotEdep * fNbins / fMaxEnergy);
0156
0157   if (idx < 0) {
0158     fEgas[0] += 1.0;
0159   }
0160   if (idx >= fNbins) {
0161     fOverflow += 1.0;
0162   }
0163   else {
0164     fEgas[idx] += 1.0;
0165   }
0166
0167   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0168
0169   // fill histo
0170
0171   analysisManager->FillH1(1, fTotEdep / GeV, 1.0);
0172   fEdep.fill(fTotEdep, 1.0);
0173 }
0174
0175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0176 //
0177 // Allows one to accumulate data from different threads
0178
0179 void Run::Merge(const G4Run* run)
0180 {
0181   const Run* localRun = static_cast<const Run*>(run);
0182
0183   fTotStepGas += localRun->fTotStepGas;
0184   fOverflow += localRun->fOverflow;
0185
0186   G4StatDouble* stat = const_cast<G4StatDouble*>(localRun->GetStat());
0187
0188   fEdep.add(stat);
0189
0190   for (G4int j = 0; j < fNbins; ++j) {
0191     fEgas[j] += localRun->fEgas[j];
0192   }
0193
0194   G4Run::Merge(run);
0195 }
0196
0197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0198 //
0199 // Is called from TargetSD
0200
0201 void Run::AddEnergy(G4double edep, const G4Step* step)
0202 {
0203   if (1 < fVerbose) {
0204     G4cout << "Run::AddEnergy: e(keV)= " << edep / keV << G4endl;
0205   }
0206   fTotEdep += edep;
0207
0208   if (step) {
0209     if (1 == step->GetTrack()->GetTrackID()) {
0210       fStepGas += 1.0;
0211     }
0212   }
0213 }
0214
0215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......