RandomConeAna/macro/Fun4All.C

0001 #pragma once
0002 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
0003
0004 #include <fun4all/SubsysReco.h>
0005 #include <fun4all/Fun4AllServer.h>
0006 #include <fun4all/Fun4AllInputManager.h>
0007 #include <fun4all/Fun4AllDstInputManager.h>
0008
0009 #include <fun4all/Fun4AllDstOutputManager.h>
0010 #include <fun4all/Fun4AllOutputManager.h>
0011
0012 #include <phool/PHRandomSeed.h>
0013 #include <phool/recoConsts.h>
0014
0015 #include <g4centrality/PHG4CentralityReco.h>
0016
0017 #include <jetbackground/RetowerCEMC.h>
0018 #include <jetbackground/FastJetAlgoSub.h>
0019
0020 #include <jetbase/JetReco.h>
0021 #include <jetbase/TowerJetInput.h>
0022 #include <jetbase/FastJetAlgo.h>
0023 #include <g4jets/TruthJetInput.h>
0024
0025 #include <HIJetReco.C>
0026
0027 #include <randomconetree/RandomConeTree.h>
0028
0029 #include <towerrho/DetermineTowerRho.h>
0030 #include <towerrho/TowerRho.h>
0031
0032 #include <randomcones/RandomConeReco.h>
0033
0034
0035
0036 // #include <randomconeana/RandomConeAna.h>
0037
0038 R__LOAD_LIBRARY(libfun4all.so)
0039 R__LOAD_LIBRARY(libg4jets.so)
0040 R__LOAD_LIBRARY(libjetbackground.so)
0041 R__LOAD_LIBRARY(libjetbase.so)
0042 R__LOAD_LIBRARY(libg4centrality.so)
0043 R__LOAD_LIBRARY(libg4dst.so)
0044 R__LOAD_LIBRARY(libRandomConeTree.so)
0045 R__LOAD_LIBRARY(libtowerrho.so)
0046 R__LOAD_LIBRARY(librandomcones.so)
0047
0048 // R__LOAD_LIBRARY(libRandomConeAna.so)
0049
0050 #endif
0051
0052
0053 void Fun4All(
0054   const char *filelistcalo = "dst_calo_cluster.list",
0055     const char *filelistglobal = "dst_global.list",
0056   const char * outputfile = "output.root",
0057   const char * outputfile_debug = "output_debug.root"
0058 )
0059 {
0060
0061   // analysis parameters
0062   //-----------------------------------
0063   // global
0064   bool do_centrality = true;
0065   double eta_max = 1.1;
0066   double cone_jet_eta_max = 0.65;
0067
0068   double tower_pT_min = 0.0;
0069   double tower_threshold = 0.05;
0070   bool do_tower_cut = true;
0071
0072   Enable::HIJETS_TRUTH=false;
0073
0074   double zvrtx_abs_cut = 10.0;
0075
0076   //-----------------------------------
0077   // rho estimation
0078   bool do_rho_Area = true;
0079   bool do_rho_Mult = true;
0080
0081   Jet::ALGO background_algo = Jet::ALGO::KT;
0082   double background_R = 0.4;
0083   int omit_nhardest = 2;
0084   double ghost_area = 0.01;
0085   double jet_min_pT = 0.0;
0086
0087   std::vector<TowerJetInput*> background_inputs =
0088   {
0089     new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER),
0090     new TowerJetInput(Jet::HCALIN_TOWERINFO),
0091     new TowerJetInput(Jet::HCALOUT_TOWERINFO)
0092   };
0093   std::vector<std::string> rho_node_names;
0094   if(do_rho_Area) rho_node_names.push_back("TowerRho_AREA");
0095   if(do_rho_Mult) rho_node_names.push_back("TowerRho_MULT");
0096
0097   //-----------------------------------
0098   // random cone reco
0099   double cone_radius = 0.4;
0100   int random_seed = 42;
0101
0102   bool do_basic_reconstruction = true;
0103   bool do_randomize_etaphi = true;
0104
0105   bool do_avoid_leading_jet = false;
0106   double lead_jet_dR = 1.4;
0107   double lead_jet_pT_threshold = 5.0; // GeV
0108
0109   std::string output_node_name = "RandomCone_Tower";
0110   std::vector<TowerJetInput*> random_cone_inputs =
0111   {
0112     new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER),
0113     new TowerJetInput(Jet::HCALIN_TOWERINFO),
0114     new TowerJetInput(Jet::HCALOUT_TOWERINFO)
0115   };
0116   std::vector<std::string> random_cone_node_names;
0117   if(do_basic_reconstruction) random_cone_node_names.push_back( output_node_name + "_basic_r0" + std::to_string(int(cone_radius*10)) );
0118   if(do_randomize_etaphi) random_cone_node_names.push_back( output_node_name + "_randomize_etaphi_r0" + std::to_string(int(cone_radius*10)) );
0119   if(do_avoid_leading_jet) random_cone_node_names.push_back( output_node_name + "_avoid_lead_jet_r0" + std::to_string(int(cone_radius*10)) );
0120
0121   std::string output_node_name_sub1 = "RandomCone_Tower_Sub1";
0122   std::vector<TowerJetInput*> random_cone_inputs_sub1 =
0123   {
0124     new TowerJetInput(Jet::CEMC_TOWERINFO_SUB1),
0125     new TowerJetInput(Jet::HCALIN_TOWERINFO_SUB1),
0126     new TowerJetInput(Jet::HCALOUT_TOWERINFO_SUB1)
0127   };
0128   if(do_basic_reconstruction) random_cone_node_names.push_back( output_node_name_sub1 + "_basic_r0" + std::to_string(int(cone_radius*10)) );
0129   if(do_randomize_etaphi) random_cone_node_names.push_back( output_node_name_sub1 + "_randomize_etaphi_r0" + std::to_string(int(cone_radius*10)) );
0130   if(do_avoid_leading_jet) random_cone_node_names.push_back( output_node_name_sub1 + "_avoid_lead_jet_r0" + std::to_string(int(cone_radius*10)) );
0131
0132   //-----------------------------------
0133
0134
0135   //-----------------------------------
0136   // Fun4All server initialization
0137   //-----------------------------------
0138
0139   // create fun4all server
0140   Fun4AllServer *se = Fun4AllServer::instance();
0141   int verbosity = 0;
0142   se->Verbosity(verbosity);
0143   recoConsts *rc = recoConsts::instance();
0144
0145   //-----------------------------------
0146   //  Centrality
0147   //-----------------------------------
0148
0149   if(do_centrality)
0150   {
0151     PHG4CentralityReco *cent = new PHG4CentralityReco();
0152     cent->Verbosity(0);
0153     cent->GetCalibrationParameters().ReadFromFile("centrality", "xml", 0, 0, string(getenv("CALIBRATIONROOT")) + string("/Centrality/"));
0154     se->registerSubsystem( cent );
0155   }
0156
0157
0158   //-----------------------------------
0159   // Jet reco
0160   //-----------------------------------
0161
0162   HIJetReco();
0163
0164   //-----------------------------------
0165   // Estimating rho
0166   //-----------------------------------
0167   DetermineTowerRho *dtr = new DetermineTowerRho();
0168   if(do_rho_Area) dtr->add_method(TowerRho::Method::AREA); // do area method
0169   if(do_rho_Mult) dtr->add_method(TowerRho::Method::MULT); // do multiplicy method
0170   dtr->set_algo(background_algo); // default = KT
0171   dtr->set_par(background_R); // default = 0.4
0172   dtr->set_tower_abs_eta(eta_max); // default=1.1
0173   dtr->set_jet_abs_eta(cone_jet_eta_max); //default is tower_abs_eta - R
0174   dtr->set_omit_nhardest(omit_nhardest); // default=2
0175   dtr->set_ghost_area(ghost_area); // default=0.01
0176   dtr->set_jet_min_pT(jet_min_pT); // default=0.0
0177   dtr->set_tower_min_pT(tower_pT_min); // default=0.0
0178   dtr->do_tower_cut(do_tower_cut); // default=false
0179   dtr->set_tower_threshold(tower_threshold); // default=0.00
0180   for(auto input : background_inputs) dtr->add_tower_input(input); // add inputs to estimate rho
0181   dtr->Verbosity(verbosity);
0182   se->registerSubsystem(dtr);
0183
0184
0185   //-----------------------------------
0186   // Random Cone Reco
0187   //-----------------------------------
0188
0189   RandomConeReco *rcr = new RandomConeReco();
0190   rcr->set_cone_radius(cone_radius); // default = 0.4
0191   rcr->set_seed(random_seed); // seed for random number generator (default is 0 (time ))
0192   rcr->set_output_node_name(output_node_name); // sets the output node prefix (default is "RandomCone")
0193   rcr->set_input_max_abs_eta(eta_max); // default (max eta for input towers)
0194   rcr->set_cone_max_abs_eta(cone_jet_eta_max); //max eta for random cone (default is input_max_abs_eta - cone_radius)
0195   rcr->set_input_min_pT(tower_pT_min); // default (min pT for input towers)
0196   rcr->do_tower_cut(do_tower_cut); // (default is false)
0197   rcr->set_tower_threshold(tower_threshold); // default (threshold for input towers)
0198   rcr->do_basic_reconstruction(do_basic_reconstruction); // (default is true)
0199   rcr->do_randomize_etaphi(do_randomize_etaphi); //  (default is false)
0200   if(do_avoid_leading_jet)
0201   {
0202     rcr->do_avoid_leading_jet(do_avoid_leading_jet); // (default is false)
0203     rcr->set_lead_jet_dR(lead_jet_dR); // defaults to 1.0 + cone_radius
0204     rcr->set_lead_jet_pT_threshold(lead_jet_pT_threshold); // defaults to 10.0
0205
0206   }
0207   for(auto input : random_cone_inputs) rcr->add_input(input); // add inputs for random cone reco
0208   rcr->Verbosity(1);
0209   se->registerSubsystem(rcr);
0210
0211   //-----------------------------------
0212   // same thing but for sub1
0213   rcr = new RandomConeReco();
0214   rcr->set_cone_radius(cone_radius); // default = 0.4
0215   rcr->set_seed(random_seed); // seed for random number generator (default is 0 (time ))
0216   rcr->set_output_node_name(output_node_name_sub1); // sets the output node prefix (default is "RandomCone")
0217   rcr->set_input_max_abs_eta(eta_max); // default (max eta for input towers)
0218   rcr->set_cone_max_abs_eta(cone_jet_eta_max); //max eta for random cone (default is input_max_abs_eta - cone_radius)
0219   rcr->set_input_min_pT(tower_pT_min); // default (min pT for input towers)
0220   rcr->do_tower_cut(false); // (default is false)
0221   rcr->set_tower_threshold(0.0); // default (threshold for input towers)
0222   rcr->do_basic_reconstruction(do_basic_reconstruction); // (default is true)
0223   rcr->do_randomize_etaphi(do_randomize_etaphi); //  (default is false)
0224   if(do_avoid_leading_jet)
0225   {
0226     rcr->do_avoid_leading_jet(do_avoid_leading_jet); // (default is false)
0227     rcr->set_lead_jet_dR(lead_jet_dR); // defaults to 1.0 + cone_radius
0228     rcr->set_lead_jet_pT_threshold(lead_jet_pT_threshold); // defaults to 10.0
0229   }
0230   for(auto input : random_cone_inputs_sub1) rcr->add_input(input); // add inputs for random cone reco
0231   rcr->Verbosity(1);
0232   se->registerSubsystem(rcr);
0233
0234
0235
0236
0237   //-----------------------------------
0238
0239
0240   // Random Cone Analysis
0241   //-----------------------------------
0242   RandomConeTree * rca = new RandomConeTree(outputfile);
0243   for(auto node_name : random_cone_node_names) rca->add_random_cone_node(node_name); // add random cone node to output tree
0244   for(auto rho_node_name : rho_node_names) rca->add_rho_node(rho_node_name); // add rho node to output tree
0245   rca->add_gvtx_cut(-zvrtx_abs_cut, zvrtx_abs_cut); // default is no cut
0246   rca->do_centrality_info(do_centrality); // default is false (centrality info)
0247   // rca->do_event_selection_on_leading_truth_jet(30,1000); // default is false (event selection on leading truth jet) MC
0248   // rca->add_weight_to_ttree(1.0) // add weight to output tree
0249   rca->do_data(false); // default is false (data)
0250   rca->Verbosity(1);
0251   se->registerSubsystem(rca);
0252
0253
0254   // manual reco
0255   // RandomConeAna *myJetTree = new RandomConeAna(outputfile_debug);
0256   // myJetTree->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER));
0257   // myJetTree->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO));
0258   // myJetTree->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO));
0259   // myJetTree->add_iter_input(new TowerJetInput(Jet::CEMC_TOWERINFO_SUB1));
0260   // myJetTree->add_iter_input(new TowerJetInput(Jet::HCALIN_TOWERINFO_SUB1));
0261   // myJetTree->add_iter_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO_SUB1));
0262   // myJetTree->set_user_seed(42);
0263   // // myJetTree->setEventSelection(30,1000);
0264   // // myJetTree->addWeight(2.178e-9);
0265   // myJetTree->Verbosity(verbosity);
0266   // se->registerSubsystem(myJetTree);
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276   //-----------------------------------
0277   // Input managers
0278   //-----------------------------------
0279
0280   Fun4AllInputManager *in2 = new Fun4AllDstInputManager("DSTcalo");
0281   in2->AddListFile(filelistcalo,1);
0282   se->registerInputManager(in2);
0283
0284   Fun4AllInputManager *in3 = new Fun4AllDstInputManager("DSTglobal");
0285   in3->AddListFile(filelistglobal,1);
0286   se->registerInputManager(in3);
0287
0288   //-----------------------------------
0289   // Run the analysis
0290   //-----------------------------------
0291
0292   se->run(10);
0293   se->End();
0294
0295   gSystem->Exit(0);
0296   return 0;
0297
0298 }