GGSPenny.cpp

/*
 * GGSPenny.cpp
 *
 *  Created on: 2010-09-29
 *     Authors: Emiliano Mocchiutti and Cecilia Pizzolotto
 */
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// C++ and STL headers
#include <sys/stat.h>
#include <vector>

// ROOT headers
#include "TROOT.h"

// GEANT4 headers
#include "G4PhysListFactory.hh"
#include "G4UImanager.hh"
#include "G4Version.hh"
#include "Randomize.hh"

#ifdef G4VIS_USE
#include "G4VisExecutive.hh"
#endif

#ifdef G4UI_USE
#include "G4UIExecutive.hh"
#endif

// GGS headers
#include "GGSVersion.h" //Automatically generated by the CMake configuration system
#include "montecarlo/mccore/GGSDetectorConstruction.h"
#include "montecarlo/mccore/GGSMTRunManager.h"
#include "montecarlo/mccore/GGSMultiUserAction.h"
#include "montecarlo/mccore/GGSRunManager.h"
#include "montecarlo/mccore/GGSUserActionInitialization.h"
#include "montecarlo/mccore/GGSUserWorkerInitialization.h"
#include "montecarlo/mccore/GGSUserWorkerThreadInitialization.h"
#include "montecarlo/pluginmanagers/GGSMCPluginManager.h"
#include "montecarlo/services/GGSRootFileService.h"
#include "utils/GGSSmartLog.h"

#ifdef CRMC
#include "montecarlo/physicslists/GGSCRMCPhysicsList.h"
#endif

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// clang-format off
void usage() {
  std::cout << "Usage: GGSPenny [options] \n";
  std::cout << "This is the GGS main MonteCarlo program\n\n";

  std::cout << "Run type selection\n";
  std::cout << "------------------\n";
  std::cout << " -X, --interactive                interactive mode [default: false]\n\n";

  std::cout << " -d, --data-card FILE             data card filename [default: not set]\n\n";

  std::cout << "If only -d is specified  the program will run in batch mode.\n";
  std::cout << "If both options are specified, the datacard will be processed before entering the interactive mode.\n\n";

  std::cout << "Configuration options\n";
  std::cout << "---------------------\n";
  std::cout << " -g, --geometry FILE              full path and name of the geometry; can be a .so or a .gdml.\n";
  std::cout << " -gd, --geo-data-card FILE        data card filename for geometry configuration [default: not set]\n\n";
  std::cout << " -pl, --physics-list NAME         physics list (standard list name or full path of physics list .so library) [default: QGSP_BERT]\n\n";
  std::cout << " -pd, --physlist-data-card FILE   data card filename for physics list configuration [default: not set]\n\n";
  std::cout << " -ap, --action-plugin FILE        full path and name of the plugin user and generator action library\n";
  std::cout << "                                  May be repeated to load many libraries: -ap libUserAction.so -ap libGeneratorAction.so\n\n";
  std::cout << " -st, --single-threaded           Force single-threaded run mode for multi-threaded build\n\n";
  std::cout << " -nt, --n-threads                 Set the number of threads for multi-threaded run [default: 2]\n\n";
  std::cout << " -ro, --root-output-file FILE     default ROOT output filename (may be overridden in datacard for each action)\n\n";
  std::cout << " -os, --output-suffix SUFFIX      output file suffix [default: ]\n\n";
  std::cout << " --seed1 INT                      first random seed [default: 1]\n\n"; // do not put default 0!
  std::cout << " --seed2 INT                      second random seed [default: 2]\n\n";
  std::cout << " -s, --input-seed-file FILE       file containing input random seeds [default: not set]\n\n";
  std::cout << "Configuration options can either be given as command line parameters or placed in an input file.\n";
  std::cout << "File syntax is the same as the above specified commands, one command per line. Set the input file via command line as:\n\n";
  std::cout << " -i, --input-file FILE            file containing input arguments [default: not set]\n\n";

  std::cout << "Miscellaneous\n";
  std::cout << "-------------\n";
  std::cout << " -v, --verbose                    be verbose [default: false]\n\n";
  std::cout << " -V, --version                    displays software version and exit\n\n";
  std::cout << " -h, --help                       display this help and exit\n\n";

  std::cout << "\nTo run GGSPenny either the -d (data card) or the -X (interactive) option MUST be set!\n";
  std::cout << "Report bugs to: https://baltig.infn.it/mori/GGSSoftware/-/issues \n";
}
// clang-format on

