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 // 27 /// \file Analysis.cc 28 /// \brief Implementation of the Analysis clas 29 /// \file Analysis.cc 30 /// \brief Implementation of the Analysis clas 31 32 #include "Analysis.hh" 33 #include "DetectorConstruction.hh" 34 #include "G4AutoDelete.hh" 35 #include "G4Filesystem.hh" 36 #include "G4DNAMolecule.hh" 37 38 39 #include <sstream> 40 #include <fstream> 41 #ifdef USE_MPI 42 #include "G4MPImanager.hh" 43 #endif 44 45 G4ThreadLocal Analysis* the_analysis = nullptr 46 47 //....oooOO0OOooo........oooOO0OOooo........oo 48 49 Analysis* Analysis::GetAnalysis() 50 { 51 if (!the_analysis) { 52 the_analysis = new Analysis(); 53 G4AutoDelete::Register(the_analysis); 54 } 55 56 return the_analysis; 57 } 58 59 //....oooOO0OOooo........oooOO0OOooo........oo 60 61 void Analysis::OpenFile(const G4String outFold 62 { 63 G4AnalysisManager* anManager = G4AnalysisM 64 anManager->SetDefaultFileType("root"); 65 G4String fullFileName = fFileName; 66 if (gRunMode == RunningMode::Chem) { 67 G4String slash = ""; 68 #if defined(_WIN32) || defined(WIN32) 69 slash= "\\"; 70 #else 71 G4String slashu = "/"; 72 slash = slashu; 73 #endif 74 fullFileName = outFolder+slash+fFileNa 75 } 76 77 // open output file 78 G4bool fileOpen = anManager->OpenFile(full 79 if (!fileOpen) { 80 G4cout << "\n---> HistoManager::book() 81 << G4endl; 82 return; 83 } 84 } 85 86 //....oooOO0OOooo........oooOO0OOooo........oo 87 88 void Analysis::Save() 89 { 90 G4AnalysisManager* anManager = G4AnalysisM 91 anManager->Write(); 92 } 93 94 //....oooOO0OOooo........oooOO0OOooo........oo 95 96 void Analysis::Close(G4bool reset) 97 { 98 G4AnalysisManager* anManager = G4AnalysisM 99 anManager->CloseFile(reset); 100 } 101 102 //....oooOO0OOooo........oooOO0OOooo........oo 103 104 void Analysis::Book() 105 { 106 G4AnalysisManager* anManager = G4AnalysisM 107 anManager->SetVerboseLevel(1); 108 if (anManager->GetFirstNtupleId() != 1) an 109 110 if (gRunMode == RunningMode::Phys) { 111 #ifdef G4MULTITHREADED 112 // MT ntuple merging 113 anManager->SetNtupleMerging(true);//fo 114 #endif 115 anManager->SetNtupleDirectoryName("ntu 116 // create ntuple 117 anManager->CreateNtuple("ntuple_1","ph 118 anManager->CreateNtupleIColumn("flagPa 119 anManager->CreateNtupleIColumn("flagPa 120 anManager->CreateNtupleIColumn("flagPr 121 anManager->CreateNtupleDColumn("x"); 122 anManager->CreateNtupleDColumn("y"); 123 anManager->CreateNtupleDColumn("z"); 124 anManager->CreateNtupleDColumn("edep") 125 anManager->CreateNtupleIColumn("eventN 126 anManager->CreateNtupleIColumn("volume 127 anManager->CreateNtupleIColumn("copyNu 128 anManager->CreateNtupleIColumn("lastMe 129 anManager->FinishNtuple(1); 130 131 // For total edep 132 anManager->CreateNtuple("ntuple_3","to 133 anManager->CreateNtupleIColumn("eventN 134 anManager->CreateNtupleDColumn("edep") 135 anManager->FinishNtuple(2); 136 } 137 138 if (gRunMode == RunningMode::Chem) { 139 // Create directories 140 const G4String directoryName = "ntuple 141 if (anManager->GetNtupleDirectoryName( 142 anManager->SetNtupleDirectoryName( 143 } 144 145 // DBScan 146 147 anManager->CreateNtuple("ntuple_2","DB 148 anManager->CreateNtupleIColumn(1,"stra 149 anManager->CreateNtupleIColumn(1,"copy 150 anManager->CreateNtupleDColumn(1,"xp") 151 anManager->CreateNtupleDColumn(1,"yp") 152 anManager->CreateNtupleDColumn(1,"zp") 153 anManager->CreateNtupleDColumn(1,"time 154 anManager->CreateNtupleIColumn(1,"base 155 anManager->FinishNtuple(1); 156 } 157 } 158 159 //....oooOO0OOooo........oooOO0OOooo........oo 160 161 G4AnalysisManager* Analysis::GetAnalysisManage 162 { 163 return G4AnalysisManager::Instance(); 164 } 165 166 //....oooOO0OOooo........oooOO0OOooo........oo 167 168 void Analysis::AddInfoForChemGeo(InfoForChemGe 169 { 170 fInfoForChemGeoVector.push_back(b); 171 } 172 173 //....oooOO0OOooo........oooOO0OOooo........oo 174 175 void Analysis::ClearVector() 176 { 177 fInfoForChemGeoVector.clear(); 178 fInfoInPhysStageVector.clear(); 179 fOutputFiles.clear(); 180 } 181 182 //....oooOO0OOooo........oooOO0OOooo........oo 183 184 void Analysis::AddInfoInPhysStage(InfoInPhysSt 185 { 186 fInfoInPhysStageVector.push_back(b); 187 } 188 189 //....oooOO0OOooo........oooOO0OOooo........oo 190 191 void Analysis::UpdateChemInputDataAndFillNtupl 192 { 193 //if (fInfoForChemGeoVector.size() == 0) r 194 195 // We will loop on all the "element" of th 196 // - New file for each event/voxel couple 197 // Create a map to define all the output f 198 199 for (auto const& binfo : fInfoForChemGeoVe 200 if( binfo.fVolume == 161 // voxelStra 201 || binfo.fVolume == 162 // voxelRi 202 || binfo.fVolume == 163 // voxelLe 203 || binfo.fVolume == 164 // voxelUp 204 || binfo.fVolume == 165 // voxelDo 205 || binfo.fVolume == 261 // voxelSt 206 || binfo.fVolume == 262 // voxelRi 207 || binfo.fVolume == 263 // voxelLe 208 || binfo.fVolume == 264 // voxelUp 209 || binfo.fVolume == 265) // voxelD 210 { 211 // We are in a voxel 212 213 std::string voxelName(""); 214 215 if(binfo.fVolume==161) voxelName = 216 else if(binfo.fVolume==162) voxelN 217 else if(binfo.fVolume==163) voxelN 218 else if(binfo.fVolume==164) voxelN 219 else if(binfo.fVolume==165) voxelN 220 else if(binfo.fVolume==261) voxelN 221 else if(binfo.fVolume==262) voxelN 222 else if(binfo.fVolume==263) voxelN 223 else if(binfo.fVolume==264) voxelN 224 else if(binfo.fVolume==265) voxelN 225 226 // Get the event number 227 int eventNum = int(binfo.fEventNum 228 // Here we have all the informatio 229 // We want to know if we have to c 230 // or if we already have one. 231 232 // Check if the event has already 233 if(fOutputFiles.find(eventNum)==fO 234 { 235 // If not then create the even 236 fOutputFiles[eventNum][binfo.f 237 CreateChemInputFile(eventNum, 238 } 239 // If the event has been registere 240 // registered. 241 else 242 { 243 if(fOutputFiles[eventNum].find 244 { 245 // We register the event-v 246 fOutputFiles[eventNum][bin 247 CreateChemInputFile(eventN 248 } 249 } 250 251 // Write in the file the informati 252 // the input water molecule or sol 253 UpdatingChemInputFile(binfo); 254 } 255 } 256 257 if (fInfoInPhysStageVector.size() == 0) re 258 for (auto const& binfo : fInfoInPhysStageV 259 // Check if the process is an ionisati 260 // Only ionisation should trigger the 261 if(binfo.fFlagProcess == 13 262 || binfo.fFlagProcess == 113 263 || binfo.fFlagProcess == 18 264 || binfo.fFlagProcess == 21 265 || binfo.fFlagProcess == 24 266 || binfo.fFlagProcess == 27 267 || binfo.fFlagProcess == 31) { 268 // Check the interaction happe 269 if( binfo.fVolumeName == 1 // 270 || binfo.fVolumeName == 11 271 || binfo.fVolumeName == 2 272 || binfo.fVolumeName == 22 273 || binfo.fVolumeName == 3 274 || binfo.fVolumeName == 4 275 || binfo.fVolumeName == 5 276 || binfo.fVolumeName == 6 277 || binfo.fVolumeName == 7 278 || binfo.fVolumeName == 71 279 || binfo.fVolumeName == 8 280 || binfo.fVolumeName == 81 281 || binfo.fVolumeName == 9 282 || binfo.fVolumeName == 10 283 || binfo.fVolumeName == 13 284 || binfo.fVolumeName == 12 285 { 286 // Retrieve the voxel copy 287 double voxelCopyNumber = 288 // Get the event number 289 double eventNum = binfo.fE 290 291 // Check if the event has 292 if(fOutputFiles.find(event 293 { 294 // Check if the volume 295 if(fOutputFiles.at(eve 296 != fOutputFiles.at 297 { 298 // If we are here 299 //file in which we 300 UpdatingChemInputF 301 } 302 } 303 } 304 } 305 306 // fill ntuple 307 auto analysisManager =G4AnalysisManage 308 analysisManager->FillNtupleIColumn(1, 309 analysisManager->FillNtupleIColumn(1, 310 analysisManager->FillNtupleIColumn(1, 311 analysisManager->FillNtupleDColumn(1, 312 analysisManager->FillNtupleDColumn(1, 313 analysisManager->FillNtupleDColumn(1, 314 analysisManager->FillNtupleDColumn(1, 315 analysisManager->FillNtupleIColumn(1, 316 analysisManager->FillNtupleIColumn(1, 317 analysisManager->FillNtupleIColumn(1, 318 analysisManager->FillNtupleIColumn(1, 319 analysisManager->AddNtupleRow(1); 320 } 321 322 } 323 324 //....oooOO0OOooo........oooOO0OOooo........oo 325 326 G4String Analysis::CreateChemInputFile(G4int e 327 { 328 std::stringstream sstream; 329 sstream<<"./"<<fChemInputFolderName<<"/eve 330 std::ofstream oFile; 331 oFile.open(sstream.str().c_str() ); 332 333 oFile<<"_eventNum"<<"\t\t"<<eventNum<<std: 334 oFile<<"_voxelType"<<"\t\t"<<voxelName<<st 335 oFile<<"_voxelCopyNumber"<<"\t\t"<<volumeC 336 oFile<<std::endl; 337 338 oFile 339 <<"# Chemistry input informations\ 340 <<"# "<<"_input, type, state, elec 341 oFile<<"# type=1 -> water molecule && type 342 oFile<<std::endl; 343 344 oFile.close(); 345 346 return sstream.str().c_str(); 347 } 348 349 //....oooOO0OOooo........oooOO0OOooo........oo 350 351 void Analysis::UpdatingChemInputFile(InfoForCh 352 { 353 std::ofstream out; 354 out.open( fOutputFiles[b.fEventNumber][b.f 355 356 if (b.fType==1) 357 { 358 out<<"_input"<<"\t" 359 << 360 << 361 <<4-b. 362 <<b.fR 363 <<b.fR 364 <<b.fR 365 <<b.fParen 366 <<"\n"; 367 } 368 if (b.fType==2) 369 { 370 out<<"_input"<<"\t" 371 << 372 << 373 <<b.fE 374 // If we are here then the event and volume co 375 //generated file in which we should add a dna 376 <<b.fR 377 <<b.fR 378 <<b.fR 379 <<b.fParen 380 <<"\n"; 381 } 382 out.close(); 383 } 384 385 //....oooOO0OOooo........oooOO0OOooo........oo 386 387 void Analysis::UpdatingChemInputFile(InfoInPhy 388 { 389 G4String name; 390 G4double volumeName = b.fVolumeName; 391 if(volumeName == 1 || volumeName == 2 || v 392 volumeName == 8) name = G4Deoxyri 393 else if(volumeName == 11 || volumeName == 394 volumeName == 71 || volumeName == 395 else if(volumeName == 6 || volumeName == 9 396 else if(volumeName == 4 || volumeName == 1 397 else if(volumeName == 5 || volumeName == 1 398 else if(volumeName == 3 || volumeName == 1 399 else 400 { 401 G4ExceptionDescription msg; 402 msg <<"Volume number "<<volumeName<<" 403 G4Exception("Analysis::UpdatingChemInp 404 "", FatalException, msg); 405 } 406 407 G4double strand (-1); 408 409 // Determine the strand 410 if(volumeName==1 411 || volumeName==11 412 || volumeName==7 413 || volumeName==71 414 || volumeName==6 // ade 415 || volumeName==9 // ade 416 || volumeName==4 // gua 417 || volumeName==10) // gua 418 strand = 1; 419 else if(volumeName==2 420 || volumeName==22 421 || volumeName==8 422 || volumeName==81 423 || volumeName==5 // thy 424 || volumeName==12 // thy 425 || volumeName==3 // cyto 426 || volumeName==13) // cyto 427 strand = 2; 428 std::ofstream out; 429 out.open(fOutputFiles[b.fEventNumber][b.fL 430 out<<"_remove"<<"\t"<<name<<"\t\t"<<b.fCop 431 out.close(); 432 } 433 434 //....oooOO0OOooo........oooOO0OOooo........oo 435 436 void Analysis::WritePhysGeo() 437 { 438 //Create a file containing geopath and Phy 439 std::string fname = "imp.info"; 440 std::ofstream fout(fname); 441 fout<<"====> Auto_generated file. Do not d 442 fout<<"====> The file conveys some informa 443 fout<<"_geocellpath "<<fCellDefFilePath<<" 444 for (const auto & entry : fVoxelDefFilesLi 445 fout<<"_geovolxelpath "<<std::string(e 446 } 447 fout<<"_numberOfBasepairs "<<fTotalNbBpPla 448 fout<<"_numberOfHistones "<<fTotalNbHiston 449 fout<<"_nucleusVolume "<<fNucleusVolume<<" 450 fout<<"_nucleusMassDensity "<<fNucleusMass 451 fout<<"_nucleusMass "<<fNucleusVolume*fNuc 452 fout.close(); 453 } 454 455 //....oooOO0OOooo........oooOO0OOooo........oo 456 457 void Analysis::DefineCommands() 458 { 459 fMessenger = std::make_unique<G4GenericMes 460 "/dsbandrepair/ou 461 "cmd controld"); 462 auto & fChemOutFolderCmd = fMessenger->Dec 463 464 fChemOutFolderCmd.SetParameterName("outFol 465 fChemOutFolderCmd.SetDefaultValue("chem_ou 466 } 467 468 //....oooOO0OOooo........oooOO0OOooo........oo 469 //....oooOO0OOooo........oooOO0OOooo........oo 470 471 void Analysis::CheckAndCreateNewFolderInChemSt 472 { 473 G4fs::file_status myFs = G4fs::file_status 474 auto folderPath = G4fs::path{fChemOutFolde 475 auto isExist = G4fs::status_known(myFs) ? 476 if (isExist && !G4fs::is_empty(folderPath) 477 G4ExceptionDescription msg; 478 msg <<"==>> Chem output folder "<<fChe 479 <<" is already existing and not em 480 msg<<"==>> Please delete or rename it, 481 <<"Chem output folder in macro file!!! 482 G4Exception("Analysis::CheckAndCreateN 483 } 484 G4fs::create_directory(folderPath); 485 } 486 487 //....oooOO0OOooo........oooOO0OOooo........oo 488 489 void Analysis::CheckAndCreateNewFolderInPhysSt 490 { 491 G4fs::file_status myFs = G4fs::file_status 492 const G4fs::path folderPath{fPhysOutFolder 493 auto isExist = G4fs::status_known(myFs) ? 494 if (isExist && !G4fs::is_empty(folderPath) 495 G4fs::remove_all(folderPath); 496 } 497 G4fs::create_directory(folderPath); 498 499 const G4fs::path folderPathPC{fChemInputFo 500 isExist = G4fs::status_known(myFs) ? G4fs: 501 if (isExist && !G4fs::is_empty(folderPathP 502 G4fs::remove_all(folderPathPC); 503 } 504 G4fs::create_directory(folderPathPC); 505 } 506 507 //....oooOO0OOooo........oooOO0OOooo........oo