Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JetRhoMedian.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file JetRhoMedian.cc
1 #include "JetRhoMedian.h"
2 #include "fastjet/AreaDefinition.hh"
3 #include "fastjet/ClusterSequenceArea.hh"
4 #include "fastjet/Selector.hh"
5 #include "fastjet/tools/BackgroundEstimatorBase.hh"
6 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
7 
8 
9 #include <TFile.h>
10 #include <TH1F.h>
11 #include <TH2F.h>
12 #include <TRandom3.h>
13 #include <TString.h>
14 #include <TTree.h>
15 #include <TH2D.h>
16 #include <TVector3.h>
17 #include <algorithm>
18 #include <cassert>
19 #include <centrality/CentralityInfo.h>
20 #include <cmath>
21 #include <fastjet/ClusterSequence.hh>
22 #include <fastjet/JetDefinition.hh>
23 #include <fastjet/PseudoJet.hh>
24 #include <fun4all/Fun4AllReturnCodes.h>
25 #include <fun4all/Fun4AllServer.h>
26 #include <fun4all/PHTFileServer.h>
27 #include <g4eval/JetEvalStack.h>
28 #include <jetbase/JetInput.h>
29 #include <jetbase/JetMapv1.h>
30 #include <iostream>
31 #include <limits>
32 #include <phool/PHCompositeNode.h>
33 #include <phool/getClass.h>
34 #include <stdexcept>
35 #include <trackbase_historic/SvtxTrackMap.h>
36 #include <globalvertex/GlobalVertexMap.h>
37 
38 // using namespace fastjet;
39 
40 JetRhoMedian::JetRhoMedian
41  ( const std::string& _outputfilename
42  , const float _jet_R
43  , const std::string& truthjetname
44  , const std::string& sub1jetname
45  , const float min_lead_truth_pt
46  )
47  : SubsysReco("JetRhoMedian")
48  , m_outputfilename { _outputfilename }
49  , jet_R { _jet_R }
50  , m_truthJetName { truthjetname }
51  , m_sub1JetName { sub1jetname }
52  , m_min_lead_truth_pt { min_lead_truth_pt }
53  , m_T { nullptr }
54  , m_inputs {}
55 { }
56 
57 JetRhoMedian::~JetRhoMedian()
58 {
59  for (unsigned int i = 0; i < m_inputs.size(); ++i) delete m_inputs[i];
60  m_inputs.clear();
61 }
62 
63 int JetRhoMedian::Init(PHCompositeNode* topNode)
64 {
65  if (Verbosity() >= JetRhoMedian::VERBOSITY_SOME)
66  std::cout << "JetRhoMedian::Init - Outoput to " << m_outputfilename << std::endl;
67 
68  std::cout << " File out: " << m_outputfilename << std::endl;
69  PHTFileServer::get().open(m_outputfilename, "RECREATE");
70 
71  //Tree
72  m_T = new TTree("T", "JetRhoMedian Tree");
73 
74  // int m_event;
75  /* m_T->Branch("id", &m_id); */
76  m_T->Branch("rho", &m_rho);
77  m_T->Branch("rho_sigma", &m_rho_sigma);
78  m_T->Branch("cent", &m_cent);
79  m_T->Branch("cent_mdb", &m_cent_mdb);
80  m_T->Branch("cent_epd", &m_cent_epd);
81  m_T->Branch("impactparam", &m_impactparam);
82 
83 
84  //Truth Jets
85  m_T->Branch("TruthJetEta", &m_TruthJetEta);
86  m_T->Branch("TruthJetPhi", &m_TruthJetPhi);
87  m_T->Branch("TruthJetPt", &m_TruthJetPt);
88 
89  //Sub1 Jets
90  m_T->Branch("sub1JetEta", &m_Sub1JetEta);
91  m_T->Branch("sub1JetPhi", &m_Sub1JetPhi);
92  m_T->Branch("sub1JetPt", &m_Sub1JetPt);
93 
94  m_T->Branch("rhoAJetEta", &m_CaloJetEta);
95  m_T->Branch("rhoAJetPhi", &m_CaloJetPhi);
96  m_T->Branch("rhoAJetPt", &m_CaloJetPt);
97  m_T->Branch("rhoAJetPtLessRhoA", &m_CaloJetPtLessRhoA);
98  m_T->Branch("rhoAJetArea", &m_CaloJetArea);
99 
100  return Fun4AllReturnCodes::EVENT_OK;
101 }
102 
103 int JetRhoMedian::End(PHCompositeNode* topNode)
104 {
105  std::cout << "JetRhoMedian::End - Output to " << m_outputfilename << std::endl;
106  PHTFileServer::get().cd(m_outputfilename);
107 
108  m_T->Write();
109  return Fun4AllReturnCodes::EVENT_OK;
110 }
111 
112 int JetRhoMedian::InitRun(PHCompositeNode* topNode)
113 {
114  topNode->print();
115  std::cout << " Input Selections:" << std::endl;
116  for (unsigned int i = 0; i < m_inputs.size(); ++i) m_inputs[i]->identify();
117  return Fun4AllReturnCodes::EVENT_OK;
118 }
119 
120 int JetRhoMedian::process_event(PHCompositeNode* topNode)
121 {
122  clear_vectors();
123  ++nevent;
124 
125  // cut for only event with good vertices
126  GlobalVertexMap *mapVtx = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
127  if (mapVtx->size() == 0) return Fun4AllReturnCodes::EVENT_OK;
128 
129  // centrality
130  CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
131  if (cent_node)
132  {
133  m_cent_epd = cent_node->get_centile(CentralityInfo::PROP::epd_NS);
134  m_cent_mdb = cent_node->get_centile(CentralityInfo::PROP::mbd_NS);
135  m_impactparam = cent_node->get_quantity(CentralityInfo::PROP::bimp);
136  m_cent = cent_node->get_centile(CentralityInfo::PROP::bimp);
137  }
138 
139  const double ghost_max_rap = 1.1;
140  const double ghost_R = 0.01;
141  const double max_jet_rap = 1.1; // - jet_R;
142 
143 
144  /* fastjet::JetDefinition jet_def(fastjet::kt_algorithm, 0.4); // JET DEFINITION */
145  /* const double ghost_max_rap { 1.0 }; */
146 
147  /* fastjet::AreaDefinition area_def_bkgd( */
148  /* fastjet::active_area_explicit_ghosts, fastjet::GhostedAreaSpec(ghost_max_rap, 1, ghost_R)); */
149  /* fastjet::Selector selector_rm2 = jetrap * (!fastjet::SelectorNHardest(2)); // <-- */
150 
151  //interface to truth jets
152  JetMap* jetsMC = findNode::getClass<JetMap>(topNode, m_truthJetName);
153  if (!jetsMC )
154  {
155  std::cout
156  << "MyJetAnalysis::process_event - Error can not find DST Truth JetMap node "
157  << m_truthJetName << std::endl;
158  exit(-1);
159  }
160 
161  bool isfirst = true;
162  int fixmek = 0;
163  auto vjets = jetsMC->vec();
164  std::sort(vjets.begin(), vjets.end(), [](Jet* a, Jet* b){return a->get_pt() > b->get_pt();});
165  if (Verbosity() > 100) {
166  if (m_cent_mdb<10) {
167  std::cout << " NEW EVENT[" << (outevent++) << "]" << std::endl;
168  std::cout << " --- Truth Jets [set(" << outevent << ")]" << std::endl;
169  }
170  }
171  for (auto& jet : vjets) {
172  float pt = jet->get_pt();
173  if (pt < 5.) continue;
174  float eta = jet->get_eta();
175  if (eta < -max_jet_rap || eta > max_jet_rap) continue;
176  if (Verbosity()>100 && m_cent_mdb<10) {
177  std::cout << Form("k[%i] pt:phi:eta[%5.2f:%5.2f:%5.2f]", fixmek++, jet->get_pt(), jet->get_phi(), jet->get_eta()) << std::endl;
178  }
179 
180  if (isfirst) {
181  isfirst = false;
182  if (pt < m_min_lead_truth_pt) return Fun4AllReturnCodes::EVENT_OK;
183  }
184  m_TruthJetPt .push_back(jet->get_pt());
185  m_TruthJetEta.push_back(jet->get_eta());
186  m_TruthJetPhi.push_back(jet->get_phi());
187  }
188 
189  //get the SUB1 jets
190  if (Verbosity() > 100 && m_cent_mdb<10) {
191  std::cout << " --- Sub1 Jets [set(" << outevent << ")]" << std::endl;
192  fixmek = 0;
193  }
194  JetMap* jets_sub1 = findNode::getClass<JetMap>(topNode, m_sub1JetName);
195  if (jets_sub1 != nullptr) {
196  vjets = jets_sub1->vec();
197  std::sort(vjets.begin(), vjets.end(), [](Jet* a, Jet* b){return a->get_pt() > b->get_pt();});
198  fixmek = 0;
199  if (Verbosity()>100 && m_cent_mdb<10) std::cout << "SUB1 jets" << std::endl;
200  for (auto& jet : jets_sub1->vec()) {
201  // will return jets in order of descending pT
202  float pt = jet->get_pt();
203  float eta = jet->get_eta();
204  if (eta < -max_jet_rap || eta > max_jet_rap) continue;
205  if (pt < 5.) continue;
206  if (Verbosity()>100 && m_cent_mdb<10) std::cout << Form("k[%i] pt:phi:eta[%5.2f:%5.2f:%5.2f]", fixmek++, jet->get_pt(), jet->get_phi(), jet->get_eta()) << std::endl;
207  m_Sub1JetPt .push_back(jet->get_pt());
208  m_Sub1JetEta.push_back(jet->get_eta());
209  m_Sub1JetPhi.push_back(jet->get_phi());
210  }
211  }
212 
213  const int TOW_PRINT_INT = 8000;
214 
215  int n_inputs {0};
216  std::vector<fastjet::PseudoJet> calo_pseudojets;
217  for (auto inp : m_inputs) {
218  if (Verbosity()>TOW_PRINT_INT) std::cout << " Input Tower Set("<<n_inputs++<<"): ";
219  float sum_pt {0.};
220  float n_towers {0};
221  std::vector<Jet *> particles = inp->get_input(topNode);
222  for (unsigned int i = 0; i < particles.size(); ++i)
223  {
224  auto& part = particles[i];
225  float this_e = part->get_e();
226  if (this_e == 0.) continue;
227  float this_px = particles[i]->get_px();
228  float this_py = particles[i]->get_py();
229  float this_pz = particles[i]->get_pz();
230  if (this_e < 0)
231  {
232  // make energy = +1 MeV for purposes of clustering
233  float e_ratio = 0.001 / this_e;
234  this_e = this_e * e_ratio;
235  this_px = this_px * e_ratio;
236  this_py = this_py * e_ratio;
237  this_pz = this_pz * e_ratio;
238  }
239  fastjet::PseudoJet pseudojet(this_px, this_py, this_pz, this_e);
240  /* if (pseudojet.pt() < 0.2) continue; */
241  calo_pseudojets .push_back(pseudojet);
242  if (Verbosity() > TOW_PRINT_INT) {
243  sum_pt += pseudojet.perp();
244  ++n_towers;
245  }
246  }
247  if (Verbosity()>TOW_PRINT_INT) std::cout << " sumpt(" << sum_pt <<") from ntowers("<<n_towers <<")" << std::endl;
248  for (auto &p : particles) delete p;
249  }
250 
251 
252  fastjet::Selector jetrap = fastjet::SelectorAbsEtaMax(max_jet_rap);
253  fastjet::Selector not_pure_ghost = !fastjet::SelectorIsPureGhost();
254 
255  fastjet::Selector selection = jetrap && not_pure_ghost;
256  fastjet::AreaDefinition area_def(
257  fastjet::active_area_explicit_ghosts, fastjet::GhostedAreaSpec(ghost_max_rap, 1, ghost_R));
258 
259  fastjet::JetDefinition jet_def_antikt(fastjet::antikt_algorithm, jet_R);
260 
261 
262  fastjet::JetDefinition jet_def_bkgd(fastjet::kt_algorithm, jet_R); // <--
263  fastjet::Selector selector_rm2 = jetrap * (!fastjet::SelectorNHardest(2)); // <--
264  //
265  fastjet::JetMedianBackgroundEstimator bge_rm2 {selector_rm2, jet_def_bkgd, area_def};
266  bge_rm2.set_particles(calo_pseudojets);
267 
268  m_rho = bge_rm2.rho();
269  m_rho_sigma = bge_rm2.sigma();
270 
271  mBGE_mean_area = bge_rm2 .mean_area();
272  mBGE_empty_area = bge_rm2 .empty_area();
273  mBGE_n_empty_jets = bge_rm2 .n_empty_jets();
274  mBGE_n_jets_used = bge_rm2 .n_jets_used();
275 
276 
277  fastjet::ClusterSequenceArea clustSeq(calo_pseudojets, jet_def_antikt, area_def);
278  std::vector<fastjet::PseudoJet> jets
279  = sorted_by_pt( selection( clustSeq.inclusive_jets(5.) ));
280 
281  if (Verbosity()>100 && m_cent_mdb<10) {
282  fixmek = 0;
283  std::cout << " --- rhoA jets [set("<<outevent<<")]" << std::endl;
284  for (auto jet : jets) {
285  if (jet.perp()-m_rho*jet.area() > 5.) {
286  auto phi = jet.phi();
287  while (phi > M_PI) phi -= 2*M_PI;
288  std::cout << Form("k[%i] pt:phi:eta[%5.2f:%5.2f:%5.2f]", fixmek++, (jet.perp()-m_rho*jet.area()), phi, jet.eta()) << std::endl;
289  }
290  }
291  }
292 
293  for (auto jet : jets) {
294  if (false && Verbosity() >= JetRhoMedian::VERBOSITY_SOME) {
295  auto cst = fastjet::sorted_by_pt(jet.constituents());
296  int i =0;
297  for (auto p : cst) {
298  if (p.is_pure_ghost()) continue;
299  std::cout << Form(" %3i|%4.2f",i++,p.perp());
300  if (i%7==0) std::cout << std::endl;
301  }
302  std::cout << std::endl;
303  }
304 
305  m_CaloJetEta .push_back( jet.eta() );
306  m_CaloJetPhi .push_back( jet.phi_std() );
307  m_CaloJetPt .push_back( jet.pt() );
308  m_CaloJetPtLessRhoA .push_back( jet.pt() - m_rho * jet.area());
309  m_CaloJetArea .push_back( jet.area());
310  }
311 
312 
313  m_T->Fill();
314  return Fun4AllReturnCodes::EVENT_OK;
315 }
316 
317 void JetRhoMedian::clear_vectors() {
318  /* if (nevent % 1 == 0) cout << " z0 " << m_CaloJetArea.size(); */
319  m_CaloJetEta .clear();
320  m_CaloJetPhi .clear();
321  m_CaloJetPt .clear();
322  m_CaloJetPtLessRhoA .clear();
323  m_CaloJetArea .clear();
324 
325  /* if (nevent % 1 == 0) cout << " " << m_TruthJetPt.size(); */
326  m_TruthJetEta .clear();
327  m_TruthJetPhi .clear();
328  m_TruthJetPt .clear();
329 
330  m_Sub1JetEta .clear();
331  m_Sub1JetPhi .clear();
332  m_Sub1JetPt .clear();
333 
334 }