Fun4All_RunNewTruthMatcher.C

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_RunNewTruthMatcher.C

#ifndef MACRO_FUN4ALLG4SPHENIX_C
#define MACRO_FUN4ALLG4SPHENIX_C

#include <GlobalVariables.C>

#include <DisplayOn.C>
#include <G4Setup_sPHENIX.C>
#include <G4_Mbd.C>
#include <G4_CaloTrigger.C>
#include <G4_Centrality.C>
#include <G4_DSTReader.C>
#include <G4_Global.C>
#include <G4_HIJetReco.C>
#include <G4_Input.C>
#include <G4_Jets.C>
#include <G4_KFParticle.C>
#include <G4_ParticleFlow.C>
#include <G4_Production.C>
#include <G4_TopoClusterReco.C>

#include <Trkr_RecoInit.C>
#include <Trkr_Clustering.C>
#include <Trkr_LaserClustering.C>
#include <Trkr_Reco.C>
#include <Trkr_Eval.C>
#include <Trkr_QA.C>

#include <Trkr_Diagnostics.C>
#include <G4_User.C>
#include <QA.C>

#include <ffamodules/FlagHandler.h>
#include <ffamodules/HeadReco.h>
#include <ffamodules/SyncReco.h>
#include <ffamodules/CDBInterface.h>

#include <fun4all/Fun4AllDstOutputManager.h>
#include <fun4all/Fun4AllOutputManager.h>
#include <fun4all/Fun4AllServer.h>

#include <phool/PHRandomSeed.h>
#include <phool/recoConsts.h>

R__LOAD_LIBRARY(libfun4all.so)
R__LOAD_LIBRARY(libffamodules.so)

// For HepMC Hijing
// try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat

