Back to home page

EIC code displayed by LXR

 
 

    


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

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 #include "Scorer.hh"
0027 
0028 #include "PrimaryGeneratorAction.hh"
0029 #include "TimeStepAction.hh"
0030 
0031 #include "G4AnalysisManager.hh"
0032 #include "G4DNAEventScheduler.hh"
0033 #include "G4DNAScavengerMaterial.hh"
0034 #include "G4Event.hh"
0035 #include "G4MoleculeTable.hh"
0036 #include "G4PhysicalConstants.hh"
0037 #include "G4RunManager.hh"
0038 #include "G4Scheduler.hh"
0039 #include "G4UIcmdWithADoubleAndUnit.hh"
0040 #include "G4UIcmdWithAnInteger.hh"
0041 #include "G4VChemistryWorld.hh"
0042 
0043 #include <G4EventManager.hh>
0044 #include <G4MolecularConfiguration.hh>
0045 #include <G4MoleculeCounter.hh>
0046 #include <G4SystemOfUnits.hh>
0047 #include <globals.hh>
0048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0049 
0050 // Dose
0051 Dose::Dose()
0052   : G4UImessenger(),
0053     fpDoseDir(new G4UIdirectory("/scorer/Dose/")),
0054     fpAddDoseCutOff(new G4UIcmdWithADoubleAndUnit("/scorer/Dose/cutoff", this)),
0055     fpAddDoseToAbort(new G4UIcmdWithADoubleAndUnit("/scorer/Dose/abortedDose", this))
0056 {
0057   fpDoseDir->SetGuidance("Dose scorer commands");
0058 }
0059 
0060 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0061 
0062 void Dose::SetNewValue(G4UIcommand* command, G4String newValue)
0063 {
0064   if (command == fpAddDoseCutOff.get()) {
0065     fDosesCutOff = G4UIcmdWithADoubleAndUnit::GetNewDoubleValue(newValue);
0066   }
0067   if (command == fpAddDoseToAbort.get()) {
0068     fDosesToAbort = G4UIcmdWithADoubleAndUnit::GetNewDoubleValue(newValue);
0069   }
0070 }
0071 
0072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0073 
0074 template<>
0075 void Scorer<Dose>::SetChemistryWorld(G4VChemistryWorld* pChemistryWorld)
0076 {
0077   fpChemistryWorld = pChemistryWorld;
0078 }
0079 
0080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0081 
0082 template<>
0083 G4VChemistryWorld* Scorer<Dose>::GetChemistryWorld() const
0084 {
0085   return fpChemistryWorld;
0086 }
0087 
0088 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0089 
0090 template<>
0091 void Scorer<Dose>::clear()
0092 {
0093   fpScorer->fCumulatedDose = 0.;
0094 }
0095 
0096 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0097 
0098 template<>
0099 void Scorer<Dose>::Initialize(G4HCofThisEvent* HCE)
0100 {
0101   clear();
0102   fpEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
0103   if (fHCID < 0) {
0104     fHCID = GetCollectionID(0);
0105   }
0106   HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fpEvtMap);
0107 }
0108 
0109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0110 
0111 template<>
0112 void Scorer<Dose>::EndOfEvent(G4HCofThisEvent*)
0113 {
0114   if (!G4RunManager::GetRunManager()->GetCurrentEvent()->IsAborted()) {
0115     fpEvtMap->add(0, fpScorer->fDosesCutOff);
0116   }
0117   fpScorer->fCumulatedDose = 0.;
0118 }
0119 
0120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0121 
0122 template<>
0123 G4bool Scorer<Dose>::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0124 {
0125   auto currentEvent = G4EventManager::GetEventManager();
0126   const G4Track* track = aStep->GetTrack();
0127   auto boundingBox = fpChemistryWorld->GetChemistryBoundary();
0128   G4double V = boundingBox->Volume() / cm3;
0129   G4double edep = aStep->GetTotalEnergyDeposit();
0130   if (edep == 0.) {
0131     return false;
0132   }
0133   (fpScorer->fCumulatedDose) += edep;
0134   if (track->GetParentID() == 0 && aStep->IsFirstStepInVolume()) {
0135     G4double DoseInGray = ((fpScorer->fCumulatedDose) / eV) / (0.001 * V * 6.242e+18);
0136     if (DoseInGray > fpScorer->fDosesCutOff / gray) {
0137       G4cout << "_____________________________________________________________________________"
0138              << G4endl;
0139       auto name = currentEvent->GetConstCurrentEvent()
0140                     ->GetPrimaryVertex()
0141                     ->GetPrimary()
0142                     ->GetParticleDefinition()
0143                     ->GetParticleName();
0144       auto energy =
0145         currentEvent->GetConstCurrentEvent()->GetPrimaryVertex()->GetPrimary()->GetKineticEnergy();
0146       G4cout << "Stop this beam line (" << name << ", " << energy
0147              << " MeV) at actual dose: " << DoseInGray
0148              << " Gy. Cut-off dose: " << fpScorer->fDosesCutOff / gray << " Gy" << G4endl;
0149       G4cout << "The beam of " << 1000000 - currentEvent->GetStackManager()->GetNUrgentTrack() - 1
0150              << " tracks"
0151              << " in a volume of " << V * 1e+12  // convert cm3 to um3
0152              << " um3. Total deposit energy: " << fpScorer->fCumulatedDose / eV << " eV. "
0153              << G4endl;
0154       if (DoseInGray > fpScorer->fDosesToAbort / gray) {
0155         G4cout << "Abort this beam line (" << name << ", " << energy
0156                << " MeV) at actual dose: " << DoseInGray << " Gy." << G4endl;
0157         G4RunManager::GetRunManager()->AbortEvent();
0158       }
0159       G4cout << "_____________________________________________________________________________"
0160              << G4endl;
0161       auto myTrack = ((G4Track*)track);
0162       myTrack->SetTrackStatus(fStopAndKill);
0163       auto secondaries = track->GetStep()->GetSecondaryInCurrentStep();
0164       if (!secondaries->empty()) {
0165         for (auto it : *(secondaries)) {
0166           if (it != nullptr) {
0167             ((G4Track*)it)->SetTrackStatus(fStopAndKill);
0168           }
0169         }
0170       }
0171       currentEvent->GetStackManager()->ClearUrgentStack();
0172     }
0173   }
0174   return true;
0175 }
0176 
0177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0178 
0179 // Gvalues
0180 Gvalues::Gvalues()
0181   : G4UImessenger(),
0182     fTimeLimit(G4Scheduler::Instance()->GetEndTime()),
0183     fSpeciesdir(new G4UIdirectory("/scorer/Gvalues/")),
0184     fTimeBincmd(new G4UIcmdWithAnInteger("/scorer/Gvalues/nOfTimeBins", this)),
0185     fAddTimeToRecordcmd(new G4UIcmdWithADoubleAndUnit("/scorer/Gvalues/addTimeToRecord", this))
0186 {
0187   fSpeciesdir->SetGuidance("ScoreSpecies commands");
0188   G4MoleculeCounter::Instance()->ResetCounter();
0189 }
0190 
0191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0192 
0193 void Gvalues::SetNewValue(G4UIcommand* command, G4String newValue)
0194 {
0195   if (command == fAddTimeToRecordcmd.get()) {
0196     G4double cmdTime = fAddTimeToRecordcmd->GetNewDoubleValue(newValue);
0197     if (fTimeLimit >= cmdTime) {
0198       AddTimeToRecord(cmdTime);
0199     }
0200     else {
0201       AddTimeToRecord(fTimeLimit);
0202     }
0203   }
0204   if (command == fTimeBincmd.get()) {
0205     G4int cmdBins = fTimeBincmd->GetNewIntValue(newValue);
0206     G4double timeMin = 1 * ps;
0207     G4double timeMax = G4Scheduler::Instance()->GetEndTime() - timeMin;
0208     G4double timeLogMin = std::log10(timeMin);
0209     G4double timeLogMax = std::log10(timeMax);
0210     for (int i = 0; i <= cmdBins; i++) {
0211       AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / cmdBins));
0212     }
0213   }
0214 }
0215 
0216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0217 
0218 void Gvalues::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager, const std::string& out)
0219 {
0220   analysisManager->CreateNtuple(out, out);
0221   G4cout << "NtupleID : " << fRunID << " name : " << out << G4endl;
0222   analysisManager->CreateNtupleIColumn(fRunID, "speciesID");
0223   analysisManager->CreateNtupleIColumn(fRunID, "number");
0224   analysisManager->CreateNtupleIColumn(fRunID, "nEvent");
0225   analysisManager->CreateNtupleSColumn(fRunID, "speciesName");
0226   analysisManager->CreateNtupleDColumn(fRunID, "time");
0227   analysisManager->CreateNtupleDColumn(fRunID, "sumG");
0228   analysisManager->CreateNtupleDColumn(fRunID, "sumG2");
0229   analysisManager->FinishNtuple(fRunID);
0230 
0231   for (const auto& it_map1 : fSpeciesInfoPerTime) {
0232     const InnerSpeciesMap& map2 = it_map1.second;
0233     for (auto it_map2 : map2) {
0234       double time = it_map1.first;
0235       auto species = it_map2.first;
0236       const G4String& name = species->GetName();
0237       int molID = it_map2.first->GetMoleculeID();
0238       auto number = it_map2.second.fNumber;
0239       double G = it_map2.second.fG;
0240       double G2 = it_map2.second.fG2;
0241 
0242       analysisManager->FillNtupleIColumn(fRunID, 0, molID);  // MolID
0243       analysisManager->FillNtupleIColumn(fRunID, 1, number);  // Number
0244       analysisManager->FillNtupleIColumn(fRunID, 2, fNEvent);  // Total nb events
0245       analysisManager->FillNtupleSColumn(fRunID, 3, name);  // molName
0246       analysisManager->FillNtupleDColumn(fRunID, 4, time);  // time
0247       analysisManager->FillNtupleDColumn(fRunID, 5, G);  // G
0248       analysisManager->FillNtupleDColumn(fRunID, 6, G2);  // G2
0249       analysisManager->AddNtupleRow(fRunID);
0250     }
0251   }
0252   fRunID++;
0253 }
0254 
0255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0256 
0257 template<>
0258 void Scorer<Gvalues>::clear()
0259 {
0260   fpEvtMap->clear();
0261   fpScorer->fNEvent = 0;
0262   fpScorer->fEdep = 0;
0263   fpScorer->fSpeciesInfoPerTime.clear();
0264 }
0265 
0266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0267 
0268 template<>
0269 void Scorer<Gvalues>::Initialize(G4HCofThisEvent* HCE)
0270 {
0271   fpEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
0272   if (fHCID < 0) {
0273     fHCID = GetCollectionID(0);
0274   }
0275   HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fpEvtMap);
0276 }
0277 
0278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0279 
0280 template<>
0281 void Scorer<Gvalues>::SetChemistryWorld(G4VChemistryWorld* pChemistryWorld)
0282 {
0283   fpChemistryWorld = pChemistryWorld;
0284 }
0285 
0286 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0287 
0288 template<>
0289 G4VChemistryWorld* Scorer<Gvalues>::GetChemistryWorld() const
0290 {
0291   return fpChemistryWorld;
0292 }
0293 
0294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0295 
0296 template<>
0297 G4bool Scorer<Gvalues>::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0298 {
0299   G4double edep = aStep->GetTotalEnergyDeposit();
0300   if (edep == 0.) {
0301     return FALSE;
0302   }
0303   edep *= aStep->GetPreStepPoint()->GetWeight();  // (Particle Weight)
0304   G4int index = GetIndex(aStep);
0305   fpEvtMap->add(index, edep);
0306   (fpScorer->fEdep) += edep;
0307   return TRUE;
0308 }
0309 
0310 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0311 
0312 template<>
0313 void Scorer<Gvalues>::SaveScavengerChange()
0314 {
0315   auto pScavengerMaterial =
0316     dynamic_cast<G4DNAScavengerMaterial*>(G4Scheduler::Instance()->GetScavengerMaterial());
0317   if (pScavengerMaterial == nullptr) {
0318     G4ExceptionDescription errMsg;
0319     errMsg << "pScavengerMaterial == nullptr";
0320     G4Exception("Scorer<Gvalues>::SaveScavengerChange()", "SaveScavengerChange",
0321                 FatalErrorInArgument, errMsg);
0322   }
0323   auto scavengerList = pScavengerMaterial->GetScavengerList();
0324   auto V = fpChemistryWorld->GetChemistryBoundary()->Volume();
0325 
0326   for (const auto& it : scavengerList) {
0327     if (it == G4MoleculeTable::Instance()->GetConfiguration("H2O")
0328         || G4MoleculeTable::Instance()->GetConfiguration("H3Op(B)") == it
0329         || G4MoleculeTable::Instance()->GetConfiguration("OHm(B)") == it)
0330     {
0331       continue;
0332     }
0333     for (auto time_mol : fpScorer->fTimeToRecord) {
0334       int64_t n_mol = pScavengerMaterial->GetNMoleculesAtTime(it, time_mol);
0335       if (n_mol < 0) {
0336         G4ExceptionDescription errMsg;
0337         errMsg << "SaveScavengerChange()::N molecules not valid < 0 : " << it->GetName()
0338                << " N : " << n_mol << G4endl;
0339         G4Exception("", "N<0", FatalException, errMsg);
0340       }
0341 
0342       Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[time_mol][it];
0343       molInfo.fNumber += n_mol;
0344       G4double gValue = n_mol / (Avogadro * V * 1.0e-6 /*mm3 to L*/);
0345       molInfo.fG += gValue;
0346       molInfo.fG2 += gValue * gValue;
0347     }
0348   }
0349 }
0350 
0351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0352 
0353 template<>
0354 void Scorer<Gvalues>::SaveMoleculeCounter()
0355 {
0356   if (fpEventScheduler == nullptr) {
0357     G4Exception("fpEventScheduler == nullptr", "Scorer<Gvalues>::SaveMoleculeCounter()",
0358                 FatalException, "fpEventScheduler == nullptr");
0359   }
0360   else {
0361     auto counterMap = fpEventScheduler->GetCounterMap();
0362     if (counterMap.empty()) {
0363       if (!G4MoleculeCounter::Instance()->InUse()) {
0364         G4Exception("No counter", "Scorer<Gvalues>::SaveMoleculeCounter()", JustWarning,
0365                     "G4MoleculeCounter::Instance() is not used");
0366         return;
0367       }
0368 
0369       G4MoleculeCounter::RecordedMolecules species;
0370       species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
0371       if (species.get() == nullptr) {
0372         return;
0373       }
0374       else if (species->empty()) {
0375         G4cout << "No molecule recorded, energy deposited" << G4endl;
0376         ++(fpScorer->fNEvent);
0377         fpScorer->fEdep = 0.;
0378         G4MoleculeCounter::Instance()->ResetCounter();
0379         return;
0380       }
0381       for (auto molecule : *species) {
0382         if (molecule == G4MoleculeTable::Instance()->GetConfiguration("O2")) {
0383           continue;
0384         }
0385         for (auto time_mol : fpScorer->fTimeToRecord) {
0386           int n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, time_mol);
0387 
0388           if (n_mol < 0) {
0389             G4ExceptionDescription errMsg;
0390             errMsg << "N molecules not valid < 0 " << G4endl;
0391             G4Exception("", "N<0", FatalException, errMsg);
0392           }
0393 
0394           Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[time_mol][molecule];
0395           molInfo.fNumber += n_mol;
0396           G4double gValue = (n_mol / (fpScorer->fEdep / eV)) * 100.;
0397           molInfo.fG += gValue;
0398           molInfo.fG2 += gValue * gValue;
0399         }
0400       }
0401     }
0402     else {
0403       for (const auto& map_mol : counterMap) {
0404         auto time_mol = map_mol.first;
0405         for (auto it_mol : map_mol.second) {
0406           auto molecule = it_mol.first;
0407           if (molecule == G4MoleculeTable::Instance()->GetConfiguration("O2")) {
0408             continue;
0409           }
0410           int n_mol = it_mol.second;
0411 
0412           if (n_mol < 0) {
0413             G4ExceptionDescription errMsg;
0414             errMsg << "N molecules not valid < 0 "
0415                    << " molecule : " << it_mol.first->GetName() << " N : " << n_mol << G4endl;
0416             G4Exception("", "N<0", FatalException, errMsg);
0417           }
0418 
0419           Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[time_mol][molecule];
0420           molInfo.fNumber += n_mol;
0421           G4double gValue = (n_mol / (fpScorer->fEdep / eV)) * 100.;
0422           molInfo.fG += gValue;
0423           molInfo.fG2 += gValue * gValue;
0424         }
0425       }
0426     }
0427   }
0428 }
0429 
0430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0431 
0432 template<>
0433 void Scorer<Gvalues>::EndOfEvent(G4HCofThisEvent*)
0434 {
0435   if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
0436     fpScorer->fEdep = 0.;
0437     G4MoleculeCounter::Instance()->ResetCounter();
0438     return;
0439   }
0440 
0441   SaveScavengerChange();
0442   SaveMoleculeCounter();
0443 
0444   ++(fpScorer->fNEvent);
0445   fpScorer->fEdep = 0.;
0446 
0447   G4MoleculeCounter::Instance()->ResetCounter();
0448   G4MoleculeCounter::Instance()->Use(true);
0449 }
0450 
0451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0452 
0453 template<>
0454 void Scorer<Gvalues>::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0455 {
0456   auto right = dynamic_cast<Scorer<Gvalues>*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
0457 
0458   if (right == nullptr) {
0459     return;
0460   }
0461   if (right == this) {
0462     return;
0463   }
0464   auto it_map1 = right->fpScorer->fSpeciesInfoPerTime.begin();
0465   auto end_map1 = right->fpScorer->fSpeciesInfoPerTime.end();
0466 
0467   for (; it_map1 != end_map1; ++it_map1) {
0468     Gvalues::InnerSpeciesMap& map2 = it_map1->second;
0469     auto it_map2 = map2.begin();
0470     auto end_map2 = map2.end();
0471 
0472     for (; it_map2 != end_map2; ++it_map2) {
0473       Gvalues::SpeciesInfo& molInfo = fpScorer->fSpeciesInfoPerTime[it_map1->first][it_map2->first];
0474       molInfo.fNumber += it_map2->second.fNumber;
0475       molInfo.fG += it_map2->second.fG;
0476       molInfo.fG2 += it_map2->second.fG2;
0477     }
0478   }
0479   right->fpScorer->fSpeciesInfoPerTime.clear();
0480   fpScorer->fNEvent += right->fpScorer->fNEvent;
0481   right->fpScorer->fNEvent = 0;
0482   right->fpScorer->fEdep = 0.;
0483 }
0484 
0485 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0486 
0487 template<>
0488 void Scorer<Gvalues>::OutputAndClear(const std::string& dose)
0489 {
0490   if (G4Threading::IsWorkerThread()) {
0491     return;
0492   }
0493   G4VAnalysisManager* analysisManager = G4AnalysisManager::Instance();
0494   if (analysisManager != nullptr) {
0495     this->fpScorer->WriteWithAnalysisManager(analysisManager, dose);
0496   }
0497   fpScorer->fNEvent = 0;
0498   fpScorer->fSpeciesInfoPerTime.clear();
0499 }
0500 
0501 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....