EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
eASTRunAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file eASTRunAction.cc
1 // ********************************************************************
2 //
3 // eASTRunAction.cc
4 // eAST Run Action class that takes care of defining and handling
5 // histograms and n-tuple.
6 // Filling histograms is taken care by eASTRun class.
7 //
8 // History
9 // May 8th, 2021 : first implementation
10 // Dec.22.2021 : migration to Geant4 version 11.0 - Makoto Asai (JLab)
11 //
12 // ********************************************************************
13 
14 #include "eASTRunAction.hh"
15 #include "eASTRunActionMessenger.hh"
16 #include "eASTAnalysis.hh"
17 #include "eASTUserActionDispatcher.hh"
18 
19 #include "G4Run.hh"
20 #include "G4RunManager.hh"
21 #include "G4UnitsTable.hh"
22 #include "G4SystemOfUnits.hh"
23 #include "G4Version.hh"
24 
25 eASTRunAction::eASTRunAction()
26 {
27  messenger = new eASTRunActionMessenger(this);
28  auto analysisManager = G4AnalysisManager::Instance();
29 #if G4VERSION_NUMBER >= 1100
30  analysisManager->SetDefaultFileType("csv");
31 #endif
32  G4cout << "Using " << analysisManager->GetType() << G4endl;
33 
34  analysisManager->SetVerboseLevel(verbose);
35  //analysisManager->SetNtupleMerging(true);
36 }
37 
38 eASTRunAction::~eASTRunAction()
39 {
40 #if G4VERSION_NUMBER < 1100
41  delete G4AnalysisManager::Instance();
42 #endif
43  delete messenger;
44 }
45 
46 void eASTRunAction::BeginOfRunAction(const G4Run* run)
47 {
48  // Open an output file
49  //
50  OpenFile();
51 
52  // Define nTuple column if needed
53  //
54  DefineNTColumn();
55 
56  // Invoke component-specific run actions
57  //
58  auto runActions = eASTUserActionDispatcher::Instance()->GetRunActions();
59  if((runActions!=nullptr) && !(runActions->empty()))
60  {
61  for(auto ua : *runActions)
62  { ua.second->BeginOfRunAction(run); }
63  }
64 }
65 
66 void eASTRunAction::OpenFile()
67 {
68  if(!fileOpen)
69  {
70  auto analysisManager = G4AnalysisManager::Instance();
71  analysisManager->OpenFile(fileName);
72  if(verbose>0) G4cout << "eASTRunAction::BeginOfRunAction ### <" << fileName << "> is opened." << G4endl;
73  fileOpen = true;
74  }
75 }
76 
77 void eASTRunAction::EndOfRunAction(const G4Run* run)
78 {
79  if(!ifCarry) Flush();
80 
81  // Invoke component-specific run actions
82  //
83  auto runActions = eASTUserActionDispatcher::Instance()->GetRunActions();
84  if((runActions!=nullptr) && !(runActions->empty()))
85  {
86  for(auto ua : *runActions)
87  { ua.second->EndOfRunAction(run); }
88  }
89 }
90 
91 void eASTRunAction::Flush()
92 {
93  auto analysisManager = G4AnalysisManager::Instance();
94 
95  analysisManager->Write();
96  analysisManager->CloseFile();
97  if(verbose>0) G4cout << "eASTRunAction::Flush ### <" << fileName << "> is closed." << G4endl;
98 
99  fileOpen = false;
100 
101  if(IsMaster()) MergeNtuple();
102 }
103 
104 void eASTRunAction::SetVerbose(G4int val)
105 {
106  verbose = val;
107  auto analysisManager = G4AnalysisManager::Instance();
108  analysisManager->SetVerboseLevel(verbose);
109 }
110 
111 void eASTRunAction::ListHistograms()
112 {
113  G4cout << "################## registered histograms/plots" << G4endl;
114  G4cout << "id\thistID\thistType\tdetName-X\tpsName-X\tcollID-X\tcopyNo-X\tdetName-Y\tpsName-Y\tcollID-Y\tcopyNo-Y" << G4endl;
115  for(auto itr : IDMap)
116  {
117  G4cout << itr.first << "\t" << itr.second->histID << "\t";
118  if(itr.second->histType==1) // 1D histogram
119  { G4cout << "1-D hist\t" << itr.second->meshName << "\t" << itr.second->primName << "\t" << itr.second->collID << "\t" << itr.second->idx; }
120  else if(itr.second->histType==2) // 1D profile
121  { G4cout << "1-D prof\t" << itr.second->meshName << "\t" << itr.second->primName << "\t" << itr.second->collID; }
122  G4cout << G4endl;
123  }
124 }
125 
126 G4bool eASTRunAction::Open(G4int id)
127 {
128  auto hItr = IDMap.find(id);
129  return (hItr!=IDMap.end());
130 }
131 
132 #include "G4SDManager.hh"
133 using namespace G4Analysis;
134 
135 G4bool eASTRunAction::SetAllPlotting(G4bool val)
136 {
137  G4bool valid = true;
138  for(auto hItr : IDMap)
139  {
140  valid = SetPlotting(hItr.first,val);
141  if(!valid) break;
142  }
143  return valid;
144 }
145 
146 G4bool eASTRunAction::SetPlotting(G4int id,G4bool val)
147 {
148  auto hItr = IDMap.find(id);
149  if(hItr==IDMap.end()) return false;
150  auto ht = hItr->second;
151  auto hTyp = ht->histType;
152  auto analysisManager = G4AnalysisManager::Instance();
153  if(hTyp==1) // 1D histogram
154  { analysisManager->SetH1Plotting(ht->histID,val); }
155  else if(hTyp==2) // 1D profile
156  { analysisManager->SetP1Plotting(ht->histID,val); }
157  else
158  { return false; }
159  return true;
160 }
161 
162 // ------------- 1D histogram
163 
164 G4int eASTRunAction::Create1D(G4String& mName,G4String& pName,G4int cn)
165 {
166  G4String collName = mName;
167  collName += "/";
168  collName += pName;
169  auto cID = G4SDManager::GetSDMpointer()->GetCollectionID(collName);
170  if(cID<0) return cID;
171 
172  G4int id = (cID+id_offset)*id_factor+cn+1;
173  auto histoTypeItr = IDMap.find(id);
174  if(histoTypeItr!=IDMap.end()) return false;
175  if(verbose) G4cout << "eASTRunAction::Create1D for <" << collName
176  << ", copyNo=" << cn << "> is registered for hitCollectionID "
177  << cID << G4endl;
178 
179  auto histTyp = new eASTHistoType;
180  histTyp->collID = cID;
181  histTyp->histType = 1; // 1D histogram
182  histTyp->meshName = mName;
183  histTyp->primName = pName;
184  histTyp->idx = cn;
185  IDMap[id] = histTyp;
186  return id;
187 }
188 
189 #include "G4SDManager.hh"
190 #include "G4ScoringManager.hh"
191 #include "G4VScoringMesh.hh"
192 #include "G4VPrimitivePlotter.hh"
193 G4int eASTRunAction::Create1DForPlotter(G4String& mName,G4String& pName,G4bool /*wgt*/)
194 {
195  using MeshShape = G4VScoringMesh::MeshShape;
196 
197  G4String collName = mName;
198  collName += "/";
199  collName += pName;
200  auto cID = G4SDManager::GetSDMpointer()->GetCollectionID(collName);
201  if(cID<0) return cID;
202 
203  auto sm = G4ScoringManager::GetScoringManagerIfExist();
204  assert(sm!=nullptr);
205  auto mesh = sm->FindMesh(mName);
206  if(mesh==nullptr)
207  { return -2; }
208  auto shape = mesh->GetShape();
209  if(shape!=MeshShape::realWorldLogVol && shape!=MeshShape::probe)
210  { return -3; }
211  G4int nBin[3];
212  mesh->GetNumberOfSegments(nBin);
213 
214  auto prim = mesh->GetPrimitiveScorer(pName);
215  if(prim==nullptr)
216  { return -3; }
217  auto pp = dynamic_cast<G4VPrimitivePlotter*>(prim);
218  if(pp==nullptr)
219  { return -4; }
220 
221  G4int id0 = (cID+id_offset)*id_factor+1;
222  for(G4int cn=0; cn<nBin[0]; cn++)
223  {
224  G4int id = id0+cn;
225  auto histoTypeItr = IDMap.find(id);
226  if(histoTypeItr!=IDMap.end())
227  { return -5; }
228  if(verbose) G4cout << "eASTRunAction::Create1D for <" << collName
229  << ", copyNo=" << cn << "> is registered for hitCollectionID "
230  << cID << G4endl;
231 
232  auto histTyp = new eASTHistoType;
233  histTyp->collID = cID;
234  histTyp->histType = 1; // 1D histogram
235  histTyp->histDup = nBin[0];
236  histTyp->meshName = mName;
237  histTyp->primName = pName;
238  histTyp->idx = cn;
239  histTyp->pplotter = pp;
240  IDMap[id] = histTyp;
241  }
242  return id0;
243 }
244 
245 #include "G4UIcommand.hh"
246 G4bool eASTRunAction::Set1D(G4int id0,G4int nBin,G4double valMin,G4double valMax,G4String& unit,
247  G4String& schem, G4bool logVal)
248 {
249  OpenFile();
250 
251  auto hIt = IDMap.find(id0);
252  if(hIt==IDMap.end()) return false;
253 
254  auto analysisManager = G4AnalysisManager::Instance();
255  auto dup = (hIt->second)->histDup;
256  for(G4int ii=0;ii<dup;ii++)
257  {
258  G4int id = id0 + ii;
259  auto hItr = IDMap.find(id);
260  auto ht = hItr->second;
261  G4String mNam = ht->primName;
262  G4String nam = ht->meshName + "_" + ht->primName;
263  if(ht->idx>-1)
264  {
265  mNam += "_";
266  mNam += G4UIcommand::ConvertToString(ht->idx);
267  nam += "_";
268  nam += G4UIcommand::ConvertToString(ht->idx);
269  }
270  G4int hid = -1;
271  if(schem=="linear")
272  { hid = analysisManager->CreateH1(mNam,nam,nBin,valMin,valMax,unit,"none","linear"); }
273  else
274  {
275  if(logVal)
276  { hid = analysisManager->CreateH1(mNam,nam,nBin,valMin,valMax,unit,"log10","linear"); }
277  else
278  {
279  hid = analysisManager->CreateH1(mNam,nam,nBin,valMin,valMax,unit,"none","log");
280  analysisManager->SetH1XAxisIsLog(hid,true);
281  }
282  }
283 
284  if(verbose) G4cout << "eASTRunAction::Set1D for " << mNam << " / " << nam
285  << " has the histogram ID " << hid << G4endl;
286  ht->histID = hid;
287  auto pp = ht->pplotter;
288  if(pp!=nullptr) pp->Plot(ht->idx,hid);
289  }
290  return true;
291 }
292 
293 G4bool eASTRunAction::Set1DTitle(G4int id,G4String& title,G4String& x_axis,G4String&y_axis)
294 {
295  auto hItr = IDMap.find(id);
296  if(hItr==IDMap.end()) return false;
297 
298  auto analysisManager = G4AnalysisManager::Instance();
299  auto hid = hItr->second->histID;
300  analysisManager->SetH1Title(hid,title);
301  analysisManager->SetH1XAxisTitle(hid,x_axis);
302  analysisManager->SetH1YAxisTitle(hid,y_axis);
303  return true;
304 }
305 
306 G4bool eASTRunAction::Set1DYAxisLog(G4int id0,G4bool val)
307 {
308  auto hIt = IDMap.find(id0);
309  if(hIt==IDMap.end()) return false;
310  auto analysisManager = G4AnalysisManager::Instance();
311  auto dup = (hIt->second)->histDup;
312  for(G4int ii=0;ii<dup;ii++)
313  {
314  G4int id = id0 + ii;
315  auto hItr = IDMap.find(id);
316  analysisManager->SetH1YAxisIsLog(hItr->second->histID,val);
317  }
318  return true;
319 }
320 
321 // ------------- 1D profile
322 
323 G4int eASTRunAction::Create1P(G4String& mName,G4String& pName,G4int cn)
324 {
325  G4String collName = mName;
326  collName += "/";
327  collName += pName;
328  auto cID = G4SDManager::GetSDMpointer()->GetCollectionID(collName);
329  if(cID<0) return cID;
330 
331  G4int id = (cID+2*id_offset)*id_factor;
332  auto histoTypeItr = IDMap.find(id);
333  if(histoTypeItr!=IDMap.end()) return false;
334  if(verbose) G4cout << "eASTRunAction::Create1P for <" << collName
335  << "> is registered for hitCollectionID "
336  << cID << G4endl;
337 
338  auto histTyp = new eASTHistoType;
339  histTyp->collID = cID;
340  histTyp->histType = 2; // 1D profile
341  histTyp->meshName = mName;
342  histTyp->primName = pName;
343  histTyp->idx = cn;
344  IDMap[id] = histTyp;
345  return id;
346 }
347 
348 G4bool eASTRunAction::Set1P(G4int id,G4double valYMin,G4double valYMax,G4String& unit,
349  G4String& funcX,G4String& funcY,G4String& schem)
350 {
351  OpenFile();
352 
353  if(verbose) G4cout << "eASTRunAction::Set1P for id = " << id << G4endl;
354  auto hItr = IDMap.find(id);
355  if(hItr==IDMap.end()) return false;
356 
357  auto ht = hItr->second;
358  if(verbose) G4cout << "eASTRunAction::Set1P for " << ht->meshName << " / " << ht->primName << G4endl;
359  auto analysisManager = G4AnalysisManager::Instance();
360  auto nBin = ht->idx;
361  G4double valMin = -0.5;
362  G4double valMax = G4double(nBin) - 0.5;
363  G4String nam = ht->meshName + "_" + ht->primName;
364  auto hid = analysisManager->CreateP1(nam,ht->primName,nBin,
365  valMin,valMax,valYMin,valYMax,"none",unit,funcX,funcY,schem);
366 
367  if(verbose) G4cout << "eASTRunAction::Set1P for " << ht->meshName << " / " << ht->primName
368  << " has the histogram ID " << hid << G4endl;
369  ht->histID = hid;
370  return true;
371 }
372 
373 G4bool eASTRunAction::Set1PTitle(G4int id,G4String& title,G4String& x_axis,G4String&y_axis)
374 {
375  auto hItr = IDMap.find(id);
376  if(hItr==IDMap.end()) return false;
377 
378  auto analysisManager = G4AnalysisManager::Instance();
379  auto hid = hItr->second->histID;
380  analysisManager->SetP1Title(hid,title);
381  analysisManager->SetP1XAxisTitle(hid,x_axis);
382  analysisManager->SetP1YAxisTitle(hid,y_axis);
383  return true;
384 }
385 
386 // ------------- Ntuple
387 
388 G4int eASTRunAction::NtupleColumn(G4String& mName,G4String& pName,G4String& unit,G4int cn)
389 {
390  G4String collName = mName;
391  collName += "/";
392  collName += pName;
393  auto cID = G4SDManager::GetSDMpointer()->GetCollectionID(collName);
394  if(cID<0) return cID;
395 
396  G4int id = NTMap.size();
397  if(verbose) G4cout << "eASTRunAction::NtupleColumn : <" << collName
398  << ", copyNo=" << cn << "> is registered for nTuple column "
399  << id << G4endl;
400 
401  auto histTyp = new eASTHistoType;
402  histTyp->collID = cID;
403  histTyp->meshName = mName;
404  histTyp->primName = pName;
405  histTyp->meshName2 = unit;
406  if(unit!="none")
407  { histTyp->fuct = 1./(G4UnitDefinition::GetValueOf(unit)); }
408  histTyp->idx = cn;
409  NTMap[id] = histTyp;
410  return id;
411 }
412 
413 #include "G4UIcommand.hh"
414 
415 void eASTRunAction::DefineNTColumn()
416 {
417  if(NTMap.size()==0) return;
418 
419  auto analysisManager = G4AnalysisManager::Instance();
420  analysisManager->CreateNtuple("eASTNtuple","Scores for each event");
421 
422  for(auto itr : NTMap)
423  {
424  G4String colNam = itr.second->meshName;
425  colNam += "_";
426  colNam += itr.second->primName;
427  if(itr.second->idx != -1)
428  { colNam += "_"; colNam += G4UIcommand::ConvertToString(itr.second->idx); }
429  if(itr.second->meshName2 != "none")
430  { colNam += "["; colNam += itr.second->meshName2; colNam += "]"; }
431  analysisManager->CreateNtupleDColumn(colNam);
432  }
433 
434  analysisManager->FinishNtuple();
435 }
436 
437 #include <fstream>
438 #include "G4Threading.hh"
439 #include "G4UImanager.hh"
440 
441 void eASTRunAction::MergeNtuple()
442 {
443  if(NTMap.size()==0) return;
444  if(!(G4Threading::IsMultithreadedApplication())) return;
445 
446  auto analysisManager = G4AnalysisManager::Instance();
447 
448  // This MergeNtuple() method is valid only for CSV file format
449  if(analysisManager->GetType()!="Csv") return;
450 
451  std::fstream target;
452  G4String targetFN = fileName + "_nt_eASTNtuple_total.csv";
453  target.open(targetFN,std::ios::out);
454 
455  enum { BUFSIZE = 4096 };
456  char* line = new char[BUFSIZE];
457 
458  G4String titleFN = fileName + "_nt_eASTNtuple.csv";
459  std::ifstream title;
460  title.open(titleFN,std::ios::in);
461  while(1)
462  {
463  title.getline(line,BUFSIZE);
464  if(title.eof()) break;
465  G4cout << line << G4endl;
466  target << line << G4endl;
467  }
468  title.close();
469 
470  auto nWorker = G4Threading::GetNumberOfRunningWorkerThreads();
471  G4String sourceFNBase = fileName + "_nt_eASTNtuple_t";
472  for(G4int i = 0; i < nWorker; i++)
473  {
474  G4String sourceFN = sourceFNBase;
475  sourceFN += G4UIcommand::ConvertToString(i);
476  sourceFN += ".csv";
477  std::ifstream source;
478  source.open(sourceFN,std::ios::in);
479  if(!source)
480  {
481  G4ExceptionDescription ed; ed << "Source file <" << sourceFN << "> is not found.";
482  G4Exception("eASTRunAction::MergeNtuple()","eASTim12345",FatalException,ed);
483  }
484  while(1)
485  {
486  source.getline(line,BUFSIZE);
487  if(line[0]=='#') continue;
488  if(source.eof()) break;
489  target << line << G4endl;
490  }
491  source.close();
492  G4String scmd = "rm -f ";
493  scmd += sourceFN;
494  auto rc = system(scmd);
495  if(rc<0)
496  {
497  G4ExceptionDescription ed;
498  ed << "File <" << sourceFN << "> could not be deleted, thought it is merged.";
499  G4Exception("eASTRunAction::MergeNtuple()","eASTim12345",JustWarning,ed);
500  }
501  }
502 
503  target.close();
504 
505  G4String cmd = "mv ";
506  cmd += targetFN;
507  cmd += " ";
508  cmd += titleFN;
509  auto rcd = system(cmd);
510  if(rcd<0)
511  {
512  G4ExceptionDescription ed;
513  ed << "File <" << targetFN << "> could not be renamed.";
514  G4Exception("eASTRunAction::MergeNtuple()","eASTim12345",JustWarning,ed);
515  }
516 }
517 
518