File indexing completed on 2025-11-04 09:30:00
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029 #include "ScoreSpecies.hh"
0030
0031 #include "G4AnalysisManager.hh"
0032 #include "G4Event.hh"
0033 #include "G4Scheduler.hh"
0034 #include "G4UImessenger.hh"
0035 #include "G4UnitsTable.hh"
0036
0037 #include <G4EventManager.hh>
0038 #include <G4MolecularConfiguration.hh>
0039 #include <G4MoleculeCounter.hh>
0040 #include <G4SystemOfUnits.hh>
0041 #include <globals.hh>
0042
0043 namespace scavenger
0044 {
0045
0046
0047
0048 ScoreSpecies::ScoreSpecies(const G4String& name, const G4int& depth)
0049 : G4VPrimitiveScorer(name, depth),
0050 G4UImessenger(),
0051 fpSpeciesdir(new G4UIdirectory("/scorer/species/")),
0052 fpTimeBincmd(new G4UIcmdWithAnInteger("/scorer/species/nOfTimeBins", this)),
0053 fpAddTimeToRecordcmd(new G4UIcmdWithADoubleAndUnit("/scorer/species/addTimeToRecord", this)),
0054 fpSetResultsFileNameCmd(new G4UIcmdWithAString("/scorer/species/setRootFileName", this))
0055 {
0056 fpSpeciesdir->SetGuidance("ScoreSpecies commands");
0057 fpSetResultsFileNameCmd->SetGuidance("Set the root file name");
0058 }
0059
0060
0061
0062 void ScoreSpecies::SetNewValue(G4UIcommand* command, G4String newValue)
0063 {
0064 if (command == fpAddTimeToRecordcmd.get()) {
0065 G4double cmdTime = fpAddTimeToRecordcmd->GetNewDoubleValue(newValue);
0066 AddTimeToRecord(cmdTime);
0067 }
0068 if (command == fpTimeBincmd.get()) {
0069 ClearTimeToRecord();
0070 G4int cmdBins = fpTimeBincmd->GetNewIntValue(newValue);
0071 G4double timeMin = 1 * ps;
0072 G4double timeMax = G4Scheduler::Instance()->GetEndTime() - timeMin;
0073 G4double timeLogMin = std::log10(timeMin);
0074 G4double timeLogMax = std::log10(timeMax);
0075 for (G4int i = 0; i < cmdBins; i++) {
0076 AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / (cmdBins - 1)));
0077 }
0078 }
0079 if (command == fpSetResultsFileNameCmd.get()) {
0080 fRootFileName = newValue;
0081 }
0082 }
0083
0084
0085
0086 G4bool ScoreSpecies::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0087 {
0088 G4double edep = aStep->GetTotalEnergyDeposit();
0089 if (edep == 0.) {
0090 return FALSE;
0091 }
0092 edep *= aStep->GetPreStepPoint()->GetWeight();
0093 auto index = GetIndex(aStep);
0094 fEvtMap->add(index, edep);
0095 fEdep += edep;
0096 return TRUE;
0097 }
0098
0099
0100
0101 void ScoreSpecies::Initialize(G4HCofThisEvent* HCE)
0102 {
0103 fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
0104 if (fHCID < 0) {
0105 fHCID = GetCollectionID(0);
0106 }
0107 HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fEvtMap);
0108 }
0109
0110
0111
0112 void ScoreSpecies::EndOfEvent(G4HCofThisEvent*)
0113 {
0114 if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
0115 fEdep = 0.;
0116 return;
0117 }
0118
0119
0120 auto counter = G4MoleculeCounterManager::Instance()->GetMoleculeCounter<G4MoleculeCounter>(0);
0121 if (counter == nullptr) {
0122 G4Exception("ScoreSpecies::EndOfEvent", "BAD_REFERENCE", FatalException,
0123 "The molecule counter could not be received!");
0124 }
0125
0126 auto indices = counter->GetMapIndices();
0127
0128 if (indices.empty()) {
0129 G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
0130 ++fNEvent;
0131 fEdep = 0.;
0132 return;
0133 }
0134 for (auto idx : indices) {
0135 for (auto time_mol : fTimeToRecord) {
0136 double n_mol = counter->GetNbMoleculesAtTime(idx, time_mol);
0137 if (n_mol < 0) {
0138 G4cerr << "N molecules not valid < 0 " << G4endl;
0139 G4Exception("", "N<0", FatalException, "");
0140 }
0141 SpeciesInfo& molInfo = fSpeciesInfoPerTime[time_mol][idx.Molecule];
0142 molInfo.fNumber += n_mol;
0143 G4double gValue = (n_mol / (fEdep / eV)) * 100.;
0144 molInfo.fG += gValue;
0145 molInfo.fG2 += gValue * gValue;
0146 }
0147 }
0148 ++fNEvent;
0149 fEdep = 0.;
0150 }
0151
0152
0153
0154 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0155 {
0156 auto right = dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
0157
0158 if (right == nullptr) {
0159 return;
0160 }
0161 if (right == this) {
0162 return;
0163 }
0164
0165 auto it_map1 = right->fSpeciesInfoPerTime.begin();
0166 auto end_map1 = right->fSpeciesInfoPerTime.end();
0167
0168 for (; it_map1 != end_map1; ++it_map1) {
0169 InnerSpeciesMap& map2 = it_map1->second;
0170 auto it_map2 = map2.begin();
0171 auto end_map2 = map2.end();
0172
0173 for (; it_map2 != end_map2; ++it_map2) {
0174 SpeciesInfo& molInfo = fSpeciesInfoPerTime[it_map1->first][it_map2->first];
0175 molInfo.fNumber += it_map2->second.fNumber;
0176 molInfo.fG += it_map2->second.fG;
0177 molInfo.fG2 += it_map2->second.fG2;
0178 }
0179 }
0180 right->fSpeciesInfoPerTime.clear();
0181 fNEvent += right->fNEvent;
0182 right->fNEvent = 0;
0183 right->fEdep = 0.;
0184 }
0185
0186
0187
0188 void ScoreSpecies::PrintAll()
0189 {
0190 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
0191 G4cout << " PrimitiveScorer " << GetName() << G4endl;
0192 G4cout << " Number of events " << fNEvent << G4endl;
0193 G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
0194
0195 for (auto itr : *fEvtMap->GetMap()) {
0196 G4cout << " copy no.: " << itr.first << " energy deposit: " << *(itr.second) / GetUnitValue()
0197 << " [" << GetUnit() << "]" << G4endl;
0198 }
0199 }
0200
0201
0202
0203 void ScoreSpecies::OutputAndClear()
0204 {
0205 if (G4Threading::IsWorkerThread()) {
0206 return;
0207 }
0208
0209 auto analysisManager = G4AnalysisManager::Instance();
0210 analysisManager->SetDefaultFileType(fOutputType);
0211 this->WriteWithAnalysisManager(analysisManager);
0212 fNEvent = 0;
0213 fSpeciesInfoPerTime.clear();
0214 }
0215
0216
0217
0218 void ScoreSpecies::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager)
0219 {
0220 analysisManager->OpenFile(fRootFileName);
0221 G4int fNtupleID = analysisManager->CreateNtuple("species", "species");
0222 analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
0223 analysisManager->CreateNtupleIColumn(fNtupleID, "number");
0224 analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
0225 analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
0226 analysisManager->CreateNtupleDColumn(fNtupleID, "time");
0227 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG");
0228 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG2");
0229 analysisManager->FinishNtuple(fNtupleID);
0230
0231 for (const auto& it_map1 : fSpeciesInfoPerTime) {
0232 const InnerSpeciesMap& map2 = it_map1.second;
0233
0234 for (const auto& it_map2 : map2) {
0235 G4double time = it_map1.first;
0236 auto species = it_map2.first;
0237 const G4String& name = species->GetName();
0238 auto molID = it_map2.first->GetMoleculeID();
0239 auto number = it_map2.second.fNumber;
0240 auto G = it_map2.second.fG;
0241 auto G2 = it_map2.second.fG2;
0242 analysisManager->FillNtupleIColumn(fNtupleID, 0, molID);
0243 analysisManager->FillNtupleIColumn(fNtupleID, 1, number);
0244 analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent);
0245 analysisManager->FillNtupleSColumn(fNtupleID, 3, name);
0246 analysisManager->FillNtupleDColumn(fNtupleID, 4, time);
0247 analysisManager->FillNtupleDColumn(fNtupleID, 5, G);
0248 analysisManager->FillNtupleDColumn(fNtupleID, 6, G2);
0249 analysisManager->AddNtupleRow(fNtupleID);
0250 }
0251 }
0252
0253 analysisManager->Write();
0254 analysisManager->CloseFile();
0255 }
0256
0257
0258
0259 }
