gorad/src/GRRunAction.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //  Gorad (Geant4 Open-source Radiation Analysis and Design)
0027 //
0028 //  Author : Makoto Asai (SLAC National Accelerator Laboratory)
0029 //
0030 //  Development of Gorad is funded by NASA Johnson Space Center (JSC)
0031 //  under the contract NNJ15HK11B.
0032 //
0033 // ********************************************************************
0034 //
0035 // GRRunAction.cc
0036 //   Gorad Run Action class that takes care of defining and handling
0037 //   histograms and n-tuple.
0038 //   Filling histograms is taken care by GRRun class.
0039 //
0040 // History
0041 //   September 8th, 2020 : first implementation
0042 //
0043 // ********************************************************************
0044
0045 #include "GRRunAction.hh"
0046 #include "GRRunActionMessenger.hh"
0047
0048 #include "G4AnalysisManager.hh"
0049 #include "G4Run.hh"
0050 #include "G4RunManager.hh"
0051 #include "G4UnitsTable.hh"
0052 #include "G4SystemOfUnits.hh"
0053
0054 GRRunAction::GRRunAction()
0055 :fileName("GoradOut"), fileOpen(false), verbose(0), ifCarry(false),
0056  id_offset(100), id_factor(100)
0057 {
0058   messenger = new GRRunActionMessenger(this);
0059   auto analysisManager = G4AnalysisManager::Instance();
0060   analysisManager->SetDefaultFileType("csv");
0061   G4cout << "Using " << analysisManager->GetType() << G4endl;
0062
0063   analysisManager->SetVerboseLevel(verbose);
0064   //analysisManager->SetNtupleMerging(true);
0065 }
0066
0067 GRRunAction::~GRRunAction()
0068 {
0069   delete messenger;
0070 }
0071
0072 void GRRunAction::BeginOfRunAction(const G4Run* /*run*/)
0073 {
0074   // Open an output file
0075   //
0076   OpenFile();
0077
0078   // Define nTuple column if needed
0079   //
0080   DefineNTColumn();
0081 }
0082
0083 void GRRunAction::OpenFile()
0084 {
0085   if(!fileOpen)
0086   {
0087     auto analysisManager = G4AnalysisManager::Instance();
0088     analysisManager->OpenFile(fileName);
0089     if(verbose>0) G4cout << "GRRunAction::BeginOfRunAction ### <" << fileName << "> is opened." << G4endl;
0090     fileOpen = true;
0091   }
0092 }
0093
0094 void GRRunAction::EndOfRunAction(const G4Run* /*run*/)
0095 {
0096   // print histogram statistics
0097   //
0098   if(!ifCarry) Flush();
0099 }
0100
0101 void GRRunAction::Flush()
0102 {
0103   auto analysisManager = G4AnalysisManager::Instance();
0104
0105   analysisManager->Write();
0106   analysisManager->CloseFile();
0107   if(verbose>0) G4cout << "GRRunAction::Flush ### <" << fileName << "> is closed." << G4endl;
0108
0109   fileOpen = false;
0110
0111   if(IsMaster()) MergeNtuple();
0112 }
0113
0114 void GRRunAction::SetVerbose(G4int val)
0115 {
0116   verbose = val;
0117   auto analysisManager = G4AnalysisManager::Instance();
0118   analysisManager->SetVerboseLevel(verbose);
0119 }
0120
0121 void GRRunAction::ListHistograms()
0122 {
0123   G4cout << "################## registered histograms/plots" << G4endl;
0124   G4cout << "id\thistID\thistType\tdetName-X\tpsName-X\tcollID-X\tcopyNo-X\tdetName-Y\tpsName-Y\tcollID-Y\tcopyNo-Y" << G4endl;
0125   for(auto itr : IDMap)
0126   {
0127     G4cout << itr.first << "\t" << itr.second->histID << "\t";
0128     if(itr.second->histType==1) // 1D histogram
0129     { G4cout << "1-D hist\t" << itr.second->meshName << "\t" << itr.second->primName << "\t" << itr.second->collID << "\t" << itr.second->idx; }
0130     else if(itr.second->histType==2) // 1D profile
0131     { G4cout << "1-D prof\t" << itr.second->meshName << "\t" << itr.second->primName << "\t" << itr.second->collID; }
0132     G4cout << G4endl;
0133   }
0134 }
0135
0136 G4bool GRRunAction::Open(G4int id)
0137 {
0138   auto hItr = IDMap.find(id);
0139   return (hItr!=IDMap.end());
0140 }
0141
0142 #include "G4SDManager.hh"
0143 using namespace G4Analysis;
0144
0145 G4bool GRRunAction::SetAllPlotting(G4bool val)
0146 {
0147   G4bool valid = true;
0148   for(auto hItr : IDMap)
0149   {
0150     valid = SetPlotting(hItr.first,val);
0151     if(!valid) break;
0152   }
0153   return valid;
0154 }
0155
0156 G4bool GRRunAction::SetPlotting(G4int id,G4bool val)
0157 {
0158   auto hItr = IDMap.find(id);
0159   if(hItr==IDMap.end()) return false;
0160   auto ht = hItr->second;
0161   auto hTyp = ht->histType;
0162   auto analysisManager = G4AnalysisManager::Instance();
0163   if(hTyp==1) // 1D histogram
0164   { analysisManager->SetH1Plotting(ht->histID,val); }
0165   else if(hTyp==2) // 1D profile
0166   { analysisManager->SetP1Plotting(ht->histID,val); }
0167   else
0168   { return false; }
0169   return true;
0170 }
0171
0172 // ------------- 1D histogram
0173
0174 G4int GRRunAction::Create1D(G4String& mName,G4String& pName,G4int cn)
0175 {
0176   G4String collName = mName;
0177   collName += "/";
0178   collName += pName;
0179   auto cID = G4SDManager::GetSDMpointer()->GetCollectionID(collName);
0180   if(cID<0) return cID;
0181
0182   G4int id = (cID+id_offset)*id_factor+cn+1;
0183   auto histoTypeItr = IDMap.find(id);
0184   if(histoTypeItr!=IDMap.end()) return false;
0185   if(verbose) G4cout << "GRRunAction::Create1D for <" << collName
0186                      << ", copyNo=" << cn << "> is registered for hitCollectionID "
0187                      << cID << G4endl;
0188
0189   auto histTyp = new GRHistoType;
0190   histTyp->collID = cID;
0191   histTyp->histType = 1; // 1D histogram
0192   histTyp->meshName = mName;
0193   histTyp->primName = pName;
0194   histTyp->idx = cn;
0195   IDMap[id] = histTyp;
0196   return id;
0197 }
0198
0199 G4int GRRunAction::Create1DForPrimary(G4String& mName,G4bool wgt)
0200 {
0201   G4int cn = wgt ? 1 : 0;
0202
0203   G4int id = 99999 - cn;
0204   auto histoTypeItr = IDMap.find(id);
0205   if(histoTypeItr!=IDMap.end()) return false;
0206   if(verbose) G4cout << "GRRunAction::Create1D for <" << mName
0207                      << "(weighted : " << cn << ")> is registered " << G4endl;
0208
0209   auto histTyp = new GRHistoType;
0210   histTyp->collID = -999;
0211   histTyp->histType = 1; // 1D histogram
0212   histTyp->meshName = "PrimPEnergy";
0213   histTyp->primName = mName;
0214   histTyp->biasf = cn;
0215   IDMap[id] = histTyp;
0216   return id;
0217 }
0218
0219 #include "G4SDManager.hh"
0220 #include "G4ScoringManager.hh"
0221 #include "G4VScoringMesh.hh"
0222 #include "G4VPrimitivePlotter.hh"
0223 G4int GRRunAction::Create1DForPlotter(G4String& mName,G4String& pName,G4bool /*wgt*/)
0224 {
0225   using MeshShape = G4VScoringMesh::MeshShape;
0226
0227   G4String collName = mName;
0228   collName += "/";
0229   collName += pName;
0230   auto cID = G4SDManager::GetSDMpointer()->GetCollectionID(collName);
0231   if(cID<0) return cID;
0232
0233   auto sm = G4ScoringManager::GetScoringManagerIfExist();
0234   assert(sm!=nullptr);
0235   auto mesh = sm->FindMesh(mName);
0236   if(mesh==nullptr)
0237   { return -2; }
0238   auto shape = mesh->GetShape();
0239   if(shape!=MeshShape::realWorldLogVol && shape!=MeshShape::probe)
0240   { return -3; }
0241   G4int nBin[3];
0242   mesh->GetNumberOfSegments(nBin);
0243
0244   auto prim = mesh->GetPrimitiveScorer(pName);
0245   if(prim==nullptr)
0246   { return -3; }
0247   auto pp = dynamic_cast<G4VPrimitivePlotter*>(prim);
0248   if(pp==nullptr)
0249   { return -4; }
0250
0251   G4int id0 = (cID+id_offset)*id_factor+1;
0252   for(G4int cn=0; cn<nBin[0]; cn++)
0253   {
0254     G4int id = id0+cn;
0255     auto histoTypeItr = IDMap.find(id);
0256     if(histoTypeItr!=IDMap.end())
0257     { return -5; }
0258     if(verbose) G4cout << "GRRunAction::Create1D for <" << collName
0259                      << ", copyNo=" << cn << "> is registered for hitCollectionID "
0260                      << cID << G4endl;
0261
0262     auto histTyp = new GRHistoType;
0263     histTyp->collID = cID;
0264     histTyp->histType = 1; // 1D histogram
0265     histTyp->histDup = nBin[0];
0266     histTyp->meshName = mName;
0267     histTyp->primName = pName;
0268     histTyp->idx = cn;
0269     histTyp->pplotter = pp;
0270     IDMap[id] = histTyp;
0271   }
0272   return id0;
0273 }
0274
0275 #include "G4UIcommand.hh"
0276 G4bool GRRunAction::Set1D(G4int id0,G4int nBin,G4double valMin,G4double valMax,G4String& unit,
0277                           G4String& schem, G4bool logVal)
0278 {
0279   OpenFile();
0280
0281   auto hIt = IDMap.find(id0);
0282   if(hIt==IDMap.end()) return false;
0283
0284   auto analysisManager = G4AnalysisManager::Instance();
0285   auto dup = (hIt->second)->histDup;
0286   for(G4int ii=0;ii<dup;ii++)
0287   {
0288     G4int id = id0 + ii;
0289     auto hItr = IDMap.find(id);
0290     auto ht = hItr->second;
0291     G4String mNam = ht->primName;
0292     G4String nam = ht->meshName + "_" + ht->primName;
0293     if(ht->idx>-1)
0294     {
0295       mNam += "_";
0296       mNam += G4UIcommand::ConvertToString(ht->idx);
0297       nam += "_";
0298       nam += G4UIcommand::ConvertToString(ht->idx);
0299     }
0300     G4int hid = -1;
0301     if(schem=="linear")
0302     { hid = analysisManager->CreateH1(mNam,nam,nBin,valMin,valMax,unit,"none","linear"); }
0303     else
0304     {
0305       if(logVal)
0306       { hid = analysisManager->CreateH1(mNam,nam,nBin,valMin,valMax,unit,"log10","linear"); }
0307       else
0308       {
0309         hid = analysisManager->CreateH1(mNam,nam,nBin,valMin,valMax,unit,"none","log");
0310         analysisManager->SetH1XAxisIsLog(hid,true);
0311       }
0312     }
0313
0314     if(verbose) G4cout << "GRRunAction::Set1D for " << mNam << " / " << nam
0315                        << " has the histogram ID " << hid << G4endl;
0316     ht->histID = hid;
0317     auto pp = ht->pplotter;
0318     if(pp!=nullptr) pp->Plot(ht->idx,hid);
0319   }
0320   return true;
0321 }
0322
0323 G4bool GRRunAction::Set1DTitle(G4int id,G4String& title,G4String& x_axis,G4String&y_axis)
0324 {
0325   auto hItr = IDMap.find(id);
0326   if(hItr==IDMap.end()) return false;
0327
0328   auto analysisManager = G4AnalysisManager::Instance();
0329   auto hid = hItr->second->histID;
0330   analysisManager->SetH1Title(hid,title);
0331   analysisManager->SetH1XAxisTitle(hid,x_axis);
0332   analysisManager->SetH1YAxisTitle(hid,y_axis);
0333   return true;
0334 }
0335
0336 G4bool GRRunAction::Set1DYAxisLog(G4int id0,G4bool val)
0337 {
0338   auto hIt = IDMap.find(id0);
0339   if(hIt==IDMap.end()) return false;
0340   auto analysisManager = G4AnalysisManager::Instance();
0341   auto dup = (hIt->second)->histDup;
0342   for(G4int ii=0;ii<dup;ii++)
0343   {
0344     G4int id = id0 + ii;
0345     auto hItr = IDMap.find(id);
0346     analysisManager->SetH1YAxisIsLog(hItr->second->histID,val);
0347   }
0348   return true;
0349 }
0350
0351 // ------------- 1D profile
0352
0353 G4int GRRunAction::Create1P(G4String& mName,G4String& pName,G4int cn)
0354 {
0355   G4String collName = mName;
0356   collName += "/";
0357   collName += pName;
0358   auto cID = G4SDManager::GetSDMpointer()->GetCollectionID(collName);
0359   if(cID<0) return cID;
0360
0361   G4int id = (cID+2*id_offset)*id_factor;
0362   auto histoTypeItr = IDMap.find(id);
0363   if(histoTypeItr!=IDMap.end()) return false;
0364   if(verbose) G4cout << "GRRunAction::Create1P for <" << collName
0365                      << "> is registered for hitCollectionID "
0366                      << cID << G4endl;
0367
0368   auto histTyp = new GRHistoType;
0369   histTyp->collID = cID;
0370   histTyp->histType = 2; // 1D profile
0371   histTyp->meshName = mName;
0372   histTyp->primName = pName;
0373   histTyp->idx = cn;
0374   IDMap[id] = histTyp;
0375   return id;
0376 }
0377
0378 G4bool GRRunAction::Set1P(G4int id,G4double valYMin,G4double valYMax,G4String& unit,
0379         G4String& funcX,G4String& funcY,G4String& schem)
0380 {
0381   OpenFile();
0382
0383   if(verbose) G4cout << "GRRunAction::Set1P for id = " << id << G4endl;
0384   auto hItr = IDMap.find(id);
0385   if(hItr==IDMap.end()) return false;
0386
0387   auto ht = hItr->second;
0388   if(verbose) G4cout << "GRRunAction::Set1P for " << ht->meshName << " / " << ht->primName << G4endl;
0389   auto analysisManager = G4AnalysisManager::Instance();
0390   auto nBin = ht->idx;
0391   G4double valMin = -0.5;
0392   G4double valMax = G4double(nBin) - 0.5;
0393   G4String nam = ht->meshName + "_" + ht->primName;
0394   auto hid = analysisManager->CreateP1(nam,ht->primName,nBin,
0395               valMin,valMax,valYMin,valYMax,"none",unit,funcX,funcY,schem);
0396
0397   if(verbose) G4cout << "GRRunAction::Set1P for " << ht->meshName << " / " << ht->primName
0398                      << " has the histogram ID " << hid << G4endl;
0399   ht->histID = hid;
0400   return true;
0401 }
0402
0403 G4bool GRRunAction::Set1PTitle(G4int id,G4String& title,G4String& x_axis,G4String&y_axis)
0404 {
0405   auto hItr = IDMap.find(id);
0406   if(hItr==IDMap.end()) return false;
0407
0408   auto analysisManager = G4AnalysisManager::Instance();
0409   auto hid = hItr->second->histID;
0410   analysisManager->SetP1Title(hid,title);
0411   analysisManager->SetP1XAxisTitle(hid,x_axis);
0412   analysisManager->SetP1YAxisTitle(hid,y_axis);
0413   return true;
0414 }
0415
0416 // ------------- Ntuple
0417
0418 G4int GRRunAction::NtupleColumn(G4String& mName,G4String& pName,G4String& unit,G4int cn)
0419 {
0420   G4String collName = mName;
0421   collName += "/";
0422   collName += pName;
0423   auto cID = G4SDManager::GetSDMpointer()->GetCollectionID(collName);
0424   if(cID<0) return cID;
0425
0426   G4int id = NTMap.size();
0427   if(verbose) G4cout << "GRRunAction::NtupleColumn : <" << collName
0428                      << ", copyNo=" << cn << "> is registered for nTuple column "
0429                      << id << G4endl;
0430
0431   auto histTyp = new GRHistoType;
0432   histTyp->collID = cID;
0433   histTyp->meshName = mName;
0434   histTyp->primName = pName;
0435   histTyp->meshName2 = unit;
0436   if(unit!="none")
0437   { histTyp->fuct = 1./(G4UnitDefinition::GetValueOf(unit)); }
0438   histTyp->idx = cn;
0439   NTMap[id] = histTyp;
0440   return id;
0441 }
0442
0443 #include "G4UIcommand.hh"
0444
0445 void GRRunAction::DefineNTColumn()
0446 {
0447   if(NTMap.size()==0) return;
0448
0449   auto analysisManager = G4AnalysisManager::Instance();
0450   analysisManager->CreateNtuple("GRimNtuple","Scores for each event");
0451
0452   for(auto itr : NTMap)
0453   {
0454     G4String colNam = itr.second->meshName;
0455     colNam += "_";
0456     colNam += itr.second->primName;
0457     if(itr.second->idx != -1)
0458     { colNam += "_"; colNam += G4UIcommand::ConvertToString(itr.second->idx); }
0459     if(itr.second->meshName2 != "none")
0460     { colNam += "["; colNam += itr.second->meshName2; colNam += "]"; }
0461     analysisManager->CreateNtupleDColumn(colNam);
0462   }
0463
0464   analysisManager->FinishNtuple();
0465 }
0466
0467 #include <fstream>
0468 #include "G4Threading.hh"
0469 #include "G4UImanager.hh"
0470
0471 void GRRunAction::MergeNtuple()
0472 {
0473   if(NTMap.size()==0) return;
0474   if(!(G4Threading::IsMultithreadedApplication())) return;
0475
0476   auto analysisManager = G4AnalysisManager::Instance();
0477
0478   // This MergeNtuple() method is valid only for CSV file format
0479   if(analysisManager->GetType()!="Csv") return;
0480
0481   std::fstream target;
0482   G4String targetFN = "GRimOut_nt_GRimNtuple_total.csv";
0483   target.open(targetFN,std::ios::out);
0484
0485   enum { BUFSIZE = 4096 };
0486   char* line = new char[BUFSIZE];
0487
0488   G4String titleFN = "GRimOut_nt_GRimNtuple.csv";
0489   std::ifstream title;
0490   title.open(titleFN,std::ios::in);
0491   while(1)
0492   {
0493     title.getline(line,BUFSIZE);
0494     if(title.eof()) break;
0495     G4cout << line << G4endl;
0496     target << line << G4endl;
0497   }
0498   title.close();
0499
0500   auto nWorker = G4Threading::GetNumberOfRunningWorkerThreads();
0501   G4String sourceFNBase = "GRimOut_nt_GRimNtuple_t";
0502   for(G4int i = 0; i < nWorker; i++)
0503   {
0504     G4String sourceFN = sourceFNBase;
0505     sourceFN += G4UIcommand::ConvertToString(i);
0506     sourceFN += ".csv";
0507     std::ifstream source;
0508     source.open(sourceFN,std::ios::in);
0509     if(!source)
0510     {
0511       G4ExceptionDescription ed; ed << "Source file <" << sourceFN << "> is not found.";
0512       G4Exception("GRRunAction::MergeNtuple()","GRim12345",FatalException,ed);
0513     }
0514     while(1)
0515     {
0516       source.getline(line,BUFSIZE);
0517       if(line[0]=='#') continue;
0518       if(source.eof()) break;
0519       target << line << G4endl;
0520     }
0521     source.close();
0522     G4String scmd = "rm -f ";
0523     scmd += sourceFN;
0524     auto rc = system(scmd);
0525     if(rc<0)
0526     {
0527       G4ExceptionDescription ed;
0528       ed << "File <" << sourceFN << "> could not be deleted, thought it is merged.";
0529       G4Exception("GRRunAction::MergeNtuple()","GRim12345",JustWarning,ed);
0530     }
0531   }
0532
0533   target.close();
0534
0535   G4String cmd = "mv ";
0536   cmd += targetFN;
0537   cmd += " ";
0538   cmd += titleFN;
0539   auto rcd = system(cmd);
0540   if(rcd<0)
0541   {
0542     G4ExceptionDescription ed;
0543     ed << "File <" << targetFN << "> could not be renamed.";
0544     G4Exception("GRRunAction::MergeNtuple()","GRim12345",JustWarning,ed);
0545   }
0546 }
0547
0548