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