Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Fun4All_JetBkgd_Embed.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Fun4All_JetBkgd_Embed.C
1 #pragma once
2 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
3 #include <fun4all/SubsysReco.h>
4 #include <fun4all/Fun4AllServer.h>
5 #include <fun4all/Fun4AllInputManager.h>
6 #include <fun4all/Fun4AllSyncManager.h>
7 #include <fun4all/Fun4AllDstInputManager.h>
8 #include <fun4all/Fun4AllNoSyncDstInputManager.h>
9 
10 #include <fun4all/Fun4AllDstOutputManager.h>
11 #include <fun4all/Fun4AllOutputManager.h>
12 
13 #include <phool/PHRandomSeed.h>
14 #include <phool/recoConsts.h>
15 
16 #include <centrality/CentralityReco.h>
17 
18 #include <caloreco/CaloTowerStatus.h>
19 
20 #include <jetbackground/RetowerCEMC.h>
21 #include <jetbackground/FastJetAlgoSub.h>
22 
23 #include <jetbase/JetReco.h>
24 #include <jetbase/TowerJetInput.h>
25 #include <jetbase/FastJetAlgo.h>
26 #include <g4jets/TruthJetInput.h>
27 
28 #include <caloembedding/caloTowerEmbed.h>
29 
30 #include <jetbkgdsub/JetBkgdSub.h>
31 
32 R__LOAD_LIBRARY(libfun4all.so)
33 R__LOAD_LIBRARY(libg4jets.so)
34 R__LOAD_LIBRARY(libjetbackground.so)
35 R__LOAD_LIBRARY(libjetbase.so)
36 R__LOAD_LIBRARY(libg4dst.so)
37 R__LOAD_LIBRARY(libJetBkgdSub.so)
38 R__LOAD_LIBRARY(libcentrality_io.so)
39 R__LOAD_LIBRARY(libcentrality.so)
40 R__LOAD_LIBRARY(libCaloEmbedding.so)
41 R__LOAD_LIBRARY(libcalo_reco.so)
42 
43 
44 #endif
45 
46 
47 void Fun4All_JetBkgd_Embed(
48  const char *filelistdata = "dst_data.list",
49  const char *filelistsim = "dst_sim.list",
50  const char * output = "output",
51  const double jet_parameter = 0.4
52 )
53 {
54 
55  std::string outputbase = output;
56 
57  bool doCEMC = true;
58  bool doIHCAL = true;
59  bool doOHCAL = true;
60 
61  //-----------------------------------
62  // Fun4All server initialization
63  //-----------------------------------
64 
65  // create fun4all server
66  Fun4AllServer *se = Fun4AllServer::instance();
67  int verbosity = 1;
68  se->Verbosity(verbosity);
69  recoConsts *rc = recoConsts::instance();
70 
71  rc->set_StringFlag("CDB_GLOBALTAG","ProdA_2023");
72  rc->set_uint64Flag("TIMESTAMP",23745);
73 
74 
75  CentralityReco *cent = new CentralityReco("CentralityReco");
76  se->registerSubsystem(cent);
77 
78 
79  /*
80  // retower CEMC
81  RetowerCEMC *rcemc = new RetowerCEMC();
82  rcemc->set_towerinfo(true);
83  rcemc->Verbosity(verbosity);
84  se->registerSubsystem(rcemc);
85  */
86 
87  /*
88  //-----------------------------------
89  // Jet reco
90  //-----------------------------------
91  Enable::HIJETS_TRUTH=false;
92  Enable::HIJETS_VERBOSITY=0;
93  HIJetReco();
94 
95  // tower jets
96  // create jetreco and jettruth node names
97  string rawname = "AntiKt_Tower_r0" + to_string(int(jet_parameter * 10));
98  JetReco *towerjetreco = new JetReco();
99  towerjetreco->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER));
100  towerjetreco->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO));
101  towerjetreco->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO));
102  towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, jet_parameter), rawname);
103  towerjetreco->set_algo_node("ANTIKT");
104  towerjetreco->set_input_node("TOWER");
105  towerjetreco->Verbosity(verbosity);
106  se->registerSubsystem(towerjetreco);
107 
108  */
109 
110 
111  if(doCEMC){
112  CaloTowerStatus *statusEMC = new CaloTowerStatus("CEMCSTATUS");
113  statusEMC->set_detector_type(CaloTowerDefs::CEMC);
114  statusEMC->set_time_cut(1);
115  statusEMC->set_inputNodePrefix("TOWERINFO_CALIB_");
116  se->registerSubsystem(statusEMC);
117 
118  caloTowerEmbed *embedder_CEMC = new caloTowerEmbed("embedder_CEMC");
119  embedder_CEMC->set_detector_type(CaloTowerDefs::CEMC);
120  embedder_CEMC->set_removeBadTowers(true);
121  embedder_CEMC->set_inputNodePrefix("TOWERINFO_CALIB_");
122  se->registerSubsystem(embedder_CEMC);
123  }
124 
125  if(doIHCAL){
126  CaloTowerStatus *statusIHCAL = new CaloTowerStatus("IHCALSTATUS");
127  statusIHCAL->set_detector_type(CaloTowerDefs::HCALIN);
128  statusIHCAL->set_time_cut(2);
129  statusIHCAL->set_inputNodePrefix("TOWERINFO_CALIB_");
130  se->registerSubsystem(statusIHCAL);
131 
132  caloTowerEmbed *embedder_IHCAL = new caloTowerEmbed("embedder_IHCAL");
133  embedder_IHCAL->set_detector_type(CaloTowerDefs::HCALIN);
134  embedder_IHCAL->set_removeBadTowers(true);
135  embedder_IHCAL->set_inputNodePrefix("TOWERINFO_CALIB_");
136  se->registerSubsystem(embedder_IHCAL);
137  }
138 
139  if(doOHCAL){
140  caloTowerEmbed *embedder_OHCAL = new caloTowerEmbed("embedder_OHCal");
141  embedder_OHCAL->set_detector_type(CaloTowerDefs::HCALOUT);
142  embedder_OHCAL->set_removeBadTowers(true);
143  embedder_OHCAL->set_inputNodePrefix("TOWERINFO_CALIB_");
144  se->registerSubsystem(embedder_OHCAL);
145 
146  CaloTowerStatus *statusOHCAL = new CaloTowerStatus("OHCALSTATUS");
147  statusOHCAL->set_detector_type(CaloTowerDefs::HCALOUT);
148  statusOHCAL->set_time_cut(2);
149  statusOHCAL->set_inputNodePrefix("TOWERINFO_CALIB_");
150  se->registerSubsystem(statusOHCAL);
151  }
152 
153  // ==============
154  // Jet Bkgd Sub
155  // ==============
156  double etamin = -1.1 + jet_parameter;
157  double etamax = 1.1 - jet_parameter;
158  //data
159  JetBkgdSub *myJetTree = new JetBkgdSub(jet_parameter,outputbase + ".root");
160  if(doCEMC) myJetTree->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO));
161  if(doIHCAL) myJetTree->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO));
162  if(doOHCAL) myJetTree->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO));
163  myJetTree->doIterative(false);
164  myJetTree->doAreaSub(true);
165  myJetTree->doMultSub(false);
166  myJetTree->doTruth(false);
167  myJetTree->doData(true);
168  myJetTree->doEmbed(false);
169  myJetTree->doTowerECut(true);
170  myJetTree->setTowerThreshold(0.05);
171  myJetTree->setMinRecoPt(5.0); // only sets range for reco jets
172  myJetTree->setEtaRange(etamin, etamax);
173  myJetTree->setPtRange(10, 100); // only sets range for truth jets
174  myJetTree->Verbosity(verbosity);
175  se->registerSubsystem(myJetTree);
176 
177 
178  //embed
179  JetBkgdSub *myJetTreeEmbed = new JetBkgdSub(jet_parameter,outputbase + "_embed.root");
180  if(doCEMC) myJetTreeEmbed->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO_EMBED));
181  if(doIHCAL) myJetTreeEmbed->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO_EMBED));
182  if(doOHCAL) myJetTreeEmbed->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO_EMBED));
183  myJetTreeEmbed->doIterative(false);
184  myJetTreeEmbed->doAreaSub(true);
185  myJetTreeEmbed->doMultSub(false);
186  myJetTreeEmbed->doTruth(false);
187  myJetTreeEmbed->doData(true);
188  myJetTreeEmbed->doEmbed(true);
189  myJetTreeEmbed->doTowerECut(true);
190  myJetTreeEmbed->setTowerThreshold(0.05);
191  myJetTreeEmbed->setMinRecoPt(5.0); // only sets range for reco jets
192  myJetTreeEmbed->setEtaRange(etamin, etamax);
193  myJetTreeEmbed->setPtRange(10, 100); // only sets range for truth jets
194  myJetTreeEmbed->Verbosity(verbosity);
195  se->registerSubsystem(myJetTreeEmbed);
196 
197  //sim
198  JetBkgdSub *myJetTreeSim = new JetBkgdSub(jet_parameter,outputbase + "_sim.root");
199  if(doCEMC) myJetTreeSim->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO_SIM));
200  if(doIHCAL) myJetTreeSim->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO_SIM));
201  if(doOHCAL) myJetTreeSim->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO_SIM));
202  myJetTreeSim->doIterative(false);
203  myJetTreeSim->doAreaSub(true);
204  myJetTreeSim->doMultSub(false);
205  myJetTreeSim->doTruth(false);
206  myJetTreeSim->doData(false);
207  myJetTreeSim->doEmbed(false);
208  myJetTreeSim->doTowerECut(true);
209  myJetTreeSim->setTowerThreshold(0.05);
210  myJetTreeSim->setMinRecoPt(5.0); // only sets range for reco jets
211  myJetTreeSim->setEtaRange(etamin, etamax);
212  myJetTreeSim->setPtRange(10, 100); // only sets range for truth jets
213  myJetTreeSim->Verbosity(verbosity);
214  se->registerSubsystem(myJetTreeSim);
215 
216  //-----------------------------------
217  // Input managers
218  //-----------------------------------
219 
220  Fun4AllSyncManager *sync = se->getSyncManager();
221  sync->MixRunsOk(true);
222 
223  Fun4AllInputManager *inData = new Fun4AllDstInputManager("DSTData","DST");
224  inData->AddListFile(filelistdata);
225  se->registerInputManager(inData);
226 
227  Fun4AllInputManager *inSim = new Fun4AllNoSyncDstInputManager("DSTSim","DST","TOPSim");
228  inSim->AddListFile(filelistsim);
229  se->registerInputManager(inSim);
230  //-----------------------------------
231  // Run the analysis
232  //-----------------------------------
233 
234  se->run(100);
235  se->End();
236 
237  gSystem->Exit(0);
238  return 0;
239 
240 }