RICHAnalysis_macros/macros_fun4all/Fun4All_G4_EICDetector_RICH.C

0001 int Fun4All_G4_EICDetector_RICH(
0002                                 const int nEvents = 1,
0003                 const double set_momentum = 50.0,
0004                                 const char * inputFile = "/gpfs/mnt/gpfs04/sphenix/user/nfeege/data/pythiaeRHIC/TREES/pythia.ep.20x250.RadCor=0.Q2gt1.5kevts.root",
0005                                 const char * outputFile = "G4EICDetector.root"
0006                                 )
0007 {
0008   // Set the number of TPC layer
0009   const int n_TPC_layers = 40;  // use 60 for backward compatibility only
0010
0011   //===============
0012   // Input options
0013   //===============
0014
0015   // Either:
0016   // read previously generated g4-hits files, in this case it opens a DST and skips
0017   // the simulations step completely. The G4Setup macro is only loaded to get information
0018   // about the number of layers used for the cell reco code
0019   //
0020   // In case reading production output, please double check your G4Setup_sPHENIX.C and G4_*.C consistent with those in the production macro folder
0021   // E.g. /sphenix/sim//sim01/production/2016-07-21/single_particle/spacal2d/
0022   const bool readhits = false;
0023   // Or:
0024   // read files in HepMC format (typically output from event generators like hijing or pythia)
0025   const bool readhepmc = false; // read HepMC files
0026   // Or:
0027   // read files in EICTree format generated by eicsmear package
0028   const bool readeictree = true;
0029   // Or:
0030   // Use Pythia 8
0031   const bool runpythia8 = false;
0032   // Or:
0033   // Use Pythia 6
0034   const bool runpythia6 = false;
0035   // Or:
0036   // Use HEPGen
0037   const bool runhepgen = false;
0038   // Or:
0039   // Use Sartre
0040   const bool runsartre = false;
0041
0042
0043
0044   // Besides the above flags. One can further choose to further put in following particles in Geant4 simulation
0045   // Use multi-particle generator (PHG4SimpleEventGenerator), see the code block below to choose particle species and kinematics
0046   const bool particles = false && !readhits;
0047   // or gun/ very simple single particle gun generator
0048   const bool usegun = false && !readhits;
0049   // Throw single Upsilons, may be embedded in Hijing by setting readhepmc flag also  (note, careful to set Z vertex equal to Hijing events)
0050   const bool upsilons = false && !readhits;
0051
0052   //======================
0053   // What to run
0054   //======================
0055
0056   // sPHENIX barrel
0057   bool do_bbc = false;
0058
0059   bool do_pipe = false;
0060
0061   bool do_svtx = true;
0062   bool do_svtx_cell = do_svtx && true;
0063   bool do_svtx_track = do_svtx_cell && true;
0064   bool do_svtx_eval = do_svtx_track && true;
0065
0066   bool do_pstof = false;
0067
0068   bool do_cemc = false;
0069   bool do_cemc_cell = do_cemc && true;
0070   bool do_cemc_twr = do_cemc_cell && true;
0071   bool do_cemc_cluster = do_cemc_twr && true;
0072   bool do_cemc_eval = do_cemc_cluster && false;
0073
0074   bool do_hcalin = false;
0075   bool do_hcalin_cell = do_hcalin && false;
0076   bool do_hcalin_twr = do_hcalin_cell && false;
0077   bool do_hcalin_cluster = do_hcalin_twr && false;
0078   bool do_hcalin_eval = do_hcalin_cluster && false;
0079
0080   bool do_magnet = true;
0081
0082   bool do_hcalout = false;
0083   bool do_hcalout_cell = do_hcalout && false;
0084   bool do_hcalout_twr = do_hcalout_cell && false;
0085   bool do_hcalout_cluster = do_hcalout_twr && false;
0086   bool do_hcalout_eval = do_hcalout_cluster && false;
0087
0088   // EICDetector geometry - barrel
0089   bool do_DIRC = false;
0090
0091   // EICDetector geometry - 'hadron' direction
0092   bool do_FGEM = true;
0093   bool do_FGEM_track = do_FGEM &&  true;
0094
0095   bool do_RICH = true;
0096   bool do_Aerogel = false;
0097
0098   bool do_FEMC = false;
0099   bool do_FEMC_cell = do_FEMC && true;
0100   bool do_FEMC_twr = do_FEMC_cell && true;
0101   bool do_FEMC_cluster = do_FEMC_twr && true;
0102   bool do_FEMC_eval = do_FEMC_cluster && false;
0103
0104   bool do_FHCAL = false;
0105   bool do_FHCAL_cell = do_FHCAL && true;
0106   bool do_FHCAL_twr = do_FHCAL_cell && true;
0107   bool do_FHCAL_cluster = do_FHCAL_twr && true;
0108   bool do_FHCAL_eval = do_FHCAL_cluster && false;
0109
0110   // EICDetector geometry - 'electron' direction
0111   bool do_EGEM = false;
0112   bool do_EGEM_track = do_EGEM &&  false;
0113
0114   bool do_EEMC = false;
0115   bool do_EEMC_cell = do_EEMC && false;
0116   bool do_EEMC_twr = do_EEMC_cell && false;
0117   bool do_EEMC_cluster = do_EEMC_twr && false;
0118   bool do_EEMC_eval = do_EEMC_cluster && false;
0119
0120   // Other options
0121   bool do_global = false;
0122   bool do_global_fastsim = false;
0123
0124   bool do_calotrigger = false && do_cemc_twr && do_hcalin_twr && do_hcalout_twr;
0125
0126   bool do_jet_reco = false;
0127   bool do_jet_eval = do_jet_reco && false;
0128
0129   bool do_fwd_jet_reco = false;
0130   bool do_fwd_jet_eval = do_fwd_jet_reco && false;
0131
0132   // HI Jet Reco for jet simulations in Au+Au (default is false for
0133   // single particle / p+p simulations, or for Au+Au simulations which
0134   // don't care about jets)
0135   bool do_HIjetreco = false && do_jet_reco && do_cemc_twr && do_hcalin_twr && do_hcalout_twr;
0136
0137   // Compress DST files
0138   bool do_dst_compress = false;
0139
0140   //Option to convert DST to human command readable TTree for quick poke around the outputs
0141   bool do_DSTReader = false;
0142
0143   //---------------
0144   // Load libraries
0145   //---------------
0146
0147   gSystem->Load("libfun4all.so");
0148   gSystem->Load("libg4detectors.so");
0149   gSystem->Load("libphhepmc.so");
0150   gSystem->Load("libg4testbench.so");
0151   gSystem->Load("libg4hough.so");
0152   gSystem->Load("libg4calo.so");
0153   gSystem->Load("libg4eval.so");
0154
0155   // establish the geometry and reconstruction setup
0156   gROOT->LoadMacro("G4Setup_EICDetector.C");
0157   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);
0158
0159   int absorberactive = 0; // set to 1 to make all absorbers active volumes
0160   //  const string magfield = "1.5"; // if like float -> solenoidal field in T, if string use as fieldmap name (including path)
0161   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)
0162   //const float magfield_rescale = 1.4/1.5.; //1.4/1.5; // scale the map to a 1.4 T field
0163   const float magfield_rescale = 0; // scale the map to a 0 T field
0164
0165   //---------------
0166   // Fun4All server
0167   //---------------
0168
0169   Fun4AllServer *se = Fun4AllServer::instance();
0170   se->Verbosity(1); // uncomment for batch production running with minimal output messages
0171   // se->Verbosity(Fun4AllServer::VERBOSITY_SOME); // uncomment for some info for interactive running
0172
0173   // just if we set some flags somewhere in this macro
0174   recoConsts *rc = recoConsts::instance();
0175   // By default every random number generator uses
0176   // PHRandomSeed() which reads /dev/urandom to get its seed
0177   // if the RANDOMSEED flag is set its value is taken as seed
0178   // You can either set this to a random value using PHRandomSeed()
0179   // which will make all seeds identical (not sure what the point of
0180   // this would be:
0181   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0182   // or set it to a fixed value so you can debug your code
0183   // rc->set_IntFlag("RANDOMSEED", 12345);
0184
0185   //-----------------
0186   // Event generation
0187   //-----------------
0188
0189   if (readhits)
0190     {
0191       // Get the hits from a file
0192       // The input manager is declared later
0193     }
0194   else if (readhepmc)
0195     {
0196       // this module is needed to read the HepMC records into our G4 sims
0197       // but only if you read HepMC input files
0198       HepMCNodeReader *hr = new HepMCNodeReader();
0199       se->registerSubsystem(hr);
0200     }
0201   else if (readeictree)
0202     {
0203       // this module is needed to read the EICTree style records into our G4 sims
0204       ReadEICFiles *eicr = new ReadEICFiles();
0205       eicr->OpenInputFile(inputFile);
0206       eicr->SetFirstEntry(0); // first entry to be processed in input file; start event counting at 0
0207
0208       se->registerSubsystem(eicr);
0209
0210       HepMCNodeReader *hr = new HepMCNodeReader();
0211       se->registerSubsystem(hr);
0212     }
0213   else if (runpythia8)
0214     {
0215       gSystem->Load("libPHPythia8.so");
0216
0217       PHPythia8* pythia8 = new PHPythia8();
0218       // see coresoftware/generators/PHPythia8 for example config
0219       pythia8->set_config_file("phpythia8.cfg");
0220       se->registerSubsystem(pythia8);
0221
0222       HepMCNodeReader *hr = new HepMCNodeReader();
0223       se->registerSubsystem(hr);
0224     }
0225   else if (runpythia6)
0226     {
0227       gSystem->Load("libPHPythia6.so");
0228
0229       PHPythia6 *pythia6 = new PHPythia6();
0230       // see coresoftware/generators/PHPythia6 for example config
0231       pythia6->set_config_file("phpythia6_ep.cfg");
0232       se->registerSubsystem(pythia6);
0233
0234       HepMCNodeReader *hr = new HepMCNodeReader();
0235       se->registerSubsystem(hr);
0236     }
0237   else if (runhepgen)
0238     {
0239       gSystem->Load("libsHEPGen.so");
0240
0241       sHEPGen *hepgen = new sHEPGen();
0242       // see HEPGen source directory/share/vggdata for required .dat files
0243       // see HEPGen source directory/share/datacards for required datacard files
0244       hepgen->set_datacard_file("hepgen_dvcs.data");
0245       hepgen->set_momentum_electron(-20);
0246       hepgen->set_momentum_hadron(250);
0247       se->registerSubsystem(hepgen);
0248
0249       HepMCNodeReader *hr = new HepMCNodeReader();
0250       se->registerSubsystem(hr);
0251     }
0252   else if (runsartre)
0253     {
0254       // see coresoftware/generators/PHSartre/README for setup instructions
0255       // before running:
0256       // setenv SARTRE_DIR /opt/sphenix/core/sartre-1.20_root-5.34.36
0257       gSystem->Load("libPHSartre.so");
0258
0259       PHSartre* mysartre = new PHSartre();
0260       // see coresoftware/generators/PHSartre for example config
0261       mysartre->set_config_file("sartre.cfg");
0262
0263       // particle trigger to enhance forward J/Psi -> ee
0264       PHSartreParticleTrigger* pTrig = new PHSartreParticleTrigger("MySartreTrigger");
0265       pTrig->AddParticles(-11);
0266       //pTrig->SetEtaHighLow(4.0,1.4);
0267       pTrig->SetEtaHighLow(1.0,-1.1);  // central arm
0268       pTrig->PrintConfig();
0269       mysartre->register_trigger((PHSartreGenTrigger *)pTrig);
0270       se->registerSubsystem(mysartre);
0271
0272       HepMCNodeReader *hr = new HepMCNodeReader();
0273       se->registerSubsystem(hr);
0274     }
0275
0276   // If "readhepMC" is also set, the particles will be embedded in Hijing events
0277   if(particles)
0278     {
0279       // toss low multiplicity dummy events
0280       PHG4SimpleEventGenerator *gen = new PHG4SimpleEventGenerator();
0281       gen->add_particles("pi-",1); // mu+,e+,proton,pi+,Upsilon
0282       //gen->add_particles("pi+",100); // 100 pion option
0283       if (readhepmc)
0284         {
0285           gen->set_reuse_existing_vertex(true);
0286           gen->set_existing_vertex_offset_vector(0.0, 0.0, 0.0);
0287         }
0288       else
0289         {
0290           gen->set_vertex_distribution_function(PHG4SimpleEventGenerator::Uniform,
0291                                                 PHG4SimpleEventGenerator::Uniform,
0292                                                 PHG4SimpleEventGenerator::Uniform);
0293           gen->set_vertex_distribution_mean(0.0, 0.0, 0.0);
0294           gen->set_vertex_distribution_width(0.0, 0.0, 0.0);
0295         }
0296       gen->set_vertex_size_function(PHG4SimpleEventGenerator::Uniform);
0297       gen->set_vertex_size_parameters(0.0, 0.0);
0298       gen->set_eta_range(2.0, 2.0);
0299       gen->set_phi_range(-1.0 * TMath::Pi(), 1.0 * TMath::Pi());
0300       //gen->set_pt_range(0.1, 50.0);
0301       gen->set_pt_range(0.1, 20.0);
0302       gen->Embed(1);
0303       gen->Verbosity(0);
0304
0305       se->registerSubsystem(gen);
0306     }
0307   if (usegun)
0308     {
0309       // PHG4ParticleGun *gun = new PHG4ParticleGun();
0310       // gun->set_name("anti_proton");
0311       // gun->set_name("geantino");
0312       // gun->set_vtx(0, 0, 0);
0313       // gun->set_mom(10, 0, 0.01);
0314       // gun->AddParticle("geantino",1.7776,-0.4335,0.);
0315       // gun->AddParticle("geantino",1.7709,-0.4598,0.);
0316       // gun->AddParticle("geantino",2.5621,0.60964,0.);
0317       // gun->AddParticle("geantino",1.8121,0.253,0.);
0318       // se->registerSubsystem(gun);
0319       PHG4ParticleGenerator *pgen = new PHG4ParticleGenerator();
0320       pgen->set_name("proton");
0321       pgen->set_z_range(0,0);
0322       pgen->set_eta_range(2,2);
0323       pgen->set_mom_range(set_momentum,set_momentum);
0324       //      pgen->set_mom_range(20,20);
0325       pgen->set_phi_range( -1.0 * TMath::Pi(), 1.0 * TMath::Pi());
0326       se->registerSubsystem(pgen);
0327     }
0328
0329   // 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
0330   if(upsilons)
0331     {
0332       // run upsilons for momentum, dca performance, alone or embedded in Hijing
0333
0334       PHG4ParticleGeneratorVectorMeson *vgen = new PHG4ParticleGeneratorVectorMeson();
0335       vgen->add_decay_particles("e+","e-",0); // i = decay id
0336       // event vertex
0337       if (readhepmc || particles)
0338         {
0339           vgen->set_reuse_existing_vertex(true);
0340         }
0341       else
0342         {
0343           vgen->set_vtx_zrange(-10.0, +10.0);
0344         }
0345
0346       // Note: this rapidity range completely fills the acceptance of eta = +/- 1 unit
0347       vgen->set_rapidity_range(-1.0, +1.0);
0348       vgen->set_pt_range(0.0, 10.0);
0349
0350       int istate = 1;
0351
0352       if(istate == 1)
0353         {
0354           // Upsilon(1S)
0355           vgen->set_mass(9.46);
0356           vgen->set_width(54.02e-6);
0357         }
0358       else if (istate == 2)
0359         {
0360           // Upsilon(2S)
0361           vgen->set_mass(10.0233);
0362           vgen->set_width(31.98e-6);
0363         }
0364       else
0365         {
0366           // Upsilon(3S)
0367           vgen->set_mass(10.3552);
0368           vgen->set_width(20.32e-6);
0369         }
0370
0371       vgen->Verbosity(0);
0372       vgen->Embed(2);
0373       se->registerSubsystem(vgen);
0374
0375       cout << "Upsilon generator for istate = " << istate << " created and registered "  << endl;
0376     }
0377
0378   if (!readhits)
0379     {
0380       //---------------------
0381       // Detector description
0382       //---------------------
0383
0384       G4Setup(absorberactive, magfield, TPythia6Decayer::kAll,
0385               do_svtx,do_cemc,do_hcalin,do_magnet,do_hcalout,do_pipe,
0386               do_FGEM,do_EGEM,do_FEMC,do_FHCAL,do_EEMC,do_DIRC,do_RICH,do_Aerogel,
0387               magfield_rescale);
0388
0389     }
0390
0391   //---------
0392   // BBC Reco
0393   //---------
0394
0395   if (do_bbc)
0396     {
0397       gROOT->LoadMacro("G4_Bbc.C");
0398       BbcInit();
0399       Bbc_Reco();
0400     }
0401
0402   //------------------
0403   // Detector Division
0404   //------------------
0405
0406   if (do_svtx_cell) Svtx_Cells();
0407
0408   if (do_cemc_cell) CEMC_Cells();
0409
0410   if (do_hcalin_cell) HCALInner_Cells();
0411
0412   if (do_hcalout_cell) HCALOuter_Cells();
0413
0414   if (do_FEMC_cell) FEMC_Cells();
0415
0416   if (do_FHCAL_cell) FHCAL_Cells();
0417
0418   if (do_EEMC_cell) EEMC_Cells();
0419
0420   //-----------------------------
0421   // CEMC towering and clustering
0422   //-----------------------------
0423
0424   if (do_cemc_twr) CEMC_Towers();
0425   if (do_cemc_cluster) CEMC_Clusters();
0426
0427   //-----------------------------
0428   // HCAL towering and clustering
0429   //-----------------------------
0430
0431   if (do_hcalin_twr) HCALInner_Towers();
0432   if (do_hcalin_cluster) HCALInner_Clusters();
0433
0434   if (do_hcalout_twr) HCALOuter_Towers();
0435   if (do_hcalout_cluster) HCALOuter_Clusters();
0436
0437   //-----------------------------
0438   // e, h direction Calorimeter  towering and clustering
0439   //-----------------------------
0440
0441   if (do_FEMC_twr) FEMC_Towers();
0442   if (do_FEMC_cluster) FEMC_Clusters();
0443
0444   if (do_FHCAL_twr) FHCAL_Towers();
0445   if (do_FHCAL_cluster) FHCAL_Clusters();
0446
0447   if (do_EEMC_twr) EEMC_Towers();
0448   if (do_EEMC_cluster) EEMC_Clusters();
0449
0450   if (do_dst_compress) ShowerCompress();
0451
0452   //--------------
0453   // SVTX tracking
0454   //--------------
0455
0456   if (do_svtx_track) Svtx_Reco();
0457
0458   //--------------
0459   // FGEM tracking
0460   //--------------
0461
0462   if(do_FGEM_track) FGEM_FastSim_Reco();
0463
0464   //--------------
0465   // EGEM tracking
0466   //--------------
0467
0468   if(do_EGEM_track) EGEM_FastSim_Reco();
0469
0470   //-----------------
0471   // Global Vertexing
0472   //-----------------
0473
0474   if (do_global)
0475     {
0476       gROOT->LoadMacro("G4_Global.C");
0477       Global_Reco();
0478     }
0479
0480   else if (do_global_fastsim)
0481     {
0482       gROOT->LoadMacro("G4_Global.C");
0483       Global_FastSim();
0484     }
0485
0486   //-----------------
0487   // Calo Trigger Simulation
0488   //-----------------
0489
0490   if (do_calotrigger)
0491     {
0492       gROOT->LoadMacro("G4_CaloTrigger.C");
0493       CaloTrigger_Sim();
0494     }
0495
0496   //---------
0497   // Jet reco
0498   //---------
0499
0500   if (do_jet_reco)
0501     {
0502       gROOT->LoadMacro("G4_Jets.C");
0503       Jet_Reco();
0504     }
0505
0506   if (do_HIjetreco) {
0507     gROOT->LoadMacro("G4_HIJetReco.C");
0508     HIJetReco();
0509   }
0510
0511   if (do_fwd_jet_reco)
0512     {
0513       gROOT->LoadMacro("G4_FwdJets.C");
0514       Jet_FwdReco();
0515     }
0516
0517   //----------------------
0518   // Simulation evaluation
0519   //----------------------
0520
0521   if (do_svtx_eval) Svtx_Eval(string(outputFile) + "_g4svtx_eval.root");
0522
0523   if (do_cemc_eval) CEMC_Eval(string(outputFile) + "_g4cemc_eval.root");
0524
0525   if (do_hcalin_eval) HCALInner_Eval(string(outputFile) + "_g4hcalin_eval.root");
0526
0527   if (do_hcalout_eval) HCALOuter_Eval(string(outputFile) + "_g4hcalout_eval.root");
0528
0529   if (do_FEMC_eval) FEMC_Eval(string(outputFile) + "_g4femc_eval.root");
0530
0531   if (do_FHCAL_eval) FHCAL_Eval(string(outputFile) + "_g4fhcal_eval.root");
0532
0533   if (do_EEMC_eval) EEMC_Eval(string(outputFile) + "_g4eemc_eval.root");
0534
0535   if (do_jet_eval) Jet_Eval(string(outputFile) + "_g4jet_eval.root");
0536
0537   if (do_fwd_jet_eval) Jet_FwdEval(string(outputFile) + "_g4fwdjet_eval.root");
0538
0539   //--------------
0540   // IO management
0541   //--------------
0542
0543   if (readhits)
0544     {
0545       // Hits file
0546       Fun4AllInputManager *hitsin = new Fun4AllDstInputManager("DSTin");
0547       hitsin->fileopen(inputFile);
0548       se->registerInputManager(hitsin);
0549     }
0550   if (readhepmc)
0551     {
0552       Fun4AllInputManager *in = new Fun4AllHepMCInputManager( "DSTIN");
0553       se->registerInputManager( in );
0554       se->fileopen( in->Name().c_str(), inputFile );
0555     }
0556   else
0557     {
0558       // for single particle generators we just need something which drives
0559       // the event loop, the Dummy Input Mgr does just that
0560       Fun4AllInputManager *in = new Fun4AllDummyInputManager( "JADE");
0561       se->registerInputManager( in );
0562     }
0563
0564   if (do_DSTReader)
0565     {
0566       //Convert DST to human command readable TTree for quick poke around the outputs
0567       gROOT->LoadMacro("G4_DSTReader_EICDetector.C");
0568
0569       G4DSTreader_EICDetector( outputFile, //
0570                                /*int*/ absorberactive ,
0571                                /*bool*/ do_svtx ,
0572                                /*bool*/ do_cemc ,
0573                                /*bool*/ do_hcalin ,
0574                                /*bool*/ do_magnet ,
0575                                /*bool*/ do_hcalout ,
0576                                /*bool*/ do_cemc_twr ,
0577                                /*bool*/ do_hcalin_twr ,
0578                                /*bool*/ do_magnet  ,
0579                                /*bool*/ do_hcalout_twr,
0580                                /*bool*/ do_FGEM,
0581                                /*bool*/ do_EGEM,
0582                                /*bool*/ do_FHCAL,
0583                                /*bool*/ do_FHCAL_twr,
0584                                /*bool*/ do_FEMC,
0585                                /*bool*/ do_FEMC_twr,
0586                                /*bool*/ do_EEMC,
0587                                /*bool*/ do_EEMC_twr
0588                                );
0589     }
0590
0591   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", outputFile);
0592   if (do_dst_compress) DstCompress(out);
0593   se->registerOutputManager(out);
0594
0595
0596   //-----------------
0597   // Event processing
0598   //-----------------
0599   if (nEvents < 0)
0600     {
0601       return;
0602     }
0603   // if we run the particle generator and use 0 it'll run forever
0604   if (nEvents == 0 && !readhits && !readhepmc)
0605     {
0606       cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0607       cout << "it will run forever, so I just return without running anything" << endl;
0608       return;
0609     }
0610
0611   se->run(nEvents);
0612
0613   //-----
0614   // Exit
0615   //-----
0616
0617   se->End();
0618   std::cout << "All done" << std::endl;
0619   delete se;
0620   gSystem->Exit(0);
0621 }