File indexing completed on 2026-05-18 07:54:41
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
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043 #include "ScoreSpecies.hh"
0044
0045 #include "G4AnalysisManager.hh"
0046 #include "G4Event.hh"
0047 #include "G4Scheduler.hh"
0048 #include "G4TScoreNtupleWriter.hh"
0049 #include "G4UImessenger.hh"
0050 #include "G4UnitsTable.hh"
0051 #include "G4SystemOfUnits.hh"
0052
0053 #include <G4EventManager.hh>
0054 #include <G4MolecularConfiguration.hh>
0055 #include <G4MoleculeCounter.hh>
0056 #include <globals.hh>
0057
0058
0059
0060
0061
0062
0063
0064
0065 extern std::ofstream out;
0066
0067
0068
0069 ScoreSpecies::ScoreSpecies(const G4String &name, const G4int depth)
0070 : G4VPrimitiveScorer(name, depth), G4UImessenger() {
0071 fSpeciesdir = std::make_unique<G4UIdirectory>("/scorer/species/");
0072 fSpeciesdir->SetGuidance("ScoreSpecies commands");
0073
0074 fAddTimeToRecordcmd = std::make_unique<G4UIcmdWithADoubleAndUnit>
0075 ("/scorer/species/addTimeToRecord", this);
0076 fTimeBincmd = std::make_unique<G4UIcmdWithAnInteger>("/scorer/species/nOfTimeBins", this);
0077 fNameCmd = std::make_unique<G4UIcmdWithAString>("/scorer/species/fileName", this);
0078 }
0079
0080
0081
0082 void ScoreSpecies::SetNewValue(G4UIcommand *command, const G4String newValue) {
0083 if (command == fAddTimeToRecordcmd.get()) {
0084 const G4double cmdTime = fAddTimeToRecordcmd->GetNewDoubleValue(newValue);
0085 AddTimeToRecord(cmdTime);
0086 }
0087 if (command == fTimeBincmd.get()) {
0088 ClearTimeToRecord();
0089 const G4int cmdBins = fTimeBincmd->GetNewIntValue(newValue);
0090 constexpr G4double timeMin = 1 * ps;
0091 const G4double timeMax = G4Scheduler::Instance()->GetEndTime() - 1 * ps;
0092 if (cmdBins <= 1) {
0093
0094 AddTimeToRecord(timeMax);
0095 } else {
0096 const G4double timeLogMin = std::log10(timeMin);
0097 const G4double timeLogMax = std::log10(timeMax);
0098 for (G4int i = 0; i < cmdBins; i++) {
0099 AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / (cmdBins - 1)));
0100 }
0101 }
0102 }
0103 if (command == fNameCmd.get()) {
0104 fFileName = newValue;
0105 }
0106 }
0107
0108
0109
0110 G4bool ScoreSpecies::ProcessHits(G4Step *aStep, G4TouchableHistory *) {
0111 G4double edep = aStep->GetTotalEnergyDeposit();
0112
0113 if (edep == 0.) { return false; }
0114
0115 edep *= aStep->GetPreStepPoint()->GetWeight();
0116 const G4int index = GetIndex(aStep);
0117 fEvtMap->add(index, edep);
0118 fEdep += edep;
0119
0120 return true;
0121 }
0122
0123
0124
0125 void ScoreSpecies::Initialize(G4HCofThisEvent *HCE) {
0126 fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
0127
0128 if (fHCID < 0) {
0129 fHCID = GetCollectionID(0);
0130 }
0131
0132 HCE->AddHitsCollection(fHCID, fEvtMap);
0133 }
0134
0135
0136
0137 void ScoreSpecies::EndOfEvent(G4HCofThisEvent *) {
0138 if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
0139 fEdep = 0.;
0140 return;
0141 }
0142
0143
0144 const auto counter = G4MoleculeCounterManager::Instance()->GetMoleculeCounter<
0145 G4MoleculeCounter>(0);
0146 if (counter == nullptr) {
0147 G4Exception("ScoreSpecies::EndOfEvent", "BAD_REFERENCE", FatalException,
0148 "The molecule counter could not be received!");
0149 }
0150
0151 const auto indices = counter->GetMapIndices();
0152
0153 if (indices.empty()) {
0154 G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
0155 ++fNEvent;
0156 fEdep = 0.;
0157 return;
0158 }
0159 for (const auto &idx: indices) {
0160 for (const auto time_mol: fTimeToRecord) {
0161 const G4double n_mol = counter->GetNbMoleculesAtTime(idx, time_mol);
0162
0163 if (n_mol < 0) {
0164 G4cerr << "N molecules not valid < 0 " << G4endl;
0165 G4Exception("", "N<0", FatalException, "");
0166 }
0167
0168 SpeciesInfo &molInfo = fSpeciesInfoPerTime[time_mol][idx.Molecule];
0169 molInfo.fNumber += n_mol;
0170
0171 if (fEdep <= 0.) {
0172
0173 continue;
0174 }
0175 const G4double gValue = (n_mol / (fEdep / eV)) * 100.;
0176 molInfo.fG += gValue;
0177 molInfo.fG2 += gValue * gValue;
0178 }
0179 }
0180
0181 ++fNEvent;
0182
0183 fEdep = 0.;
0184 }
0185
0186
0187
0188 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer *workerScorer) {
0189 auto *right = dynamic_cast<ScoreSpecies *>(workerScorer);
0190 if (right == nullptr || right == this) {
0191 return;
0192 }
0193
0194 for (const auto &[timeKey, innerMap]: right->fSpeciesInfoPerTime) {
0195 for (const auto &[speciesPtr, info]: innerMap) {
0196 SpeciesInfo &molInfo = fSpeciesInfoPerTime[timeKey][speciesPtr];
0197 molInfo.fNumber += info.fNumber;
0198 molInfo.fG += info.fG;
0199 molInfo.fG2 += info.fG2;
0200 }
0201 }
0202 right->fSpeciesInfoPerTime.clear();
0203
0204 fNEvent += right->fNEvent;
0205 right->fNEvent = 0;
0206 right->fEdep = 0.;
0207 }
0208
0209
0210
0211 void ScoreSpecies::PrintAll() {
0212 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
0213 G4cout << " PrimitiveScorer " << GetName() << G4endl;
0214 G4cout << " Number of events " << fNEvent << G4endl;
0215 G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
0216
0217 for (auto [fst, snd]: *fEvtMap->GetMap()) {
0218 G4cout << " copy no.: " << fst << " energy deposit: " << *(snd) / GetUnitValue()
0219 << " [" << GetUnit() << "]" << G4endl;
0220 }
0221 }
0222
0223
0224
0225 void ScoreSpecies::OutputAndClear() {
0226 if (G4Threading::IsWorkerThread()) { return; }
0227
0228
0229
0230
0231 G4AnalysisManager *analysisManager = G4AnalysisManager::Instance();
0232 analysisManager->SetDefaultFileType(fOutputType);
0233
0234 if (analysisManager) {
0235 this->WriteWithAnalysisManager(analysisManager);
0236 }
0237
0238 fNEvent = 0;
0239 fSpeciesInfoPerTime.clear();
0240 }
0241
0242
0243
0244 void ScoreSpecies::WriteWithAnalysisManager(G4VAnalysisManager *analysisManager) {
0245 const G4String fileN = fFileName + G4UIcommand::ConvertToString(fRunID);
0246 analysisManager->OpenFile(fileN);
0247 const int fNtupleID = analysisManager->CreateNtuple("species", "species");
0248 analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
0249 analysisManager->CreateNtupleIColumn(fNtupleID, "number");
0250 analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
0251 analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
0252 analysisManager->CreateNtupleDColumn(fNtupleID, "time");
0253 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG");
0254 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG2");
0255 analysisManager->FinishNtuple(fNtupleID);
0256
0257 for (auto it_map1: fSpeciesInfoPerTime) {
0258 InnerSpeciesMap &map2 = it_map1.second;
0259
0260 if (it_map1.first == fSpeciesInfoPerTime.begin()->first) {
0261 for (auto it_molname: map2) {
0262 const auto tmp_species = it_molname.first;
0263 out << std::setw(12) << tmp_species->GetName() << std::setw(12)
0264 << tmp_species->GetMoleculeID();
0265 }
0266 out << '\n';
0267 }
0268
0269 for (auto it_map2: map2) {
0270 const G4double time = it_map1.first;
0271 const auto species = it_map2.first;
0272 const G4String &name = species->GetName();
0273 const G4int molID = it_map2.first->GetMoleculeID();
0274 const G4int number = it_map2.second.fNumber;
0275 const G4double G = it_map2.second.fG;
0276 const G4double G2 = it_map2.second.fG2;
0277 const G4int N = fNEvent;
0278
0279 if (time == *fTimeToRecord.rbegin()) {
0280 if (N > 1) {
0281 out << std::setw(12) << G / N << std::setw(12)
0282 << std::sqrt(((G2 / N) - std::pow(G / N, 2)) / (N - 1));
0283 } else {
0284 out << std::setw(12) << G / N << std::setw(12)
0285 << std::sqrt(((G2 / N) - std::pow(G / N, 2)) / N);
0286 }
0287 }
0288
0289 analysisManager->FillNtupleIColumn(fNtupleID, 0, molID);
0290 analysisManager->FillNtupleIColumn(fNtupleID, 1, number);
0291 analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent);
0292 analysisManager->FillNtupleSColumn(fNtupleID, 3, name);
0293 analysisManager->FillNtupleDColumn(fNtupleID, 4, time);
0294 analysisManager->FillNtupleDColumn(fNtupleID, 5, G);
0295 analysisManager->FillNtupleDColumn(fNtupleID, 6, G2);
0296 analysisManager->AddNtupleRow(fNtupleID);
0297 }
0298 }
0299
0300 analysisManager->Write();
0301 analysisManager->CloseFile();
0302 fRunID++;
0303 }
0304
0305
