electronScattering/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 "HistoManager.hh"
0033 #include "PrimaryGeneratorAction.hh"
0034
0035 #include "G4PhysicalConstants.hh"
0036 #include "G4Run.hh"
0037 #include "G4RunManager.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4UnitsTable.hh"
0040 #include "Randomize.hh"
0041
0042 #include <iomanip>
0043 #include <stdexcept>  // std::out_of_range
0044
0045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0046
0047 Run::Run(DetectorConstruction* det) : G4Run(), fDetector(det), fParticle(0), fEnergy(0)
0048 {
0049   InitFluence();
0050 }
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053
0054 Run::~Run() {}
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057
0058 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0059 {
0060   fParticle = particle;
0061   fEnergy = energy;
0062 }
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0065
0066 void Run::Merge(const G4Run* run)
0067 {
0068   const Run* localRun = static_cast<const Run*>(run);
0069
0070   // pass information about primary particle
0071   fParticle = localRun->fParticle;
0072   fEnergy = localRun->fEnergy;
0073
0074   // In MT all threads have the same fNbBins and fDr value
0075   fNbBins = localRun->fNbBins;
0076   fDr = localRun->fDr;
0077
0078   // Some value by value
0079
0080   for (unsigned int i = 0; i < localRun->fluence.size(); ++i) {
0081     try {
0082       fluence[i] += localRun->fluence[i];
0083     }
0084     catch (const std::out_of_range&) {
0085       fluence.push_back(localRun->fluence[i]);
0086     }
0087   }
0088
0089   for (unsigned int i = 0; i < localRun->fluence1.size(); ++i) {
0090     try {
0091       fluence1[i] += localRun->fluence1[i];
0092     }
0093     catch (const std::out_of_range&) {
0094       fluence1.push_back(localRun->fluence1[i]);
0095     }
0096   }
0097
0098   for (unsigned int i = 0; i < localRun->fluence2.size(); ++i) {
0099     try {
0100       fluence2[i] += localRun->fluence2[i];
0101     }
0102     catch (const std::out_of_range&) {
0103       fluence2.push_back(localRun->fluence2[i]);
0104     }
0105   }
0106
0107   for (unsigned int i = 0; i < localRun->fNbEntries.size(); ++i) {
0108     try {
0109       fNbEntries[i] += localRun->fNbEntries[i];
0110     }
0111     catch (const std::out_of_range&) {
0112       fNbEntries.push_back(localRun->fNbEntries[i]);
0113     }
0114   }
0115
0116   G4Run::Merge(run);
0117 }
0118
0119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0120
0121 void Run::EndOfRun()
0122 {
0123   if (numberOfEvent == 0) return;
0124
0125   // Scatter foil
0126
0127   G4Material* material = fDetector->GetMaterialScatter();
0128   G4double length = fDetector->GetThicknessScatter();
0129   G4double density = material->GetDensity();
0130   G4String partName = fParticle->GetParticleName();
0131
0132   G4cout << "\n ======================== run summary ======================\n";
0133
0134   G4int prec = G4cout.precision(3);
0135
0136   G4cout << "\n The run was " << numberOfEvent << " " << partName << " of "
0137          << G4BestUnit(fEnergy, "Energy") << " through " << G4BestUnit(length, "Length") << " of "
0138          << material->GetName() << " (density: " << G4BestUnit(density, "Volumic Mass") << ")"
0139          << G4endl;
0140
0141   G4cout.precision(prec);
0142
0143   // normalize histograms
0144   //
0145   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0146   G4double fac = 1. / (double(numberOfEvent));
0147   analysisManager->ScaleH1(1, fac);
0148   analysisManager->ScaleH1(2, fac);
0149   analysisManager->ScaleH1(3, fac);
0150   analysisManager->ScaleH1(5, fac);
0151   analysisManager->ScaleH1(6, fac);
0152
0153   ComputeFluenceError();
0154   PrintFluence(numberOfEvent);
0155 }
0156
0157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0158
0159 void Run::InitFluence()
0160 {
0161   // create Fluence histo in any case
0162
0163   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0164   G4int ih = 4;
0165   analysisManager->SetH1(ih, 120, 0 * mm, 240 * mm, "mm");
0166
0167   // construct vectors for fluence distribution
0168
0169   fNbBins = analysisManager->GetH1Nbins(ih);
0170   fDr = analysisManager->GetH1Width(ih);
0171   fluence.resize(fNbBins, 0.);
0172   fluence1.resize(fNbBins, 0.);
0173   fluence2.resize(fNbBins, 0.);
0174   fNbEntries.resize(fNbBins, 0);
0175 }
0176
0177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0178
0179 void Run::SumFluence(G4double r, G4double fl)
0180 {
0181   G4int ibin = (int)(r / fDr);
0182   if (ibin >= fNbBins) return;
0183   fNbEntries[ibin]++;
0184   fluence[ibin] += fl;
0185   fluence2[ibin] += fl * fl;
0186 }
0187
0188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0189
0190 void Run::ComputeFluenceError()
0191 {
0192   // compute rms
0193
0194   G4double ds, variance, rms;
0195   G4double rmean = -0.5 * fDr;
0196
0197   for (G4int bin = 0; bin < fNbBins; bin++) {
0198     rmean += fDr;
0199     ds = twopi * rmean * fDr;
0200     fluence[bin] /= ds;
0201     fluence2[bin] /= (ds * ds);
0202     variance = 0.;
0203     if (fNbEntries[bin] > 0)
0204       variance = fluence2[bin] - (fluence[bin] * fluence[bin]) / fNbEntries[bin];
0205     rms = 0.;
0206     if (variance > 0.) rms = std::sqrt(variance);
0207     fluence2[bin] = rms;
0208   }
0209
0210   // normalize to first bins, compute error and fill histo
0211
0212   G4double rnorm(4 * mm), radius(0.), fnorm(0.), fnorm2(0.);
0213   G4int inorm = -1;
0214   do {
0215     inorm++;
0216     radius += fDr;
0217     fnorm += fluence[inorm];
0218     fnorm2 += fluence2[inorm];
0219   } while (radius < rnorm);
0220   fnorm /= (inorm + 1);
0221   fnorm2 /= (inorm + 1);
0222
0223   G4double ratio, error;
0224   G4double scale = 1. / fnorm;
0225   G4double err0 = fnorm2 / fnorm, err1 = 0.;
0226
0227   rmean = -0.5 * fDr;
0228
0229   for (G4int bin = 0; bin < fNbBins; bin++) {
0230     ratio = fluence[bin] * scale;
0231     error = 0.;
0232     if (ratio > 0.) {
0233       err1 = fluence2[bin] / fluence[bin];
0234       error = ratio * std::sqrt(err1 * err1 + err0 * err0);
0235     }
0236     fluence1[bin] = ratio;
0237     fluence2[bin] = error;
0238     rmean += fDr;
0239     G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0240     analysisManager->FillH1(4, rmean, ratio);
0241   }
0242 }
0243
0244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0245
0246 #include <fstream>
0247
0248 void Run::PrintFluence(G4int TotEvents)
0249 {
0250   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0251   G4String name = analysisManager->GetFileName();
0252   G4String fileName = name + ".ascii";
0253   std::ofstream File(fileName, std::ios::out);
0254
0255   std::ios::fmtflags mode = File.flags();
0256   File.setf(std::ios::scientific, std::ios::floatfield);
0257   G4int prec = File.precision(3);
0258
0259   File << "  Fluence density distribution \n "
0260        << "\n  ibin \t radius (mm) \t Nb \t fluence\t norma fl\t rms/nfl (%) \n"
0261        << G4endl;
0262
0263   G4double rmean = -0.5 * fDr;
0264   for (G4int bin = 0; bin < fNbBins; bin++) {
0265     rmean += fDr;
0266     G4double error = 0.;
0267     if (fluence1[bin] > 0.) error = 100 * fluence2[bin] / fluence1[bin];
0268     File << "  " << bin << "\t " << rmean / mm << "\t " << fNbEntries[bin] << "\t "
0269          << fluence[bin] / double(TotEvents) << "\t " << fluence1[bin] << "\t " << error << G4endl;
0270   }
0271
0272   // restaure default formats
0273   File.setf(mode, std::ios::floatfield);
0274   File.precision(prec);
0275 }
0276
0277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......