int Fun4All_RunNewTruthMatcher(
    const int nEvents = 1,
    const int skip = 0,
    const string &outputFile = "G4sPHENIX.root",
    const string &inputFile = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
    const string &embed_input_file = "https://www.phenix.bnl.gov/WWW/publish/phnxbld/sPHENIX/files/sPHENIX_G4Hits_sHijing_9-11fm_00000_00010.root",
    const string &outdir = ".")
{
  Fun4AllServer *se = Fun4AllServer::instance();
  se->Verbosity(0);

  //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
  PHRandomSeed::Verbosity(1);

  // just if we set some flags somewhere in this macro
  recoConsts *rc = recoConsts::instance();
  // By default every random number generator uses
  // PHRandomSeed() which reads /dev/urandom to get its seed
  // if the RANDOMSEED flag is set its value is taken as seed
  // You can either set this to a random value using PHRandomSeed()
  // which will make all seeds identical (not sure what the point of
  // this would be:
  //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
  // or set it to a fixed value so you can debug your code
  //  rc->set_IntFlag("RANDOMSEED", 12345);


  //===============
  // Input options
  //===============
  // verbosity setting (applies to all input managers)
  Input::VERBOSITY = 0;
  // First enable the input generators
  // Either:
  // read previously generated g4-hits files, in this case it opens a DST and skips
  // the simulations step completely. The G4Setup macro is only loaded to get information
  // about the number of layers used for the cell reco code
  //  Input::READHITS = true;
  INPUTREADHITS::filename[0] = inputFile;
  // if you use a filelist
  // INPUTREADHITS::listfile[0] = inputFile;
  // Or:
  // Use particle generator
  // And
  // Further choose to embed newly simulated events to a previous simulation. Not compatible with `readhits = true`
  // In case embedding into a production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
  // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
  //  Input::EMBED = true;
  INPUTEMBED::filename[0] = embed_input_file;
  // if you use a filelist
  //INPUTEMBED::listfile[0] = embed_input_file;

  Input::SIMPLE = true;
  // Input::SIMPLE_NUMBER = 2; // if you need 2 of them
  // Input::SIMPLE_VERBOSITY = 1;

  // Enable this is emulating the nominal pp/pA/AA collision vertex distribution
  // Input::BEAM_CONFIGURATION = Input::AA_COLLISION; // Input::AA_COLLISION (default), Input::pA_COLLISION, Input::pp_COLLISION

  //  Input::PYTHIA6 = true;

  // Input::PYTHIA8 = true;

  //  Input::GUN = true;
  //  Input::GUN_NUMBER = 3; // if you need 3 of them
  // Input::GUN_VERBOSITY = 1;

  // Input::COSMIC = true;

  //D0 generator
  //Input::DZERO = false;
  //Input::DZERO_VERBOSITY = 0;
  //Lambda_c generator //Not ready yet
  //Input::LAMBDAC = false;
  //Input::LAMBDAC_VERBOSITY = 0;
  // Upsilon generator
  //Input::UPSILON = true;
  //Input::UPSILON_NUMBER = 3; // if you need 3 of them
  //Input::UPSILON_VERBOSITY = 0;

  //  Input::HEPMC = true;
  INPUTHEPMC::filename = inputFile;
  //-----------------
  // Hijing options (symmetrize hijing, add flow, add fermi motion)
  //-----------------
  //  INPUTHEPMC::HIJINGFLIP = true;
  //  INPUTHEPMC::FLOW = true;
  //  INPUTHEPMC::FLOW_VERBOSITY = 3;
  //  INPUTHEPMC::FERMIMOTION = true;


  // Event pile up simulation with collision rate in Hz MB collisions.
  //Input::PILEUPRATE = 50e3; // 50 kHz for AuAu
  //Input::PILEUPRATE = 3e6; // 3MHz for pp

  // Enable this is emulating the nominal pp/pA/AA collision vertex distribution
  // for HepMC records (hijing, pythia8)
  //  Input::BEAM_CONFIGURATION = Input::AA_COLLISION; // for 2023 sims we want the AA geometry for no pileup sims
  //  Input::BEAM_CONFIGURATION = Input::pp_COLLISION; // for 2024 sims we want the pp geometry for no pileup sims
  //  Input::BEAM_CONFIGURATION = Input::pA_COLLISION; // for pAu sims we want the pA geometry for no pileup sims

  //-----------------
  // Initialize the selected Input/Event generation
  //-----------------
  // This creates the input generator(s)
  InputInit();

  //--------------
  // Set generator specific options
  //--------------
  // can only be set after InputInit() is called

  // Simple Input generator:
  // if you run more than one of these Input::SIMPLE_NUMBER > 1
  // add the settings for other with [1], next with [2]...

  // TEST [Derek, 09.05.2023]
  const UInt_t   nPerEvt(1);
  const TString  sParticle("pi-");
  const Double_t ptParMin(10.);
  const Double_t ptParMax(10.);

  // TEST [Derek, 02.06.2024]
  const Bool_t   doGausVertex(false);
  const Double_t meanVtxXG(0.);
  const Double_t meanVtxYG(0.);
  const Double_t meanVtxZG(0.);
  const Double_t meanVtxXU(0.);
  const Double_t meanVtxYU(0.);
  const Double_t meanVtxZU(0.);
  const Double_t widthVtxXG(0.01);
  const Double_t widthVtxYG(0.01);
  const Double_t widthVtxZG(5.);
  const Double_t widthVtxXU(0.);
  const Double_t widthVtxYU(0.);
  const Double_t widthVtxZU(0.);

  if (Input::SIMPLE)
  {
    INPUTGENERATOR::SimpleEventGenerator[0]->add_particles(sParticle.Data(), nPerEvt);
    if (Input::HEPMC || Input::EMBED)
    {
      INPUTGENERATOR::SimpleEventGenerator[0]->set_reuse_existing_vertex(true);
      INPUTGENERATOR::SimpleEventGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
    }
    else
    {
      if (doGausVertex) {
        INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_function(
          PHG4SimpleEventGenerator::Gaus,
          PHG4SimpleEventGenerator::Gaus,
          PHG4SimpleEventGenerator::Gaus
        );
        INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_mean(meanVtxXG, meanVtxYG, meanVtxZG);
        INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_width(widthVtxXG, widthVtxYG, widthVtxZG);
      } else {
        INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_function(
          PHG4SimpleEventGenerator::Uniform,
          PHG4SimpleEventGenerator::Uniform,
          PHG4SimpleEventGenerator::Uniform
        );
        INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_mean(meanVtxXU, meanVtxYU, meanVtxZU);
        INPUTGENERATOR::SimpleEventGenerator[0] -> set_vertex_distribution_width(widthVtxXU, widthVtxYU, widthVtxZU);
      }
    }
    INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-1, 1);
    INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI, M_PI);
    INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(ptParMin, ptParMax);
  }
  // Upsilons
  // if you run more than one of these Input::UPSILON_NUMBER > 1
  // add the settings for other with [1], next with [2]...
  if (Input::UPSILON)
  {
    INPUTGENERATOR::VectorMesonGenerator[0]->add_decay_particles("e", 0);
    INPUTGENERATOR::VectorMesonGenerator[0]->set_rapidity_range(-1, 1);
    INPUTGENERATOR::VectorMesonGenerator[0]->set_pt_range(0., 10.);
    // Y species - select only one, last one wins
    INPUTGENERATOR::VectorMesonGenerator[0]->set_upsilon_1s();
    if (Input::HEPMC || Input::EMBED)
    {
      INPUTGENERATOR::VectorMesonGenerator[0]->set_reuse_existing_vertex(true);
      INPUTGENERATOR::VectorMesonGenerator[0]->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
    }
  }
  // particle gun
  // if you run more than one of these Input::GUN_NUMBER > 1
  // add the settings for other with [1], next with [2]...
  if (Input::GUN)
  {
    INPUTGENERATOR::Gun[0]->AddParticle("pi-", 0, 1, 0);
    INPUTGENERATOR::Gun[0]->set_vtx(0, 0, 0);
  }

  // pythia6
  if (Input::PYTHIA6)
  {
    Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia6);
  }
  // pythia8
  if (Input::PYTHIA8)
  {
    Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8);
  }

  //--------------
  // Set Input Manager specific options
  //--------------
  // can only be set after InputInit() is called

  if (Input::HEPMC)
  {
    Input::ApplysPHENIXBeamParameter(INPUTMANAGER::HepMCInputManager);

    // optional overriding beam parameters
    //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_width(100e-4, 100e-4, 8, 0);  //optional collision smear in space, time
    //    INPUTMANAGER::HepMCInputManager->set_vertex_distribution_mean(0,0,0,0);//optional collision central position shift in space, time
    // //optional choice of vertex distribution function in space, time
    //INPUTMANAGER::HepMCInputManager->set_vertex_distribution_function(PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus, PHHepMCGenHelper::Gaus);
    //INPUTMANAGER::HepMCInputManager->set_embedding_id(Input::EmbedID);
    if (Input::PILEUPRATE > 0)
    {
      // Copy vertex settings from foreground hepmc input
      INPUTMANAGER::HepMCPileupInputManager->CopyHelperSettings(INPUTMANAGER::HepMCInputManager);
      // and then modify the ones you want to be different
      // INPUTMANAGER::HepMCPileupInputManager->set_vertex_distribution_width(100e-4,100e-4,8,0);
    }
  }
  if (Input::PILEUPRATE > 0)
  {
    Input::ApplysPHENIXBeamParameter(INPUTMANAGER::HepMCPileupInputManager);
  }
  // register all input generators with Fun4All
  InputRegister();

  if (! Input::READHITS)
  {
    rc->set_IntFlag("RUNNUMBER",1);

    SyncReco *sync = new SyncReco();
    se->registerSubsystem(sync);

    HeadReco *head = new HeadReco();
    se->registerSubsystem(head);
  }
