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