gemc.cc

/// \mainpage
/// \htmlonly <center><img src="gemc_logo.gif" width="230"></center>\endhtmlonly
/// \section overview Overview
/// gemc (<b>GE</b>ant4 <b>M</b>onte<b>C</b>arlo) GEMC is a C++ framework
/// based on <a href="http://geant4.web.cern.ch/geant4/"> Geant4 </a>
/// Libraries to simulate the passage of particles through matter.\n
/// The simulation parameters are external to the software:
/// Geometry, Materials, Fields, Banks definitions are stored in
/// external databases in various format and are loaded at run
/// time using factory methods.\n
/// \section databases Databases
/// gemc currently supports <i> mysql, gdml, TEXT </i> to build the detector systems: \n
/// - Geometry.
/// - Sensitive Detectors.
/// - Hit Process.
/// - Banks Format.
/// - Materials.
/// \section platforms Platforms Supported:
/// - Linux (32, 64)
/// - Mac OS X
/// \section docs Documentation:
/// - <a href="http://gemc.jlab.org">  gemc website </a>
/// \n\n
/// \author \n &copy; Maurizio Ungaro
/// \author e-mail: ungaro@jlab.org\n\n\n
/// \file gemc.cc
/// Defines the gemc main( int argc, char **argv )
/// \author \n &copy; Maurizio Ungaro
/// \author e-mail: ungaro@jlab.org\n\n\n

const char *GEMC_VERSION = "gemc 2.12";

// G4 headers
#include "G4RunManager.hh"
#include "G4UImanager.hh"
#include "G4UIterminal.hh"
#include "G4VisExecutive.hh"
#include "G4VModularPhysicsList.hh"
#include "G4PropagatorInField.hh"
#include "G4TransportationManager.hh"
#include "G4UIQt.hh"
#include "G4Qt.hh"

// Qt headers
#include <QApplication>
#include <QSplashScreen>

// gemc headers
#include "HitProcess_MapRegister.h"
#include "detector_factory.h"
#include "gemc_MainGui.h"
#include "gbank.h"
#include "MDetectorConstruction.h"
#include "outputFactory.h"
#include "HitProcess.h"
#include "PhysicsList.h"
#include "gemcOptions.h"
#include "dmesg_init.h"
#include "run_conditions.h"
#include "fieldFactory.h"
#include "material_factory.h"
#include "mirrors_factory.h"
#include "parameter_factory.h"
#include "string_utilities.h"
#include "gemcUtils.h"
#include "ActionInitialization.h"

// c++ headers
#include <unistd.h>  // needed for get_pid

/////////////////////////
/// <b> Main Program </b>
/////////////////////////
///  -# Sets the goptions\n
///  -# Starts QT application if USE_GUI=1 or 2
///  -# Starts the CLHEP random engine
///  -# Instantiates the Geant4 Run Manager
///  -# Builds detector map object from database
///  -# Builds Processes Routines Map
///  -# Builds Materials Map
///  -# Builds G4 Physical Volumes
///  -# Initialize Physics List
///  -# Initialize Generator, Event and Stepping Actions
///  -# Initialize G4Qt User Interface if USE_GUI>0
///  -# Initialize Visualization Manager if USE_GUI>0


// get_pid is useful only on the farm to set the seed
// can set to zero in Windows environment
// ideally we'd want __get_pid();
#ifdef _MSC_VER
#include <stdio.h>
#include <process.h>
//int get_pid(){return __get_pid();}
int get_pid(){return 0;}
#endif

// distinguishing between graphical and batch mode
QCoreApplication *createApplication(int &argc, char *argv[], double use_gui) {
    if (!use_gui)
        return new QCoreApplication(argc, argv);
    return new QApplication(argc, argv);
}


