test/DijetRatio/MassAnaBackgrd.C

0001 # root script for making Djiet Ratio plots from
0002 # output of dijet ratio analysis.  Used by Manoj Jha.
0003
0004 #include <TH1.h>
0005 #include <TFile.h>
0006 #include <TBranch.h>
0007 #include <TTree.h>
0008 #include <TLeaf.h>
0009
0010 #include <iostream>
0011 using namespace std;
0012
0013 void  scale(TH1* &h1) {
0014
0015             Int_t nbin = h1->GetNbinsX();
0016         TH1* h2 = (TH1F*) h1->Clone();
0017 //      h1->Sumw2();
0018         for (Int_t i =0; i < nbin; ++i){
0019             float binWidth = h2->GetBinWidth(i);
0020             float binContent = h2->GetBinContent(i);
0021             float value = binContent/binWidth;
0022             float error = h2->GetBinError(i);
0023             float errVal = error/binWidth;
0024             h1->SetBinContent(i, value);
0025             h1->SetBinError(i, errVal);
0026         }
0027 }
0028
0029 void Addition(TString hName, TH1F* &h0, TH1F* &h0Var){
0030
0031     //Reading cross-sections
0032     ifstream in ("datasetBackgrd.txt");
0033     if (!in){
0034         cout << "Not able to open the file \n" << endl;
0035         return 1;
0036     }
0037
0038
0039     double n1;
0040     //X-section
0041     double n2;
0042
0043     vector<double> tried ;  // Vectors for number of events
0044     vector<double> xSection;
0045     vector<double> weight;
0046
0047     tried.clear();
0048     xSection.clear();
0049     weight.clear();
0050
0051     while (in){
0052     in >> n1 >> n2 ;
0053     tried.push_back(n1);
0054     xSection.push_back(n2);
0055     weight.push_back(n2*pow(10.0,9.0)/n1);
0056     }
0057
0058     in.close();
0059
0060
0061         double nEvents = 0.;
0062     double crossSec = 0.;
0063
0064     //Number of files
0065     TFile* f[21];
0066     TString s1 = "QcdBackgrd_";
0067     TString s2 = ".root";
0068
0069     TF1 *f1 = new TF1("f1","1.0",0,10);
0070
0071     for (int i =0; i < 21; i++){
0072         int j = i+1 ;
0073         TString fileName = "file25/" + s1 + j + s2;
0074         f[i] = new TFile(fileName);
0075         double wt = weight[i];
0076         double wt2 = wt*wt;
0077
0078             if (i==0) {
0079          h0 = (TH1F*)f[i]->Get(hName)->Clone(fileName + hName + "0");
0080          h0Var = (TH1F*)f[i]->Get(hName)->Clone(fileName + hName + "Var");
0081
0082         h0->Multiply(f1, wt);
0083         h0Var->Multiply(f1, wt2);
0084         }
0085         if (i >=1){
0086         TH1F* hA0 = (TH1F*)f[i]->Get(hName)->Clone(fileName + "A0");
0087
0088         h0->Add(hA0,wt);
0089         h0Var->Add(hA0, wt2);
0090
0091         delete hA0;
0092         } // i>=1
0093
0094     } //for (i =0 ; i <21; i++)
0095
0096     int nBin = h0->GetNbinsX();
0097     for (int i = 0; i < nBin; i++){
0098 //      double valueVar1 = h0Var->GetBinError(i);
0099         double valueVar1 = h0Var->GetBinContent(i);
0100         double error = TMath::Sqrt(valueVar1);
0101         h0->SetBinError(i,error);
0102     }
0103
0104     h0 = (TH1F*)h0->Clone(hName);
0105
0106 }//Addition
0107
0108 void f(){
0109     cout << "program is working fine " << endl;
0110 }
0111
0112
0113 TH1F* Ratio(TH1F* h1, TH1F* h2){
0114
0115     TH1F* ratio = (TH1F*) h1->Clone("ratio");
0116     ratio->Divide(h2);
0117     int nbin = h1->GetNbinsX();
0118
0119     for (Int_t i = 0; i < nbin; i++){
0120         double sigma1 = h1->GetBinError(i);
0121         double sigma2 = h2->GetBinError(i);
0122
0123         double x1 = h1->GetBinContent(i);
0124         double x2 = h2->GetBinContent(i);
0125
0126         double r = ratio->GetBinContent(i);
0127         double error =0. ;
0128         if (x1 !=0. && x2 !=0.){
0129              error = r * TMath::Sqrt(TMath::Power(sigma1/x1, 2.0) + TMath::Power(sigma2/x2, 2.0));
0130         }
0131         ratio->SetBinError(i, error);
0132     }
0133
0134     return ratio;
0135 }//Ratio
0136
0137 TCanvas* fitting(TCanvas* c1, TH1F* ratio, TString jetFlavor, TFile* f){
0138
0139     char labelCons[150];
0140     char labelLine[150];
0141     char labelPoly[150];
0142
0143     ratio->SetTitle("DiJet Ratio");
0144     ratio->GetXaxis()->SetTitle("DiJet Mass (in GeV)");
0145     ratio->GetXaxis()->SetRangeUser(330.,7600.);
0146         if (!strcmp(jetFlavor,"ratioGen") || !strcmp(jetFlavor,"ratioCalo") || !strcmp(jetFlavor,"ratioCor") )
0147         ratio->GetYaxis()->SetRangeUser(0.,1.2);
0148     else
0149         ratio->GetYaxis()->SetRangeUser(0.6,1.6);
0150         if (!strcmp(jetFlavor,"ratioGen") || !strcmp(jetFlavor,"ratioCalo") || !strcmp(jetFlavor,"ratioCor") )
0151         ratio->GetYaxis()->SetTitle("Ratio=N(|#eta|<0.5)/N(0.5<|#eta|<1)");
0152     else
0153         ratio->GetYaxis()->SetTitle("Ratio(CorrectedJets/GenJets)");
0154
0155     ratio->GetYaxis()->SetTitleOffset(1.3);
0156     ratio->SetStats(kFALSE);
0157     ratio->Draw("E1P");
0158
0159     TF1* f1 = new TF1("f1", "pol0");
0160     f1->SetLineColor(1);
0161     ratio->Fit("f1");
0162     TF1* fit = ratio->GetFunction("f1");
0163     float constChi2 = fit->GetChisquare();
0164     float constP0 = fit->GetParameter(0);
0165     float constP0Er = fit->GetParError(0);
0166     int constNDF = fit->GetNDF();
0167     float constChi = constChi2/constNDF;
0168     sprintf(labelCons, "y = %5.3f #pm %5.3f, #chi^{2} = %5.3f, NDF = %d, #chi^{2}/NDF = %5.3f", constP0, constP0Er, constChi2, constNDF, constChi);
0169
0170
0171     TF1* f2 = new TF1("f2", "pol1");
0172     f2->SetLineColor(3);
0173     ratio->Fit("f2", "+");
0174     f2->Draw("SAME");
0175     TF1* fit = ratio->GetFunction("f2");
0176     float lineChi2 = fit->GetChisquare();
0177     float lineP0 = fit->GetParameter(0);
0178     float lineP1 = fit->GetParameter(1);
0179     int lineNDF = fit->GetNDF();
0180     float lineChi = lineChi2/lineNDF;
0181     sprintf(labelLine, "y = %5.3e*x + %5.3f, #chi^{2} = %5.3f, NDF = %d, #chi^{2}/NDF = %5.3f", lineP1, lineP0, lineChi2, lineNDF, lineChi);
0182
0183     TF1* f3 = new TF1("f3", "pol2");
0184     f3->SetLineColor(4);
0185     ratio->Fit("f3", "+");
0186     TF1* fit = ratio->GetFunction("f3");
0187     float polyChi2 = fit->GetChisquare();
0188     float polyP0 = fit->GetParameter(0);
0189     float polyP1 = fit->GetParameter(1);
0190     float polyP2 = fit->GetParameter(2);
0191     int polyNDF = fit->GetNDF();
0192     float polyChi = polyChi2/polyNDF;
0193     sprintf(labelPoly, "y = %5.3e*x^{2} + %5.3e*x + %5.3f, #chi^{2} = %5.3f, NDF = %d, #chi^{2}/NDF = %5.3f",polyP2, polyP1, polyP0, polyChi2, polyNDF, polyChi);
0194
0195         if (!strcmp(jetFlavor,"ratioGen"))
0196     TLegend *leg = new TLegend(0.25,0.7,0.9,0.9, "Generated Jets");
0197     else if (!strcmp(jetFlavor,"ratioCalo"))
0198     TLegend *leg = new TLegend(0.25,0.7,0.9,0.9, "Calorimetry Jets");
0199     else if (!strcmp(jetFlavor,"ratioCor"))
0200     TLegend *leg = new TLegend(0.25,0.7,0.9,0.9, "Corrected Jets");
0201     else if (!strcmp(jetFlavor,"R0_CorGen"))
0202     TLegend *leg = new TLegend(0.25,0.7,0.9,0.9, "|#eta| < 0.5");
0203     else
0204     TLegend *leg = new TLegend(0.25,0.7,0.9,0.9, "|#eta| < 1.0");
0205
0206     leg->AddEntry(f1,labelCons,"L");
0207     leg->AddEntry(f2,labelLine,"L");
0208     leg->AddEntry(f3,labelPoly,"L");
0209     leg->Draw();
0210     //save ratio histo into the root file
0211         if (!strcmp(jetFlavor,"ratioGen")){
0212         TH1F* ratioGen = (TH1F*)ratio->Clone("ratioGen");
0213         f->cd();
0214         ratioGen->Write();
0215     }
0216         else if (!strcmp(jetFlavor,"ratioCalo")){
0217         TH1F* ratioCalo = (TH1F*)ratio->Clone("ratioCalo");
0218         f->cd();
0219         ratioCalo->Write();
0220     }
0221         if (!strcmp(jetFlavor,"ratioCor")){
0222         TH1F* ratioCor = (TH1F*)ratio->Clone("ratioCor");
0223         f->cd();
0224         ratioCor->Write();
0225     }
0226
0227     c1->Print("c1.ps");
0228
0229     return c1;
0230 }// fitting
0231
0232 TCanvas* plot(TCanvas* c1, TString id, TFile* f){
0233
0234     c1 = new TCanvas("c1", "c1",3,25,999,799);
0235     c1->SetLogy(); c1->SetGridx(); c1->SetGridy();
0236
0237     TH1F* hGen; TH1F* hGenVar;
0238     TH1F* hCalo;    TH1F* hCaloVar;
0239     TH1F* hCor; TH1F* hCorVar;
0240
0241     TString sGen = "hGen" + id;
0242     TString sCalo = "hCalo" + id;
0243     TString sCor = "hCor" + id;
0244
0245     //title
0246     TString title = "DiJet Mass Distribution";
0247
0248     Addition(sGen, hGen, hGenVar);
0249     f->cd();
0250     hGen->Write();
0251     Addition(sCalo, hCalo, hCaloVar);
0252     f->cd();
0253     hCalo->Write();
0254     Addition(sCor, hCor, hCorVar);
0255     f->cd();
0256     hCor->Write();
0257
0258     float GenMax = hGen->GetMaximumStored();
0259     float CaloMax = hCalo->GetMaximumStored();
0260     float CorMax = hCor->GetMaximumStored();
0261
0262     float max1 = TMath::Max(GenMax, CaloMax);
0263     float max = TMath::Max(max1, CorMax);
0264     float maxY = max + 1000.;
0265
0266     hGen->GetXaxis()->SetTitle("DiJet Mass (in GeV)");
0267     hGen->GetXaxis()->SetRangeUser(330.,7600.);
0268 //  hGen->GetYaxis()->SetRangeUser(0.,maxY);
0269     hGen->GetYaxis()->SetTitle("d#sigma/dM_{diJet} (pb/GeV)");
0270     hGen->GetYaxis()->SetTitleOffset(1.3);
0271     hGen->SetMarkerStyle(20);
0272     hGen->SetMarkerColor(1);
0273     hGen->SetLineColor(1);
0274     hGen->SetStats(kFALSE);
0275     scale(hGen);
0276     hGen->Draw("EP");
0277
0278     hCalo->SetMarkerStyle(21);
0279     hCalo->SetMarkerColor(2);
0280     hCalo->SetLineColor(2);
0281     hCalo->SetStats(kFALSE);
0282     scale(hCalo);
0283     hCalo->Draw("EPSAME");
0284
0285     hCor->SetMarkerStyle(29);
0286     hCor->SetMarkerColor(3);
0287     hCor->SetLineColor(3);
0288     hCor->SetStats(kFALSE);
0289     scale(hCor);
0290     hCor->Draw("EPSAME");
0291
0292     if (!strcmp(id,"0")){
0293     TLegend *leg = new TLegend(0.65,0.75,0.9,0.9, "-0.5<#eta_{1}<0.5,  -0.5<#eta_{2}<0.5");
0294     leg->AddEntry(hGen,"Generated Jets","P");
0295     leg->AddEntry(hCalo,"Calorimetry Jets","P");
0296     leg->AddEntry(hCor,"Corrected Jets","P");}
0297
0298     else if (!strcmp(id,"1")){
0299     TLegend *leg = new TLegend(0.65,0.75,0.9,0.9, "-1.0<#eta_{1}<1.0,  -1.0<#eta_{2}<1.0");
0300     leg->AddEntry(hGen,"Generated Jets","P");
0301     leg->AddEntry(hCalo,"Calorimetry Jets","P");
0302     leg->AddEntry(hCor,"Corrected Jets","P");}
0303
0304     else if (!strcmp(id,"2")){
0305     TLegend *leg = new TLegend(0.65,0.75,0.9,0.9, "0.5<#eta_{1}<1.0,   0.5<#eta_{2}<1.0 ");
0306     leg->AddEntry(hGen,"Generated Jets","P");
0307     leg->AddEntry(hCalo,"Calorimetry Jets","P");
0308     leg->AddEntry(hCor,"Corrected Jets","P");}
0309
0310     else if (!strcmp(id,"3")){
0311     TLegend *leg = new TLegend(0.65,0.75,0.9,0.9, "-1.0<#eta_{1}<-0.5,  -1.0<#eta_{2}<-0.5");
0312     leg->AddEntry(hGen,"Generated Jets","P");
0313     leg->AddEntry(hCalo,"Calorimetry Jets","P");
0314     leg->AddEntry(hCor,"Corrected Jets","P");}
0315
0316     else if (!strcmp(id,"4")){
0317     TLegend *leg = new TLegend(0.65,0.75,0.9,0.9, "0.5<#eta_{1}<1.0,  -1.0<#eta_{2}<-0.5");
0318     leg->AddEntry(hGen,"Generated Jets","P");
0319     leg->AddEntry(hCalo,"Calorimetry Jets","P");
0320     leg->AddEntry(hCor,"Corrected Jets","P");}
0321
0322     else if (!strcmp(id,"5")){
0323     TLegend *leg = new TLegend(0.65,0.75,0.9,0.9, "-1.0<#eta_{1}<-0.5,  0.5<#eta_{2}<1.0");
0324     leg->AddEntry(hGen,"Generated Jets","P");
0325     leg->AddEntry(hCalo,"Calorimetry Jets","P");
0326     leg->AddEntry(hCor,"Corrected Jets","P");}
0327
0328     else  (!strcmp(id,"6")){
0329     TLegend *leg = new TLegend(0.65,0.75,0.9,0.9, "0.5 < |#eta| < 1.0");
0330     leg->AddEntry(hGen,"Generated Jets","P");
0331     leg->AddEntry(hCalo,"Calorimetry Jets","P");
0332     leg->AddEntry(hCor,"Corrected Jets","P");}
0333
0334     leg->Draw();
0335     c1->Print("c1.ps");
0336
0337
0338     if (!strcmp(id, "6")){
0339
0340         TH1F* h0Gen;    TH1F* h0GenVar;
0341         TH1F* h0Calo;   TH1F* h0CaloVar;
0342         TH1F* h0Cor;    TH1F* h0CorVar;
0343
0344         TString sGen = "hGen0";
0345         TString sCalo = "hCalo0";
0346         TString sCor = "hCor0";
0347
0348         //title
0349         TString ratioTitle = "Dijet Ratio";
0350
0351         Addition(sGen, h0Gen, h0GenVar);
0352         Addition(sCalo, h0Calo, h0CaloVar);
0353         Addition(sCor, h0Cor, h0CorVar);
0354
0355         //scaling
0356         scale(h0Gen);
0357         scale(h0Calo);
0358         scale(h0Cor);
0359
0360         //Ratio
0361         TCanvas *c1 = new TCanvas("c1", "c1",3,25,999,799);
0362         TH1F* ratioGen = Ratio(h0Gen,hGen);
0363         TH1F* ratioCalo = Ratio(h0Calo,hCalo);
0364         TH1F* ratioCor = Ratio(h0Cor,hCor);
0365
0366
0367
0368     cout << " Ratio = " << ratioCor->GetMean(2) << "+/-" << ratioCor->GetMeanError(2) << endl;
0369
0370     ratioGen->SetTitle("DiJet Ratio");
0371     ratioGen->GetXaxis()->SetTitle("DiJet Mass (in GeV)");
0372     ratioGen->GetXaxis()->SetRangeUser(330.,7600.);
0373     ratioGen->GetYaxis()->SetRangeUser(0.,1.2);
0374     ratioGen->GetYaxis()->SetTitle("Ratio=N(|#eta|<0.5)/N(0.5<|#eta|<1)");
0375     ratioGen->GetYaxis()->SetTitleOffset(1.3);
0376     ratioGen->SetMarkerStyle(20);
0377     ratioGen->SetMarkerColor(1);
0378     ratioGen->SetStats(kFALSE);
0379     ratioGen->Draw("E1P");
0380
0381
0382     ratioCalo->SetMarkerStyle(21);
0383     ratioCalo->SetMarkerColor(2);
0384     ratioCalo->SetStats(kFALSE);
0385     ratioCalo->Draw("E1PSAME");
0386
0387     ratioCor->SetMarkerStyle(29);
0388     ratioCor->SetMarkerColor(3);
0389     ratioCor->SetStats(kFALSE);
0390     ratioCor->Draw("E1PSAME");
0391
0392     TLegend *leg = new TLegend(0.65,0.75,0.9,0.9, "Ratio");
0393     leg->AddEntry(ratioGen,"Generated Jets","P");
0394     leg->AddEntry(ratioCalo,"Calorimetry Jets","P");
0395     leg->AddEntry(ratioCor,"Corrected Jets","P");
0396     leg->Draw();
0397
0398     c1->Print("c1.ps");
0399
0400     c1 = fitting(c1, ratioGen, "ratioGen", f);
0401     c1 = fitting(c1, ratioCalo,"ratioCalo", f);
0402     c1 = fitting(c1, ratioCor, "ratioCor", f);
0403
0404
0405
0406 }// id ==6
0407
0408     if (!strcmp(id, "0")){
0409
0410         TH1F* h1Gen;    TH1F* h1GenVar;
0411         TH1F* h1Calo;   TH1F* h1CaloVar;
0412         TH1F* h1Cor;    TH1F* h1CorVar;
0413
0414         TString sGen = "hGen1";
0415         TString sCalo = "hCalo1";
0416         TString sCor = "hCor1";
0417
0418         //title
0419         Addition(sGen, h1Gen, h1GenVar);
0420         Addition(sCalo, h1Calo, h1CaloVar);
0421         Addition(sCor, h1Cor, h1CorVar);
0422
0423         //scaling
0424         scale(h1Gen);
0425         scale(h1Calo);
0426         scale(h1Cor);
0427
0428         // ration CorJets/GenJets
0429         TCanvas *c1 = new TCanvas("c1", "c1",3,25,999,799);
0430         TH1F* R0_CorGen = Ratio(hCor,hGen);
0431         TH1F* R1_CorGen = Ratio(h1Cor,h1Gen);
0432
0433         R0_CorGen->SetTitle("DiJet Ratio");
0434         R0_CorGen->GetXaxis()->SetTitle("DiJet Mass (in GeV)");
0435         R0_CorGen->GetXaxis()->SetRangeUser(330.,7600.);
0436         R0_CorGen->GetYaxis()->SetRangeUser(0.6,1.6);
0437         R0_CorGen->GetYaxis()->SetTitle("Ratio (CorrectedJet/GenJet)");
0438         R0_CorGen->GetYaxis()->SetTitleOffset(1.3);
0439         R0_CorGen->SetMarkerStyle(20);
0440         R0_CorGen->SetMarkerColor(1);
0441         R0_CorGen->SetLineColor(1);
0442         R0_CorGen->SetStats(kFALSE);
0443         R0_CorGen->Draw("E1P");
0444
0445         R1_CorGen->SetMarkerStyle(21);
0446         R1_CorGen->SetMarkerColor(2);
0447         R1_CorGen->SetLineColor(2);
0448         R1_CorGen->SetStats(kFALSE);
0449         R1_CorGen->Draw("E1PSAME");
0450
0451         TLegend *leg = new TLegend(0.65,0.75,0.9,0.9, "Ratio: CorrectedJets/GenJets");
0452         leg->AddEntry(R0_CorGen,"|#eta| < 0.5","P");
0453         leg->AddEntry(R1_CorGen,"|#eta| < 1.0","P");
0454         leg->Draw();
0455         c1->Print("c1.ps");
0456             c1 = fitting(c1, R0_CorGen, "R0_CorGen", f);
0457             c1 = fitting(c1, R1_CorGen, "R1_CorGen", f);
0458     } //id == 0
0459
0460     return c1;
0461     delete hGen, hCalo, hCor;
0462     delete hGenVar, hCaloVar, hCorVar;
0463
0464 }//plot
0465
0466 int MassAnaBackgrd(){
0467
0468     TCanvas *c1 = new TCanvas("c1", "c1",3,25,999,799);
0469     c1->SetLogy(); c1->SetGridx(); c1->SetGridy();
0470     gSystem->Exec("rm -f backgrd.root");
0471     TFile* f = new TFile("backgrd.root", "new");
0472
0473     c1->Print("c1.ps(");
0474
0475     c1 = plot(c1, "0", f);
0476     c1 = plot(c1, "1", f);
0477     c1 = plot(c1, "2", f);
0478     c1 = plot(c1, "3", f);
0479     c1 = plot(c1, "4", f);
0480     c1 = plot(c1, "5", f);
0481     c1 = plot(c1, "6", f);
0482
0483     c1->Print("c1.ps)");
0484     f->Close();
0485
0486 //  delete c1;
0487     return 0;
0488 }
0489