Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
likelihood_topology.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file likelihood_topology.C
1 #include "tau_commons.h"
2 
10 double Average(vector<double> v)
11 {
12  double sum=0;
13  for(int i=0;(unsigned)i<v.size();i++)
14  sum+=v[i];
15  return sum/v.size();
16 }
17 
18 
19 int likelihood_topology()
20 {
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[2] = {"3pion","DIScharged"};
29  string seed[10] = {"1","2","3","4","5","6","7","8","9","10"};
30  string R_max = "5";
31  string nevents;
32  string file;
33 
34  /* open inout files and merge trees */
35  TChain chain_CC("event");
36  TChain chain_tau("event");
37 
38 
39  for (int b = 0; b<2; b++){
40  for (int g = 0; g<10; g++){
41  // Set some labels //
42  if(b==0) nevents = "1000";
43  if(b==1 || b==2) nevents = "100";
44 
45  if(b==0) file = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/data/LeptoAna_p250_e20_"+nevents+"events_"+seed[g]+"seed_"+type[b]+"_r0"+R_max+".root";
46  if(b==1) file = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/data/LeptoAna_p250_e20_"+nevents+"events_"+seed[g]+"seed_"+type[b]+"_tau_r0"+R_max+".root";
47 
48  if(b==0) chain_tau.Add(file.c_str());
49  if(b==1) chain_CC.Add(file.c_str());
50  }
51  }
52 
53  // Add file with selected tau event topologies //
54  string inputFile = "./data/Event_topology.csv";
55 
56  TTree *t=new TTree("MyTree","MyTree");
57  t->ReadFile(inputFile.c_str(),"ptmiss:acoplan:is_lq");
58 
59  // Create temporary canvas
60  TCanvas* ctemp = new TCanvas();
61 
62  vector< TString > observables;
63  vector< TString > observables_name;
64 
65  vector< float > plots_xmin;
66  vector< float > plots_xmax;
67 
68  // Missing energy
69  observables.push_back( "Et_miss" );
70  observables_name.push_back( "p_{T}^{miss} (GeV)" );
71  plots_xmin.push_back(0);
72  plots_xmax.push_back(50);
73 
74  // Acoplanarity
75  observables.push_back( "Et_miss_phi" );
76  observables_name.push_back( "#Delta#phi_{miss - #tau jet}" );
77  plots_xmin.push_back(0);
78  plots_xmax.push_back(180);
79 
80  TString name_CC_base("h_CC_");
81  TString name_tau_base("h_tau_");
82 
83  const int nvar = observables.size();
84 
85  TH1F* h_CC[nvar];
86  TH1F* h_tau[nvar];
87 
88  TString name_CC_i = name_CC_base;
89  TString name_tau_i = name_tau_base;
90 
91  // create histograms
92  for ( int l = 0; l < observables.size(); l++ ){
93  name_CC_i = name_CC_base;
94  name_tau_i = name_tau_base;
95 
96  name_CC_i.Append(l);
97  name_tau_i.Append(l);
98 
99  h_CC[l] = new TH1F( name_CC_i, "", 50, plots_xmin.at(l), plots_xmax.at(l));
100  h_tau[l] = new TH1F( name_tau_i, "", 50, plots_xmin.at(l), plots_xmax.at(l));
101  }
102 
103  // Fill Histograms //
104  for(int i = 0;i < observables.size(); i++){
105 
106  ctemp->cd();
107 
108  name_CC_i = name_CC_base;
109  name_tau_i = name_tau_base;
110 
111  name_CC_i.Append(i);
112  name_tau_i.Append(i);
113 
114  // Fill Et_miss straight from tree //
115  if(i==0){
116  // Fill CC histogram //
117  int n_CC = chain_CC.Draw( observables.at(i) + " >> " + name_CC_i, "", "goff");
118  h_CC[i]->Scale(1./h_CC[i]->Integral());
119  h_CC[i]->GetXaxis()->SetTitle( observables_name.at(i) );
120  h_CC[i]->GetYaxis()->SetRangeUser(0, .13 );
121  h_CC[i]->SetLineColor(col1);
122  h_CC[i]->SetFillColor(col1);
123  h_CC[i]->SetFillStyle(0);
124  h_CC[i]->GetYaxis()->SetTitle("# entries / #Sigma entries");
125  h_CC[i]->GetYaxis()->SetTitleOffset(1.5);
126 
127  // Fill LQ histogram //
128  int n_tau = chain_tau.Draw( observables.at(i) + " >> " + name_tau_i, "", "goff");
129  h_tau[i]->Scale(1./h_tau[i]->Integral());
130  h_tau[i]->SetLineColor(col2);
131  h_tau[i]->SetFillColor(col2);
132  h_tau[i]->SetFillStyle(0);
133 
134  TCanvas *c1 = new TCanvas();
135 
136  h_CC[i]->DrawClone("");
137  h_tau[i]->DrawClone("sames");
138  gPad->RedrawAxis();
139 
140  TLegend *legend = new TLegend(0.4,0.6,0.7,0.89);
141  legend->AddEntry(h_CC[i],"Charged Current","l");
142  legend->AddEntry(h_tau[i],"Leptoquark","l");
143  legend->Draw();
144  }
145 
146  // Calculate acoplanarity and then fill histogram //
147  if(i==1){
148 
149  float Et_miss_phi_tau, Et_miss_phi_CC;
150  float tau_phi_tau, tau_phi_CC;
151 
152  // Read tau phi and missing momentum phi from tree //
153  chain_CC.SetBranchAddress("Et_miss_phi",&Et_miss_phi_CC);
154  chain_CC.SetBranchAddress("reco_tau_phi",&tau_phi_CC);
155 
156  chain_tau.SetBranchAddress("Et_miss_phi",&Et_miss_phi_tau);
157  chain_tau.SetBranchAddress("reco_tau_phi",&tau_phi_tau);
158 
159  // Phi angle between reco tau and miss mom //
160  double delta_phi;
161 
162  // Calculate delta phi for each event type //
163  for(int j = 0; j<n_CC;j++){
164  chain_CC.GetEntry(j);
165  delta_phi = fabs(Et_miss_phi_CC - tau_phi_CC) * 180 / TMath::Pi();
166  if(delta_phi > 180) delta_phi = 360 - delta_phi;
167  h_CC[i]->Fill(delta_phi);
168  }
169 
170  for(int k = 0; k<n_tau;k++){
171  chain_tau.GetEntry(k);
172  delta_phi = fabs(Et_miss_phi_tau - tau_phi_tau) * 180 / TMath::Pi();
173  if(delta_phi > 180) delta_phi = 360 - delta_phi;
174  h_tau[i]->Fill(delta_phi);
175  }
176 
177  //Fill histograms //
178  h_CC[i]->Scale(1./h_CC[i]->Integral());
179  h_CC[i]->GetXaxis()->SetTitle( observables_name.at(i) );
180  h_CC[i]->GetYaxis()->SetRangeUser(0, 0.27 );
181  h_CC[i]->SetLineColor(col1);
182  h_CC[i]->SetFillColor(col1);
183  h_CC[i]->SetFillStyle(0);
184  h_CC[i]->GetYaxis()->SetTitle("# entries / #Sigma entries");
185  h_CC[i]->GetYaxis()->SetTitleOffset(1.5);
186 
187  h_tau[i]->Scale(1./h_tau[i]->Integral());
188  h_tau[i]->SetLineColor(col2);
189  h_tau[i]->SetFillColor(col2);
190  h_tau[i]->SetFillStyle(0);
191 
192  TCanvas *c1 = new TCanvas();
193 
194  h_CC[i]->DrawClone("");
195  h_tau[i]->DrawClone("sames");
196  gPad->RedrawAxis();
197 
198  TLegend *legend = new TLegend(0.4,0.6,0.7,0.89);
199  legend->AddEntry(h_CC[i],"Charged Current","l");
200  legend->AddEntry(h_tau[i],"Leptoquark","l");
201  legend->Draw();
202 
203  }
204 
205  }
206 
207 
208  float acoplan, ptmiss, is_lq;
209 
210  // Get variables from tree //
211  t->SetBranchAddress("acoplan",&acoplan);
212  t->SetBranchAddress("ptmiss",&ptmiss);
213  t->SetBranchAddress("is_lq",&is_lq);
214 
215  // identification counts //
216  int tau_as_tau = 0;
217  int tau_as_DIS = 0;
218  int DIS_as_DIS = 0;
219  int DIS_as_tau = 0;
220 
221  // Define likelihood cut //
222  double lik_cut = 0.60;
223 
224  int n_tot = t->GetEntries();
225 
226  // Loop over all events //
227  for( int k = 0; k < n_tot; k++){
228  t->GetEntry(k);
229 
230  vector<double> likelihood;
231 
232  int x_bin;
233 
234  // Find likelihood for each observable //
235  x_bin = h_tau[0]->GetXaxis()->FindBin(ptmiss);
236  likelihood.push_back( h_tau[0]->GetBinContent(x_bin) / (h_tau[0]->GetBinContent(x_bin) + h_CC[0]->GetBinContent(x_bin)));
237 
238  x_bin = h_tau[1]->GetXaxis()->FindBin(acoplan);
239  likelihood.push_back( h_tau[1]->GetBinContent(x_bin) / (h_tau[1]->GetBinContent(x_bin) + h_CC[1]->GetBinContent(x_bin)));
240 
241  // Average the likelihoods //
242  double avg = Average(likelihood);
243 
244  // Make sure likelihood is real number //
245  if(avg != avg) continue;
246 
247  // Use likelihood cut for identification //
248  if(avg>lik_cut){
249  // Call event LQ and check if true //
250  if(is_lq == 1) tau_as_tau++;
251  else DIS_as_tau++;
252  }
253  else{
254  // Call event SM and check if true //
255  if(is_lq == 1) tau_as_DIS++;
256  else DIS_as_DIS++;
257  }
258 
259 
260  }
261 
262 
263  // Define efficiency and false positive rate //
264  double efficiency = tau_as_tau*1.0/(n_tot*1.0);
265  double FP;
266 
267  if(tau_as_tau == 0 && DIS_as_tau == 0) FP = 0;
268  else FP = DIS_as_tau*1.0/((tau_as_tau+DIS_as_tau)*1.0);
269 
270  // Confusion Matrix //
271  cout<<" | ID SM ID LQ"<<endl;
272  cout<<"---------------------"<<endl;
273  cout<<"Act SM| "<<DIS_as_DIS<<" "<<DIS_as_tau<<endl;
274  cout<<"Act LQ| "<<tau_as_DIS<<" "<<tau_as_tau<<endl;
275 
276  cout<<endl;
277 
278  cout<<"LQ Recovery Rate: "<<efficiency<<endl;
279  cout<<"False Positive Rate: "<<FP<<endl;
280 
281 
282 
283  return 0;
284 }