// Flag Handler is always needed to read flags from input (if used)
// and update our rc flags with them. At the end it saves all flags
// again on the DST in the Flags node under the RUN node
  FlagHandler *flag = new FlagHandler();
  se->registerSubsystem(flag);

  // set up production relatedstuff
  //   Enable::PRODUCTION = true;

  //======================
  // Write the DST
  //======================

  //Enable::DSTOUT = true;
  Enable::DSTOUT_COMPRESS = false;
  DstOut::OutputDir = outdir;
  DstOut::OutputFile = outputFile;

  //Option to convert DST to human command readable TTree for quick poke around the outputs
  //  Enable::DSTREADER = true;

  // turn the display on (default off)
   //Enable::DISPLAY = true;

  //======================
  // What to run
  //======================

  // QA, main switch
  Enable::QA = false;

  // Global options (enabled for all enables subsystems - if implemented)
  //  Enable::ABSORBER = true;
  //  Enable::OVERLAPCHECK = true;
  //  Enable::VERBOSITY = 1;

  // Enable::MBD = true;
  // Enable::MBD_SUPPORT = true; // save hist in MBD/BBC support structure
  // Enable::MBDRECO = Enable::MBD && true;
  Enable::MBDFAKE = true;  // Smeared vtx and t0, use if you don't want real MBD/BBC in simulation

  Enable::PIPE = true;
  Enable::PIPE_ABSORBER = true;

  // central tracking
  Enable::MVTX = true;
  Enable::MVTX_CELL = Enable::MVTX && true;
  Enable::MVTX_CLUSTER = Enable::MVTX_CELL && true;
  Enable::MVTX_QA = Enable::MVTX_CLUSTER && Enable::QA && true;

  Enable::INTT = true;
