Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:57

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