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