Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MyJetAnalysis.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MyJetAnalysis.cc
1 #include "MyJetAnalysis.h"
2 #include <vector>
3 #include <fun4all/Fun4AllReturnCodes.h>
4 #include <fun4all/Fun4AllServer.h>
5 #include <fun4all/PHTFileServer.h>
6 //#include <sPHAnalysis.h>
7 
8 #include <HepMC/GenVertex.h>
9 #include <HepMC/IteratorRange.h>
10 #include <HepMC/SimpleVector.h>
11 
12 #include <phool/PHCompositeNode.h>
13 #include <phool/getClass.h>
14 
15 //#include <g4eval/JetEvalStack.h>
16 
17 #include <trackbase_historic/SvtxTrackMap.h>
18 #include <g4jets/FastJetAlgo.h>
19 #include <jetbackground/FastJetAlgoSub.h>
20 #include <g4jets/JetMap.h>
21 #include <g4jets/Jetv1.h>
22 
23 #include <g4vertex/GlobalVertex.h>
24 #include <g4vertex/GlobalVertexMap.h>
25 
26 #include <calobase/RawTowerContainer.h>
27 #include <calobase/RawTower.h>
28 #include <calobase/RawTowerGeomContainer.h>
29 #include <calobase/RawTowerGeom.h>
30 #include <calobase/TowerInfoContainer.h>
31 #include <calobase/TowerInfo.h>
32 
33 #include <g4main/PHG4Particle.h>
34 #include <g4main/PHG4TruthInfoContainer.h>
35 
36 #include <g4eval/CaloRawTowerEval.h>
37 #include <g4eval/CaloEvalStack.h>
38 
39 #include <ffaobjects/EventHeader.h>
40 #include <ffaobjects/EventHeaderv1.h>
41 #include <centrality/CentralityInfo.h>
42 // fastjet includes
43 #include <fastjet/ClusterSequence.hh>
44 #include <fastjet/JetDefinition.hh>
45 #include <fastjet/PseudoJet.hh>
46 #include <jetbackground/TowerBackground.h>
47 #include <phhepmc/PHHepMCGenEvent.h>
48 #include <phhepmc/PHHepMCGenEventMap.h>
49 
50 #include <TFile.h>
51 #include <TH1F.h>
52 #include <TH2F.h>
53 #include <TString.h>
54 #include <TTree.h>
55 #include <TVector3.h>
56 #include <TLorentzVector.h>
57 #include <TArrayF.h>
58 
59 #include <fstream>
60 #include <algorithm>
61 #include <cassert>
62 #include <cmath>
63 #include <iostream>
64 #include <limits>
65 #include <stdexcept>
66 //#include <memory>
67 
68 double DeltaPhi(double phi1, double phi2){
69  double_t dPhi_temp = phi1 - phi2;
70  if(dPhi_temp<-TMath::Pi()) dPhi_temp += TMath::TwoPi();
71  if(dPhi_temp>=TMath::Pi()) dPhi_temp -= TMath::TwoPi();
72  return dPhi_temp;
73  }
74 
75 //#include "/sphenix/user/stapiaar/jet_correlator/EnergyEnergyCorrelators/eec/include/EEC.hh"
76 
77 using namespace std;
78 
79 MyJetAnalysis::MyJetAnalysis(const std::string& recojetname, const std::string& truthjetname, const std::string& outputfilename)
80  : SubsysReco("MyJetAnalysis_" + recojetname + "_" + truthjetname)
81  , m_recoJetName(recojetname)
82  , m_truthJetName(truthjetname)
83  , m_outputFileName(outputfilename)
84  , m_etaRange(-1.1, 1.1)
85  , m_ptRange(5, 100)
86  , m_trackJetMatchingRadius(.7)
87  , m_hInclusiveE(nullptr)
88  , m_hInclusiveEta(nullptr)
89  , m_hInclusivePhi(nullptr)
90  , m_T(nullptr)
91  , m_event(-1)
92  , _caloevalstackHCALIN(nullptr)
93  , _caloevalstackHCALOUT(nullptr)
94  , _caloevalstackCEMC(nullptr)
95  /* , m_id(-1)
96  , m_nComponent(-1)
97  , m_eta(numeric_limits<float>::signaling_NaN())
98  , m_phi(numeric_limits<float>::signaling_NaN())
99  // , m_e(numeric_limits<vector>::signaling_NaN())
100  , m_pt(numeric_limits<float>::signaling_NaN())
101  // , m_px(numeric_limits<vector>::signaling_NaN())
102  // , m_py(numeric_limits<vector>::signaling_NaN())
103  // , m_pz(numeric_limits<vector>::signaling_NaN())
104  //m_v4(numeric_limits<float>::signaling_NaN())
105  , //m_truthID(-1)
106  , m_truthNComponent(-1)
107  , m_truthEta(numeric_limits<float>::signaling_NaN())
108  , m_truthPhi(numeric_limits<float>::signaling_NaN())
109  // , m_truthE(numeric_limits<vector>::signaling_NaN())
110  , m_truthPt(numeric_limits<float>::signaling_NaN())
111  // , m_truthPx(numeric_limits<vector>::signaling_NaN())
112  // , m_truthPy(numeric_limits<vector>::signaling_NaN())
113  / , m_truthPz(numeric_limits<vecto>::signaling_NaN())
114  */
115  , m_nMatchedTrack(-1)
116 
117 {
118  m_trackdR.fill(numeric_limits<float>::signaling_NaN());
119  m_trackpT.fill(numeric_limits<float>::signaling_NaN());
120  m_trackPID.fill(numeric_limits<float>::signaling_NaN());
121 }
122 
123 MyJetAnalysis::~MyJetAnalysis()
124 {
125 }
126 
127 int MyJetAnalysis::Init(PHCompositeNode* topNode)
128 {
129  if (Verbosity() >= MyJetAnalysis::VERBOSITY_SOME)
130  cout << "MyJetAnalysis::Init - Outoput to " << m_outputFileName << endl;
131 
132  PHTFileServer::get().open(m_outputFileName, "RECREATE");
133 
134  cout << "MyJetAnalysis::Init - Outoput to " << m_outputFileName << endl;
135 
136 
137  // Histograms
138  m_hInclusiveE = new TH1F(
139  "hInclusive_E", //
140  TString(m_recoJetName) + " inclusive jet E;Total jet energy (GeV)", 100, 0, 100);
141 
142  m_hInclusiveEta =
143  new TH1F(
144  "hInclusive_eta", //
145  TString(m_recoJetName) + " inclusive jet #eta;#eta;Jet energy density", 50, -1, 1);
146  m_hInclusivePhi =
147  new TH1F(
148  "hInclusive_phi", //
149  TString(m_recoJetName) + " inclusive jet #phi;#phi;Jet energy density", 50, -M_PI, M_PI);
150 
151  initializeTrees();
152  return Fun4AllReturnCodes::EVENT_OK;
153 }
154 
155 int MyJetAnalysis::End(PHCompositeNode* topNode)
156 {
157  PHTFileServer::get().cd(m_outputFileName);
158  m_T->Write();
159  cout << "MyJetAnalysis::End - Output to " << m_outputFileName << endl;
160  return Fun4AllReturnCodes::EVENT_OK;
161 }
162 
163 int MyJetAnalysis::InitRun(PHCompositeNode* topNode)
164 {
165 return Fun4AllReturnCodes::EVENT_OK;
166 }
167 
168 int MyJetAnalysis::process_event(PHCompositeNode* topNode)
169 {
170  //cout << "MyJetAnalysis::process_event" << endl;
171  if (Verbosity() >= MyJetAnalysis::VERBOSITY_SOME)
172  cout << "MyJetAnalysis::process_event() entered" << endl;
173 
174  ++m_event;
175  if( (m_event % 10) == 0 ) cout << "Event number = "<< m_event << endl;
176 
177  // Making Jets (Needed if not using JetMap)
178  GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
179 
180  // interface to jets (JetMap)
181  //JetMap* jets = findNode::getClass<JetMap>(topNode, m_recoJetName);
182  JetMap* jetsMC = findNode::getClass<JetMap>(topNode, m_truthJetName);
183 
184  // Making Jets (Non JetMap)
185 
186  RawTowerContainer *ihcal_towers_sub = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN_SUB1");
187  RawTowerContainer *ihcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALIN");
188  RawTowerContainer *cemc_towers_sub = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC_RETOWER_SUB1");
189  RawTowerContainer *cemc_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_CEMC");
190  RawTowerContainer *ohcal_towers_sub = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT_SUB1");
191  RawTowerContainer *ohcal_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_HCALOUT");
192 
193  //TowerInfoContainer *ihcal_towers_sub = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
194  //TowerInfoContainer *ihcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
195  //TowerInfoContainer *cemc_towers_sub = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
196  //TowerInfoContainer *cemc_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
197  //TowerInfoContainer *ohcal_towers_sub = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
198  //TowerInfoContainer *ohcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
199 
200  RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
201  RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
202  RawTowerGeomContainer *geomEC = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
203 
204  if (!geomIH || !geomOH ){
205  cout << "Cannot find Tower Geom Nodes" << endl;
206  exit(-1);
207  }
208 
209  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
210  if (!truthinfo){
211  cout << "Error cant find G4 Truth Info" << endl;
212  exit(-1);
213  }
214 
215  CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
216  if (!cent_node){
217  cout << "Error can't find CentralityInfo" << endl;
218  exit(-1);
219  }
220 
221  float cent_mbd = cent_node->get_centile(CentralityInfo::PROP::mbd_NS);
222  m_cent.push_back(cent_mbd);
223 
224  /*
225  std::vector<fastjet::PseudoJet> pseudojet_truth;
226  PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
227  for(PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter){
228  PHG4Particle *part = iter->second;
229  if ((abs(part->get_pid()) >= 12) && (abs(part->get_pid()) <= 16)) continue;
230  if ((part->get_px() == 0.0) && (part->get_py() == 0.0)) continue;
231  float eta = asinh(part->get_pz()/sqrt(pow(part->get_px(),2) + pow(part->get_py(),2)));
232  if (eta < -1.1) continue;
233  if (eta > 1.1) continue;
234  float Part_ID = truthinfo->isEmbeded(part->get_track_id());
235  if (Part_ID != 2) continue;
236  m_truthConstitID.push_back(Part_ID);
237  m_truthConstitE.push_back(part->get_e());
238  //m_truthConstitPt.push_back(part->get_pt());
239  if (Part_ID != 2){
240  cout << "TruthJet has a particle with ID = " << Part_ID << endl;
241  }
242  fastjet::PseudoJet jet_truth(part->get_px(), part->get_py(), part->get_pz(), part->get_e());
243  //jet_truth.set_user_index(Part_ID);
244  pseudojet_truth.push_back(jet_truth);
245  }
246  */
247  GlobalVertex *vtx = vertexmap->begin()->second;
248  float vtxz = NAN;
249  if (vtx) vtxz = vtx->get_z();
250 
251  std::vector<fastjet::PseudoJet> pseudojet_reco;
252  std::vector<fastjet::PseudoJet> pseudojet_reco_truth;
253  CaloRawTowerEval* towerevalHCALIN = new CaloRawTowerEval(topNode, "HCALIN");
254  CaloRawTowerEval* towerevalHCALOUT = new CaloRawTowerEval(topNode, "HCALOUT");
255  CaloRawTowerEval* towerevalCEMC = new CaloRawTowerEval(topNode, "CEMC");
256 
257  RawTowerContainer::ConstRange ICB_begin_end = ihcal_towers->getTowers();
258  RawTowerContainer::ConstIterator ICB_rtiter;
259  for (ICB_rtiter = ICB_begin_end.first; ICB_rtiter != ICB_begin_end.second; ++ICB_rtiter){
260  float Tower_Embedded = 0;
261  RawTower *towerB = ICB_rtiter->second;
262  if (towerB->get_energy() <= 0) continue;
263  PHG4Particle* primary = towerevalHCALIN->max_truth_primary_particle_by_energy(towerB);
264  if (primary){
265  if (truthinfo->isEmbeded(primary->get_track_id()) == 2){
266  Tower_Embedded = 0;
267  }
268  else{
269  Tower_Embedded = 1;
270  }
271  }
272  else Tower_Embedded = -10;
273  towers_id.push_back(towerB->get_id());
274  towers_primary.push_back(Tower_Embedded);
275  }
276 
277  RawTowerContainer::ConstRange IC_begin_end = ihcal_towers_sub->getTowers();
278  RawTowerContainer::ConstIterator IC_rtiter;
279  for (IC_rtiter = IC_begin_end.first; IC_rtiter != IC_begin_end.second; ++IC_rtiter){
280  RawTower *tower = IC_rtiter->second;
281  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
282  assert(tower_geom);
283  if(tower->get_energy() <= 0) continue;
284  float HCALIN_Embedded = -2;
285  //for (ICB_rtiter = ICB_begin_end.first; ICB_rtiter != ICB_begin_end.second; ++ICB_rtiter){
286  //RawTower *towerB = ICB_rtiter->second;
287  //if (towerB->get_id() != tower->get_id()) continue;
288  for (int j = 0; j < (int)towers_id.size(); ++j){
289  if (tower->get_id() != towers_id.at(j)) continue;
290  HCALIN_Embedded = towers_primary.at(j);
291  }
292 
293  double r = tower_geom->get_center_radius();
294  double phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
295  double z0 = tower_geom->get_center_z();
296  double z = z0 - vtxz;
297 
298  double eta = asinh(z/r);
299  double pt = tower->get_energy()/cosh(eta);
300  double px = pt*cos(phi);
301  double py = pt*sin(phi);
302  double pz = pt*sinh(eta);
303 
304  m_IHCAL_Tower_Energy.push_back(tower->get_energy());
305  m_IHCAL_Cent.push_back(cent_mbd);
306  //if (HCALIN_Embedded == -10 || HCALIN_Embedded == 1) continue;
307  fastjet::PseudoJet pseudojet(px, py, pz, tower->get_energy());
308  pseudojet.set_user_index(HCALIN_Embedded);
309  pseudojet_reco.push_back(pseudojet);
310  if (HCALIN_Embedded == -10 || HCALIN_Embedded == 1) continue;
311  pseudojet_reco_truth.push_back(pseudojet);
312  }
313 
314  towers_id.clear();
315  towers_primary.clear();
316 
317  RawTowerContainer::ConstRange EC_begin_end = cemc_towers_sub->getTowers();
318  RawTowerContainer::ConstRange ECB_begin_end = cemc_towers->getTowers();
319  RawTowerContainer::ConstIterator EC_rtiter;
320  RawTowerContainer::ConstIterator ECB_rtiter;
321  for (EC_rtiter = EC_begin_end.first; EC_rtiter != EC_begin_end.second; ++EC_rtiter){
322  RawTower *tower = EC_rtiter->second;
323  RawTowerGeom *tower_geom = geomIH->get_tower_geometry(tower->get_key());
324  assert(tower_geom);
325  if(tower->get_energy() <= 0) continue;
326  int this_ECetabin = geomIH->get_etabin(tower_geom->get_eta());
327  int this_ECphibin = geomIH->get_phibin(tower_geom->get_phi());
328  float HCEMC_Embedded = -2;
329  float Tower_Energy = 0;
330  for (ECB_rtiter = ECB_begin_end.first; ECB_rtiter != ECB_begin_end.second; ++ECB_rtiter){
331  RawTower *towerB = ECB_rtiter->second;
332  RawTowerGeom *towerB_geom = geomEC->get_tower_geometry(towerB->get_key());
333  int this_ECBetabin = geomIH->get_etabin(towerB_geom->get_eta());
334  int this_ECBphibin = geomIH->get_phibin(towerB_geom->get_phi());
335  if (this_ECBetabin != this_ECetabin && this_ECBphibin != this_ECphibin) continue;
336  if(towerB->get_energy() < Tower_Energy) continue;
337  PHG4Particle* primaryEC = towerevalCEMC->max_truth_primary_particle_by_energy(towerB);
338  if (primaryEC){
339  if (truthinfo->isEmbeded(primaryEC->get_track_id()) == 2){
340  HCEMC_Embedded = 2;
341  }
342  else{
343  HCEMC_Embedded = 3;
344  }
345  }
346  else{
347  HCEMC_Embedded = -11;
348  }
349  Tower_Energy = towerB->get_energy();
350  }
351 
352 
353  double r = tower_geom->get_center_radius();
354  double phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
355  double z0 = tower_geom->get_center_z();
356  double z = z0 - vtxz;
357 
358  double eta = asinh(z/r);
359  double pt = tower->get_energy()/cosh(eta);
360  double px = pt*cos(phi);
361  double py = pt*sin(phi);
362  double pz = pt*sinh(eta);
363 
364  m_EMCAL_Tower_Energy.push_back(tower->get_energy());
365  m_EMCAL_Cent.push_back(cent_mbd);
366  //if(HCEMC_Embedded == -11 || HCEMC_Embedded == 3) continue;
367  fastjet::PseudoJet pseudojet(px, py, pz, tower->get_energy());
368  pseudojet.set_user_index(HCEMC_Embedded);
369  pseudojet_reco.push_back(pseudojet);
370  if(HCEMC_Embedded == -11 || HCEMC_Embedded == 3) continue;
371  pseudojet_reco_truth.push_back(pseudojet);
372  }
373 
374  RawTowerContainer::ConstRange OCB_begin_end = ohcal_towers->getTowers();
375  RawTowerContainer::ConstIterator OCB_rtiter;
376  for (OCB_rtiter = OCB_begin_end.first; OCB_rtiter != OCB_begin_end.second; ++OCB_rtiter){
377  float Tower_Embedded = 0;
378  RawTower *towerB = OCB_rtiter->second;
379  if (towerB->get_energy() <= 0) continue;
380  PHG4Particle* primary = towerevalHCALOUT->max_truth_primary_particle_by_energy(towerB);
381  if (primary){
382  if (truthinfo->isEmbeded(primary->get_track_id()) == 2){
383  Tower_Embedded = 4;
384  }
385  else{
386  Tower_Embedded = 5;
387  }
388  }
389  else Tower_Embedded = -12;
390  towers_id.push_back(towerB->get_id());
391  towers_primary.push_back(Tower_Embedded);
392  }
393 
394  RawTowerContainer::ConstRange OC_begin_end = ohcal_towers_sub->getTowers();
395  RawTowerContainer::ConstIterator OC_rtiter;
396  for (OC_rtiter = OC_begin_end.first; OC_rtiter != OC_begin_end.second; ++OC_rtiter){
397  RawTower *tower = OC_rtiter->second;
398  RawTowerGeom *tower_geom = geomOH->get_tower_geometry(tower->get_key());
399  assert(tower_geom);
400  if(tower->get_energy() <= 0) continue;
401  float HCALOUT_Embedded = -1;
402  for (int j = 0; j < (int)towers_id.size(); ++j){
403  if (tower->get_id() != towers_id.at(j)) continue;
404  HCALOUT_Embedded = towers_primary.at(j);
405  }
406  double r = tower_geom->get_center_radius();
407  double phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
408  double z0 = tower_geom->get_center_z();
409  double z = z0 - vtxz;
410 
411  double eta = asinh(z/r);
412  double pt = tower->get_energy()/cosh(eta);
413  double px = pt*cos(phi);
414  double py = pt*sin(phi);
415  double pz = pt*sinh(eta);
416 
417  m_OHCAL_Tower_Energy.push_back(tower->get_energy());
418  m_OHCAL_Cent.push_back(cent_mbd);
419  //if (HCALOUT_Embedded == -12 || HCALOUT_Embedded == 5) continue;
420  fastjet::PseudoJet pseudojet(px, py, pz, tower->get_energy());
421  pseudojet.set_user_index(HCALOUT_Embedded);
422  pseudojet_reco.push_back(pseudojet);
423  if (HCALOUT_Embedded == -12 || HCALOUT_Embedded == 5) continue;
424  pseudojet_reco_truth.push_back(pseudojet);
425  }
426 
427  towers_id.clear();
428  towers_primary.clear();
429 
430  fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, 0.2, fastjet::E_scheme,fastjet::Best);
431  fastjet::ClusterSequence jetFinder_reco_truth(pseudojet_reco_truth, jetDef);
432  fastjet::ClusterSequence jetFinder_reco(pseudojet_reco, jetDef);
433  std::vector<fastjet::PseudoJet> fastjets_reco_truth = jetFinder_reco_truth.inclusive_jets();
434  std::vector<fastjet::PseudoJet> fastjets_reco = jetFinder_reco.inclusive_jets();
435 
436  pseudojet_reco_truth.clear();
437  pseudojet_reco.clear();
438 
439  //Saves Truth Jet Info
440  for (JetMap::Iter iter = jetsMC->begin(); iter != jetsMC->end(); ++iter){
441  //for (int j = 0; j < (int)fastjets.size(); ++j){
442  Jet* truthjet = iter->second;
443  assert(truthjet);
444  bool eta_cut = (truthjet->get_eta() >= m_etaRange.first) and (truthjet->get_eta() <= m_etaRange.second);
445  if (not eta_cut) continue;
446  //cout << "Jet Truth = " << j << endl;
447  /*
448  vector<fastjet::PseudoJet> consts = jetFinder.constituents(fastjets[j]);
449  //cout << "Size of Truth Jet" << (int)consts.size() << endl;
450  for (int i = 0; i < (int)consts.size(); i++){
451  if (consts[i].user_index() != -10){
452  //cout << consts[i].user_index() << endl;
453  }
454  }
455 
456  for(Jet::ConstIter comp = truthjet->begin_comp(); comp != truthjet->end_comp(); ++comp){
457  const int this_embed_id = truthinfo->isEmbeded((*comp).second);
458  if (this_embed_id != 2){
459  cout << "Truth jet has non embedded particle with id = " << this_embed_id << endl;
460  }
461  }
462  */
463  m_truthEta.push_back(truthjet->get_eta());
464  m_truthPhi.push_back(truthjet->get_phi());
465  m_truthE.push_back(truthjet->get_e());
466  m_truthPt.push_back(truthjet->get_pt());
467  m_truthPx.push_back(truthjet->get_px());
468  m_truthPy.push_back(truthjet->get_py());
469  m_truthPz.push_back(truthjet->get_pz());
470  m_truthNComponent.push_back(truthjet->size_comp());
471 
472  }
473  //Saving Reco Truth Jets
474  for (int j = 0; j < (int)fastjets_reco_truth.size(); ++j){
475 
476  vector<fastjet::PseudoJet> consts_reco_truth = jetFinder_reco_truth.constituents(fastjets_reco_truth[j]);
477 
478  m_recotruthnComponent.push_back((int)consts_reco_truth.size());
479  m_recotruthEta.push_back(fastjets_reco_truth[j].eta());
480  m_recotruthPhi.push_back(fastjets_reco_truth[j].phi_std());
481  m_recotruthPt.push_back(fastjets_reco_truth[j].perp());
482  m_recotruthPx.push_back(fastjets_reco_truth[j].px());
483  m_recotruthPy.push_back(fastjets_reco_truth[j].py());
484  m_recotruthPz.push_back(fastjets_reco_truth[j].pz());
485  m_recotruthE.push_back(fastjets_reco_truth[j].e());
486 
487  }
488  //saving closest recojets
489  //for (JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter)
490  for (int j = 0; j < (int)fastjets_reco.size(); ++j)
491  {
492  vector<fastjet::PseudoJet> consts = jetFinder_reco.constituents(fastjets_reco[j]);
493 
494  //Jet* jet = iter->second;
495 
496  //Needed For Making Jets
497  /*
498  float nEmbedded = 0;
499  float nTotal = 0;
500  float nUnidentified_IHCAL = 0;
501  float nUnidentified_OHCAL = 0;
502  float nIHCAL_EM = 0;
503  float nIHCAL_BC = 0;
504  float nCEMC_EM = 0;
505  float nCEMC_BC = 0;
506  float nOHCAL_EM = 0;
507  float nOHCAL_BC = 0;
508  */
509  //ID ordered by ID = 0 embedded Inner, ID = 1 nonembedded Inner, continues like this for EMCal then Outer
510  for (int i = 0; i < (int)consts.size(); i++){
511  m_CAL_ID.push_back(consts[i].user_index());
512  m_Constit_E.push_back(consts[i].e());
513  m_Constit_Cent.push_back(cent_mbd);
514  /*
515  if (consts[i].user_index() == 0 || consts[i].user_index() == 2 || consts[i].user_index() == 4){
516  nEmbedded = nEmbedded + 1;
517  if (consts[i].user_index() == 0){
518  nIHCAL_EM = nIHCAL_EM + 1;
519  }
520  if (consts[i].user_index() == 2){
521  nCEMC_EM = nCEMC_EM + 1;
522  }
523  if (consts[i].user_index() == 4){
524  nOHCAL_EM = nOHCAL_EM + 1;
525  }
526  }
527  else{
528  if (consts[i].user_index() == 1){
529  nIHCAL_BC = nIHCAL_BC + 1;
530  }
531  if (consts[i].user_index() == 3){
532  nCEMC_BC = nCEMC_BC + 1;
533  }
534  if (consts[i].user_index() == 5){
535  nOHCAL_BC = nOHCAL_BC + 1;
536  }
537  }
538 
539  if (consts[i].user_index() == -11){
540  nUnidentified_IHCAL = nUnidentified_IHCAL + 1;
541  }
542  if (consts[i].user_index() == -10){
543  nUnidentified_OHCAL = nUnidentified_OHCAL + 1;
544  }
545  nTotal = nTotal + 1;
546  */
547  }
548  /*
549  m_Embedded_IHCAL.push_back(nIHCAL_EM);
550  m_Embedded_CEMC.push_back(nCEMC_EM);
551  m_Embedded_OHCAL.push_back(nOHCAL_EM);
552  m_Background_IHCAL.push_back(nIHCAL_BC);
553  m_Background_CEMC.push_back(nCEMC_BC);
554  m_Background_OHCAL.push_back(nOHCAL_BC);
555  m_Embedded_Count.push_back(nEmbedded);
556  m_Unidentified_IHCAL.push_back(nUnidentified_IHCAL);
557  m_Unidentified_OHCAL.push_back(nUnidentified_OHCAL);
558  m_Total_Count.push_back(nTotal);
559  */
560  // Save reco jet information (In code jet reconstruction)
561  m_eta.push_back(fastjets_reco[j].eta());
562  m_phi.push_back(fastjets_reco[j].phi_std());
563  m_pt.push_back(fastjets_reco[j].perp());
564  m_px.push_back(fastjets_reco[j].px());
565  m_py.push_back(fastjets_reco[j].py());
566  m_pz.push_back(fastjets_reco[j].pz());
567  m_e.push_back(fastjets_reco[j].e());
568 
569  //Save the reco jet Information (Jet Map)
570 
571  //m_id.push_back(jet->get_id());
572  //m_nComponent.push_back(jet->size_comp());
573  //m_eta.push_back(jet->get_eta());
574  //m_phi.push_back(jet->get_phi());
575  //m_pt.push_back(jet->get_pt());
576  //m_px.push_back(jet->get_px());
577  //m_py.push_back(jet->get_py());
578  //m_pz.push_back(jet->get_pz());
579  //m_e.push_back(jet->get_e());
580 
581  float m_Loop_dR = 100000;// dR used for Loop to save closer truth jets
582 
583  //for (int i = 0; i < (int)fastjets.size(); ++i){
584  for(std::size_t i = 0, max = (int)m_truthPhi.size(); i != max; ++i){
585  //if(m_truthE.at(i) < 0) continue;
586 
587  double dEta_temp = fastjets_reco[j].eta() - m_truthEta.at(i);
588  double dPhi_temp = DeltaPhi(fastjets_reco[j].phi_std(), m_truthPhi.at(i));
589 
590  //double dEta_temp = jet->get_eta() - m_truthEta.at(i);
591  //double dPhi_temp = DeltaPhi(jet->get_phi(), m_truthPhi.at(i));
592 
593  float dR_temp = sqrt(dEta_temp * dEta_temp + dPhi_temp * dPhi_temp);
594  if(dR_temp > m_Loop_dR) continue;
595 
596  // Saving truth R = 0.2 matched Info
597  temp_TE_Matched = m_truthE.at(i);
598  temp_TPhi_Matched = m_truthPhi.at(i);
599  temp_TEta_Matched = m_truthEta.at(i);
600 
601 
602  m_Loop_dR = dR_temp;//Resets cutoff dR to mae sure to save only closer truth jets
603  }
604 
605  m_E_Matched.push_back(fastjets_reco[j].e());
606  m_Phi_Matched.push_back(fastjets_reco[j].phi_std());
607  m_Eta_Matched.push_back(fastjets_reco[j].eta());
608 
609  m_TE_Matched.push_back(temp_TE_Matched);
610  m_TPhi_Matched.push_back(temp_TPhi_Matched);
611  m_TEta_Matched.push_back(temp_TEta_Matched);
612  m_dR.push_back(m_Loop_dR);
613  }
614 
615  m_T->Fill();
616  m_IHCAL_Tower_Energy.clear();
617  m_IHCAL_Cent.clear();
618  m_EMCAL_Tower_Energy.clear();
619  m_EMCAL_Cent.clear();
620  m_OHCAL_Tower_Energy.clear();
621  m_OHCAL_Cent.clear();
622  //clearing truth info
623  //m_truthID.clear();
624  //m_truthNComponent.clear();
625  m_truthEta.clear();
626  m_truthPhi.clear();
627  m_truthE.clear();
628  m_truthPt.clear();
629  m_truthPx.clear();
630  m_truthPy.clear();
631  m_truthPz.clear();
632  //m_truthConstitE.clear();
633  //m_truthConstitID.clear();
634  //m_truthConstitPt.clear();
635  m_cent.clear();
636 
637  //clearing reco info
638  //m_id.clear();
639  m_nComponent.clear();
640  m_eta.clear();
641  m_phi.clear();
642  m_pt.clear();
643  m_e.clear();
644  m_px.clear();
645  m_py.clear();
646  m_pz.clear();
647 
648  m_recotruthnComponent.clear();
649  m_recotruthEta.clear();
650  m_recotruthPhi.clear();
651  m_recotruthPt.clear();
652  m_recotruthE.clear();
653  m_recotruthPx.clear();
654  m_recotruthPy.clear();
655  m_recotruthPz.clear();
656 
657  //MAtched infor clear
658  m_TE_Matched.clear();
659  m_TPhi_Matched.clear();
660  m_TEta_Matched.clear();
661 
662  m_E_Matched.clear();
663  m_Phi_Matched.clear();
664  m_Eta_Matched.clear();
665 
666  m_dR.clear();
667  m_cent.clear();
668  m_CAL_ID.clear();
669  m_Constit_E.clear();
670  m_Constit_Cent.clear();
671  /*
672  m_Embedded_IHCAL.clear();
673  m_Embedded_CEMC.clear();
674  m_Embedded_OHCAL.clear();
675  m_Background_IHCAL.clear();
676  m_Background_CEMC.clear();
677  m_Background_OHCAL.clear();
678  m_Embedded_Count.clear();
679  m_Total_Count.clear();
680  m_Unidentified_IHCAL.clear();
681  m_Unidentified_OHCAL.clear();
682  */
683  Clean();
684 
685  return Fun4AllReturnCodes::EVENT_OK;
686 }
687 
688 void MyJetAnalysis::initializeTrees()
689 {
690 
691 m_T = new TTree("T", "MyJetAnalysis Tree");
692 
693 //Storing All Reco Info
694 m_T->Branch("m_event", &m_event, "event/I");
695 m_T->Branch("id", &m_id);
696 m_T->Branch("nComponent", &m_nComponent);
697 m_T->Branch("eta", &m_eta);
698 m_T->Branch("phi", &m_phi);
699 m_T->Branch("e", &m_e);
700 m_T->Branch("pt", &m_pt);
701 m_T->Branch("px", &m_px);
702 m_T->Branch("py", &m_py);
703 m_T->Branch("pz", &m_pz);
704 
705 m_T->Branch("IHCAL_Tower_Energy",&m_IHCAL_Tower_Energy);
706 m_T->Branch("IHCAL_Cent",&m_IHCAL_Cent);
707 m_T->Branch("EMCAL_Tower_Energy",&m_EMCAL_Tower_Energy);
708 m_T->Branch("EMCAL_Cent",&m_EMCAL_Cent);
709 m_T->Branch("OHCAL_Tower_Energy",&m_OHCAL_Tower_Energy);
710 m_T->Branch("OHCAL_Cent",&m_OHCAL_Cent);
711 
712 
713 //Storing All Reco Truth Info
714 m_T->Branch("recotruthnComponent", &m_recotruthnComponent);
715 m_T->Branch("recotruthEta", &m_recotruthEta);
716 m_T->Branch("recotruthPhi", &m_recotruthPhi);
717 m_T->Branch("recotruthPt", &m_recotruthPt);
718 m_T->Branch("recotruthPx", &m_recotruthPx);
719 m_T->Branch("recotruthPy", &m_recotruthPy);
720 m_T->Branch("recotruthPz", &m_recotruthPz);
721 m_T->Branch("recotruthE", &m_recotruthE);
722 
723 //Storing Matched Jets Info
724 m_T->Branch("dR", &m_dR);
725 
726 //Storing Reco Jet Info
727 m_T->Branch("cent", &m_cent);
728 m_T->Branch("CAL_ID", &m_CAL_ID);
729 //m_T->Branch("CAL_ID", &m_CAL_ID_RT);
730 m_T->Branch("Constit_Energy", &m_Constit_E);
731 //m_T->Branch("Constit_Energy_RT", &m_Constit_RT);
732 m_T->Branch("Constit_Cent", &m_Constit_Cent);
733 //m_T->Branch("Constit_Cent_RT", &m_Constit_Cent_RT);
734 //m_T->Branch("Embedded_Count", &m_Embedded_Count);
735 //m_T->Branch("Embedded_nIHCAL", &m_Embedded_IHCAL);
736 //m_T->Branch("Embedded_nCEMC", &m_Embedded_CEMC);
737 //m_T->Branch("Embedded_nOHCAL", &m_Embedded_OHCAL);
738 //m_T->Branch("Background_nIHCAL", &m_Background_IHCAL);
739 //m_T->Branch("Background_nCEMC", &m_Background_CEMC);
740 //m_T->Branch("Background_nOHCAL", &m_Background_OHCAL);
741 //m_T->Branch("Unidentified_IHCAL", &m_Unidentified_IHCAL);
742 //m_T->Branch("Unidentified_OHCAL", &m_Unidentified_OHCAL);
743 //m_T->Branch("Total_Count", &m_Total_Count);
744 //m_T->Branch("truthID", &m_truthID);
745 //m_T->Branch("truthNComponent", &m_truthNComponent);
746 //m_T->Branch("truthConstitE", &m_truthConstitE);
747 //m_T->Branch("truthConstitID", &m_truthConstitID);
748 m_T->Branch("truthEta", &m_truthEta);
749 m_T->Branch("truthPhi", &m_truthPhi);
750 m_T->Branch("truthE", &m_truthE);
751 m_T->Branch("truthPt", &m_truthPt);
752 m_T->Branch("truthPx", &m_truthPx);
753 m_T->Branch("truthPy", &m_truthPy);
754 m_T->Branch("truthPz", &m_truthPz);
755 //m_T->Branch("truthdR", &m_truthdR, "truthdR/F");
756 
757 }
758 
759  void MyJetAnalysis::Clean()
760 {
761 //m_cent.clear();
762 }
763 
764 
765 
766 
767 
768 
769