int main(int argc, char **argv) {
    clock_t startTime = clock();
    cout << endl;

    goptions gemcOpt;
    gemcOpt.setGoptions();
    gemcOpt.setOptMap(argc, argv, GEMC_VERSION);

    double use_gui = gemcOpt.optMap["USE_GUI"].arg;

    cout << endl << "  > Initializing GEant4 MonteCarlo:  " << GEMC_VERSION << endl << endl;

    QScopedPointer <QCoreApplication> app(createApplication(argc, argv, use_gui));

    // Initializing gemc splash class
    // This class will log initialization messages
    // This class will show a splashscreen if use_gui is non zero
    // The screen log verbosity is controlled by LOG_VERBOSITY
    gui_splash gemc_splash(gemcOpt);
    gemc_splash.message(" Initializing GEant4 MonteCarlo version " + string(GEMC_VERSION));


    // random seed initialization
    G4Random::setTheEngine(new CLHEP::MixMaxRng);

    G4int seed;

    if (gemcOpt.optMap["RANDOM"].args == "TIME") {
        gemc_splash.message(" Initializing CLHEP Random Engine from local time " \
 + stringify((double) time(nullptr)) \
 + ", cpu clock "        \
 + stringify((double) clock())    \
 + " and process id "    \
 + stringify(getpid()) + ".");
        seed = (G4int)((double) time(nullptr) - (double) clock() - getpid());
    } else {
        seed = atoi(gemcOpt.optMap["RANDOM"].args.c_str());
        gemc_splash.message(" Initializing CLHEP Random Engine from user defined seed.");
    }

    CLHEP::HepRandom::setTheSeed(seed);
    gemc_splash.message(" Seed initialized to: " + stringify(seed));

    // Construct the default G4 run manager
    gemc_splash.message(" Instantiating Run Manager...");
    G4RunManager *runManager = new G4RunManager;

    // Initializing run_condition class
    gemc_splash.message(" Instantiating Run Conditions...");
    runConditions runConds(gemcOpt);

    // GEMC Detector Map
    gemc_splash.message(" Registering Detectors Factories...");
    // Initializing Detector Factory
    map <string, detectorFactoryInMap> detectorFactoryMap = registerDetectorFactory();

    // Building detector with factories
    map <string, detector> hallMap = buildDetector(detectorFactoryMap, gemcOpt, runConds);

    // Initialize Materials Map Factory
    gemc_splash.message(" Initializing Material Factories...");
    map <string, materialFactory> materialFactoriesMap = registerMaterialFactories();
    // Build all materials
    map < string, G4Material * > mats = buildMaterials(materialFactoriesMap, gemcOpt, runConds);

    // Initialize Mirrors Map Factory
    gemc_splash.message(" Initializing Mirrors Factories...");
    map <string, mirrorFactory> mirrorFactoriesMap = registerMirrorFactories();
    // Build all mirrors
    map < string, mirror * > mirs = buildMirrors(mirrorFactoriesMap, gemcOpt, runConds);

    // Initialize Parameters Map Factory
    gemc_splash.message(" Registering Parameters Factories...");
    map <string, parameterFactoryInMap> parameterFactoriesMap = registerParameterFactories();
    // All Parameters with factories
    map<string, double> gParameters = loadAllParameters(parameterFactoriesMap, gemcOpt, runConds);

    // Process Hit Map
    gemc_splash.message(" Building gemc Process Hit Factory...");
    map <string, HitProcess_Factory> hitProcessMap = HitProcess_Map(gemcOpt.optMap["HIT_PROCESS_LIST"].args);

    ///< magnetic Field Map
    gemc_splash.message(" Creating fields Map...");
    map <string, fieldFactoryInMap> fieldFactoryMap = registerFieldFactories();
    map <string, gfield> fieldsMap = loadAllFields(fieldFactoryMap, gemcOpt);

    // Build the detector
    gemc_splash.message(" Building Detector Map...");
    MDetectorConstruction *ExpHall = new MDetectorConstruction(gemcOpt);
    ExpHall->hallMap = &hallMap;
    ExpHall->mirs = &mirs;
    ExpHall->mats = &mats;
    ExpHall->fieldsMap = &fieldsMap;
    // this is what calls Construct inside MDetectorConstruction
    runManager->SetUserInitialization(ExpHall);


    ///< Physics List
    string phys_list = gemcOpt.optMap["PHYSICS"].args;
    gemc_splash.message(" Initializing Physics List " + phys_list + "...");
    runManager->SetUserInitialization(new PhysicsList(gemcOpt));

    // Setting Max step for all the simulation.
    // Notice: on the forum:
    // http://hypernews.slac.stanford.edu/HyperNews/geant4/get/emfields/183/1.html
    // it is mentioned that going through volumes of different materials could create problems.
    // this is verified in clas12 when going from target to "hall" - even when vacuum was not involved.
    // the solution to that was to create a transitional tube from target to the vacuum line
    // This solution allowed to avoid setting MAX_FIELD_STEP to a value that would slow down the
    // simulation by a factor of 5

    double max_step = gemcOpt.optMap["MAX_FIELD_STEP"].arg;
    if (max_step != 0) {
        G4TransportationManager::GetTransportationManager()->GetPropagatorInField()->SetLargestAcceptableStep(max_step);
    }

    // User action initialization
    gemc_splash.message(" Initializing User Actions...");
    ActionInitialization *gActions = new ActionInitialization(&gemcOpt, &gParameters);
    runManager->SetUserInitialization(gActions);

    ///< User Interface manager
    gemc_splash.message(" Initializing User Interface...");

    ///< Vis Manager
    G4VisManager *visManager = nullptr;
    if (use_gui) {
        // G4VisExecutive can take a verbosity argument - see /vis/verbose guidance.
        visManager = new G4VisExecutive("Quiet");
        visManager->Initialize();
    }

    G4UIsession *session = nullptr;
    if (use_gui) {
        session = new G4UIQt(1, argv);
    }


    // Output File: registering output type, output process factory,
    // sensitive detectors into Event Action
    gemc_splash.message(" Initializing Output Action...");
    outputContainer outContainer(gemcOpt);
    map <string, outputFactoryInMap> outputFactoryMap = registerOutputFactories();

    // Initialize G4 kernel
    gemc_splash.message(" Initializing Run Manager...\n");
    // physical volumes, sensitive detectors are built here
    runManager->Initialize();

    // registering activated field in the option so they're written out
    if (ExpHall->activeFields.size()) {
        gemcOpt.optMap["ACTIVEFIELDS"].args = "";
    }

    for (set<string>::iterator fit = ExpHall->activeFields.begin(); fit != ExpHall->activeFields.end(); fit++) {
        gemcOpt.optMap["ACTIVEFIELDS"].args = gemcOpt.optMap["ACTIVEFIELDS"].args + *fit + " ";
    }

    // Bank Map, derived from sensitive detector map
    gemc_splash.message(" Creating gemc Banks Map...");
    map <string, gBank> banksMap = read_banks(gemcOpt, runConds.get_systems());

    // Getting UI manager, restoring G4Out to cout
    G4UImanager *UImanager = G4UImanager::GetUIpointer();
    UImanager->SetCoutDestination(nullptr);

    // saving simulation condition in the output file
    if (outContainer.outType != "no") {
        // Creating the sim_condition map to save to the output
        gemc_splash.message(" Writing simulation parameters in the output...");

        // filling gcard option content
        map <string, string> sim_condition = gemcOpt.getOptMap();

        // adding detectors conditions to sim_condition
        mergeMaps(sim_condition, runConds.getDetectorConditionsMap());

        // adding parameters value to sim_condition
        mergeMaps(sim_condition, getParametersMap(gParameters));

        sim_condition["JSON"] = gemcOpt.jSonOptions();

        outputFactory *processOutputFactory = getOutputFactory(&outputFactoryMap, outContainer.outType);
        processOutputFactory->recordSimConditions(&outContainer, sim_condition);

        // then deleting process output pointer, not needed anymore
        delete processOutputFactory;
    }

    gActions->evtAction->outContainer = &outContainer;
    gActions->evtAction->outputFactoryMap = &outputFactoryMap;
    gActions->evtAction->hitProcessMap = &hitProcessMap;
    gActions->evtAction->SeDe_Map = ExpHall->SeDe_Map;
    gActions->evtAction->banksMap = &banksMap;
    gActions->evtAction->gen_action = gActions->genAction;

    ///< passing output process factory to sensitive detectors
    map<string, sensitiveDetector *>::iterator it;
    for (it = ExpHall->SeDe_Map.begin(); it != ExpHall->SeDe_Map.end(); it++) {
        it->second->hitProcessMap = &hitProcessMap;
    }

    gemc_splash.message(" Executing initial directives...\n");
    vector <string> init_commands = init_dmesg(gemcOpt);
    for (unsigned int i = 0; i < init_commands.size(); i++)
        UImanager->ApplyCommand(init_commands[i].c_str());
    string exec_macro = "/control/execute " + gemcOpt.optMap["EXEC_MACRO"].args;

    clock_t start_events;

    long int nEventsToProcess = gemcOpt.optMap["N"].arg;

    // if it is not set explicitely, and it is a file input, then run all the event in the file by default
    // only in batch mode
    if (use_gui == 0 && gemcOpt.optMap["N"].arg == 0 && gemcOpt.optMap["INPUT_GEN_FILE"].args != "gemc_internal") {
        nEventsToProcess = 1000000000;
    }

    if (use_gui) {
        gemc_splash.message("Starting GUI...");
        qApp->processEvents();

        gemcMainWidget gemcW(&gemcOpt, ExpHall->SeDe_Map, &hallMap, mats);
        gemcW.setWindowTitle(GEMC_VERSION);
        vector <string> wpos = get_info(gemcOpt.optMap["GUIPOS"].args);

        gemcW.move(get_number(wpos[0]), get_number(wpos[1]));
        gemcW.show();

        // splash can finish once gemcW is up
        gemc_splash.splash->finish(&gemcW);

        gemc_splash.message(" Executing initial visual directives...\n");
        vector <string> init_vcommands = init_dvmesg(gemcOpt, visManager);
        for (unsigned int i = 0; i < init_vcommands.size(); i++) {
            gemc_splash.message(" Now executing: " + init_vcommands[i]);

            UImanager->ApplyCommand(init_vcommands[i].c_str());
        }

        if (exec_macro != "/control/execute no") {
            UImanager->ApplyCommand(exec_macro.c_str());
        }

        if (nEventsToProcess > 0) {
            start_events = clock();
            char command[100];
            // starting clock after the first event is much more precise
            if (nEventsToProcess > 10) {
                snprintf(command, 100, "/run/beamOn 1");
                UImanager->ApplyCommand(command);
                start_events = clock();
                nEventsToProcess--;
            }
            snprintf(command, 100, "/run/beamOn %ld", nEventsToProcess);
            UImanager->ApplyCommand(command);
        }

        return qApp->exec();
        // deleting runManager is now taken care
        // in the gemc_quit slot
        delete visManager;
        if (session != nullptr) delete session;
    } else {
        if (exec_macro != "/control/execute no") UImanager->ApplyCommand(exec_macro.c_str());
        start_events = clock();
        if (nEventsToProcess > 0) {
            start_events = clock();
            char command[100];
            // starting clock after the first event is much more precise
            if (nEventsToProcess > 10) {
                snprintf(command, 100, "/run/beamOn 1");
                UImanager->ApplyCommand(command);
                start_events = clock();
                nEventsToProcess--;
            }
            snprintf(command, 100, "/run/beamOn %ld", nEventsToProcess);
            UImanager->ApplyCommand(command);
        }
    }

    clock_t endTime = clock();
    clock_t clockAllTaken = endTime - startTime;
    clock_t clockEventTaken = endTime - start_events;

    if (nEventsToProcess > 10) {
        clockEventTaken = clockEventTaken * (nEventsToProcess + 1) / nEventsToProcess;
    }

    cout << " > Total gemc time: " << clockAllTaken / (double) CLOCKS_PER_SEC << " seconds. "
         << " Events only time: " << clockEventTaken / (double) CLOCKS_PER_SEC << " seconds. " << endl;

    // closing db connection
    closeGdb();

    delete runManager;
    return 0;
}


// introducing OPTICALPHOTONPID here to be semi-transparent to G4 changes
// this pid changed from 0 to -22 with geant4 10.7
int MHit::OPTICALPHOTONPID = -22;