//  Enable::INTT_ABSORBER = true; // enables layerwise support structure readout
//  Enable::INTT_SUPPORT = true; // enable global support structure readout
  Enable::INTT_CELL = Enable::INTT && true;
  Enable::INTT_CLUSTER = Enable::INTT_CELL && true;
  Enable::INTT_QA = Enable::INTT_CLUSTER && Enable::QA && true;

  Enable::TPC = true;
  Enable::TPC_ABSORBER = true;
  Enable::TPC_CELL = Enable::TPC && true;
  Enable::TPC_CLUSTER = Enable::TPC_CELL && true;
  Enable::TPC_QA = Enable::TPC_CLUSTER && Enable::QA && true;

  Enable::MICROMEGAS = true;
  Enable::MICROMEGAS_CELL = Enable::MICROMEGAS && true;
  Enable::MICROMEGAS_CLUSTER = Enable::MICROMEGAS_CELL && true;
  Enable::MICROMEGAS_QA = Enable::MICROMEGAS_CLUSTER && Enable::QA && true;

  Enable::TRACKING_TRACK = (Enable::MICROMEGAS_CLUSTER && Enable::TPC_CLUSTER && Enable::INTT_CLUSTER && Enable::MVTX_CLUSTER) && true;
  Enable::TRACKING_EVAL = Enable::TRACKING_TRACK && Enable::GLOBAL_RECO && true;
  Enable::TRACKING_QA = Enable::TRACKING_TRACK && Enable::QA && true;

  // only do track matching if TRACKINGTRACK is also used
  Enable::TRACK_MATCHING = Enable::TRACKING_TRACK && true;
  Enable::TRACK_MATCHING_TREE = Enable::TRACK_MATCHING && true;
  Enable::TRACK_MATCHING_TREE_CLUSTERS = Enable::TRACK_MATCHING_TREE && true;

  //Additional tracking tools 
  //Enable::TRACKING_DIAGNOSTICS = Enable::TRACKING_TRACK && true;
  //G4TRACKING::filter_conversion_electrons = true;


  //  cemc electronics + thin layer of W-epoxy to get albedo from cemc
  //  into the tracking, cannot run together with CEMC
  //  Enable::CEMCALBEDO = true;

  Enable::CEMC = true;
  Enable::CEMC_ABSORBER = true;
  Enable::CEMC_CELL = Enable::CEMC && true;
  Enable::CEMC_TOWER = Enable::CEMC_CELL && true;
  Enable::CEMC_CLUSTER = Enable::CEMC_TOWER && true;
  Enable::CEMC_EVAL = Enable::CEMC_G4Hit && Enable::CEMC_CLUSTER && false;
  Enable::CEMC_QA = Enable::CEMC_CLUSTER && Enable::QA && false;

  Enable::HCALIN = true;
  Enable::HCALIN_ABSORBER = true;
  Enable::HCALIN_CELL = Enable::HCALIN && true;
  Enable::HCALIN_TOWER = Enable::HCALIN_CELL && true;
  Enable::HCALIN_CLUSTER = Enable::HCALIN_TOWER && true;
  Enable::HCALIN_EVAL = Enable::HCALIN_G4Hit && Enable::HCALIN_CLUSTER && false;
  Enable::HCALIN_QA = Enable::HCALIN_CLUSTER && Enable::QA && false;

  Enable::MAGNET = true;
  Enable::MAGNET_ABSORBER = true;

  Enable::HCALOUT = true;
  Enable::HCALOUT_ABSORBER = true;
  Enable::HCALOUT_CELL = Enable::HCALOUT && true;
  Enable::HCALOUT_TOWER = Enable::HCALOUT_CELL && true;
  Enable::HCALOUT_CLUSTER = Enable::HCALOUT_TOWER && true;
  Enable::HCALOUT_EVAL = Enable::HCALOUT_G4Hit && Enable::HCALOUT_CLUSTER && false;
  Enable::HCALOUT_QA = Enable::HCALOUT_CLUSTER && Enable::QA && false;

  Enable::EPD = true;
  Enable::EPD_TILE = Enable::EPD && true;

  Enable::BEAMLINE = true;
