ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Hammad/macros/EfficiencyPlot.h
Revision: 1.1
Committed: Mon Aug 31 11:53:54 2009 UTC (15 years, 8 months ago) by ghammad
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Log Message:
personal code : first version

File Contents

# User Rev Content
1 ghammad 1.1 #include "/home/gregory/Analysis/TtSemiLeptonic/tdrstyle.C"
2     #include "map.h"
3    
4     void EfficiencyPlot(){
5    
6     //using namespace std;
7    
8     //
9     // load framework libraries
10     //
11     //gSystem->Load( "libFWCoreFWLite" );
12     //AutoLibraryLoader::enable();
13     //
14     // Set P-TDR plots style
15     //
16    
17     setTDRStyle();
18    
19     //TFile *InputFile = TFile::Open("../rootfiles/multijets_EvtSelect_NonIso_4jetsPt30_JetComb_MediumBtagging_obs.root");
20     //TFile *InputFile = TFile::Open("../rootfiles/multijets_EvtSelect_NonIso_VetoIso_4jetsPt30_JetComb_MediumBtagging_obs.root");
21     TFile *InputFile = TFile::Open("../rootfiles/allprocesses_EvtSelect_NonIso_VetoIso_4jetsPt30_JetComb_MediumBtagging_obs.root");
22    
23     //TString ControlHistoFile = "multijets_EvtSelect_NonIso_4jetsPt30_JetComb_MediumBtagging_eff.ps";
24     //TString ControlHistoFile = "multijets_EvtSelect_NonIso_VetoIso_4jetsPt30_JetComb_MediumBtagging_eff.ps";
25     TString ControlHistoFile = "allprocesses_EvtSelect_NonIso_VetoIso_4jetsPt30_JetComb_MediumBtagging_eff.ps";
26    
27     //TFile *OutputFile = new TFile("multijets_EvtSelect_NonIso_VetoIso_4jetsPt30_JetComb_MediumBtagging_eff.root","RECREATE");
28     TFile *OutputFile = new TFile("allprocesses_EvtSelect_NonIso_VetoIso_4jetsPt30_JetComb_MediumBtagging_eff.root","RECREATE");
29    
30     bool Print = false;
31    
32     const int NbOfObs = 19;
33    
34     TString Obs[NbOfObs]={"p_{T}^{#mu}[GeV/c]","H_{T}[GeV]","H_{T}[GeV]","CaloIso[GeV]","TrackIso[GeV/c]","CaloIso/p_{T}","TrackIso/p_{T}","RelIso","|d_{0}|/#sigma","M^{l#nu}_{T}[GeV/c^{2}]","p_{T}^{rel}[GeV/c]","#DeltaR(#mu,jet)","Nb of jets (p_{T}>30GeV/c)","#eta^{#mu}","Nb of b-tagged jets ","m_{qq'} [GeV/c^{2}]","|#eta_{lead. jet}|","#Chi^{2}","M^{bl#nu}_{T}[GeV/c^{2}]"};
35    
36     TString ObsName[NbOfObs]={"MuonPt","HT","HTall","CaloIso","TrackIso","RelCaloIso","RelTrackIso","RelIso","d0Sig","MT","PtRel","DRMuJet","Njets","MuonEta","Nbjets","Mqq","EtaLeadJet","Chi2","LeptMT"};
37    
38     TString EffName[NbOfObs]={"N_{p_{T}^{#mu}<","N_{H_{T}<","N_{H_{T}^{all}<","N_{CaloIso<","N_{TrackIso<","N_{RelCaloIso<","N_{RelTrackIso<","N_{RelIso<","N_{|d0|/#sigma<","N_{M^{l#nu}_{T}<","N_{p_{T}^{rel}<","N_{#DeltaR(#mu,jet)<","N_{Njets<","N_{#eta^{#mu}<","N_{Nbjets<","N_{m_{qq'}<","N_{|#eta_{lead.jet}|<","N_{#Chi^{2}<","N_{M^{bl#nu}_{T}<"};
39    
40     TString EffNameEnd = "}/N";
41    
42     const Double_t Cuts[NbOfObs] = {30,250,300,3,2,0.1,0.1,0.1,3,20,10,0.35,4,2.0,1,150,1.8,5,125}; //Cuts on MuonPt, HT, HTall, CaloIso, TrackIso, RelCaloIso, RelTrackIso, RelIso, MET, d0...
43     TString CutsValue[NbOfObs] = {"30","250","300","3","2","0.1","0.1","0.1","3","20","10","0.35","4","2.0","1","150","1.8","5","125"}; //Cuts on MuonPt, HT, HTall, CaloIso, TrackIso, RelCaloIso, RelTrackIso, RelIso, MET, d0...
44    
45     const Int_t Nbins[NbOfObs] = {13,12,13,12,12,8,8,12,20,9,7,10,10,8,10,10,8,10,10};
46    
47     Double_t MuonPt_bins[Nbins[0]+1] = {0,15,20,22,25,30,40,50,65,80,100,130,165,200};
48     Double_t HT_bins[Nbins[1]+1] = {0,150,200,220,250,270,300,350,400,500,600,700,1000};
49     Double_t HTall_bins[Nbins[2]+1] = {0,150,200,220,250,270,300,350,400,500,600,700,900,1200};
50     Double_t CaloIso_bins[Nbins[3]+1] = {0,1,2,3,4,8,12,16,20,24,28,32,40};
51     Double_t TrackIso_bins[Nbins[4]+1] = {0,1,2,3,4,8,12,16,20,24,28,32,40};
52     Double_t RelCaloIso_bins[Nbins[5]+1] = {0,0.1,0.2,0.3,0.5,0.8,1.3,2,3};
53     Double_t RelTrackIso_bins[Nbins[6]+1] = {0,0.1,0.2,0.3,0.5,0.8,1.3,2,3};
54     Double_t RelIso_bins[Nbins[7]+1] = {0,0.05,0.08,0.12,0.2,0.4,0.7,1.0,1.5,2.0,3.0,4.0,5.0};
55     Double_t d0Sig_bins[Nbins[8]+1] = {0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5,6,7,8,9,10,11,12,13,14,15};
56     Double_t MT_bins[Nbins[9]+1] = {0,10,20,30,40,50,70,100,150,250};
57     Double_t PtRel_bins[Nbins[10]+1] = {0,5,10,15,30,70,200,300};
58     Double_t DRMuJet_bins[Nbins[11]+1] = {0,0.1,0.2,0.3,0.4,0.5,0.8,1.5,3,4,5};
59     Double_t Njets_bins[Nbins[12]+1] = {0,1,2,3,4,5,6,7,8,9,10};
60     Double_t MuonEta_bins[Nbins[13]+1] = {-2.4,-2.1,-1.5,-0.7,0,0.7,1.5,2.1,2.4};
61     Double_t Nbjets_bins[Nbins[14]+1] = {0,1,2,3,4,5,6,7,8,9,10};
62     Double_t Nmqq_bins[Nbins[15]+1] = {0,25,50,100,125,175,200,250,300,350,400};
63     Double_t LeadJetEta_bins[Nbins[16]+1] = {-2.4,-2.1,-1.5,-0.7,0,0.7,1.5,2.1,2.4};
64     Double_t NChi2_bins[Nbins[17]+1] = {0,2.5,5,10,12,17,20,25,30,35,40};
65     Double_t Nmblv_bins[Nbins[18]+1] = {0,25,50,100,125,175,200,250,300,350,400};
66    
67     Double_t **bins = new Double_t*[NbOfObs];
68     for(Int_t i=0;i<NbOfObs;i++)
69     {
70     bins[i] = new Double_t[Nbins[i]+1];
71     }
72    
73     bins[0] = MuonPt_bins;
74     bins[1] = HT_bins;
75     bins[2] = HTall_bins;
76     bins[3] = CaloIso_bins;
77     bins[4] = TrackIso_bins;
78     bins[5] = RelCaloIso_bins;
79     bins[6] = RelTrackIso_bins;
80     bins[7] = RelIso_bins;
81     bins[8] = d0Sig_bins;
82     bins[9] = MT_bins;
83     bins[10] = PtRel_bins;
84     bins[11]= DRMuJet_bins;
85     bins[12]= Njets_bins;
86     bins[13]= MuonEta_bins;
87     bins[14]= Nbjets_bins;
88     bins[15]= Nmqq_bins;
89     bins[16]= LeadJetEta_bins;
90     bins[17]= NChi2_bins;
91     bins[18]= Nmblv_bins;
92     //bins[12]= Pt4thJet_bins;
93    
94     std::cout<<"Initialization done"<<std::endl;
95    
96     // --------------------------------------------------------------
97     // Pointers to histograms
98     // --------------------------------------------------------------
99    
100     TDirectory *dir = InputFile;
101     //(TDirectory*) InputFile->Get("TSAnalyser");
102    
103     TString ObsToGet = "";
104     TString ObsToGet_EffX = "";
105     TString ObsToGet_EffY = "";
106     Int_t Index = 0;
107    
108     TH2F *hHisto2D[NbOfObs*(NbOfObs-1)/2];
109     TH1F *hHistoEff[NbOfObs*(NbOfObs-1)/2][2];
110    
111     Double_t Eff_X = 0;
112     Double_t Eff_Y = 0;
113     Double_t S = 0;
114     Double_t Ntot = 0;
115    
116     Double_t Xmin = 0;
117     Double_t Ymin = 0;
118     Double_t Xmax = 0;
119     Double_t Ymax = 0;
120    
121     Int_t NbinsX = 0;
122     Int_t NbinsY = 0;
123    
124     Int_t X0 = 0;
125     Int_t Y0 = 0;
126    
127     Double_t error = 0;
128    
129     TString YaxisName = "";
130    
131     OutputFile->cd();
132    
133     TCanvas c("dummy","",1);
134     c.Print(ControlHistoFile + "[");
135    
136     for(Int_t i = 0; i<NbOfObs ; i++)
137     {
138     for(Int_t j = i+1; j<NbOfObs ; j++)
139     {
140     ObsToGet = "Obs_"; ObsToGet += ObsName[i]; ObsToGet +="_vs_Obs_"; ObsToGet += ObsName[j];
141    
142     ObsToGet_EffX = ObsToGet; ObsToGet_EffX += "_EffX";
143     ObsToGet_EffY = ObsToGet; ObsToGet_EffY += "_EffY";
144     cout<<ObsToGet<<endl;
145    
146     hHisto2D[Index] = (TH2F *) dir->Get(ObsToGet);
147    
148     ////////////////// Efficiency along the X axis /////////////////////////////
149     hHistoEff[Index][0] = new TH1F(ObsToGet_EffX,ObsToGet_EffX,Nbins[i],bins[i]);
150     hHistoEff[Index][0]->SetTitle("");
151     hHistoEff[Index][0]->GetXaxis()->SetTitle(Obs[i]);
152     YaxisName = EffName[j]; YaxisName += CutsValue[j]; YaxisName += EffNameEnd;
153     hHistoEff[Index][0]->GetYaxis()->SetTitle(YaxisName);
154     hHistoEff[Index][0]->GetYaxis()->SetRangeUser(0,1.1);
155     //hHistoEff[Index][0]->GetYaxis()->SetRangeUser(0.5*hHistoEff[Index][0]->GetMinimum(),1.5*hHistoEff[Index][0]->GetMaximum());
156    
157     ////////////////// Efficiency along the Y axis /////////////////////////////
158     hHistoEff[Index][1] = new TH1F(ObsToGet_EffY,ObsToGet_EffY,Nbins[j],bins[j]);
159     hHistoEff[Index][1]->SetTitle("");
160     hHistoEff[Index][1]->GetXaxis()->SetTitle(Obs[j]);
161     YaxisName = EffName[i]; YaxisName += CutsValue[i]; YaxisName += EffNameEnd;
162     hHistoEff[Index][1]->GetYaxis()->SetTitle(YaxisName);
163     hHistoEff[Index][1]->GetYaxis()->SetRangeUser(0,1.1);
164     //hHistoEff[Index][1]->GetYaxis()->SetRangeUser(0.5*hHistoEff[Index][1]->GetMinimum(),1.5*hHistoEff[Index][1]->GetMaximum());
165    
166     NbinsX = hHisto2D[Index]->GetNbinsX();
167     NbinsY = hHisto2D[Index]->GetNbinsY();
168     Xmin = hHisto2D[Index]->GetXaxis()->GetXmin();Xmax = hHisto2D[Index]->GetXaxis()->GetXmax();
169     Ymin = hHisto2D[Index]->GetYaxis()->GetXmin();Ymax = hHisto2D[Index]->GetYaxis()->GetXmax();
170     X0 = NbinsX*(Cuts[i]-Xmin)/(Xmax-Xmin);
171     Y0 = NbinsY*(Cuts[j]-Ymin)/(Ymax-Ymin);
172    
173     if(Print){
174     std::cout<<"----------------------------------------------------------------------------------"<<std::endl;
175     std::cout<<"------------------------------"<<ObsName[i]<<" vs "<<ObsName[j]<<"------------------------------"<<std::endl;
176     std::cout<<"----------------------------------------------------------------------------------"<<std::endl;
177     std::cout<<"--------"<<ObsName[i]<<"range = "<<Xmin<<"-"<<Xmax<<" / Cut = "<<Cuts[i]<<" / Nb of bins = "<<NbinsX<<" / X0 = "<<X0<<" --------"<<std::endl;
178     std::cout<<"--------"<<ObsName[j]<<"range = "<<Ymin<<"-"<<Ymax<<" / Cut = "<<Cuts[j]<<" / Nb of bins = "<<NbinsY<<" / Y0 = "<<Y0<<" --------"<<std::endl;
179     }
180    
181     if(Print) cout<<"X Binning = ";
182     for(Int_t k=1; k <= Nbins[i]; k++)
183     {
184     if(Print) cout<<bins[i][k]<<"/";
185     if(hHisto2D[Index]->Integral((bins[i][k-1]*NbinsX/Xmax)+1,bins[i][k]*NbinsX/Xmax,1,NbinsY) > 0 )
186     {
187     S = hHisto2D[Index]->Integral((bins[i][k-1]*NbinsX/Xmax)+1,bins[i][k]*NbinsX/Xmax,1,Y0);
188     Ntot = hHisto2D[Index]->Integral((bins[i][k-1]*NbinsX/Xmax)+1,bins[i][k]*NbinsX/Xmax,1,NbinsY);
189    
190     hHistoEff[Index][0]->SetBinContent(k,S/Ntot);
191     hHistoEff[Index][0]->SetBinError(k,EffError(S,Ntot));
192     }
193     else hHistoEff[Index][0]->SetBinContent(k,0);
194     }
195     if(Print) cout<<endl;
196     if(Print) cout<<"Y Binning = ";
197     for(Int_t k=1; k <= Nbins[j]; k++)
198     {
199     if(Print) cout<<bins[j][k]<<"/";
200     if( hHisto2D[Index]->Integral(1 ,NbinsX,(bins[j][k-1]*NbinsY/Ymax)+1,bins[j][k]*NbinsY/Ymax) > 0 )
201     {
202     //S = hHisto2D[Index]->Integral(X0,NbinsX,(bins[j][k-1]*NbinsY/Ymax)+1,bins[j][k]*NbinsY/Ymax);
203     S = hHisto2D[Index]->Integral(1, X0,(bins[j][k-1]*NbinsY/Ymax)+1,bins[j][k]*NbinsY/Ymax);
204     Ntot = hHisto2D[Index]->Integral(1 ,NbinsX,(bins[j][k-1]*NbinsY/Ymax)+1,bins[j][k]*NbinsY/Ymax);
205    
206     hHistoEff[Index][1]->SetBinContent(k,S/Ntot);
207     hHistoEff[Index][1]->SetBinError(k,EffError(S,Ntot));
208     }
209     else hHistoEff[Index][1]->SetBinContent(k,0);
210     }
211     hHistoEff[Index][0]->Draw();
212     hHistoEff[Index][0]->Write(ObsToGet_EffX);
213     c.Print(ControlHistoFile);
214     hHistoEff[Index][1]->Draw();
215     hHistoEff[Index][1]->Write(ObsToGet_EffY);
216     c.Print(ControlHistoFile);
217     Index++;
218     }
219     }
220    
221    
222     //--------------------------------------------------------------------------------------------------------------------------------------------
223     //--------------------------------------------------------------------------------------------------------------------------------------------
224     //--------------------------------------------------------------------------------------------------------------------------------------------
225    
226     OutputFile->Close();
227     c.Print(ControlHistoFile + "]");
228     gROOT->ProcessLine(".q");
229    
230     }
231    
232     Double_t EffError(Double_t S, Double_t Total)
233     {
234     Double_t Eff = S/Total;
235     return (sqrt(S*(1-Eff))/Total);
236     }
237    
238