Fun4All_G4_EICDetector_LQ_reference.C

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

int Fun4All_G4_EICDetector_LQ_reference(
                              string n="1093",
                              string ebeam="20",
                              string pbeam="250",
                              //string inputFile,
                              string inputFile="/direct/phenix+u/spjeffas/LQGENEP/TestOut.1093event.root",
                              string output="",
                              const char * outputFile = "G4EICDetector.root"
                           )
{
  // Set the number of TPC layer
  const int n_TPC_layers = 40;  // use 60 for backward compatibility only
  
  //Get parameter variables from parameter file
  
  int nEvents;
  stringstream geek(n);
  geek>>nEvents;



  string directory = "/direct/phenix+u/spjeffas/leptoquark/output/"+output+"/";


  //===============
  // Input options
  //===============
  // 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
  //
  // In case reading 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/
  const bool readhits = false;
  // Or:
  // read files in HepMC format (typically output from event generators like hijing or pythia)
  const bool readhepmc = false; // read HepMC files
  // Or:
  // read files in EICTree format generated by eicsmear package
  const bool readeictree = true;
  // Or:
  // Use Pythia 8
  const bool runpythia8 = false;
  // Or:
  // Use Pythia 6
  const bool runpythia6 = false;
  // Or:
  // Use HEPGen
  const bool runhepgen = false;
  // Or:
  // Use Sartre
  const bool runsartre = false;



  // Besides the above flags. One can further choose to further put in following particles in Geant4 simulation
  // Use multi-particle generator (PHG4SimpleEventGenerator), see the code block below to choose particle species and kinematics
  const bool particles = false && !readhits;
  // or gun/ very simple single particle gun generator
  const bool usegun = false && !readhits;
  // Throw single Upsilons, may be embedded in Hijing by setting readhepmc flag also  (note, careful to set Z vertex equal to Hijing events)
  const bool upsilons = false && !readhits;

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

  // sPHENIX barrel
  bool do_bbc = true;

  bool do_pipe = true;

  bool do_svtx = true;
  bool do_svtx_cell = do_svtx && true;
  bool do_svtx_track = do_svtx_cell && true;
  bool do_svtx_eval = do_svtx_track && true;

  bool do_pstof = false;

  bool do_cemc = true;

  bool do_cemc_cell = true;
  bool do_cemc_twr = true;
  bool do_cemc_cluster = true;
  bool do_cemc_eval = true;

  bool do_hcalin = true;
  bool do_hcalin_cell = true;
  bool do_hcalin_twr = true;
  bool do_hcalin_cluster = true;
  bool do_hcalin_eval = true;

  bool do_cemc_cell = do_cemc && true;
  bool do_cemc_twr = do_cemc_cell && true;
  bool do_cemc_cluster = do_cemc_twr && true;
  bool do_cemc_eval = do_cemc_cluster && true;

  bool do_hcalin = true;
  bool do_hcalin_cell = do_hcalin && true;
  bool do_hcalin_twr = do_hcalin_cell && true;
  bool do_hcalin_cluster = do_hcalin_twr && true;
  bool do_hcalin_eval = do_hcalin_cluster && true;


  bool do_magnet = true;

  bool do_hcalout = true;

  bool do_hcalout_cell = true;
  bool do_hcalout_twr = true;
  bool do_hcalout_cluster = true;
  bool do_hcalout_eval = true;

  bool do_global = true;
  bool do_global_fastsim = false;

  bool do_jet_reco = true;
  bool do_jet_eval = true;

  bool do_fwd_jet_reco = true;
  bool do_fwd_jet_eval = false;

  bool do_hcalout_cell = do_hcalout && true;
  bool do_hcalout_twr = do_hcalout_cell && true;
  bool do_hcalout_cluster = do_hcalout_twr && true;
  bool do_hcalout_eval = do_hcalout_cluster && true;


  // EICDetector geometry - barrel
  bool do_DIRC = true;

  // EICDetector geometry - 'hadron' direction
  bool do_FGEM = true;
  bool do_FGEM_track = do_FGEM &&  false;

  bool do_RICH = true;
  bool do_Aerogel = true;

  bool do_FEMC = true;
  bool do_FEMC_cell = do_FEMC && true;
  bool do_FEMC_twr = do_FEMC_cell && true;
  bool do_FEMC_cluster = do_FEMC_twr && true;
  bool do_FEMC_eval = do_FEMC_cluster && true;

  bool do_FHCAL = true;
  bool do_FHCAL_cell = do_FHCAL && true;
  bool do_FHCAL_twr = do_FHCAL_cell && true;
  bool do_FHCAL_cluster = do_FHCAL_twr && true;
  bool do_FHCAL_eval = do_FHCAL_cluster && true;

  // EICDetector geometry - 'electron' direction
  bool do_EGEM = true;
  bool do_EGEM_track = do_EGEM &&  false;

  bool do_EEMC = true;
  bool do_EEMC_cell = do_EEMC && true;
  bool do_EEMC_twr = do_EEMC_cell && true;
  bool do_EEMC_cluster = do_EEMC_twr && true;
  bool do_EEMC_eval = do_EEMC_cluster && true;

  //do leptoquark analysis modules
  bool do_lepto_analysis = true;

  // Other options
  bool do_global = true;
  bool do_global_fastsim = false;

  bool do_calotrigger = false && do_cemc_twr && do_hcalin_twr && do_hcalout_twr;

  bool do_jet_reco = true;
  bool do_jet_eval = do_jet_reco && true;

  bool do_fwd_jet_reco = true;
  bool do_fwd_jet_eval = do_fwd_jet_reco && true;

  // HI Jet Reco for jet simulations in Au+Au (default is false for
  // single particle / p+p simulations, or for Au+Au simulations which
  // don't care about jets)
  bool do_HIjetreco = false && do_jet_reco && do_cemc_twr && do_hcalin_twr && do_hcalout_twr;

  // Compress DST files
  bool do_dst_compress = false;

  //Option to convert DST to human command readable TTree for quick poke around the outputs
  bool do_DSTReader = false;

  //---------------
  // Load libraries
  //---------------

  gSystem->Load("libfun4all.so");
  gSystem->Load("libg4detectors.so");
  gSystem->Load("libphhepmc.so");
  gSystem->Load("libg4testbench.so");
  gSystem->Load("libg4hough.so");
  gSystem->Load("libg4calo.so");
  gSystem->Load("libg4eval.so");
  gSystem->Load("libeicana.so");

  // establish the geometry and reconstruction setup
  gROOT->LoadMacro("G4Setup_EICDetector.C");
  G4Init(do_svtx,do_cemc,do_hcalin,do_magnet,do_hcalout,do_pipe,do_FGEM,do_EGEM,do_FEMC,do_FHCAL,do_EEMC,do_DIRC,do_RICH,do_Aerogel,n_TPC_layers);

  int absorberactive = 0; // set to 1 to make all absorbers active volumes
  //  const string magfield = "1.5"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
  const string magfield = "/phenix/upgrades/decadal/fieldmaps/sPHENIX.2d.root"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
  const float magfield_rescale = 1.4/1.5; // scale the map to a 1.4 T field

  //---------------
  // Fun4All server
  //---------------

  Fun4AllServer *se = Fun4AllServer::instance();
  se->Verbosity(0); // uncomment for batch production running with minimal output messages
  // se->Verbosity(Fun4AllServer::VERBOSITY_SOME); // uncomment for some info for interactive running

  // 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);

  //-----------------
  // Event generation
  //-----------------

  if (readhits)
    {
      // Get the hits from a file
      // The input manager is declared later
    }
  else if (readhepmc)
    {
      // this module is needed to read the HepMC records into our G4 sims
      // but only if you read HepMC input files
      HepMCNodeReader *hr = new HepMCNodeReader();
      se->registerSubsystem(hr);
    }
  else if (readeictree)
    {
      // this module is needed to read the EICTree style records into our G4 sims
      ReadEICFiles *eicr = new ReadEICFiles();
      eicr->OpenInputFile(inputFile);

      se->registerSubsystem(eicr);
    }
  else if (runpythia8)
    {
      gSystem->Load("libPHPythia8.so");

      PHPythia8* pythia8 = new PHPythia8();
      // see coresoftware/generators/PHPythia8 for example config
      pythia8->set_config_file("/direct/phenix+u/spjeffas/coresoftware/generators/PHPythia8/phpythia8.cfg");
      se->registerSubsystem(pythia8);

      HepMCNodeReader *hr = new HepMCNodeReader();
      se->registerSubsystem(hr);
    }
  else if (runpythia6)
    {
      gSystem->Load("libPHPythia6.so");

      PHPythia6 *pythia6 = new PHPythia6();
      // see coresoftware/generators/PHPythia6 for example config
      pythia6->set_config_file("/direct/phenix+u/spjeffas/coresoftware/generators/PHPythia6/phpythia6_ep.cfg");
      se->registerSubsystem(pythia6);

      HepMCNodeReader *hr = new HepMCNodeReader();
      se->registerSubsystem(hr);
    }
  else if (runhepgen)
    {
      gSystem->Load("libsHEPGen.so");

      sHEPGen *hepgen = new sHEPGen();
      // see HEPGen source directory/share/vggdata for required .dat files
      // see HEPGen source directory/share/datacards for required datacard files
      hepgen->set_datacard_file("hepgen_dvcs.data");
      hepgen->set_momentum_electron(-20);
      hepgen->set_momentum_hadron(250);
      se->registerSubsystem(hepgen);

      HepMCNodeReader *hr = new HepMCNodeReader();
      se->registerSubsystem(hr);
    }
  else if (runsartre)
    {
      // see coresoftware/generators/PHSartre/README for setup instructions
      // before running:
      // setenv SARTRE_DIR /opt/sphenix/core/sartre-1.20_root-5.34.36
      gSystem->Load("libPHSartre.so");

      PHSartre* mysartre = new PHSartre();
      // see coresoftware/generators/PHSartre for example config
      mysartre->set_config_file("sartre.cfg");

      // particle trigger to enhance forward J/Psi -> ee
      PHSartreParticleTrigger* pTrig = new PHSartreParticleTrigger("MySartreTrigger");
      pTrig->AddParticles(-11);
      //pTrig->SetEtaHighLow(4.0,1.4);
      pTrig->SetEtaHighLow(1.0,-1.1);  // central arm
      pTrig->PrintConfig();
      mysartre->register_trigger((PHSartreGenTrigger *)pTrig);
      se->registerSubsystem(mysartre);

      HepMCNodeReader *hr = new HepMCNodeReader();
      se->registerSubsystem(hr);
    }

  // If "readhepMC" is also set, the particles will be embedded in Hijing events
  if(particles)
    {
      // toss low multiplicity dummy events
      PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();

      //gen->add_particles("e-",5); // mu+,e+,proton,pi+,Upsilon
      //gen->add_particles("e+",5); // mu-,e-,anti_proton,pi-
      gen->add_particles("tau-",1); // mu-,e-,anti_proton,pi-
      if (readhepmc) {
        gen->set_reuse_existing_vertex(true);
        gen->set_existing_vertex_offset_vector(0.0,0.0,0.0);
      } else {
        gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
                                              PHG4SimpleEventGenerator::Uniform,
                                              PHG4SimpleEventGenerator::Uniform);
        gen->set_vertex_distribution_mean(0.0,0.0,0.0);
        gen->set_vertex_distribution_width(0.0,0.0,5.0);
      }
      gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
      gen->set_vertex_size_parameters(0.0,0.0);
      gen->set_eta_range(0.1, 0.1);
      //gen->set_eta_range(3.0, 3.0); //EICDetector FWD
      gen->set_phi_range(TMath::Pi()/2-0.1, TMath::Pi()/2-0.1);
      //gen->set_phi_range(TMath::Pi()/2-0.1, TMath::Pi()/2-0.1);
      gen->set_p_range(30.0, 30.0);

      //gen->add_particles("pi-",1); // mu+,e+,proton,pi+,Upsilon
      //gen->add_particles("pi+",100); // 100 pion option
      if (readhepmc)
        {
          gen->set_reuse_existing_vertex(true);
          gen->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
        }
      else
        {
          gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
                                                PHG4SimpleEventGenerator::Uniform,
                                                PHG4SimpleEventGenerator::Uniform);
          gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
          gen->set_vertex_distribution_width(0.0, 0.0, 0.0);
        }
      gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
      gen->set_vertex_size_parameters(0.0, 0.0);
      gen->set_eta_range(-1.0, 1.0);
      gen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
      //gen->set_pt_range(0.1, 50.0);
      gen->set_pt_range(0.1, 20.0);

      gen->Embed(1);
      gen->Verbosity(0);

      se->registerSubsystem(gen);
    }
  if (usegun)
    {
      // PHG4ParticleGun *gun = new PHG4ParticleGun();
      // gun->set_name("anti_proton");
      // gun->set_name("geantino");
      // gun->set_vtx(0, 0, 0);
      // gun->set_mom(10, 0, 0.01);
      // gun->AddParticle("geantino",1.7776,-0.4335,0.);
      // gun->AddParticle("geantino",1.7709,-0.4598,0.);
      // gun->AddParticle("geantino",2.5621,0.60964,0.);
      // gun->AddParticle("geantino",1.8121,0.253,0.);
      // se->registerSubsystem(gun);
      PHG4ParticleGenerator *pgen = new PHG4ParticleGenerator();
      pgen->set_name("e-");
      pgen->set_z_range(0,0);
      pgen->set_eta_range(0.01,0.01);
      pgen->set_mom_range(10,10);
      pgen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
      se->registerSubsystem(pgen);
    }

  // If "readhepMC" is also set, the Upsilons will be embedded in Hijing events, if 'particles" is set, the Upsilons will be embedded in whatever particles are thrown
  if(upsilons)
    {
      // run upsilons for momentum, dca performance, alone or embedded in Hijing

      PHG4ParticleGeneratorVectorMeson *vgen = new PHG4ParticleGeneratorVectorMeson();
      vgen->add_decay_particles("e+","e-",0); // i = decay id
      // event vertex
      if (readhepmc || particles)
        {
          vgen->set_reuse_existing_vertex(true);
        }
      else
        {
          vgen->set_vtx_zrange(-10.0, +10.0);
        }

      // Note: this rapidity range completely fills the acceptance of eta = +/- 1 unit
      vgen->set_rapidity_range(-1.0, +1.0);
      vgen->set_pt_range(0.0, 10.0);

      int istate = 1;

      if(istate == 1)
        {
          // Upsilon(1S)
          vgen->set_mass(9.46);
          vgen->set_width(54.02e-6);
        }
      else if (istate == 2)
        {
          // Upsilon(2S)
          vgen->set_mass(10.0233);
          vgen->set_width(31.98e-6);
        }
      else
        {
          // Upsilon(3S)
          vgen->set_mass(10.3552);
          vgen->set_width(20.32e-6);
        }

      vgen->Verbosity(0);
      vgen->Embed(2);
      se->registerSubsystem(vgen);

      cout << "Upsilon generator for istate = " << istate << " created and registered "  << endl;
    }

  if (!readhits)
    {
      //---------------------
      // Detector description
      //---------------------

      G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
              do_svtx,do_cemc,do_hcalin,do_magnet,do_hcalout,do_pipe,
              do_FGEM,do_EGEM,do_FEMC,do_FHCAL,do_EEMC,do_DIRC,do_RICH,do_Aerogel,
              magfield_rescale);
   
    }

  //---------
  // BBC Reco
  //---------

  if (do_bbc)
    {
      gROOT->LoadMacro("G4_Bbc.C");
      BbcInit();
      Bbc_Reco();
    }

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

  if (do_svtx_cell) Svtx_Cells();

  if (do_cemc_cell) CEMC_Cells();

  if (do_hcalin_cell) HCALInner_Cells();

  if (do_hcalout_cell) HCALOuter_Cells();

  if (do_FEMC_cell) FEMC_Cells();

  if (do_FHCAL_cell) FHCAL_Cells();

  if (do_EEMC_cell) EEMC_Cells();

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

  if (do_cemc_twr) CEMC_Towers();
  if (do_cemc_cluster) CEMC_Clusters();

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

  if (do_hcalin_twr) HCALInner_Towers();
  if (do_hcalin_cluster) HCALInner_Clusters();

  if (do_hcalout_twr) HCALOuter_Towers();
  if (do_hcalout_cluster) HCALOuter_Clusters();

  //-----------------------------
  // e, h direction Calorimeter  towering and clustering
  //-----------------------------

  if (do_FEMC_twr) FEMC_Towers();
  if (do_FEMC_cluster) FEMC_Clusters();

  if (do_FHCAL_twr) FHCAL_Towers();
  if (do_FHCAL_cluster) FHCAL_Clusters();

  if (do_EEMC_twr) EEMC_Towers();
  if (do_EEMC_cluster) EEMC_Clusters();

  if (do_dst_compress) ShowerCompress();

  //--------------
  // SVTX tracking
  //--------------

  if (do_svtx_track) Svtx_Reco();

  //--------------
  // FGEM tracking
  //--------------

  if(do_FGEM_track) FGEM_FastSim_Reco();

  //--------------
  // EGEM tracking
  //--------------

  if(do_EGEM_track) EGEM_FastSim_Reco();

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

  if (do_global)
    {
      gROOT->LoadMacro("G4_Global.C");
      Global_Reco();
    }

  else if (do_global_fastsim)
    {
      gROOT->LoadMacro("G4_Global.C");
      Global_FastSim();
    }

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

  if (do_calotrigger)
    {
      gROOT->LoadMacro("G4_CaloTrigger.C");
      CaloTrigger_Sim();
    }

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

  if (do_jet_reco)
    {
      gROOT->LoadMacro("G4_Jets.C");
      Jet_Reco();
    }

  if (do_HIjetreco) {
    gROOT->LoadMacro("G4_HIJetReco.C");
    HIJetReco();
  }

  if (do_fwd_jet_reco)
    {
      gROOT->LoadMacro("G4_FwdJets.C");
      Jet_FwdReco();
    }

  //----------------------
  // Simulation evaluation
  //----------------------
  
  
  if (do_svtx_eval) Svtx_Eval(directory+"g4svtx_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");
  
  if (do_cemc_eval) CEMC_Eval(directory+"g4cemc_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");

  if (do_hcalin_eval) HCALInner_Eval(directory+"g4hcalin_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");

  if (do_hcalout_eval) HCALOuter_Eval(directory+"g4hcalout_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");

  if (do_jet_eval) Jet_Eval(directory+"g4jet_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");

  if (do_fwd_jet_eval) Jet_FwdEval(directory+"g4fwdjet_p"+pbeam+"_e"+ebeam+"_"+n+"events_eval.root");
  
  if(do_lepto_analysis){
    gROOT->LoadMacro("G4_Lepto.C");
    G4_Lepto(directory+"LeptoAna_p"+pbeam+"_e"+ebeam+"_"+n+"events");
  }
  

  
 
  
  //--------------
  // IO management
  //--------------

  if (readhits)
    {
      // Hits file
      Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
      hitsin->fileopen(inputFile);
      se->registerInputManager(hitsin);
    }
  if (readhepmc)
    {
      Fun4AllInputManager *in = new Fun4AllHepMCInputManager( "DSTIN");
      se->registerInputManager( in );
      se->fileopen( in->Name().c_str(), inputFile );
    }
  else
    {
      // for single particle generators we just need something which drives
      // the event loop, the Dummy Input Mgr does just that
      Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
      se->registerInputManager( in );
    }

  if (do_DSTReader)
    {
      //Convert DST to human command readable TTree for quick poke around the outputs
      gROOT->LoadMacro("G4_DSTReader_EICDetector.C");

      G4DSTreader_EICDetector( outputFile, //
                               /*int*/ absorberactive ,
                               /*bool*/ do_svtx ,
                               /*bool*/ do_cemc ,
                               /*bool*/ do_hcalin ,
                               /*bool*/ do_magnet ,
                               /*bool*/ do_hcalout ,
                               /*bool*/ do_cemc_twr ,
                               /*bool*/ do_hcalin_twr ,
                               /*bool*/ do_magnet  ,
                               /*bool*/ do_hcalout_twr,
                               /*bool*/ do_FGEM,
                               /*bool*/ do_EGEM,
                               /*bool*/ do_FHCAL,
                               /*bool*/ do_FHCAL_twr,
                               /*bool*/ do_FEMC,
                               /*bool*/ do_FEMC_twr,
                               /*bool*/ do_EEMC,
                               /*bool*/ do_EEMC_twr
                               );
    }

  Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
  if (do_dst_compress) DstCompress(out);
  se->registerOutputManager(out);

  //-----------------
  // Event processing
  //-----------------
  if (nEvents < 0)
    {
      return;
    }
  // if we run the particle generator and use 0 it'll run forever
  if (nEvents == 0 && !readhits && !readhepmc)
    {
      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;
    }

  se->run(nEvents);

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

  se->End();
  std::cout << "All done" << std::endl;
  delete se;
  gSystem->Exit(0);
}