chem4/src/ScoreSpecies.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // This example is provided by the Geant4-DNA collaboration
0027 // Any report or published results obtained using the Geant4-DNA software
0028 // shall cite the following Geant4-DNA collaboration publication:
0029 // Med. Phys. 37 (2010) 4692-4708
0030 // J. Comput. Phys. 274 (2014) 841-882
0031 // The Geant4-DNA web site is available at http://geant4-dna.org
0032 //
0033 // ScoreSpecies.cc
0034 //
0035 #include "ScoreSpecies.hh"
0036
0037 #include "G4Event.hh"
0038 #include "G4UnitsTable.hh"
0039
0040 #include <G4AnalysisManager.hh>
0041 #include <G4EventManager.hh>
0042 #include <G4MolecularConfiguration.hh>
0043 #include <G4MoleculeCounter.hh>
0044 #include <G4SystemOfUnits.hh>
0045 #include <globals.hh>
0046
0047 /**
0048  \file ScoreSpecies.cc
0049  \class ScoreSpecies
0050   This is a primitive scorer class for molecular species.
0051   The number of species is recorded for all times (predetermined or
0052   user chosen). It also scores the energy deposition in order to compute the
0053   radiochemical yields.
0054 */
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057
0058 ScoreSpecies::ScoreSpecies(G4String name, G4int depth)
0059   : G4VPrimitiveScorer(name, depth),
0060     fEdep(0),
0061     fOutputType("root"),  // other options: "csv", "hdf5", "xml"
0062     fHCID(-1),
0063     fEvtMap(0)
0064 {
0065   fNEvent = 0;
0066   AddTimeToRecord(1 * CLHEP::picosecond);
0067   AddTimeToRecord(10 * CLHEP::picosecond);
0068   AddTimeToRecord(100 * CLHEP::picosecond);
0069   AddTimeToRecord(1000 * CLHEP::picosecond);
0070   AddTimeToRecord(10000 * CLHEP::picosecond);
0071   AddTimeToRecord(100000 * CLHEP::picosecond);
0072   AddTimeToRecord(999999 * CLHEP::picosecond);
0073   fEdep = 0;
0074 }
0075
0076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0077
0078 ScoreSpecies::~ScoreSpecies()
0079 {
0080   ;
0081 }
0082
0083 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0084
0085 G4bool ScoreSpecies::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0086 {
0087   G4double edep = aStep->GetTotalEnergyDeposit();
0088
0089   if (edep == 0.) return FALSE;
0090
0091   edep *= aStep->GetPreStepPoint()->GetWeight();  // (Particle Weight)
0092   G4int index = GetIndex(aStep);
0093   fEvtMap->add(index, edep);
0094   fEdep += edep;
0095
0096   return TRUE;
0097 }
0098
0099 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0100
0101 void ScoreSpecies::Initialize(G4HCofThisEvent* HCE)
0102 {
0103   fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
0104
0105   if (fHCID < 0) {
0106     fHCID = GetCollectionID(0);
0107   }
0108
0109   HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fEvtMap);
0110 }
0111
0112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0113
0114 void ScoreSpecies::EndOfEvent(G4HCofThisEvent*)
0115 {
0116   if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
0117     fEdep = 0.;
0118     return;
0119   }
0120
0121   // get the first, and in this case only, counter
0122   auto counter = G4MoleculeCounterManager::Instance()->GetMoleculeCounter<G4MoleculeCounter>(0);
0123   if (counter == nullptr) {
0124     G4Exception("ScoreSpecies::EndOfEvent", "BAD_REFERENCE", FatalException,
0125                 "The molecule counter could not be received!");
0126   }
0127
0128   auto indices = counter->GetMapIndices();
0129
0130   if (indices.empty()) {
0131     G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
0132     ++fNEvent;
0133     fEdep = 0.;
0134     return;
0135   }
0136
0137   //  G4cout << "ScoreSpecies::EndOfEvent"<<G4endl;
0138   //  G4cout << "End of event, deposited energy: "
0139   //  << G4BestUnit(fEdep, "Energy") << G4endl;
0140
0141 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
0142   int eventID = G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID();
0143 #endif
0144
0145   for (auto idx : indices) {
0146     for (auto time_mol : fTimeToRecord) {
0147       double n_mol = counter->GetNbMoleculesAtTime(idx, time_mol);
0148
0149       if (n_mol < 0) {
0150         G4cerr << "N molecules not valid < 0 " << G4endl;
0151         G4Exception("", "N<0", FatalException, "");
0152       }
0153
0154       SpeciesInfo& molInfo = fSpeciesInfoPerTime[time_mol][idx.Molecule];
0155       molInfo.fNumber += n_mol;
0156       double gValue = (n_mol / (fEdep / eV)) * 100.;
0157       molInfo.fG += gValue;
0158       molInfo.fG2 += gValue * gValue;
0159
0160 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
0161       SpeciesInfoSOA& molInfoPerEvent = fSpeciesInfoPerEvent[time_mol][molecule];
0162       molInfoPerEvent.fNumber.push_back(n_mol);
0163       molInfoPerEvent.fG.push_back(gValue);
0164       molInfoPerEvent.fG2.push_back(gValue * gValue);
0165       molInfoPerEvent.fEventID.push_back(eventID);
0166 #endif
0167       //      G4cout << "In Save molucule: fNumber " << molInfo.fNumber
0168       //            << " fG " << molInfo.fG
0169       //            << " fEdep " << fEdep/eV
0170       //            << G4endl;
0171     }
0172   }
0173
0174   ++fNEvent;
0175
0176   //  G4cout << "End of event " << fNEvent
0177   //         << ", energy deposited=" << G4BestUnit(fEdep, "Energy") << G4endl;
0178
0179   fEdep = 0.;
0180 }
0181
0182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0183
0184 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0185 {
0186   ScoreSpecies* right =
0187     dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
0188
0189   if (right == 0) {
0190     return;
0191   }
0192   if (right == this) {
0193     return;
0194   }
0195
0196   // G4cout<<"ScoreSpecies::AbsorbResultsFromWorkerScorer"<<G4endl;
0197   {
0198     SpeciesMap::iterator it_map1 = right->fSpeciesInfoPerTime.begin();
0199     SpeciesMap::iterator end_map1 = right->fSpeciesInfoPerTime.end();
0200
0201     for (; it_map1 != end_map1; ++it_map1) {
0202       InnerSpeciesMap& map2 = it_map1->second;
0203       InnerSpeciesMap::iterator it_map2 = map2.begin();
0204       InnerSpeciesMap::iterator end_map2 = map2.end();
0205
0206       for (; it_map2 != end_map2; ++it_map2) {
0207         SpeciesInfo& molInfo = fSpeciesInfoPerTime[it_map1->first][it_map2->first];
0208         molInfo.fNumber += it_map2->second.fNumber;
0209         molInfo.fG += it_map2->second.fG;
0210         molInfo.fG2 += it_map2->second.fG2;
0211
0212         //      G4cout << "In AbsorbeResultsFromWorkerScorer: fNumber "
0213         //             << molInfo.fNumber
0214         //             << " fG "
0215         //             << molInfo.fG
0216         //             << G4endl;
0217       }
0218     }
0219   }
0220   //---------------------------------------------------------
0221 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
0222   {
0223     SpeciesMapPerEvent::iterator it_map1 = right->fSpeciesInfoPerEvent.begin();
0224     SpeciesMapPerEvent::iterator end_map1 = right->fSpeciesInfoPerEvent.end();
0225
0226     for (; it_map1 != end_map1; ++it_map1) {
0227       auto& map2 = it_map1->second;
0228       InnerSpeciesMapPerEvent::iterator it_map2 = map2.begin();
0229       InnerSpeciesMapPerEvent::iterator end_map2 = map2.end();
0230
0231       for (; it_map2 != end_map2; ++it_map2) {
0232         SpeciesInfoSOA& molInfo = fSpeciesInfoPerEvent[it_map1->first][it_map2->first];
0233         molInfo.fNumber.insert(molInfo.fNumber.end(), it_map2->second.fNumber.begin(),
0234                                it_map2->second.fNumber.end());
0235         molInfo.fG.insert(molInfo.fG.end(), it_map2->second.fG.begin(), it_map2->second.fG.end());
0236         molInfo.fG2.insert(molInfo.fG2.end(), it_map2->second.fG2.begin(),
0237                            it_map2->second.fG2.end());
0238         molInfo.fEventID.insert(molInfo.fEventID.end(), it_map2->second.fEventID.begin(),
0239                                 it_map2->second.fEventID.end());
0240         // G4cout << "In AbsorbeResultsFromWorkerScorer: fNumber "
0241         //        << molInfo.fNumber
0242         //        << " fG "
0243         //        << molInfo.fG
0244         //        << G4endl;
0245       }
0246     }
0247     right->fSpeciesInfoPerEvent.clear();
0248   }
0249 #endif
0250
0251   fNEvent += right->fNEvent;
0252   right->fNEvent = 0;
0253   right->fEdep = 0.;
0254 }
0255
0256 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0257
0258 void ScoreSpecies::DrawAll()
0259 {
0260   ;
0261 }
0262
0263 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0264
0265 void ScoreSpecies::PrintAll()
0266 {
0267   G4cout << " MultiFunctionalDet  " << detector->GetName() << G4endl;
0268   G4cout << " PrimitiveScorer " << GetName() << G4endl;
0269   G4cout << " Number of events " << fNEvent << G4endl;
0270   G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
0271
0272   for (auto itr : *fEvtMap->GetMap()) {
0273     G4cout << "  copy no.: " << itr.first << "  energy deposit: " << *(itr.second) / GetUnitValue()
0274            << " [" << GetUnit() << "]" << G4endl;
0275   }
0276 }
0277
0278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0279
0280 void ScoreSpecies::ASCII()
0281 {
0282   std::ofstream out("Species.Txt");
0283   if (!out) return;
0284
0285   out << "Time is in ns" << G4endl;
0286
0287   for (auto it_map1 : fSpeciesInfoPerTime) {
0288     InnerSpeciesMap& map2 = it_map1.second;
0289
0290     out << it_map1.first << G4endl;
0291
0292     for (auto it_map2 : map2) {
0293       out << it_map2.first->GetName() << " " << it_map2.second.fNumber << G4endl;
0294     }
0295   }
0296
0297   out.close();
0298 }
0299
0300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0301
0302 void ScoreSpecies::OutputAndClear()
0303 {
0304   if (G4Threading::IsWorkerThread()) return;
0305
0306   //----------------------------------------------------------------------------
0307   // Save results
0308
0309   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0310   analysisManager->SetDefaultFileType(fOutputType);
0311
0312   if (analysisManager) {
0313     this->WriteWithAnalysisManager(analysisManager);
0314   }
0315
0316   fNEvent = 0;
0317   fSpeciesInfoPerTime.clear();
0318 }
0319
0320 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0321
0322 void ScoreSpecies::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager)
0323 {
0324   //  G4cout << "ScoreSpecies::WriteWithAnalysisManager" << G4endl;
0325   analysisManager->OpenFile("Species.root");
0326   int fNtupleID = analysisManager->CreateNtuple("species", "species");
0327   analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
0328   analysisManager->CreateNtupleIColumn(fNtupleID, "number");
0329   analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
0330   analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
0331   analysisManager->CreateNtupleDColumn(fNtupleID, "time");
0332   analysisManager->CreateNtupleDColumn(fNtupleID, "sumG");
0333   analysisManager->CreateNtupleDColumn(fNtupleID, "sumG2");
0334   analysisManager->FinishNtuple(fNtupleID);
0335
0336   for (auto it_map1 : fSpeciesInfoPerTime) {
0337     InnerSpeciesMap& map2 = it_map1.second;
0338
0339     for (auto it_map2 : map2) {
0340       double time = it_map1.first;
0341       auto species = it_map2.first;
0342       const G4String& name = species->GetName();
0343       int molID = it_map2.first->GetMoleculeID();
0344       int number = it_map2.second.fNumber;
0345       double G = it_map2.second.fG;
0346       double G2 = it_map2.second.fG2;
0347
0348       analysisManager->FillNtupleIColumn(fNtupleID, 0, molID);  // MolID
0349       analysisManager->FillNtupleIColumn(fNtupleID, 1, number);  // Number
0350       analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent);  // Total nb events
0351       analysisManager->FillNtupleSColumn(fNtupleID, 3, name);  // molName
0352       analysisManager->FillNtupleDColumn(fNtupleID, 4, time);  // time
0353       analysisManager->FillNtupleDColumn(fNtupleID, 5, G);  // G
0354       analysisManager->FillNtupleDColumn(fNtupleID, 6, G2);  // G2
0355       analysisManager->AddNtupleRow(fNtupleID);
0356     }
0357   }
0358
0359   //----------------------------------------------------------------------------
0360
0361 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
0362   fNtupleID = analysisManager->CreateNtuple("species_all", "species_all");
0363   analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
0364   analysisManager->CreateNtupleIColumn(fNtupleID, "number");
0365   analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
0366   analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
0367   analysisManager->CreateNtupleDColumn(fNtupleID, "time");
0368   analysisManager->CreateNtupleDColumn(fNtupleID, "G");
0369   analysisManager->CreateNtupleDColumn(fNtupleID, "G2");
0370   analysisManager->CreateNtupleIColumn(fNtupleID, "eventID");
0371   analysisManager->FinishNtuple(fNtupleID);
0372
0373   for (auto it_map1 : fSpeciesInfoPerEvent) {
0374     InnerSpeciesMapPerEvent& map2 = it_map1.second;
0375
0376     for (auto it_map2 : map2) {
0377       double time = it_map1.first;
0378       const Species& species = it_map2.first;
0379       const G4String& name = species->GetName();
0380       int molID = it_map2.first->GetMoleculeID();
0381
0382       size_t nG = it_map2.second.fG.size();
0383
0384       for (size_t i = 0; i < nG; ++i) {
0385         int number = it_map2.second.fNumber[i];
0386         double G = it_map2.second.fG[i];
0387         double G2 = it_map2.second.fG2[i];
0388         int eventID = it_map2.second.fEventID[i];
0389
0390         analysisManager->FillNtupleIColumn(fNtupleID, 0, molID);  // MolID
0391         analysisManager->FillNtupleIColumn(fNtupleID, 1, number);  // Number
0392         analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent);  // Total nb events
0393         analysisManager->FillNtupleSColumn(fNtupleID, 3, name);  // molName
0394         analysisManager->FillNtupleDColumn(fNtupleID, 4, time);  // time
0395         analysisManager->FillNtupleDColumn(fNtupleID, 5, G);  // G
0396         analysisManager->FillNtupleDColumn(fNtupleID, 6, G2);  // G2
0397         analysisManager->FillNtupleIColumn(fNtupleID, 7, eventID);  // EventID
0398         analysisManager->AddNtupleRow(fNtupleID);
0399       }
0400     }
0401   }
0402 #endif
0403
0404   analysisManager->Write();
0405   analysisManager->CloseFile();
0406 }