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