Back to home page

sPhenix code displayed by LXR

 
 

    


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 // Truth binning
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   static const int nbins_truth = 40;
0040   static const double boundaries_truth[nbins_truth+1] = {
0041   100, 110, 120, 130, 140, 150, 160, 170,
0042   180, 190, 200, 210, 220, 230, 240, 250, 260, 270, 280, 290, 300, 310,
0043   320, 330, 340, 350, 360, 370, 380, 390, 400, 410, 420, 430, 440, 450,
0044   460, 470, 480, 490, 500
0045   };
0046 */
0047 // Measurement binning
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   static const int nbins_rec = 40;
0059   static const double boundaries_rec[nbins_rec+1] = {
0060   100, 110, 120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 220,
0061   230, 240, 250, 260, 270, 280, 290, 300, 310, 320, 330, 340, 350, 360,
0062   370, 380, 390, 400, 410, 420, 430, 440, 450, 460, 470, 480, 490, 500
0063   };
0064 */
0065 
0066 /*
0067   static const int nbins_rec = 20;
0068   static const double boundaries_rec[nbins_rec+1] = {
0069   0,5,10,15,20,25,30,35,40,45,50,55,60,70,80,90,110,130,170,250,300
0070   };
0071 
0072   static const int nbins_truth = 20;
0073   static const double boundaries_truth[nbins_truth+1] = {
0074   0,5,10,15,20,25,30,35,40,45,50,55,60,70,80,90,110,130,170,250,300
0075   };
0076 */
0077 //changed it for yaxian's comparison
0078 
0079 //static const int nbins_recrebin = 20;
0080 //static const double boundaries_recrebin[nbins_recrebin+1] = {
0081 //  0,5,10,15,20,25,30,35,40,45,50,55,60,70,80,90,110,130,170,250,330
0082 //};
0083 
0084 //temporary nbins_recrebin to accomodate the old binning style for the NLO comparison
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 //the original bin values. 
0096 //static const int nbins_recrebin = 18;
0097 //static const double boundaries_recrebin[nbins_recrebin+1] = {
0098 //30, 40, 50, 60, 70, 80, 90,100, 110, 120, 130, 140, 150, 160,
0099 //170, 180, 200, 240, 300
0100 //};
0101 
0102 //static const int nbins_recrebinM = 18;
0103 //static const double boundaries_recrebinM[nbins_recrebinM+1] = {
0104 //30, 40, 50, 60, 70, 80, 90,100, 110, 120, 130, 140, 150, 160,
0105 //170, 180, 200, 240, 300
0106 //};
0107 
0108 static const int nbins_recrebin_Npart = 1;
0109 static const double boundaries_recrebin_Npart[nbins_recrebin_Npart+1] = {
0110   100, 300//changed from 100 to 30
0111 };
0112 
0113 static const int nColor = 5;
0114 static const int colorCode[nColor] = {
0115   1, 2, kGreen+1, 4, 6
0116 };
0117 
0118 // Algos, following Ying's convension
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 // Centrality binning
0127 const int nbins_cent=     6;
0128 Double_t boundaries_cent[nbins_cent+1] = {   0,2,4,12,20,28,36   };// multiply by 2.5 to get your actual centrality % (old 2011 data) 
0129 Double_t ncoll[nbins_cent] = { 1660, 1310, 745, 251, 62.8, 10.8 };
0130 //the following is due to request form Yaxian 
0131 //const int nbins_cent = 1;
0132 //Double_t boundaries_cent[nbins_cent+1] = {0,40};
0133 //Double_t ncoll[nbins_cent] = {362.24}; //use taa instead of ncoll. 
0134 //Double_t TAA[nbins_cent] = {};
0135 
0136 // TAA 
0137 
0138 TGraphErrors *tTAAerr[6]={0};
0139 TGraphErrors *tTAAerrNpart=0;
0140 
0141 
0142 //static const double trigEffInc[nbins_recrebin] = {
0143 //  0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.966051, 0.99591, 0.999665, 1., 1., 1., 1.
0144 //};
0145 
0146 
0147 // make a histogram from TF1 function
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 // draw envelope using a systematic uncertainty histogram
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   //   hMirror->SetLineStyle(2);
0183   //   h->SetLineStyle(2);
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       //b->SetFillColor(kGray);
0318       b->SetFillStyle(1001);
0319       //b->SetLineColor(kGray);
0320       
0321       
0322       //***********For Gunther's Color Systematics Band Peference
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       //cout <<"value" <<val<<" "<<err<<" "<<hSysGeneral[i]->GetBinLowEdge(j+1)<<" "<<val*(1-err)<<" "<<hSysGeneral[i]->GetBinLowEdge(j+2)<<" "<<val*(1+err)<<endl;
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       //cout <<"value" << val<<" "<<err<<" "<<h->GetBinLowEdge(j)<<" "<<val*(1-err)<<" "<<h->GetBinLowEdge(j+1)<<" "<<val*(1+err)<<endl;
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     //cout << "value " << yvalue<<" err   "<<err<<" xvalue  "<<xvalue<<" "<<yvalue*(1-err)<<" "<<yvalue*(1+err)<<endl;
0363     b->SetFillColor(kGray);
0364     b->SetFillStyle(1001);
0365     b->SetLineColor(kGray);
0366     
0367     //***********For Gunther's Color Systematics Band Peference
0368     //b->SetFillColor(5);
0369     //b->SetLineColor(5);
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  // TH1F* sysIter = drawEnvelope(hSysIter[i],"same",hSysIter[i]->GetLineColor(),3004,hSysIter[i]->GetLineColor(),-1);
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); //Gen
0453     hResMeas       = new TH1F(Form("Resmeas_cent%d",i), Form(" RooUnfold ResMeasured_cent1%d",i),    nbins_rec, boundaries_rec); //Reco
0454     hGen           = new TH1F(Form("hGen_cent%d",i) , Form(" Generator Truth_cent%d",i)    ,    nbins_truth, boundaries_truth); //Gen
0455     hRecoMC           = new TH1F(Form("hRecoMC_cent%d",i) , Form(" Generator Reco_cent%d",i)    ,    nbins_truth, boundaries_truth); //RecoMC
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); //Matrix
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); //Matrix Rebin
0458     hMeas          = new TH1F(Form("hMeas_cent%d",i)       , Form("Measured jet p_{T} spectra_cent%d",i)       ,nbins_rec,boundaries_rec);
0459     //hCombined      = new TH1F(Form("hCombined_cent%d",i)   , Form("Combined measured jet p_{T} spectra_cent%d",i)       ,nbins_rec,boundaries_rec);
0460     //hMeas_80        = new TH1F(Form("hMeas_80_cent%d",i)       , Form("Measured jet p_{T} from HLT80 spectra_cent%d",i)       ,nbins_rec,boundaries_rec);
0461     //hMeas_65        = new TH1F(Form("hMeas_65_cent%d",i)       , Form("Measured jet p_{T} from HLT65 spectra_cent%d",i)       ,nbins_rec,boundaries_rec);
0462     //hMeas_55        = new TH1F(Form("hMeas_55_cent%d",i)       , Form("Measured jet p_{T} from HLT55 spectra_cent%d",i)       ,nbins_rec,boundaries_rec);
0463     //hMeas80JECSys    = new TH1F(Form("hMeas80JECSys_cent%d",i) , Form(" Measured jet p_{T} 80 sectra for JEC systematics_cent%d",i)       ,nbins_rec,boundaries_rec);
0464     //hCombinedJECSys    = new TH1F(Form("hCombinedJECSys_cent%d",i) , Form(" Combined Measured jet p_{T} sectra for JEC systematics_cent%d",i)       ,nbins_rec,boundaries_rec);
0465     //hMeas80SmearSys    = new TH1F(Form("hMeas80SmearSys_cent%d",i) , Form(" Measured jet p_{T}80 sectra for Smear systematics_cent%d",i)       ,nbins_rec,boundaries_rec);
0466     //hCombinedSmearSys    = new TH1F(Form("hCombinedSmearSys_cent%d",i) , Form(" Combined Measured jet p_{T} sectra for Smear systematics_cent%d",i)       ,nbins_rec,boundaries_rec);
0467     //hMeas65JECSys    = new TH1F(Form("hMeas65JECSys_cent%d",i) , Form(" Measured jet p_{T}65 sectra for JEC systematics_cent%d",i)       ,nbins_rec,boundaries_rec);
0468     //hMeas65SmearSys    = new TH1F(Form("hMeas65SmearSys_cent%d",i) , Form(" Measured jet p_{T}65 sectra for Smear systematics_cent%d",i)       ,nbins_rec,boundaries_rec);
0469     //hMeas55JECSys    = new TH1F(Form("hMeas55JECSys_cent%d",i) , Form(" Measured jet p_{T}55 sectra for JEC systematics_cent%d",i)       ,nbins_rec,boundaries_rec);
0470     //hMeas55SmearSys    = new TH1F(Form("hMeas55SmearSys_cent%d",i) , Form(" Measured jet p_{T}55 sectra for Smear systematics_cent%d",i)       ,nbins_rec,boundaries_rec);
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       hMeas_80        ->Sumw2();
0476       hMeas_65        ->Sumw2();
0477       hMeas_55        ->Sumw2();
0478       hCombined      ->Sumw2();
0479       hCombinedJECSys->Sumw2();
0480       hCombinedSmearSys->Sumw2();
0481     */
0482     hMeasJECSys    ->Sumw2();
0483     hMeasSmearSys  ->Sumw2();
0484     /*
0485       hMeas80JECSys    ->Sumw2();
0486       hMeas80SmearSys  ->Sumw2();
0487       hMeas65JECSys    ->Sumw2();
0488       hMeas65SmearSys  ->Sumw2();
0489       hMeas55JECSys    ->Sumw2();
0490       hMeas55SmearSys  ->Sumw2();
0491     */
0492     hResTrue       ->Sumw2();
0493     hGen           ->Sumw2();
0494     hRecoMC        ->Sumw2();
0495     hResMeas       ->Sumw2();
0496     hMatrix        ->Sumw2();
0497     hMatrixRebin   ->Sumw2();
0498   };
0499   TH1F* hGen;
0500   //start adding by raghav
0501   TH1F* hRecoMC;
0502   /*
0503     TH1F* hTrue;//to check with the pp data distribution
0504     TH1F* hUnfold;// make sure that unfolding works properly
0505     TH1F* hMC;
0506     //Try to define the different histograms to take in the entries for different HLT datasets here
0507     /*
0508     TH1F* hMeas_80;
0509     TH1F* hMeas_65;
0510     TH1F* hMeas_55;
0511     TH1F* hCombined;
0512     TH1F* hCombinedJECSys;
0513     TH1F* hCombinedSmearSys;
0514     TH1F* hMeas80JECSys;
0515     TH1F* hMeas80SmearSys;
0516     TH1F* hMeas65JECSys;
0517     TH1F* hMeas65SmearSys;
0518     TH1F* hMeas55JECSys;
0519     TH1F* hMeas55SmearSys;
0520   */
0521   //end adding by raghav
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 // Remove bins with error > central value
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 // Remove error 
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 // Remove Zero
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 // rebin the spectra
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 // rebin the spectra
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 // Npart rebin the spectra  
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 // divide by bin width
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 // make systematic histogram
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       //cout <<i<<" "<<val<<" "<<hSys->GetBinContent(i)<<" "<<h->GetBinContent(i)<<endl;
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   //PadHeight = 0.5*PadHeight;
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       //Yup[i] = 0.2*Yup[i];
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       //pad[i][j] = new TPad(padName.Data(),padName.Data(),
0743       //Xlow[i],Ylow[j],Xup[i],Yup[j]);
0744       // this is hacked version to create aysmmetric pads around low 
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       //else pad[i][j]->SetTopMargin(0.01);
0762       else pad[i][j]->SetTopMargin(0.02);
0763       
0764       //if(j==0) pad[i][j]->SetTopMargin(edge*3.5);
0765       //else pad[i][j]->SetTopMargin(0.0);
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   //tex->SetTextFont(42);
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     // TAA
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     // add 6% error 
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   // int ci = 29;
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   void DrawNpartTAABand(){
0864   double xvalueNpart[6];
0865   double yerrorNpart[6];
0866   xvalueNpart[0] = 381.29; xvalueNpart[1] = 329.41; xvalueNpart[2] = 224.28;
0867   xvalueNpart[3] = 108.12; xvalueNpart[4] = 42.04;  xvalueNpart[5] = 11.43;
0868   yerrorNpart[0]=0.0409, yerrorNpart[1]=0.0459, yerrorNpart[2]=0.0578, yerrorNpart[3]=0.0944, yerrorNpart[4]=0.143, yerrorNpart[5]=0.176;
0869 
0870 int ci = 30;
0871 
0872 for (int i=0;i<6;i++) {
0873 
0874 TBox *b = new TBox(xvalueNpart[i]-5,1.-yerrorNpart[i]/2,xvalueNpart[i]+5,1.+yerrorNpart[i]/2);
0875 b->SetFillColor(ci);
0876 b->SetFillStyle(3001);
0877 b->SetLineColor(ci);
0878 b->Draw();
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   // shifts a graph by the absolute value in the argument
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 }