#if !defined(__CINT__) && !defined(__CLING__)
#include "TString.h"
#include "TStopwatch.h"
#include "TROOT.h"
#include "TSystem.h"
#include "TMath.h"
#include "TRandom.h"
#include "TDatabasePDG.h"

#include "FairRunSim.h"
#include "FairRuntimeDb.h"
#include "FairParRootFileIo.h"
#include "FairTrajFilter.h"
#include "FairUrqmdGenerator.h"
#include "FairPrimaryGenerator.h"
#include "MpdLAQGSMGenerator.h"
#include "FairCave.h"
#include "FairPipe.h"
#include "FairMagnet.h"
#include "FairBoxGenerator.h"
#include "MpdUrqmdGenerator.h"

#include "TpcDetector.h"
#include "MpdEmc.h"
#include "MpdTof.h"
#include "MpdZdc.h"
#include "MpdFfd.h"
#include "MpdMultiField.h"
#include "MpdMultiFieldPar.h"
#include "MpdConstField.h"
#include "MpdFieldMap.h"
#include "MpdGetNumEvents.h"

#include <iostream>
using namespace std;
#endif

R__ADD_INCLUDE_PATH($VMCWORKDIR)
#include "macro/mpd/mpdloadlibs.C"
#include "macro/mpd/geometry_stage1.C"
//#include "macro/mpd/geometry_v2.C"

#define URQMD  // Choose generator: URQMD VHLLE FLUID PART ION BOX HSD LAQGSM HADGEN
#define GEANT3  // Choose: GEANT3 GEANT4

