emc_pi0_tbt/towInfoClustering/rundst_spiNo.C

0001 #ifndef MACRO_FUN4ALLG4CALO_C
0002 #define MACRO_FUN4ALLG4CALO_C
0003
0004 #include <GlobalVariables.C>
0005
0006 #include <DisplayOn.C>
0007 #include <G4Setup_sPHENIX.C>
0008 #include <G4_Mbd.C>
0009 #include <G4_CaloTrigger.C>
0010 #include <G4_DSTReader.C>
0011 #include <G4_Global.C>
0012 #include <G4_HIJetReco.C>
0013 #include <G4_Input.C>
0014 #include <G4_Jets.C>
0015 #include <G4_ParticleFlow.C>
0016 #include <G4_Production.C>
0017 #include <G4_TopoClusterReco.C>
0018 #include <G4_Tracking.C>
0019 #include <G4_User.C>
0020
0021 #include <fun4all/Fun4AllDstOutputManager.h>
0022 #include <fun4all/Fun4AllOutputManager.h>
0023 #include <fun4all/Fun4AllServer.h>
0024
0025 //#include <litecaloeval/LiteCaloEval.h>
0026
0027
0028 #include <phool/PHRandomSeed.h>
0029 #include <phool/recoConsts.h>
0030
0031 #include <litecaloeval/LiteCaloEval.h>
0032 #include <calib_emc_pi0/CaloCalibEmc_Pi0.h>
0033
0034 R__LOAD_LIBRARY(libfun4all.so)
0035 R__LOAD_LIBRARY(libcaloCalibDBFile.so)
0036 //R__LOAD_LIBRARY(libcalibCaloEmc_pi0.so)
0037 R__LOAD_LIBRARY(libLiteCaloEvalTowSlope.so)
0038
0039
0040 // For HepMC Hijing
0041 // try inputFile = /sphenix/sim/sim01/sphnxpro/sHijing_HepMC/sHijing_0-12fm.dat
0042
0043 int rundst_spiNo(
0044     const int nEvents = 5,
0045     //    const string &inputFile0 = "DST_CALO_G4HIT_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000004-10000.root",
0046     //const string &inputFile1 = "DST_VERTEX_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000004-10000.root",
0047     const int mdc2_4_file_num = 20,
0048     const string &outputFile = "dstnewoutput5_calo5",
0049     //      const string &inputFile0 = "/gpfs/mnt/gpfs02/sphenix/sim/sim01/sphnxpro/MDC2/sHijing_HepMC/fm_0_20/pileup/DST_CALO_G4HIT_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000003-00001.root",  //"DST_CALO_G4HIT_sHijing_0_12fm-0000000001-00000.root",
0050     //      const string &inputFile1 = "/gpfs/mnt/gpfs02/sphenix/sim/sim01/sphnxpro/MDC2/sHijing_HepMC/fm_0_20/pileup/DST_VERTEX_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000003-00001.root", //"DST_VERTEX_sHijing_0_12fm-0000000001-00000.root",
0051     //      const string &outputFile = "G4sPHENIX_calo1",
0052       const string &outdir = ".")
0053 {
0054    Fun4AllServer *se = Fun4AllServer::instance();
0055    se->Verbosity(0);
0056
0057
0058    //   string inputFile0 = "DST_CALO_G4HIT_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000004-10000.root";
0059    //   string inputFile1 = "DST_VERTEX_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000004-10000.root";
0060
0061    //   string inputFile0 = "DST_CALO_G4HIT_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-";
0062    //string inputFile0 = "DST_CALO_CLUSTER_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-00000.root";
0063    //   string inputFile1 = "DST_VERTEX_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-";
0064
0065    string inputFile0 = "DST_CALO_G4HIT_single_pi0_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-";
0066
0067    string inputFile1 =     "DST_VERTEX_single_pi0_sHijing_0_20fm_50kHz_bkg_0_20fm-0000000062-";
0068
0069
0070    int ynum_int = 100000+ mdc2_4_file_num;
0071    TString yn_tstr = "";
0072    yn_tstr += ynum_int;
0073    yn_tstr.Remove(0,1);
0074    inputFile0 += yn_tstr.Data();
0075    inputFile1 += yn_tstr.Data();
0076
0077    inputFile0 += ".root";
0078    inputFile1 += ".root";
0079
0080    cout << "running over these files" << endl;
0081    cout << inputFile0 << endl;
0082    cout << inputFile1 << endl;
0083
0084    //const string &outputFile = "newoutput1_calo1",
0085
0086
0087   //Opt to print all random seed used for debugging reproducibility. Comment out to reduce stdout prints.
0088   PHRandomSeed::Verbosity(1);
0089   /*
0090   long mtime = gSystem->Now();
0091   TString fileOut = outputFile.c_str();
0092   fileOut += (mtime / 100000 - 8542922)/3;
0093   string outputFile2 = fileOut.Data();
0094   */
0095   string outputFile2 = outputFile.c_str();
0096   outputFile2 = outputFile2 + ".root";
0097
0098   // just if we set some flags somewhere in this macro
0099   recoConsts *rc = recoConsts::instance();
0100   // By default every random number generator uses
0101   // PHRandomSeed() which reads /dev/urandom to get its seed
0102   // if the RANDOMSEED flag is set its value is taken as seed
0103   // You can either set this to a random value using PHRandomSeed()
0104   // which will make all seeds identical (not sure what the point of
0105   // this would be:
0106   //  rc->set_IntFlag("RANDOMSEED",PHRandomSeed());
0107   // or set it to a fixed value so you can debug your code
0108   //  rc->set_IntFlag("RANDOMSEED", 12345);
0109
0110   //===============
0111   // Input options
0112   //===============
0113   // verbosity setting (applies to all input managers)
0114   Input::VERBOSITY = 1;
0115   // First enable the input generators
0116   // Either:
0117   // read previously generated g4-hits files, in this case it opens a DST and skips
0118   // the simulations step completely. The G4Setup macro is only loaded to get information
0119   // about the number of layers used for the cell reco code
0120   Input::READHITS = true;
0121   INPUTREADHITS::filename[0] = inputFile0;
0122   INPUTREADHITS::filename[1] = inputFile1;
0123
0124   //-----------------
0125   // Initialize the selected Input/Event generation
0126   //-----------------
0127   // This creates the input generator(s)
0128   InputInit();
0129
0130   // register all input generators with Fun4All
0131   InputRegister();
0132
0133   // set up production relatedstuff
0134    Enable::PRODUCTION = true;
0135
0136   //======================
0137   // Write the DST
0138   //======================
0139
0140   Enable::DSTOUT = true;
0141   Enable::DSTOUT_COMPRESS = false;
0142   DstOut::OutputDir = outdir;
0143   //DstOut::OutputFile = outputFile;
0144   DstOut::OutputFile = outputFile2;
0145
0146
0147   //Option to convert DST to human command readable TTree for quick poke around the outputs
0148   //  Enable::DSTREADER = true;
0149
0150   // turn the display on (default off)
0151   Enable::DISPLAY = false;
0152
0153   //======================
0154   // What to run
0155   //======================
0156   // Global options (enabled for all enables subsystems - if implemented)
0157   //  Enable::ABSORBER = true;
0158   //  Enable::OVERLAPCHECK = true;
0159   //  Enable::VERBOSITY = 1;
0160
0161   Enable::CEMC = true;
0162   Enable::CEMC_CELL = Enable::CEMC && true;
0163   Enable::CEMC_TOWER = Enable::CEMC_CELL && true;
0164   //  Enable::CEMC_CLUSTER = Enable::CEMC_TOWER && true;
0165 //  Enable::CEMC_EVAL = Enable::CEMC_CLUSTER && true;
0166
0167   Enable::HCALIN = true;
0168   Enable::HCALIN_CELL = Enable::HCALIN && true;
0169   Enable::HCALIN_TOWER = Enable::HCALIN_CELL && true;
0170   //  Enable::HCALIN_CLUSTER = Enable::HCALIN_TOWER && true;
0171 //  Enable::HCALIN_EVAL = Enable::HCALIN_CLUSTER && true;
0172
0173   Enable::HCALOUT = true;
0174   Enable::HCALOUT_CELL = Enable::HCALOUT && true;
0175   Enable::HCALOUT_TOWER = Enable::HCALOUT_CELL && true;
0176   // Enable::HCALOUT_CLUSTER = Enable::HCALOUT_TOWER && true;
0177 //  Enable::HCALOUT_EVAL = Enable::HCALOUT_CLUSTER && true;
0178
0179
0180 //  Enable::EPD = false;
0181
0182 //  Enable::GLOBAL_RECO = true;
0183   //  Enable::GLOBAL_FASTSIM = true;
0184
0185 //  Enable::CALOTRIGGER = Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER && false;
0186
0187 //  Enable::JETS = true;
0188 //  Enable::JETS_EVAL = Enable::JETS && true;
0189
0190   // HI Jet Reco for p+Au / Au+Au collisions (default is false for
0191   // single particle / p+p-only simulations, or for p+Au / Au+Au
0192   // simulations which don't particularly care about jets)
0193 //  Enable::HIJETS = false && Enable::JETS && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER;
0194
0195   // 3-D topoCluster reconstruction, potentially in all calorimeter layers
0196 //  Enable::TOPOCLUSTER = true && Enable::CEMC_TOWER && Enable::HCALIN_TOWER && Enable::HCALOUT_TOWER;
0197   // particle flow jet reconstruction - needs topoClusters!
0198 //  Enable::PARTICLEFLOW = true && Enable::TOPOCLUSTER;
0199
0200   //---------------
0201   // Magnet Settings
0202   //---------------
0203
0204   //  const string 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)
0205   //  G4MAGNET::magfield = string(getenv("CALIBRATIONROOT")) + string("/Field/Map/sPHENIX.2d.root");  // default map from the calibration database
0206 //  G4MAGNET::magfield_rescale = 1;  // make consistent with expected Babar field strength of 1.4T
0207
0208   //---------------
0209   // Pythia Decayer
0210   //---------------
0211   // list of decay types in
0212   // $OFFLINE_MAIN/include/g4decayer/EDecayType.hh
0213   // default is All:
0214   // G4P6DECAYER::decayType = EDecayType::kAll;
0215
0216   // Initialize the selected subsystems
0217   G4Init();
0218
0219   //---------------------
0220   // GEANT4 Detector description
0221   //---------------------
0222   if (!Input::READHITS)
0223   {
0224     G4Setup();
0225   }
0226
0227   //------------------
0228   // Detector Division
0229   //------------------
0230
0231   if (Enable::MBD || Enable::MBDFAKE) Mbd_Reco();
0232
0233   if (Enable::MVTX_CELL) Mvtx_Cells();
0234   if (Enable::INTT_CELL) Intt_Cells();
0235   if (Enable::TPC_CELL) TPC_Cells();
0236   if (Enable::MICROMEGAS_CELL) Micromegas_Cells();
0237
0238   if (Enable::CEMC_CELL) CEMC_Cells();
0239
0240   if (Enable::HCALIN_CELL) HCALInner_Cells();
0241
0242   if (Enable::HCALOUT_CELL) HCALOuter_Cells();
0243
0244   //-----------------------------
0245   // CEMC towering and clustering
0246   //-----------------------------
0247
0248   if (Enable::CEMC_TOWER) CEMC_Towers();
0249   //  if (Enable::CEMC_CLUSTER) CEMC_Clusters();
0250
0251   //-----------------------------
0252   // HCAL towering and clustering
0253   //-----------------------------
0254
0255   if (Enable::HCALIN_TOWER) HCALInner_Towers();
0256   if (Enable::HCALIN_CLUSTER) HCALInner_Clusters();
0257
0258   if (Enable::HCALOUT_TOWER) HCALOuter_Towers();
0259   if (Enable::HCALOUT_CLUSTER) HCALOuter_Clusters();
0260
0261   // if enabled, do topoClustering early, upstream of any possible jet reconstruction
0262   if (Enable::TOPOCLUSTER) TopoClusterReco();
0263
0264   if (Enable::DSTOUT_COMPRESS) ShowerCompress();
0265
0266   //--------------
0267   // SVTX tracking
0268   //--------------
0269   if (Enable::MVTX_CLUSTER) Mvtx_Clustering();
0270   if (Enable::INTT_CLUSTER) Intt_Clustering();
0271   if (Enable::TPC_CLUSTER) TPC_Clustering();
0272   if (Enable::MICROMEGAS_CLUSTER) Micromegas_Clustering();
0273
0274   if (Enable::TRACKING_TRACK)
0275   {
0276     TrackingInit();
0277     Tracking_Reco();
0278   }
0279   //-----------------
0280   // Global Vertexing
0281   //-----------------
0282
0283   if (Enable::GLOBAL_RECO && Enable::GLOBAL_FASTSIM)
0284   {
0285     cout << "You can only enable Enable::GLOBAL_RECO or Enable::GLOBAL_FASTSIM, not both" << endl;
0286     gSystem->Exit(1);
0287   }
0288   if (Enable::GLOBAL_RECO)
0289   {
0290     Global_Reco();
0291   }
0292   else if (Enable::GLOBAL_FASTSIM)
0293   {
0294     Global_FastSim();
0295   }
0296
0297   //-----------------
0298   // Calo Trigger Simulation
0299   //-----------------
0300
0301   if (Enable::CALOTRIGGER)
0302   {
0303     CaloTrigger_Sim();
0304   }
0305
0306   //---------
0307   // Jet reco
0308   //---------
0309
0310   if (Enable::JETS) Jet_Reco();
0311   if (Enable::HIJETS) HIJetReco();
0312
0313   if (Enable::PARTICLEFLOW) ParticleFlow();
0314
0315   //----------------------
0316   // Simulation evaluation
0317   //----------------------
0318   string outputroot = outputFile;
0319   string remove_this = ".root";
0320   size_t pos = outputroot.find(remove_this);
0321   if (pos != string::npos)
0322   {
0323     outputroot.erase(pos, remove_this.length());
0324   }
0325
0326   //--------------
0327   // Set up Input Managers
0328   //--------------
0329
0330   InputManagers();
0331
0332   if (Enable::PRODUCTION)
0333   {
0334     Production_CreateOutputDir();
0335   }
0336
0337   if (Enable::DSTOUT)
0338   {
0339     string FullOutFile = DstOut::OutputDir + "/" + DstOut::OutputFile;
0340     Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
0341     out->AddNode("Sync");
0342     out->AddNode("EventHeader");
0343     out->AddNode("TOWER_SIM_HCALIN");
0344     out->AddNode("TOWER_RAW_HCALIN");
0345     out->AddNode("TOWER_CALIB_HCALIN");
0346     out->AddNode("CLUSTER_HCALIN");
0347     out->AddNode("TOWER_SIM_HCALOUT");
0348     out->AddNode("TOWER_RAW_HCALOUT");
0349     out->AddNode("TOWER_CALIB_HCALOUT");
0350     out->AddNode("CLUSTER_HCALOUT");
0351     out->AddNode("TOWER_SIM_CEMC");
0352     out->AddNode("TOWER_RAW_CEMC");
0353     out->AddNode("TOWER_CALIB_CEMC");
0354     out->AddNode("TOWERINFO_SIM_CEMC");
0355     out->AddNode("TOWERINFO_RAW_CEMC");
0356     out->AddNode("TOWERINFO_CALIB_CEMC");
0357     //    out->AddNode("CLUSTER_CEMC");
0358     //out->AddNode("CLUSTER_POS_COR_CEMC");
0359 // leave the topo cluster here in case we run this during pass3
0360     //out->AddNode("TOPOCLUSTER_ALLCALO");
0361     //out->AddNode("TOPOCLUSTER_EMCAL");
0362     //out->AddNode("TOPOCLUSTER_HCAL");
0363     out->AddNode("GlobalVertexMap");
0364     if (Enable::DSTOUT_COMPRESS) DstCompress(out);
0365     se->registerOutputManager(out);
0366   }
0367   //-----------------
0368   // Event processing
0369   //-----------------
0370   if (Enable::DISPLAY)
0371   {
0372     DisplayOn();
0373
0374     gROOT->ProcessLine("Fun4AllServer *se = Fun4AllServer::instance();");
0375     gROOT->ProcessLine("PHG4Reco *g4 = (PHG4Reco *) se->getSubsysReco(\"PHG4RECO\");");
0376
0377     cout << "-------------------------------------------------" << endl;
0378     cout << "You are in event display mode. Run one event with" << endl;
0379     cout << "se->run(1)" << endl;
0380     cout << "Run Geant4 command with following examples" << endl;
0381     gROOT->ProcessLine("displaycmd()");
0382
0383     return 0;
0384   }
0385
0386   string outputfile = outputFile2 + "_ho4ho_eval.root";
0387   string outputfile2 = outputFile2 + "_piemc.root";
0388   string outputfile3 = outputFile2 + "_towslopemc.root";
0389   string outputfile4 = outputFile2 + "_hcalin.root";
0390   string outputfile5 = outputFile2 + "_hcalin_mod.root";
0391   string outputfile6 = outputFile2 + "_emcal.root";
0392
0393
0394   LiteCaloEval *eval3a = new LiteCaloEval("CEMCEVALUATOR", "CEMC", outputfile3.c_str());
0395   //  LiteCaloEval *eval3a = new LiteCaloEval("HOCEMCEVALUATOR", "HCALOUT", outputfile3.c_str());
0396   //  eval->Verbosity(verbosity);
0397   eval3a->CaloType(LiteCaloEval::CEMC);
0398   //eval3a->CaloType(LiteCaloEval::HCALOUT);
0399   //eval3a->set_mode(1);
0400   se->registerSubsystem(eval3a);
0401
0402   /*
0403   //  LiteCaloEval *eval = new LiteCaloEval("CEMCEVALUATOR2", "CEMC", outputfile.c_str());
0404   LiteCaloEval *eval = new LiteCaloEval("HOCEMCEVALUATOR2", "HCALOUT", outputfile.c_str());
0405   //  eval->Verbosity(verbosity);
0406   //  eval->CaloType(LiteCaloEval::CEMC);
0407   eval->CaloType(LiteCaloEval::HCALOUT);
0408   se->registerSubsystem(eval);
0409
0410   LiteCaloEval *eval7 = new LiteCaloEval("CEMCEVALUATOR", "CEMC", outputfile6.c_str());
0411   //  LiteCaloEval *eval = new LiteCaloEval("HOCEMCEVALUATOR2", "HCALOUT", outputfile.c_str());
0412   //  eval->Verbosity(verbosity);
0413   eval7->CaloType(LiteCaloEval::CEMC);
0414   //eval->CaloType(LiteCaloEval::HCALOUT);
0415   se->registerSubsystem(eval7);
0416
0417
0418   */
0419   /*
0420   CaloCalibEmc_Pi0 *eval_pi1 = new CaloCalibEmc_Pi0("CEMC_CALIB_PI0", outputfile2);
0421   //  eval_pi1->set_mode(1);
0422   //  eval->Verbosity(verbosity);
0423   se->registerSubsystem(eval_pi1);
0424   */
0425
0426
0427   /*
0428   CaloCalibEmc_Pi0 *eval_pi2 = new CaloCalibEmc_Pi0("CEMC_CALIB_PI02", outputfile2);
0429   //  eval_pi2->set_mode(1);
0430   //  eval->Verbosity(verbosity);
0431   se->registerSubsystem(eval_pi2);
0432   cout << "successful registration of pi0 " << endl;
0433   */
0434
0435   /*
0436   LiteCaloEval *eval3 = new LiteCaloEval("HcalIn_towslope", "HCALIN", outputfile4.c_str());
0437   //  eval->Verbosity(verbosity);
0438   eval3->CaloType(LiteCaloEval::HCALIN);
0439   se->registerSubsystem(eval3);
0440
0441   LiteCaloEval *eval4 = new LiteCaloEval("HIN_towslope2", "HCALIN", outputfile5.c_str());
0442   //  eval->Verbosity(verbosity);
0443   eval4->CaloType(LiteCaloEval::HCALIN);
0444   eval4->set_mode(1);
0445   se->registerSubsystem(eval4);
0446
0447   */
0448
0449   // if we use a negative number of events we go back to the command line here
0450   if (nEvents < 0)
0451   {
0452     return 0;
0453   }
0454   // if we run the particle generator and use 0 it'll run forever
0455   if (nEvents == 0 && !Input::HEPMC && !Input::READHITS)
0456   {
0457     cout << "using 0 for number of events is a bad idea when using particle generators" << endl;
0458     cout << "it will run forever, so I just return without running anything" << endl;
0459     return 0;
0460   }
0461
0462   se->run(nEvents);
0463
0464   //-----
0465   // Exit
0466   //-----
0467
0468   se->End();
0469   //  se->PrintTimer();
0470   //se->PrintMemoryTracker();
0471   std::cout << "All done" << std::endl;
0472   gSystem->Exit(0);
0473   delete se;
0474   /*
0475   if (Enable::PRODUCTION)
0476   {
0477     Production_MoveOutput();
0478   }
0479   */
0480
0481   return 0;
0482 }
0483 #endif