Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fill_jet_tree.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file fill_jet_tree.C
1 #include "fastjet/ClusterSequence.hh"
2 
3 #include <vector>
4 #include <iostream>
5 #include <fstream>
6 #include <cstdio>
7 #include <cstdlib>
8 #include <cassert>
9 
10 #include "TTree.h"
11 #include "TFile.h"
12 
13 #include "eicsmear/erhic/EventBase.h"
14 #include "eicsmear/erhic/EventPythia.h"
15 #include "eicsmear/erhic/Particle.h"
16 #include "eicsmear/erhic/ParticleMC.h"
17 #include "eicsmear/smear/EventSmear.h"
18 #include "eicsmear/smear/ParticleMCS.h"
19 
20 using namespace fastjet;
21 using namespace std;
22 
23 
26 PseudoJet* find_matching_jet( const PseudoJet* refjet, vector<PseudoJet>* vjets )
27 {
28  /* PidCandidate with eta, phi closest to reference */
29  PseudoJet* best_candidate = NULL;
30 
31  float eta_ref = refjet->eta();
32  float phi_ref = refjet->phi();
33 
34  /* find jet with smallest delta_R */
35  float min_delta_R = 100;
36  vector<PseudoJet>::iterator min_delta_R_iter = vjets->end();
37 
38  vector<PseudoJet>::iterator ijet;
39 
40  for ( ijet = vjets->begin(); ijet != vjets->end(); ijet++ )
41  {
42  float eta = ijet->eta();
43  float phi = ijet->phi();
44 
45  float delta_R = sqrt( pow( eta - eta_ref, 2 ) + pow( phi - phi_ref, 2 ) );
46 
47  if ( delta_R < min_delta_R )
48  {
49  min_delta_R_iter = ijet; ;
50  min_delta_R = delta_R;
51  }
52  }
53 
54  /* set best_candidate to PidCandiate with smallest delta_R within reasonable range*/
55  if ( min_delta_R_iter != vjets->end() && min_delta_R < 0.5 )
56  best_candidate = &(*min_delta_R_iter);
57 
58  return best_candidate;
59 }
60 
61 
64 int main(int argc, char* argv[]) {
65 
66  if ( argc != 6 )
67  {
68  cerr << "Wrong number of arguments! Exit." << endl;
69  return -1;
70  }
71 
72  const char* truthFileName = argv[1];
73  const char* smearFileName = argv[2];
74  const char* outFileName = argv[3];
75 
76  /* R value for fastjet algorithm */
77  const float fastjetR = atof(argv[4]);
78 
79  /* minimum pt for jets that are kept for analysis */
80  const float ptmin = atof(argv[5]);
81 
82  cout << "truthFileName = " << truthFileName << endl;
83  cout << "smearFileName = " << smearFileName << endl;
84  cout << "outFileName = " << outFileName << endl;
85  cout << "fastjetR = " << fastjetR << endl;
86  cout << "ptmin = " << ptmin << endl;
87 
88  /* Open files and retrieve trees */
89  TFile *file_mc_truth = new TFile(truthFileName, "OPEN");
90  TFile *file_mc_smeared = new TFile(smearFileName, "OPEN");
91 
92  TTree* tree_truth = (TTree*)file_mc_truth->Get("EICTree");
93  TTree *tree_smeared = (TTree*)file_mc_smeared->Get("Smeared");
94 
95  tree_truth->AddFriend(tree_smeared);
96 
97  /* Setup Input Event Buffer */
98  erhic::EventPythia* event_truth(NULL);
99  Smear::Event* event_smear(NULL);
100 
101  tree_truth->SetBranchAddress("event",&event_truth);
102  tree_truth->SetBranchAddress("eventS", &event_smear);
103 
104  /* Set more buffers */
105  Double_t event_x = NAN;
106  Double_t event_q2 = NAN;
107 
108  tree_truth->SetBranchAddress("trueX",&event_x);
109  tree_truth->SetBranchAddress("trueQ2",&event_q2);
110 
111  /* Open Output File */
112  TFile *ofile = TFile::Open(outFileName,"recreate");
113  assert(ofile);
114 
115  /* Create Jet Tree */
116  Int_t _event_id = -999;
117  Int_t _event_njets = -999;
118  Double_t _event_truth_x = NAN;
119  Double_t _event_truth_q2 = NAN;
120 
121  Int_t _jet_truth_id = -999;
122  Int_t _jet_truth_ncomp = -999;
123  Int_t _jet_truth_ncharged = -999;
124  Double_t _jet_truth_e = NAN;
125  Double_t _jet_truth_et = NAN;
126  Double_t _jet_truth_eta = NAN;
127  Double_t _jet_truth_phi = NAN;
128  Double_t _jet_truth_minv = NAN;
129  Double_t _jet_truth_eem = NAN;
130  Double_t _jet_truth_rvtx = NAN;
131  Double_t _jet_truth_rmean = NAN;
132  Double_t _jet_truth_rstdev = NAN;
133 
134  Int_t _jet_smear_id = -999;
135  Int_t _jet_smear_ncomp = -999;
136  Int_t _jet_smear_ncharged = -999;
137  Double_t _jet_smear_e = NAN;
138  Double_t _jet_smear_et = NAN;
139  Double_t _jet_smear_eta = NAN;
140  Double_t _jet_smear_phi = NAN;
141  Double_t _jet_smear_minv = NAN;
142  Double_t _jet_smear_eem = NAN;
143  Double_t _jet_smear_rvtx = NAN;
144  Double_t _jet_smear_rmean = NAN;
145  Double_t _jet_smear_rstdev = NAN;
146 
147  /* Create jet tree */
148  TTree *mTree = new TTree("jets","Jet Tree");
149  mTree->Branch("event_id", &_event_id);
150  mTree->Branch("event_njets", &_event_njets);
151  mTree->Branch("event_truth_x", &_event_truth_x);
152  mTree->Branch("event_truth_q2", &_event_truth_q2);
153 
154  mTree->Branch("jet_truth_id", &_jet_truth_id);
155  mTree->Branch("jet_truth_ncomp", &_jet_truth_ncomp);
156  mTree->Branch("jet_truth_ncharged", &_jet_truth_ncharged);
157  mTree->Branch("jet_truth_e", &_jet_truth_e);
158  mTree->Branch("jet_truth_et", &_jet_truth_et);
159  mTree->Branch("jet_truth_eta", &_jet_truth_eta);
160  mTree->Branch("jet_truth_phi", &_jet_truth_phi);
161  mTree->Branch("jet_truth_minv", &_jet_truth_minv);
162  mTree->Branch("jet_truth_eem", &_jet_truth_eem);
163  mTree->Branch("jet_truth_rvtx", &_jet_truth_rvtx);
164  mTree->Branch("jet_truth_rmean", &_jet_truth_rmean);
165  mTree->Branch("jet_truth_rstdev", &_jet_truth_rstdev);
166 
167  mTree->Branch("jet_smear_id", &_jet_smear_id);
168  mTree->Branch("jet_smear_ncomp", &_jet_smear_ncomp);
169  mTree->Branch("jet_smear_ncharged", &_jet_smear_ncharged);
170  mTree->Branch("jet_smear_e", &_jet_smear_e);
171  mTree->Branch("jet_smear_et", &_jet_smear_et);
172  mTree->Branch("jet_smear_eta", &_jet_smear_eta);
173  mTree->Branch("jet_smear_phi", &_jet_smear_phi);
174  mTree->Branch("jet_smear_minv", &_jet_smear_minv);
175  mTree->Branch("jet_smear_eem", &_jet_smear_eem);
176  mTree->Branch("jet_smear_rvtx", &_jet_smear_rvtx);
177  mTree->Branch("jet_smear_rmean", &_jet_smear_rmean);
178  mTree->Branch("jet_smear_rstdev", &_jet_smear_rstdev);
179 
180  /* Loop Over Events in Simu Trees */
181  for(Int_t iEvent=0; iEvent<tree_truth->GetEntries(); iEvent++)
182  {
183  /* Read Next Event */
184  if(tree_truth->GetEntry(iEvent) <=0) break;
185 
186  if(iEvent%10000 == 0) cout << "Event " << iEvent << endl;
187 
188  /* reset tree variables */
189  _event_id = iEvent;
190  _event_njets = 0;
191  _event_truth_x = event_x;
192  _event_truth_q2 = event_q2;
193 
194  /* Create laboratory frame jet vectors */
195  vector<PseudoJet> jetcomponent_truth;
196  vector<PseudoJet> jetcomponent_smear;
197 
198  /* Loop over truth Particles, fill jet component vector */
199  for(Int_t j=0; j<event_truth->GetNTracks(); j++)
200  {
201  const Particle* inParticle = event_truth->GetTrack(j);
202 
203  if ( !inParticle)
204  continue;
205 
206  /* Select Particles for Jets */
207  if(j>10 && inParticle->GetStatus() == 1 && inParticle->GetParentIndex() != 3)
208  {
209  if( abs(inParticle->GetEta()) <= 5 && inParticle->GetE() >= 0.250 )
210  {
211  /* Truth particle: Get all information directly from particle */
212  Double_t px = inParticle->GetPx();
213  Double_t py = inParticle->GetPy();
214  Double_t pz = inParticle->GetPz();
215  Double_t E = inParticle->GetE();
216 
217  fastjet::PseudoJet pPt(px,py,pz,E);
218  pPt.set_user_index(inParticle->GetIndex());
219  jetcomponent_truth.push_back(pPt);
220  }
221  }
222  }
223 
224  /* Loop over SMEARED jet components: Energy in Calorimeter (like 'tower jets'). Fill jet component vector. */
225  for(Int_t js=0; js<event_smear->GetNTracks(); js++)
226  {
227  const Smear::ParticleMCS* inParticle = event_smear->GetTrack(js);
228 
229  if ( !inParticle)
230  continue;
231 
232  /* Select Particles for Jets */
233  if(js>10 && inParticle->GetStatus() == 1 && inParticle->GetParentIndex() != 3)
234  {
235  if( inParticle->GetE() >= 0.250 )
236  {
237  /* Calorimeter: Get energy from smeared particle, and ... */
238  Double_t E = inParticle->GetE();
239 
240 // if ( E == 0 )
241 // cout << "E == 0 found! PID = " << event_truth->GetTrack(js)->GetPdgCode() << " , E_true = " << event_truth->GetTrack(js)->GetE() << " , Eta_true = " << event_truth->GetTrack(js)->GetEta() << endl;
242 
243  /* ... get theta, phi from truth particle */
244  Double_t phi = event_truth->GetTrack(js)->GetPhi();
245  Double_t eta = event_truth->GetTrack(js)->GetEta();
246 
247  /* Recalculate px, py, pz based on calorimeter enrgy and truth direction */
248  Double_t pT = E / cosh( eta );
249  Double_t px = pT * cos( phi );
250  Double_t py = pT * sin( phi );
251  Double_t pz = pT * sinh( eta );
252 
253  fastjet::PseudoJet pPt(px,py,pz,E);
254  pPt.set_user_index(event_truth->GetTrack(js)->GetIndex());
255  jetcomponent_smear.push_back(pPt);
256  }
257  }
258  }
259 
260  /* Set Jet Definitions */
261  JetDefinition jet_def_antikt(antikt_algorithm,fastjetR);
262 
263  /* Run Clustering and Extract the Jets */
264 
265  /* Lab Frame Cluster */
266  ClusterSequence cluster_truth_antikt(jetcomponent_truth, jet_def_antikt);
267  ClusterSequence cluster_smear_antikt(jetcomponent_smear, jet_def_antikt);
268 
269  /* Lab Frame Jets*/
270  vector<PseudoJet> jets_truth_antikt = sorted_by_pt(cluster_truth_antikt.inclusive_jets(ptmin));
271  vector<PseudoJet> jets_smear_antikt = sorted_by_pt(cluster_smear_antikt.inclusive_jets(ptmin));
272 
273  /* loop over SMEARED jets */
274  _event_njets = jets_smear_antikt.size();
275  for ( unsigned ijet = 0; ijet < _event_njets; ijet++ )
276  {
277  PseudoJet* jetMatch = find_matching_jet( &(jets_smear_antikt.at(ijet)), &jets_truth_antikt );
278 
279  /* Set SMEARED jet variables */
280  _jet_smear_id = ijet;
281  _jet_smear_ncomp = -999;
282  _jet_smear_ncharged = -999;
283  _jet_smear_e = jets_smear_antikt.at(ijet).E();
284  _jet_smear_et = jets_smear_antikt.at(ijet).Et();
285  _jet_smear_eta = jets_smear_antikt.at(ijet).eta();
286  _jet_smear_phi = jets_smear_antikt.at(ijet).phi_std();
287  _jet_smear_minv = jets_smear_antikt.at(ijet).m();
288  _jet_smear_eem = NAN;
289  _jet_smear_rvtx = NAN;
290  _jet_smear_rmean = NAN;
291  _jet_smear_rstdev = NAN;
292 
293  /* Set TRUTH jet variables */
294  _jet_truth_id = -999;
295  _jet_truth_ncomp = -999;
296  _jet_truth_ncharged = -999;
297  _jet_truth_e = NAN;
298  _jet_truth_et = NAN;
299  _jet_truth_eta = NAN;
300  _jet_truth_phi = NAN;
301  _jet_truth_minv = NAN;
302  _jet_truth_eem = NAN;
303  _jet_truth_rvtx = NAN;
304  _jet_truth_rmean = NAN;
305  _jet_truth_rstdev = NAN;
306 
307  if ( jetMatch )
308  {
309  _jet_truth_id = ijet;
310  _jet_truth_ncomp = -999;
311  _jet_truth_ncharged = -999;
312  _jet_truth_e = jetMatch->E();
313  _jet_truth_et = jetMatch->Et();
314  _jet_truth_eta = jetMatch->eta();
315  _jet_truth_phi = jetMatch->phi_std();
316  _jet_truth_minv = jetMatch->m();
317  _jet_truth_eem = NAN;
318  _jet_truth_rvtx = NAN;
319  _jet_truth_rmean = NAN;
320  _jet_truth_rstdev = NAN;
321  }
322 
323  /* Fill tree */
324  mTree->Fill();
325  }
326  }
327 
328  /* Write output tree and close file */
329  mTree->Write();
330  ofile->Close();
331 
332  return 0;
333 }