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