fanoCavity2/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 "G4Electron.hh"
0036 #include "G4EmCalculator.hh"
0037 #include "G4Run.hh"
0038 #include "G4RunManager.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4UnitsTable.hh"
0041 #include "Randomize.hh"
0042
0043 #include <iomanip>
0044
0045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0046
0047 Run::Run(DetectorConstruction* det, PrimaryGeneratorAction* kin, bool isMaster)
0048   : G4Run(),
0049     fDetector(det),
0050     fKinematic(kin),
0051     fProcCounter(0),
0052     fEdepCavity(0.),
0053     fEdepCavity2(0.),
0054     fTrkSegmCavity(0.),
0055     fNbEventCavity(0),
0056     fStepWall(0.),
0057     fStepWall2(0.),
0058     fStepCavity(0.),
0059     fStepCavity2(0.),
0060     fNbStepWall(0),
0061     fNbStepCavity(0),
0062     fEnergyGun(0.),
0063     fMassWall(0.),
0064     fMassCavity(0.),
0065     fIsMaster(isMaster)
0066 {
0067   // run conditions
0068   //
0069   G4ParticleDefinition* particleGun = fKinematic->GetParticleGun()->GetParticleDefinition();
0070   G4String partName = particleGun->GetParticleName();
0071   fEnergyGun = fKinematic->GetParticleGun()->GetParticleEnergy();
0072
0073   // geometry : effective wall volume
0074   //
0075   G4double cavityThickness = fDetector->GetCavityThickness();
0076   G4Material* mateCavity = fDetector->GetCavityMaterial();
0077   G4double densityCavity = mateCavity->GetDensity();
0078   fMassCavity = cavityThickness * densityCavity;
0079
0080   G4double wallThickness = fDetector->GetWallThickness();
0081   G4Material* mateWall = fDetector->GetWallMaterial();
0082   G4double densityWall = mateWall->GetDensity();
0083
0084   G4EmCalculator emCal;
0085   G4double RangeWall = emCal.GetCSDARange(fEnergyGun, particleGun, mateWall);
0086   G4double factor = 1.2;
0087   G4double effWallThick = factor * RangeWall;
0088   if ((effWallThick > wallThickness) || (effWallThick <= 0.)) effWallThick = wallThickness;
0089   fMassWall = 2 * effWallThick * densityWall;
0090
0091   G4double massTotal = fMassWall + fMassCavity;
0092   G4double fMassWallRatio = fMassWall / massTotal;
0093   fKinematic->RunInitialisation(effWallThick, fMassWallRatio);
0094
0095   G4double massRatio = fMassCavity / fMassWall;
0096
0097   // check radius
0098   //
0099   G4double worldRadius = fDetector->GetWallRadius();
0100   G4double RangeCavity = emCal.GetCSDARange(fEnergyGun, particleGun, mateCavity);
0101
0102   // G4String partName    = particleGun->GetParticleName();
0103
0104   std::ios::fmtflags mode = G4cout.flags();
0105   G4cout.setf(std::ios::fixed, std::ios::floatfield);
0106   G4int prec = G4cout.precision(3);
0107
0108   G4cout << "\n ===================== run conditions =====================\n";
0109
0110   G4cout << "\n The run will be " << numberOfEvent << " " << partName << " of "
0111          << G4BestUnit(fEnergyGun, "Energy") << " through 2*" << G4BestUnit(effWallThick, "Length")
0112          << " of " << mateWall->GetName()
0113          << " (density: " << G4BestUnit(densityWall, "Volumic Mass")
0114          << "); Mass/cm2 = " << G4BestUnit(fMassWall * cm2, "Mass")
0115          << "\n csdaRange: " << G4BestUnit(RangeWall, "Length") << G4endl;
0116
0117   G4cout << "\n the cavity is " << G4BestUnit(cavityThickness, "Length") << " of "
0118          << mateCavity->GetName() << " (density: " << G4BestUnit(densityCavity, "Volumic Mass")
0119          << "); Mass/cm2 = " << G4BestUnit(fMassCavity * cm2, "Mass")
0120          << " --> massRatio = " << std::setprecision(6) << massRatio << G4endl;
0121
0122   G4cout.precision(3);
0123   G4cout << " Wall radius: " << G4BestUnit(worldRadius, "Length")
0124          << "; range in cavity: " << G4BestUnit(RangeCavity, "Length") << G4endl;
0125
0126   G4cout << "\n ==========================================================\n";
0127
0128   // stopping power from EmCalculator
0129   //
0130   G4double dedxWall = emCal.GetDEDX(fEnergyGun, G4Electron::Electron(), mateWall);
0131   dedxWall /= densityWall;
0132   G4double dedxCavity = emCal.GetDEDX(fEnergyGun, G4Electron::Electron(), mateCavity);
0133   dedxCavity /= densityCavity;
0134
0135   G4cout << std::setprecision(4)
0136          << "\n StoppingPower in wall   = " << G4BestUnit(dedxWall, "Energy*Surface/Mass")
0137          << "\n               in cavity = " << G4BestUnit(dedxCavity, "Energy*Surface/Mass")
0138          << G4endl;
0139
0140   // process counter
0141   //
0142   fProcCounter = new ProcessesCount;
0143
0144   // charged particles and energy flow in cavity
0145   //
0146   fPartFlowCavity[0] = fPartFlowCavity[1] = 0;
0147   fEnerFlowCavity[0] = fEnerFlowCavity[1] = 0.;
0148
0149   // total energy deposit and charged track segment in cavity
0150   //
0151   fEdepCavity = fEdepCavity2 = fTrkSegmCavity = 0.;
0152   fNbEventCavity = 0;
0153
0154   // stepLenth of charged particles
0155   //
0156   fStepWall = fStepWall2 = fStepCavity = fStepCavity2 = 0.;
0157   fNbStepWall = fNbStepCavity = 0;
0158
0159   // reset default formats
0160   G4cout.setf(mode, std::ios::floatfield);
0161   G4cout.precision(prec);
0162
0163   // histograms
0164   //
0165   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0166   if (analysisManager->IsActive()) {
0167     analysisManager->OpenFile();
0168   }
0169 }
0170
0171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0172
0173 Run::~Run() {}
0174
0175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0176
0177 void Run::CountProcesses(G4String procName)
0178 {
0179   // does the process  already encounted ?
0180   size_t nbProc = fProcCounter->size();
0181   size_t i = 0;
0182   while ((i < nbProc) && ((*fProcCounter)[i]->GetName() != procName))
0183     i++;
0184   if (i == nbProc) fProcCounter->push_back(new OneProcessCount(procName));
0185
0186   (*fProcCounter)[i]->Count();
0187 }
0188
0189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0190
0191 void Run::SurveyConvergence(G4int NbofEvents)
0192 {
0193   if (NbofEvents == 0) return;
0194
0195   // beam fluence
0196   //
0197   G4int Nwall = fKinematic->GetWallCount();
0198   G4int Ncavity = fKinematic->GetCavityCount();
0199   G4double Iwall = Nwall / fMassWall;
0200   G4double Icavity = Ncavity / fMassCavity;
0201   G4double Iratio = Icavity / Iwall;
0202   G4double Itot = NbofEvents / (fMassWall + fMassCavity);
0203   G4double energyFluence = fEnergyGun * Itot;
0204
0205   // total dose in cavity
0206   //
0207   G4double doseCavity = fEdepCavity / fMassCavity;
0208   G4double ratio = doseCavity / energyFluence;
0209   G4double err = 100 * (ratio - 1.);
0210
0211   std::ios::fmtflags mode = G4cout.flags();
0212   G4cout.setf(std::ios::fixed, std::ios::floatfield);
0213   G4int prec = G4cout.precision(5);
0214
0215   G4cout << "--->evntNb= " << NbofEvents << " Nwall= " << Nwall << " Ncav= " << Ncavity
0216          << " Ic/Iw= " << Iratio << " Ne-_cav= " << fPartFlowCavity[0]
0217          << " doseCavity/Ebeam= " << ratio << "  (100*(ratio-1) = " << err << " %) \n"
0218          << G4endl;
0219
0220   // reset default formats
0221   G4cout.setf(mode, std::ios::floatfield);
0222   G4cout.precision(prec);
0223 }
0224
0225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0226
0227 void Run::EndOfRun()
0228 {  // Only call by Master thread
0229
0230   std::ios::fmtflags mode = G4cout.flags();
0231   G4cout.setf(std::ios::fixed, std::ios::floatfield);
0232   G4int prec = G4cout.precision(3);
0233
0234   if (numberOfEvent == 0) return;
0235
0236   // frequency of processes
0237   //
0238   G4cout << "\n Process calls frequency --->";
0239   for (size_t i = 0; i < fProcCounter->size(); i++) {
0240     G4String procName = (*fProcCounter)[i]->GetName();
0241     G4int count = (*fProcCounter)[i]->GetCounter();
0242     G4cout << "  " << procName << "= " << count;
0243   }
0244   G4cout << G4endl;
0245
0246   // charged particle flow in cavity
0247   //
0248   G4cout << "\n Charged particle flow in cavity :"
0249          << "\n      Enter --> nbParticles = " << fPartFlowCavity[0]
0250          << "\t Energy = " << G4BestUnit(fEnerFlowCavity[0], "Energy")
0251          << "\n      Exit  --> nbParticles = " << fPartFlowCavity[1]
0252          << "\t Energy = " << G4BestUnit(fEnerFlowCavity[1], "Energy") << G4endl;
0253
0254   if (fPartFlowCavity[0] == 0) return;
0255
0256   G4int Nwall = fKinematic->GetWallCount();
0257   G4int Ncavity = fKinematic->GetCavityCount();
0258
0259   G4double Iwall = Nwall / fMassWall;
0260   G4double Icavity = Ncavity / fMassCavity;
0261   G4double Iratio = Icavity / Iwall;
0262   G4double Itot = numberOfEvent / (fMassWall + fMassCavity);
0263   G4double energyFluence = fEnergyGun * Itot;
0264
0265   G4cout.precision(5);
0266   G4cout << "\n beamFluence in wall = " << Nwall << "\t in cavity = " << Ncavity
0267          << "\t Icav/Iwall = " << Iratio
0268          << "\t energyFluence = " << energyFluence / (MeV * cm2 / mg) << " MeV*cm2/mg" << G4endl;
0269
0270   // error on Edep in cavity
0271   //
0272   if (fNbEventCavity == 0) return;
0273   G4double meanEdep = fEdepCavity / fNbEventCavity;
0274   G4double meanEdep2 = fEdepCavity2 / fNbEventCavity;
0275   G4double varianceEdep = meanEdep2 - meanEdep * meanEdep;
0276   G4double dEoverE = 0.;
0277   if (varianceEdep > 0.) dEoverE = std::sqrt(varianceEdep / fNbEventCavity) / meanEdep;
0278
0279   // total dose in cavity
0280   //
0281   G4double doseCavity = fEdepCavity / fMassCavity;
0282
0283   G4double ratio = doseCavity / energyFluence, error = ratio * dEoverE;
0284
0285   G4cout << "\n Total edep in cavity = " << G4BestUnit(fEdepCavity, "Energy") << " +- "
0286          << 100 * dEoverE << " %"
0287          << "\n Total dose in cavity = " << doseCavity / (MeV * cm2 / mg) << " MeV*cm2/mg"
0288          << " +- " << 100 * dEoverE << " %"
0289          << "\n\n DoseCavity/EnergyFluence = " << ratio << " +- " << error << G4endl;
0290
0291   // track length in cavity
0292   G4double meantrack = fTrkSegmCavity / fPartFlowCavity[0];
0293
0294   G4cout.precision(4);
0295   G4cout << "\n Total charged trackLength in cavity = " << G4BestUnit(fTrkSegmCavity, "Length")
0296          << "   (mean value = " << G4BestUnit(meantrack, "Length") << ")" << G4endl;
0297
0298   // compute mean step size of charged particles
0299   //
0300   fStepWall /= fNbStepWall;
0301   fStepWall2 /= fNbStepWall;
0302   G4double rms = fStepWall2 - fStepWall * fStepWall;
0303   if (rms > 0.)
0304     rms = std::sqrt(rms);
0305   else
0306     rms = 0.;
0307   G4double nbTrackWall = fKinematic->GetWallCount();
0308
0309   G4cout << "\n StepSize of ch. tracks in wall   = " << G4BestUnit(fStepWall, "Length") << " +- "
0310          << G4BestUnit(rms, "Length") << "\t (nbSteps/track = " << double(fNbStepWall) / nbTrackWall
0311          << ")";
0312
0313   fStepCavity /= fNbStepCavity;
0314   fStepCavity2 /= fNbStepCavity;
0315   rms = fStepCavity2 - fStepCavity * fStepCavity;
0316   if (rms > 0.)
0317     rms = std::sqrt(rms);
0318   else
0319     rms = 0.;
0320
0321   G4cout << "\n StepSize of ch. tracks in cavity = " << G4BestUnit(fStepCavity, "Length") << " +- "
0322          << G4BestUnit(rms, "Length")
0323          << "\t (nbSteps/track = " << double(fNbStepCavity) / fPartFlowCavity[0] << ")";
0324
0325   G4cout << G4endl;
0326
0327   // reset default formats
0328   G4cout.setf(mode, std::ios::floatfield);
0329   G4cout.precision(prec);
0330
0331   // delete and remove all contents in fProcCounter
0332   while (fProcCounter->size() > 0) {
0333     OneProcessCount* aProcCount = fProcCounter->back();
0334     fProcCounter->pop_back();
0335     delete aProcCount;
0336   }
0337   delete fProcCounter;
0338
0339   // show Rndm status
0340   CLHEP::HepRandom::showEngineStatus();
0341 }
0342
0343 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0344
0345 void Run::Merge(const G4Run* run)
0346 {
0347   const Run* localRun = static_cast<const Run*>(run);
0348
0349   // Merge Run variables
0350   fPartFlowCavity[0] += localRun->fPartFlowCavity[0];
0351   fPartFlowCavity[1] += localRun->fPartFlowCavity[1];
0352   fEnerFlowCavity[0] += localRun->fEnerFlowCavity[0];
0353   fEnerFlowCavity[1] += localRun->fEnerFlowCavity[1];
0354   fEdepCavity += localRun->fEdepCavity;
0355   fEdepCavity2 += localRun->fEdepCavity2;
0356   fTrkSegmCavity += localRun->fTrkSegmCavity;
0357   fNbEventCavity += localRun->fNbEventCavity;
0358   fStepWall += localRun->fStepWall;
0359   fStepWall2 += localRun->fStepWall2;
0360   fStepCavity += localRun->fStepCavity;
0361   fStepCavity2 += localRun->fStepCavity2;
0362   fNbStepWall += localRun->fNbStepWall;
0363   fNbStepCavity += localRun->fNbStepCavity;
0364
0365   // Merge PrimaryGenerator variables
0366   fKinematic->AddWallCount(localRun->fKinematic->GetWallCount());
0367   fKinematic->AddCavityCount(localRun->fKinematic->GetCavityCount());
0368
0369   // Merge ProcessCount varaibles
0370   std::vector<OneProcessCount*>::iterator it;
0371   for (it = localRun->fProcCounter->begin(); it != localRun->fProcCounter->end(); it++) {
0372     OneProcessCount* process = *it;
0373     for (G4int i = 0; i < process->GetCounter(); i++)
0374       this->CountProcesses(process->GetName());
0375   }
0376
0377   G4Run::Merge(run);
0378 }
0379
0380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......