FullJetFinder/macro/Fun4All_FullJetFinder.C

0001 #include <fun4all/SubsysReco.h>
0002 #include <fun4all/Fun4AllServer.h>
0003 #include <fun4all/Fun4AllInputManager.h>
0004 #include <fun4all/Fun4AllDstInputManager.h>
0005
0006 #include <fun4all/Fun4AllDstOutputManager.h>
0007 #include <fun4all/Fun4AllOutputManager.h>
0008 #include <fun4all/Fun4AllServer.h>
0009
0010 #include <jetbase/FastJetAlgo.h>
0011 #include <jetbase/JetReco.h>
0012 #include <g4jets/TruthJetInput.h>
0013 #include <jetbackground/RetowerCEMC.h>
0014 #include <jetbase/TrackJetInput.h>
0015
0016 #include <phool/PHRandomSeed.h>
0017 #include <phool/recoConsts.h>
0018
0019 #include <TSystem.h>
0020
0021 #include <g4centrality/PHG4CentralityReco.h>
0022 #include <string>
0023
0024 //#include "HIJetReco.C"
0025
0026 #include <caloreco/RawClusterBuilderTopo.h>
0027 #include <particleflowreco/ParticleFlowReco.h>
0028 #include <particleflowreco/ParticleFlowJetInput.h>
0029
0030
0031 //#include <jetvertextagging/JetVertexTagging.h>
0032 #pragma GCC diagnostic push
0033 #pragma GCC diagnostic ignored "-Wundefined-internal"
0034 #include <FullJetFinder.h>
0035 #pragma GCC diagnostic pop
0036
0037 #include <fstream>
0038
0039 #include <globalvertex/GlobalVertexReco.h>
0040
0041 R__LOAD_LIBRARY(libfun4all.so)
0042 R__LOAD_LIBRARY(libjetbase.so)
0043 R__LOAD_LIBRARY(libjetbackground.so)
0044 R__LOAD_LIBRARY(libFullJetFinder.so)
0045 R__LOAD_LIBRARY(libg4centrality.so)
0046 R__LOAD_LIBRARY(libg4dst.so)
0047 R__LOAD_LIBRARY(libcalo_reco.so)
0048 R__LOAD_LIBRARY(libparticleflow.so)
0049 R__LOAD_LIBRARY(libglobalvertex.so)
0050 R__LOAD_LIBRARY(libg4jets.so)
0051
0052
0053 //void Fun4All_JetVal(const char *filelisttruth = "dst_truth_jet.list",
0054 //                     const char *filelistcalo = "dst_calo_cluster.list",
0055 //                     const char *outname = "outputest.root")
0056 void Fun4All_FullJetFinder(std::string outDir = "./", std::vector<std::string> myInputLists = {"condorJob/fileLists/productionFiles-CHARM-dst_tracks-00000.LIST"}, const int nEvents = 1){
0057
0058
0059   Fun4AllServer *se = Fun4AllServer::instance();
0060   int verbosity = 0;
0061
0062   bool whichR[7] = {false, false, true,false, false, false, false};
0063
0064   FullJetFinder::TYPE jet_type = FullJetFinder::TYPE::CHARGEJET;
0065
0066   //std::string outDir = "./";
0067
0068   //std::string outDir = "/sphenix/tg/tg01/hf/jkvapil/JET30_r8/";
0069   if (outDir.substr(outDir.size() - 1, 1) != "/") outDir += "/";
0070   outDir += "myTestJets";
0071
0072
0073
0074   std::string fileNumber = myInputLists[0];
0075   size_t findLastDash = fileNumber.find_last_of("-");
0076   if (findLastDash != std::string::npos) fileNumber.erase(0, findLastDash + 1);
0077   std::string remove_this = ".list";
0078   size_t pos = fileNumber.find(remove_this);
0079   if (pos != std::string::npos) fileNumber.erase(pos, remove_this.length());
0080   std::string outputFileName = "outputData_" + fileNumber + ".root";
0081
0082   std::string outputRecoDir = outDir + "/inReconstruction/";
0083   std::string makeDirectory = "mkdir -p " + outputRecoDir;
0084   system(makeDirectory.c_str());
0085   std::string outputRecoFile = outputRecoDir + outputFileName;
0086
0087   //Create the server
0088   se->Verbosity(verbosity);
0089
0090   //Add all required input files
0091   for (unsigned int i = 0; i < myInputLists.size(); ++i)
0092   {
0093     Fun4AllInputManager *infile = new Fun4AllDstInputManager("DSTin_" + std::to_string(i));
0094     infile->AddListFile(myInputLists[i]);
0095     se->registerInputManager(infile);
0096   }
0097   //recoConsts *rc = recoConsts::instance();
0098
0099   PHG4CentralityReco *cent = new PHG4CentralityReco();
0100   cent->Verbosity(0);
0101   cent->GetCalibrationParameters().ReadFromFile("centrality", "xml", 0, 0, std::string(getenv("CALIBRATIONROOT")) + std::string("/Centrality/"));
0102   se->registerSubsystem( cent );
0103
0104   /* GlobalVertexReco* gblvertex = new GlobalVertexReco();
0105   gblvertex->Verbosity(0);
0106   se->registerSubsystem(gblvertex);*/
0107
0108   //HIJetReco();
0109   //JetReco *truthjetreco = new JetReco("JetReco_truth",JetReco::TRANSITION::JET_MAP);
0110   JetReco *truthjetreco = new JetReco("JetReco_truth");
0111       TruthJetInput *tji = new TruthJetInput(Jet::PARTICLE);
0112       //tji->Verbosity(10);
0113       tji->add_embedding_flag(1);  // changes depending on signal vs. embedded
0114       truthjetreco->add_input(tji);
0115      /* if(whichR[0])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.2), "AntiKt_Truth_r02");  actually you cannot rename it, or you get empty trees...sad
0116       if(whichR[1])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.3), "AntiKt_Truth_r03");
0117       if(whichR[2])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.4), "AntiKt_Truth_r04");
0118       if(whichR[3])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.5), "AntiKt_Truth_r05");
0119       if(whichR[4])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.6), "AntiKt_Truth_r06");
0120       if(whichR[5])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.7), "AntiKt_Truth_r07");
0121       if(whichR[6])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.8), "AntiKt_Truth_r08");*/
0122       if(whichR[0])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.2, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r02");
0123       if(whichR[1])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.3, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r03");
0124       if(whichR[2])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.4, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r04");
0125       if(whichR[3])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.5, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r05");
0126       if(whichR[4])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.6, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r06");
0127       if(whichR[5])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.7, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r07");
0128       if(whichR[6])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.8, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r08");
0129       truthjetreco->set_algo_node("ANTIKT");
0130       truthjetreco->set_input_node("TRUTH");
0131       truthjetreco->Verbosity(0);
0132       se->registerSubsystem(truthjetreco);
0133
0134   if(jet_type == FullJetFinder::TYPE::FULLJET){
0135
0136         RetowerCEMC *rcemc = new RetowerCEMC();
0137   rcemc->Verbosity(verbosity);
0138   rcemc->set_towerinfo(true);
0139  se->registerSubsystem(rcemc);
0140
0141
0142   RawClusterBuilderTopo* ClusterBuilder1 = new RawClusterBuilderTopo("HcalRawClusterBuilderTopo1");
0143   ClusterBuilder1->Verbosity(verbosity);
0144   ClusterBuilder1->set_nodename("TOPOCLUSTER_EMCAL");
0145   ClusterBuilder1->set_enable_HCal(false);
0146   ClusterBuilder1->set_enable_EMCal(true);
0147   ClusterBuilder1->set_noise(0.0025, 0.006, 0.03);
0148   ClusterBuilder1->set_significance(4.0, 2.0, 0.0);
0149   ClusterBuilder1->allow_corner_neighbor(true);
0150   ClusterBuilder1->set_do_split(true);
0151   ClusterBuilder1->set_minE_local_max(1.0, 2.0, 0.5);
0152   ClusterBuilder1->set_R_shower(0.025);
0153   se->registerSubsystem(ClusterBuilder1);
0154
0155   RawClusterBuilderTopo* ClusterBuilder2 = new RawClusterBuilderTopo("HcalRawClusterBuilderTopo2");
0156   ClusterBuilder2->Verbosity(verbosity);
0157   ClusterBuilder2->set_nodename("TOPOCLUSTER_HCAL");
0158   ClusterBuilder2->set_enable_HCal(true);
0159   ClusterBuilder2->set_enable_EMCal(false);
0160   ClusterBuilder2->set_noise(0.0025, 0.006, 0.03);
0161   ClusterBuilder2->set_significance(4.0, 2.0, 0.0);
0162   ClusterBuilder2->allow_corner_neighbor(true);
0163   ClusterBuilder2->set_do_split(true);
0164   ClusterBuilder2->set_minE_local_max(1.0, 2.0, 0.5);
0165   ClusterBuilder2->set_R_shower(0.025);
0166   se->registerSubsystem(ClusterBuilder2);
0167
0168   ParticleFlowReco *pfr = new ParticleFlowReco();
0169   pfr->set_energy_match_Nsigma(1.5);
0170   pfr->Verbosity(verbosity);
0171   se->registerSubsystem(pfr);
0172
0173   }
0174
0175   JetReco *towerjetreco = new JetReco();
0176   if(jet_type == FullJetFinder::TYPE::FULLJET)towerjetreco->add_input(new ParticleFlowJetInput());
0177   if(jet_type == FullJetFinder::TYPE::CHARGEJET)towerjetreco->add_input(new TrackJetInput(Jet::TRACK));
0178   /*if(whichR[0])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.2, verbosity), "AntiKt_reco_r02");
0179   if(whichR[1])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.3, verbosity), "AntiKt_reco_r03");
0180   if(whichR[2])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.4, verbosity), "AntiKt_reco_r04");
0181   if(whichR[3])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.5, verbosity), "AntiKt_reco_r05");
0182   if(whichR[4])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.6, verbosity), "AntiKt_reco_r06");
0183   if(whichR[5])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.7, verbosity), "AntiKt_reco_r07");
0184   if(whichR[6])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.8, verbosity), "AntiKt_reco_r08");*/
0185   if(whichR[0])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.2, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r02");
0186   if(whichR[1])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.3, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r03");
0187   if(whichR[2])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.4, VERBOSITY, verbosity, CALC_AREA, CONSTITUENT_MIN_PT, 0.3}}), "AntiKt_reco_r04");
0188   if(whichR[3])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.5, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r05");
0189   if(whichR[4])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.6, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r06");
0190   if(whichR[5])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.7, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r07");
0191   if(whichR[6])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.8, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r08");
0192   towerjetreco->set_algo_node("ANTIKT");
0193   towerjetreco->set_input_node("INCLUSIVE_RECO");
0194   towerjetreco->Verbosity(verbosity);
0195   se->registerSubsystem(towerjetreco);
0196
0197   FullJetFinder *myJetVal = new FullJetFinder(outputRecoFile,jet_type);//{"AntiKt_r02","AntiKt_r03","AntiKt_r04","AntiKt_r05","AntiKt_r06"});
0198   if(whichR[0])myJetVal->add_input("AntiKt_reco_r02", "C_AntiKt_Truth_r02","AntiKt_r02");
0199   if(whichR[1])myJetVal->add_input("AntiKt_reco_r03", "C_AntiKt_Truth_r03","AntiKt_r03");
0200   if(whichR[2])myJetVal->add_input("AntiKt_reco_r04", "C_AntiKt_Truth_r04","AntiKt_r04");
0201   if(whichR[3])myJetVal->add_input("AntiKt_reco_r05", "C_AntiKt_Truth_r05","AntiKt_r05");
0202   if(whichR[4])myJetVal->add_input("AntiKt_reco_r06", "C_AntiKt_Truth_r06","AntiKt_r06");
0203   if(whichR[5])myJetVal->add_input("AntiKt_reco_r07", "C_AntiKt_Truth_r07","AntiKt_r07");
0204   if(whichR[6])myJetVal->add_input("AntiKt_reco_r08", "C_AntiKt_Truth_r08","AntiKt_r08");
0205   myJetVal->doFiducialAcceptance(true);
0206   myJetVal->setPtRangeReco(5, 100);
0207   myJetVal->setPtRangeTruth(10, 100);
0208   myJetVal->doTruth(true);
0209   se->registerSubsystem(myJetVal);
0210
0211
0212
0213   se->run(nEvents);
0214   se->End();
0215
0216   std::ifstream file(outputRecoFile.c_str());
0217   if (file.good())
0218   {
0219     std::string moveOutput = "mv " + outputRecoFile + " " + outDir;
0220     system(moveOutput.c_str());
0221   }
0222
0223   gSystem->Exit(0);
0224   return;
0225
0226 }