chem5/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 // Phys. Med. Biol. 63(10) (2018) 105014-12pp
0032 // The Geant4-DNA web site is available at http://geant4-dna.org
0033 //
0034 // ScoreSpecies.cc
0035 //
0036 #include "ScoreSpecies.hh"
0037
0038 #include "G4Event.hh"
0039 #include "G4UnitsTable.hh"
0040
0041 #include <G4EventManager.hh>
0042 #include <G4MolecularConfiguration.hh>
0043 #include <G4MoleculeCounter.hh>
0044 #include <G4SystemOfUnits.hh>
0045 #include <globals.hh>
0046 #include <iomanip>
0047
0048 /**
0049  \file ScoreSpecies.cc
0050  \class ScoreSpecies
0051   This is a primitive scorer class for molecular species.
0052   The number of species is recorded for all times (log spaced). It
0053   also scores the energy deposition in order to compute the
0054   radiochemical yields.
0055 */
0056
0057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0058
0059 ScoreSpecies::ScoreSpecies(G4String name, G4int depth)
0060   : G4VPrimitiveScorer(name, depth)
0061 {
0062   auto tMin = 1.0 * CLHEP::picosecond;
0063   auto tMax = 999999 * CLHEP::picosecond;
0064   auto tLogMin = std::log10(tMin);
0065   auto tLogMax = std::log10(tMax);
0066   auto tBins = 50;
0067   for (G4int i = 0; i <= tBins; i++)
0068     AddTimeToRecord(std::pow(10., tLogMin + i * (tLogMax - tLogMin) / tBins));
0069 }
0070
0071 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0072
0073 G4bool ScoreSpecies::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0074 {
0075   G4double edep = aStep->GetTotalEnergyDeposit();
0076
0077   if (edep == 0.) return FALSE;
0078
0079   edep *= aStep->GetPreStepPoint()->GetWeight();
0080   G4int index = GetIndex(aStep);
0081   fEvtMap->add(index, edep);
0082   fEdep += edep;
0083
0084   return TRUE;
0085 }
0086
0087 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0088
0089 void ScoreSpecies::Initialize(G4HCofThisEvent* HCE)
0090 {
0091   fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
0092
0093   if (fHCID < 0) {
0094     fHCID = GetCollectionID(0);
0095   }
0096
0097   HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fEvtMap);
0098 }
0099
0100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0101
0102 void ScoreSpecies::EndOfEvent(G4HCofThisEvent*)
0103 {
0104   if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
0105     fEdep = 0.;
0106     return;
0107   }
0108
0109   // get the first, and in this case only, counter
0110   auto counter = G4MoleculeCounterManager::Instance()->GetMoleculeCounter<G4MoleculeCounter>(0);
0111   if (counter == nullptr) {
0112     G4Exception("ScoreSpecies::EndOfEvent", "BAD_REFERENCE", FatalException,
0113                 "The molecule counter could not be received!");
0114   }
0115
0116   auto indices = counter->GetMapIndices();
0117
0118   if (indices.empty()) {
0119     G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
0120     ++fNEvent;
0121     fEdep = 0.;
0122     return;
0123   }
0124
0125   for (const auto& idx : indices) {
0126     for (auto time_mol : fTimeToRecord) {
0127       double n_mol = counter->GetNbMoleculesAtTime(idx, time_mol);
0128
0129       if (n_mol < 0) {
0130         G4cerr << "N molecules not valid < 0 " << G4endl;
0131         G4Exception("", "N<0", FatalException, "");
0132       }
0133
0134       SpeciesInfo& molInfo = fSpeciesInfoPerTime[time_mol][idx.Molecule];
0135       molInfo.fNumber += n_mol;
0136       G4double gValue = (n_mol / (fEdep / eV)) * 100.;
0137       molInfo.fG += gValue;
0138       molInfo.fG2 += gValue * gValue;
0139     }
0140   }
0141   ++fNEvent;
0142
0143   fEdep = 0.;
0144 }
0145
0146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0147
0148 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0149 {
0150   auto right =
0151     dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
0152
0153   if (right == nullptr) {
0154     return;
0155   }
0156   if (right == this) {
0157     return;
0158   }
0159
0160   auto it_map1 = right->fSpeciesInfoPerTime.begin();
0161   auto end_map1 = right->fSpeciesInfoPerTime.end();
0162
0163   for (; it_map1 != end_map1; ++it_map1) {
0164     InnerSpeciesMap& map2 = it_map1->second;
0165     auto it_map2 = map2.begin();
0166     auto end_map2 = map2.end();
0167
0168     for (; it_map2 != end_map2; ++it_map2) {
0169       SpeciesInfo& molInfo = fSpeciesInfoPerTime[it_map1->first][it_map2->first];
0170       molInfo.fNumber += it_map2->second.fNumber;
0171       molInfo.fG += it_map2->second.fG;
0172       molInfo.fG2 += it_map2->second.fG2;
0173     }
0174   }
0175   right->fSpeciesInfoPerTime.clear();
0176
0177   fNEvent += right->fNEvent;
0178   right->fNEvent = 0;
0179   right->fEdep = 0.;
0180 }
0181
0182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0183
0184 void ScoreSpecies::PrintAll()
0185 {
0186   G4cout << " MultiFunctionalDet  " << detector->GetName() << G4endl;
0187   G4cout << " PrimitiveScorer " << GetName() << G4endl;
0188   G4cout << " Number of events " << fNEvent << G4endl;
0189   G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
0190
0191   for (auto itr : *fEvtMap->GetMap()) {
0192     G4cout << "  copy no.: " << itr.first << "  energy deposit: " << *(itr.second) / GetUnitValue()
0193            << " [" << GetUnit() << "]" << G4endl;
0194   }
0195 }
0196
0197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0198
0199 void ScoreSpecies::ASCII()
0200 {
0201   std::ofstream out("Species.txt");
0202   if (!out) return;
0203
0204   out << "# Time [ps]  G-value (/100 eV)  RMS   Molecule" << G4endl;
0205
0206   std::map<G4String, std::map<G4double, std::pair<G4double, G4double>>> mol;
0207
0208   for (auto& it_map1 : fSpeciesInfoPerTime) {
0209     InnerSpeciesMap& map2 = it_map1.second;
0210     G4double time = it_map1.first / ps;
0211     for (auto& it_map2 : map2) {
0212       G4double G = it_map2.second.fG;
0213       G4double G2 = it_map2.second.fG2;
0214       G4double N = fNEvent;
0215       G /= N;
0216       G2 = std::sqrt(N / (N - 1) * (G2 / N - G * G));
0217       mol[it_map2.first->GetName()][time] = std::make_pair(G, G2);
0218     }
0219   }
0220
0221   for (const auto& it1 : mol)
0222     for (auto it2 : it1.second)
0223       out << std::setw(12) << it2.first << std::setw(12) << it2.second.first << std::setw(12)
0224           << it2.second.second << std::setw(12) << std::setw(12) << it1.first << G4endl;
0225
0226   out.close();
0227 }
0228
0229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0230
0231 void ScoreSpecies::OutputAndClear()
0232 {
0233   if (G4Threading::IsWorkerThread()) return;
0234
0235   //----------------------------------------------------------------------------
0236   // Save results in ASCII format
0237   ASCII();
0238
0239   fNEvent = 0;
0240   fSpeciesInfoPerTime.clear();
0241 }