ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ScoreSpecies.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ScoreSpecies.cc
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 // ScoreSpecies.cc
34 //
35 #include "ScoreSpecies.hh"
36 
37 #include "G4UnitsTable.hh"
38 #include <G4MolecularConfiguration.hh>
39 #include <G4MoleculeCounter.hh>
40 #include "G4Event.hh"
41 #include <G4SystemOfUnits.hh>
42 #include <globals.hh>
43 #include <G4EventManager.hh>
44 #include "g4analysis.hh"
45 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
57 ScoreSpecies::ScoreSpecies(G4String name, G4int depth)
58 : G4VPrimitiveScorer(name,depth),
59  fEdep(0),
60  fOutputToRoot(true),
61  fOutputToXml(false),
62  fOutputToCsv(false),
63  fHCID(-1),
64  fEvtMap(0)
65 {
66  fNEvent = 0;
67  AddTimeToRecord(1 * CLHEP::picosecond);
68  AddTimeToRecord(10 * CLHEP::picosecond);
69  AddTimeToRecord(100 * CLHEP::picosecond);
70  AddTimeToRecord(1000 * CLHEP::picosecond);
71  AddTimeToRecord(10000 * CLHEP::picosecond);
72  AddTimeToRecord(100000 * CLHEP::picosecond);
73  AddTimeToRecord(999999 * CLHEP::picosecond);
74  fEdep = 0;
75 }
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78 
79 ScoreSpecies::~ScoreSpecies()
80 {;}
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83 
84 G4bool ScoreSpecies::ProcessHits(G4Step* aStep,G4TouchableHistory*)
85 {
86  G4double edep = aStep->GetTotalEnergyDeposit();
87 
88  if ( edep == 0. ) return FALSE;
89 
90  edep *= aStep->GetPreStepPoint()->GetWeight(); // (Particle Weight)
91  G4int index = GetIndex(aStep);
92  fEvtMap->add(index,edep);
93  fEdep+=edep;
94 
95  return TRUE;
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99 
100 void ScoreSpecies::Initialize(G4HCofThisEvent* HCE)
101 {
102  fEvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(),
103  GetName());
104 
105  if(fHCID < 0)
106  {
107  fHCID = GetCollectionID(0);
108  }
109 
110  HCE->AddHitsCollection(fHCID, (G4VHitsCollection*)fEvtMap);
111 }
112 
113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114 
115 void ScoreSpecies::EndOfEvent(G4HCofThisEvent*)
116 {
117  if(G4EventManager::GetEventManager()->GetConstCurrentEvent()->IsAborted())
118  {
119  fEdep = 0.;
120  G4MoleculeCounter::Instance()->ResetCounter();
121  return;
122  }
123 
124  auto species = G4MoleculeCounter::Instance()->GetRecordedMolecules();
125 
126  if(species.get() == 0 || species->size() == 0)
127  {
128  G4cout << "No molecule recorded, energy deposited= "
129  << G4BestUnit(fEdep, "Energy") << G4endl;
130  ++fNEvent;
131  fEdep = 0.;
132  G4MoleculeCounter::Instance()->ResetCounter();
133  return;
134  }
135 
136  // G4cout << "ScoreSpecies::EndOfEvent"<<G4endl;
137 // G4cout << "End of event, deposited energy: "
138 // << G4BestUnit(fEdep, "Energy") << G4endl;
139 
140 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
141  int eventID=
142  G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetEventID();
143 #endif
144 
145  for(auto molecule: *species)
146  {
147  for(auto time_mol: fTimeToRecord)
148  {
149  double n_mol =
150  G4MoleculeCounter::Instance()->GetNMoleculesAtTime(molecule,
151  time_mol);
152 
153  if(n_mol < 0)
154  {
155  G4cerr << "N molecules not valid < 0 " << G4endl;
156  G4Exception("","N<0",FatalException,"");
157  }
158 
159  SpeciesInfo& molInfo = fSpeciesInfoPerTime[time_mol][molecule];
160  molInfo.fNumber += n_mol;
161  double gValue = (n_mol/(fEdep/eV)) * 100.;
162  molInfo.fG += gValue;
163  molInfo.fG2 += gValue*gValue;
164 
165 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
166  SpeciesInfoSOA& molInfoPerEvent =
167  fSpeciesInfoPerEvent[time_mol][molecule];
168  molInfoPerEvent.fNumber.push_back(n_mol);
169  molInfoPerEvent.fG.push_back(gValue);
170  molInfoPerEvent.fG2.push_back(gValue*gValue);
171  molInfoPerEvent.fEventID.push_back(eventID);
172 #endif
173  // G4cout << "In Save molucule: fNumber " << molInfo.fNumber
174  // << " fG " << molInfo.fG
175  // << " fEdep " << fEdep/eV
176  // << G4endl;
177  }
178  }
179 
180  ++fNEvent;
181 
182 // G4cout << "End of event " << fNEvent
183 // << ", energy deposited=" << G4BestUnit(fEdep, "Energy") << G4endl;
184 
185  fEdep = 0.;
186  G4MoleculeCounter::Instance()->ResetCounter();
187 }
188 
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
190 
191 void
192 ScoreSpecies::AbsorbResultsFromWorkerScorer(G4VPrimitiveScorer* workerScorer)
193 {
194  ScoreSpecies* right =
195  dynamic_cast<ScoreSpecies*>(dynamic_cast<G4VPrimitiveScorer*>(workerScorer));
196 
197  if(right == 0)
198  {
199  return;
200  }
201  if(right == this)
202  {
203  return;
204  }
205 
206  // G4cout<<"ScoreSpecies::AbsorbResultsFromWorkerScorer"<<G4endl;
207  {
208  SpeciesMap::iterator it_map1 = right->fSpeciesInfoPerTime.begin();
209  SpeciesMap::iterator end_map1 = right->fSpeciesInfoPerTime.end();
210 
211  for(; it_map1 != end_map1; ++it_map1)
212  {
213  InnerSpeciesMap& map2 = it_map1->second;
214  InnerSpeciesMap::iterator it_map2 = map2.begin();
215  InnerSpeciesMap::iterator end_map2 = map2.end();
216 
217  for(; it_map2 != end_map2; ++it_map2)
218  {
219  SpeciesInfo& molInfo =
220  fSpeciesInfoPerTime[it_map1->first][it_map2->first] ;
221  molInfo.fNumber += it_map2->second.fNumber;
222  molInfo.fG += it_map2->second.fG;
223  molInfo.fG2 += it_map2->second.fG2;
224 
225  // G4cout << "In AbsorbeResultsFromWorkerScorer: fNumber "
226  // << molInfo.fNumber
227  // << " fG "
228  // << molInfo.fG
229  // << G4endl;
230  }
231  }
232  }
233  //---------------------------------------------------------
234 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
235  {
236  SpeciesMapPerEvent::iterator it_map1 = right->fSpeciesInfoPerEvent.begin();
237  SpeciesMapPerEvent::iterator end_map1 = right->fSpeciesInfoPerEvent.end();
238 
239  for(; it_map1 != end_map1; ++it_map1)
240  {
241  auto& map2 = it_map1->second;
242  InnerSpeciesMapPerEvent::iterator it_map2 = map2.begin();
243  InnerSpeciesMapPerEvent::iterator end_map2 = map2.end();
244 
245  for(; it_map2 != end_map2; ++it_map2)
246  {
247  SpeciesInfoSOA& molInfo =
248  fSpeciesInfoPerEvent[it_map1->first][it_map2->first] ;
249  molInfo.fNumber.insert(molInfo.fNumber.end(),
250  it_map2->second.fNumber.begin(),
251  it_map2->second.fNumber.end());
252  molInfo.fG.insert(molInfo.fG.end(),
253  it_map2->second.fG.begin(),
254  it_map2->second.fG.end());
255  molInfo.fG2.insert(molInfo.fG2.end(),
256  it_map2->second.fG2.begin(),
257  it_map2->second.fG2.end());
258  molInfo.fEventID.insert(molInfo.fEventID.end(),
259  it_map2->second.fEventID.begin(),
260  it_map2->second.fEventID.end());
261  // G4cout << "In AbsorbeResultsFromWorkerScorer: fNumber "
262  // << molInfo.fNumber
263  // << " fG "
264  // << molInfo.fG
265  // << G4endl;
266  }
267  }
268  right->fSpeciesInfoPerEvent.clear();
269  }
270 #endif
271 
272  fNEvent += right->fNEvent;
273  right->fNEvent = 0;
274  right->fEdep = 0.;
275 }
276 
277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
278 
279 void ScoreSpecies::clear()
280 {
281  fEvtMap->clear();
282  fNEvent = 0;
283  fEdep = 0;
284  fSpeciesInfoPerTime.clear();
285 }
286 
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 
289 void ScoreSpecies::DrawAll()
290 {;}
291 
292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
293 
294 void ScoreSpecies::PrintAll()
295 {
296  G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl;
297  G4cout << " PrimitiveScorer " << GetName() << G4endl;
298  G4cout << " Number of events " << fNEvent << G4endl;
299  G4cout << " Number of energy deposition recorded "
300  << fEvtMap->entries() << G4endl;
301 
302  for(auto itr : *fEvtMap->GetMap()) {
303  G4cout << " copy no.: " << itr.first
304  << " energy deposit: "
305  << *(itr.second)/GetUnitValue()
306  << " [" << GetUnit()<<"]"
307  << G4endl;
308  }
309 }
310 
311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
312 
313 void ScoreSpecies::ASCII()
314 {
315  std::ofstream out("Species.Txt");
316  if(!out) return;
317 
318  out << "Time is in ns" << G4endl;
319 
320  for(auto it_map1: fSpeciesInfoPerTime)
321  {
322  InnerSpeciesMap& map2 = it_map1.second;
323 
324  out << it_map1.first << G4endl;
325 
326  for(auto it_map2: map2)
327  {
328  out << it_map2.first->GetName()<< " "
329  << it_map2.second.fNumber << G4endl;
330  }
331  }
332 
333  out.close();
334 }
335 
336 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
337 
338 void ScoreSpecies::OutputAndClear()
339 {
340  if(G4Threading::IsWorkerThread()) return;
341 
342  //----------------------------------------------------------------------------
343  // Save results
344 
345  G4VAnalysisManager* analysisManager(0);
346 
347  if(fOutputToCsv)
348  {
349  analysisManager = G4Analysis::ManagerInstance("csv");
350  // this->ASCII(); // useful ?
351  }
352  else if (fOutputToRoot)
353  {
354  analysisManager = G4Analysis::ManagerInstance("root");
355  }
356  else if(fOutputToXml)
357  {
358  analysisManager = G4Analysis::ManagerInstance("xml");
359 
360  }
361  if(analysisManager)
362  {
363  this->WriteWithAnalysisManager(analysisManager);
364  }
365 
366  fNEvent = 0;
367  fSpeciesInfoPerTime.clear();
368 }
369 
370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
371 
372 void
373 ScoreSpecies::WriteWithAnalysisManager(G4VAnalysisManager* analysisManager)
374 {
375  // G4cout << "ScoreSpecies::WriteWithAnalysisManager" << G4endl;
376  analysisManager->OpenFile("Species.root");
377  int fNtupleID = analysisManager->CreateNtuple("species", "species");
378  analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
379  analysisManager->CreateNtupleIColumn(fNtupleID, "number");
380  analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
381  analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
382  analysisManager->CreateNtupleDColumn(fNtupleID, "time");
383  analysisManager->CreateNtupleDColumn(fNtupleID, "sumG");
384  analysisManager->CreateNtupleDColumn(fNtupleID, "sumG2");
385  analysisManager->FinishNtuple(fNtupleID);
386 
387  for(auto it_map1: fSpeciesInfoPerTime)
388  {
389  InnerSpeciesMap& map2 = it_map1.second;
390 
391  for(auto it_map2 : map2)
392  {
393  double time = it_map1.first;
394  auto species = it_map2.first;
395  const G4String& name = species->GetName();
396  int molID = it_map2.first->GetMoleculeID();
397  int number = it_map2.second.fNumber;
398  double G = it_map2.second.fG;
399  double G2 = it_map2.second.fG2;
400 
401  analysisManager->FillNtupleIColumn(fNtupleID, 0, molID); // MolID
402  analysisManager->FillNtupleIColumn(fNtupleID, 1, number); // Number
403  analysisManager->FillNtupleIColumn(fNtupleID,
404  2, fNEvent); // Total nb events
405  analysisManager->FillNtupleSColumn(fNtupleID, 3, name); // molName
406  analysisManager->FillNtupleDColumn(fNtupleID, 4, time); // time
407  analysisManager->FillNtupleDColumn(fNtupleID, 5, G); // G
408  analysisManager->FillNtupleDColumn(fNtupleID, 6, G2); // G2
409  analysisManager->AddNtupleRow(fNtupleID);
410  }
411  }
412 
413  //----------------------------------------------------------------------------
414 
415 #ifdef _ScoreSpecies_FOR_ALL_EVENTS
416  fNtupleID = analysisManager->CreateNtuple("species_all", "species_all");
417  analysisManager->CreateNtupleIColumn(fNtupleID, "speciesID");
418  analysisManager->CreateNtupleIColumn(fNtupleID, "number");
419  analysisManager->CreateNtupleIColumn(fNtupleID, "nEvent");
420  analysisManager->CreateNtupleSColumn(fNtupleID, "speciesName");
421  analysisManager->CreateNtupleDColumn(fNtupleID, "time");
422  analysisManager->CreateNtupleDColumn(fNtupleID, "G");
423  analysisManager->CreateNtupleDColumn(fNtupleID, "G2");
424  analysisManager->CreateNtupleIColumn(fNtupleID, "eventID");
425  analysisManager->FinishNtuple(fNtupleID);
426 
427  for(auto it_map1: fSpeciesInfoPerEvent)
428  {
429  InnerSpeciesMapPerEvent& map2 = it_map1.second;
430 
431  for(auto it_map2 : map2)
432  {
433  double time = it_map1.first;
434  const Species& species = it_map2.first;
435  const G4String& name = species->GetName();
436  int molID = it_map2.first->GetMoleculeID();
437 
438  size_t nG = it_map2.second.fG.size();
439 
440  for(size_t i=0; i<nG;++i){
441  int number = it_map2.second.fNumber[i];
442  double G = it_map2.second.fG[i];
443  double G2 = it_map2.second.fG2[i];
444  int eventID = it_map2.second.fEventID[i];
445 
446  analysisManager->FillNtupleIColumn(fNtupleID, 0, molID); // MolID
447  analysisManager->FillNtupleIColumn(fNtupleID, 1, number); // Number
448  analysisManager->FillNtupleIColumn(fNtupleID,
449  2, fNEvent); // Total nb events
450  analysisManager->FillNtupleSColumn(fNtupleID, 3, name); // molName
451  analysisManager->FillNtupleDColumn(fNtupleID, 4, time); // time
452  analysisManager->FillNtupleDColumn(fNtupleID, 5, G); // G
453  analysisManager->FillNtupleDColumn(fNtupleID, 6, G2); // G2
454  analysisManager->FillNtupleIColumn(fNtupleID, 7, eventID); // EventID
455  analysisManager->AddNtupleRow(fNtupleID);
456  }
457  }
458  }
459 #endif
460 
461  analysisManager->Write();
462  analysisManager->CloseFile();
463 }