File indexing completed on 2025-02-23 09:21:54
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 G4MoleculeCounter::Instance()->ResetCounter();
0084 }
0085
0086
0087
0088 ScoreSpecies::~ScoreSpecies()
0089 {
0090 delete fSpeciesdir;
0091 delete fAddTimeToRecordcmd;
0092 delete fTimeBincmd;
0093 delete fOutputTypeUI;
0094 delete fOutputFileUI;
0095 }
0096
0097
0098
0099 void ScoreSpecies::SetNewValue(G4UIcommand* command, G4String newValue)
0100 {
0101 if (command == fAddTimeToRecordcmd) {
0102 G4double cmdTime = fAddTimeToRecordcmd->GetNewDoubleValue(newValue);
0103 AddTimeToRecord(cmdTime);
0104 }
0105 if (command == fTimeBincmd) {
0106 ClearTimeToRecord();
0107 G4int cmdBins = fTimeBincmd->GetNewIntValue(newValue);
0108 G4double timeMin = 1 * ps;
0109 G4double timeMax = G4Scheduler::Instance()->GetEndTime() - 1 * ps;
0110 G4double timeLogMin = std::log10(timeMin);
0111 G4double timeLogMax = std::log10(timeMax);
0112 for (G4int i = 0; i < cmdBins; i++) {
0113 AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / (cmdBins - 1)));
0114 }
0115 }
0116
0117 if (command == fOutputTypeUI) {
0118 fOutputType = newValue;
0119 G4StrUtil::to_lower(fOutputType);
0120 }
0121
0122 if (command == fOutputFileUI) fOutputFile = newValue;
0123 }
0124
0125
0126
0127 G4bool ScoreSpecies::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0128 {
0129 G4double edep = aStep->GetTotalEnergyDeposit();
0130
0131 if (edep == 0.) return FALSE;
0132
0133 edep *= aStep->GetPreStepPoint()->GetWeight();
0134 G4int index = GetIndex(aStep);
0135 fEvtMap->add(index, edep);
0136 fEdep += edep;
0137
0138 return TRUE;
0139 }
0140
0141
0142
0143 void ScoreSpecies::Initialize(G4HCofThisEvent* HCE)
0144 {
0145 fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
0146
0147 if (fHCID < 0) {
0148 fHCID = GetCollectionID(0);
0149 }
0150
0151 HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fEvtMap);
0152 G4MoleculeCounter::Instance()->ResetCounter();
0153 }
0154
0155
0156
0157 void ScoreSpecies::EndOfEvent(G4HCofThisEvent*)
0158 {
0159 if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
0160 fEdep = 0.;
0161 G4MoleculeCounter::Instance()->ResetCounter();
0162 return;
0163 }
0164
0165 auto species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
0166
0167 if (species.get() == 0 || species->size() == 0) {
0168 G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
0169 ++fNEvent;
0170 fEdep = 0.;
0171 G4MoleculeCounter::Instance()->ResetCounter();
0172 return;
0173 }
0174
0175 for (auto molecule : *species) {
0176 for (auto time_mol : fTimeToRecord) {
0177 double n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, 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][molecule];
0185 molInfo.fNumber += n_mol;
0186 molInfo.fNumber2 += n_mol * n_mol;
0187 double gValue = (n_mol / (fEdep / eV)) * 100.;
0188 molInfo.fG += gValue;
0189 molInfo.fG2 += gValue * gValue;
0190 }
0191 }
0192
0193 ++fNEvent;
0194
0195 fEdep = 0.;
0196 G4MoleculeCounter::Instance()->ResetCounter();
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) {
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