//  Enable::BEAMLINE_ABSORBER = true;  // makes the beam line magnets sensitive volumes
//  Enable::BEAMLINE_BLACKHOLE = true; // turns the beamline magnets into black holes
  Enable::ZDC = true;
//  Enable::ZDC_ABSORBER = true;
//  Enable::ZDC_SUPPORT = true;
  Enable::ZDC_TOWER = Enable::ZDC && true;
  Enable::ZDC_EVAL = Enable::ZDC_TOWER && true;

  //Enable::PLUGDOOR = true;
  Enable::PLUGDOOR_ABSORBER = true;

  Enable::GLOBAL_RECO = (Enable::MBDFAKE || Enable::TRACKING_TRACK) && true;
  //Enable::GLOBAL_FASTSIM = true;

  //Enable::KFPARTICLE = true;
  //Enable::KFPARTICLE_VERBOSITY = 1;
  //Enable::KFPARTICLE_TRUTH_MATCH = true;
  //Enable::KFPARTICLE_SAVE_NTUPLE = true;

  Enable::CALOTRIGGER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;

  Enable::JETS = (Enable::GLOBAL_RECO || Enable::GLOBAL_FASTSIM) && false;
  Enable::JETS_EVAL = Enable::JETS && false;
  Enable::JETS_QA = Enable::JETS && Enable::QA && false;

  // HI Jet Reco for p+Au / Au+Au collisions (default is false for
  // single particle / p+p-only simulations, or for p+Au / Au+Au
  // simulations which don't particularly care about jets)
  Enable::HIJETS = Enable::JETS && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;

  // 3-D topoCluster reconstruction, potentially in all calorimeter layers
  Enable::TOPOCLUSTER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
  // particle flow jet reconstruction - needs topoClusters!
  Enable::PARTICLEFLOW = Enable::TOPOCLUSTER && false;
  // centrality reconstruction
  Enable::CENTRALITY = true;

  // new settings using Enable namespace in GlobalVariables.C
  Enable::BLACKHOLE = true;
  //Enable::BLACKHOLE_SAVEHITS = false; // turn off saving of bh hits
  //Enable::BLACKHOLE_FORWARD_SAVEHITS = false; // disable forward/backward hits
  //BlackHoleGeometry::visible = true;

  // run user provided code (from local G4_User.C)
  //Enable::USER = true;

  //===============
  // conditions DB flags
  //===============
  Enable::CDB = true;
  // global tag
  rc->set_StringFlag("CDB_GLOBALTAG",CDB::global_tag);
  // 64 bit timestamp
  rc->set_uint64Flag("TIMESTAMP",CDB::timestamp);
  //---------------
  // World Settings
  //---------------
  //  G4WORLD::PhysicsList = "FTFP_BERT"; //FTFP_BERT_HP best for calo
  //  G4WORLD::WorldMaterial = "G4_AIR"; // set to G4_GALACTIC for material scans

  //---------------
  // Magnet Settings
  //---------------

  //  G4MAGNET::magfield =  string(getenv("CALIBRATIONROOT"))+ string("/Field/Map/sphenix3dbigmapxyz.root");  // default map from the calibration database
  //  G4MAGNET::magfield = "1.5"; // alternatively to specify a constant magnetic field, give a float number, which will be translated to solenoidal field in T, if string use as fieldmap name (including path)
