Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-02-23 09:22:42 

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 //
0030 
0031 #include "Run.hh"
0032 
0033 #include "DetectorConstruction.hh"
0034 #include "HistoManager.hh"
0035 #include "PrimaryGeneratorAction.hh"
0036 
0037 #include "G4ProcessTable.hh"
0038 #include "G4Radioactivation.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4TwoVector.hh"
0041 #include "G4UnitsTable.hh"
0042 
0043 #include <fstream>
0044 
0045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0046 
0047 Run::Run(DetectorConstruction* det) : fDetector(det) {}
0048 
0049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0050 
0051 void Run::SetPrimary(G4ParticleDefinition* particle, G4double energy)
0052 {
0053   fParticle = particle;
0054   fEkin = energy;
0055 }
0056 
0057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0058 
0059 void Run::CountProcesses(const G4VProcess* process, G4int iVol)
0060 {
0061   if (process == nullptr) return;
0062   G4String procName = process->GetProcessName();
0063   if (iVol == 1) {
0064     std::map<G4String, G4int>::iterator it1 = fProcCounter1.find(procName);
0065     if (it1 == fProcCounter1.end()) {
0066       fProcCounter1[procName] = 1;
0067     }
0068     else {
0069       fProcCounter1[procName]++;
0070     }
0071   }
0072 
0073   if (iVol == 2) {
0074     std::map<G4String, G4int>::iterator it2 = fProcCounter2.find(procName);
0075     if (it2 == fProcCounter2.end()) {
0076       fProcCounter2[procName] = 1;
0077     }
0078     else {
0079       fProcCounter2[procName]++;
0080     }
0081   }
0082 }
0083 
0084 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0085 
0086 void Run::ParticleCount(G4String name, G4double Ekin, G4int iVol)
0087 {
0088   if (iVol == 1) {
0089     std::map<G4String, ParticleData>::iterator it = fParticleDataMap1.find(name);
0090     if (it == fParticleDataMap1.end()) {
0091       fParticleDataMap1[name] = ParticleData(1, Ekin, Ekin, Ekin);
0092     }
0093     else {
0094       ParticleData& data = it->second;
0095       data.fCount++;
0096       data.fEmean += Ekin;
0097       // update min max
0098       G4double emin = data.fEmin;
0099       if (Ekin < emin) data.fEmin = Ekin;
0100       G4double emax = data.fEmax;
0101       if (Ekin > emax) data.fEmax = Ekin;
0102     }
0103   }
0104 
0105   if (iVol == 2) {
0106     std::map<G4String, ParticleData>::iterator it = fParticleDataMap2.find(name);
0107     if (it == fParticleDataMap2.end()) {
0108       fParticleDataMap2[name] = ParticleData(1, Ekin, Ekin, Ekin);
0109     }
0110     else {
0111       ParticleData& data = it->second;
0112       data.fCount++;
0113       data.fEmean += Ekin;
0114       // update min max
0115       G4double emin = data.fEmin;
0116       if (Ekin < emin) data.fEmin = Ekin;
0117       G4double emax = data.fEmax;
0118       if (Ekin > emax) data.fEmax = Ekin;
0119     }
0120   }
0121 }
0122 
0123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0124 
0125 void Run::AddEdep(G4double edep1, G4double edep2)
0126 {
0127   fEdepTarget += edep1;
0128   fEdepTarget2 += edep1 * edep1;
0129   fEdepDetect += edep2;
0130   fEdepDetect2 += edep2 * edep2;
0131 }
0132 
0133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0134 
0135 void Run::Merge(const G4Run* run)
0136 {
0137   const Run* localRun = static_cast<const Run*>(run);
0138 
0139   // primary particle info
0140   //
0141   fParticle = localRun->fParticle;
0142   fEkin = localRun->fEkin;
0143 
0144   // accumulate sums
0145   //
0146   fEdepTarget += localRun->fEdepTarget;
0147   fEdepTarget2 += localRun->fEdepTarget2;
0148   fEdepDetect += localRun->fEdepDetect;
0149   fEdepDetect2 += localRun->fEdepDetect2;
0150 
0151   // map: processes count in target
0152 
0153   std::map<G4String, G4int>::const_iterator itp1;
0154   for (itp1 = localRun->fProcCounter1.begin(); itp1 != localRun->fProcCounter1.end(); ++itp1) {
0155     G4String procName = itp1->first;
0156     G4int localCount = itp1->second;
0157     if (fProcCounter1.find(procName) == fProcCounter1.end()) {
0158       fProcCounter1[procName] = localCount;
0159     }
0160     else {
0161       fProcCounter1[procName] += localCount;
0162     }
0163   }
0164 
0165   // map: processes count in detector
0166 
0167   std::map<G4String, G4int>::const_iterator itp2;
0168   for (itp2 = localRun->fProcCounter2.begin(); itp2 != localRun->fProcCounter2.end(); ++itp2) {
0169     G4String procName = itp2->first;
0170     G4int localCount = itp2->second;
0171     if (fProcCounter2.find(procName) == fProcCounter2.end()) {
0172       fProcCounter2[procName] = localCount;
0173     }
0174     else {
0175       fProcCounter2[procName] += localCount;
0176     }
0177   }
0178 
0179   // map: created particles in target
0180   std::map<G4String, ParticleData>::const_iterator itc;
0181   for (itc = localRun->fParticleDataMap1.begin(); itc != localRun->fParticleDataMap1.end(); ++itc) {
0182     G4String name = itc->first;
0183     const ParticleData& localData = itc->second;
0184     if (fParticleDataMap1.find(name) == fParticleDataMap1.end()) {
0185       fParticleDataMap1[name] =
0186         ParticleData(localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax);
0187     }
0188     else {
0189       ParticleData& data = fParticleDataMap1[name];
0190       data.fCount += localData.fCount;
0191       data.fEmean += localData.fEmean;
0192       G4double emin = localData.fEmin;
0193       if (emin < data.fEmin) data.fEmin = emin;
0194       G4double emax = localData.fEmax;
0195       if (emax > data.fEmax) data.fEmax = emax;
0196     }
0197   }
0198 
0199   // map: created particle in detector
0200   std::map<G4String, ParticleData>::const_iterator itn;
0201   for (itn = localRun->fParticleDataMap2.begin(); itn != localRun->fParticleDataMap2.end(); ++itn) {
0202     G4String name = itn->first;
0203     const ParticleData& localData = itn->second;
0204     if (fParticleDataMap2.find(name) == fParticleDataMap2.end()) {
0205       fParticleDataMap2[name] =
0206         ParticleData(localData.fCount, localData.fEmean, localData.fEmin, localData.fEmax);
0207     }
0208     else {
0209       ParticleData& data = fParticleDataMap2[name];
0210       data.fCount += localData.fCount;
0211       data.fEmean += localData.fEmean;
0212       G4double emin = localData.fEmin;
0213       if (emin < data.fEmin) data.fEmin = emin;
0214       G4double emax = localData.fEmax;
0215       if (emax > data.fEmax) data.fEmax = emax;
0216     }
0217   }
0218 
0219   G4Run::Merge(run);
0220 }
0221 
0222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0223 
0224 void Run::EndOfRun()
0225 {
0226   G4int prec = 5, wid = prec + 2;
0227   G4int dfprec = G4cout.precision(prec);
0228 
0229   // run condition
0230   //
0231   G4String Particle = fParticle->GetParticleName();
0232   G4cout << "\n The run is " << numberOfEvent << " " << Particle << " of "
0233          << G4BestUnit(fEkin, "Energy") << " through : ";
0234 
0235   G4cout << "\n Target   : Length = " << G4BestUnit(fDetector->GetTargetLength(), "Length")
0236          << " Radius    = " << G4BestUnit(fDetector->GetTargetRadius(), "Length")
0237          << " Material = " << fDetector->GetTargetMaterial()->GetName();
0238   G4cout << "\n Detector : Length = " << G4BestUnit(fDetector->GetDetectorLength(), "Length")
0239          << " Thickness = " << G4BestUnit(fDetector->GetDetectorThickness(), "Length")
0240          << " Material = " << fDetector->GetDetectorMaterial()->GetName() << G4endl;
0241 
0242   if (numberOfEvent == 0) {
0243     G4cout.precision(dfprec);
0244     return;
0245   }
0246 
0247   // compute mean Energy deposited and rms in target
0248   //
0249   G4int TotNbofEvents = numberOfEvent;
0250   fEdepTarget /= TotNbofEvents;
0251   fEdepTarget2 /= TotNbofEvents;
0252   G4double rmsEdep = fEdepTarget2 - fEdepTarget * fEdepTarget;
0253   if (rmsEdep > 0.)
0254     rmsEdep = std::sqrt(rmsEdep);
0255   else
0256     rmsEdep = 0.;
0257 
0258   G4cout << "\n Mean energy deposit in target,   in time window = "
0259          << G4BestUnit(fEdepTarget, "Energy") << ";  rms = " << G4BestUnit(rmsEdep, "Energy")
0260          << G4endl;
0261 
0262   // compute mean Energy deposited and rms in detector
0263   //
0264   fEdepDetect /= TotNbofEvents;
0265   fEdepDetect2 /= TotNbofEvents;
0266   rmsEdep = fEdepDetect2 - fEdepDetect * fEdepDetect;
0267   if (rmsEdep > 0.)
0268     rmsEdep = std::sqrt(rmsEdep);
0269   else
0270     rmsEdep = 0.;
0271 
0272   G4cout << " Mean energy deposit in detector, in time window = "
0273          << G4BestUnit(fEdepDetect, "Energy") << ";  rms = " << G4BestUnit(rmsEdep, "Energy")
0274          << G4endl;
0275 
0276   // frequency of processes in target
0277   //
0278   G4cout << "\n Process calls frequency in target :" << G4endl;
0279   G4int index = 0;
0280   std::map<G4String, G4int>::iterator it1;
0281   for (it1 = fProcCounter1.begin(); it1 != fProcCounter1.end(); it1++) {
0282     G4String procName = it1->first;
0283     G4int count = it1->second;
0284     G4String space = " ";
0285     if (++index % 3 == 0) space = "\n";
0286     G4cout << " " << std::setw(20) << procName << "=" << std::setw(7) << count << space;
0287   }
0288   G4cout << G4endl;
0289 
0290   // frequency of processes in detector
0291   //
0292   G4cout << "\n Process calls frequency in detector:" << G4endl;
0293   index = 0;
0294   std::map<G4String, G4int>::iterator it2;
0295   for (it2 = fProcCounter2.begin(); it2 != fProcCounter2.end(); it2++) {
0296     G4String procName = it2->first;
0297     G4int count = it2->second;
0298     G4String space = " ";
0299     if (++index % 3 == 0) space = "\n";
0300     G4cout << " " << std::setw(20) << procName << "=" << std::setw(7) << count << space;
0301   }
0302   G4cout << G4endl;
0303 
0304   // particles count in target
0305   //
0306   G4cout << "\n List of generated particles in target:" << G4endl;
0307 
0308   std::map<G4String, ParticleData>::iterator itc;
0309   for (itc = fParticleDataMap1.begin(); itc != fParticleDataMap1.end(); itc++) {
0310     G4String name = itc->first;
0311     ParticleData data = itc->second;
0312     G4int count = data.fCount;
0313     G4double eMean = data.fEmean / count;
0314     G4double eMin = data.fEmin;
0315     G4double eMax = data.fEmax;
0316 
0317     G4cout << "  " << std::setw(13) << name << ": " << std::setw(7) << count
0318            << "  Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0319            << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")" << G4endl;
0320   }
0321 
0322   // particles count in detector
0323   //
0324   G4cout << "\n List of generated particles in detector:" << G4endl;
0325 
0326   std::map<G4String, ParticleData>::iterator itn;
0327   for (itn = fParticleDataMap2.begin(); itn != fParticleDataMap2.end(); itn++) {
0328     G4String name = itn->first;
0329     ParticleData data = itn->second;
0330     G4int count = data.fCount;
0331     G4double eMean = data.fEmean / count;
0332     G4double eMin = data.fEmin;
0333     G4double eMax = data.fEmax;
0334 
0335     G4cout << "  " << std::setw(13) << name << ": " << std::setw(7) << count
0336            << "  Emean = " << std::setw(wid) << G4BestUnit(eMean, "Energy") << "\t( "
0337            << G4BestUnit(eMin, "Energy") << " --> " << G4BestUnit(eMax, "Energy") << ")" << G4endl;
0338   }
0339   G4cout << G4endl;
0340 
0341   // activities in VR mode
0342   //
0343   WriteActivity(numberOfEvent);
0344 
0345   // remove all contents in fProcCounter, fCount
0346   fProcCounter1.clear();
0347   fProcCounter2.clear();
0348   fParticleDataMap1.clear();
0349   fParticleDataMap2.clear();
0350 
0351   // restore default format
0352   G4cout.precision(dfprec);
0353 }
0354 
0355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0356 
0357 void Run::WriteActivity(G4int nevent)
0358 {
0359   G4ProcessTable* pTable = G4ProcessTable::GetProcessTable();
0360   G4Radioactivation* rDecay =
0361     (G4Radioactivation*)pTable->FindProcess("Radioactivation", "GenericIon");
0362 
0363   // output the induced radioactivities (in VR mode only)
0364   //
0365   if ((rDecay == 0) || (rDecay->IsAnalogueMonteCarlo())) return;
0366 
0367   G4String fileName = G4AnalysisManager::Instance()->GetFileName() + ".activity";
0368   std::ofstream outfile(fileName, std::ios::out);
0369 
0370   std::vector<G4RadioactivityTable*> theTables = rDecay->GetTheRadioactivityTables();
0371 
0372   for (size_t i = 0; i < theTables.size(); i++) {
0373     G4double rate, error;
0374     outfile << "Radioactivities in decay window no. " << i << G4endl;
0375     outfile << "Z \tA \tE \tActivity (decays/window) \tError (decays/window) " << G4endl;
0376 
0377     map<G4ThreeVector, G4TwoVector>* aMap = theTables[i]->GetTheMap();
0378     map<G4ThreeVector, G4TwoVector>::iterator iter;
0379     for (iter = aMap->begin(); iter != aMap->end(); iter++) {
0380       rate = iter->second.x() / nevent;
0381       error = std::sqrt(iter->second.y()) / nevent;
0382       if (rate < 0.) rate = 0.;  // statically it can be < 0.
0383       outfile << iter->first.x() << "\t" << iter->first.y() << "\t" << iter->first.z() << "\t"
0384               << rate << "\t" << error << G4endl;
0385     }
0386     outfile << G4endl;
0387   }
0388   outfile.close();
0389 }
0390 
0391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!