sdcounters/src/ScoreBasicReactionCounts.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 "ScoreBasicReactionCounts.hh"
0040
0041 #include "G4AnalysisManager.hh"
0042 #include "G4Event.hh"
0043 #include "G4EventManager.hh"
0044 #include "G4MolecularConfiguration.hh"
0045 #include "G4MoleculeCounterManager.hh"
0046 #include "G4MoleculeReactionCounter.hh"
0047 #include "G4Scheduler.hh"
0048 #include "G4SystemOfUnits.hh"
0049 #include "G4TScoreNtupleWriter.hh"
0050 #include "G4UImessenger.hh"
0051 #include "G4UnitsTable.hh"
0052 #include "globals.hh"
0053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0054 ScoreBasicReactionCounts::ScoreBasicReactionCounts(G4String name, G4int depth,
0055                                                    G4String moleculeCounterName)
0056   : G4VPrimitiveScorer(name, depth),
0057     G4UImessenger(),
0058     fMoleculeCounterName(moleculeCounterName),
0059     fRunID(),
0060     fNbOfScoredEvents(),
0061     fTimesToRecord(),
0062     fReactionCountPerIndexPerTime()
0063 {
0064   fAddTimeToRecordcmd =
0065     new G4UIcmdWithADoubleAndUnit("/scorer/basicreactioncounts/addTimeToRecord", this);
0066   fTimeBincmd = new G4UIcmdWithAnInteger("/scorer/basicreactioncounts/nOfTimeBins", this);
0067 }
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0069 ScoreBasicReactionCounts::~ScoreBasicReactionCounts()
0070 {
0071   delete fAddTimeToRecordcmd;//cannot be unique_ptr ?
0072   delete fTimeBincmd;
0073 }
0074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0075 void ScoreBasicReactionCounts::SetNewValue(G4UIcommand* command, G4String newValue)
0076 {
0077   if (command == fAddTimeToRecordcmd) {
0078     G4double cmdTime = fAddTimeToRecordcmd->GetNewDoubleValue(newValue);
0079     AddTimeToRecord(cmdTime);
0080   }
0081   if (command == fTimeBincmd) {
0082     ClearTimesToRecord();
0083     G4int cmdBins = fTimeBincmd->GetNewIntValue(newValue);
0084     G4double timeMin = 1 * ps;
0085     G4double timeMax = G4Scheduler::Instance()->GetEndTime();
0086     G4double timeLogMin = std::log10(timeMin);
0087     G4double timeLogMax = std::log10(timeMax);
0088     for (G4int i = 0; i < cmdBins; i++) {
0089       AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / (cmdBins - 1)));
0090     }
0091   }
0092 }
0093 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0094 G4bool ScoreBasicReactionCounts::ProcessHits(G4Step*, G4TouchableHistory*)
0095 {
0096   return true;
0097 }
0098 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0099 void ScoreBasicReactionCounts::EndOfEvent(G4HCofThisEvent*)
0100 {
0101   if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) return;
0102
0103   auto counters =
0104     G4MoleculeCounterManager::Instance()->GetMoleculeReactionCounters(fMoleculeCounterName);
0105
0106   G4MoleculeReactionCounter *counter = nullptr;
0107   if (counters.rbegin() != counters.rend()) {
0108     counter = const_cast<G4MoleculeReactionCounter*>(
0109       dynamic_cast<const G4MoleculeReactionCounter*>(*counters.begin()));
0110     if (counter == nullptr)
0111       G4Exception("ScoreReactionCounter::EndOfEvent", "SCOREMOLCOUNT", FatalException,
0112                   "Molecule counter has wrong type!");
0113   }
0114   else {
0115     G4Exception("ScoreReactionCounter::EndOfEvent", "SCOREMOLCOUNT", FatalException,
0116                 "No molecule counter with given name found!");
0117   }
0118
0119   for (auto const& entry : counter->GetCounterMap()) {
0120     auto key = entry.first.GetInfo();
0121     for (auto const& time : fTimesToRecord) {
0122       G4int nbOfMoleculesAtTime = counter->GetNbReactionsAtTime(entry.first, time);
0123       fReactionCountPerIndexPerTime[key][time] += nbOfMoleculesAtTime;
0124     }
0125   }
0126
0127   ++fNbOfScoredEvents;
0128 }
0129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0130 void ScoreBasicReactionCounts::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0131 {
0132   auto worker = dynamic_cast<ScoreBasicReactionCounts*>(workerScorer);
0133   if (worker == nullptr || worker == this) return;
0134
0135   fNbOfScoredEvents += worker->fNbOfScoredEvents;
0136
0137   for (auto const& workerEntry : worker->fReactionCountPerIndexPerTime) {
0138     auto emplacedPair =
0139       fReactionCountPerIndexPerTime.emplace(workerEntry.first, workerEntry.second);
0140     if (!emplacedPair.second) {  // workerEntry.first already exists, emplacedPair.first == it
0141       auto it = emplacedPair.first;
0142
0143       for (auto const& workerMoleculeCountWithTimeIter : workerEntry.second) {
0144         auto masterMoleculeCountWithTimeIter =
0145           it->second.find(workerMoleculeCountWithTimeIter.first);
0146         masterMoleculeCountWithTimeIter->second += workerMoleculeCountWithTimeIter.second;
0147       }
0148     }
0149   }
0150
0151   worker->clear();
0152 }
0153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0154 void ScoreBasicReactionCounts::clear()
0155 {
0156   G4VPrimitiveScorer::clear();
0157   fNbOfScoredEvents = 0;
0158   fReactionCountPerIndexPerTime.clear();
0159 }
0160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0161 void ScoreBasicReactionCounts::OutputAndClear()
0162 {
0163   if (G4Threading::IsWorkerThread()) return;
0164
0165   G4cout << "\n\nWriting Scorers: " << GetName() << G4endl;
0166
0167   WriteWithAnalysisManager();
0168
0169   ++fRunID;
0170   clear();
0171 }
0172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0173 void ScoreBasicReactionCounts::WriteWithAnalysisManager()
0174 {
0175   auto analysisManager = G4AnalysisManager::Instance();
0176
0177   int fNtupleID = analysisManager->CreateNtuple("BasicReactionCount_" + fMoleculeCounterName,
0178                                                 "BasicReactionCount_" + fMoleculeCounterName);
0179   analysisManager->CreateNtupleDColumn("Time__ps_");
0180   analysisManager->CreateNtupleDColumn("Reaction_Count");
0181   analysisManager->CreateNtupleSColumn("Reaction_Name");
0182   analysisManager->FinishNtuple();
0183
0184   for (auto const& entry : fReactionCountPerIndexPerTime) {
0185     G4cout << entry.first << ":" << G4endl;
0186
0187     for (auto const& reactionCountWithTime : entry.second) {
0188       auto reactionCount = reactionCountWithTime.second;
0189       auto time = reactionCountWithTime.first;
0190       auto reactionName = entry.first;
0191
0192       G4cout << "  t=" << G4BestUnit(time, "Time") << ", n=" << reactionCount << G4endl;
0193
0194       analysisManager->FillNtupleDColumn(fNtupleID, 0, time / ps);  // time
0195       analysisManager->FillNtupleDColumn(fNtupleID, 1, reactionCount);  // Number
0196       analysisManager->FillNtupleSColumn(fNtupleID, 2, reactionName);  // molecule name
0197       analysisManager->AddNtupleRow(fNtupleID);
0198     }
0199   }
0200 }
0201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......