macros/leptoquarks_analysis/likelihood_topology.C

0001 #include "tau_commons.h"
0002
0003 /**
0004  * Plot missing trasnverse energy and other event-wise observables
0005  *
0006  * Written by nils.feege@stonybrook.edu & sean.jeffas@stonybrook.edu
0007  */
0008
0009
0010 double Average(vector<double> v)
0011 {
0012   double sum=0;
0013   for(int i=0;(unsigned)i<v.size();i++)
0014     sum+=v[i];
0015   return sum/v.size();
0016 }
0017
0018
0019 int likelihood_topology()
0020 {
0021
0022   gStyle->SetOptStat(0);
0023
0024   unsigned col1 = kOrange+7;
0025   unsigned col2 = kBlue+2;
0026
0027   // labels for geant files //
0028   string type[2] = {"3pion","DIScharged"};
0029   string seed[10] = {"1","2","3","4","5","6","7","8","9","10"};
0030   string R_max = "5";
0031   string nevents;
0032   string file;
0033
0034   /* open inout files and merge trees */
0035   TChain chain_CC("event");
0036   TChain chain_tau("event");
0037
0038
0039   for (int b = 0; b<2; b++){
0040     for (int g = 0; g<10; g++){
0041       // Set some labels //
0042       if(b==0) nevents = "1000";
0043       if(b==1 || b==2) nevents = "100";
0044
0045       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";
0046       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";
0047
0048       if(b==0) chain_tau.Add(file.c_str());
0049       if(b==1) chain_CC.Add(file.c_str());
0050     }
0051   }
0052
0053   // Add file with selected tau event topologies //
0054   string inputFile = "./data/Event_topology.csv";
0055
0056   TTree *t=new TTree("MyTree","MyTree");
0057   t->ReadFile(inputFile.c_str(),"ptmiss:acoplan:is_lq");
0058
0059   // Create temporary canvas
0060   TCanvas* ctemp = new TCanvas();
0061
0062   vector< TString > observables;
0063   vector< TString > observables_name;
0064
0065   vector< float > plots_xmin;
0066   vector< float > plots_xmax;
0067
0068   // Missing energy
0069   observables.push_back( "Et_miss" );
0070   observables_name.push_back( "p_{T}^{miss} (GeV)" );
0071   plots_xmin.push_back(0);
0072   plots_xmax.push_back(50);
0073
0074   // Acoplanarity
0075   observables.push_back( "Et_miss_phi" );
0076   observables_name.push_back( "#Delta#phi_{miss - #tau jet}" );
0077   plots_xmin.push_back(0);
0078   plots_xmax.push_back(180);
0079
0080   TString name_CC_base("h_CC_");
0081   TString name_tau_base("h_tau_");
0082
0083   const int nvar = observables.size();
0084
0085   TH1F* h_CC[nvar];
0086   TH1F* h_tau[nvar];
0087
0088   TString name_CC_i = name_CC_base;
0089   TString name_tau_i = name_tau_base;
0090
0091   // create histograms
0092   for ( int l = 0; l < observables.size(); l++ ){
0093     name_CC_i = name_CC_base;
0094     name_tau_i = name_tau_base;
0095
0096     name_CC_i.Append(l);
0097     name_tau_i.Append(l);
0098
0099     h_CC[l] = new TH1F( name_CC_i, "", 50, plots_xmin.at(l), plots_xmax.at(l));
0100     h_tau[l] = new TH1F( name_tau_i, "", 50, plots_xmin.at(l), plots_xmax.at(l));
0101   }
0102
0103   // Fill Histograms //
0104   for(int i = 0;i < observables.size(); i++){
0105
0106     ctemp->cd();
0107
0108     name_CC_i = name_CC_base;
0109     name_tau_i = name_tau_base;
0110
0111     name_CC_i.Append(i);
0112     name_tau_i.Append(i);
0113
0114     // Fill Et_miss straight from tree //
0115     if(i==0){
0116       // Fill CC histogram //
0117       int n_CC = chain_CC.Draw( observables.at(i) + " >> " + name_CC_i, "", "goff");
0118       h_CC[i]->Scale(1./h_CC[i]->Integral());
0119       h_CC[i]->GetXaxis()->SetTitle( observables_name.at(i) );
0120       h_CC[i]->GetYaxis()->SetRangeUser(0, .13 );
0121       h_CC[i]->SetLineColor(col1);
0122       h_CC[i]->SetFillColor(col1);
0123       h_CC[i]->SetFillStyle(0);
0124       h_CC[i]->GetYaxis()->SetTitle("# entries / #Sigma entries");
0125       h_CC[i]->GetYaxis()->SetTitleOffset(1.5);
0126
0127       // Fill LQ histogram //
0128       int n_tau = chain_tau.Draw( observables.at(i) + " >> " + name_tau_i, "", "goff");
0129       h_tau[i]->Scale(1./h_tau[i]->Integral());
0130       h_tau[i]->SetLineColor(col2);
0131       h_tau[i]->SetFillColor(col2);
0132       h_tau[i]->SetFillStyle(0);
0133
0134       TCanvas *c1 = new TCanvas();
0135
0136       h_CC[i]->DrawClone("");
0137       h_tau[i]->DrawClone("sames");
0138       gPad->RedrawAxis();
0139
0140       TLegend *legend = new TLegend(0.4,0.6,0.7,0.89);
0141       legend->AddEntry(h_CC[i],"Charged Current","l");
0142       legend->AddEntry(h_tau[i],"Leptoquark","l");
0143       legend->Draw();
0144     }
0145
0146     // Calculate acoplanarity and then fill histogram //
0147     if(i==1){
0148
0149       float Et_miss_phi_tau, Et_miss_phi_CC;
0150       float tau_phi_tau, tau_phi_CC;
0151
0152       // Read tau phi and missing momentum phi from tree //
0153       chain_CC.SetBranchAddress("Et_miss_phi",&Et_miss_phi_CC);
0154       chain_CC.SetBranchAddress("reco_tau_phi",&tau_phi_CC);
0155
0156       chain_tau.SetBranchAddress("Et_miss_phi",&Et_miss_phi_tau);
0157       chain_tau.SetBranchAddress("reco_tau_phi",&tau_phi_tau);
0158
0159       // Phi angle between reco tau and miss mom   //
0160       double delta_phi;
0161
0162       // Calculate delta phi for each event type //
0163       for(int j = 0; j<n_CC;j++){
0164     chain_CC.GetEntry(j);
0165     delta_phi = fabs(Et_miss_phi_CC - tau_phi_CC) * 180 / TMath::Pi();
0166     if(delta_phi > 180) delta_phi = 360 - delta_phi;
0167     h_CC[i]->Fill(delta_phi);
0168       }
0169
0170       for(int k = 0; k<n_tau;k++){
0171         chain_tau.GetEntry(k);
0172         delta_phi = fabs(Et_miss_phi_tau - tau_phi_tau) * 180 / TMath::Pi();
0173         if(delta_phi > 180) delta_phi = 360 - delta_phi;
0174         h_tau[i]->Fill(delta_phi);
0175       }
0176
0177       //Fill histograms //
0178       h_CC[i]->Scale(1./h_CC[i]->Integral());
0179       h_CC[i]->GetXaxis()->SetTitle( observables_name.at(i) );
0180       h_CC[i]->GetYaxis()->SetRangeUser(0, 0.27 );
0181       h_CC[i]->SetLineColor(col1);
0182       h_CC[i]->SetFillColor(col1);
0183       h_CC[i]->SetFillStyle(0);
0184       h_CC[i]->GetYaxis()->SetTitle("# entries / #Sigma entries");
0185       h_CC[i]->GetYaxis()->SetTitleOffset(1.5);
0186
0187       h_tau[i]->Scale(1./h_tau[i]->Integral());
0188       h_tau[i]->SetLineColor(col2);
0189       h_tau[i]->SetFillColor(col2);
0190       h_tau[i]->SetFillStyle(0);
0191
0192       TCanvas *c1 = new TCanvas();
0193
0194       h_CC[i]->DrawClone("");
0195       h_tau[i]->DrawClone("sames");
0196       gPad->RedrawAxis();
0197
0198       TLegend *legend = new TLegend(0.4,0.6,0.7,0.89);
0199       legend->AddEntry(h_CC[i],"Charged Current","l");
0200       legend->AddEntry(h_tau[i],"Leptoquark","l");
0201       legend->Draw();
0202
0203     }
0204
0205   }
0206
0207
0208   float acoplan, ptmiss, is_lq;
0209
0210   // Get variables from tree //
0211   t->SetBranchAddress("acoplan",&acoplan);
0212   t->SetBranchAddress("ptmiss",&ptmiss);
0213   t->SetBranchAddress("is_lq",&is_lq);
0214
0215   // identification counts  //
0216   int tau_as_tau = 0;
0217   int tau_as_DIS = 0;
0218   int DIS_as_DIS = 0;
0219   int DIS_as_tau = 0;
0220
0221   // Define likelihood cut //
0222   double lik_cut = 0.60;
0223
0224   int n_tot = t->GetEntries();
0225
0226   // Loop over all events //
0227   for( int k = 0; k < n_tot; k++){
0228     t->GetEntry(k);
0229
0230     vector<double> likelihood;
0231
0232     int x_bin;
0233
0234     // Find likelihood for each observable //
0235     x_bin = h_tau[0]->GetXaxis()->FindBin(ptmiss);
0236     likelihood.push_back( h_tau[0]->GetBinContent(x_bin) / (h_tau[0]->GetBinContent(x_bin) + h_CC[0]->GetBinContent(x_bin)));
0237
0238     x_bin = h_tau[1]->GetXaxis()->FindBin(acoplan);
0239     likelihood.push_back( h_tau[1]->GetBinContent(x_bin) / (h_tau[1]->GetBinContent(x_bin) + h_CC[1]->GetBinContent(x_bin)));
0240
0241     // Average the likelihoods //
0242     double avg = Average(likelihood);
0243
0244     // Make sure likelihood is real number //
0245     if(avg != avg) continue;
0246
0247     // Use likelihood cut for identification //
0248     if(avg>lik_cut){
0249       // Call event LQ and check if true //
0250       if(is_lq == 1) tau_as_tau++;
0251       else DIS_as_tau++;
0252     }
0253     else{
0254       // Call event SM and check if true //
0255       if(is_lq == 1) tau_as_DIS++;
0256       else DIS_as_DIS++;
0257     }
0258
0259
0260   }
0261
0262
0263   // Define efficiency and false positive rate //
0264   double efficiency = tau_as_tau*1.0/(n_tot*1.0);
0265   double FP;
0266
0267   if(tau_as_tau == 0 && DIS_as_tau == 0) FP = 0;
0268   else FP = DIS_as_tau*1.0/((tau_as_tau+DIS_as_tau)*1.0);
0269
0270   // Confusion Matrix //
0271   cout<<"      | ID SM   ID LQ"<<endl;
0272   cout<<"---------------------"<<endl;
0273   cout<<"Act SM| "<<DIS_as_DIS<<"       "<<DIS_as_tau<<endl;
0274   cout<<"Act LQ| "<<tau_as_DIS<<"       "<<tau_as_tau<<endl;
0275
0276   cout<<endl;
0277
0278   cout<<"LQ Recovery Rate: "<<efficiency<<endl;
0279   cout<<"False Positive Rate: "<<FP<<endl;
0280
0281
0282
0283   return 0;
0284 }