Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MatchJets.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MatchJets.C
1 #include <iostream>
2 #include <vector>
3 
4 using namespace std;
5 
6 
7 Int_t MatchJets(TString infile="../macro/output.root", TString outfile="matched.root", Float_t JetParameter=0.4, Float_t pTHard=30.0)
8 {
9 
10  // Set event selection cuts
11  double dRMax = 0.75*JetParameter;
12  float pTmin_truth = 10;
13  float pTmin_reco = 5;
14  float pThard_min, pThard_max;
15  float weight;
16  if(pTHard == 10){
17  pThard_min = 10;
18  pThard_max = 30;
19  weight = 3.210e-6;
20  }
21  else{
22  pThard_min = 30;
23  pThard_max = 999;
24  weight = 2.178e-9;
25  }
26 
27 
28  //=======================================================
29  //================ Input File ===========================
30  //=======================================================
31  std::cout << "Input file: " << infile.Data() << std::endl;
32 
33  // open input file
34  TFile fin(infile.Data(), "READ");
35  if(!fin.IsOpen())
36  {
37  std::cout << "Error: could not open input file" << std::endl;
38  exit(-1);
39  }
40 
41  // get tree
42  TTree *tree = (TTree*)fin.Get("T");
43  if(!tree)
44  {
45  std::cout << "Error: could not find tree" << std::endl;
46  exit(-1);
47  }
48 
49  int event = 0;
50  int centrality = 0;
51  float event_leading_truth_pt = 0;
52  double rho_area = 0;
53  double rho_area_sigma = 0;
54  double rho_mult = 0;
55 
56  // reco jet variables
57  std::vector<float> * iter_eta = 0;
58  std::vector<float> * iter_phi = 0;
59  std::vector<float> * iter_pt = 0;
60  std::vector<float> * iter_pt_unsub = 0;
61 
62  // mult jet variables
63  std::vector<int> * mult_ncomponent = 0;
64  std::vector<float> * mult_eta = 0;
65  std::vector<float> * mult_phi = 0;
66  std::vector<float> * mult_pt = 0;
67  std::vector<float> * mult_pt_unsub = 0;
68  std::vector<int> * mult_nsignal = 0;
69 
70  // truth jet variables
71  std::vector<int> * truth_ncomponent = 0;
72  std::vector<float> * truth_eta = 0;
73  std::vector<float> * truth_phi = 0;
74  std::vector<float> * truth_pt = 0;
75 
76  //rhoA jet variables
77  std::vector<float> * rhoA_eta = 0;
78  std::vector<float> * rhoA_phi = 0;
79  std::vector<float> * rhoA_pt = 0;
80  std::vector<float> * rhoA_pt_unsub = 0;
81  std::vector<float> * rhoA_area = 0;
82 
83  tree->SetBranchAddress("event", &event);
84  tree->SetBranchAddress("event_leading_truth_pt", &event_leading_truth_pt);
85  tree->SetBranchAddress("centrality", &centrality);
86  tree->SetBranchAddress("truth_ncomponent", &truth_ncomponent);
87  tree->SetBranchAddress("truth_eta", &truth_eta);
88  tree->SetBranchAddress("truth_phi", &truth_phi);
89  tree->SetBranchAddress("truth_pt", &truth_pt);
90  tree->SetBranchAddress("rho_area", &rho_area);
91  tree->SetBranchAddress("rho_area_sigma", &rho_area_sigma);
92  tree->SetBranchAddress("rhoA_eta", &rhoA_eta);
93  tree->SetBranchAddress("rhoA_phi", &rhoA_phi);
94  tree->SetBranchAddress("rhoA_pt", &rhoA_pt);
95  tree->SetBranchAddress("rhoA_area", &rhoA_area);
96  tree->SetBranchAddress("rhoA_pt_unsub", &rhoA_pt_unsub);
97  tree->SetBranchAddress("rho_mult", &rho_mult);
98  tree->SetBranchAddress("mult_ncomponent", &mult_ncomponent);
99  tree->SetBranchAddress("mult_nsignal", &mult_nsignal);
100  tree->SetBranchAddress("mult_eta", &mult_eta);
101  tree->SetBranchAddress("mult_phi", &mult_phi);
102  tree->SetBranchAddress("mult_pt", &mult_pt);
103  tree->SetBranchAddress("mult_pt_unsub", &mult_pt_unsub);
104  tree->SetBranchAddress("iter_eta", &iter_eta);
105  tree->SetBranchAddress("iter_phi", &iter_phi);
106  tree->SetBranchAddress("iter_pt", &iter_pt);
107  tree->SetBranchAddress("iter_pt_unsub", &iter_pt_unsub);
108 
109 
110  //=======================================================
111 
112  // output file
113  std::cout << "Output file: " << outfile.Data() << std::endl;
114  TFile *fout = new TFile(outfile.Data(), "RECREATE");
115  TTree *out_tree = new TTree("T", "MatchedJets");
116 
117  // output variables
118  int out_centrality = 0;
119  double out_rho_area = 0;
120  double out_rho_mult = 0;
121 
122  std::vector<float> matched_pt_iter, matched_pt_area, matched_pt_mult;
123  std::vector<float> matched_pt_iter_unsub, matched_pt_area_unsub, matched_pt_mult_unsub;
124  std::vector<float> matched_truth_pt_iter, matched_truth_pt_area, matched_truth_pt_mult;
125  std::vector<float> unmatched_pt_iter, unmatched_pt_area, unmatched_pt_mult;
126  std::vector<float> unmatched_pt_iter_unsub, unmatched_pt_area_unsub, unmatched_pt_mult_unsub;
127  std::vector<float> missed_truth_pt_iter, missed_truth_pt_area, missed_truth_pt_mult;
128 
129  // set branch addresses
130  out_tree->Branch("weight", &weight);
131  out_tree->Branch("centrality", &out_centrality);
132  out_tree->Branch("rho_area", &out_rho_area);
133  out_tree->Branch("rho_mult", &out_rho_mult);
134  out_tree->Branch("matched_pt_iter", &matched_pt_iter);
135  out_tree->Branch("matched_pt_area", &matched_pt_area);
136  out_tree->Branch("matched_pt_mult", &matched_pt_mult);
137  out_tree->Branch("matched_pt_iter_unsub", &matched_pt_iter_unsub);
138  out_tree->Branch("matched_pt_area_unsub", &matched_pt_area_unsub);
139  out_tree->Branch("matched_pt_mult_unsub", &matched_pt_mult_unsub);
140  out_tree->Branch("matched_truth_pt_iter", &matched_truth_pt_iter);
141  out_tree->Branch("matched_truth_pt_area", &matched_truth_pt_area);
142  out_tree->Branch("matched_truth_pt_mult", &matched_truth_pt_mult);
143  out_tree->Branch("unmatched_pt_iter", &unmatched_pt_iter);
144  out_tree->Branch("unmatched_pt_area", &unmatched_pt_area);
145  out_tree->Branch("unmatched_pt_mult", &unmatched_pt_mult);
146  out_tree->Branch("unmatched_pt_iter_unsub", &unmatched_pt_iter_unsub);
147  out_tree->Branch("unmatched_pt_area_unsub", &unmatched_pt_area_unsub);
148  out_tree->Branch("unmatched_pt_mult_unsub", &unmatched_pt_mult_unsub);
149  out_tree->Branch("missed_truth_pt_iter", &missed_truth_pt_iter);
150  out_tree->Branch("missed_truth_pt_area", &missed_truth_pt_area);
151  out_tree->Branch("missed_truth_pt_mult", &missed_truth_pt_mult);
152 
153  // get number of entries
154  int nEntries = tree->GetEntries();
155  for (int iEvent = 0; iEvent < nEntries; iEvent++){
156  // get entry
157  tree->GetEntry(iEvent);
158 
159  // event selection
160  if(event_leading_truth_pt < pThard_min || event_leading_truth_pt > pThard_max) continue;
161 
162  // clear output variables
163  out_centrality = centrality;
164  out_rho_area = rho_area;
165  out_rho_mult = rho_mult;
166  // clear vectors
167  matched_pt_iter.clear();
168  matched_pt_area.clear();
169  matched_pt_mult.clear();
170  matched_pt_iter_unsub.clear();
171  matched_pt_area_unsub.clear();
172  matched_pt_mult_unsub.clear();
173  matched_truth_pt_iter.clear();
174  matched_truth_pt_area.clear();
175  matched_truth_pt_mult.clear();
176  unmatched_pt_iter.clear();
177  unmatched_pt_area.clear();
178  unmatched_pt_mult.clear();
179  unmatched_pt_iter_unsub.clear();
180  unmatched_pt_area_unsub.clear();
181  unmatched_pt_mult_unsub.clear();
182  missed_truth_pt_iter.clear();
183  missed_truth_pt_area.clear();
184  missed_truth_pt_mult.clear();
185 
186  // vector of matched reco jets for each method
187  std::vector<unsigned int> matched_area_jets, matched_mult_jets, matched_iter_jets;
188 
189  // loop over truth jets
190  for(unsigned int ijet =0; ijet < truth_pt->size(); ijet++)
191  {
192  // apply truth cut
193  if(truth_pt->at(ijet) < pTmin_truth) continue;
194 
195  float dr = 999;
196  float matched_area_pt, matched_area_unsubpt;
197  float matched_mult_pt, matched_mult_unsubpt;
198  float matched_iter_pt, matched_iter_unsubpt;
199  unsigned int matched_idx_area, matched_idx_mult, matched_idx_iter;
200 
201  // match area jets to truth jets
202  for(unsigned int jjet = 0; jjet < rhoA_pt->size(); jjet++)
203  {
204  // apply reco jet cut
205  if(rhoA_pt->at(jjet) < pTmin_reco) continue;
206 
207  // dr calculation
208  float dphi = rhoA_phi->at(jjet) - truth_phi->at(ijet);
209  float deta = rhoA_eta->at(jjet) - truth_eta->at(ijet);
210  if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
211  if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
212  float dr_tmp = TMath::Sqrt(dphi*dphi + deta*deta);
213  if(dr_tmp < dr && dr_tmp < dRMax)
214  {
215  dr = dr_tmp;
216  matched_area_pt = rhoA_pt->at(jjet);
217  matched_area_unsubpt = rhoA_pt_unsub->at(jjet);
218  matched_idx_area = jjet;
219  }
220  }
221  if(dr > dRMax){
222  missed_truth_pt_area.push_back(truth_pt->at(ijet));
223  }
224  else{
225  matched_truth_pt_area.push_back(truth_pt->at(ijet));
226  matched_pt_area.push_back(matched_area_pt);
227  matched_pt_area_unsub.push_back(matched_area_unsubpt);
228  matched_area_jets.push_back(matched_idx_area);
229  }
230 
231  // reset dr
232  dr = 999;
233  // match mult jets to truth jets
234  for(unsigned int jjet = 0; jjet < mult_pt->size(); jjet++)
235  {
236  // apply reco jet cut
237  if(mult_pt->at(jjet) < pTmin_reco) continue;
238 
239  // dr calculation
240  float dphi = mult_phi->at(jjet) - truth_phi->at(ijet);
241  float deta = mult_eta->at(jjet) - truth_eta->at(ijet);
242  if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
243  if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
244  float dr_tmp = TMath::Sqrt(dphi*dphi + deta*deta);
245  if(dr_tmp < dr && dr_tmp < dRMax)
246  {
247  dr = dr_tmp;
248  matched_mult_pt = mult_pt->at(jjet);
249  matched_mult_unsubpt = mult_pt_unsub->at(jjet);
250  matched_idx_mult = jjet;
251  }
252  }
253  if (dr > dRMax){
254  missed_truth_pt_mult.push_back(truth_pt->at(ijet));
255  }
256  else{
257  matched_truth_pt_mult.push_back(truth_pt->at(ijet));
258  matched_pt_mult.push_back(matched_mult_pt);
259  matched_pt_mult_unsub.push_back(matched_mult_unsubpt);
260  matched_mult_jets.push_back(matched_idx_mult);
261  }
262 
263  // // reset dr
264  dr = 999;
265  // match iter jets to truth jets
266  for (unsigned int jjet = 0; jjet < iter_pt->size(); jjet++)
267  {
268  // apply reco jet cut
269  if(iter_pt->at(jjet) < pTmin_reco) continue;
270 
271  // dr calculation
272  float dphi = iter_phi->at(jjet) - truth_phi->at(ijet);
273  float deta = iter_eta->at(jjet) - truth_eta->at(ijet);
274  if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
275  if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
276  float dr_tmp = TMath::Sqrt(dphi*dphi + deta*deta);
277  if(dr_tmp < dr && dr_tmp < dRMax)
278  {
279  dr = dr_tmp;
280  matched_iter_pt = iter_pt->at(jjet);
281  matched_iter_unsubpt = iter_pt_unsub->at(jjet);
282  matched_idx_iter = jjet;
283  }
284  }
285  if (dr > dRMax){
286  missed_truth_pt_iter.push_back(truth_pt->at(ijet));
287  }
288  else{
289  matched_truth_pt_iter.push_back(truth_pt->at(ijet));
290  matched_pt_iter.push_back(matched_iter_pt);
291  matched_pt_iter_unsub.push_back(matched_iter_unsubpt);
292  matched_iter_jets.push_back(matched_idx_iter);
293  }
294 
295 
296  }
297  // fill unmatched area jets
298  for(unsigned int jjet = 0; jjet < rhoA_pt->size(); jjet++)
299  {
300  if(std::find(matched_area_jets.begin(), matched_area_jets.end(), jjet) != matched_area_jets.end()) continue;
301  if(rhoA_pt->at(jjet) < pTmin_reco) continue;
302  unmatched_pt_area.push_back(rhoA_pt->at(jjet));
303  unmatched_pt_area_unsub.push_back(rhoA_pt_unsub->at(jjet));
304  }
305  // fill unmatched reco mult jets
306  for(unsigned int jjet = 0; jjet < mult_pt->size(); jjet++)
307  {
308  if(std::find(matched_mult_jets.begin(), matched_mult_jets.end(), jjet) != matched_mult_jets.end()) continue;
309  if(mult_pt->at(jjet) < pTmin_reco) continue;
310  unmatched_pt_mult.push_back(mult_pt->at(jjet));
311  unmatched_pt_mult_unsub.push_back(mult_pt_unsub->at(jjet));
312  }
313  // fill unmatched reco iter jets
314  for(unsigned int jjet = 0; jjet < iter_pt->size(); jjet++)
315  {
316  if(std::find(matched_iter_jets.begin(), matched_iter_jets.end(), jjet) != matched_iter_jets.end()) continue;
317  if(iter_pt->at(jjet) < pTmin_reco) continue;
318  unmatched_pt_iter.push_back(iter_pt->at(jjet));
319  unmatched_pt_iter_unsub.push_back(iter_pt_unsub->at(jjet));
320  }
321 
322  // fill tree
323  out_tree->Fill();
324  }
325 
326  // write output file
327  fout->Write();
328  fout->Close();
329  // close input file
330  fin.Close();
331  return 0;
332 }
333