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