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