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