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), fEdep(0), fHCID(-1), fEvtMap(0)
0061 {
0062   fNEvent = 0;
0063   G4double tMin = 1.0 * CLHEP::picosecond;
0064   G4double tMax = 999999 * CLHEP::picosecond;
0065   G4double tLogMin = std::log10(tMin);
0066   G4double tLogMax = std::log10(tMax);
0067   G4int tBins = 50;
0068   for (int i = 0; i <= tBins; i++)
0069     AddTimeToRecord(std::pow(10., tLogMin + i * (tLogMax - tLogMin) / tBins));
0070
0071   fEdep = 0;
0072 }
0073
0074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0075
0076 ScoreSpecies::~ScoreSpecies()
0077 {
0078   ;
0079 }
0080
0081 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0082
0083 G4bool ScoreSpecies::ProcessHits(G4Step* aStep, G4TouchableHistory*)
0084 {
0085   G4double edep = aStep->GetTotalEnergyDeposit();
0086
0087   if (edep == 0.) return FALSE;
0088
0089   edep *= aStep->GetPreStepPoint()->GetWeight();
0090   G4int index = GetIndex(aStep);
0091   fEvtMap->add(index, edep);
0092   fEdep += edep;
0093
0094   return TRUE;
0095 }
0096
0097 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0098
0099 void ScoreSpecies::Initialize(G4HCofThisEvent* HCE)
0100 {
0101   fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName());
0102
0103   if (fHCID < 0) {
0104     fHCID = GetCollectionID(0);
0105   }
0106
0107   HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fEvtMap);
0108 }
0109
0110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0111
0112 void ScoreSpecies::EndOfEvent(G4HCofThisEvent*)
0113 {
0114   if (G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted()) {
0115     fEdep = 0.;
0116     G4MoleculeCounter::Instance()->ResetCounter();
0117     return;
0118   }
0119
0120   auto species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
0121
0122   if (species.get() == 0 || species->size() == 0) {
0123     G4cout << "No molecule recorded, energy deposited= " << G4BestUnit(fEdep, "Energy") << G4endl;
0124     ++fNEvent;
0125     fEdep = 0.;
0126     G4MoleculeCounter::Instance()->ResetCounter();
0127     return;
0128   }
0129
0130   for (auto molecule : *species) {
0131     for (auto time_mol : fTimeToRecord) {
0132       double n_mol = G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule, time_mol);
0133
0134       if (n_mol < 0) {
0135         G4cerr << "N molecules not valid < 0 " << G4endl;
0136         G4Exception("", "N<0", FatalException, "");
0137       }
0138
0139       SpeciesInfo& molInfo = fSpeciesInfoPerTime[time_mol][molecule];
0140       molInfo.fNumber += n_mol;
0141       double gValue = (n_mol / (fEdep / eV)) * 100.;
0142       molInfo.fG += gValue;
0143       molInfo.fG2 += gValue * gValue;
0144     }
0145   }
0146   ++fNEvent;
0147
0148   fEdep = 0.;
0149   G4MoleculeCounter::Instance()->ResetCounter();
0150 }
0151
0152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0153
0154 void ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
0155 {
0156   ScoreSpecies* right =
0157     dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
0158
0159   if (right == 0) {
0160     return;
0161   }
0162   if (right == this) {
0163     return;
0164   }
0165
0166   SpeciesMap::iterator it_map1 = right->fSpeciesInfoPerTime.begin();
0167   SpeciesMap::iterator end_map1 = right->fSpeciesInfoPerTime.end();
0168
0169   for (; it_map1 != end_map1; ++it_map1) {
0170     InnerSpeciesMap& map2 = it_map1->second;
0171     InnerSpeciesMap::iterator it_map2 = map2.begin();
0172     InnerSpeciesMap::iterator end_map2 = map2.end();
0173
0174     for (; it_map2 != end_map2; ++it_map2) {
0175       SpeciesInfo& molInfo = fSpeciesInfoPerTime[it_map1->first][it_map2->first];
0176       molInfo.fNumber += it_map2->second.fNumber;
0177       molInfo.fG += it_map2->second.fG;
0178       molInfo.fG2 += it_map2->second.fG2;
0179     }
0180   }
0181
0182   fNEvent += right->fNEvent;
0183   right->fNEvent = 0;
0184   right->fEdep = 0.;
0185 }
0186
0187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0188
0189 void ScoreSpecies::DrawAll()
0190 {
0191   ;
0192 }
0193
0194 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0195
0196 void ScoreSpecies::PrintAll()
0197 {
0198   G4cout << " MultiFunctionalDet  " << detector->GetName() << G4endl;
0199   G4cout << " PrimitiveScorer " << GetName() << G4endl;
0200   G4cout << " Number of events " << fNEvent << G4endl;
0201   G4cout << " Number of energy deposition recorded " << fEvtMap->entries() << G4endl;
0202
0203   for (auto itr : *fEvtMap->GetMap()) {
0204     G4cout << "  copy no.: " << itr.first << "  energy deposit: " << *(itr.second) / GetUnitValue()
0205            << " [" << GetUnit() << "]" << G4endl;
0206   }
0207 }
0208
0209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0210
0211 void ScoreSpecies::ASCII()
0212 {
0213   std::ofstream out("Species.txt");
0214   if (!out) return;
0215
0216   out << "# Time [ps]  G-value (/100 eV)  RMS   Molecule" << G4endl;
0217
0218   std::map<G4String, std::map<G4double, std::pair<G4double, G4double>>> mol;
0219
0220   for (auto it_map1 : fSpeciesInfoPerTime) {
0221     InnerSpeciesMap& map2 = it_map1.second;
0222     G4double time = it_map1.first / ps;
0223     for (auto it_map2 : map2) {
0224       G4double G = it_map2.second.fG;
0225       G4double G2 = it_map2.second.fG2;
0226       G4double N = fNEvent;
0227       G /= N;
0228       G2 = std::sqrt(N / (N - 1) * (G2 / N - G * G));
0229       mol[it_map2.first->GetName()][time] = std::make_pair(G, G2);
0230     }
0231   }
0232
0233   for (auto it1 : mol)
0234     for (auto it2 : it1.second)
0235       out << std::setw(12) << it2.first << std::setw(12) << it2.second.first << std::setw(12)
0236           << it2.second.second << std::setw(12) << std::setw(12) << it1.first << G4endl;
0237
0238   out.close();
0239 }
0240
0241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0242
0243 void ScoreSpecies::OutputAndClear()
0244 {
0245   if (G4Threading::IsWorkerThread()) return;
0246
0247   //----------------------------------------------------------------------------
0248   // Save results in ASCII format
0249   ASCII();
0250
0251   fNEvent = 0;
0252   fSpeciesInfoPerTime.clear();
0253 }