/**
 * \file global_sim.C
 * \brief Macro runs simulation for "electron" or "muon" setups of CBM.
 * \author Andrey Lebedev <andrey.lebedev@gsi.de>
 * \date 2010
 **/
/*
#include "CbmCave.h"
#include "CbmEcal.h"
#include "CbmFieldMapSym2.h"
#include "CbmFieldPar.h"
#include "CbmMagnet.h"
#include "CbmMuch.h"
#include "CbmMvd.h"
#include "CbmPipe.h"
#include "CbmRich.h"
#include "CbmShield.h"
#include "CbmSts.h"
#include "CbmTarget.h"
#include "CbmTof.h"
#include "CbmTrd.h"
#include "CbmUnigenGenerator.h"

#include "FairDetector.h"
#include "FairModule.h"
#include "FairParRootFileIo.h"
#include "FairPlutoGenerator.h"
#include "FairPrimaryGenerator.h"
#include "FairRunSim.h"
#include "FairRuntimeDb.h"
#include "FairUrqmdGenerator.h"

#include "TRandom.h"
#include "TROOT.h"
#include "TStopwatch.h"
#include "TString.h"
*/
#include <iostream>
using std::cout;
using std::endl;

void global_sim(Int_t nEvents = 100)
{
    //gSystem->Setenv("G3_CONFIG", "/u/ekrebs/workspace/cbmroot/JUL13/gconfig_cus");
    //TTree::SetMaxTreeSize(90000000000);
    TString script = TString(gSystem->Getenv("LIT_SCRIPT"));

    // Specify "electron" or "muon" setup of CBM
    //TString setup = "muon";
    TString setup = "electron";

    TString urqmd  = "yes"; // If "yes" than UrQMD will be used as background
    TString unigen = "yes"; // If "yes" than CbmUnigenGenerator will be used instead of FairUrqmdGenerator
    TString pluto  = "yes"; // If "yes" PLUTO generator will be used

    // Files
    TString urqmdFile  = "/hera/cbm/prod/gen/urqmd/auau/25gev/centr/urqmd.auau.25gev.centr.00001.root"; // input UrQMD file
    TString dir        = "/hera/cbm/users/ekrebs/analysis/dielectron/auau/centr/25gev/1.0field_50.0z_mvd3_noms/rho0/epem"; // Directory for output simulation files
    TString mcFile     = dir + "/sim/mc.0001.root"; //MC file name
    TString parFile    = dir + "/par/param.0001.root"; //Parameter file name
    TString plutoFile  = "/hera/cbm/users/ekrebs/pluto/cktA/25gev/rho0/epem/pluto.auau.25gev.rho0.epem.0001.root";

    // Geometries
    TString caveGeom = "", targetGeom = "", pipeGeom = "", shieldGeom = "",
        mvdGeom = "", stsGeom = "", muchGeom = "", richGeom = "", trdGeom = "",
    tofGeom = "", ecalGeom = "", fieldMap = "", magnetGeom = "";
    if (setup == "muon") {
        caveGeom   = "cave.geo";
        targetGeom = "target/target_au_025mu.geo";
        pipeGeom   = "pipe/pipe_much.geo";
        shieldGeom = "shield_segmented_CPbFe_wo2last.geo";//"shield_standard.geo";
        mvdGeom    = "mvd/mvd_3_Nov_2012.geo";//"mvd/mvd_v07a.geo";
        stsGeom    = "sts/sts_v12b.geo.root";
        muchGeom   = "much/much/much_v12c.geo";
        trdGeom    = "";//"trd_muon_setup_new.geo";
        tofGeom    = "tof/tof_v13b.root";
        fieldMap   = "field_v12a";
        magnetGeom = "magnet/magnet_v12a.geo";
    } else if (setup == "electron") {
        caveGeom   = "cave.geo";
        targetGeom = "target/target_au_025mu.geo";
        pipeGeom   = "pipe/pipe_standard_short.geo";
        mvdGeom    = "mvd/mvd_3_Nov_2012.geo";//"mvd/mvd_v07a.geo";
        stsGeom    = "sts/sts_v12b.geo.root";
        richGeom   = "rich/rich_v08a.geo";
        trdGeom    = "trd/trd_v13g.geo.root";
        tofGeom    = "tof/tof_v13b.root";
        ecalGeom   = "";//"ecal_FastMC.geo";
        fieldMap   = "field_v12b";
        magnetGeom = "magnet/magnet_v12a.geo";
    }

    // -----   Magnetic field   -----------------------------------------------
    Double_t fieldZ = 50.; // field center z position
    Double_t fieldScale = 1.; // field scaling factor
    TString fieldCut("no");

    // If SCRIPT environment variable is set to "yes", i.e. macro is run via script
    if (script == "yes") {
        urqmdFile  = TString(gSystem->Getenv("LIT_URQMD_FILE"));
        mcFile     = TString(gSystem->Getenv("LIT_MC_FILE"));
        parFile    = TString(gSystem->Getenv("LIT_PAR_FILE"));
        plutoFile  = TString(gSystem->Getenv("LIT_PLUTO_FILE"));

        urqmd      = TString(gSystem->Getenv("LIT_URQMD"));
        unigen     = TString(gSystem->Getenv("LIT_UNIGEN"));
        pluto      = TString(gSystem->Getenv("LIT_PLUTO"));

        caveGeom   = TString(gSystem->Getenv("LIT_CAVE_GEOM"));
        targetGeom = TString(gSystem->Getenv("LIT_TARGET_GEOM"));
        pipeGeom   = TString(gSystem->Getenv("LIT_PIPE_GEOM"));
        shieldGeom = TString(gSystem->Getenv("LIT_SHIELD_GEOM"));
        mvdGeom    = TString(gSystem->Getenv("LIT_MVD_GEOM"));
        stsGeom    = TString(gSystem->Getenv("LIT_STS_GEOM"));
        muchGeom   = TString(gSystem->Getenv("LIT_MUCH_GEOM"));
        richGeom   = TString(gSystem->Getenv("LIT_RICH_GEOM"));
        trdGeom    = TString(gSystem->Getenv("LIT_TRD_GEOM"));
        tofGeom    = TString(gSystem->Getenv("LIT_TOF_GEOM"));
        ecalGeom   = TString(gSystem->Getenv("LIT_ECAL_GEOM"));
        fieldMap   = TString(gSystem->Getenv("LIT_FIELD_MAP"));
	magnetGeom = TString(gSystem->Getenv("LIT_MAGNET_GEOM"));

	fieldScale = TString(gSystem->Getenv("LIT_FIELD_MAP_SCALE")).Atof();
        fieldZ     = TString(gSystem->Getenv("LIT_FIELD_Z")).Atof();
        fieldCut   = TString(gSystem->Getenv("LIT_FIELD_CUT"));
    }

    TStopwatch timer;
    timer.Start();

    gROOT->LoadMacro("$VMCWORKDIR/macro/littrack/loadlibs.C");
    loadlibs();

    FairRunSim* run = new FairRunSim();
    run->SetName("TGeant3"); // Transport engine
    run->SetOutputFile(mcFile); // Output file
    FairRuntimeDb* rtdb = run->GetRuntimeDb();

    run->SetMaterials("media.geo"); // Materials
    //run->SetStoreTraj(kTRUE);

    if ( caveGeom != "" ) {
        FairModule* cave = new CbmCave("CAVE");
        cave->SetGeometryFileName(caveGeom);
        run->AddModule(cave);
    }

    if ( pipeGeom != "" ) {
        FairModule* pipe = new CbmPipe("PIPE");
        pipe->SetGeometryFileName(pipeGeom);
        run->AddModule(pipe);
    }

    if ( shieldGeom != "" ) {
        FairModule* shield = new CbmShield("SHIELD");
        shield->SetGeometryFileName(shieldGeom);
        run->AddModule(shield);
    }

    if ( targetGeom != "" ) {
        FairModule* target = new CbmTarget("Target");
        target->SetGeometryFileName(targetGeom);
        run->AddModule(target);
    }

    if ( magnetGeom != "" ) {
        FairModule* magnet = new CbmMagnet("MAGNET");
        magnet->SetGeometryFileName(magnetGeom);
        run->AddModule(magnet);
    }

    if ( mvdGeom != "" ) {
        FairDetector* mvd = new CbmMvd("MVD", kTRUE);
        mvd->SetGeometryFileName(mvdGeom);
        run->AddModule(mvd);
    }

    if ( stsGeom != "" ) {
        FairDetector* sts = new CbmSts("STS", kTRUE);
        sts->SetGeometryFileName(stsGeom);
        run->AddModule(sts);
    }

    if ( richGeom != "" ) {
        FairDetector* rich = new CbmRich("RICH", kTRUE);
        rich->SetGeometryFileName(richGeom);
        run->AddModule(rich);
    }

    if ( muchGeom != "" ) {
        FairDetector* much = new CbmMuch("MUCH", kTRUE);
        much->SetGeometryFileName(muchGeom);
        run->AddModule(much);
    }

    if ( trdGeom != "" ) {
        FairDetector* trd = new CbmTrd("TRD",kTRUE );
        trd->SetGeometryFileName(trdGeom);
        run->AddModule(trd);
    }

    if ( tofGeom != "" ) {
        FairDetector* tof = new CbmTof("TOF", kTRUE);
        tof->SetGeometryFileName(tofGeom);
        run->AddModule(tof);
    }

    if ( ecalGeom != "" ) {
        FairDetector* ecal = new CbmEcal("ECAL", kTRUE, ecalGeom.Data());
        run->AddModule(ecal);
    }
    // ------------------------------------------------------------------------

    // -----   Create magnetic field   ----------------------------------------
    CbmFieldMap* magField = new CbmFieldMapSym3(fieldMap);
    magField->SetPosition(0., 0., fieldZ);
    magField->SetScale(fieldScale);
    if(fieldCut.Contains("yes")) {
        cout<<"magField->SetFieldRegionNew(-1000., 1000., -1000., 1000., 13., 1000.)  "<<endl;
        magField->SetFieldRegionNew(-1000., 1000., -1000., 1000., 13., 1000.);
    }
    magField->Print();
    run->SetField(magField);
    // ------------------------------------------------------------------------

    // ------------------------------------------------------------------------
    FairPrimaryGenerator* primGen = new FairPrimaryGenerator();

    if (urqmd == "yes" && unigen == "yes") {
        CbmUnigenGenerator*  urqmdGen = new CbmUnigenGenerator(urqmdFile);
        primGen->AddGenerator(urqmdGen);
    }

    if (urqmd == "yes" && unigen != "yes") {
        FairUrqmdGenerator* urqmdGen = new FairUrqmdGenerator(urqmdFile);
        primGen->AddGenerator(urqmdGen);
    }

    if (pluto == "yes") {
        FairPlutoGenerator* plutoGen = new FairPlutoGenerator(plutoFile);
        primGen->AddGenerator(plutoGen);
    }


    run->SetGenerator(primGen);
    run->Init();

    if (pluto == "yes" && urqmd == "yes") {
        Float_t bratioEta[6];
        Int_t modeEta[6];

        TGeant3* gMC3 = (TGeant3*) gMC;

        for (Int_t kz = 0; kz < 6; ++kz) {
            bratioEta[kz] = 0.;
            modeEta[kz]   = 0;
        }

        Int_t ipa    = 17;
        bratioEta[0] = 39.38;  //2gamma
        bratioEta[1] = 32.20;  //3pi0
        bratioEta[2] = 22.70;  //pi+pi-pi0
        bratioEta[3] = 4.69;   //pi+pi-gamma
        bratioEta[4] = 0.60;   //e+e-gamma
        bratioEta[5] = 4.4e-2; //pi02gamma

        modeEta[0] = 101;    //2gamma
        modeEta[1] = 70707;  //3pi0
        modeEta[2] = 80907;  //pi+pi-pi0
        modeEta[3] = 80901;  //pi+pi-gamma
        modeEta[4] = 30201;  //e+e-gamma
        modeEta[5] = 10107;  //pi02gamma
        gMC3->Gsdk(ipa, bratioEta, modeEta);

        Float_t bratioPi0[6];
        Int_t modePi0[6];

        for (Int_t kz = 0; kz < 6; ++kz) {
            bratioPi0[kz] = 0.;
            modePi0[kz] = 0;
        }

        ipa = 7;
        bratioPi0[0] = 98.798;
        bratioPi0[1] = 1.198;

        modePi0[0] = 101;
        modePi0[1] = 30201;

        gMC3->Gsdk(ipa, bratioPi0, modePi0);

        Int_t t = time(NULL);
        TRandom *rnd = new TRandom(t);
        gMC->SetRandom(rnd);
    }

    // -----   Runtime database   ---------------------------------------------
    CbmFieldPar* fieldPar = (CbmFieldPar*) rtdb->getContainer("CbmFieldPar");
    fieldPar->SetParameters(magField);
    fieldPar->setChanged();
    fieldPar->setInputVersion(run->GetRunId(),1);
    Bool_t kParameterMerged = kTRUE;
    FairParRootFileIo* parOut = new FairParRootFileIo(kParameterMerged);
    parOut->open(parFile.Data());
    rtdb->setOutput(parOut);
    rtdb->saveOutput();
    rtdb->print();
    // ------------------------------------------------------------------------

    // -----   Start run   ----------------------------------------------------
    run->Run(nEvents);
    // ------------------------------------------------------------------------

    // -----   Finish   -------------------------------------------------------
    timer.Stop();
    Double_t rtime = timer.RealTime();
    Double_t ctime = timer.CpuTime();
    cout << endl << endl;
    cout << "Macro finished successfully." << endl;
    cout << "Test passed"<< endl;
    cout << " All ok " << endl;
    cout << "Output file is " << mcFile << endl;
    cout << "Real time used: " << rtime << "s " << endl;
    cout << "CPU time used : " << ctime << "s " << endl << endl << endl;
    // ------------------------------------------------------------------------
}