void GGSPennyVersion() {
  std::cout << "GGSPenny simulation software";
#ifdef GGSMULTITHREADED
  std::cout << " (MT build)\n";
#else
  std::cout << " (ST build)\n";
#endif
  std::cout << "Using GGS " << GGSVERSION << std::endl;
}

namespace GGSPenny {
G4String rootOutputFile = "";
G4String outputSuffix = "";
G4String dataCard = "";
G4String geoDataCard = "";
G4String physListDataCard = "";
G4String inputSeedFile = "";
G4String geometry = "";
G4String physicsListName = "FTFP_BERT";
std::vector<G4String> actionPluginFiles;
G4bool interactiveMode = false;
#ifdef GGSMULTITHREADED
G4bool multiThreadedMode = true;
G4int nThreads = 2;
#else
G4bool multiThreadedMode = false;
G4int nThreads = 1;
#endif
G4long seed1 = 1;
G4long seed2 = 2;
} // namespace GGSPenny
using namespace GGSPenny;

void DumpGlobals() {
  std::cout << "==============> MonteCarlo globals used for this run <===============\n";
  std::cout << "dataCard          : " << dataCard.c_str() << "\n";
  std::cout << "geometry          : " << geometry.c_str() << "\n";
  std::cout << "geoDataCard       : " << geoDataCard.c_str() << "\n";
  std::cout << "physicsListName   : " << physicsListName.c_str() << "\n";
  std::cout << "physListDataCard  : " << physListDataCard.c_str() << "\n";
  std::cout << "rootOutputFile    : " << rootOutputFile.c_str() << "\n";
  std::cout << "outputSuffix      : " << outputSuffix.c_str() << "\n";
  std::cout << "inputSeedFile     : " << inputSeedFile.c_str() << "\n";
  std::cout << "seed1             : " << seed1 << "\n";
  std::cout << "seed2             : " << seed2 << "\n";
  std::cout << "interactiveMode   : " << interactiveMode << "\n";
  std::cout << "multiThreadedMode : " << multiThreadedMode << "\n";
  if (multiThreadedMode) {
    std::cout << "nThreads          : " << nThreads << "\n";
  }
  std::cout << "======================================================================" << std::endl;
}

