macros/leptoquarks_analysis/event_topology_reco.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 event_topology_reco()
0020 {
0021
0022   gStyle->SetOptStat(0);
0023
0024   unsigned col1 = kOrange+7;
0025   unsigned col2 = kBlue+2;
0026   unsigned col3 = kGreen+2;
0027
0028   // labels for geant files //
0029   string type[3] = {"3pion","DIScharged","DISneutral"};
0030   string seed[10] = {"1","2","3","4","5","6","7","8","9","10"};
0031   string R_max = "5";
0032   string nevents;
0033   string file;
0034
0035   /* open inout files and merge trees */
0036   TChain chain_CC("event");
0037   TChain chain_NC("event");
0038   TChain chain_tau("event");
0039
0040
0041   for (int b = 0; b<3; b++){
0042     for (int g = 0; g<10; g++){
0043       // Skip file that is used to test identification //
0044       if(b==0 & g==0) continue;
0045       // Set some labels //
0046       if(b==0) nevents = "1000";
0047       if(b==1 || b==2) nevents = "100";
0048
0049       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";
0050       if(b==1 || b==2) file = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/data/LeptoAna_p250_e20_"+nevents+"events_"+seed[g]+"seed_"+type[b]+"_tau_r0"+R_max+".root";
0051
0052       if(b==0) chain_tau.Add(file.c_str());
0053       if(b==1) chain_CC.Add(file.c_str());
0054       if(b==2) chain_NC.Add(file.c_str());
0055     }
0056   }
0057
0058   // Add file that will be used to test //
0059   string inputFile = "/gpfs/mnt/gpfs02/phenix/scratch/spjeffas/data/LeptoAna_p250_e20_1000events_1seed_"+type[0]+"_r0"+R_max+".root";
0060
0061   TFile *f = TFile::Open(inputFile.c_str());
0062   TTree *t = (TTree*)f->Get("event");
0063
0064   // Create temporary canvas
0065   TCanvas* ctemp = new TCanvas();
0066
0067   vector< TString > observables;
0068   vector< TString > observables_name;
0069
0070   vector< float > plots_xmin;
0071   vector< float > plots_xmax;
0072
0073   // Missing energy
0074   observables.push_back( "Et_miss" );
0075   observables_name.push_back( "p_{T}^{miss} (GeV)" );
0076   plots_xmin.push_back(0);
0077   plots_xmax.push_back(50);
0078
0079   // Acoplanarity
0080   observables.push_back( "Et_miss_phi" );
0081   observables_name.push_back( "#Delta#phi_{miss - #tau jet}" );
0082   plots_xmin.push_back(0);
0083   plots_xmax.push_back(180);
0084
0085   TString name_CC_base("h_CC_");
0086   TString name_NC_base("h_NC_");
0087   TString name_tau_base("h_tau_");
0088
0089   const int nvar = observables.size();
0090
0091   TH1F* h_CC[nvar];
0092   TH1F* h_NC[nvar];
0093   TH1F* h_tau[nvar];
0094
0095   TString name_CC_i = name_CC_base;
0096   TString name_NC_i = name_NC_base;
0097   TString name_tau_i = name_tau_base;
0098
0099   // create histograms
0100   for ( int l = 0; l < observables.size(); l++ ){
0101     name_CC_i = name_CC_base;
0102     name_NC_i = name_NC_base;
0103     name_tau_i = name_tau_base;
0104
0105     name_CC_i.Append(l);
0106     name_NC_i.Append(l);
0107     name_tau_i.Append(l);
0108
0109     h_CC[l] = new TH1F( name_CC_i, "", 50, plots_xmin.at(l), plots_xmax.at(l));
0110     h_NC[l] = new TH1F( name_NC_i, "", 50, plots_xmin.at(l), plots_xmax.at(l));
0111     h_tau[l] = new TH1F( name_tau_i, "", 50, plots_xmin.at(l), plots_xmax.at(l));
0112   }
0113
0114   // Fill Histograms //
0115   for(int i = 0;i < observables.size(); i++){
0116
0117     ctemp->cd();
0118
0119     name_CC_i = name_CC_base;
0120     name_NC_i = name_NC_base;
0121     name_tau_i = name_tau_base;
0122
0123     name_CC_i.Append(i);
0124     name_NC_i.Append(i);
0125     name_tau_i.Append(i);
0126
0127     // Fill Et_miss straight from tree //
0128     if(i==0){
0129       // Fill CC histogram //
0130       int n_CC = chain_CC.Draw( observables.at(i) + " >> " + name_CC_i, "", "goff");
0131       h_CC[i]->Scale(1./h_CC[i]->Integral());
0132       h_CC[i]->GetXaxis()->SetTitle( observables_name.at(i) );
0133       h_CC[i]->GetYaxis()->SetRangeUser(0, .13 );
0134       h_CC[i]->SetLineColor(col1);
0135       h_CC[i]->SetFillColor(col1);
0136       h_CC[i]->SetFillStyle(0);
0137       h_CC[i]->GetYaxis()->SetTitle("# entries / #Sigma entries");
0138       h_CC[i]->GetYaxis()->SetTitleOffset(1.5);
0139
0140       // Fill NC histogram //
0141       int n_NC = chain_NC.Draw( observables.at(i) + " >> " + name_NC_i, "", "goff");
0142       h_NC[i]->Scale(1./h_NC[i]->Integral());
0143       h_NC[i]->SetLineColor(col3);
0144       h_NC[i]->SetFillColor(col3);
0145       h_NC[i]->SetFillStyle(0);
0146
0147       // Fill LQ histogram //
0148       int n_tau = chain_tau.Draw( observables.at(i) + " >> " + name_tau_i, "", "goff");
0149       h_tau[i]->Scale(1./h_tau[i]->Integral());
0150       h_tau[i]->SetLineColor(col2);
0151       h_tau[i]->SetFillColor(col2);
0152       h_tau[i]->SetFillStyle(0);
0153
0154       TCanvas *c1 = new TCanvas();
0155
0156       h_CC[i]->DrawClone("");
0157       h_NC[i]->DrawClone("sames");
0158       h_tau[i]->DrawClone("sames");
0159       gPad->RedrawAxis();
0160
0161       TLegend *legend = new TLegend(0.4,0.6,0.7,0.89);
0162       legend->AddEntry(h_CC[i],"Neutral Current","l");
0163       legend->AddEntry(h_NC[i],"Charged Current","l");
0164       legend->AddEntry(h_tau[i],"Leptoquark","l");
0165       legend->Draw();
0166     }
0167
0168     // Calculate acoplanarity and then fill histogram //
0169     if(i==1){
0170
0171       float Et_miss_phi_tau, Et_miss_phi_CC, Et_miss_phi_NC;
0172       float tau_phi_tau, tau_phi_CC, tau_phi_NC;
0173
0174       // Read tau phi and missing momentum phi from tree //
0175       chain_CC.SetBranchAddress("Et_miss_phi",&Et_miss_phi_CC);
0176       chain_CC.SetBranchAddress("reco_tau_phi",&tau_phi_CC);
0177
0178       chain_NC.SetBranchAddress("Et_miss_phi",&Et_miss_phi_NC);
0179       chain_NC.SetBranchAddress("reco_tau_phi",&tau_phi_NC);
0180
0181       chain_tau.SetBranchAddress("Et_miss_phi",&Et_miss_phi_tau);
0182       chain_tau.SetBranchAddress("reco_tau_phi",&tau_phi_tau);
0183
0184       // Phi angle between reco tau and miss mom   //
0185       double delta_phi;
0186
0187       // Calculate delta phi for each event type //
0188       for(int j = 0; j<n_CC;j++){
0189     chain_CC.GetEntry(j);
0190     delta_phi = fabs(Et_miss_phi_CC - tau_phi_CC) * 180 / TMath::Pi();
0191     if(delta_phi > 180) delta_phi = 360 - delta_phi;
0192     h_CC[i]->Fill(delta_phi);
0193       }
0194
0195       for(int l = 0; l<n_NC;l++){
0196         chain_NC.GetEntry(l);
0197         delta_phi = fabs(Et_miss_phi_NC - tau_phi_NC) * 180 / TMath::Pi();
0198         if(delta_phi > 180) delta_phi = 360 - delta_phi;
0199         h_NC[i]->Fill(delta_phi);
0200       }
0201
0202       for(int k = 0; k<n_tau;k++){
0203         chain_tau.GetEntry(k);
0204         delta_phi = fabs(Et_miss_phi_tau - tau_phi_tau) * 180 / TMath::Pi();
0205         if(delta_phi > 180) delta_phi = 360 - delta_phi;
0206         h_tau[i]->Fill(delta_phi);
0207       }
0208
0209       //Fill histograms //
0210       h_CC[i]->Scale(1./h_CC[i]->Integral());
0211       h_CC[i]->GetXaxis()->SetTitle( observables_name.at(i) );
0212       h_CC[i]->GetYaxis()->SetRangeUser(0, 0.27 );
0213       h_CC[i]->SetLineColor(col1);
0214       h_CC[i]->SetFillColor(col1);
0215       h_CC[i]->SetFillStyle(0);
0216       h_CC[i]->GetYaxis()->SetTitle("# entries / #Sigma entries");
0217       h_CC[i]->GetYaxis()->SetTitleOffset(1.5);
0218
0219       h_NC[i]->Scale(1./h_NC[i]->Integral());
0220       h_NC[i]->SetLineColor(col3);
0221       h_NC[i]->SetFillColor(col3);
0222       h_NC[i]->SetFillStyle(0);
0223
0224
0225       h_tau[i]->Scale(1./h_tau[i]->Integral());
0226       h_tau[i]->SetLineColor(col2);
0227       h_tau[i]->SetFillColor(col2);
0228       h_tau[i]->SetFillStyle(0);
0229
0230       TCanvas *c1 = new TCanvas();
0231
0232       h_CC[i]->DrawClone("");
0233       h_NC[i]->DrawClone("sames");
0234       h_tau[i]->DrawClone("sames");
0235       gPad->RedrawAxis();
0236
0237       TLegend *legend = new TLegend(0.4,0.6,0.7,0.89);
0238       legend->AddEntry(h_NC[i],"Neutral Current","l");
0239       legend->AddEntry(h_CC[i],"Charged Current","l");
0240       legend->AddEntry(h_tau[i],"Leptoquark","l");
0241       legend->Draw();
0242
0243     }
0244
0245   }
0246
0247   return 0;
0248 }