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