//  G4MAGNET::magfield_rescale = 1.;  // make consistent with expected Babar field strength of 1.4T

  //---------------
  // Pythia Decayer
  //---------------
  // list of decay types in
  // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
  // default is All:
  // G4P6DECAYER::decayType = EDecayType::kAll;

  // Initialize the selected subsystems
  G4Init();

  //---------------------
  // GEANT4 Detector description
  //---------------------
  if (!Input::READHITS)
  {
    G4Setup();
  }

  //------------------
  // Detector Division
  //------------------

  if ((Enable::MBD && Enable::MBDRECO) || Enable::MBDFAKE) Mbd_Reco();

  if (Enable::MVTX_CELL) Mvtx_Cells();
  if (Enable::INTT_CELL) Intt_Cells();
  if (Enable::TPC_CELL) TPC_Cells();
  if (Enable::MICROMEGAS_CELL) Micromegas_Cells();

  if (Enable::CEMC_CELL) CEMC_Cells();

  if (Enable::HCALIN_CELL) HCALInner_Cells();

  if (Enable::HCALOUT_CELL) HCALOuter_Cells();

  //-----------------------------
  // CEMC towering and clustering
  //-----------------------------

  if (Enable::CEMC_TOWER) CEMC_Towers();
  if (Enable::CEMC_CLUSTER) CEMC_Clusters();

  //--------------
  // EPD tile reconstruction
  //--------------

  if (Enable::EPD_TILE) EPD_Tiles();

  //-----------------------------
  // HCAL towering and clustering
  //-----------------------------

  if (Enable::HCALIN_TOWER) HCALInner_Towers();
  if (Enable::HCALIN_CLUSTER) HCALInner_Clusters();

  if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
  if (Enable::HCALOUT_CLUSTER) HCALOuter_Clusters();

  // if enabled, do topoClustering early, upstream of any possible jet reconstruction
  if (Enable::TOPOCLUSTER) TopoClusterReco();

  //--------------
  // SVTX tracking
  //--------------
  if(Enable::TRACKING_TRACK)
    {
      TrackingInit();
    }
  if (Enable::MVTX_CLUSTER) Mvtx_Clustering();
  if (Enable::INTT_CLUSTER) Intt_Clustering();
  if (Enable::TPC_CLUSTER)
    {
      if(G4TPC::ENABLE_DIRECT_LASER_HITS || G4TPC::ENABLE_CENTRAL_MEMBRANE_HITS)
        {
          TPC_LaserClustering();
        }
      else
        {
          TPC_Clustering();
        }
    }
  if (Enable::MICROMEGAS_CLUSTER) Micromegas_Clustering();

  if (Enable::TRACKING_TRACK)
  {
    Tracking_Reco();
  }

  

  if(Enable::TRACKING_DIAGNOSTICS)
    {
      const std::string kshortFile = "./kshort_" + outputFile;
      const std::string residualsFile = "./residuals_" + outputFile;
 
      G4KshortReconstruction(kshortFile);
      seedResiduals(residualsFile);
    }


  //-----------------
  // Global Vertexing
  //-----------------

  if (Enable::GLOBAL_RECO && Enable::GLOBAL_FASTSIM)
  {
    cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
    gSystem->Exit(1);
  }
  if (Enable::GLOBAL_RECO)
  {
    Global_Reco();
  }
  else if (Enable::GLOBAL_FASTSIM)
  {
    Global_FastSim();
  }

  //-----------------
  // Centrality Determination
  //-----------------

  if (Enable::CENTRALITY)
  {
      Centrality();
  }

  //-----------------
  // Calo Trigger Simulation
  //-----------------

  if (Enable::CALOTRIGGER)
  {
    CaloTrigger_Sim();
  }

  //---------
  // Jet reco
  //---------

  if (Enable::JETS) Jet_Reco();
  if (Enable::HIJETS) HIJetReco();

  if (Enable::PARTICLEFLOW) ParticleFlow();

  //----------------------
  // Simulation evaluation
  //----------------------
  string outputroot = outputFile;
  string remove_this = ".root";
  size_t pos = outputroot.find(remove_this);
  if (pos != string::npos)
  {
    outputroot.erase(pos, remove_this.length());
  }

  //if (Enable::TRACKING_EVAL) {
  if (true) {

    build_truthreco_tables(); 

    std::cout << "CHECK I am inside the svtx evaluator block!" << std::endl;

    SvtxEvaluator* eval = new SvtxEvaluator(
      "SVTXEVALUATOR",
      outputroot + "_g4svtx_eval.root",
      "SvtxTrackMap",
      G4MVTX::n_maps_layer,
      G4INTT::n_intt_layer,
      G4TPC::n_gas_layer,
      G4MICROMEGAS::n_micromegas_layer
    );
    eval->do_cluster_eval(true);
    eval->do_g4hit_eval(false);
    eval->do_hit_eval(false);  // enable to see the hits that includes the chamber physics...
    eval->do_gpoint_eval(true);
    eval->do_vtx_eval_light(true);
    eval->do_eval_light(true);
    eval->do_track_eval(true);
    eval->do_gtrack_eval(true);
    eval->do_track_match(true);
    eval->set_use_initial_vertex(G4TRACKING::g4eval_use_initial_vertex);

    bool embed_scan = true;
    if(TRACKING::pp_mode) {
      embed_scan = false;
    }
    eval->scan_for_embedded(embed_scan);   // take all tracks if false - take only embedded tracks if true
    eval->scan_for_primaries(embed_scan);  // defaults to only thrown particles for ntp_gtrack
    std::cout << "SvtxEvaluator: pp_mode set to " << TRACKING::pp_mode << " and scan_for_embedded set to " << embed_scan << std::endl;
    eval->Verbosity(Enable::TRACKING_VERBOSITY);
    se->registerSubsystem(eval);
  }  // end if (Enable::TRACKING_EVAL)

  if (Enable::CEMC_EVAL) CEMC_Eval(outputroot + "_g4cemc_eval.root");

  if (Enable::HCALIN_EVAL) HCALInner_Eval(outputroot + "_g4hcalin_eval.root");

  if (Enable::HCALOUT_EVAL) HCALOuter_Eval(outputroot + "_g4hcalout_eval.root");

  if (Enable::JETS_EVAL) Jet_Eval(outputroot + "_g4jet_eval.root");

  if (Enable::DSTREADER) G4DSTreader(outputroot + "_DSTReader.root");



  if (Enable::USER) UserAnalysisInit();

  // Writes electrons from conversions to a new track map on the node tree
  // the ntuple file is for diagnostics, it is produced only if the flag is set in G4_Tracking.C
  if(G4TRACKING::filter_conversion_electrons) Filter_Conversion_Electrons(outputroot + "_secvert_ntuple.root");


  //======================
  // Run KFParticle on evt
  //======================
  if (Enable::KFPARTICLE && Input::UPSILON) KFParticle_Upsilon_Reco();
  if (Enable::KFPARTICLE && Input::DZERO) KFParticle_D0_Reco();

  //----------------------
  // Standard QAs
  //----------------------

  if (Enable::CEMC_QA) CEMC_QA();
  if (Enable::HCALIN_QA) HCALInner_QA();
  if (Enable::HCALOUT_QA) HCALOuter_QA();

  if (Enable::JETS_QA) Jet_QA();

  if (Enable::MVTX_QA) Mvtx_QA();
  if (Enable::INTT_QA) Intt_QA();
  if (Enable::TPC_QA) TPC_QA();
  if (Enable::MICROMEGAS_QA) Micromegas_QA();
  if (Enable::TRACKING_QA) Tracking_QA();

  if (Enable::TRACKING_QA && Enable::CEMC_QA && Enable::HCALIN_QA && Enable::HCALOUT_QA) QA_G4CaloTracking();

  //if (Enable::TRACK_MATCHING) {
  if (true) {

    // initialize track matcher
    TrkrClusterIsMatcher* ismatcher = new TrkrClusterIsMatcher();

    // These are the default values -- uncomment and change as desired
    //ismatcher->single_pixel_phi_MVTX = false ; // default to pitch*max(N_pixels_M,N_pixels_T)*tol_MVTX
    //ismatcher->single_pixel_phi_INTT = false ; // ... same as for MVTX
    //ismatcher->single_bin_phi_TPC    = true  ;  // default to pitch*tol_phi_TPC
    //ismatcher->single_pixel_z_MVTX = false ; // default to pitch*max(N_pixels_M,N_pixels_T)*tol_z_MVTX
    //ismatcher->single_pixel_z_INTT = false ; // ... same as for MVTX
    //ismatcher->single_bin_t_TPC    = true  ;  // default to pitch*tol_t_TPC
    //ismatcher-> tol_phi_MVTX = 0.5;
    //ismatcher-> tol_phi_INTT = 0.5;
    //ismatcher-> tol_phi_TPC  = 1.0;
    //ismatcher-> tol_z_MVTX  = 0.5;
    //ismatcher-> tol_t_TPC   = 1.0;

    auto trackmatcher = new TruthRecoTrackMatching(ismatcher);
    trackmatcher->set_min_cl_match     (5);    // minimum number of matched clusters to make a matched track
    trackmatcher->set_min_cl_ratio     (0.1);  // at least 10% of truth clusters must be matched
    trackmatcher->set_cutoff_deta      (0.3);  // won't compare tracks with |Δeta|>0.3 away
    trackmatcher->set_cutoff_dphi      (0.3);  // won't compare tracks with |Δphi|>0.3 away
    trackmatcher->set_smallsearch_deta (0.05); // will first compare tracks within this |Δphi|
    trackmatcher->set_smallsearch_dphi (0.05); // will first compare tracks within this |Δeta|
    //trackmatcher->set_max_nreco_per_truth (4); // maximum reco tracks matched for any truth track
    //trackmatcher->set_max_ntruth_per_reco (4); // maximum truth tracks matched for any reco track
    trackmatcher->set_max_nreco_per_truth (1); // maximum reco tracks matched for any truth track // TEST [Derek, 10.13.2023]
    trackmatcher->set_max_ntruth_per_reco (1); // maximum truth tracks matched for any reco track // TEST [Derek, 10.13.2023]

    // set verbosity
    int verbosity = std::max(Enable::VERBOSITY, Enable::TRACKING_VERBOSITY);
    trackmatcher->Verbosity(verbosity);
    se->registerSubsystem(trackmatcher);

    if (Enable::TRACK_MATCHING_TREE) {
      auto treefiller = new FillClusMatchTree(ismatcher, outputroot + "_g4trackmatching.root");//, true, true, true, outputFile);
      treefiller->Verbosity(verbosity);
      if (Enable::TRACK_MATCHING_TREE_CLUSTERS) {
        treefiller->m_fill_clusters     = true;
        treefiller->m_fill_SvUnmatched  = true;
      } else {
        treefiller->m_fill_clusters     = false;
        treefiller->m_fill_SvUnmatched  = false;
      }
      treefiller->m_fill_clusverbose = false;
      se->registerSubsystem(treefiller);
    }
  }  // end if (ENABLE::TRACK_MATCHING)

  //--------------
  // Set up Input Managers
  //--------------

  InputManagers();

  if (Enable::PRODUCTION)
  {
    Production_CreateOutputDir();
  }

  if (Enable::DSTOUT)
  {
    string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
    Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
    if (Enable::DSTOUT_COMPRESS)
    {
      ShowerCompress();
      DstCompress(out);
    }
    se->registerOutputManager(out);
  }
  //-----------------
  // Event processing
  //-----------------
  if (Enable::DISPLAY)
  {
    DisplayOn();

    gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
    gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");

    cout << "-------------------------------------------------" << endl;
    cout << "You are in event display mode. Run one event with" << endl;
    cout << "se->run(1)" << endl;
    cout << "Run Geant4 command with following examples" << endl;
    gROOT->ProcessLine("displaycmd()");

    return 0;
  }

  // if we use a negative number of events we go back to the command line here
  if (nEvents < 0)
  {
    return 0;
  }
  // if we run the particle generator and use 0 it'll run forever
  // for embedding it runs forever if the repeat flag is set
  if (nEvents == 0 && !Input::HEPMC && !Input::READHITS && INPUTEMBED::REPEAT)
  {
    cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
    cout << "it will run forever, so I just return without running anything" << endl;
    return 0;
  }

  se->skip(skip);
  se->run(nEvents);

  //-----
  // QA output
  //-----

  if (Enable::QA) QA_Output(outputroot + "_qa.root");

  //-----
  // Exit
  //-----

//  CDBInterface::instance()->Print(); // print used DB files
  se->End();
  std::cout << "All done" << std::endl;
  delete se;
  if (Enable::PRODUCTION)
  {
    Production_MoveOutput();
  }

  gSystem->Exit(0);
  return 0;
}
#endif