sdcounters/src/ScoreBasicMoleculeCounts.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 // The `molcounters` example(s) are provided as part of Geant4-DNA
0027 // and any report or published result obtained using it shall cite
0028 // the respective Geant4-DNA collaboration publications.
0029 //
0030 // Reports or results obtained using the spatially-aware `MoleculeCounter`
0031 // provided in this example, shall further cite:
0032 //
0033 // Velten & Tomé, Radiation Physics and Chemistry, 2023 (10.1016/j.radphyschem.2023.111194)
0034 //
0035 //
0036 // Author: Christian Velten (2025)
0037 //
0038
0039 #include "ScoreBasicMoleculeCounts.hh"
0040
0041 #include "G4AnalysisManager.hh"
0042 #include "G4Event.hh"
0043 #include "G4EventManager.hh"
0044 #include "G4MolecularConfiguration.hh"
0045 #include "G4MoleculeCounterManager.hh"
0046 #include "G4Scheduler.hh"
0047 #include "G4SystemOfUnits.hh"
0048 #include "G4TScoreNtupleWriter.hh"
0049 #include "G4UImessenger.hh"
0050 #include "G4UnitsTable.hh"
0051 #include "globals.hh"
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053 ScoreBasicMoleculeCounts::ScoreBasicMoleculeCounts(G4String name, G4int depth,
0054                                                    G4String moleculeCounterName)
0055   : G4VPrimitiveScorer(name, depth),
0056     G4UImessenger(),
0057     fMoleculeCounterName(moleculeCounterName),
0058     fRunID(),
0059     fNbOfScoredEvents(),
0060     fTimesToRecord(),
0061     fMoleculeCountPerIndexPerTime()
0062 {
0063   fAddTimeToRecordcmd = new G4UIcmdWithADoubleAndUnit(ScorerCommand("addTimeToRecord"), this);
0064   fTimeBincmd = new G4UIcmdWithAnInteger(ScorerCommand("nOfTimeBins"), this);
0065 }
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0067 const char* ScoreBasicMoleculeCounts::ScorerCommand(G4String command)
0068 {
0069   command = "/scorer/" + GetName() + "/" + command;
0070   return command.c_str();
0071 }
0072 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0073 ScoreBasicMoleculeCounts::~ScoreBasicMoleculeCounts()
0074 {
0075   delete fAddTimeToRecordcmd;//can be smart ptr ?
0076   delete fTimeBincmd;
0077 }
0078 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0079 void ScoreBasicMoleculeCounts::SetNewValue(G4UIcommand* command, G4String newValue)
0080 {
0081   if (command == fAddTimeToRecordcmd) {
0082     G4double cmdTime = fAddTimeToRecordcmd->GetNewDoubleValue(newValue);
0083     AddTimeToRecord(cmdTime);
0084   }
0085   if (command == fTimeBincmd) {
0086     ClearTimesToRecord();
0087     G4int cmdBins = fTimeBincmd->GetNewIntValue(newValue);
0088     G4double timeMin = 1 * ps;
0089     G4double timeMax = G4Scheduler::Instance()->GetEndTime();
0090     G4double timeLogMin = std::log10(timeMin);
0091     G4double timeLogMax = std::log10(timeMax);
0092     for (G4int i = 0; i < cmdBins; i++) {
0093       AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / (cmdBins - 1)));
0094     }
0095   }
0096 }
0097 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0098 G4bool ScoreBasicMoleculeCounts::ProcessHits(G4Step*, G4TouchableHistory*)
0099 {
0100   return true;
0101 }
0102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0103 void ScoreBasicMoleculeCounts::EndOfEvent(G4HCofThisEvent*)
0104 {
0105   if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) return;
0106
0107   auto counters = G4MoleculeCounterManager::Instance()->GetMoleculeCounters(fMoleculeCounterName);
0108
0109   G4MoleculeCounter* counter = nullptr;
0110   if (counters.rbegin() != counters.rend()) {
0111     counter =
0112       const_cast<G4MoleculeCounter*>(dynamic_cast<const G4MoleculeCounter*>(*counters.begin()));
0113     if (counter == nullptr)
0114       G4Exception("ScoreMoleculeCounter::EndOfEvent", "SCOREMOLCOUNT", FatalException,
0115                   "Molecule counter has wrong type!");
0116   }
0117   else {
0118     G4Exception("ScoreMoleculeCounter::EndOfEvent", "SCOREMOLCOUNT", FatalException,
0119                 "No molecule counter with given name found!");
0120   }
0121
0122   for (auto const& entry : counter->GetCounterMap()) {
0123     auto key = entry.first.GetMolecule()->GetName();
0124     for (auto const& time : fTimesToRecord) {
0125       G4int nbOfMoleculesAtTime = counter->GetNbMoleculesAtTime(entry.first, time);
0126       fMoleculeCountPerIndexPerTime[key][time] += nbOfMoleculesAtTime;
0127     }
0128   }
0129
0130   ++fNbOfScoredEvents;
0131 }
0132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0133 void ScoreBasicMoleculeCounts::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0134 {
0135   ScoreBasicMoleculeCounts* worker = dynamic_cast<ScoreBasicMoleculeCounts*>(workerScorer);
0136   if (worker == nullptr || worker == this) return;
0137
0138   fNbOfScoredEvents += worker->fNbOfScoredEvents;
0139
0140   for (auto const& workerEntry : worker->fMoleculeCountPerIndexPerTime) {
0141     auto emplacedPair =
0142       fMoleculeCountPerIndexPerTime.emplace(workerEntry.first, workerEntry.second);
0143     if (!emplacedPair.second) {  // workerEntry.first already exists, emplacedPair.first == it
0144       auto it = emplacedPair.first;
0145
0146       for (auto const& workerMoleculeCountWithTimeIter : workerEntry.second) {
0147         auto masterMoleculeCountWithTimeIter =
0148           it->second.find(workerMoleculeCountWithTimeIter.first);
0149         masterMoleculeCountWithTimeIter->second += workerMoleculeCountWithTimeIter.second;
0150       }
0151     }
0152   }
0153
0154   worker->clear();
0155 }
0156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0157 void ScoreBasicMoleculeCounts::clear()
0158 {
0159   G4VPrimitiveScorer::clear();
0160   fNbOfScoredEvents = 0;
0161   fMoleculeCountPerIndexPerTime.clear();
0162 }
0163
0164 void ScoreBasicMoleculeCounts::OutputAndClear()
0165 {
0166   if (G4Threading::IsWorkerThread()) return;
0167
0168   G4cout << "\n\nWriting Scorers: " << GetName() << G4endl;
0169
0170   WriteWithAnalysisManager();
0171
0172   ++fRunID;
0173   clear();
0174 }
0175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0176 void ScoreBasicMoleculeCounts::WriteWithAnalysisManager()
0177 {
0178   auto analysisManager = G4AnalysisManager::Instance();
0179
0180   int fNtupleID = analysisManager->CreateNtuple("BasicMoleculeCount_" + fMoleculeCounterName,
0181                                                 "BasicMoleculeCount_" + fMoleculeCounterName);
0182
0183   analysisManager->CreateNtupleDColumn("Time__ps_");
0184   analysisManager->CreateNtupleDColumn("Molecule_Count");
0185   analysisManager->CreateNtupleSColumn("Molecule_Name");
0186   analysisManager->FinishNtuple();
0187
0188   for (auto const& entry : fMoleculeCountPerIndexPerTime) {
0189     G4cout << entry.first << ":" << G4endl;
0190
0191     for (auto const& moleculeCountWithTime : entry.second) {
0192       auto moleculeCount = moleculeCountWithTime.second;
0193       //auto moleculeCountPerEvent = moleculeCount / fNbOfScoredEvents;//not used
0194
0195       auto time = moleculeCountWithTime.first;
0196       auto moleculeName = entry.first;
0197
0198       G4cout << "  t=" << G4BestUnit(time, "Time") << ", n=" << moleculeCount << G4endl;
0199
0200       analysisManager->FillNtupleDColumn(fNtupleID, 0, time / ps);  // time
0201       analysisManager->FillNtupleDColumn(fNtupleID, 1, moleculeCount);  // Number
0202       analysisManager->FillNtupleSColumn(fNtupleID, 2, moleculeName);  // molecule name
0203       analysisManager->AddNtupleRow(fNtupleID);
0204     }
0205   }
0206 }
0207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......