File indexing completed on 2025-02-23 09:21:50
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 G4MoleculeCounter::Instance()->ResetCounter();
0091 }
0092
0093
0094
0095 ScoreSpecies::~ScoreSpecies()
0096 {
0097 delete fSpeciesdir;
0098 delete fAddTimeToRecordcmd;
0099 delete fTimeBincmd;
0100 }
0101
0102
0103
0104 void ScoreSpecies::SetNewValue(G4UIcommand* command, G4String newValue)
0105 {
0106 if (command == fAddTimeToRecordcmd) {
0107 G4double cmdTime = fAddTimeToRecordcmd->GetNewDoubleValue(newValue);
0108 AddTimeToRecord(cmdTime);
0109 }
0110 if (command == fTimeBincmd) {
0111 ClearTimeToRecord();
0112 G4int cmdBins = fTimeBincmd->GetNewIntValue(newValue);
0113 G4double timeMin = 1 * ps;
0114 G4double timeMax = G4Scheduler::Instance()->GetEndTime() - 1 * ps;
0115 G4double timeLogMin = std::log10(timeMin);
0116 G4double timeLogMax = std::log10(timeMax);
0117 for (G4int i = 0; i < cmdBins; i++) {
0118 AddTimeToRecord(std::pow(10, timeLogMin + i * (timeLogMax - timeLogMin) / (cmdBins - 1)));
0119 }
0120 }
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 G4MoleculeCounter::Instance()->ResetCounter();
0151 }
0152
0153
0154
0155 void ScoreSpecies::EndOfEvent(G4HCofThisEvent*)
0156 {
0157 if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
0158 fEdep = 0.;
0159 G4MoleculeCounter::Instance()->ResetCounter();
0160 return;
0161 }
0162
0163 auto species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
0164
0165 if (species.get() == 0 || species->size() == 0) {
0166 G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
0167 ++fNEvent;
0168 fEdep = 0.;
0169 G4MoleculeCounter::Instance()->ResetCounter();
0170 return;
0171 }
0172
0173
0174
0175
0176
0177 for (auto molecule : *species) {
0178 for (auto time_mol : fTimeToRecord) {
0179 double n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, 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][molecule];
0187 molInfo.fNumber += n_mol;
0188 double gValue = (n_mol / (fEdep / eV)) * 100.;
0189 molInfo.fG += gValue;
0190 molInfo.fG2 += gValue * gValue;
0191
0192
0193
0194
0195
0196 }
0197 }
0198
0199 ++fNEvent;
0200
0201
0202
0203
0204 fEdep = 0.;
0205 G4MoleculeCounter::Instance()->ResetCounter();
0206 }
0207
0208
0209
0210 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0211 {
0212 ScoreSpecies* right =
0213 dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
0214
0215 if (right == 0) {
0216 return;
0217 }
0218 if (right == this) {
0219 return;
0220 }
0221
0222
0223 SpeciesMap::iterator it_map1 = right->fSpeciesInfoPerTime.begin();
0224 SpeciesMap::iterator end_map1 = right->fSpeciesInfoPerTime.end();
0225
0226 for (; it_map1 != end_map1; ++it_map1) {
0227 InnerSpeciesMap& map2 = it_map1->second;
0228 InnerSpeciesMap::iterator it_map2 = map2.begin();
0229 InnerSpeciesMap::iterator end_map2 = map2.end();
0230
0231 for (; it_map2 != end_map2; ++it_map2) {
0232 SpeciesInfo& molInfo = fSpeciesInfoPerTime[it_map1->first][it_map2->first];
0233 molInfo.fNumber += it_map2->second.fNumber;
0234 molInfo.fG += it_map2->second.fG;
0235 molInfo.fG2 += it_map2->second.fG2;
0236
0237
0238
0239
0240
0241
0242 }
0243 }
0244 right->fSpeciesInfoPerTime.clear();
0245
0246 fNEvent += right->fNEvent;
0247 right->fNEvent = 0;
0248 right->fEdep = 0.;
0249 }
0250
0251
0252
0253 void ScoreSpecies::DrawAll()
0254 {
0255 ;
0256 }
0257
0258
0259
0260 void ScoreSpecies::PrintAll()
0261 {
0262 G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
0263 G4cout << " PrimitiveScorer " << GetName() << G4endl;
0264 G4cout << " Number of events " << fNEvent << G4endl;
0265 G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
0266
0267 for (auto itr : *fEvtMap->GetMap()) {
0268 G4cout << " copy no.: " << itr.first << " energy deposit: " << *(itr.second) / GetUnitValue()
0269 << " [" << GetUnit() << "]" << G4endl;
0270 }
0271 }
0272
0273
0274
0275 void ScoreSpecies::OutputAndClear()
0276 {
0277 if (G4Threading::IsWorkerThread()) return;
0278
0279
0280
0281
0282 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0283 analysisManager->SetDefaultFileType(fOutputType);
0284
0285 if (analysisManager) {
0286 this->WriteWithAnalysisManager(analysisManager);
0287 }
0288
0289 fNEvent = 0;
0290 fSpeciesInfoPerTime.clear();
0291 }
0292
0293
0294
0295 void ScoreSpecies::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager)
0296 {
0297 G4String fileN = "Species" + G4UIcommand::ConvertToString(fRunID);
0298 analysisManager->OpenFile(fileN);
0299 int fNtupleID = analysisManager->CreateNtuple("species", "species");
0300 analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
0301 analysisManager->CreateNtupleIColumn(fNtupleID, "number");
0302 analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
0303 analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
0304 analysisManager->CreateNtupleDColumn(fNtupleID, "time");
0305 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG");
0306 analysisManager->CreateNtupleDColumn(fNtupleID, "sumG2");
0307 analysisManager->FinishNtuple(fNtupleID);
0308
0309 for (auto it_map1 : fSpeciesInfoPerTime) {
0310 InnerSpeciesMap& map2 = it_map1.second;
0311
0312 if (it_map1.first == fSpeciesInfoPerTime.begin()->first) {
0313 for (auto it_molname : map2) {
0314 auto tmp_species = it_molname.first;
0315 out << std::setw(12) << tmp_species->GetName() << std::setw(12)
0316 << tmp_species->GetMoleculeID();
0317 }
0318 out << '\n';
0319 }
0320
0321 for (auto it_map2 : map2) {
0322 double time = it_map1.first;
0323 auto species = it_map2.first;
0324 const G4String& name = species->GetName();
0325 int molID = it_map2.first->GetMoleculeID();
0326 int number = it_map2.second.fNumber;
0327 double G = it_map2.second.fG;
0328 double G2 = it_map2.second.fG2;
0329 G4int N = fNEvent;
0330
0331 if (time == *fTimeToRecord.rbegin()) {
0332 if (N > 1) {
0333 out << std::setw(12) << G / N << std::setw(12)
0334 << std::sqrt(((G2 / N) - std::pow(G / N, 2)) / (N - 1));
0335 }
0336 else {
0337 out << std::setw(12) << G / N << std::setw(12)
0338 << std::sqrt(((G2 / N) - std::pow(G / N, 2)) / N);
0339 }
0340 }
0341
0342 analysisManager->FillNtupleIColumn(fNtupleID, 0, molID);
0343 analysisManager->FillNtupleIColumn(fNtupleID, 1, number);
0344 analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent);
0345 analysisManager->FillNtupleSColumn(fNtupleID, 3, name);
0346 analysisManager->FillNtupleDColumn(fNtupleID, 4, time);
0347 analysisManager->FillNtupleDColumn(fNtupleID, 5, G);
0348 analysisManager->FillNtupleDColumn(fNtupleID, 6, G2);
0349 analysisManager->AddNtupleRow(fNtupleID);
0350 }
0351 }
0352
0353 analysisManager->Write();
0354 analysisManager->CloseFile();
0355 fRunID++;
0356 }
0357
0358
