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     G4MoleculeCounter::Instance()->ResetCounter();
0119     return;
0120   }
0121
0122   auto species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
0123
0124   if (species.get() == 0 || species->size() == 0) {
0125     G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
0126     ++fNEvent;
0127     fEdep = 0.;
0128     G4MoleculeCounter::Instance()->ResetCounter();
0129     return;
0130   }
0131
0132   //  G4cout << "ScoreSpecies::EndOfEvent"<<G4endl;
0133   //  G4cout << "End of event, deposited energy: "
0134   //  << G4BestUnit(fEdep, "Energy") << G4endl;
0135
0136 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
0137   int eventID = G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID();
0138 #endif
0139
0140   for (auto molecule : *species) {
0141     for (auto time_mol : fTimeToRecord) {
0142       double n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, time_mol);
0143
0144       if (n_mol < 0) {
0145         G4cerr << "N molecules not valid < 0 " << G4endl;
0146         G4Exception("", "N<0", FatalException, "");
0147       }
0148
0149       SpeciesInfo& molInfo = fSpeciesInfoPerTime[time_mol][molecule];
0150       molInfo.fNumber += n_mol;
0151       double gValue = (n_mol / (fEdep / eV)) * 100.;
0152       molInfo.fG += gValue;
0153       molInfo.fG2 += gValue * gValue;
0154
0155 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
0156       SpeciesInfoSOA& molInfoPerEvent = fSpeciesInfoPerEvent[time_mol][molecule];
0157       molInfoPerEvent.fNumber.push_back(n_mol);
0158       molInfoPerEvent.fG.push_back(gValue);
0159       molInfoPerEvent.fG2.push_back(gValue * gValue);
0160       molInfoPerEvent.fEventID.push_back(eventID);
0161 #endif
0162       //      G4cout << "In Save molucule: fNumber " << molInfo.fNumber
0163       //            << " fG " << molInfo.fG
0164       //            << " fEdep " << fEdep/eV
0165       //            << G4endl;
0166     }
0167   }
0168
0169   ++fNEvent;
0170
0171   //  G4cout << "End of event " << fNEvent
0172   //         << ", energy deposited=" << G4BestUnit(fEdep, "Energy") << G4endl;
0173
0174   fEdep = 0.;
0175   G4MoleculeCounter::Instance()->ResetCounter();
0176 }
0177
0178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0179
0180 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0181 {
0182   ScoreSpecies* right =
0183     dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
0184
0185   if (right == 0) {
0186     return;
0187   }
0188   if (right == this) {
0189     return;
0190   }
0191
0192   // G4cout<<"ScoreSpecies::AbsorbResultsFromWorkerScorer"<<G4endl;
0193   {
0194     SpeciesMap::iterator it_map1 = right->fSpeciesInfoPerTime.begin();
0195     SpeciesMap::iterator end_map1 = right->fSpeciesInfoPerTime.end();
0196
0197     for (; it_map1 != end_map1; ++it_map1) {
0198       InnerSpeciesMap& map2 = it_map1->second;
0199       InnerSpeciesMap::iterator it_map2 = map2.begin();
0200       InnerSpeciesMap::iterator end_map2 = map2.end();
0201
0202       for (; it_map2 != end_map2; ++it_map2) {
0203         SpeciesInfo& molInfo = fSpeciesInfoPerTime[it_map1->first][it_map2->first];
0204         molInfo.fNumber += it_map2->second.fNumber;
0205         molInfo.fG += it_map2->second.fG;
0206         molInfo.fG2 += it_map2->second.fG2;
0207
0208         //      G4cout << "In AbsorbeResultsFromWorkerScorer: fNumber "
0209         //             << molInfo.fNumber
0210         //             << " fG "
0211         //             << molInfo.fG
0212         //             << G4endl;
0213       }
0214     }
0215   }
0216   //---------------------------------------------------------
0217 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
0218   {
0219     SpeciesMapPerEvent::iterator it_map1 = right->fSpeciesInfoPerEvent.begin();
0220     SpeciesMapPerEvent::iterator end_map1 = right->fSpeciesInfoPerEvent.end();
0221
0222     for (; it_map1 != end_map1; ++it_map1) {
0223       auto& map2 = it_map1->second;
0224       InnerSpeciesMapPerEvent::iterator it_map2 = map2.begin();
0225       InnerSpeciesMapPerEvent::iterator end_map2 = map2.end();
0226
0227       for (; it_map2 != end_map2; ++it_map2) {
0228         SpeciesInfoSOA& molInfo = fSpeciesInfoPerEvent[it_map1->first][it_map2->first];
0229         molInfo.fNumber.insert(molInfo.fNumber.end(), it_map2->second.fNumber.begin(),
0230                                it_map2->second.fNumber.end());
0231         molInfo.fG.insert(molInfo.fG.end(), it_map2->second.fG.begin(), it_map2->second.fG.end());
0232         molInfo.fG2.insert(molInfo.fG2.end(), it_map2->second.fG2.begin(),
0233                            it_map2->second.fG2.end());
0234         molInfo.fEventID.insert(molInfo.fEventID.end(), it_map2->second.fEventID.begin(),
0235                                 it_map2->second.fEventID.end());
0236         // G4cout << "In AbsorbeResultsFromWorkerScorer: fNumber "
0237         //        << molInfo.fNumber
0238         //        << " fG "
0239         //        << molInfo.fG
0240         //        << G4endl;
0241       }
0242     }
0243     right->fSpeciesInfoPerEvent.clear();
0244   }
0245 #endif
0246
0247   fNEvent += right->fNEvent;
0248   right->fNEvent = 0;
0249   right->fEdep = 0.;
0250 }
0251
0252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0253
0254 void ScoreSpecies::DrawAll()
0255 {
0256   ;
0257 }
0258
0259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0260
0261 void ScoreSpecies::PrintAll()
0262 {
0263   G4cout << " MultiFunctionalDet  " << detector->GetName() << G4endl;
0264   G4cout << " PrimitiveScorer " << GetName() << G4endl;
0265   G4cout << " Number of events " << fNEvent << G4endl;
0266   G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
0267
0268   for (auto itr : *fEvtMap->GetMap()) {
0269     G4cout << "  copy no.: " << itr.first << "  energy deposit: " << *(itr.second) / GetUnitValue()
0270            << " [" << GetUnit() << "]" << G4endl;
0271   }
0272 }
0273
0274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0275
0276 void ScoreSpecies::ASCII()
0277 {
0278   std::ofstream out("Species.Txt");
0279   if (!out) return;
0280
0281   out << "Time is in ns" << G4endl;
0282
0283   for (auto it_map1 : fSpeciesInfoPerTime) {
0284     InnerSpeciesMap& map2 = it_map1.second;
0285
0286     out << it_map1.first << G4endl;
0287
0288     for (auto it_map2 : map2) {
0289       out << it_map2.first->GetName() << " " << it_map2.second.fNumber << G4endl;
0290     }
0291   }
0292
0293   out.close();
0294 }
0295
0296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0297
0298 void ScoreSpecies::OutputAndClear()
0299 {
0300   if (G4Threading::IsWorkerThread()) return;
0301
0302   //----------------------------------------------------------------------------
0303   // Save results
0304
0305   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0306   analysisManager->SetDefaultFileType(fOutputType);
0307
0308   if (analysisManager) {
0309     this->WriteWithAnalysisManager(analysisManager);
0310   }
0311
0312   fNEvent = 0;
0313   fSpeciesInfoPerTime.clear();
0314 }
0315
0316 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0317
0318 void ScoreSpecies::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager)
0319 {
0320   //  G4cout << "ScoreSpecies::WriteWithAnalysisManager" << G4endl;
0321   analysisManager->OpenFile("Species.root");
0322   int fNtupleID = analysisManager->CreateNtuple("species", "species");
0323   analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
0324   analysisManager->CreateNtupleIColumn(fNtupleID, "number");
0325   analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
0326   analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
0327   analysisManager->CreateNtupleDColumn(fNtupleID, "time");
0328   analysisManager->CreateNtupleDColumn(fNtupleID, "sumG");
0329   analysisManager->CreateNtupleDColumn(fNtupleID, "sumG2");
0330   analysisManager->FinishNtuple(fNtupleID);
0331
0332   for (auto it_map1 : fSpeciesInfoPerTime) {
0333     InnerSpeciesMap& map2 = it_map1.second;
0334
0335     for (auto it_map2 : map2) {
0336       double time = it_map1.first;
0337       auto species = it_map2.first;
0338       const G4String& name = species->GetName();
0339       int molID = it_map2.first->GetMoleculeID();
0340       int number = it_map2.second.fNumber;
0341       double G = it_map2.second.fG;
0342       double G2 = it_map2.second.fG2;
0343
0344       analysisManager->FillNtupleIColumn(fNtupleID, 0, molID);  // MolID
0345       analysisManager->FillNtupleIColumn(fNtupleID, 1, number);  // Number
0346       analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent);  // Total nb events
0347       analysisManager->FillNtupleSColumn(fNtupleID, 3, name);  // molName
0348       analysisManager->FillNtupleDColumn(fNtupleID, 4, time);  // time
0349       analysisManager->FillNtupleDColumn(fNtupleID, 5, G);  // G
0350       analysisManager->FillNtupleDColumn(fNtupleID, 6, G2);  // G2
0351       analysisManager->AddNtupleRow(fNtupleID);
0352     }
0353   }
0354
0355   //----------------------------------------------------------------------------
0356
0357 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
0358   fNtupleID = analysisManager->CreateNtuple("species_all", "species_all");
0359   analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
0360   analysisManager->CreateNtupleIColumn(fNtupleID, "number");
0361   analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
0362   analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
0363   analysisManager->CreateNtupleDColumn(fNtupleID, "time");
0364   analysisManager->CreateNtupleDColumn(fNtupleID, "G");
0365   analysisManager->CreateNtupleDColumn(fNtupleID, "G2");
0366   analysisManager->CreateNtupleIColumn(fNtupleID, "eventID");
0367   analysisManager->FinishNtuple(fNtupleID);
0368
0369   for (auto it_map1 : fSpeciesInfoPerEvent) {
0370     InnerSpeciesMapPerEvent& map2 = it_map1.second;
0371
0372     for (auto it_map2 : map2) {
0373       double time = it_map1.first;
0374       const Species& species = it_map2.first;
0375       const G4String& name = species->GetName();
0376       int molID = it_map2.first->GetMoleculeID();
0377
0378       size_t nG = it_map2.second.fG.size();
0379
0380       for (size_t i = 0; i < nG; ++i) {
0381         int number = it_map2.second.fNumber[i];
0382         double G = it_map2.second.fG[i];
0383         double G2 = it_map2.second.fG2[i];
0384         int eventID = it_map2.second.fEventID[i];
0385
0386         analysisManager->FillNtupleIColumn(fNtupleID, 0, molID);  // MolID
0387         analysisManager->FillNtupleIColumn(fNtupleID, 1, number);  // Number
0388         analysisManager->FillNtupleIColumn(fNtupleID, 2, fNEvent);  // Total nb events
0389         analysisManager->FillNtupleSColumn(fNtupleID, 3, name);  // molName
0390         analysisManager->FillNtupleDColumn(fNtupleID, 4, time);  // time
0391         analysisManager->FillNtupleDColumn(fNtupleID, 5, G);  // G
0392         analysisManager->FillNtupleDColumn(fNtupleID, 6, G2);  // G2
0393         analysisManager->FillNtupleIColumn(fNtupleID, 7, eventID);  // EventID
0394         analysisManager->AddNtupleRow(fNtupleID);
0395       }
0396     }
0397   }
0398 #endif
0399
0400   analysisManager->Write();
0401   analysisManager->CloseFile();
0402 }