monopole/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/monopole/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 "DetectorConstruction.hh"
0036 #include "PrimaryGeneratorAction.hh"
0037
0038 #include "G4EmCalculator.hh"
0039 #include "G4Proton.hh"
0040 #include "G4SystemOfUnits.hh"
0041 #include "G4UnitsTable.hh"
0042
0043 #include <iomanip>
0044
0045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0046
0047 Run::Run(DetectorConstruction* det, PrimaryGeneratorAction* prim) : fDetector(det), fPrimary(prim)
0048 {
0049   fAnalysisManager = G4AnalysisManager::Instance();
0050
0051   G4double length = fDetector->GetAbsorSizeX();
0052   fOffsetX = -0.5 * length;
0053
0054   fVerboseLevel = 1;
0055   fNevt = 0;
0056   fProjRange = fProjRange2 = 0.;
0057 }
0058
0059 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0060
0061 Run::~Run() {}
0062
0063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0064
0065 void Run::Merge(const G4Run* run)
0066 {
0067   const Run* localRun = static_cast<const Run*>(run);
0068
0069   fNevt += localRun->GetNumberOfEvent();
0070   fProjRange += localRun->fProjRange;
0071   fProjRange2 += localRun->fProjRange2;
0072
0073   G4Run::Merge(run);
0074 }
0075
0076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0077
0078 void Run::EndOfRun(G4double binLength)
0079 {
0080   if (!G4Threading::IsMultithreadedApplication()) {
0081     fNevt += this->GetNumberOfEvent();
0082   }
0083
0084   G4int nEvents = fNevt;
0085   if (nEvents == 0) {
0086     return;
0087   }
0088
0089   // run conditions
0090   //
0091   const G4Material* material = fDetector->GetAbsorMaterial();
0092   G4double density = material->GetDensity();
0093   G4String matName = material->GetName();
0094
0095   const G4ParticleDefinition* part = fPrimary->GetParticleGun()->GetParticleDefinition();
0096   G4String particle = part->GetParticleName();
0097   const G4ParticleDefinition* proton = G4Proton::Proton();
0098
0099   G4double energy = fPrimary->GetParticleGun()->GetParticleEnergy();
0100
0101   if (GetVerbose() > 0) {
0102     G4cout << "\n The run consists of " << nEvents << " " << particle << " of "
0103            << G4BestUnit(energy, "Energy") << "\n through "
0104            << G4BestUnit(fDetector->GetAbsorSizeX(), "Length") << " of " << matName
0105            << " (density: " << G4BestUnit(density, "Volumic Mass") << ")" << G4endl;
0106     // G4cout<<"Proj "<<fProjRange<<" "<<fProjRange2<<G4endl;
0107   };
0108
0109   // compute projected range and straggling
0110   fProjRange /= nEvents;
0111   fProjRange2 /= nEvents;
0112   G4double rms = fProjRange2 - fProjRange * fProjRange;
0113   if (rms > 0.) {
0114     rms = std::sqrt(rms);
0115   }
0116   else {
0117     rms = 0.;
0118   }
0119
0120   if (GetVerbose() > 0) {
0121     G4cout.precision(5);
0122     G4cout << " Projected Range= " << G4BestUnit(fProjRange, "Length")
0123            << "   rms= " << G4BestUnit(rms, "Length") << "\n"
0124            << G4endl;
0125   };
0126
0127   G4double ekin[100], dedxp[100], dedxmp[100], tdedxp[100], tdedxmp[100], xsp[100], xsmp[100];
0128   G4EmCalculator calc;
0129   // calc.SetVerbose(2);
0130   G4int i;
0131   for (i = 0; i < 100; ++i) {
0132     ekin[i] = std::pow(10., 0.1 * G4double(i)) * keV;
0133     dedxp[i] = calc.GetDEDX(ekin[i], proton, material);
0134     xsp[i] = calc.GetCrossSectionPerVolume(ekin[i], proton, "hIoni", material);
0135     tdedxp[i] = calc.ComputeElectronicDEDX(ekin[i], proton, material);
0136     dedxmp[i] = calc.GetDEDX(ekin[i], part, material);
0137     xsmp[i] = calc.GetCrossSectionPerVolume(ekin[i], part, "mplIoni", material);
0138     tdedxmp[i] = calc.ComputeElectronicDEDX(ekin[i], part, material);
0139   }
0140
0141   if (GetVerbose() > 0) {
0142     G4int prec = G4cout.precision(3);
0143     G4cout << "##################################################################" << G4endl;
0144     G4cout << "### Stopping Powers and Cross Sections" << G4endl;
0145     G4cout << "##################################################################" << G4endl;
0146
0147     G4cout << "# N   E(MeV)  p_dEdx(MeV/mm) mpl_dEdx(MeV/mm)           xs(1/mm)" << G4endl;
0148     G4cout << "               restr    tot      restr    tot          p      mpl" << G4endl;
0149     G4cout << "##################################################################" << G4endl;
0150     for (i = 0; i < 100; ++i) {
0151       G4cout << std::setw(2) << i << "." << std::setw(9) << ekin[i] << std::setw(8) << dedxp[i]
0152              << std::setw(8) << tdedxp[i] << std::setw(9) << dedxmp[i] << std::setw(9) << tdedxmp[i]
0153              << std::setw(10) << xsp[i] << std::setw(10) << xsmp[i] << G4endl;
0154     }
0155     G4cout.precision(prec);
0156     G4cout << "##################################################################" << G4endl;
0157   }
0158
0159   // normalize histogram
0160   G4double fac = (mm / MeV) / (nEvents * binLength);
0161   fAnalysisManager->ScaleH1(1, fac);
0162
0163   for (i = 0; i < 100; ++i) {
0164     G4double e = std::log10(ekin[i] / MeV) + 0.05;
0165     fAnalysisManager->FillH1(2, e, tdedxp[i]);
0166     fAnalysisManager->FillH1(3, e, tdedxmp[i]);
0167     fAnalysisManager->FillH1(4, e, std::log10(calc.GetRange(ekin[i], "proton", matName) / mm));
0168     fAnalysisManager->FillH1(5, e, std::log10(calc.GetRange(ekin[i], "monopole", matName) / mm));
0169     fAnalysisManager->FillH1(6, e, dedxp[i]);
0170     fAnalysisManager->FillH1(7, e, dedxmp[i]);
0171     fAnalysisManager->FillH1(8, e, xsp[i]);
0172     fAnalysisManager->FillH1(9, e, xsmp[i]);
0173   }
0174 }
0175
0176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0177
0178 void Run::FillHisto(G4int histoId, G4double v1, G4double v2)
0179 {
0180   fAnalysisManager->FillH1(histoId, v1, v2);
0181 }
0182
0183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......