// inFile - input file with generator data, default: auau.09gev.mbias.98k.ftn14
// nStartEvent - for compatibility, any number
// nEvents - number of events to transport, default: 1
// flag_store_FairRadLenPoint
// FieldSwitcher: 0 - corresponds to the ConstantField (0, 0, 5) kG (It is used by default); 1 - corresponds to the FieldMap ($VMCWORKDIR/input/B-field_v2.dat)
void runMC(TString inFile = "auau.04gev.0_3fm.10k.f14.gz", TString outFile = "evetest.root", Int_t nStartEvent = 0, Int_t nEvents = 10,
        Bool_t flag_store_FairRadLenPoint = kFALSE, Int_t FieldSwitcher = 0)
{
    TStopwatch timer;
    timer.Start();
    gDebug = 0;

    FairRunSim* fRun = new FairRunSim();
    // Choose the Geant Navigation System
#ifdef GEANT3
    fRun->SetName("TGeant3");
#else
    fRun->SetName("TGeant4");
#endif

    geometry_stage1(fRun); // load mpd geometry
    //geometry_v2(fRun); // load mpd geometry

    // Use extended MC Event header instead of the default one.
    //MpdMCEventHeader* mcHeader = new MpdMCEventHeader();
    //fRun->SetMCEventHeader(mcHeader);

    // Create and Set Event Generator
    FairPrimaryGenerator* primGen = new FairPrimaryGenerator();
    fRun->SetGenerator(primGen);

    // smearing of beam interaction point
    primGen->SetBeam(0.0,0.0,0.1,0.1);
    primGen->SetTarget(0.0,24.0);
    primGen->SmearGausVertexZ(kTRUE);
    primGen->SmearVertexXY(kTRUE);

    // Use user defined decays https://fairroot.gsi.de/?q=node/57
    fRun->SetUserDecay(kTRUE);

#ifdef URQMD // <---- Urqmd  Generator

//V
    gRandom->SetSeed(0);

//Add Rho(770)0 (pid = 113)
    FairBoxGenerator* boxGen1 = new FairBoxGenerator(113, 1);
    boxGen1->SetPtRange(0.0, 3.0); // GeV/c, setPRange vs setPtRange; SetEkinRange    (Double32_t kmin=0   , Double32_t kmax=10)
    boxGen1->SetPhiRange(0, 360); // Azimuth angle range [degree]
    boxGen1->SetEtaRange(-1.0, 1.0); // Polar angle in lab system range [degree]
    boxGen1->SetXYZ(0., 0., 0.); // mm o cm ??
    primGen->AddGenerator(boxGen1);

//Add Phi(1020) (pid = 333)
    FairBoxGenerator* boxGen2 = new FairBoxGenerator(333, 1);
    boxGen2->SetPtRange(0.0, 3.0); // GeV/c, setPRange vs setPtRange; SetEkinRange    (Double32_t kmin=0   , Double32_t kmax=10)
    boxGen2->SetPhiRange(0, 360); // Azimuth angle range [degree]
    boxGen2->SetEtaRange(-1.0, 1.0); // Polar angle in lab system range [degree]
    boxGen2->SetXYZ(0., 0., 0.); // mm o cm ??
    primGen->AddGenerator(boxGen2);

//Add K*(892)0 (pid = 313)
    FairBoxGenerator* boxGen3 = new FairBoxGenerator(313, 1);
    boxGen3->SetPtRange(0.0, 3.0); // GeV/c, setPRange vs setPtRange; SetEkinRange    (Double32_t kmin=0   , Double32_t kmax=10)
    boxGen3->SetPhiRange(0, 360); // Azimuth angle range [degree]
    boxGen3->SetEtaRange(-1.0, 1.0); // Polar angle in lab system range [degree]
    boxGen3->SetXYZ(0., 0., 0.); // mm o cm ??
    primGen->AddGenerator(boxGen3);

//Add anti-K*(892)0 (pid = -313)
    FairBoxGenerator* boxGen4 = new FairBoxGenerator(-313, 1);
    boxGen4->SetPtRange(0.0, 3.0); // GeV/c, setPRange vs setPtRange; SetEkinRange    (Double32_t kmin=0   , Double32_t kmax=10)
    boxGen4->SetPhiRange(0, 360); // Azimuth angle range [degree]
    boxGen4->SetEtaRange(-1.0, 1.0); // Polar angle in lab system range [degree]
    boxGen4->SetXYZ(0., 0., 0.); // mm o cm ??
    primGen->AddGenerator(boxGen4);

//Add K*(892)+ (pid = 323)
    FairBoxGenerator* boxGen5 = new FairBoxGenerator(323, 1);
    boxGen5->SetPtRange(0.0, 3.0); // GeV/c, setPRange vs setPtRange; SetEkinRange    (Double32_t kmin=0   , Double32_t kmax=10)
    boxGen5->SetPhiRange(0, 360); // Azimuth angle range [degree]
    boxGen5->SetEtaRange(-1.0, 1.0); // Polar angle in lab system range [degree]
    boxGen5->SetXYZ(0., 0., 0.); // mm o cm ??
    primGen->AddGenerator(boxGen5);

//Add K*(892)- (pid = -323)
    FairBoxGenerator* boxGen6 = new FairBoxGenerator(-323, 1);
    boxGen6->SetPtRange(0.0, 3.0); // GeV/c, setPRange vs setPtRange; SetEkinRange    (Double32_t kmin=0   , Double32_t kmax=10)
    boxGen6->SetPhiRange(0, 360); // Azimuth angle range [degree]
    boxGen6->SetEtaRange(-1.0, 1.0); // Polar angle in lab system range [degree]
    boxGen6->SetXYZ(0., 0., 0.); // mm o cm ??
    primGen->AddGenerator(boxGen6);

//Add Lambda(1520)- (pid = 3124)
    FairBoxGenerator* boxGen7 = new FairBoxGenerator(3124, 1);
    boxGen7->SetPtRange(0.0, 3.0); // GeV/c, setPRange vs setPtRange; SetEkinRange    (Double32_t kmin=0   , Double32_t kmax=10)
    boxGen7->SetPhiRange(0, 360); // Azimuth angle range [degree]
    boxGen7->SetEtaRange(-1.0, 1.0); // Polar angle in lab system range [degree]
    boxGen7->SetXYZ(0., 0., 0.); // mm o cm ??
    primGen->AddGenerator(boxGen7);

//Add Sigma(1385)+ (pid = 3224)
    FairBoxGenerator* boxGen8 = new FairBoxGenerator(3224, 1);
    boxGen8->SetPtRange(0.0, 3.0); // GeV/c, setPRange vs setPtRange; SetEkinRange    (Double32_t kmin=0   , Double32_t kmax=10)
    boxGen8->SetPhiRange(0, 360); // Azimuth angle range [degree]
    boxGen8->SetEtaRange(-1.0, 1.0); // Polar angle in lab system range [degree]
    boxGen8->SetXYZ(0., 0., 0.); // mm o cm ??
    primGen->AddGenerator(boxGen8);

//Add Sigma(1385)- (pid = 3114)
    FairBoxGenerator* boxGen9 = new FairBoxGenerator(3114, 1);
    boxGen9->SetPtRange(0.0, 3.0); // GeV/c, setPRange vs setPtRange; SetEkinRange    (Double32_t kmin=0   , Double32_t kmax=10)
    boxGen9->SetPhiRange(0, 360); // Azimuth angle range [degree]
    boxGen9->SetEtaRange(-1.0, 1.0); // Polar angle in lab system range [degree]
    boxGen9->SetXYZ(0., 0., 0.); // mm o cm ??
    primGen->AddGenerator(boxGen9);

//Add Xi(1530)0 (pid = 3324)
    FairBoxGenerator* boxGen10 = new FairBoxGenerator(3324, 2);
    boxGen10->SetPtRange(0.0, 3.0); // GeV/c, setPRange vs setPtRange; SetEkinRange    (Double32_t kmin=0   , Double32_t kmax=10)
    boxGen10->SetPhiRange(0, 360); // Azimuth angle range [degree]
    boxGen10->SetEtaRange(-1.0, 1.0); // Polar angle in lab system range [degree]
    boxGen10->SetXYZ(0., 0., 0.); // mm o cm ??
    primGen->AddGenerator(boxGen10);
//V

    if (!CheckFileExist(inFile)) return;

    MpdUrqmdGenerator* urqmdGen = new MpdUrqmdGenerator(inFile);
    // Event plane angle (in degrees) will be generated by uniform distribution from min to max
    Float_t min = 0.0, max = 360.0;
    urqmdGen->SetEventPlane(min * TMath::DegToRad(), max * TMath::DegToRad());
    primGen->AddGenerator(urqmdGen);
    if (nStartEvent > 0) urqmdGen->SkipEvents(nStartEvent);

    // if nEvents is equal 0 then all events (start with nStartEvent) of the given file should be processed
    if (nEvents == 0)
        nEvents = MpdGetNumEvents::GetNumURQMDEvents(inFile.Data()) - nStartEvent;



#else
#ifdef VHLLE
    if (!CheckFileExist(inFile)) return;

    MpdVHLLEGenerator* vhlleGen = new MpdVHLLEGenerator(inFile, kTRUE); // kTRUE corresponds to hydro + cascade, kFALSE -- hydro only
    vhlleGen->SkipEvents(0);
    primGen->AddGenerator(vhlleGen);

#else 
#ifdef FLUID
    if (!CheckFileExist(inFile)) return;

    Mpd3fdGenerator* fluidGen = new Mpd3fdGenerator(inFile);
    fluidGen->SkipEvents(0);
    primGen->AddGenerator(fluidGen);
    //if (nStartEvent > 0) urqmdGen->SkipEvents(nStartEvent);

    // if nEvents is equal 0 then all events (start with nStartEvent) of the given file should be processed
    // if (nEvents == 0)
    //      nEvents = MpdGetNumEvents::GetNumURQMDEvents(inFile.Data()) - nStartEvent;

#else
#ifdef PART // <---- Particle Generator
    FairParticleGenerator* partGen =
            new FairParticleGenerator(211, 10, 1, 0, 3, 1, 0, 0);
    primGen->AddGenerator(partGen);

#else
#ifdef ION // <---- Ion Generator
    FairIonGenerator *fIongen =
            new FairIonGenerator(79, 197, 79, 1, 0., 0., 25, 0., 0., -1.);
    primGen->AddGenerator(fIongen);

#else
#ifdef BOX // <---- Box Generator
    gRandom->SetSeed(0);
    FairBoxGenerator* boxGen = new FairBoxGenerator(22, 1); // 22 = photon; 1 = multipl.
    boxGen->SetPRange(0.5, 0.5); // GeV/c, setPRange vs setPtRange; SetEkinRange    (Double32_t kmin=0   , Double32_t kmax=10)
    boxGen->SetPhiRange(0, 360); // Azimuth angle range [degree]
//    boxGen->SetThetaRange(0, 180); // Polar angle in lab system range [degree]
//    boxGen->SetThetaRange(90, 90); // Polar angle in lab system range [degree]
//    boxGen->SetThetaRange(45, 135); // Polar angle in lab system range [degree]
    boxGen->SetEtaRange(-1.2, 1.2); // Polar angle in lab system range [degree]
    boxGen->SetXYZ(0., 0., 0.); // mm o cm ??
    primGen->AddGenerator(boxGen);

#else
#ifdef HSD // <---- HSD/PHSD Generator
    if (!CheckFileExist(inFile)) return;

    MpdPHSDGenerator *hsdGen = new MpdPHSDGenerator(inFile.Data());
    //hsdGen->SetPsiRP(0.); // set fixed Reaction Plane angle [rad] instead of random
    primGen->AddGenerator(hsdGen);
    if (nStartEvent > 0) hsdGen->SkipEvents(nStartEvent);

    // if nEvents is equal 0 then all events (start with nStartEvent) of the given file should be processed
    if (nEvents == 0)
        nEvents = MpdGetNumEvents::GetNumPHSDEvents(inFile.Data()) - nStartEvent;

#else
#ifdef LAQGSM // <---- LAQGSM Generator
    if (!CheckFileExist(inFile)) return;

    MpdLAQGSMGenerator* guGen = new MpdLAQGSMGenerator(inFile.Data(), kTRUE, 0, 1+nStartEvent+nEvents);
    // kTRUE - for NICA/MPD, 1+nStartEvent+nEvents - search ions in selected part of file.
    primGen->AddGenerator(guGen);
    if (nStartEvent > 0) guGen->SkipEvents(nStartEvent);

    // if nEvents is equal 0 then all events (start with nStartEvent) of the given file should be processed
    if (nEvents == 0)
        nEvents = MpdGetNumEvents::GetNumQGSMEvents(inFile.Data()) - nStartEvent;

#else
#ifdef HADGEN
    THadgen* hadGen = new THadgen();
    hadGen->SetRandomSeed(clock() + time(0));
    hadGen->SetParticleFromPdgCode(0, 196.9665, 79);
    hadGen->SetEnergy(6.5E3);
    MpdGeneralGenerator* generalHad = new MpdGeneralGenerator(hadGen);
    primGen->AddGenerator(generalHad);

#endif
#endif
#endif
#endif
#endif
#endif
#endif
#endif
#endif 

    fRun->SetOutputFile(outFile.Data());

    // Magnetic Field Map - for proper use in the analysis MultiField is necessary here
    MpdMultiField* fField = new MpdMultiField();

    if (FieldSwitcher == 0) {
        MpdConstField* fMagField = new MpdConstField();
        fMagField->SetField(0., 0., 5.); // values are in kG:  1T = 10kG
        fMagField->SetFieldRegion(-230, 230, -230, 230, -375, 375);
        fField->AddField(fMagField);
        fRun->SetField(fField);
        cout << "FIELD at (0., 0., 0.) = (" <<
                fMagField->GetBx(0., 0., 0.) << "; " << fMagField->GetBy(0., 0., 0.) << "; " << fMagField->GetBz(0., 0., 0.) << ")" << endl;
    }
    else if (FieldSwitcher == 1) {
        MpdFieldMap* fMagField = new MpdFieldMap("B-field_v2", "A");
        fMagField->Init();
        fField->AddField(fMagField);
        fRun->SetField(fField);
        cout << "FIELD at (0., 0., 0.) = (" <<
                fMagField->GetBx(0., 0., 0.) << "; " << fMagField->GetBy(0., 0., 0.) << "; " << fMagField->GetBz(0., 0., 0.) << ")" << endl;
    }

    fRun->SetStoreTraj(kTRUE);
    fRun->SetRadLenRegister(flag_store_FairRadLenPoint); // radiation length manager

    //  MpdTpcDigitizerTask* tpcDigitizer = new MpdTpcDigitizerTask();
    //  tpcDigitizer->SetOnlyPrimary(kTRUE); /// Digitize only primary track
    //  tpcDigitizer->SetMakeQA(kTRUE);  /// SetMakeQA(kTRUE) prepares Quality Assurance Histograms  
    //  tpcDigitizer->SetDiffuse(kFALSE);
    //  tpcDigitizer->SetDebug(kFALSE);
    //  tpcDigitizer->SetDistort(kFALSE);
    //  tpcDigitizer->SetResponse(kFALSE);
    //  tpcDigitizer->SetDistribute(kFALSE);
    //  fRun->AddTask(tpcDigitizer);

    fRun->Init();

    // -Trajectories Visualization (TGeoManager Only)
    // Set cuts for storing the trajectories
    FairTrajFilter* trajFilter = FairTrajFilter::Instance();
    trajFilter->SetStepSizeCut(0.01); // 1 cm
    //  trajFilter->SetVertexCut(-2000., -2000., 4., 2000., 2000., 100.);
    trajFilter->SetMomentumCutP(.50); // p_lab > 500 MeV
    //  trajFilter->SetEnergyCut(.2, 3.02); // 0 < Etot < 1.04 GeV

    trajFilter->SetStorePrimaries(kTRUE);
    trajFilter->SetStoreSecondaries(kFALSE);

    // Fill the Parameter containers for this run
    FairRuntimeDb* rtdb = fRun->GetRuntimeDb();

    Bool_t kParameterMerged = kTRUE;
    FairParRootFileIo* output = new FairParRootFileIo(kParameterMerged);
    //AZ output->open(parFile.Data());
    output->open(gFile);
    rtdb->setOutput(output);

    MpdMultiFieldPar* Par = (MpdMultiFieldPar*) rtdb->getContainer("MpdMultiFieldPar");
    if (fField)
        Par->SetParameters(fField);
    Par->setInputVersion(fRun->GetRunId(), 1);
    Par->setChanged();
    // Par->printParams();

    rtdb->saveOutput();
    rtdb->print();

    // Transport nEvents
    fRun->Run(nEvents);

#ifdef LAQGSM
    TString Pdg_table_name = TString::Format("%s%s%c%s", gSystem->BaseName(inFile.Data()), ".g", (fRun->GetName())[6], ".pdg_table.dat");
    (TDatabasePDG::Instance())->WritePDGTable(Pdg_table_name.Data());
#endif

    timer.Stop();
    Double_t rtime = timer.RealTime(), ctime = timer.CpuTime();
    printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime);
    cout << "Macro finished successfully." << endl;     // marker of successful execution for CDASH
}