int HandleInputPar(int argc, char **argv) {
  static const std::string routineName("HandleInputPar");
  // Set verbose level to print warnings about unrecognized options
  GGSSmartLog::verboseLevel = GGSSmartLog::WARNING;
  int vLevel = GGSSmartLog::INFO;
  if (argc > 1) {
    int i = 1;
    while (i < argc) {
      if (!strcmp(argv[i], "-h") || !strcmp(argv[i], "--help")) {
        usage();
        return 1;
      }
      // -----------------------------------------------------//
      if (!strcmp(argv[i], "-V") || !strcmp(argv[i], "--version")) {
        GGSPennyVersion();
        return 1;
      }
      // -----------------------------------------------------//
      if (!strcmp(argv[i], "--")) {
        return 0;
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-d") || !strcmp(argv[i], "--data-card")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        dataCard = argv[i];
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-gd") || !strcmp(argv[i], "--geo-data-card")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        geoDataCard = argv[i];
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-g") || !strcmp(argv[i], "--geometry")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        geometry = argv[i];
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-pl") || !strcmp(argv[i], "--physics-list")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        physicsListName = argv[i];
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-pd") || !strcmp(argv[i], "--physlist-data-card")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        physListDataCard = argv[i];
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-ap") || !strcmp(argv[i], "--action-plugin")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        actionPluginFiles.push_back(G4String(argv[i]));
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-st") || !strcmp(argv[i], "--single-threaded")) {
        multiThreadedMode = false;
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-nt") || !strcmp(argv[i], "--n-threads")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        nThreads = std::stoi(argv[i]);
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-X") || !strcmp(argv[i], "--interactive")) {
        interactiveMode = true;
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-ro") || !strcmp(argv[i], "--root-output-file")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        rootOutputFile = argv[i];
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-os") || !strcmp(argv[i], "--output-suffix")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        outputSuffix = argv[i];
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-s") || !strcmp(argv[i], "--input-seed-file")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        inputSeedFile = argv[i];
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "--seed1")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        seed1 = std::stol(argv[i]);
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "--seed2")) {
        if (++i >= argc) {
          usage();
          return 1;
        }
        seed2 = std::stol(argv[i]);
      }
      // -----------------------------------------------------//
      else if (!strcmp(argv[i], "-v") || !strcmp(argv[i], "--verbose")) {
        std::string vLevStr(argv[++i]);
        if (vLevStr == "error")
          vLevel = GGSSmartLog::ERROR;
        else if (vLevStr == "warning")
          vLevel = GGSSmartLog::WARNING;
        else if (vLevStr == "normal" || vLevStr == "info")
          vLevel = GGSSmartLog::INFO;
        else if (vLevStr == "debug")
          vLevel = GGSSmartLog::DEBUG;
        else if (vLevStr == "deepdeb")
          vLevel = GGSSmartLog::DEEPDEB;
        else {
          GGSCOUT(WARNING) << "Unrecognized verbosity level: " << vLevStr << ". Ignoring." << GGSENDL;
        }
      }
      // -----------------------------------------------------//
      else {
        GGSCOUT(WARNING) << "Unidentified input argument \"" << argv[i] << "\". Ignoring." << std::endl;
      }
      i++;
    }
  } else {
    usage();
    return 1;
  }
  GGSSmartLog::verboseLevel = vLevel;
  return 0;
}

