amsEcal/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/TestEm11/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 #include "DetectorConstruction.hh"
0035 #include "HistoManager.hh"
0036
0037 #include "G4ParticleDefinition.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4UnitsTable.hh"
0040
0041 #include <iomanip>
0042
0043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0044
0045 Run::Run(DetectorConstruction* det)
0046 : G4Run(),
0047   fDetector(det),
0048   fParticle(0), fEkin(0.),
0049   nbOfModules(0), nbOfLayers(0), kLayerMax(0),
0050   EtotCalor(0.), Etot2Calor(0.), EvisCalor(0.), Evis2Calor(0.),
0051   Eleak(0.), Eleak2(0.)
0052 {
0053   nbOfModules = fDetector->GetNbModules();
0054   nbOfLayers  = fDetector->GetNbLayers();
0055   kLayerMax = nbOfModules*nbOfLayers + 1;
0056
0057   //initialize vectors
0058   //
0059   EtotLayer.resize(kLayerMax); Etot2Layer.resize(kLayerMax);
0060   EvisLayer.resize(kLayerMax); Evis2Layer.resize(kLayerMax);
0061   for (G4int k=0; k<kLayerMax; k++) {
0062     EtotLayer[k] = Etot2Layer[k] = EvisLayer[k] = Evis2Layer[k] = 0.0;
0063   }
0064
0065   EtotCalor = Etot2Calor = EvisCalor = Evis2Calor = Eleak = Eleak2 = 0.;
0066   EdLeak[0] = EdLeak[1] = EdLeak[2] = 0.;
0067 }
0068
0069 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0070
0071 Run::~Run()
0072 { }
0073
0074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0075
0076 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0077 {
0078   fParticle = particle;
0079   fEkin = energy;
0080 }
0081
0082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0083
0084 void Run::SumEvents_1(G4int layer, G4double Etot, G4double Evis)
0085 {
0086   //accumulate statistic per layer
0087   //
0088   EtotLayer[layer] += Etot;  Etot2Layer[layer] += Etot*Etot;
0089   EvisLayer[layer] += Evis;  Evis2Layer[layer] += Evis*Evis;
0090 }
0091
0092 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0093
0094 void Run::SumEvents_2(G4double etot, G4double evis, G4double eleak)
0095 {
0096   //accumulate statistic for full calorimeter
0097   //
0098   EtotCalor += etot;  Etot2Calor += etot*etot;
0099   EvisCalor += evis;  Evis2Calor += evis*evis;
0100   Eleak += eleak;  Eleak2 += eleak*eleak;
0101 }
0102
0103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0104
0105 void Run::DetailedLeakage(G4int icase, G4double energy)
0106 {
0107   //forward, backward, lateral leakage
0108   //
0109   EdLeak[icase] += energy;
0110 }
0111
0112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0113
0114 void Run::Merge(const G4Run* run)
0115 {
0116   const Run* localRun = static_cast<const Run*>(run);
0117
0118   // pass information about primary particle
0119   fParticle = localRun->fParticle;
0120   fEkin     = localRun->fEkin;
0121
0122   // accumulate sums
0123   //
0124   for (G4int k=0; k<kLayerMax; k++) {
0125     EtotLayer[k]  += localRun->EtotLayer[k];
0126     Etot2Layer[k] += localRun->Etot2Layer[k];
0127     EvisLayer[k]  += localRun->EvisLayer[k];
0128     Evis2Layer[k] += localRun->Evis2Layer[k];
0129   }
0130
0131   EtotCalor  += localRun->EtotCalor;
0132   Etot2Calor += localRun->Etot2Calor;
0133   EvisCalor  += localRun->EvisCalor;
0134   Evis2Calor += localRun->Evis2Calor;
0135   Eleak      += localRun->Eleak;
0136   Eleak2     += localRun->Eleak2;
0137   EdLeak[0]  += localRun->EdLeak[0];
0138   EdLeak[1]  += localRun->EdLeak[1];
0139   EdLeak[2]  += localRun->EdLeak[2];
0140
0141   G4Run::Merge(run);
0142 }
0143
0144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0145
0146 void Run::EndOfRun()
0147 {
0148   //calorimeter
0149   //
0150   fDetector->PrintCalorParameters();
0151
0152   //run conditions
0153   //
0154   G4String partName = fParticle->GetParticleName();
0155   G4int nbEvents = numberOfEvent;
0156
0157   G4int prec = G4cout.precision(3);
0158
0159   G4cout << " The run was " << nbEvents << " " << partName << " of "
0160          << G4BestUnit(fEkin,"Energy") << " through the calorimeter" << G4endl;
0161
0162   G4cout << "------------------------------------------------------------"
0163          << G4endl;
0164
0165   //if no events, return
0166   //
0167   if (nbEvents == 0) return;
0168
0169   //compute and print statistic
0170   //
0171   std::ios::fmtflags mode = G4cout.flags();
0172
0173   // energy in layers
0174   //
0175   G4cout.precision(prec);
0176   G4cout << "\n             "
0177          << "total Energy          (rms/mean)      "
0178          << "visible Energy        (rms/mean)" << G4endl;
0179
0180   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0181
0182   G4double meanEtot,meanEtot2,varianceEtot,rmsEtot,resEtot;
0183   G4double meanEvis,meanEvis2,varianceEvis,rmsEvis,resEvis;
0184
0185   for (G4int i1=1; i1<kLayerMax; i1++) {
0186     //total energy
0187     meanEtot  = EtotLayer[i1] /nbEvents;
0188     meanEtot2 = Etot2Layer[i1]/nbEvents;
0189     varianceEtot = meanEtot2 - meanEtot*meanEtot;
0190     resEtot = rmsEtot = 0.;
0191     if (varianceEtot > 0.) rmsEtot = std::sqrt(varianceEtot);
0192     if (meanEtot > 0.) resEtot = 100*rmsEtot/meanEtot;
0193     analysisManager->FillH1(3, i1+0.5, meanEtot);
0194
0195     //visible energy
0196     meanEvis  = EvisLayer[i1] /nbEvents;
0197     meanEvis2 = Evis2Layer[i1]/nbEvents;
0198     varianceEvis = meanEvis2 - meanEvis*meanEvis;
0199     resEvis = rmsEvis = 0.;
0200     if (varianceEvis > 0.) rmsEvis = std::sqrt(varianceEvis);
0201     if (meanEvis > 0.) resEvis = 100*rmsEvis/meanEvis;
0202     analysisManager->FillH1(4, i1+0.5, meanEvis);
0203
0204     //print
0205     //
0206     G4cout
0207       << "\n   layer " << i1 << ": "
0208       << std::setprecision(5)
0209       << std::setw(6) << G4BestUnit(meanEtot,"Energy") << " +- "
0210       << std::setprecision(4)
0211       << std::setw(5) << G4BestUnit( rmsEtot,"Energy") << "  ("
0212       << std::setprecision(2)
0213       << std::setw(3) << resEtot  << " %)"
0214       << "     "
0215       << std::setprecision(5)
0216       << std::setw(6) << G4BestUnit(meanEvis,"Energy") << " +- "
0217       << std::setprecision(4)
0218       << std::setw(5) << G4BestUnit( rmsEvis,"Energy") << "  ("
0219       << std::setprecision(2)
0220       << std::setw(3) << resEvis  << " %)";
0221   }
0222   G4cout << G4endl;
0223
0224   //calorimeter: total energy
0225   meanEtot  = EtotCalor /nbEvents;
0226   meanEtot2 = Etot2Calor/nbEvents;
0227   varianceEtot = meanEtot2 - meanEtot*meanEtot;
0228   resEtot = rmsEtot = 0.;
0229   if (varianceEtot > 0.) rmsEtot = std::sqrt(varianceEtot);
0230   if (meanEtot > 0.) resEtot = 100*rmsEtot/meanEtot;
0231
0232   //calorimeter: visible energy
0233   meanEvis  = EvisCalor /nbEvents;
0234   meanEvis2 = Evis2Calor/nbEvents;
0235   varianceEvis = meanEvis2 - meanEvis*meanEvis;
0236   resEvis = rmsEvis = 0.;
0237   if (varianceEvis > 0.) rmsEvis = std::sqrt(varianceEvis);
0238   if (meanEvis > 0.) resEvis = 100*rmsEvis/meanEvis;
0239
0240   //print
0241   //
0242   G4cout
0243     << "\n   total calor : "
0244     << std::setprecision(5)
0245     << std::setw(6) << G4BestUnit(meanEtot,"Energy") << " +- "
0246     << std::setprecision(4)
0247     << std::setw(5) << G4BestUnit( rmsEtot,"Energy") << "  ("
0248     << std::setprecision(2)
0249     << std::setw(3) << resEtot  << " %)"
0250     << "     "
0251     << std::setprecision(5)
0252     << std::setw(6) << G4BestUnit(meanEvis,"Energy") << " +- "
0253     << std::setprecision(4)
0254     << std::setw(5) << G4BestUnit( rmsEvis,"Energy") << "  ("
0255     << std::setprecision(2)
0256     << std::setw(3) << resEvis  << " %)";
0257
0258   G4cout << "\n------------------------------------------------------------"
0259          << G4endl;
0260
0261   //leakage
0262   G4double meanEleak,meanEleak2,varianceEleak,rmsEleak,ratio;
0263   meanEleak  = Eleak /nbEvents;
0264   meanEleak2 = Eleak2/nbEvents;
0265   varianceEleak = meanEleak2 - meanEleak*meanEleak;
0266   rmsEleak = 0.;
0267   if (varianceEleak > 0.) rmsEleak = std::sqrt(varianceEleak);
0268   ratio = 100*meanEleak/fEkin;
0269
0270   G4double forward = 100*EdLeak[0]/(nbEvents*fEkin);
0271   G4double bakward = 100*EdLeak[1]/(nbEvents*fEkin);
0272   G4double lateral = 100*EdLeak[2]/(nbEvents*fEkin);
0273
0274   //print
0275   //
0276   G4cout
0277     << "\n   Leakage : "
0278     << std::setprecision(5)
0279     << std::setw(6) << G4BestUnit(meanEleak,"Energy") << " +- "
0280     << std::setprecision(4)
0281     << std::setw(5) << G4BestUnit( rmsEleak,"Energy")
0282     << "\n   Eleak/Ebeam ="
0283     << std::setprecision(3)
0284     << std::setw(4) << ratio  << " %  ( forward ="
0285     << std::setw(4) << forward  << " %;   backward ="
0286     << std::setw(4) << bakward  << " %;   lateral ="
0287     << std::setw(4) << lateral  << " %)"
0288     << G4endl;
0289
0290   G4cout.setf(mode,std::ios::floatfield);
0291   G4cout.precision(prec);
0292
0293   //normalize histograms
0294   G4double factor = 1./nbEvents;
0295   analysisManager->ScaleH1(5,factor);
0296 }
0297
0298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......