| 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 |
}
|
