dna/dsbandrepair/dsbandrepair.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 //  author: Le Tuan Anh, 20/10/2023
0027 /// \file main.cc
0028 /// \brief Main program of the dsbandrepair
0029
0030 #include "G4UImanager.hh"
0031 #include "G4UIterminal.hh"
0032 #include "G4UItcsh.hh"
0033 #include "G4UIExecutive.hh"
0034
0035 #include "G4RunManagerFactory.hh"
0036
0037 #ifdef G4VIS_USE
0038 #include "G4VisExecutive.hh"
0039 #endif
0040
0041 #include "G4Timer.hh"
0042 #include "G4ExceptionSeverity.hh"
0043 #include "G4DNAChemistryManager.hh"
0044 #include "G4VisExecutive.hh"
0045 #include "G4Filesystem.hh"
0046
0047 #include "ActionInitialization.hh"
0048 #include "DetectorConstruction.hh"
0049 #include "PhysicsList.hh"
0050
0051 #include "Analysis.hh"
0052
0053 #ifdef USE_MPI
0054 #include "G4MPImanager.hh"
0055 #include "G4MPIsession.hh"
0056 #include "G4MPIextraWorker.hh"
0057 #endif
0058 #include <ctime>
0059
0060 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0061
0062 G4String ExtractChemListNameFromMacroFile(G4String);
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0065
0066 int main(int argc,char** argv)
0067 {
0068 #ifdef USE_MPI
0069         G4MPImanager* g4MPI = new G4MPImanager(argc, argv, 0);
0070         g4MPI->SetVerbose(1);
0071         G4MPIsession* session = g4MPI-> GetMPIsession();
0072         G4String prompt = "[40;01;33m";
0073         prompt += "G4MPI";
0074         prompt += "[40;31m(%s)[40;36m[%/][00;30m:";
0075         session-> SetPrompt(prompt);
0076 #else
0077         G4UIExecutive* ui = nullptr;
0078         if ( argc == 1 ) { ui = new G4UIExecutive(argc, argv); }
0079 #endif // USE_MPI
0080     if (argc < 2) {
0081         G4cerr<<"====>> Wrong input. To run Physgeo, type : ./dsbandrepair macrofile\n"
0082             <<"To run Chem_geo, type : ./dsbandrepair macrofile chem"<<G4endl;
0083         #ifdef USE_MPI
0084             delete g4MPI;
0085         #endif // USE_MPI
0086         return EXIT_SUCCESS;
0087     }
0088     G4String stgstr = "";
0089     G4String macrofileName = argv[1];
0090     if (argc > 2) {
0091         const G4String rmode = argv[2];
0092         if (rmode == "phys") gRunMode = RunningMode::Phys;
0093         else if (rmode == "chem") gRunMode = RunningMode::Chem;
0094         else {
0095             G4cout<<"Undefined Running Mode = "<<rmode<<" ;dsbansrepair will quit now. See you!\n";
0096             #ifdef USE_MPI
0097                 delete g4MPI;
0098             #endif // USE_MPI
0099             return  EXIT_SUCCESS;
0100         }
0101     }
0102     // Choose the Random engine
0103     time_t timeStart;
0104     time(&timeStart);
0105     unsigned long seed = timeStart;
0106 #ifdef USE_MPI
0107     // Le Tuan Anh: add rankID to get different seeds for multi parallel processes
0108     seed += g4MPI->GetRank();
0109 #endif // USE_MPI
0110     G4cout<<"Initial Seed for random engine: "<<seed<<G4endl;
0111     CLHEP::HepRandom::setTheEngine(new CLHEP::MTwistEngine);
0112     CLHEP::HepRandom::setTheSeed(seed);
0113     G4RunManager* runManager{nullptr};
0114     auto analysis = Analysis::GetAnalysis();
0115     if (gRunMode == RunningMode::Phys) {
0116         stgstr = "physical stage";
0117         runManager = G4RunManagerFactory::CreateRunManager(G4RunManagerType::Default);
0118 #ifdef G4MULTITHREADED
0119         G4int threadNumber= 1;
0120         runManager-> SetNumberOfThreads(threadNumber);
0121 #endif  // G4MULTITHREADED
0122 #ifdef USE_MPI
0123         stgstr += " in #rank"+std::to_string(g4MPI->GetRank());
0124         if (g4MPI->IsMaster()) analysis->CheckAndCreateNewFolderInPhysStage();
0125 #else
0126         analysis->CheckAndCreateNewFolderInPhysStage();
0127 #endif
0128     } else if (gRunMode == RunningMode::Chem) {
0129         stgstr = "chemical stage";
0130         runManager = G4RunManagerFactory::CreateRunManager(G4RunManagerType::Serial);
0131         G4DNAChemistryManager::Instance()->SetChemistryActivation(true);
0132         G4Scheduler::Instance();
0133     }
0134
0135     DetectorConstruction* detector = new DetectorConstruction(1.,0,false);
0136     runManager->SetUserInitialization(detector);
0137     PhysicsList* physList = new PhysicsList;
0138     ActionInitialization* actionIni = new ActionInitialization();
0139
0140     if (gRunMode == RunningMode::Phys) {
0141         runManager->SetUserInitialization(physList);
0142         runManager->SetUserInitialization(actionIni);
0143 #ifdef USE_MPI
0144         session-> SessionStart();
0145         if (g4MPI->GetRank() == 0 ){
0146             analysis->WritePhysGeo();
0147         }
0148 #else
0149         // Get the pointer to the User Interface manager
0150         G4UImanager* UImanager = G4UImanager::GetUIpointer();
0151         // Process macro or start UI session
0152         if ( ! ui ) {
0153             // batch mode
0154             G4String command = "/control/execute ";
0155             UImanager->ApplyCommand(command+macrofileName);
0156         }
0157         analysis->WritePhysGeo();
0158 #endif // USE_MPI
0159     }
0160
0161
0162     if (gRunMode == RunningMode::Chem) {
0163         //get the pointer to the User Interface manager
0164         G4UImanager* UI = G4UImanager::GetUIpointer();
0165 #ifdef USE_MPI
0166         session->SessionStart();
0167         stgstr += " in #rank"+std::to_string(g4MPI->GetRank());
0168         if (g4MPI->IsMaster()) analysis->CheckAndCreateNewFolderInChemStage();
0169 #else
0170         G4String command = "/control/execute ";
0171         UI->ApplyCommand(command+macrofileName);
0172         analysis->CheckAndCreateNewFolderInChemStage();
0173 #endif
0174         //------------------------------------------
0175         // Prepare input file
0176         //------------------------------------------
0177         std::string inputFileorFolder = "chem_input";
0178         if (argc == 4) inputFileorFolder = argv[3];
0179         G4fs::path p{inputFileorFolder};
0180         G4String outputFileName = "test";
0181         std::vector<G4String> totalNumberofFilesVector, numberOfFilesTobeProcessedVector;
0182         if (G4fs::is_directory(p)) {
0183             for (const auto& entry : G4fs::directory_iterator(p)) {
0184                 if (entry.path().extension() == ".dat") {
0185                     totalNumberofFilesVector.push_back(entry.path().string());
0186                 }
0187             }
0188             std::sort(totalNumberofFilesVector.begin(),totalNumberofFilesVector.end());
0189 #ifdef USE_MPI
0190             G4int numberofRanks = g4MPI->GetActiveSize();
0191             size_t filesTobeProcessedSlave = (size_t)(
0192                 std::floor(G4double(totalNumberofFilesVector.size())/G4double(numberofRanks)));
0193             // note: should not use "std::ceil"
0194             size_t filesTobeProcessedMaster =
0195                 totalNumberofFilesVector.size() - (numberofRanks-1)*filesTobeProcessedSlave;
0196             if (g4MPI->IsMaster()) {
0197                 for (size_t ii=0; ii< filesTobeProcessedMaster; ii++) {
0198                     numberOfFilesTobeProcessedVector.push_back(totalNumberofFilesVector.at(ii));
0199                 }
0200             } else {
0201                 for (size_t ii=0; ii< filesTobeProcessedSlave; ii++) {
0202                     auto rankID = g4MPI->GetRank();
0203                     size_t kk = filesTobeProcessedMaster + (rankID-1)*filesTobeProcessedSlave + ii;
0204                     numberOfFilesTobeProcessedVector.push_back(totalNumberofFilesVector.at(kk));
0205                 }
0206             }
0207             G4cout<<"-----> "<<numberOfFilesTobeProcessedVector.size()
0208                     <<" files will be processed on rank #"<<g4MPI->GetRank()<<G4endl;
0209 #else
0210             numberOfFilesTobeProcessedVector = totalNumberofFilesVector;
0211 #endif
0212             if (totalNumberofFilesVector.size() == 0) {
0213                 G4cout<<"===>> There is no files found in "<<inputFileorFolder
0214                 <<". You have to run Phys_geo first!!!"<<G4endl;
0215                 #ifdef USE_MPI
0216                 delete g4MPI;
0217                 #endif // USE_MPI
0218                 delete runManager;
0219                 return EXIT_SUCCESS;
0220             } else {
0221                 G4cout<<"===>> Total files found in "<<inputFileorFolder
0222                 <<" : "<<totalNumberofFilesVector.size()<<G4endl;
0223             }
0224         } else if (G4fs::is_regular_file(p)) {
0225             numberOfFilesTobeProcessedVector.push_back(inputFileorFolder);
0226             if (p.has_stem()) {
0227                 outputFileName = p.stem().string();
0228             } else outputFileName = inputFileorFolder;
0229         }
0230         else G4cout<<"===>>dsbandrepair: "<<p.string()<<" is Not Directory or file !!!"<<G4endl;
0231         G4String firstFileForInit="";
0232         if (numberOfFilesTobeProcessedVector.size()>0) {
0233             firstFileForInit=numberOfFilesTobeProcessedVector.at(0);
0234             detector->ParseGeoFileForChemMode(firstFileForInit); // read to build voxel
0235         }
0236         //------------------------------------------
0237         // Initialization classes
0238         //------------------------------------------
0239
0240         runManager->SetUserInitialization(physList);
0241         runManager->SetUserInitialization(actionIni);
0242         runManager->Initialize();
0243         if (numberOfFilesTobeProcessedVector.size()>0) {
0244             size_t nprocessedfiles{0}, ncounts{1};
0245             if (numberOfFilesTobeProcessedVector.size()>=100) ncounts=10;
0246             if (numberOfFilesTobeProcessedVector.size()>=10000) ncounts=100;
0247             if (numberOfFilesTobeProcessedVector.size()>=1000000) ncounts=500;
0248             for (auto const &fileInput : numberOfFilesTobeProcessedVector) {
0249                 G4fs::path aP{std::string(fileInput)};
0250                 if (aP.has_stem()) {
0251                     outputFileName = aP.stem().string();
0252                 } else outputFileName = fileInput;
0253                 analysis->SetFileName(outputFileName);
0254                 if (fileInput != firstFileForInit) detector->ParseGeoFileForChemMode(fileInput);
0255                 detector->InsertMoleculeInWorld();
0256                 UI->ApplyCommand("/run/beamOn 1");
0257                 nprocessedfiles++;
0258                 if (nprocessedfiles == 1 ||
0259                     nprocessedfiles == numberOfFilesTobeProcessedVector.size() ||
0260                     0 == (nprocessedfiles % ncounts)) {
0261                     G4cout<<"=====> Processed file: "<<nprocessedfiles<<"-th/("
0262                           <<numberOfFilesTobeProcessedVector.size()<<" files)"
0263 #ifdef USE_MPI
0264                           <<" in rank #"<<g4MPI->GetRank()
0265 #endif
0266                           <<"!!!"<<G4endl;
0267                 }
0268             }
0269         } else {
0270             UI->ApplyCommand("/run/beamOn 1");
0271         }
0272     }
0273
0274 #ifdef USE_MPI
0275     delete g4MPI;
0276 #endif // USE_MPI
0277     delete runManager;
0278
0279     G4cout  <<"----------------------> Finish "<<stgstr<<"!!! Good bye :) <----------------------"<<G4endl;
0280     return EXIT_SUCCESS;
0281 }
0282
0283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......