File indexing completed on 2025-04-10 08:06:46
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 G4MoleculeCounter::Instance()->ResetCounter();
0059 }
0060
0061
0062
0063 void ScoreSpecies::SetNewValue(G4UIcommand* command, G4String newValue)
0064 {
0065 if (command == fpAddTimeToRecordcmd.get()) {
0066 G4double cmdTime = fpAddTimeToRecordcmd->GetNewDoubleValue(newValue);
0067 AddTimeToRecord(cmdTime);
0068 }
0069 if (command == fpTimeBincmd.get()) {
0070 ClearTimeToRecord();
0071 G4int cmdBins = fpTimeBincmd->GetNewIntValue(newValue);
0072 G4double timeMin = 1 * ps;
0073 G4double timeMax = G4Scheduler::Instance()->GetEndTime() - timeMin;
0074 G4double timeLogMin = std::log10(timeMin);
0075 G4double timeLogMax = std::log10(timeMax);
0076 for (G4int i = 0; i < cmdBins; i++) {
0077 AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / (cmdBins - 1)));
0078 }
0079 }
0080 if (command == fpSetResultsFileNameCmd.get()) {
0081 fRootFileName = newValue;
0082 }
0083 }
0084
0085
0086
0087 G4bool ScoreSpecies::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0088 {
0089 G4double edep = aStep->GetTotalEnergyDeposit();
0090 if (edep == 0.) {
0091 return FALSE;
0092 }
0093 edep *= aStep->GetPreStepPoint()->GetWeight();
0094 auto index = GetIndex(aStep);
0095 fEvtMap->add(index, edep);
0096 fEdep += edep;
0097 return TRUE;
0098 }
0099
0100
0101
0102 void ScoreSpecies::Initialize(G4HCofThisEvent* HCE)
0103 {
0104 fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
0105 if (fHCID < 0) {
0106 fHCID = GetCollectionID(0);
0107 }
0108 HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fEvtMap);
0109 G4MoleculeCounter::Instance()->ResetCounter();
0110 }
0111
0112
0113
0114 void ScoreSpecies::EndOfEvent(G4HCofThisEvent*)
0115 {
0116 if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
0117 fEdep = 0.;
0118 G4MoleculeCounter::Instance()->ResetCounter();
0119 return;
0120 }
0121
0122 auto species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
0123
0124 if (species == nullptr || species->empty()) {
0125 G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
0126 ++fNEvent;
0127 fEdep = 0.;
0128 G4MoleculeCounter::Instance()->ResetCounter();
0129 return;
0130 }
0131 for (auto molecule : *species) {
0132 for (auto time_mol : fTimeToRecord) {
0133 G4double n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, time_mol);
0134 if (n_mol < 0) {
0135 G4cerr << "N molecules not valid < 0 " << G4endl;
0136 G4Exception("", "N<0", FatalException, "");
0137 }
0138 SpeciesInfo& molInfo = fSpeciesInfoPerTime[time_mol][molecule];
0139 molInfo.fNumber += n_mol;
0140 G4double gValue = (n_mol / (fEdep / eV)) * 100.;
0141 molInfo.fG += gValue;
0142 molInfo.fG2 += gValue * gValue;
0143 }
0144 }
0145 ++fNEvent;
0146 fEdep = 0.;
0147 G4MoleculeCounter::Instance()->ResetCounter();
0148 }
0149
0150
0151
0152 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0153 {
0154 auto right = dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
0155
0156 if (right == nullptr) {
0157 return;
0158 }
0159 if (right == this) {
0160 return;
0161 }
0162
0163 auto it_map1 = right->fSpeciesInfoPerTime.begin();
0164 auto end_map1 = right->fSpeciesInfoPerTime.end();
0165
0166 for (; it_map1 != end_map1; ++it_map1) {
0167 InnerSpeciesMap& map2 = it_map1->second;
0168 auto it_map2 = map2.begin();
0169 auto end_map2 = map2.end();
0170
0171 for (; it_map2 != end_map2; ++it_map2) {
0172 SpeciesInfo& molInfo = fSpeciesInfoPerTime[it_map1->first][it_map2->first];
0173 molInfo.fNumber += it_map2->second.fNumber;
0174 molInfo.fG += it_map2->second.fG;
0175 molInfo.fG2 += it_map2->second.fG2;
0176 }
0177 }
0178 right->fSpeciesInfoPerTime.clear();
0179 fNEvent += right->fNEvent;
0180 right->fNEvent = 0;
0181 right->fEdep = 0.;
0182 }
0183
0184
0185
0186 void ScoreSpecies::PrintAll()
0187 {
0188 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
0189 G4cout << " PrimitiveScorer " << GetName() << G4endl;
0190 G4cout << " Number of events " << fNEvent << G4endl;
0191 G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
0192
0193 for (auto itr : *fEvtMap->GetMap()) {
0194 G4cout << " copy no.: " << itr.first << " energy deposit: " << *(itr.second) / GetUnitValue()
0195 << " [" << GetUnit() << "]" << G4endl;
0196 }
0197 }
0198
0199
0200
0201 void ScoreSpecies::OutputAndClear()
0202 {
0203 if (G4Threading::IsWorkerThread()) {
0204 return;
0205 }
0206
0207 auto analysisManager = G4AnalysisManager::Instance();
0208 analysisManager->SetDefaultFileType(fOutputType);
0209 this->WriteWithAnalysisManager(analysisManager);
0210 fNEvent = 0;
0211 fSpeciesInfoPerTime.clear();
0212 }
0213
0214
0215
0216 void ScoreSpecies::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager)
0217 {
0218 analysisManager->OpenFile(fRootFileName);
0219 G4int fNtupleID = analysisManager->CreateNtuple("species", "species");
0220 analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
0221 analysisManager->CreateNtupleIColumn(fNtupleID, "number");
0222 analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
0223 analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
0224 analysisManager->CreateNtupleDColumn(fNtupleID, "time");
0225 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG");
0226 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG2");
0227 analysisManager->FinishNtuple(fNtupleID);
0228
0229 for (const auto& it_map1 : fSpeciesInfoPerTime) {
0230 const InnerSpeciesMap& map2 = it_map1.second;
0231
0232 for (const auto& it_map2 : map2) {
0233 G4double time = it_map1.first;
0234 auto species = it_map2.first;
0235 const G4String& name = species->GetName();
0236 auto molID = it_map2.first->GetMoleculeID();
0237 auto number = it_map2.second.fNumber;
0238 auto G = it_map2.second.fG;
0239 auto G2 = it_map2.second.fG2;
0240 analysisManager->FillNtupleIColumn(fNtupleID, 0, molID);
0241 analysisManager->FillNtupleIColumn(fNtupleID, 1, number);
0242 analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent);
0243 analysisManager->FillNtupleSColumn(fNtupleID, 3, name);
0244 analysisManager->FillNtupleDColumn(fNtupleID, 4, time);
0245 analysisManager->FillNtupleDColumn(fNtupleID, 5, G);
0246 analysisManager->FillNtupleDColumn(fNtupleID, 6, G2);
0247 analysisManager->AddNtupleRow(fNtupleID);
0248 }
0249 }
0250
0251 analysisManager->Write();
0252 analysisManager->CloseFile();
0253 }
0254
0255
0256
0257 }
