File indexing completed on 2025-08-05 08:13:20
0001
0002 #include <iostream>
0003 #include <iomanip>
0004 #include <fstream>
0005 #include <utility>
0006 #include <TRandom.h>
0007 #include <TH1D.h>
0008 #include <TH1F.h>
0009 #include <TH2F.h>
0010 #include <TROOT.h>
0011 #include <TStyle.h>
0012 #include "TCanvas.h"
0013 #include <TFile.h>
0014 #include <TTree.h>
0015 #include <TCanvas.h>
0016 #include <TBox.h>
0017 #include "TF1.h"
0018 #include "TDirectory.h"
0019 #include "TDirectoryFile.h"
0020 #include "TGraph.h"
0021 #include "TMath.h"
0022 #include "TLegend.h"
0023 #include "TLine.h"
0024 #include "TLatex.h"
0025
0026 using namespace std;
0027
0028
0029
0030 static const int nbins_truth = 50;
0031 static const double boundaries_truth[nbins_truth+1] = {
0032 0.,10.,20.,30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170,
0033 180, 190, 200, 210, 220, 230, 240, 250, 260, 270, 280, 290, 300, 310,
0034 320, 330, 340, 350, 360, 370, 380, 390, 400, 410, 420, 430, 440, 450,
0035 460, 470, 480, 490, 500
0036 };
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049 static const int nbins_rec = 50;
0050 static const double boundaries_rec[nbins_rec+1] = {
0051 0,10,20,30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170,
0052 180, 190, 200, 210, 220, 230, 240, 250, 260, 270, 280, 290, 300, 310,
0053 320, 330, 340, 350, 360, 370, 380, 390, 400, 410, 420, 430, 440, 450,
0054 460, 470, 480, 490, 500
0055 };
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085 static const int nbins_recrebin = 25;
0086 static const double boundaries_recrebin[nbins_recrebin+1] = {
0087 30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,200,220,240,260,280,300,340,380,420,460
0088 };
0089
0090
0091 static const int nbins_recrebinM = 19;
0092 static const double boundaries_recrebinM[nbins_recrebinM+1] = {
0093 0,5,10,15,20,25,30,35,40,45,50,55,60,70,80,90,110,130,170,250
0094 };
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108 static const int nbins_recrebin_Npart = 1;
0109 static const double boundaries_recrebin_Npart[nbins_recrebin_Npart+1] = {
0110 100, 300
0111 };
0112
0113 static const int nColor = 5;
0114 static const int colorCode[nColor] = {
0115 1, 2, kGreen+1, 4, 6
0116 };
0117
0118
0119 static const int nAlgos = 8;
0120 static const int BinLabelN = 11;
0121 static const char *algoName[nAlgos] = { "", "icPu5", "akPu2PF", "akPu3PF", "akPu4PF", "akPu2Calo", "akPu3Calo", "akPu4Calo" };
0122 static const char *algoNamePP[nAlgos] = { "", "icPu5", "ak2PF", "ak3PF", "ak4PF", "ak2Calo", "ak3Calo", "ak4Calo" };
0123 static const char *algoNameGen[nAlgos] = { "", "icPu5", "akPu2PF", "akPu3PF", "akPu4PF", "akPu2PF", "akPu3PF", "akPu4PF" };
0124 static const char *BinLabel[BinLabelN] = {"100-110", "110-120", "120-130", "130-140", "140-150", "150-160", "160-170", "170-180", "180-200", "200-240","240-300" };
0125
0126
0127 const int nbins_cent= 6;
0128 Double_t boundaries_cent[nbins_cent+1] = { 0,2,4,12,20,28,36 };
0129 Double_t ncoll[nbins_cent] = { 1660, 1310, 745, 251, 62.8, 10.8 };
0130
0131
0132
0133
0134
0135
0136
0137
0138 TGraphErrors *tTAAerr[6]={0};
0139 TGraphErrors *tTAAerrNpart=0;
0140
0141
0142
0143
0144
0145
0146
0147
0148 TH1F *functionHist(TF1 *f, TH1F* h,char *fHistname)
0149 {
0150 TH1F *hF = (TH1F*)h->Clone(fHistname);
0151 for (int i=1;i<=h->GetNbinsX();i++)
0152 {
0153 double var = f->Integral(h->GetBinLowEdge(i),h->GetBinLowEdge(i+1))/h->GetBinWidth(i);
0154 hF->SetBinContent(i,var);
0155 hF->SetBinError(i,0);
0156 }
0157 return hF;
0158 }
0159
0160 TLegend *myLegend(double x1,double y1,double x2, double y2)
0161 {
0162 TLegend *leg = new TLegend(x1,y1,x2,y2);
0163 leg->SetBorderSize(0);
0164 leg->SetFillStyle(0);
0165 return leg;
0166
0167 }
0168
0169
0170
0171 TH1F* drawEnvelope(TH1F *h,char *opt,int color = kGray,int fillStyle = 0, int fillColor = 0,double shift = 0)
0172 {
0173 TH1F *hClone = (TH1F*) h->Clone(Form("%s_mirror",h->GetTitle()));
0174 TH1F *hMirror = (TH1F*) h->Clone(Form("%s_mirror2",h->GetTitle()));
0175 for (int i=1;i<=h->GetNbinsX();i++)
0176 {
0177 double val = h->GetBinContent(i);
0178 hMirror->SetBinContent(i,1-fabs(val-1)+shift);
0179 hClone->SetBinContent(i,val+shift);
0180 }
0181
0182
0183
0184 hMirror->SetLineColor(color);
0185 hMirror->SetFillColor(fillColor);
0186 hMirror->SetFillStyle(fillStyle);
0187 hClone->SetLineColor(color);
0188 hClone->SetFillColor(fillColor);
0189 hClone->SetFillStyle(fillStyle);
0190 hClone->Draw(opt);
0191 hMirror->Draw(opt);
0192 return hMirror;
0193 }
0194
0195 void makeHistTitle(TH1 *h,char *title, char *xTitle, char *yTitle, int color = -1, bool centerTitle = 1)
0196 {
0197 h->SetTitle(title);
0198 h->SetXTitle(xTitle);
0199 h->SetYTitle(yTitle);
0200
0201 if (centerTitle) {
0202 h->GetXaxis()->CenterTitle();
0203 h->GetYaxis()->CenterTitle();
0204
0205 }
0206
0207 if (color!=-1) {
0208 h->SetLineColor(color);
0209 h->SetMarkerColor(color);
0210 }
0211
0212 h->GetYaxis()->SetNdivisions(610);
0213 h->GetXaxis()->SetNdivisions(505);
0214
0215
0216 h->GetYaxis()->SetLabelFont(43);
0217 h->GetYaxis()->SetTitleFont(43);
0218 h->GetYaxis()->SetLabelSize(20);
0219 h->GetYaxis()->SetTitleSize(22);
0220 h->GetYaxis()->SetTitleOffset(2.6);
0221
0222 h->GetXaxis()->SetLabelFont(43);
0223 h->GetXaxis()->SetTitleFont(43);
0224 h->GetXaxis()->SetLabelSize(20);
0225 h->GetXaxis()->SetTitleSize(22);
0226 h->GetXaxis()->SetTitleOffset(3.1);
0227
0228 h->GetXaxis()->SetNoExponent();
0229 h->GetXaxis()->SetMoreLogLabels();
0230
0231 h->GetXaxis()->SetTitleOffset(2.4);
0232
0233 }
0234
0235
0236
0237 class SysData
0238 {
0239 public:
0240 SysData() {
0241 for (int i=0;i<=nbins_cent;i++) {
0242 hSys[i] = new TH1F(Form("hSys_cent%d",i), Form("Totalsys_cent%d",i),nbins_recrebin, boundaries_recrebin);
0243 hSysGeneral[i]= new TH1F(Form("hSysGeneral_cent%d",i), Form("TotalsysGeneral_cent%d",i),nbins_recrebin, boundaries_recrebin);
0244 hSysJEC[i] = new TH1F(Form("hSysJEC_cent%d",i), Form("JECsys_cent%d",i),nbins_recrebin, boundaries_recrebin);
0245 hSysEff[i] = new TH1F(Form("hSysEff_cent%d",i), Form("Effsys_cent%d",i),nbins_recrebin, boundaries_recrebin);
0246 hSysSmear[i] = new TH1F(Form("hSysSmear_cent%d",i), Form("Smearsys_cent%d",i),nbins_recrebin, boundaries_recrebin);
0247 hSysIter[i] = new TH1F(Form("hSysIter_cent%d",i), Form("Itersys_cent%d",i),nbins_recrebin, boundaries_recrebin);
0248 hSys[i]->SetLineColor(kGray);
0249 hSysJEC[i]->SetLineColor(4);
0250 hSysSmear[i]->SetLineColor(kGreen+1);
0251 hSysIter[i]->SetLineColor(2);
0252 }
0253 }
0254 TH1F *hSys[nbins_cent+1];
0255 TH1F *hSysGeneral[nbins_cent+1];
0256 TH1F *hSysJEC[nbins_cent+1];
0257 TH1F *hSysEff[nbins_cent+1];
0258 TH1F *hSysSmear[nbins_cent+1];
0259 TH1F *hSysIter[nbins_cent+1];
0260 TH1F *hSysNoise[nbins_cent+1];
0261
0262 void calcTotalSys(int i) {
0263 TF1 *fNoise = new TF1("f","1+0.3*0.16*abs(1-([0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x))");
0264 fNoise->SetParameters(0.9521,0.001105,-9.397e-6,3.32e-8,-5.618e-11);
0265 hSysNoise[i] = functionHist(fNoise,hSys[i],Form("hSysNoise_cent%d",i));
0266 hSysNoise[i]->SetName(Form("hSysNoise_cent%d",i));
0267 hSysNoise[i]->SetLineColor(6);
0268 for (int j=1;j<=hSys[i]->GetNbinsX();j++) {
0269 double effSys = 0.01;
0270 hSysEff[i]->SetBinContent(j,1+effSys);
0271 hSysSmear[i]->SetBinContent(j,1.02);
0272 double JECSys = hSysJEC[i]->GetBinContent(j)-1;
0273 double SmearSys = hSysSmear[i]->GetBinContent(j)-1;
0274 double IterSys = hSysIter[i]->GetBinContent(j)-1;
0275 double NoiseSys = hSysNoise[i]->GetBinContent(j)-1;
0276 cout <<effSys<<" "<<JECSys<<" "<<IterSys<<endl;
0277 double totalSys = sqrt( effSys * effSys +
0278 JECSys * JECSys +
0279 SmearSys * SmearSys +
0280 NoiseSys * NoiseSys +
0281 IterSys* IterSys
0282 );
0283 hSys[i]->SetBinContent(j,totalSys+1);
0284 hSys[i]->SetLineWidth(2);
0285 }
0286 }
0287
0288 void calcTotalSysNoUnfolding(int i) {
0289 TF1 *fNoise = new TF1("f","1+0.3*0.16*abs(1-([0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x))");
0290 fNoise->SetParameters(0.9521,0.001105,-9.397e-6,3.32e-8,-5.618e-11);
0291 hSysNoise[i] = functionHist(fNoise,hSys[i],Form("hSysNoise_cent%d",i));
0292 hSysNoise[i]->SetName(Form("hSysNoise_cent%d",i));
0293 hSysNoise[i]->SetLineColor(6);
0294 for (int j=1;j<=hSysGeneral[i]->GetNbinsX();j++) {
0295 double effSys = 0.01;
0296 hSysEff[i]->SetBinContent(j,1+effSys);
0297 hSysSmear[i]->SetBinContent(j,1.02);
0298 double JECSys = hSysJEC[i]->GetBinContent(j)-1;
0299 double SmearSys = hSysSmear[i]->GetBinContent(j)-1;
0300 double NoiseSys = hSysNoise[i]->GetBinContent(j)-1;
0301 double totalSys = sqrt( effSys * effSys +
0302 JECSys * JECSys +
0303 SmearSys * SmearSys +
0304 NoiseSys * NoiseSys
0305 );
0306 hSysGeneral[i]->SetBinContent(j,totalSys+1);
0307 hSysGeneral[i]->SetLineWidth(2);
0308 }
0309 }
0310
0311 void Draw(TH1F *h,int i) {
0312 for (int j=1;j<=hSys[i]->GetNbinsX();j++) {
0313 double val = h->GetBinContent(j);
0314 double err = hSys[i]->GetBinContent(j)-1;
0315 cout << "Sys Value Check" <<val<<" "<<err<<" "<<h->GetBinLowEdge(j)<<" "<<val*(1-err)<<" "<<h->GetBinLowEdge(j+1)<<" "<<val*(1+err)<<endl;
0316 TBox *b = new TBox(h->GetBinLowEdge(j),val*(1-err),h->GetBinLowEdge(j+1),val*(1+err));
0317
0318 b->SetFillStyle(1001);
0319
0320
0321
0322
0323 b->SetFillColor(5);
0324 b->SetLineColor(5);
0325 b->Draw();
0326 }
0327 }
0328
0329 void DrawTGraph(TGraphErrors *h,int i) {
0330 double xv;
0331 double val;
0332 for(int j=0;j<h->GetN();j++){
0333 h->GetPoint(j,xv,val);
0334 double err = hSysGeneral[i]->GetBinContent(j+1)-1;
0335
0336 TBox *b_ = new TBox(hSysGeneral[i]->GetBinLowEdge(j+1),val*(1-err),hSysGeneral[i]->GetBinLowEdge(j+2),val*(1+err));
0337 b_->SetFillColor(kGray);
0338 b_->SetFillStyle(1001);
0339 b_->SetLineColor(kGray);
0340 b_->Draw();
0341 }
0342 }
0343
0344 void DrawUnfoErr(TH1F *h,int i) {
0345 for (int j=1;j<=hSysIter[i]->GetNbinsX();j++) {
0346 double val = h->GetBinContent(j);
0347 double err = hSysIter[i]->GetBinContent(j)-1;
0348
0349 TBox *b = new TBox(h->GetBinLowEdge(j),val*(1-err),h->GetBinLowEdge(j+1),val*(1+err));
0350 b->SetFillColor(29);
0351 b->SetFillStyle(1001);
0352 b->SetLineColor(29);
0353 b->Draw();
0354 }
0355 }
0356
0357
0358 void DrawNpartSys(double yvalue,int i,double xvalue) {
0359 double yerrorNpart[6]= {0.0409, 0.0459,0.0578,0.0944, 0.143, 0.176 };
0360 double err = hSys[i]->GetBinContent(1)-1;
0361 TBox *b = new TBox(xvalue-6.,yvalue*(1-err-yerrorNpart[i]),xvalue+6.,yvalue*(1+err+yerrorNpart[i]));
0362
0363 b->SetFillColor(kGray);
0364 b->SetFillStyle(1001);
0365 b->SetLineColor(kGray);
0366
0367
0368
0369
0370
0371 b->Draw();
0372
0373 }
0374
0375
0376 void DrawComponent(int i) {
0377 calcTotalSys(i);
0378 TH1D *h = new TH1D(Form("hSysTmp_cent%d",i),"",nbins_recrebin, boundaries_recrebin);
0379 makeHistTitle(h,"","Jet p_{T} (GeV/c)","Systematic uncertainty");
0380 h->SetAxisRange(-0.25,0.4,"Y");
0381 h->Draw();
0382 TH1F* sys = drawEnvelope(hSys[i],"same",hSys[i]->GetLineColor(),1001,hSys[i]->GetLineColor(),-1);
0383 TH1F* sysIter = drawEnvelope(hSysIter[i], "same", hSysIter[i]->GetLineColor(), 3004, hSysIter[i]->GetLineColor(), -1);
0384
0385 TH1F* sysJEC = drawEnvelope(hSysJEC[i],"same",hSysJEC[i]->GetLineColor(),3005,hSysJEC[i]->GetLineColor(),-1);
0386 TH1F* sysSmear = drawEnvelope(hSysSmear[i],"same",hSysSmear[i]->GetLineColor(),3001,hSysSmear[i]->GetLineColor(),-1);
0387 TH1F* sysEff = drawEnvelope(hSysEff[i],"same",hSysEff[i]->GetLineColor(),3002,hSysEff[i]->GetLineColor(),-1);
0388 TH1F* sysNoise = drawEnvelope(hSysNoise[i],"same",hSysNoise[i]->GetLineColor(),3001,hSysNoise[i]->GetLineColor(),-1);
0389 TLine *l = new TLine(h->GetBinLowEdge(1),0,h->GetBinLowEdge(h->GetNbinsX()+1),0);
0390 l->Draw();
0391 TLine *l2 = new TLine(h->GetBinLowEdge(1),-0.25,h->GetBinLowEdge(1),0.4);
0392 l2->Draw();
0393 TLegend *leg = myLegend(0.52,0.6,0.95,0.93);
0394 leg->SetTextSize(0.043);
0395 leg->AddEntry(sys,"Total Systematics","f");
0396 leg->AddEntry(sysIter,"Unfolding","f");
0397 leg->AddEntry(sysJEC,"Jet Energy Scale","f");
0398 leg->AddEntry(sysEff,"Jet Trigger Efficiency","f");
0399 leg->AddEntry(sysSmear,"UE fluctuation","f");
0400 leg->AddEntry(sysNoise,"HCAL Noise","f");
0401 if (i==nbins_cent-1)leg->Draw();
0402 }
0403
0404
0405 };
0406
0407 class JetData
0408 {
0409 public:
0410 JetData(char *fileName, char *jetTree, char *genJetTree, bool loadGenJet = 1) {
0411 cout <<"Open "<<fileName<<endl;
0412 tFile = new TFile(fileName,"read");
0413 tEvt = (TTree*)tFile->Get("hiEvtAnalyzer/HiTree");
0414 tJet = (TTree*)tFile->Get(jetTree);
0415 tJet->SetBranchAddress("jtpt" , jtpt );
0416 tJet->SetBranchAddress("trackMax" , trackMax );
0417 tJet->SetBranchAddress("chargedMax" , chargedMax );
0418 tJet->SetBranchAddress("refpt", refpt);
0419 tJet->SetBranchAddress("nref" ,&njets);
0420 tJet->SetBranchAddress("jteta", jteta);
0421 tJet->SetBranchAddress("pthat",&pthat);
0422 if (loadGenJet) tGenJet = (TTree*)tFile->Get(genJetTree);
0423 if (loadGenJet) tGenJet->SetBranchAddress("ngen" ,&ngen);
0424 if (loadGenJet) tGenJet->SetBranchAddress("genpt", genpt);
0425 if (loadGenJet) tGenJet->SetBranchAddress("gensubid", gensubid);
0426 tEvt->SetBranchAddress("hiBin",&bin );
0427 tEvt->SetBranchAddress("vz",&vz );
0428 tJet->AddFriend(tEvt);
0429 };
0430 TFile *tFile;
0431 TTree *tJet;
0432 TTree *tGenJet;
0433 TTree *tEvt;
0434 float jtpt[1000];
0435 float refpt[1000];
0436 float jteta[1000];
0437 float trackMax[1000];
0438 float chargedMax[1000];
0439 float genpt[1000];
0440 int gensubid[1000];
0441 float vz;
0442 float pthat;
0443 int njets;
0444 int ngen;
0445 int bin;
0446 };
0447
0448 class UnfoldingHistos
0449 {
0450 public:
0451 UnfoldingHistos(int i) {
0452 hResTrue = new TH1F(Form("Restrue_cent%d",i), Form(" RooUnfold Reconstructed Truth cent%d",i) , nbins_truth, boundaries_truth);
0453 hResMeas = new TH1F(Form("Resmeas_cent%d",i), Form(" RooUnfold ResMeasured_cent1%d",i), nbins_rec, boundaries_rec);
0454 hGen = new TH1F(Form("hGen_cent%d",i) , Form(" Generator Truth_cent%d",i) , nbins_truth, boundaries_truth);
0455 hRecoMC = new TH1F(Form("hRecoMC_cent%d",i) , Form(" Generator Reco_cent%d",i) , nbins_truth, boundaries_truth);
0456 hMatrix = new TH2F(Form("hMatrix_cent%d",i), Form(" Response Matrix;Genjet^{p_{T}};Recojet^{p_{T}}",i), nbins_truth, boundaries_truth,nbins_rec, boundaries_rec);
0457 hMatrixRebin = new TH2F(Form("hMatrixRebin_cent%d",i), Form(" Response Matrix;Genjet^{p_{T}};Recojet^{p_{T}}",i), nbins_recrebinM, boundaries_recrebinM,nbins_recrebinM, boundaries_recrebinM);
0458 hMeas = new TH1F(Form("hMeas_cent%d",i) , Form("Measured jet p_{T} spectra_cent%d",i) ,nbins_rec,boundaries_rec);
0459
0460
0461
0462
0463
0464
0465
0466
0467
0468
0469
0470
0471 hMeasJECSys = new TH1F(Form("hMeasJECSys_cent%d",i) , Form(" Measured jet p_{T} sectra for JEC systematics_cent%d",i) ,nbins_rec,boundaries_rec);
0472 hMeasSmearSys = new TH1F(Form("hMeasSmearSys_cent%d",i) , Form(" Measured jet p_{T} sectra for Smear systematics_cent%d",i) ,nbins_rec,boundaries_rec);
0473 hMeas ->Sumw2();
0474
0475
0476
0477
0478
0479
0480
0481
0482 hMeasJECSys ->Sumw2();
0483 hMeasSmearSys ->Sumw2();
0484
0485
0486
0487
0488
0489
0490
0491
0492 hResTrue ->Sumw2();
0493 hGen ->Sumw2();
0494 hRecoMC ->Sumw2();
0495 hResMeas ->Sumw2();
0496 hMatrix ->Sumw2();
0497 hMatrixRebin ->Sumw2();
0498 };
0499 TH1F* hGen;
0500
0501 TH1F* hRecoMC;
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520
0521
0522 TH1F* hMCMeas;
0523 TH1F* hMeas;
0524 TH1F* hMeasJECSys;
0525 TH1F* hMeasSmearSys;
0526 TH1F* hResTrue;
0527 TH1F* hResMeas;
0528 TH1F* hReco;
0529 TH1F* hRecoBinByBin;
0530 TH1F* hRecoJECSys;
0531 TH1F* hRecoSmearSys;
0532 TH1F* hRecoIterSys[10];
0533 TH2F* hMatrix;
0534 TH2F* hMatrixRebin;
0535 TH2F* hMatrixNorm;
0536 TH2F* hResponse;
0537 TH2F* hResponseNorm;
0538 TH1F *hMeasMatch;
0539 };
0540
0541
0542 void cleanup(TH1F *h)
0543 {
0544 for (int i=1;i<=h->GetNbinsX();i++)
0545 {
0546 double val1 = h->GetBinContent(i);
0547 double valErr1 = h->GetBinError(i);
0548 if (valErr1>=val1) {
0549 h->SetBinContent(i,0);
0550 h->SetBinError(i,0);
0551 }
0552 }
0553
0554 }
0555
0556
0557 void removeError(TH1F *h)
0558 {
0559 for (int i=1;i<=h->GetNbinsX();i++)
0560 {
0561 h->SetBinError(i,0);
0562 }
0563
0564 }
0565
0566
0567 void removeZero(TH1 *h)
0568 {
0569 double min = 0;
0570 for(int i = 1;i<h->GetNbinsX();i++){
0571 if(h->GetBinContent(i)>min&&h->GetBinContent(i)>0)
0572 min = h->GetBinContent(i);
0573 }
0574
0575 for(int i = 1;i<h->GetNbinsX();i++){
0576 if(h->GetBinContent(i) == 0){
0577 h->SetBinContent(i,min/10.);
0578 h->SetBinError(i,min/10.);
0579 }
0580 }
0581 }
0582
0583
0584
0585 TH1F *rebin(TH1F *h, char *histName)
0586 {
0587 TH1F *hRebin = new TH1F(Form("%s_rebin",h->GetName()),Form("rebin %s",h->GetTitle()),nbins_recrebin,boundaries_recrebin);
0588 for (int i=1;i<=h->GetNbinsX();i++)
0589 {
0590 double val=h->GetBinContent(i);
0591 double valErr=h->GetBinError(i);
0592 int binNum = hRebin->FindBin(h->GetBinCenter(i));
0593 double val1 = hRebin->GetBinContent(binNum);
0594 double valErr1 = hRebin->GetBinError(binNum);
0595 hRebin->SetBinContent(binNum,val+val1);
0596 hRebin->SetBinError(binNum,sqrt(valErr1*valErr1+valErr*valErr));
0597 }
0598 cleanup(hRebin);
0599 hRebin->SetName(histName);
0600 return hRebin;
0601 }
0602
0603
0604 TH1F *rebin2(TH1F *h, char *histName)
0605 {
0606 TH1F *hRebin = new TH1F(Form("%s_rebin",h->GetName()),Form("rebin %s",h->GetTitle()),nbins_rec,boundaries_rec);
0607 for (int i=1;i<=h->GetNbinsX();i++)
0608 {
0609 double val=h->GetBinContent(i);
0610 double valErr=h->GetBinError(i);
0611 int binNum = hRebin->FindBin(h->GetBinCenter(i));
0612 double val1 = hRebin->GetBinContent(binNum);
0613 double valErr1 = hRebin->GetBinError(binNum);
0614 hRebin->SetBinContent(binNum,val+val1);
0615 hRebin->SetBinError(binNum,sqrt(valErr1*valErr1+valErr*valErr));
0616 }
0617 cleanup(hRebin);
0618 hRebin->SetName(histName);
0619 return hRebin;
0620 }
0621
0622
0623
0624 TH1F *rebin_Npart(TH1F *h, char *histName)
0625 {
0626 TH1F *hRebin = new TH1F(Form("%s_rebin",h->GetName()),Form("rebin %s",h->GetTitle()),nbins_recrebin_Npart,boundaries_recrebin_Npart);
0627 for (int i=1;i<=h->GetNbinsX();i++)
0628 {
0629 double val=h->GetBinContent(i);
0630 double valErr=h->GetBinError(i);
0631 int binNum = hRebin->FindBin(h->GetBinCenter(i));
0632 double val1 = hRebin->GetBinContent(binNum);
0633 double valErr1 = hRebin->GetBinError(binNum);
0634 hRebin->SetBinContent(binNum,val+val1);
0635 hRebin->SetBinError(binNum,sqrt(valErr1*valErr1+valErr*valErr));
0636 }
0637 cleanup(hRebin);
0638 hRebin->SetName(histName);
0639 return hRebin;
0640 }
0641
0642
0643 void divideBinWidth(TH1 *h)
0644 {
0645 h->Sumw2();
0646 for (int i=0;i<=h->GetNbinsX();i++)
0647 {
0648 Float_t val = h->GetBinContent(i);
0649 Float_t valErr = h->GetBinError(i);
0650 val/=h->GetBinWidth(i);
0651 valErr/=h->GetBinWidth(i);
0652 h->SetBinContent(i,val);
0653 h->SetBinError(i,valErr);
0654 }
0655 h->GetXaxis()->CenterTitle();
0656 h->GetYaxis()->CenterTitle();
0657 }
0658
0659
0660
0661 void checkMaximumSys(TH1F *hSys, TH1F *h, int opt=0,double minVal = 1)
0662 {
0663 if (h->GetNbinsX()!=hSys->GetNbinsX()) {
0664 cout <<"ERROR! Different NBins in subroutine checkMaximumSys!"<<endl;
0665 } else {
0666 double val = minVal;
0667 for (int i=1;i<=h->GetNbinsX();i++) {
0668
0669 if (h->GetBinContent(i)==0) continue;
0670 if (opt==0) val=minVal;
0671 if (fabs(hSys->GetBinContent(i))>val) val = fabs(hSys->GetBinContent(i));
0672 if (fabs(h->GetBinContent(i)-1)+1>val) val=fabs(h->GetBinContent(i)-1)+1;
0673 hSys->SetBinContent(i,val);
0674 }
0675 }
0676 }
0677
0678
0679
0680 void makeMultiPanelCanvasWithGap(TCanvas*& canv,
0681 const Int_t columns,
0682 const Int_t rows,
0683 const Float_t leftOffset,
0684 const Float_t bottomOffset,
0685 const Float_t leftMargin,
0686 const Float_t bottomMargin,
0687 const Float_t edge, const Float_t asyoffset) {
0688 if (canv==0) {
0689 Error("makeMultiPanelCanvas","Got null canvas.");
0690 return;
0691 }
0692 canv->Clear();
0693
0694 TPad* pad[columns][rows];
0695
0696 Float_t Xlow[columns];
0697 Float_t Xup[columns];
0698 Float_t Ylow[rows];
0699 Float_t Yup[rows];
0700
0701 Float_t PadWidth =
0702 (1.0-leftOffset)/((1.0/(1.0-leftMargin)) +
0703 (1.0/(1.0-edge))+(Float_t)columns-2.0);
0704 Float_t PadHeight =
0705 (1.0-bottomOffset)/((1.0/(1.0-bottomMargin)) +
0706 (1.0/(1.0-edge))+(Float_t)rows-2.0);
0707
0708
0709
0710 Xlow[0] = leftOffset;
0711 Xup[0] = leftOffset + PadWidth/(1.0-leftMargin);
0712 Xup[columns-1] = 1;
0713 Xlow[columns-1] = 1.0-PadWidth/(1.0-edge);
0714
0715 Yup[0] = 1;
0716 Ylow[0] = 1.0-PadHeight/(1.0-edge);
0717 Ylow[rows-1] = bottomOffset;
0718 Yup[rows-1] = bottomOffset + PadHeight/(1.0-bottomMargin);
0719
0720 for(Int_t i=1;i<columns-1;i++) {
0721 Xlow[i] = Xup[0] + (i-1)*PadWidth;
0722 Xup[i] = Xup[0] + (i)*PadWidth;
0723 }
0724 Int_t ct = 0;
0725 for(Int_t i=rows-2;i>0;i--) {
0726 if(i==rows-2){
0727 Ylow[i] = Yup[rows-1] + ct*PadHeight;
0728 Yup[i] = Yup[rows-1] + (ct+1)*PadHeight;
0729 }else{
0730 Ylow[i] = Yup[rows-1] + ct*PadHeight;
0731 Yup[i] = Yup[rows-1] + (ct+1)*PadHeight;
0732
0733 }
0734 ct++;
0735 }
0736
0737 TString padName;
0738 for(Int_t i=0;i<columns;i++) {
0739 for(Int_t j=0;j<rows;j++) {
0740 canv->cd();
0741 padName = Form("p_%d_%d",i,j);
0742
0743
0744
0745 if(j==0){
0746 pad[i][j] = new TPad(padName.Data(),padName.Data(),
0747 Xlow[i],Ylow[j]-asyoffset,Xup[i],Yup[j]);
0748 }else{
0749 pad[i][j] = new TPad(padName.Data(),padName.Data(),
0750 Xlow[i],Ylow[j],Xup[i],Yup[j]-asyoffset);
0751 }
0752
0753
0754 if(i==0) pad[i][j]->SetLeftMargin(leftMargin);
0755 else pad[i][j]->SetLeftMargin(0);
0756
0757 if(i==(columns-1)) pad[i][j]->SetRightMargin(edge);
0758 else pad[i][j]->SetRightMargin(0);
0759
0760 if(j==0) pad[i][j]->SetTopMargin(edge);
0761
0762 else pad[i][j]->SetTopMargin(0.02);
0763
0764
0765
0766
0767 if(j==(rows-1)) pad[i][j]->SetBottomMargin(bottomMargin);
0768 else pad[i][j]->SetBottomMargin(0.15);
0769
0770 pad[i][j]->Draw();
0771 pad[i][j]->cd();
0772 pad[i][j]->SetNumber(columns*j+i+1);
0773
0774 }
0775 }
0776 }
0777
0778 void putCMSPrel(double x, double y, double size){
0779 TLatex *tex=0;
0780 tex = new TLatex(x,y,"CMS Preliminary");
0781 tex->SetTextSize(size);
0782 tex->SetLineWidth(2);
0783 tex->SetNDC();
0784 tex->Draw();
0785 }
0786 void drawText(const char *text, float xp, float yp, int size){
0787 TLatex *tex = new TLatex(xp,yp,text);
0788 tex->SetTextFont(63);
0789 tex->SetTextSize(size);
0790 tex->SetTextColor(kBlack);
0791 tex->SetLineWidth(1);
0792
0793 tex->SetNDC();
0794 tex->Draw();
0795 }
0796
0797
0798 void prepareNcollUnc(int nbins, float maxpt=300.){
0799
0800 int fillsty = 1001;
0801
0802 const int n = nbins;
0803
0804 double xvalue[n];
0805 double yvalue[n];
0806 double xerror[n];
0807 double yerror1[n], yerror2[n], yerror3[n], yerror4[n], yerror5[n], yerror6[n];
0808
0809
0810
0811 for(int i=0;i<nbins;i++){
0812
0813 xvalue[i] = 300.1 + 1.2*(double)i, yvalue[i]=1.0, xerror[i]=0.0;
0814
0815
0816 yerror1[i]=0.0409, yerror2[i]=0.0459, yerror3[i]=0.0578, yerror4[i]=0.0944, yerror5[i]=0.143, yerror6[i]=0.176;
0817
0818
0819
0820 yerror1[i]=TMath::Sqrt(yerror1[i]*yerror1[i]+0.06*0.06);
0821 yerror2[i]=TMath::Sqrt(yerror2[i]*yerror2[i]+0.06*0.06);
0822 yerror3[i]=TMath::Sqrt(yerror3[i]*yerror3[i]+0.06*0.06);
0823 yerror4[i]=TMath::Sqrt(yerror4[i]*yerror4[i]+0.06*0.06);
0824 yerror5[i]=TMath::Sqrt(yerror5[i]*yerror5[i]+0.06*0.06);
0825 yerror6[i]=TMath::Sqrt(yerror6[i]*yerror6[i]+0.06*0.06);
0826 cout<<"TAA + Lumi uncertainty = "<<yerror1[i]<<endl;
0827
0828 }
0829
0830
0831
0832 int ci = 15;
0833
0834 tTAAerr[0] = new TGraphErrors(n,xvalue,yvalue,xerror,yerror1);
0835 tTAAerr[0]->SetFillColor(ci);
0836 tTAAerr[0]->SetLineColor(ci);
0837 tTAAerr[0]->SetFillStyle(fillsty);
0838
0839 tTAAerr[1] = new TGraphErrors(n,xvalue,yvalue,xerror,yerror2);
0840 tTAAerr[1]->SetFillColor(ci);
0841 tTAAerr[1]->SetFillStyle(fillsty);
0842
0843 tTAAerr[2] = new TGraphErrors(n,xvalue,yvalue,xerror,yerror3);
0844 tTAAerr[2] ->SetFillColor(ci);
0845 tTAAerr[2] ->SetFillStyle(fillsty);
0846
0847 tTAAerr[3] = new TGraphErrors(n,xvalue,yvalue,xerror,yerror4);
0848 tTAAerr[3]->SetFillColor(ci);
0849 tTAAerr[3]->SetFillStyle(fillsty);
0850
0851 tTAAerr[4] = new TGraphErrors(n,xvalue,yvalue,xerror,yerror5);
0852 tTAAerr[4]->SetFillColor(ci);
0853 tTAAerr[4]->SetFillStyle(fillsty);
0854
0855 tTAAerr[5] = new TGraphErrors(n,xvalue,yvalue,xerror,yerror6);
0856 tTAAerr[5]->SetFillColor(ci);
0857 tTAAerr[5]->SetFillStyle(fillsty);
0858
0859
0860
0861 }
0862
0863
0864
0865
0866
0867
0868
0869
0870
0871
0872
0873
0874
0875
0876
0877
0878
0879
0880
0881
0882
0883
0884
0885
0886
0887 void dumpDatatoTxt(const char *centbin,TH1F *h, TH1F *hsys, TH1F *htotStat, const char *txtfile)
0888 {
0889 ofstream outf(txtfile,ios::out);
0890 for(int ix=1;ix<=h->GetNbinsX();ix++){
0891 double pt = h->GetBinCenter(ix);
0892 double val = h->GetBinContent(ix);
0893 double Uncorerr = h->GetBinError(ix);
0894 double syserr = hsys->GetBinContent(ix)-1;
0895 double totStaterr = htotStat->GetBinError(ix);
0896
0897 outf<<setprecision(0) << fixed <<pt<<"\t" << setprecision(3) << fixed <<val<<"\t" << setprecision(5) << fixed << totStaterr<<"\t"<< setprecision(4) << fixed << syserr*val << endl;
0898 }
0899 outf.close();
0900 }
0901
0902
0903 TGraphErrors *ShiftGraph(TGraphErrors* pGraph, Double_t pNumber){
0904
0905 TGraphErrors *pGraphtmp;
0906 for (Int_t i=0;i<pGraph->GetN();i++){
0907 Double_t x,y;
0908 Double_t yerr;
0909 pGraph->GetPoint(i,x,y);
0910 yerr = pGraph->GetErrorY(i);
0911 x = x+pNumber;
0912 pGraphtmp->SetPoint(i,x,y);
0913 pGraphtmp->SetPointError(i,yerr);
0914 }
0915 return pGraphtmp;
0916 }
0917
0918
0919 TGraphErrors* HistToTgraphShift(TH1F* hist,Double_t pNumber){
0920
0921 TGraphErrors *pGraphtmp;
0922 int nbins = hist->GetNbinsX();
0923
0924 const int nlines = nbins;
0925
0926 float pt[nlines], xsec[nlines];
0927 float pterr[nlines], xsecerr[nlines];
0928
0929 for(int i = 0; i<nbins; i++ ){
0930 pt[i] = hist->GetBinCenter(i+1);
0931 xsec[i] = hist->GetBinContent(i+1);
0932 xsecerr[i] = hist->GetBinError(i+1);
0933 pterr[i] = 0;
0934 }
0935
0936 pGraphtmp = new TGraphErrors(nlines,pt,xsec,pterr,xsecerr);
0937
0938 for (Int_t i=0;i<pGraphtmp->GetN();i++){
0939 Double_t x,y;
0940 Double_t yerr;
0941 pGraphtmp->GetPoint(i,x,y);
0942 yerr = pGraphtmp->GetErrorY(i);
0943 x = x+pNumber;
0944 pGraphtmp->SetPoint(i,x,y);
0945 pGraphtmp->SetPointError(i,yerr);
0946 }
0947 return pGraphtmp;
0948 }
0949
0950
0951 void DrawPanelLabel(int i){
0952 TLatex *tex;
0953
0954 if(i==0) tex = new TLatex(0.07,0.9,"(f)");
0955 if(i==1) tex = new TLatex(0.07,0.9,"(e)");
0956 if(i==2) tex = new TLatex(0.19,0.9,"(d)");
0957 if(i==3) tex = new TLatex(0.07,0.9,"(c)");
0958 if(i==4) tex = new TLatex(0.07,0.9,"(b)");
0959 if(i==5) tex = new TLatex(0.19,0.9,"(a)");
0960 tex->SetTextFont(63);
0961 tex->SetTextSize(18);
0962 tex->SetTextColor(kBlack);
0963 tex->SetLineWidth(1);
0964 tex->SetNDC();
0965
0966
0967 tex->Draw();
0968
0969 }