macros/leptoquarks_analysis/likelihood.C

0001 #include "tau_commons.h"
0002
0003 /**
0004  * Compare reconstructed observables between tau jets and quark jets
0005  *
0006  * Written by nils.feege@stonybrook.edu & sean.jeffas@stonybrook.edu
0007  */
0008
0009
0010
0011 double Average(vector<double> v)
0012 {
0013   double sum=0;
0014   for(int i=0;(unsigned)i<v.size();i++)
0015     sum+=v[i];
0016   return sum/v.size();
0017 }
0018
0019
0020 int likelihood()
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[4] = {"3pion","SM","DIScharged","DISneutral"};
0029   string seed[10] = {"1","2","3","4","5","6","7","8","9","10"};
0030   string R_max = "5";
0031
0032   /* open inout files and merge trees */
0033   TChain chain("event");
0034   TChain chain_tau("event");
0035   TChain t("event");
0036
0037   // File for event topology of ID tau events
0038   ofstream myfile;
0039   string filename = "./data/Event_topology.csv";
0040   myfile.open(filename.c_str());
0041
0042   // Loop over geant files and split LQ and SM events into two TChains used for jet characterization //
0043   for (int b = 0; b<2; b++){
0044     for (int g = 0; g<10; g++){
0045       // Skip file that is used to test //
0046       if(b == 1 && g==0) continue;
0047       string file = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/data/LeptoAna_p250_e20_1000events_"+seed[g]+"seed_"+type[b]+"_r0"+R_max+".root";
0048
0049       if(b==0) chain_tau.Add(file.c_str());
0050       else chain.Add(file.c_str());
0051     }
0052   }
0053
0054   for(int b = 1; b<2; b++){
0055     for(int g=0; g<1; g++){
0056       // File to test identification //
0057       string inputFile = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/data/LeptoAna_p250_e20_1000events_"+seed[g]+"seed_"+type[b]+"_r05.root";
0058
0059       t.Add(inputFile.c_str());
0060     }
0061   }
0062
0063
0064
0065   /* Create temporary canvas */
0066   TCanvas* ctemp = new TCanvas();
0067
0068   vector< TString > observables;
0069   vector< TString > observables_name;
0070
0071   vector< float > plots_ymax;
0072   vector< float > plots_xmin;
0073   vector< float > plots_xmax;
0074
0075   /* R_max */
0076   observables.push_back( "tracks_rmax_r04" );
0077   observables_name.push_back( "R_{track}^{max}" );
0078   plots_ymax.push_back(0.06);
0079   plots_xmin.push_back(0);
0080   plots_xmax.push_back(0.4);
0081
0082   /* Number of tracks */
0083   observables.push_back( "tracks_count_r04" );
0084   observables_name.push_back( "N_{tracks}" );
0085   plots_ymax.push_back(0.9);
0086   plots_xmin.push_back(0);
0087   plots_xmax.push_back(10);
0088
0089   /* Charge sum from tracks */
0090   observables.push_back( "tracks_chargesum_r04" );
0091   observables_name.push_back( "#Sigma q_{tracks}" );
0092   plots_ymax.push_back(0.9);
0093   plots_xmin.push_back(-5);
0094   plots_xmax.push_back(5);
0095
0096   /* Jet radius */
0097   observables.push_back( "jetshape_radius" );
0098   observables_name.push_back( "R_{jet}" );
0099   plots_ymax.push_back(0.08);
0100   plots_xmin.push_back(0);
0101   plots_xmax.push_back(0.5);
0102
0103   /* Jetshape */
0104   observables.push_back( "jetshape_econe_r01 / jetshape_econe_r05" );
0105   observables_name.push_back( "E_{cone}^{R<0.1} / E_{cone}^{R<0.5}" );
0106   plots_ymax.push_back(0.08);
0107   plots_xmin.push_back(0);
0108   plots_xmax.push_back(1);
0109
0110   /* Jetshape */
0111   observables.push_back( "jetshape_econe_r02 / jetshape_econe_r05" );
0112   observables_name.push_back( "E_{cone}^{R<0.2} / E_{cone}^{R<0.5}" );
0113   plots_ymax.push_back(0.08);
0114   plots_xmin.push_back(0);
0115   plots_xmax.push_back(1);
0116
0117   /* Jet eta */
0118   observables.push_back( "jet_eta" );
0119   observables_name.push_back( "#eta_{jet}" );
0120   plots_ymax.push_back(0.1);
0121   plots_xmin.push_back(-1);
0122   plots_xmax.push_back(3);
0123
0124   /* Jet mass */
0125   observables.push_back( "jet_minv" );
0126   observables_name.push_back( "Mass_{jet} (GeV)" );
0127   plots_ymax.push_back(0.1);
0128   plots_xmin.push_back(0);
0129   plots_xmax.push_back(10);
0130
0131   /* Jet energy */
0132   observables.push_back( "jet_etotal" );
0133   observables_name.push_back( "E_{jet} (GeV)" );
0134   plots_ymax.push_back(0.08);
0135   plots_xmin.push_back(0);
0136   plots_xmax.push_back(70);
0137
0138   /* Jet vertex */
0139   observables.push_back( "tracks_vertex" );
0140   observables_name.push_back( "vertex distance (cm)" );
0141   plots_ymax.push_back(0.6);
0142   plots_xmin.push_back(0);
0143   plots_xmax.push_back(1);
0144
0145
0146   /* Plot distributions */
0147   TString name_uds_base("h_uds_");
0148   TString name_tau_base("h_tau_");
0149
0150   const int nvar = observables.size();
0151
0152   TH1F* h_uds[nvar];
0153   TH1F* h_tau[nvar];
0154
0155   TString name_uds_i = name_uds_base;
0156   TString name_tau_i = name_tau_base;
0157
0158   /* create histograms */
0159   for ( unsigned j = 0; j < observables.size(); j++ ){
0160     name_uds_i = name_uds_base;
0161     name_tau_i = name_tau_base;
0162
0163     name_uds_i.Append(j);
0164     name_tau_i.Append(j);
0165
0166     h_uds[j] = new TH1F( name_uds_i, "", 100, plots_xmin.at(j), plots_xmax.at(j));
0167     h_tau[j] = new TH1F( name_tau_i, "", 100, plots_xmin.at(j), plots_xmax.at(j));
0168   }
0169
0170   /* Fill histograms*/
0171   for ( unsigned i = 0; i < observables.size(); i++ )
0172     {
0173       //cout << "Plotting " << observables.at(i) << endl;
0174
0175       ctemp->cd();
0176
0177       name_uds_i = name_uds_base;
0178       name_tau_i = name_tau_base;
0179
0180       name_uds_i.Append(i);
0181       name_tau_i.Append(i);
0182
0183       // SM parton jet observable histograms //
0184       chain.Draw( observables.at(i) + " >> " + name_uds_i, tau_commons::select_true_uds && tau_commons::select_accept_jet, "goff");
0185       h_uds[i]->Scale(1./h_uds[i]->Integral());
0186       h_uds[i]->GetXaxis()->SetTitle( observables_name.at(i) );
0187       h_uds[i]->GetYaxis()->SetRangeUser(0, plots_ymax.at(i) );
0188       h_uds[i]->SetLineColor(col1);
0189       h_uds[i]->SetFillColor(col1);
0190       h_uds[i]->SetFillStyle(1001);
0191       h_uds[i]->GetYaxis()->SetTitle("# entries / #Sigma entries");
0192       h_uds[i]->GetYaxis()->SetTitleOffset(1.5);
0193
0194       // LQ tau jet observable histograms //
0195       chain_tau.Draw( observables.at(i) + " >> " + name_tau_i, tau_commons::select_true_tau && tau_commons::select_accept_jet, "goff");
0196       h_tau[i]->Scale(1./h_tau[i]->Integral());
0197       h_tau[i]->SetLineColor(col2);
0198       h_tau[i]->SetFillColor(col2);
0199       h_tau[i]->SetFillStyle(0);
0200
0201
0202       /* create Canvas and draw histograms */
0203       TCanvas *c1 = new TCanvas();
0204
0205       h_uds[i]->DrawClone("");
0206       h_tau[i]->DrawClone("sames");
0207       gPad->RedrawAxis();
0208
0209       TLegend *legend = new TLegend(0.4,0.6,0.7,0.89);
0210       legend->AddEntry(h_uds[i],"SM parton jet","l");
0211       legend->AddEntry(h_tau[i],"LQ #tau jet","l");
0212       legend->Draw();
0213
0214       //      c1->Print( Form( "plots/compare_observables_%d.eps", i ) );
0215       //   c1->Print( Form( "plots/compare_observables_%d.png", i ) );
0216     }
0217
0218   // Define jet and //
0219   vector<float> * tracks_rmax;
0220   vector<float> * tracks_count;
0221   vector<float> * tracks_chargesum;
0222   vector<float> * tracks_vertex;
0223   vector<float> * jetshape_radius;
0224   vector<float> * jetshape_econe_1;
0225   vector<float> * jetshape_econe_2;
0226   vector<float> * jetshape_econe_5;
0227   vector<float> * jet_eta;
0228   vector<float> * jet_phi;
0229   vector<float> * jet_minv;
0230   vector<float> * jet_etotal;
0231   vector<float> * jet_ptrans;
0232   vector<float> * evtgen_pid;
0233
0234   // Define event variables //
0235   float ptmiss;
0236   float miss_phi;
0237   float is_lq;
0238
0239   // See variables from tree //
0240   t.SetBranchAddress("tracks_rmax_R",&tracks_rmax);
0241   t.SetBranchAddress("tracks_count_R",&tracks_count);
0242   t.SetBranchAddress("tracks_chargesum_R",&tracks_chargesum);
0243   t.SetBranchAddress("tracks_vertex",&tracks_vertex);
0244   t.SetBranchAddress("jetshape_radius",&jetshape_radius);
0245   t.SetBranchAddress("jetshape_econe_r01",&jetshape_econe_1);
0246   t.SetBranchAddress("jetshape_econe_r02",&jetshape_econe_2);
0247   t.SetBranchAddress("jetshape_econe_r05",&jetshape_econe_5);
0248   t.SetBranchAddress("jet_eta",&jet_eta);
0249   t.SetBranchAddress("jet_phi",&jet_phi);
0250   t.SetBranchAddress("jet_minv",&jet_minv);
0251   t.SetBranchAddress("jet_etotal",&jet_etotal);
0252   t.SetBranchAddress("evtgen_pid",&evtgen_pid);
0253   t.SetBranchAddress("jet_ptrans",&jet_ptrans);
0254   t.SetBranchAddress("Et_miss",&ptmiss);
0255   t.SetBranchAddress("Et_miss_phi",&miss_phi);
0256   t.SetBranchAddress("is_lq_event",&is_lq);
0257
0258   /* identification counts */
0259   int tau_as_tau = 0;
0260   int tau_as_DIS = 0;
0261   int DIS_as_DIS = 0;
0262   int DIS_as_tau = 0;
0263   int id_e = 0;
0264
0265   int tot_jets = 0;
0266
0267   // Define likelihood cut for tau jet identification //
0268   double lik_cut = 0.71;
0269
0270   int n_tot = t.GetEntries();
0271
0272   // Loop over all events //
0273   for( int k = 0; k < n_tot; k++){
0274     t.GetEntry(k);
0275
0276     // Define some event variables //
0277     double delta_phi;
0278     bool is_e = false;
0279     bool is_miss = false;
0280     vector<double> jet_lik;
0281     vector<double> jet_id;
0282     vector<double> miss_jet_phi;
0283
0284     // Loop over jets in event to look for electron //
0285     for(int l = 0; l < tracks_rmax->size(); l++){
0286       int count = (*tracks_count)[l];
0287       double rmax = (*tracks_rmax)[l];
0288       int chargesum = (*tracks_chargesum)[l];
0289       double radius = (*jetshape_radius)[l];
0290       double eta = (*jet_eta)[l];
0291       double phi = (*jet_phi)[l];
0292       double etotal = (*jet_etotal)[l];
0293       double vertex = (*tracks_vertex)[l];
0294       double ptrans = (*jet_ptrans)[l];
0295       int pid = (*evtgen_pid)[l];
0296
0297       // If any jets pass these cuts then this is considered an electron event //
0298       if(count == 1 && chargesum == -1 && ptrans > 5 && vertex < 0.03 && radius < 0.05) is_e = true;
0299
0300     }
0301
0302     // Loop over all jets in events //
0303     for(int l = 0; l < tracks_rmax->size(); l++){
0304
0305       // Likelihood variable //
0306       vector<double> likelihood;
0307
0308       // Get jet observables //
0309       double rmax = (*tracks_rmax)[l];
0310       int count = (*tracks_count)[l];
0311       int chargesum = (*tracks_chargesum)[l];
0312       double vertex = (*tracks_vertex)[l];
0313       double radius = (*jetshape_radius)[l];
0314       double econe_1 = (*jetshape_econe_1)[l];
0315       double econe_2 = (*jetshape_econe_2)[l];
0316       double econe_5 = (*jetshape_econe_5)[l];
0317       double eta = (*jet_eta)[l];
0318       double phi = (*jet_phi)[l];
0319       double minv = (*jet_minv)[l];
0320       double etotal = (*jet_etotal)[l];
0321       double ptrans = (*jet_ptrans)[l];
0322       int pid = (*evtgen_pid)[l];
0323       tot_jets++;
0324
0325       // If electron event then ID as parton and continue //
0326       if(is_e && pid == 15) tau_as_DIS++;
0327       if(is_e && pid != 15) DIS_as_DIS++;
0328       if(is_e) continue;
0329
0330
0331       // x bin number //
0332       int x_bin;
0333
0334       // Find where jet observable lies on characteristic histograms and use histogram heights to find likelihood //
0335       // Currently count, chargesum, eta, and vertex give best results. Rest are left commented out //
0336
0337       //x_bin = h_tau[0]->GetXaxis()->FindBin(rmax);
0338       //likelihood.push_back( h_tau[0]->GetBinContent(x_bin) / (h_tau[0]->GetBinContent(x_bin) + h_uds[0]->GetBinContent(x_bin)));
0339
0340       x_bin = h_tau[1]->GetXaxis()->FindBin(count);
0341       likelihood.push_back( h_tau[1]->GetBinContent(x_bin) / (h_tau[1]->GetBinContent(x_bin) + h_uds[1]->GetBinContent(x_bin)));
0342
0343       x_bin = h_tau[2]->GetXaxis()->FindBin(chargesum);
0344       likelihood.push_back( h_tau[2]->GetBinContent(x_bin) / (h_tau[2]->GetBinContent(x_bin) + h_uds[2]->GetBinContent(x_bin)));
0345
0346       //x_bin = h_tau[3]->GetXaxis()->FindBin(radius);
0347       //likelihood.push_back( h_tau[3]->GetBinContent(x_bin) / (h_tau[3]->GetBinContent(x_bin) + h_uds[3]->GetBinContent(x_bin)));
0348       /*
0349       x_bin = h_tau[4]->GetXaxis()->FindBin(econe_1/econe_5);
0350       likelihood.push_back( h_tau[4]->GetBinContent(x_bin) / (h_tau[4]->GetBinContent(x_bin) + h_uds[4]->GetBinContent(x_bin)));
0351
0352       x_bin = h_tau[5]->GetXaxis()->FindBin(econe_2/econe_5);
0353       likelihood.push_back( h_tau[5]->GetBinContent(x_bin) / (h_tau[5]->GetBinContent(x_bin) + h_uds[5]->GetBinContent(x_bin)));
0354       */
0355
0356       x_bin = h_tau[6]->GetXaxis()->FindBin(eta);
0357       likelihood.push_back( h_tau[6]->GetBinContent(x_bin) / (h_tau[6]->GetBinContent(x_bin) + h_uds[6]->GetBinContent(x_bin)));
0358
0359       //x_bin = h_tau[7]->GetXaxis()->FindBin(minv);
0360       //likelihood.push_back( h_tau[7]->GetBinContent(x_bin) / (h_tau[7]->GetBinContent(x_bin) + h_uds[7]->GetBinContent(x_bin)));
0361
0362       //x_bin = h_tau[8]->GetXaxis()->FindBin(etotal);
0363       //likelihood.push_back( h_tau[8]->GetBinContent(x_bin) / (h_tau[8]->GetBinContent(x_bin) + h_uds[8]->GetBinContent(x_bin)));
0364
0365       // Exclude jets where vertex = inf which causes problems //
0366       x_bin = h_tau[9]->GetXaxis()->FindBin(vertex);
0367       if(vertex == vertex) likelihood.push_back( h_tau[9]->GetBinContent(x_bin) / (h_tau[9]->GetBinContent(x_bin) + h_uds[9]->GetBinContent(x_bin)));
0368
0369
0370       // Define total likelihood as average of likelihoods of observables //
0371       double avg = Average(likelihood);
0372
0373       // Make sure final likelihood is real finite number //
0374       if(avg != avg) continue;
0375
0376       // Only use jets with transverse momentum > 5 GeV //
0377       // Write jet likelihoods to array for event //
0378       if(ptrans < 5) avg = 0;
0379       if(ptrans > 5) {
0380     jet_lik.push_back(avg);
0381     jet_id.push_back(pid);
0382     miss_jet_phi.push_back(phi);
0383       }
0384
0385       // If pt < 5 call parton jet //
0386       if(ptrans < 5 && pid == 15) tau_as_DIS++;
0387       if(ptrans < 5 && pid != 15) DIS_as_DIS++;
0388
0389       // This block uses hard cuts on variables for identification //
0390       /*
0391       if(count == 3 && chargesum == -1 && ptrans > 5 && vertex > 0.028 && eta < 0.82)) {
0392     if(pid == 15) tau_as_tau++;
0393     else DIS_as_tau++;
0394       }
0395       else {
0396     if(pid == 15) tau_as_DIS++;
0397     else DIS_as_DIS++;
0398       }
0399     */
0400
0401     }
0402
0403     // Increment if electron is in event //
0404     if(is_e) id_e++;
0405
0406     // This block is used for likelihood identification //
0407     // Define some variables for jet selection //
0408     int a = -1;
0409     int b = -1;
0410     double max_lik = 0;
0411     int max_lik_i;
0412     bool is_tau = false;
0413
0414     // If no jets with pt > 5 then skip //
0415     if(jet_lik.size() == 0) continue;
0416
0417     // Loop over all jets in event with pt > 5 //
0418     for(int m = 0; m < jet_lik.size(); m++){
0419       // Find if there is a tau in this event //
0420       if(jet_id[m] == 15) is_tau = true;
0421       // Find jet with maximum tau likelihood //
0422       if(jet_lik[m] > max_lik) {
0423     max_lik = jet_lik[m];
0424     max_lik_i = m;
0425       }
0426       // There is a problem with same jets being tagged with two different true values
0427       // Here we check if one of the jets is repeated //
0428       // Loop over all jets to check if two have exact same likelihood //
0429       for(int n = 0; n < jet_lik.size(); n++){
0430     // If different jets have same likelihood then set a and b to the jet indices //
0431     if(m!=n && jet_lik[m]==jet_lik[n]) {
0432       a = m;
0433       b = n;
0434       // If either is tau jet then set them both to a tau jet id //
0435       if(jet_id[m] == 15 || jet_id[n] 15) {
0436         jet_id[m] == 15;
0437         jet_id[n] ==15;
0438       }
0439       // Otherwise set them both to 0 id //
0440       else {
0441         jet_id[m] = 0;
0442         jet_id[n] = 0;
0443       }
0444     }
0445       }
0446     }
0447     // End loop over all jets with pt > 5 //
0448
0449     // Check if maximum likelihood jet passes cut then ID as tau //
0450     if(jet_lik[max_lik_i] > lik_cut){
0451       // Find ID tau jet as true tau //
0452       if(jet_id[max_lik_i] == 15) tau_as_tau++;
0453       // Otherwise find ID tau jet as true parton //
0454       else DIS_as_tau++;
0455       // All other jets are ID as parton, minus one that was ID as tau //
0456       DIS_as_DIS = DIS_as_DIS + jet_lik.size() - 1;
0457       // Set event topology to be recorded for later study //
0458       is_miss = true;
0459     }
0460     // If fails likelihood cut then ID as parton
0461     else{
0462       // If a tau was in the event then it is ID as parton with all other jets //
0463       if(is_tau) tau_as_DIS++;
0464       if(is_tau) DIS_as_DIS = DIS_as_DIS + jet_lik.size() - 1;
0465       // If no tau then every jet is parton and ID correctly //
0466       if(!is_tau) DIS_as_DIS = DIS_as_DIS + jet_lik.size();
0467     }
0468
0469     // Check if duplicate jets were found and subtract one from ID jets if it was //
0470     if(a > -1 && b > -1){
0471       if(jet_id[a] == 15 && jet_lik[max_lik_i] > lik_cut) DIS_as_tau--;
0472       if(jet_id[a] != 15 && jet_lik[max_lik_i] < lik_cut) DIS_as_DIS--;
0473     }
0474
0475     // If jet is ID as tau then calculate ptmiss and acoplanarity for this event //
0476     if(is_miss){
0477       delta_phi = fabs(miss_phi - miss_jet_phi[max_lik_i]) * 180 / TMath::Pi();
0478       if(delta_phi > 180) delta_phi = 360 - delta_phi;
0479       // Write variables with true LQ tagging to file for later study //
0480       myfile<<ptmiss<<","<<delta_phi<<","<<is_lq<<endl;
0481     }
0482   }
0483
0484   // Define efficiency and false positive rate //
0485   double efficiency = tau_as_tau*1.0/(n_tot*1.0);
0486   double FP;
0487
0488   if(tau_as_tau == 0 && DIS_as_tau == 0) FP = 0;
0489   else FP = DIS_as_tau*1.0/((tau_as_tau+DIS_as_tau)*1.0);
0490
0491
0492   // Confusion Matrix //
0493   cout<<"         | Act Tau   Act Parton"<<endl;
0494   cout<<"------------------------------"<<endl;
0495   cout<<"ID Tau   | "<<tau_as_tau<<"         "<<DIS_as_tau<<endl;
0496   cout<<"ID Parton| "<<tau_as_DIS<<"         "<<DIS_as_DIS<<endl;
0497
0498   cout<<endl;
0499
0500   // Write out some important variables //
0501   cout<<"Tau Recovery Rate: "<<efficiency<<endl;
0502   cout<<"False Positive Rate: "<<FP<<endl;
0503   cout<<"Identified Electron Events : "<<id_e<<endl;
0504   cout<<"Total Jets: "<<tot_jets<<endl;
0505
0506   return 0;
0507   }
0508
0509