int main(int argc, char **argv) {

  static const std::string routineName("GGSPenny");

  if (HandleInputPar(argc, argv))
    return 0;

  GGSPennyVersion();

  if (!strcmp(dataCard.c_str(), "") && !interactiveMode) {
    usage();
    return 1;
  }

  if (geometry == "") {
    GGSCOUT(ERROR) << "No geometry provided. Exit." << GGSENDL;
    return 1;
  }

  // Check the existence of given files
  {
    struct stat st;
    if (dataCard != "" && stat(dataCard.c_str(), &st) != 0) {
      GGSCOUT(ERROR) << "Datacard " << dataCard << " not found." << GGSENDL;
      return 1;
    }
    if (geoDataCard != "" && stat(geoDataCard.c_str(), &st) != 0) {
      GGSCOUT(ERROR) << "Geometry datacard " << geoDataCard << " not found." << GGSENDL;
      return 1;
    }
    if (physListDataCard != "" && stat(physListDataCard.c_str(), &st) != 0) {
      GGSCOUT(ERROR) << "Physics list datacard " << physListDataCard << " not found." << GGSENDL;
      return 1;
    }
  }

  // print on the screen input to this run
  DumpGlobals();

  if (multiThreadedMode) {
    ROOT::EnableThreadSafety();
  }

  GGSCOUT(DEBUG) << "Initialize random generator." << GGSENDL;
  // Choose the Random engine
  CLHEP::MixMaxRng rndmEngine;
  G4Random::setTheEngine(&rndmEngine);
  G4long seeds[2] = {seed1, seed2};
  G4Random::setTheSeeds(seeds, 0);
  G4Random::showEngineStatus();

  // Construct the default run manager
  //
  GGSCOUT(DEBUG) << "Create the run manager." << GGSENDL;
  G4RunManager *runManager = nullptr;
  if (multiThreadedMode) {
    runManager = new GGSMTRunManager;
    static_cast<GGSMTRunManager *>(runManager)->SetNumberOfThreads(nThreads);
  } else {
    runManager = new GGSRunManager;
  }
  // Set mandatory initialization classes
  GGSCOUT(DEBUG) << "Create the detector." << GGSENDL;
  GGSDetectorConstruction *detector = nullptr;
  if (geometry.length() > 4 && geometry.substr(geometry.length() - 5) == ".gdml") {
    detector = new GGSDetectorConstruction(geometry);
  }
#ifdef GGS_USE_DD4HEP
  else if (geometry.length() > 3 && geometry.substr(geometry.length() - 4) == ".xml") {
    detector = new GGSDetectorConstruction(geometry);
  }
#endif
  else {
    detector = new GGSDetectorConstruction(geometry, geoDataCard);
  }
  runManager->SetUserInitialization(detector);

  G4VUserPhysicsList *physicsList = nullptr;
  // Try to create reference physics list
  G4PhysListFactory *physListFactory = new G4PhysListFactory;
  if (physListFactory->IsReferencePhysList(physicsListName)) {
    physicsList = physListFactory->GetReferencePhysList(physicsListName);
  }
  // Try to create standard GGS physics lists
  if (physicsList == nullptr) {
#ifdef CRMC
    if (physicsListName == "GGSCRMC") {
      physicsList = new GGSCRMCPhysicsList;
    }
#endif
  }
  // Try to create user-defined physics list
  if (physicsList == nullptr) {
    // Try to load user defined list
    physicsList = GGSMCPluginManager::GetInstance().LoadPhysicsListPlugin(physicsListName);
    if (physicsList == nullptr) {
      GGSCOUT(WARNING) << "User physics list has NOT been loaded. Fall back to standard list." << GGSENDL;
    }
  }
  // Build standard list if no list have been selected or library loading wasn't successful
  if (physicsList == nullptr) {
    GGSCOUT(INFO) << "Load standard physics list: FTFP_BERT." << GGSENDL;
    physicsListName = "FTFP_BERT";
    physicsList = physListFactory->GetReferencePhysList(physicsListName);
  }

  // Get the pointer to the User Interface manager
  G4UImanager *UImanager = G4UImanager::GetUIpointer();

  // Execute physics list data card
  if (physListDataCard != "") {
    GGSCOUT(DEBUG) << "Execute physics list data card" << GGSENDL;
    G4String command = "/control/execute ";
    if (UImanager->ApplyCommand(command + physListDataCard) != 0) {
      GGSCOUT(ERROR) << "Error in executing physics list datacard " << physListDataCard << GGSENDL;
      return 1;
    }
  }

  runManager->SetUserInitialization(physicsList);

  // Load user and generator actions plugin libraries
  // Must be done before creating the GGSMultiUserAction
  for (unsigned int i = 0; i < actionPluginFiles.size(); i++)
    GGSMCPluginManager::GetInstance().LoadPlugin(actionPluginFiles[i]);

  // Set user action initialization
  if (multiThreadedMode) {
    GGSCOUT(DEBUG) << "Set the user action initialization." << GGSENDL;
    static_cast<GGSMTRunManager *>(runManager)->SetUserInitialization(new GGSUserActionInitialization);
  } else {
    GGSCOUT(DEBUG) << "Set the user action." << GGSENDL;
    GGSMultiUserAction *muAction = new GGSMultiUserAction;
    runManager->SetUserAction((G4UserSteppingAction *)muAction);
    runManager->SetUserAction((G4UserTrackingAction *)muAction);
    runManager->SetUserAction((G4UserEventAction *)muAction);
    runManager->SetUserAction((G4UserRunAction *)muAction);
    runManager->SetUserAction((G4UserStackingAction *)muAction);
  }

  // Set default name, suffix of ROOT files
  GGSCOUT(DEBUG) << "Initialize output." << GGSENDL;
  if (rootOutputFile != "")
    GGSRootFileService::GetInstance().SetDefaultFileBase(rootOutputFile.data());
  GGSRootFileService::GetInstance().SetSuffix(outputSuffix);
  // Set the MC info to be saved on ROOT files
  GGSTSimInfo simInfo;
  simInfo.SetName("GGSSimInfo");
  simInfo.G4Version = G4Version.data();
  simInfo.GGSVersion = GGSVERSION;
  simInfo.dataCard = dataCard.data();
  simInfo.geometry = geometry.data();
  simInfo.geoDataCard = geoDataCard.data();
  simInfo.physicsList = physicsListName.data();
  simInfo.seed1 = seed1;
  simInfo.seed2 = seed2;
  GGSRootFileService::GetInstance().SetSimInfo(simInfo);

  if (multiThreadedMode) {
    static_cast<GGSMTRunManager *>(runManager)->SetUserInitialization(new GGSUserWorkerInitialization);
    static_cast<GGSMTRunManager *>(runManager)->SetUserInitialization(new GGSUserWorkerThreadInitialization);
  }

  // Initialize G4 kernel
  GGSCOUT(DEBUG) << "Initialize run manager." << GGSENDL;
  try {
    runManager->Initialize();
  } catch (std::runtime_error &exc) {
    GGSCOUT(ERROR) << "Error while building the geometry: " << exc.what() << GGSENDL;
    delete runManager;
    delete physListFactory;
    return 1;
  } catch (...) {
    GGSCOUT(ERROR) << "Impossible to set up the geometry due to an unknown exception. Exit." << GGSENDL;
    delete runManager;
    delete physListFactory;
    return 1;
  }

#ifdef G4VIS_USE
  // Create the ROOT cling interpreter before opening the graphical UI. This is needed to avoid
  // loading LLVM symbols (e.g. when using open source Radeon driver for AMD chips) before
  // creating the interpreter, which would cause a segmentation fault.
  gROOT->GetInterpreter();
  // Initialize visualization
  //
  G4VisManager *visManager = new G4VisExecutive;
  visManager->Initialize();
#endif

  if (dataCard != "") { // batch or configured interactive mode
    if (!interactiveMode) {
      GGSCOUT(DEBUG) << "Start batch mode." << GGSENDL;
    }
    G4String command = "/control/execute ";
    if (UImanager->ApplyCommand(command + dataCard) != 0) {
      GGSCOUT(ERROR) << "Error in executing datacard " << dataCard << GGSENDL;
      return 1;
    }
  }
  if (interactiveMode) { // interactive mode : define UI session
    if (dataCard != "") {
      GGSCOUT(DEBUG) << "Start configured interactive mode." << GGSENDL;
    } else {
      GGSCOUT(DEBUG) << "Start interactive mode." << GGSENDL;
    }

#ifdef G4UI_USE
    G4UIExecutive *ui = new G4UIExecutive(argc, argv);
    struct stat st;
#ifdef G4VIS_USE
    if (stat("vis.mac", &st) == 0) {
      if (UImanager->ApplyCommand("/control/execute vis.mac") != 0) {
        GGSCOUT(ERROR) << "Error in executing vis.mac visualization macro." << GGSENDL;
        return 1;
      }
    }
#endif
    if (ui->IsGUI()) {
      if (stat("visTutor/gui.mac", &st) == 0) {
        if (UImanager->ApplyCommand("/control/execute visTutor/gui.mac") != 0) {
          GGSCOUT(ERROR) << "Error in executing visTutor/gui.mac visualization macro." << GGSENDL;
          return 1;
        }
      }
    }
    ui->SessionStart();
    delete ui;
#endif
  }

  GGSCOUT(DEBUG) << "Job termination." << GGSENDL;

#ifdef G4VIS_USE
  delete visManager;
#endif
  delete runManager;

  delete physListFactory;

  return 0;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......