Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
likelihood.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file likelihood.C
1 #include "tau_commons.h"
2 
11 double Average(vector<double> v)
12 {
13  double sum=0;
14  for(int i=0;(unsigned)i<v.size();i++)
15  sum+=v[i];
16  return sum/v.size();
17 }
18 
19 
20 int likelihood()
21 {
22  gStyle->SetOptStat(0);
23 
24  unsigned col1 = kOrange+7;
25  unsigned col2 = kBlue+2;
26 
27  // Labels for geant files //
28  string type[4] = {"3pion","SM","DIScharged","DISneutral"};
29  string seed[10] = {"1","2","3","4","5","6","7","8","9","10"};
30  string R_max = "5";
31 
32  /* open inout files and merge trees */
33  TChain chain("event");
34  TChain chain_tau("event");
35  TChain t("event");
36 
37  // File for event topology of ID tau events
38  ofstream myfile;
39  string filename = "./data/Event_topology.csv";
40  myfile.open(filename.c_str());
41 
42  // Loop over geant files and split LQ and SM events into two TChains used for jet characterization //
43  for (int b = 0; b<2; b++){
44  for (int g = 0; g<10; g++){
45  // Skip file that is used to test //
46  if(b == 1 && g==0) continue;
47  string file = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/data/LeptoAna_p250_e20_1000events_"+seed[g]+"seed_"+type[b]+"_r0"+R_max+".root";
48 
49  if(b==0) chain_tau.Add(file.c_str());
50  else chain.Add(file.c_str());
51  }
52  }
53 
54  for(int b = 1; b<2; b++){
55  for(int g=0; g<1; g++){
56  // File to test identification //
57  string inputFile = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/data/LeptoAna_p250_e20_1000events_"+seed[g]+"seed_"+type[b]+"_r05.root";
58 
59  t.Add(inputFile.c_str());
60  }
61  }
62 
63 
64 
65  /* Create temporary canvas */
66  TCanvas* ctemp = new TCanvas();
67 
68  vector< TString > observables;
69  vector< TString > observables_name;
70 
71  vector< float > plots_ymax;
72  vector< float > plots_xmin;
73  vector< float > plots_xmax;
74 
75  /* R_max */
76  observables.push_back( "tracks_rmax_r04" );
77  observables_name.push_back( "R_{track}^{max}" );
78  plots_ymax.push_back(0.06);
79  plots_xmin.push_back(0);
80  plots_xmax.push_back(0.4);
81 
82  /* Number of tracks */
83  observables.push_back( "tracks_count_r04" );
84  observables_name.push_back( "N_{tracks}" );
85  plots_ymax.push_back(0.9);
86  plots_xmin.push_back(0);
87  plots_xmax.push_back(10);
88 
89  /* Charge sum from tracks */
90  observables.push_back( "tracks_chargesum_r04" );
91  observables_name.push_back( "#Sigma q_{tracks}" );
92  plots_ymax.push_back(0.9);
93  plots_xmin.push_back(-5);
94  plots_xmax.push_back(5);
95 
96  /* Jet radius */
97  observables.push_back( "jetshape_radius" );
98  observables_name.push_back( "R_{jet}" );
99  plots_ymax.push_back(0.08);
100  plots_xmin.push_back(0);
101  plots_xmax.push_back(0.5);
102 
103  /* Jetshape */
104  observables.push_back( "jetshape_econe_r01 / jetshape_econe_r05" );
105  observables_name.push_back( "E_{cone}^{R<0.1} / E_{cone}^{R<0.5}" );
106  plots_ymax.push_back(0.08);
107  plots_xmin.push_back(0);
108  plots_xmax.push_back(1);
109 
110  /* Jetshape */
111  observables.push_back( "jetshape_econe_r02 / jetshape_econe_r05" );
112  observables_name.push_back( "E_{cone}^{R<0.2} / E_{cone}^{R<0.5}" );
113  plots_ymax.push_back(0.08);
114  plots_xmin.push_back(0);
115  plots_xmax.push_back(1);
116 
117  /* Jet eta */
118  observables.push_back( "jet_eta" );
119  observables_name.push_back( "#eta_{jet}" );
120  plots_ymax.push_back(0.1);
121  plots_xmin.push_back(-1);
122  plots_xmax.push_back(3);
123 
124  /* Jet mass */
125  observables.push_back( "jet_minv" );
126  observables_name.push_back( "Mass_{jet} (GeV)" );
127  plots_ymax.push_back(0.1);
128  plots_xmin.push_back(0);
129  plots_xmax.push_back(10);
130 
131  /* Jet energy */
132  observables.push_back( "jet_etotal" );
133  observables_name.push_back( "E_{jet} (GeV)" );
134  plots_ymax.push_back(0.08);
135  plots_xmin.push_back(0);
136  plots_xmax.push_back(70);
137 
138  /* Jet vertex */
139  observables.push_back( "tracks_vertex" );
140  observables_name.push_back( "vertex distance (cm)" );
141  plots_ymax.push_back(0.6);
142  plots_xmin.push_back(0);
143  plots_xmax.push_back(1);
144 
145 
146  /* Plot distributions */
147  TString name_uds_base("h_uds_");
148  TString name_tau_base("h_tau_");
149 
150  const int nvar = observables.size();
151 
152  TH1F* h_uds[nvar];
153  TH1F* h_tau[nvar];
154 
155  TString name_uds_i = name_uds_base;
156  TString name_tau_i = name_tau_base;
157 
158  /* create histograms */
159  for ( unsigned j = 0; j < observables.size(); j++ ){
160  name_uds_i = name_uds_base;
161  name_tau_i = name_tau_base;
162 
163  name_uds_i.Append(j);
164  name_tau_i.Append(j);
165 
166  h_uds[j] = new TH1F( name_uds_i, "", 100, plots_xmin.at(j), plots_xmax.at(j));
167  h_tau[j] = new TH1F( name_tau_i, "", 100, plots_xmin.at(j), plots_xmax.at(j));
168  }
169 
170  /* Fill histograms*/
171  for ( unsigned i = 0; i < observables.size(); i++ )
172  {
173  //cout << "Plotting " << observables.at(i) << endl;
174 
175  ctemp->cd();
176 
177  name_uds_i = name_uds_base;
178  name_tau_i = name_tau_base;
179 
180  name_uds_i.Append(i);
181  name_tau_i.Append(i);
182 
183  // SM parton jet observable histograms //
184  chain.Draw( observables.at(i) + " >> " + name_uds_i, tau_commons::select_true_uds && tau_commons::select_accept_jet, "goff");
185  h_uds[i]->Scale(1./h_uds[i]->Integral());
186  h_uds[i]->GetXaxis()->SetTitle( observables_name.at(i) );
187  h_uds[i]->GetYaxis()->SetRangeUser(0, plots_ymax.at(i) );
188  h_uds[i]->SetLineColor(col1);
189  h_uds[i]->SetFillColor(col1);
190  h_uds[i]->SetFillStyle(1001);
191  h_uds[i]->GetYaxis()->SetTitle("# entries / #Sigma entries");
192  h_uds[i]->GetYaxis()->SetTitleOffset(1.5);
193 
194  // LQ tau jet observable histograms //
195  chain_tau.Draw( observables.at(i) + " >> " + name_tau_i, tau_commons::select_true_tau && tau_commons::select_accept_jet, "goff");
196  h_tau[i]->Scale(1./h_tau[i]->Integral());
197  h_tau[i]->SetLineColor(col2);
198  h_tau[i]->SetFillColor(col2);
199  h_tau[i]->SetFillStyle(0);
200 
201 
202  /* create Canvas and draw histograms */
203  TCanvas *c1 = new TCanvas();
204 
205  h_uds[i]->DrawClone("");
206  h_tau[i]->DrawClone("sames");
207  gPad->RedrawAxis();
208 
209  TLegend *legend = new TLegend(0.4,0.6,0.7,0.89);
210  legend->AddEntry(h_uds[i],"SM parton jet","l");
211  legend->AddEntry(h_tau[i],"LQ #tau jet","l");
212  legend->Draw();
213 
214  // c1->Print( Form( "plots/compare_observables_%d.eps", i ) );
215  // c1->Print( Form( "plots/compare_observables_%d.png", i ) );
216  }
217 
218  // Define jet and //
219  vector<float> * tracks_rmax;
220  vector<float> * tracks_count;
221  vector<float> * tracks_chargesum;
222  vector<float> * tracks_vertex;
223  vector<float> * jetshape_radius;
224  vector<float> * jetshape_econe_1;
225  vector<float> * jetshape_econe_2;
226  vector<float> * jetshape_econe_5;
227  vector<float> * jet_eta;
228  vector<float> * jet_phi;
229  vector<float> * jet_minv;
230  vector<float> * jet_etotal;
231  vector<float> * jet_ptrans;
232  vector<float> * evtgen_pid;
233 
234  // Define event variables //
235  float ptmiss;
236  float miss_phi;
237  float is_lq;
238 
239  // See variables from tree //
240  t.SetBranchAddress("tracks_rmax_R",&tracks_rmax);
241  t.SetBranchAddress("tracks_count_R",&tracks_count);
242  t.SetBranchAddress("tracks_chargesum_R",&tracks_chargesum);
243  t.SetBranchAddress("tracks_vertex",&tracks_vertex);
244  t.SetBranchAddress("jetshape_radius",&jetshape_radius);
245  t.SetBranchAddress("jetshape_econe_r01",&jetshape_econe_1);
246  t.SetBranchAddress("jetshape_econe_r02",&jetshape_econe_2);
247  t.SetBranchAddress("jetshape_econe_r05",&jetshape_econe_5);
248  t.SetBranchAddress("jet_eta",&jet_eta);
249  t.SetBranchAddress("jet_phi",&jet_phi);
250  t.SetBranchAddress("jet_minv",&jet_minv);
251  t.SetBranchAddress("jet_etotal",&jet_etotal);
252  t.SetBranchAddress("evtgen_pid",&evtgen_pid);
253  t.SetBranchAddress("jet_ptrans",&jet_ptrans);
254  t.SetBranchAddress("Et_miss",&ptmiss);
255  t.SetBranchAddress("Et_miss_phi",&miss_phi);
256  t.SetBranchAddress("is_lq_event",&is_lq);
257 
258  /* identification counts */
259  int tau_as_tau = 0;
260  int tau_as_DIS = 0;
261  int DIS_as_DIS = 0;
262  int DIS_as_tau = 0;
263  int id_e = 0;
264 
265  int tot_jets = 0;
266 
267  // Define likelihood cut for tau jet identification //
268  double lik_cut = 0.71;
269 
270  int n_tot = t.GetEntries();
271 
272  // Loop over all events //
273  for( int k = 0; k < n_tot; k++){
274  t.GetEntry(k);
275 
276  // Define some event variables //
277  double delta_phi;
278  bool is_e = false;
279  bool is_miss = false;
280  vector<double> jet_lik;
281  vector<double> jet_id;
282  vector<double> miss_jet_phi;
283 
284  // Loop over jets in event to look for electron //
285  for(int l = 0; l < tracks_rmax->size(); l++){
286  int count = (*tracks_count)[l];
287  double rmax = (*tracks_rmax)[l];
288  int chargesum = (*tracks_chargesum)[l];
289  double radius = (*jetshape_radius)[l];
290  double eta = (*jet_eta)[l];
291  double phi = (*jet_phi)[l];
292  double etotal = (*jet_etotal)[l];
293  double vertex = (*tracks_vertex)[l];
294  double ptrans = (*jet_ptrans)[l];
295  int pid = (*evtgen_pid)[l];
296 
297  // If any jets pass these cuts then this is considered an electron event //
298  if(count == 1 && chargesum == -1 && ptrans > 5 && vertex < 0.03 && radius < 0.05) is_e = true;
299 
300  }
301 
302  // Loop over all jets in events //
303  for(int l = 0; l < tracks_rmax->size(); l++){
304 
305  // Likelihood variable //
306  vector<double> likelihood;
307 
308  // Get jet observables //
309  double rmax = (*tracks_rmax)[l];
310  int count = (*tracks_count)[l];
311  int chargesum = (*tracks_chargesum)[l];
312  double vertex = (*tracks_vertex)[l];
313  double radius = (*jetshape_radius)[l];
314  double econe_1 = (*jetshape_econe_1)[l];
315  double econe_2 = (*jetshape_econe_2)[l];
316  double econe_5 = (*jetshape_econe_5)[l];
317  double eta = (*jet_eta)[l];
318  double phi = (*jet_phi)[l];
319  double minv = (*jet_minv)[l];
320  double etotal = (*jet_etotal)[l];
321  double ptrans = (*jet_ptrans)[l];
322  int pid = (*evtgen_pid)[l];
323  tot_jets++;
324 
325  // If electron event then ID as parton and continue //
326  if(is_e && pid == 15) tau_as_DIS++;
327  if(is_e && pid != 15) DIS_as_DIS++;
328  if(is_e) continue;
329 
330 
331  // x bin number //
332  int x_bin;
333 
334  // Find where jet observable lies on characteristic histograms and use histogram heights to find likelihood //
335  // Currently count, chargesum, eta, and vertex give best results. Rest are left commented out //
336 
337  //x_bin = h_tau[0]->GetXaxis()->FindBin(rmax);
338  //likelihood.push_back( h_tau[0]->GetBinContent(x_bin) / (h_tau[0]->GetBinContent(x_bin) + h_uds[0]->GetBinContent(x_bin)));
339 
340  x_bin = h_tau[1]->GetXaxis()->FindBin(count);
341  likelihood.push_back( h_tau[1]->GetBinContent(x_bin) / (h_tau[1]->GetBinContent(x_bin) + h_uds[1]->GetBinContent(x_bin)));
342 
343  x_bin = h_tau[2]->GetXaxis()->FindBin(chargesum);
344  likelihood.push_back( h_tau[2]->GetBinContent(x_bin) / (h_tau[2]->GetBinContent(x_bin) + h_uds[2]->GetBinContent(x_bin)));
345 
346  //x_bin = h_tau[3]->GetXaxis()->FindBin(radius);
347  //likelihood.push_back( h_tau[3]->GetBinContent(x_bin) / (h_tau[3]->GetBinContent(x_bin) + h_uds[3]->GetBinContent(x_bin)));
348  /*
349  x_bin = h_tau[4]->GetXaxis()->FindBin(econe_1/econe_5);
350  likelihood.push_back( h_tau[4]->GetBinContent(x_bin) / (h_tau[4]->GetBinContent(x_bin) + h_uds[4]->GetBinContent(x_bin)));
351 
352  x_bin = h_tau[5]->GetXaxis()->FindBin(econe_2/econe_5);
353  likelihood.push_back( h_tau[5]->GetBinContent(x_bin) / (h_tau[5]->GetBinContent(x_bin) + h_uds[5]->GetBinContent(x_bin)));
354  */
355 
356  x_bin = h_tau[6]->GetXaxis()->FindBin(eta);
357  likelihood.push_back( h_tau[6]->GetBinContent(x_bin) / (h_tau[6]->GetBinContent(x_bin) + h_uds[6]->GetBinContent(x_bin)));
358 
359  //x_bin = h_tau[7]->GetXaxis()->FindBin(minv);
360  //likelihood.push_back( h_tau[7]->GetBinContent(x_bin) / (h_tau[7]->GetBinContent(x_bin) + h_uds[7]->GetBinContent(x_bin)));
361 
362  //x_bin = h_tau[8]->GetXaxis()->FindBin(etotal);
363  //likelihood.push_back( h_tau[8]->GetBinContent(x_bin) / (h_tau[8]->GetBinContent(x_bin) + h_uds[8]->GetBinContent(x_bin)));
364 
365  // Exclude jets where vertex = inf which causes problems //
366  x_bin = h_tau[9]->GetXaxis()->FindBin(vertex);
367  if(vertex == vertex) likelihood.push_back( h_tau[9]->GetBinContent(x_bin) / (h_tau[9]->GetBinContent(x_bin) + h_uds[9]->GetBinContent(x_bin)));
368 
369 
370  // Define total likelihood as average of likelihoods of observables //
371  double avg = Average(likelihood);
372 
373  // Make sure final likelihood is real finite number //
374  if(avg != avg) continue;
375 
376  // Only use jets with transverse momentum > 5 GeV //
377  // Write jet likelihoods to array for event //
378  if(ptrans < 5) avg = 0;
379  if(ptrans > 5) {
380  jet_lik.push_back(avg);
381  jet_id.push_back(pid);
382  miss_jet_phi.push_back(phi);
383  }
384 
385  // If pt < 5 call parton jet //
386  if(ptrans < 5 && pid == 15) tau_as_DIS++;
387  if(ptrans < 5 && pid != 15) DIS_as_DIS++;
388 
389  // This block uses hard cuts on variables for identification //
390  /*
391  if(count == 3 && chargesum == -1 && ptrans > 5 && vertex > 0.028 && eta < 0.82)) {
392  if(pid == 15) tau_as_tau++;
393  else DIS_as_tau++;
394  }
395  else {
396  if(pid == 15) tau_as_DIS++;
397  else DIS_as_DIS++;
398  }
399  */
400 
401  }
402 
403  // Increment if electron is in event //
404  if(is_e) id_e++;
405 
406  // This block is used for likelihood identification //
407  // Define some variables for jet selection //
408  int a = -1;
409  int b = -1;
410  double max_lik = 0;
411  int max_lik_i;
412  bool is_tau = false;
413 
414  // If no jets with pt > 5 then skip //
415  if(jet_lik.size() == 0) continue;
416 
417  // Loop over all jets in event with pt > 5 //
418  for(int m = 0; m < jet_lik.size(); m++){
419  // Find if there is a tau in this event //
420  if(jet_id[m] == 15) is_tau = true;
421  // Find jet with maximum tau likelihood //
422  if(jet_lik[m] > max_lik) {
423  max_lik = jet_lik[m];
424  max_lik_i = m;
425  }
426  // There is a problem with same jets being tagged with two different true values
427  // Here we check if one of the jets is repeated //
428  // Loop over all jets to check if two have exact same likelihood //
429  for(int n = 0; n < jet_lik.size(); n++){
430  // If different jets have same likelihood then set a and b to the jet indices //
431  if(m!=n && jet_lik[m]==jet_lik[n]) {
432  a = m;
433  b = n;
434  // If either is tau jet then set them both to a tau jet id //
435  if(jet_id[m] == 15 || jet_id[n] 15) {
436  jet_id[m] == 15;
437  jet_id[n] ==15;
438  }
439  // Otherwise set them both to 0 id //
440  else {
441  jet_id[m] = 0;
442  jet_id[n] = 0;
443  }
444  }
445  }
446  }
447  // End loop over all jets with pt > 5 //
448 
449  // Check if maximum likelihood jet passes cut then ID as tau //
450  if(jet_lik[max_lik_i] > lik_cut){
451  // Find ID tau jet as true tau //
452  if(jet_id[max_lik_i] == 15) tau_as_tau++;
453  // Otherwise find ID tau jet as true parton //
454  else DIS_as_tau++;
455  // All other jets are ID as parton, minus one that was ID as tau //
456  DIS_as_DIS = DIS_as_DIS + jet_lik.size() - 1;
457  // Set event topology to be recorded for later study //
458  is_miss = true;
459  }
460  // If fails likelihood cut then ID as parton
461  else{
462  // If a tau was in the event then it is ID as parton with all other jets //
463  if(is_tau) tau_as_DIS++;
464  if(is_tau) DIS_as_DIS = DIS_as_DIS + jet_lik.size() - 1;
465  // If no tau then every jet is parton and ID correctly //
466  if(!is_tau) DIS_as_DIS = DIS_as_DIS + jet_lik.size();
467  }
468 
469  // Check if duplicate jets were found and subtract one from ID jets if it was //
470  if(a > -1 && b > -1){
471  if(jet_id[a] == 15 && jet_lik[max_lik_i] > lik_cut) DIS_as_tau--;
472  if(jet_id[a] != 15 && jet_lik[max_lik_i] < lik_cut) DIS_as_DIS--;
473  }
474 
475  // If jet is ID as tau then calculate ptmiss and acoplanarity for this event //
476  if(is_miss){
477  delta_phi = fabs(miss_phi - miss_jet_phi[max_lik_i]) * 180 / TMath::Pi();
478  if(delta_phi > 180) delta_phi = 360 - delta_phi;
479  // Write variables with true LQ tagging to file for later study //
480  myfile<<ptmiss<<","<<delta_phi<<","<<is_lq<<endl;
481  }
482  }
483 
484  // Define efficiency and false positive rate //
485  double efficiency = tau_as_tau*1.0/(n_tot*1.0);
486  double FP;
487 
488  if(tau_as_tau == 0 && DIS_as_tau == 0) FP = 0;
489  else FP = DIS_as_tau*1.0/((tau_as_tau+DIS_as_tau)*1.0);
490 
491 
492  // Confusion Matrix //
493  cout<<" | Act Tau Act Parton"<<endl;
494  cout<<"------------------------------"<<endl;
495  cout<<"ID Tau | "<<tau_as_tau<<" "<<DIS_as_tau<<endl;
496  cout<<"ID Parton| "<<tau_as_DIS<<" "<<DIS_as_DIS<<endl;
497 
498  cout<<endl;
499 
500  // Write out some important variables //
501  cout<<"Tau Recovery Rate: "<<efficiency<<endl;
502  cout<<"False Positive Rate: "<<FP<<endl;
503  cout<<"Identified Electron Events : "<<id_e<<endl;
504  cout<<"Total Jets: "<<tot_jets<<endl;
505 
506  return 0;
507  }
508 
509