ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/miscFR/test.C
Revision: 1.1
Committed: Fri Jun 24 07:44:49 2011 UTC (13 years, 10 months ago) by fronga
Content type: text/plain
Branch: MAIN
CVS Tags: cbaf_4_98ifb_paper, cbaf_4p7ifb, Honeypot, cbaf_2p1ifb, HEAD
Log Message:
First commit of various scripts used in the first 2011 analyses.

File Contents

# User Rev Content
1 fronga 1.1 //________________________________________________________________________________________
2     // Compare mumu and ee JZB distributions in data
3     int jzbMuEcomp(void) {
4    
5     TString jzbType("jzb[1]");
6     TString jzbName("Run2011AB");
7     // TString jzbType("jzb[4]");
8     // TString jzbName("TC");
9    
10     TFile *_file0 = TFile::Open("/scratch/buchmann/AllDataCertified.root");
11     TTree* events = (TTree*)_file0->Get("events");
12     TCut kbase("abs(mll-91.2)<20&&pfJetGoodNum>1"); //&&passed_triggers");
13     TCut kSF("id1==id2");
14     TCut kOF("id1!=id2");
15    
16     gStyle->SetOptStat(111111110);
17    
18     // Do the fit of the peak on the full sample
19     TH1F* hJZBll = new TH1F("hJZBll","JZB ; ee+#mu#mu JZB [GeV]",40,-40,40);
20     TH1F* hJZBem = new TH1F("hJZBem","JZB ; e#mu JZB [GeV]",40,-40,40);
21    
22     hJZBll->Sumw2();
23     hJZBem->Sumw2();
24     hJZBll->SetMarkerColor(kBlue);
25    
26     events->Project("hJZBll",jzbType,kbase&&kSF);
27     events->Project("hJZBem",jzbType,kbase&&kOF);
28    
29     // Get the mean
30     TCanvas* c1 = new TCanvas(); //c1.Divide(1,2);
31     TF1* gauss = new TF1("gauss","gaus");
32     Float_t range = 20.;
33     Float_t mean = 0.;
34     int niter = 10;
35     for ( int i=0; i<niter; ++i ) {
36     gauss->SetRange(mean-range,mean+range);
37     hJZBll->Fit(gauss,"QR","");
38     mean = gauss->GetParameter(1);
39     range = 2.0*gauss->GetParameter(2);
40     std::cout << "Mean: " << mean << "+-" << range << std::endl;
41     }
42    
43     char var[30]; sprintf(var,(jzbType+"%+5.3f").Data(),-mean);
44     std::cout << var << std::endl;
45    
46    
47     hJZBll->Draw();
48     hJZBem->Draw("same");
49     c1->SaveAs("plots/jzbFitEEmm"+jzbName+".eps");
50    
51     // Project ee and mumu separately
52     Float_t xbins[] = { -100., -80., -60., -50., -40., -30., -20., -10., 0., 10, 20., 30., 40., 50., 60., 80., 100. };
53     Int_t nbins = sizeof(xbins)/4-1;
54     TH1F* hJZBee2 = new TH1F("hJZBee2","JZB ; ee JZB [GeV]",nbins,xbins);
55     TH1F* hJZBmm2 = new TH1F("hJZBmm2","JZB ; #mu#mu JZB [GeV]",nbins,xbins);
56     hJZBee2->Sumw2();
57     hJZBmm2->Sumw2();
58     TH1F* hJZBemee = new TH1F("hJZBemee","JZB ; e#mu (ee weighted) JZB [GeV]",nbins,xbins);
59     TH1F* hJZBemmm = new TH1F("hJZBemmm","JZB ; e#mu (#mu#mu weighted) JZB [GeV]",nbins,xbins);
60     hJZBemee->Sumw2();
61     hJZBemmm->Sumw2();
62    
63     TCut of(kbase&&kOF);
64     events->Project("hJZBee2",var,kbase&&kSF&&"id1==0");
65     events->Project("hJZBmm2",var,kbase&&kSF&&"id1==1");
66     events->Project("hJZBemee",var,kbase&&kOF);
67     events->Project("hJZBemmm",var,kbase&&kOF);
68    
69     // Reweight by the efficiency
70     Float_t eEff = 0.82, muEff = 0.94;
71     Float_t ratio = eEff/muEff;
72     hJZBemee->Scale(0.5*ratio);
73     hJZBemmm->Scale(0.5/ratio);
74    
75     // Subtract em from ll (but make sure there is no negative bin
76     std::cout << "hJZBemee " << hJZBemee->Integral() << std::endl;
77     std::cout << "hJZBemmm " << hJZBemmm->Integral() << std::endl;
78     std::cout << "hJZBee2 " << hJZBee2->Integral() << std::endl;
79     std::cout << "hJZBmm2 " << hJZBmm2->Integral() << std::endl;
80    
81     for ( int ibin=0; ibin<hJZBee2->GetNbinsX(); ++ibin ) {
82     Float_t diff = hJZBee2->GetBinContent(ibin)-hJZBemee->GetBinContent(ibin);
83     if (diff>0.) hJZBee2->SetBinContent(ibin,diff);
84     diff = hJZBmm2->GetBinContent(ibin)-hJZBemmm->GetBinContent(ibin);
85     if (diff>0.) hJZBmm2->SetBinContent(ibin,diff);
86     }
87     std::cout << "hJZBee2- " << hJZBee2->Integral() << std::endl;
88     std::cout << "hJZBmm2- " << hJZBmm2->Integral() << std::endl;
89    
90     // Get the ratio
91     TH1F* hRatio = (TH1F*)hJZBee2->Clone();
92     hRatio->SetTitle("JZB ratio ; JZB ee/#mu#mu");
93     hRatio->Sumw2();
94     hRatio->Divide(hJZBmm2);
95    
96     // Plot JZBs
97     TCanvas* c2 = new TCanvas();
98     c2->Divide(1,2);
99     c2->cd(1);
100     gPad->SetLogy(1);
101     hJZBee2->SetMarkerColor(kRed);
102     hJZBmm2->SetMarkerColor(kBlue);
103     hJZBee2->Draw("E1");
104     //hJZBemee->Draw("same");
105    
106     // Draw ee and mm superimposed
107     c2->cd(2);
108     gPad->SetLogy(1);
109     TH1F* hJZBee3 = (TH1F*)hJZBee2->Clone();
110     hJZBee3->Scale(hJZBmm2->Integral()/hJZBee2->Integral());
111     hJZBee3->SetFillColor(kRed);
112     hJZBmm2->Draw("E1");
113     hJZBee3->Draw("E3same");
114     //hJZBemmm->Draw("same");
115     hJZBmm2->Draw("E1same");
116     c2->SaveAs("plots/jzbCompEEmm"+jzbName+".eps");
117    
118     // Plot ratio
119     TCanvas* c3 = new TCanvas();
120     TF1* pol = new TF1("pol","pol1");
121     hRatio->SetMaximum(1.5);
122     hRatio->SetMinimum(0.);
123     hRatio->Fit("pol0","","E1");
124     c3->SaveAs("plots/jzbRatioEEmm"+jzbName+".eps");
125    
126     return 0;
127    
128     }
129    
130    
131     int compHist( TTree* t1, TTree* t2, TCut kbase, TString var, TString title, int nbins, float xmin, float xmax ) {
132    
133     TH1F* h1 = new TH1F("h1","original ; "+title,nbins,xmin,xmax);
134     TH1F* h2 = new TH1F("h2","boosted ; "+title,nbins,xmin,xmax);
135     h1->SetLineColor(kBlue);
136     h1->SetLineWidth(2);
137     h2->SetLineColor(kRed);
138     h2->SetLineWidth(2);
139    
140     // Project on histograms
141     t1->Project(h1->GetName(),var,kbase);
142     t2->Project(h2->GetName(),var,kbase);
143     // Scale largest to smallest
144     if ( h1->GetEntries()>h2->GetEntries() )
145     h1->Scale(h2->Integral()/h1->Integral());
146     else
147     h2->Scale(h1->Integral()/h2->Integral());
148     // Make sure maximum is good for both
149     if ( h1->GetMaximum()<h2->GetMaximum() )
150     h1->SetMaximum( h2->GetMaximum()*1.05 );
151    
152     h1->Draw();
153     h2->Draw("sames");
154     TVirtualPad::Pad()->Update();
155    
156     TPaveStats *s1 = (TPaveStats*)h1->GetListOfFunctions()->FindObject("stats");
157     TPaveStats *s2 = (TPaveStats*)h2->GetListOfFunctions()->FindObject("stats");
158     s1->SetTextColor(h1->GetLineColor());
159     s2->SetTextColor(h2->GetLineColor());
160     s2->SetY1NDC(s2->GetY1NDC()-0.2);
161     s2->SetY2NDC(s2->GetY2NDC()-0.2);
162     TVirtualPad::Pad()->Update();
163    
164    
165     // TLegend* legend = new TLegend();
166    
167     }
168    
169     //________________________________________________________________________________________
170     int getAsym(TString file,TCut kmll = "abs(mll-91)<20") {
171    
172     TFile *_file2 = TFile::Open("dcap://t3se01.psi.ch:22125/pnHifs/psi.ch/cms/trivcat/store/user/fronga/DYJetsToLL"+file+".root");
173     TTree* tree = (TTree*)_file2->Get("events");
174     TCut ksel("pfJetGoodNum>1&&id1==id2&&id1==1&&pt1>20&&pt2>20");
175    
176     // Get numbers, left and right
177     Int_t n1 = tree->Draw("abs(jzb[1]+1.3)",ksel&&kmll&&"jzb[1]+1.3>50","e1");
178     Int_t n2 = tree->Draw("abs(jzb[1]+1.3)",ksel&&kmll&&"jzb[1]+1.3<-50","same,hist");
179    
180     std::cout << "right: " << n1 << ", left " << n2 << std::endl;
181    
182     return 0;
183     }
184    
185    
186     //________________________________________________________________________________________
187     // Bifurcated gaussian
188     Double_t bigauss(Double_t *v, Double_t *par)
189     {
190     // par[3] = (const,mean,sigmaL,sigmaR)
191    
192     Double_t arg = v[0]-par[1];
193     Double_t coef = 0.0;
194    
195     if ( v[0]<par[1] ) {
196     if ( fabs(par[2])>0. ) coef = -0.5/(par[2]*par[2]);
197     } else {
198     if ( fabs(par[3])>0. ) coef = -0.5/(par[3]*par[3]);
199     }
200    
201     return par[0]*TMath::Exp(coef*arg*arg);
202    
203     }
204    
205    
206     //________________________________________________________________________________________
207     // Compare different types of MET measurements
208     int compMET(void) {
209    
210     TFile *_file0 = TFile::Open("/shome/theofil/JZBAnalysis/JZBsamples/PabloV13/All-Data.root");
211     TTree* events = (TTree*)_file0->Get("events");
212     TCut kbase("abs(mll-91.2)<20&&pfJetGoodNum>1&&id1==id2&&abs(jzb[1])>45");
213    
214     TCanvas* c1 = new TCanvas(); //c1.Divide(1,2);
215     TH1F* hcomp = new TH1F("hcomp","JZB comparison; (pfJZB-tcJZB)/tcJZB;",15,-0.5,1);
216    
217     gStyle->SetOptStat(111110);
218     events->Project("hcomp","(jzb[1]-jzb[4])/jzb[4]",kbase);
219     hcomp->Draw();
220    
221     c1->SaveAs("plots/metComp.eps");
222    
223    
224    
225     }
226    
227    
228     //________________________________________________________________________________________
229     // Compare mumu and ee JZB distributions in DY MC
230     int jzbMuEcompMC(void) {
231    
232     TString jzbType("jzb[1]");
233     TString jzbName("PF");
234     // TString jzbType("jzb[4]");
235     // TString jzbName("TC");
236    
237     TFile *_file0 = TFile::Open("dcap://t3se01.psi.ch:22125/pnfs/psi.ch/cms/trivcat/store/user/buchmann/MCSpring2011PU/DYJetsToLL_TuneZ2_M-50_7TeV-madgraph-tauola.root");
238     TTree* events = (TTree*)_file0->Get("events");
239     TCut kbase("abs(mll-91.2)<20&&pfJetGoodNum>1");
240     TCut kSF("id1==id2");
241     TCut kOF("id1!=id2");
242    
243     gStyle->SetOptStat(111111110);
244    
245     // Do the fit of the peak on the full sample
246     TH1F* hJZBee = new TH1F("hJZBee","JZB ; ee JZB [GeV]",40,-40,40);
247     TH1F* hJZBmm = new TH1F("hJZBmm","JZB ; #mu#mu JZB [GeV]",40,-40,40);
248    
249     hJZBee->Sumw2();
250     hJZBmm->Sumw2();
251    
252     events->Project("hJZBee",jzbType,kbase&&kSF&&"id1==0");
253     events->Project("hJZBmm",jzbType,kbase&&kSF&&"id1==1");
254    
255     // Get the mean
256     TCanvas* c1 = new TCanvas(); c1->Divide(1,2);
257     TF1* gauss = new TF1("gauss","gaus");
258     Float_t range = 10.;
259     Float_t mean = 0.;
260     int niter = 5;
261     c1->cd(1);
262     for ( int i=0; i<niter; ++i ) {
263     gauss->SetRange(mean-range,mean+range);
264     hJZBee->Fit(gauss,"QR","");
265     mean = gauss->GetParameter(1);
266     range = 2.0*gauss->GetParameter(2);
267     std::cout << "ee mean: " << mean << "+-" << range << std::endl;
268     }
269     char var1[30]; sprintf(var1,(jzbType+"%+5.3f").Data(),-mean);
270     std::cout << var1 << std::endl;
271    
272     range = 10.;
273     mean = 0.;
274     for ( int i=0; i<niter; ++i ) {
275     gauss->SetRange(mean-range,mean+range);
276     hJZBmm->Fit(gauss,"QR","");
277     mean = gauss->GetParameter(1);
278     range = 2.0*gauss->GetParameter(2);
279     std::cout << "mumu mean: " << mean << "+-" << range << std::endl;
280     }
281     char var2[30]; sprintf(var2,(jzbType+"%+5.3f").Data(),-mean);
282     std::cout << var2 << std::endl;
283    
284     c1->cd(1);
285     hJZBee->Draw();
286     c1->cd(2);
287     hJZBmm->Draw("");
288     //c1->SaveAs("plots/jzbFitEEmm"+jzbName+".eps");
289    
290     // Project ee and mumu separately
291     //Float_t xbins[] = { -100., -80., -60., -50., -40., -30., -20., -10., 0., 10, 20., 30., 40., 50., 60., 80., 100. };
292     //Int_t nbins = sizeof(xbins)/4-1;
293     Int_t nbins = 50;
294     Float_t min = -100, max = 100;
295     TH1F* hJZBee2 = new TH1F("hJZBee2","JZB ; ee JZB [GeV]",nbins,min,max);
296     TH1F* hJZBmm2 = new TH1F("hJZBmm2","JZB ; #mu#mu JZB [GeV]",nbins,min,max);
297     hJZBee2->Sumw2();
298     hJZBmm2->Sumw2();
299    
300     TCut of(kbase&&kOF);
301     events->Project("hJZBee2",var1,kbase&&kSF&&"id1==0");
302     events->Project("hJZBmm2",var2,kbase&&kSF&&"id1==1");
303    
304     // Reweight by the efficiency
305     Float_t eEff = 0.82, muEff = 0.94;
306     Float_t ratio = eEff/muEff;
307    
308     // Subtract em from ll (but make sure there is no negative bin
309     std::cout << "hJZBee2 " << hJZBee2->Integral() << std::endl;
310     std::cout << "hJZBmm2 " << hJZBmm2->Integral() << std::endl;
311    
312     // Get the ratio
313     TH1F* hRatio = (TH1F*)hJZBee2->Clone();
314     hRatio->SetTitle("JZB ratio ; JZB ee/#mu#mu");
315     hRatio->Sumw2();
316     hRatio->Divide(hJZBmm2);
317    
318     // Plot JZBs
319     TCanvas* c2 = new TCanvas();
320     c2->Divide(1,2);
321     c2->cd(1);
322     gPad->SetLogy(1);
323     hJZBee2->SetMarkerColor(kRed);
324     hJZBmm2->SetMarkerColor(kBlue);
325     hJZBee2->Draw("E1");
326     //hJZBemee->Draw("same");
327    
328     // Draw ee and mm superimposed
329     c2->cd(2);
330     gPad->SetLogy(1);
331     TH1F* hJZBee3 = (TH1F*)hJZBee2->Clone();
332     hJZBee3->Scale(hJZBmm2->Integral()/hJZBee2->Integral());
333     hJZBee3->SetFillColor(kRed);
334     hJZBmm2->Draw("E1");
335     hJZBee3->Draw("E3same");
336     hJZBmm2->Draw("E1same");
337     //c2->SaveAs("plots/jzbCompEEmm"+jzbName+".eps");
338    
339     // Plot ratio
340     TCanvas* c3 = new TCanvas();
341     TF1* pol = new TF1("pol","pol1");
342     hRatio->SetMaximum(1.5);
343     hRatio->SetMinimum(0.);
344     hRatio->Fit("pol0","","E1");
345     //c3->SaveAs("plots/jzbRatioEEmm"+jzbName+".eps");
346    
347     return 0;
348    
349     }