| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Gorad (Geant4 Open-source Radiation Analys
27 //
28 // Author : Makoto Asai (SLAC National Accele
29 //
30 // Development of Gorad is funded by NASA Joh
31 // under the contract NNJ15HK11B.
32 //
33 // *******************************************
34 //
35 // GRRunAction.cc
36 // Gorad Run Action class that takes care of
37 // histograms and n-tuple.
38 // Filling histograms is taken care by GRRun
39 //
40 // History
41 // September 8th, 2020 : first implementatio
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), verbos
56 id_offset(100), id_factor(100)
57 {
58 messenger = new GRRunActionMessenger(this);
59 auto analysisManager = G4AnalysisManager::In
60 analysisManager->SetDefaultFileType("csv");
61 G4cout << "Using " << analysisManager->GetTy
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
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::
88 analysisManager->OpenFile(fileName);
89 if(verbose>0) G4cout << "GRRunAction::Begi
90 fileOpen = true;
91 }
92 }
93
94 void GRRunAction::EndOfRunAction(const G4Run*
95 {
96 // print histogram statistics
97 //
98 if(!ifCarry) Flush();
99 }
100
101 void GRRunAction::Flush()
102 {
103 auto analysisManager = G4AnalysisManager::In
104
105 analysisManager->Write();
106 analysisManager->CloseFile();
107 if(verbose>0) G4cout << "GRRunAction::Flush
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::In
118 analysisManager->SetVerboseLevel(verbose);
119 }
120
121 void GRRunAction::ListHistograms()
122 {
123 G4cout << "################## registered his
124 G4cout << "id\thistID\thistType\tdetName-X\t
125 for(auto itr : IDMap)
126 {
127 G4cout << itr.first << "\t" << itr.second-
128 if(itr.second->histType==1) // 1D histogra
129 { G4cout << "1-D hist\t" << itr.second->me
130 else if(itr.second->histType==2) // 1D pro
131 { G4cout << "1-D prof\t" << itr.second->me
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,G4boo
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::In
163 if(hTyp==1) // 1D histogram
164 { analysisManager->SetH1Plotting(ht->histID,
165 else if(hTyp==2) // 1D profile
166 { analysisManager->SetP1Plotting(ht->histID,
167 else
168 { return false; }
169 return true;
170 }
171
172 // ------------- 1D histogram
173
174 G4int GRRunAction::Create1D(G4String& mName,G4
175 {
176 G4String collName = mName;
177 collName += "/";
178 collName += pName;
179 auto cID = G4SDManager::GetSDMpointer()->Get
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
186 << ", copyNo=" << cn << "
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
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
207 << "(weighted : " << cn <
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
224 {
225 using MeshShape = G4VScoringMesh::MeshShape;
226
227 G4String collName = mName;
228 collName += "/";
229 collName += pName;
230 auto cID = G4SDManager::GetSDMpointer()->Get
231 if(cID<0) return cID;
232
233 auto sm = G4ScoringManager::GetScoringManage
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 && shap
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*>
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::Create
259 << ", copyNo=" << cn << "
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
277 G4String& schem, G4b
278 {
279 OpenFile();
280
281 auto hIt = IDMap.find(id0);
282 if(hIt==IDMap.end()) return false;
283
284 auto analysisManager = G4AnalysisManager::In
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->pr
293 if(ht->idx>-1)
294 {
295 mNam += "_";
296 mNam += G4UIcommand::ConvertToString(ht-
297 nam += "_";
298 nam += G4UIcommand::ConvertToString(ht->
299 }
300 G4int hid = -1;
301 if(schem=="linear")
302 { hid = analysisManager->CreateH1(mNam,nam
303 else
304 {
305 if(logVal)
306 { hid = analysisManager->CreateH1(mNam,n
307 else
308 {
309 hid = analysisManager->CreateH1(mNam,n
310 analysisManager->SetH1XAxisIsLog(hid,t
311 }
312 }
313
314 if(verbose) G4cout << "GRRunAction::Set1D
315 << " has the histogram
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,G4Stri
324 {
325 auto hItr = IDMap.find(id);
326 if(hItr==IDMap.end()) return false;
327
328 auto analysisManager = G4AnalysisManager::In
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,G4
337 {
338 auto hIt = IDMap.find(id0);
339 if(hIt==IDMap.end()) return false;
340 auto analysisManager = G4AnalysisManager::In
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->sec
347 }
348 return true;
349 }
350
351 // ------------- 1D profile
352
353 G4int GRRunAction::Create1P(G4String& mName,G4
354 {
355 G4String collName = mName;
356 collName += "/";
357 collName += pName;
358 auto cID = G4SDManager::GetSDMpointer()->Get
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
365 << "> is registered for h
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 va
379 G4String& funcX,G4String& funcY,G4Stri
380 {
381 OpenFile();
382
383 if(verbose) G4cout << "GRRunAction::Set1P fo
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 fo
389 auto analysisManager = G4AnalysisManager::In
390 auto nBin = ht->idx;
391 G4double valMin = -0.5;
392 G4double valMax = G4double(nBin) - 0.5;
393 G4String nam = ht->meshName + "_" + ht->prim
394 auto hid = analysisManager->CreateP1(nam,ht-
395 valMin,valMax,valYMin,valYMax,"n
396
397 if(verbose) G4cout << "GRRunAction::Set1P fo
398 << " has the histogram ID
399 ht->histID = hid;
400 return true;
401 }
402
403 G4bool GRRunAction::Set1PTitle(G4int id,G4Stri
404 {
405 auto hItr = IDMap.find(id);
406 if(hItr==IDMap.end()) return false;
407
408 auto analysisManager = G4AnalysisManager::In
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& mNam
419 {
420 G4String collName = mName;
421 collName += "/";
422 collName += pName;
423 auto cID = G4SDManager::GetSDMpointer()->Get
424 if(cID<0) return cID;
425
426 G4int id = NTMap.size();
427 if(verbose) G4cout << "GRRunAction::NtupleCo
428 << ", copyNo=" << cn << "
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::GetV
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::In
450 analysisManager->CreateNtuple("GRimNtuple","
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::Co
459 if(itr.second->meshName2 != "none")
460 { colNam += "["; colNam += itr.second->mes
461 analysisManager->CreateNtupleDColumn(colNa
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
475
476 auto analysisManager = G4AnalysisManager::In
477
478 // This MergeNtuple() method is valid only f
479 if(analysisManager->GetType()!="Csv") return
480
481 std::fstream target;
482 G4String targetFN = "GRimOut_nt_GRimNtuple_t
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.cs
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::GetNumberOfRunni
501 G4String sourceFNBase = "GRimOut_nt_GRimNtup
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
512 G4Exception("GRRunAction::MergeNtuple()"
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 n
529 G4Exception("GRRunAction::MergeNtuple()"
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
544 G4Exception("GRRunAction::MergeNtuple()","
545 }
546 }
547
548
549