Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_FullJetFinder.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_FullJetFinder.C
1 #include <fun4all/SubsysReco.h>
2 #include <fun4all/Fun4AllServer.h>
3 #include <fun4all/Fun4AllInputManager.h>
4 #include <fun4all/Fun4AllDstInputManager.h>
5 
6 #include <fun4all/Fun4AllDstOutputManager.h>
7 #include <fun4all/Fun4AllOutputManager.h>
8 #include <fun4all/Fun4AllServer.h>
9 
10 #include <jetbase/FastJetAlgo.h>
11 #include <jetbase/JetReco.h>
12 #include <g4jets/TruthJetInput.h>
13 #include <jetbackground/RetowerCEMC.h>
14 
15 #include <phool/PHRandomSeed.h>
16 #include <phool/recoConsts.h>
17 
18 #include <TSystem.h>
19 
20 #include <g4centrality/PHG4CentralityReco.h>
21 #include <string>
22 
23 //#include "HIJetReco.C"
24 
25 #include <caloreco/RawClusterBuilderTopo.h>
26 #include <particleflowreco/ParticleFlowReco.h>
27 #include <particleflowreco/ParticleFlowJetInput.h>
28 
29 
30 //#include <jetvertextagging/JetVertexTagging.h>
31 #pragma GCC diagnostic push
32 #pragma GCC diagnostic ignored "-Wundefined-internal"
33 #include <FullJetFinder.h>
34 #pragma GCC diagnostic pop
35 
36 #include <fstream>
37 
38 #include <globalvertex/GlobalVertexReco.h>
39 
40 R__LOAD_LIBRARY(libfun4all.so)
41 R__LOAD_LIBRARY(libjetbase.so)
42 R__LOAD_LIBRARY(libjetbackground.so)
43 R__LOAD_LIBRARY(libFullJetFinder.so)
44 R__LOAD_LIBRARY(libg4centrality.so)
45 R__LOAD_LIBRARY(libg4dst.so)
46 R__LOAD_LIBRARY(libcalo_reco.so)
47 R__LOAD_LIBRARY(libparticleflow.so)
48 R__LOAD_LIBRARY(libglobalvertex.so)
49 
50 
51 //void Fun4All_JetVal(const char *filelisttruth = "dst_truth_jet.list",
52 // const char *filelistcalo = "dst_calo_cluster.list",
53 // const char *outname = "outputest.root")
54 void Fun4All_FullJetFinder(std::string outDir = "./", std::vector<std::string> myInputLists = {"condorJob/fileLists/productionFiles-CHARM-dst_tracks-00000.LIST"}, const int nEvents = 1){
55 
56 
57  Fun4AllServer *se = Fun4AllServer::instance();
58  int verbosity = 0;
59 
60  bool whichR[7] = {true, false, true,false, true, false, false};
61 
62  //std::string outDir = "./";
63 
64  //std::string outDir = "/sphenix/tg/tg01/hf/jkvapil/JET30_r8/";
65  if (outDir.substr(outDir.size() - 1, 1) != "/") outDir += "/";
66  outDir += "myTestJets";
67 
68 
69 
70  std::string fileNumber = myInputLists[0];
71  size_t findLastDash = fileNumber.find_last_of("-");
72  if (findLastDash != std::string::npos) fileNumber.erase(0, findLastDash + 1);
73  std::string remove_this = ".list";
74  size_t pos = fileNumber.find(remove_this);
75  if (pos != std::string::npos) fileNumber.erase(pos, remove_this.length());
76  std::string outputFileName = "outputData_" + fileNumber + ".root";
77 
78  std::string outputRecoDir = outDir + "/inReconstruction/";
79  std::string makeDirectory = "mkdir -p " + outputRecoDir;
80  system(makeDirectory.c_str());
81  std::string outputRecoFile = outputRecoDir + outputFileName;
82 
83  //Create the server
84  se->Verbosity(verbosity);
85 
86  //Add all required input files
87  for (unsigned int i = 0; i < myInputLists.size(); ++i)
88  {
89  Fun4AllInputManager *infile = new Fun4AllDstInputManager("DSTin_" + std::to_string(i));
90  infile->AddListFile(myInputLists[i]);
91  se->registerInputManager(infile);
92  }
93  //recoConsts *rc = recoConsts::instance();
94 
95  PHG4CentralityReco *cent = new PHG4CentralityReco();
96  cent->Verbosity(0);
97  cent->GetCalibrationParameters().ReadFromFile("centrality", "xml", 0, 0, std::string(getenv("CALIBRATIONROOT")) + std::string("/Centrality/"));
98  se->registerSubsystem( cent );
99 
100  /* GlobalVertexReco* gblvertex = new GlobalVertexReco();
101  gblvertex->Verbosity(0);
102  se->registerSubsystem(gblvertex);*/
103 
104  //HIJetReco();
105  //JetReco *truthjetreco = new JetReco("JetReco_truth",JetReco::TRANSITION::JET_MAP);
106  JetReco *truthjetreco = new JetReco("JetReco_truth");
107  TruthJetInput *tji = new TruthJetInput(Jet::PARTICLE);
108  //tji->Verbosity(10);
109  tji->add_embedding_flag(1); // changes depending on signal vs. embedded
110  truthjetreco->add_input(tji);
111  /* 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
112  if(whichR[1])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.3), "AntiKt_Truth_r03");
113  if(whichR[2])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.4), "AntiKt_Truth_r04");
114  if(whichR[3])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.5), "AntiKt_Truth_r05");
115  if(whichR[4])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.6), "AntiKt_Truth_r06");
116  if(whichR[5])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.7), "AntiKt_Truth_r07");
117  if(whichR[6])truthjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.8), "AntiKt_Truth_r08");*/
118  if(whichR[0])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.2, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r02");
119  if(whichR[1])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.3, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r03");
120  if(whichR[2])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.4, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r04");
121  if(whichR[3])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.5, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r05");
122  if(whichR[4])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.6, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r06");
123  if(whichR[5])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.7, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r07");
124  if(whichR[6])truthjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.8, VERBOSITY, verbosity, CALC_AREA}}), "C_AntiKt_Truth_r08");
125  truthjetreco->set_algo_node("ANTIKT");
126  truthjetreco->set_input_node("TRUTH");
127  truthjetreco->Verbosity(0);
128  se->registerSubsystem(truthjetreco);
129 
130  RetowerCEMC *rcemc = new RetowerCEMC();
131  rcemc->Verbosity(verbosity);
132  rcemc->set_towerinfo(true);
133  se->registerSubsystem(rcemc);
134 
135 
136  RawClusterBuilderTopo* ClusterBuilder1 = new RawClusterBuilderTopo("HcalRawClusterBuilderTopo1");
137  ClusterBuilder1->Verbosity(verbosity);
138  ClusterBuilder1->set_nodename("TOPOCLUSTER_EMCAL");
139  ClusterBuilder1->set_enable_HCal(false);
140  ClusterBuilder1->set_enable_EMCal(true);
141  ClusterBuilder1->set_noise(0.0025, 0.006, 0.03);
142  ClusterBuilder1->set_significance(4.0, 2.0, 0.0);
143  ClusterBuilder1->allow_corner_neighbor(true);
144  ClusterBuilder1->set_do_split(true);
145  ClusterBuilder1->set_minE_local_max(1.0, 2.0, 0.5);
146  ClusterBuilder1->set_R_shower(0.025);
147  se->registerSubsystem(ClusterBuilder1);
148 
149  RawClusterBuilderTopo* ClusterBuilder2 = new RawClusterBuilderTopo("HcalRawClusterBuilderTopo2");
150  ClusterBuilder2->Verbosity(verbosity);
151  ClusterBuilder2->set_nodename("TOPOCLUSTER_HCAL");
152  ClusterBuilder2->set_enable_HCal(true);
153  ClusterBuilder2->set_enable_EMCal(false);
154  ClusterBuilder2->set_noise(0.0025, 0.006, 0.03);
155  ClusterBuilder2->set_significance(4.0, 2.0, 0.0);
156  ClusterBuilder2->allow_corner_neighbor(true);
157  ClusterBuilder2->set_do_split(true);
158  ClusterBuilder2->set_minE_local_max(1.0, 2.0, 0.5);
159  ClusterBuilder2->set_R_shower(0.025);
160  se->registerSubsystem(ClusterBuilder2);
161 
162  ParticleFlowReco *pfr = new ParticleFlowReco();
163  pfr->set_energy_match_Nsigma(1.5);
164  pfr->Verbosity(verbosity);
165  se->registerSubsystem(pfr);
166 
167 
168 
169  JetReco *towerjetreco = new JetReco();
170  towerjetreco->add_input(new ParticleFlowJetInput());
171  /*if(whichR[0])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.2, verbosity), "AntiKt_reco_r02");
172  if(whichR[1])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.3, verbosity), "AntiKt_reco_r03");
173  if(whichR[2])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.4, verbosity), "AntiKt_reco_r04");
174  if(whichR[3])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.5, verbosity), "AntiKt_reco_r05");
175  if(whichR[4])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.6, verbosity), "AntiKt_reco_r06");
176  if(whichR[5])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.7, verbosity), "AntiKt_reco_r07");
177  if(whichR[6])towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, 0.8, verbosity), "AntiKt_reco_r08");*/
178  if(whichR[0])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.2, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r02");
179  if(whichR[1])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.3, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r03");
180  if(whichR[2])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.4, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r04");
181  if(whichR[3])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.5, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r05");
182  if(whichR[4])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.6, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r06");
183  if(whichR[5])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.7, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r07");
184  if(whichR[6])towerjetreco->add_algo(new FastJetAlgo({{Jet::ANTIKT, JET_R, 0.8, VERBOSITY, verbosity, CALC_AREA}}), "AntiKt_reco_r08");
185  towerjetreco->set_algo_node("ANTIKT");
186  towerjetreco->set_input_node("INCLUSIVE_RECO");
187  towerjetreco->Verbosity(verbosity);
188  se->registerSubsystem(towerjetreco);
189 
190  FullJetFinder *myJetVal = new FullJetFinder(outputRecoFile);//{"AntiKt_r02","AntiKt_r03","AntiKt_r04","AntiKt_r05","AntiKt_r06"});
191  if(whichR[0])myJetVal->add_input("AntiKt_reco_r02", "C_AntiKt_Truth_r02","AntiKt_r02");
192  if(whichR[1])myJetVal->add_input("AntiKt_reco_r03", "C_AntiKt_Truth_r03","AntiKt_r03");
193  if(whichR[2])myJetVal->add_input("AntiKt_reco_r04", "C_AntiKt_Truth_r04","AntiKt_r04");
194  if(whichR[3])myJetVal->add_input("AntiKt_reco_r05", "C_AntiKt_Truth_r05","AntiKt_r05");
195  if(whichR[4])myJetVal->add_input("AntiKt_reco_r06", "C_AntiKt_Truth_r06","AntiKt_r06");
196  if(whichR[5])myJetVal->add_input("AntiKt_reco_r07", "C_AntiKt_Truth_r07","AntiKt_r07");
197  if(whichR[6])myJetVal->add_input("AntiKt_reco_r08", "C_AntiKt_Truth_r08","AntiKt_r08");
198  myJetVal->doFiducialAcceptance(true);
199  myJetVal->setPtRangeReco(5, 100);
200  myJetVal->setPtRangeTruth(10, 100);
201  myJetVal->doTruth(true);
202  se->registerSubsystem(myJetVal);
203 
204 
205 
206  se->run(nEvents);
207  se->End();
208 
209  std::ifstream file(outputRecoFile.c_str());
210  if (file.good())
211  {
212  std::string moveOutput = "mv " + outputRecoFile + " " + outDir;
213  system(moveOutput.c_str());
214  }
215 
216  gSystem->Exit(0);
217  return;
218 
219 }