Geant4  10.03
ScoreSpecies.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708
30 // J. Comput. Phys. 274 (2014) 841-882
31 // The Geant4-DNA web site is available at http://geant4-dna.org
32 //
33 // $Id$
34 // ScoreSpecies.cc
35 //
36 #include "ScoreSpecies.hh"
37 
38 #include "G4UnitsTable.hh"
40 #include <G4Molecule.hh>
41 #include <G4MoleculeCounter.hh>
42 #include "G4Event.hh"
43 #include <G4SystemOfUnits.hh>
44 #include <globals.hh>
45 #include <G4VAnalysisManager.hh>
46 #include <G4RootAnalysisManager.hh>
47 #include <G4XmlAnalysisManager.hh>
48 #include <G4EventManager.hh>
49 #include "g4analysis_defs.hh"
50 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
63 : G4VPrimitiveScorer(name,depth),
64  fEdep(0),
65  fOutputToRoot(true),
66  fOutputToXml(false),
67  fOutputToCsv(false),
68  fHCID(-1),
69  fEvtMap(0)
70 {
71  fNEvent = 0;
79  fEdep = 0;
80 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83 
85 {;}
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 
90 {
91  G4double edep = aStep->GetTotalEnergyDeposit();
92 
93  if ( edep == 0. ) return FALSE;
94 
95  edep *= aStep->GetPreStepPoint()->GetWeight(); // (Particle Weight)
96  G4int index = GetIndex(aStep);
97  fEvtMap->add(index,edep);
98  fEdep+=edep;
99 
100  return TRUE;
101 }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104 
106 {
108  GetName());
109 
110  if(fHCID < 0)
111  {
112  fHCID = GetCollectionID(0);
113  }
114 
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
119 
121 {
122  if(G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted())
123  {
124  fEdep = 0.;
126  return;
127  }
128 
130 
131  if(species.get() == 0 || species->size() == 0)
132  {
133  G4cout << "No molecule recorded, energy deposited= "
134  << G4BestUnit(fEdep, "Energy") << G4endl;
135  ++fNEvent;
136  fEdep = 0.;
138  return;
139  }
140 
141  // G4cout << "ScoreSpecies::EndOfEvent"<<G4endl;
142 // G4cout << "End of event, deposited energy: "
143 // << G4BestUnit(fEdep, "Energy") << G4endl;
144 
145 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
146  int eventID=
148 #endif
149 
150  for(auto molecule: *species)
151  {
152  for(auto time_mol: fTimeToRecord)
153  {
154  double n_mol =
156  time_mol);
157 
158  if(n_mol < 0)
159  {
160  G4cerr << "N molecules not valid < 0 " << G4endl;
161  G4Exception("","N<0",FatalException,"");
162  }
163 
164  SpeciesInfo& molInfo =
165  fSpeciesInfoPerTime[time_mol][molecule];
166  molInfo.fNumber += n_mol;
167  double gValue = (n_mol/(fEdep/eV)) * 100.;
168  molInfo.fG += gValue;
169  molInfo.fG2 += gValue*gValue;
170 
171 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
172  SpeciesInfoSOA& molInfoPerEvent =
173  fSpeciesInfoPerEvent[time_mol][molecule];
174  molInfoPerEvent.fNumber.push_back(n_mol);
175  molInfoPerEvent.fG.push_back(gValue);
176  molInfoPerEvent.fG2.push_back(gValue*gValue);
177  molInfoPerEvent.fEventID.push_back(eventID);
178 #endif
179  // G4cout << "In Save molucule: fNumber " << molInfo.fNumber
180  // << " fG " << molInfo.fG
181  // << " fEdep " << fEdep/eV
182  // << G4endl;
183  }
184  }
185 
186  ++fNEvent;
187 
188 // G4cout << "End of event " << fNEvent
189 // << ", energy deposited=" << G4BestUnit(fEdep, "Energy") << G4endl;
190 
191  fEdep = 0.;
193 }
194 
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
196 
197 void
199 {
201  dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
202 
203  if(right == 0)
204  {
205  return;
206  }
207  if(right == this)
208  {
209  return;
210  }
211 
212  // G4cout<<"ScoreSpecies::AbsorbResultsFromWorkerScorer"<<G4endl;
213  {
214  SpeciesMap::iterator it_map1 = right->fSpeciesInfoPerTime.begin();
215  SpeciesMap::iterator end_map1 = right->fSpeciesInfoPerTime.end();
216 
217  for(; it_map1 != end_map1; ++it_map1)
218  {
219  InnerSpeciesMap& map2 = it_map1->second;
220  InnerSpeciesMap::iterator it_map2 = map2.begin();
221  InnerSpeciesMap::iterator end_map2 = map2.end();
222 
223  for(; it_map2 != end_map2; ++it_map2)
224  {
225  SpeciesInfo& molInfo =
226  fSpeciesInfoPerTime[it_map1->first][it_map2->first] ;
227  molInfo.fNumber += it_map2->second.fNumber;
228  molInfo.fG += it_map2->second.fG;
229  molInfo.fG2 += it_map2->second.fG2;
230 
231  // G4cout << "In AbsorbeResultsFromWorkerScorer: fNumber "
232  // << molInfo.fNumber
233  // << " fG "
234  // << molInfo.fG
235  // << G4endl;
236  }
237  }
238  }
239  //---------------------------------------------------------
240 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
241  {
242  SpeciesMapPerEvent::iterator it_map1 = right->fSpeciesInfoPerEvent.begin();
243  SpeciesMapPerEvent::iterator end_map1 = right->fSpeciesInfoPerEvent.end();
244 
245  for(; it_map1 != end_map1; ++it_map1)
246  {
247  auto& map2 = it_map1->second;
248  InnerSpeciesMapPerEvent::iterator it_map2 = map2.begin();
249  InnerSpeciesMapPerEvent::iterator end_map2 = map2.end();
250 
251  for(; it_map2 != end_map2; ++it_map2)
252  {
253  SpeciesInfoSOA& molInfo =
254  fSpeciesInfoPerEvent[it_map1->first][it_map2->first] ;
255  molInfo.fNumber.insert(molInfo.fNumber.end(),
256  it_map2->second.fNumber.begin(),
257  it_map2->second.fNumber.end());
258  molInfo.fG.insert(molInfo.fG.end(),
259  it_map2->second.fG.begin(),
260  it_map2->second.fG.end());
261  molInfo.fG2.insert(molInfo.fG2.end(),
262  it_map2->second.fG2.begin(),
263  it_map2->second.fG2.end());
264  molInfo.fEventID.insert(molInfo.fEventID.end(),
265  it_map2->second.fEventID.begin(),
266  it_map2->second.fEventID.end());
267  // G4cout << "In AbsorbeResultsFromWorkerScorer: fNumber "
268  // << molInfo.fNumber
269  // << " fG "
270  // << molInfo.fG
271  // << G4endl;
272  }
273  }
274  right->fSpeciesInfoPerEvent.clear();
275  }
276 #endif
277 
278  fNEvent += right->fNEvent;
279  right->fNEvent = 0;
280  right->fEdep = 0.;
281 }
282 
283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
284 
286 {
287  fEvtMap->clear();
288  fNEvent = 0;
289  fEdep = 0;
290  fSpeciesInfoPerTime.clear();
291 }
292 
293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
294 
296 {;}
297 
298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
299 
301 {
302  G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
303  G4cout << " PrimitiveScorer " << GetName() << G4endl;
304  G4cout << " Number of events " << fNEvent << G4endl;
305  G4cout << " Number of energy deposition recorded "
306  << fEvtMap->entries() << G4endl;
307 
308  for(auto itr : *fEvtMap->GetMap()) {
309  G4cout << " copy no.: " << itr.first
310  << " energy deposit: "
311  << *(itr.second)/GetUnitValue()
312  << " [" << GetUnit()<<"]"
313  << G4endl;
314  }
315 }
316 
317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
318 
320 {
321  std::ofstream out("Species.Txt");
322  if(!out) return;
323 
324  out << "Time is in ns" << G4endl;
325 
326  for(auto it_map1: fSpeciesInfoPerTime)
327  {
328  InnerSpeciesMap& map2 = it_map1.second;
329 
330  out << it_map1.first << G4endl;
331 
332  for(auto it_map2: map2)
333  {
334  out << it_map2.first->GetName()<< " "
335  << it_map2.second.fNumber << G4endl;
336  }
337  }
338 
339  out.close();
340 }
341 
342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
343 
345 {
346  if(G4Threading::IsWorkerThread()) return;
347 
348  //----------------------------------------------------------------------------
349  // Save results
350 
351  G4VAnalysisManager* analysisManager(0);
352 
353  if(fOutputToCsv)
354  {
355  analysisManager = G4CsvAnalysisManager::Instance(); // TODO?
356  // this->ASCII(); // useful ?
357  }
358  else if (fOutputToRoot)
359  {
360  analysisManager = G4Root::G4AnalysisManager::Instance();
361  }
362  else if(fOutputToXml)
363  {
364  analysisManager = G4Xml::G4AnalysisManager::Instance();
365 
366  }
367  if(analysisManager)
368  {
369  this->WriteWithAnalysisManager(analysisManager);
370  }
371 
372  fNEvent = 0;
373  fSpeciesInfoPerTime.clear();
374 }
375 
376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
377 
378 void
380 {
381  // G4cout << "ScoreSpecies::WriteWithAnalysisManager" << G4endl;
382  analysisManager->OpenFile("Species.root");
383  int fNtupleID = analysisManager->CreateNtuple("species", "species");
384  analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
385  analysisManager->CreateNtupleIColumn(fNtupleID, "number");
386  analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
387  analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
388  analysisManager->CreateNtupleDColumn(fNtupleID, "time");
389  analysisManager->CreateNtupleDColumn(fNtupleID, "sumG");
390  analysisManager->CreateNtupleDColumn(fNtupleID, "sumG2");
391  analysisManager->FinishNtuple(fNtupleID);
392 
393  for(auto it_map1: fSpeciesInfoPerTime)
394  {
395  InnerSpeciesMap& map2 = it_map1.second;
396 
397  for(auto it_map2 : map2)
398  {
399  double time = it_map1.first;
400  const Species& species = it_map2.first;
401  const G4String& name = species->GetName();
402  int molID = it_map2.first->GetMoleculeID();
403  int number = it_map2.second.fNumber;
404  double G = it_map2.second.fG;
405  double G2 = it_map2.second.fG2;
406 
407  analysisManager->FillNtupleIColumn(fNtupleID, 0, molID); // MolID
408  analysisManager->FillNtupleIColumn(fNtupleID, 1, number); // Number
409  analysisManager->FillNtupleIColumn(fNtupleID,
410  2, fNEvent); // Total nb events
411  analysisManager->FillNtupleSColumn(fNtupleID, 3, name); // molName
412  analysisManager->FillNtupleDColumn(fNtupleID, 4, time); // time
413  analysisManager->FillNtupleDColumn(fNtupleID, 5, G); // G
414  analysisManager->FillNtupleDColumn(fNtupleID, 6, G2); // G2
415  analysisManager->AddNtupleRow(fNtupleID);
416  }
417  }
418 
419  //----------------------------------------------------------------------------
420 
421 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
422  fNtupleID = analysisManager->CreateNtuple("species_all", "species_all");
423  analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
424  analysisManager->CreateNtupleIColumn(fNtupleID, "number");
425  analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
426  analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
427  analysisManager->CreateNtupleDColumn(fNtupleID, "time");
428  analysisManager->CreateNtupleDColumn(fNtupleID, "G");
429  analysisManager->CreateNtupleDColumn(fNtupleID, "G2");
430  analysisManager->CreateNtupleIColumn(fNtupleID, "eventID");
431  analysisManager->FinishNtuple(fNtupleID);
432 
433  for(auto it_map1: fSpeciesInfoPerEvent)
434  {
435  InnerSpeciesMapPerEvent& map2 = it_map1.second;
436 
437  for(auto it_map2 : map2)
438  {
439  double time = it_map1.first;
440  const Species& species = it_map2.first;
441  const G4String& name = species->GetName();
442  int molID = it_map2.first->GetMoleculeID();
443 
444  size_t nG = it_map2.second.fG.size();
445 
446  for(size_t i=0; i<nG;++i){
447  int number = it_map2.second.fNumber[i];
448  double G = it_map2.second.fG[i];
449  double G2 = it_map2.second.fG2[i];
450  int eventID = it_map2.second.fEventID[i];
451 
452  analysisManager->FillNtupleIColumn(fNtupleID, 0, molID); // MolID
453  analysisManager->FillNtupleIColumn(fNtupleID, 1, number); // Number
454  analysisManager->FillNtupleIColumn(fNtupleID,
455  2, fNEvent); // Total nb events
456  analysisManager->FillNtupleSColumn(fNtupleID, 3, name); // molName
457  analysisManager->FillNtupleDColumn(fNtupleID, 4, time); // time
458  analysisManager->FillNtupleDColumn(fNtupleID, 5, G); // G
459  analysisManager->FillNtupleDColumn(fNtupleID, 6, G2); // G2
460  analysisManager->FillNtupleIColumn(fNtupleID, 7, eventID); // EventID
461  analysisManager->AddNtupleRow(fNtupleID);
462  }
463  }
464  }
465 #endif
466 
467  analysisManager->Write();
468  analysisManager->CloseFile();
469 }
The pointer G4MolecularConfiguration will be shared by all the molecules having the same molecule def...
void WriteWithAnalysisManager(G4VAnalysisManager *)
Write results to whatever chosen file format.
void clear()
Definition: G4THitsMap.hh:267
virtual G4bool ProcessHits(G4Step *, G4TouchableHistory *)
Definition: ScoreSpecies.cc:89
const G4String & GetName() const
Returns the name of the molecule.
G4double GetWeight() const
G4bool fOutputToXml
G4int CreateNtupleIColumn(const G4String &name)
G4bool FillNtupleSColumn(G4int id, const G4String &value)
G4String GetName() const
G4int first(char) const
void ResetCounter() override
virtual void PrintAll()
G4bool fOutputToRoot
std::set< G4double > fTimeToRecord
G4MultiFunctionalDetector * GetMultiFunctionalDetector() const
virtual void Initialize(G4HCofThisEvent *)
G4int CreateNtuple(const G4String &name, const G4String &title)
static G4XmlAnalysisManager * Instance()
const char * name(G4int ptype)
void ASCII()
Write results to an text file.
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
virtual void DrawAll()
G4int CreateNtupleSColumn(const G4String &name)
int G4int
Definition: G4Types.hh:78
G4bool OpenFile(const G4String &fileName="")
static G4CsvAnalysisManager * Instance()
static constexpr double picosecond
Definition: G4SIunits.hh:161
int GetNMoleculesAtTime(G4MolecularConfiguration *molecule, double time)
SpeciesMap fSpeciesInfoPerTime
G4int GetEventID() const
Definition: G4Event.hh:151
G4StepPoint * GetPreStepPoint() const
RecordedMolecules GetRecordedMolecules()
G4bool FillNtupleIColumn(G4int id, G4int value)
This is a primitive scorer class for molecular species.
Definition: ScoreSpecies.hh:56
G4GLOB_DLL std::ostream G4cout
virtual void clear()
static G4MoleculeCounter * Instance()
G4bool FillNtupleDColumn(G4int id, G4double value)
G4double GetUnitValue() const
bool G4bool
Definition: G4Types.hh:79
#define FALSE
Definition: globals.hh:52
virtual G4int GetIndex(G4Step *)
static constexpr double eV
Definition: G4SIunits.hh:215
#define TRUE
Definition: globals.hh:55
void AddTimeToRecord(double time)
Add a time at which the number of species should be recorded.
Definition: ScoreSpecies.hh:65
Definition: G4Step.hh:76
virtual ~ScoreSpecies()
Definition: ScoreSpecies.cc:84
G4int GetCollectionID(G4int)
virtual void EndOfEvent(G4HCofThisEvent *)
G4double GetTotalEnergyDeposit() const
void AddHitsCollection(G4int HCID, G4VHitsCollection *aHC)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4bool IsWorkerThread()
Definition: G4Threading.cc:145
G4int entries() const
Definition: G4THitsMap.hh:200
virtual void OutputAndClear()
G4int add(const G4int &key, T *&aHit) const
Definition: G4THitsMap.hh:106
static G4EventManager * GetEventManager()
static G4RootAnalysisManager * Instance()
std::map< G4int, T * > * GetMap() const
Definition: G4THitsMap.hh:99
std::map< Species, SpeciesInfo > InnerSpeciesMap
#define G4endl
Definition: G4ios.hh:61
virtual void AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer *)
Method used in multithreading mode in order to merge the results.
G4MultiFunctionalDetector * detector
G4int CreateNtupleDColumn(const G4String &name)
double G4double
Definition: G4Types.hh:76
const G4Event * GetConstCurrentEvent()
const G4String & GetUnit() const
G4bool fOutputToCsv
G4THitsMap< G4double > * fEvtMap
G4GLOB_DLL std::ostream G4cerr
ScoreSpecies(G4String name, G4int depth=0)
Definition: ScoreSpecies.cc:62
std::vector< double > fG
Definition: plotG.cc:149