JetAnalyzers/test/CMSDASmacro.c

0001 #include "TFile.h"
0002 #include "TH2F.h"
0003 #include "TCanvas.h"
0004 #include "TROOT.h"
0005 #include "TF1.h"
0006 #include "TStyle.h"
0007 #include "TLine.h"
0008 #include <iostream>
0009
0010 using std::cout;
0011
0012 void CMSDASmacro() {
0013   TFile* rootfile = new TFile("histosData.root");
0014   rootfile->cd("dijetAna");
0015   gDirectory->ls();
0016
0017   //------------------Dijet mass spectrum---------------------------
0018
0019   gStyle->SetOptStat(0);
0020   gStyle->SetOptFit(1111);
0021   //Get pointers to some histograms we will want
0022   TH1D* h_CorDijetMass=(TH1D*)gROOT->FindObject("hCorDijetXsec");
0023
0024   //Divide by approx. luminosity to get cross section
0025   h_CorDijetMass->Scale(0.333333);
0026
0027   //Error bars, please.
0028   h_CorDijetMass->Sumw2();
0029
0030   //Make a home for the pretty plots
0031   TCanvas *cDijetMass=new TCanvas();
0032   // Set the current TPad to "Logy()". Not an English word.
0033   gPad->SetLogy();
0034   h_CorDijetMass->Draw();
0035
0036   // Declare the function to fit, over some range, and fit it.
0037   TF1 *massfit= new TF1("Dijet Mass Spectrum", "[0]*pow(1-x/7000.0,[1])/pow(x/7000,[2]+[3]*log(x/7000.))",489,2132);
0038   massfit->SetParameter(1,5.077);
0039   massfit->SetParameter(2,6.994);
0040   massfit->SetParameter(3,0.2658);
0041
0042   h_CorDijetMass->Fit(massfit, "R");
0043   h_CorDijetMass->SetMinimum(1E-3);
0044
0045   //Make a histogram for the values of the fit, with the binning of h_CorDijetMass
0046   TH1D* h_DataMinusFit = new TH1D(*h_CorDijetMass);
0047   h_DataMinusFit->SetTitle("(Data- Fit)/Fit;dijet mass (GeV)");
0048
0049   //Fill the histogram of the data minus the fit's values
0050   for (int bin=1; bin<=h_CorDijetMass->GetNbinsX(); bin++){
0051     double data_val = h_CorDijetMass->GetBinContent(bin);
0052     double fit_val = massfit->Eval(h_CorDijetMass->GetBinCenter(bin));
0053     double err_val  = h_CorDijetMass->GetBinError(bin);
0054     // Skip bins with no data value
0055     if (data_val != 0.0) {
0056       h_DataMinusFit->SetBinContent(bin, (data_val - fit_val)/fit_val );
0057       h_DataMinusFit->SetBinError(bin, err_val /fit_val );
0058     }
0059   }
0060
0061   //Move to the lower TPad and display the result
0062   TCanvas *cDijetMassResiduals=new TCanvas();
0063   h_DataMinusFit->SetMinimum(-1.);
0064   h_DataMinusFit->SetMaximum( 4);
0065   h_DataMinusFit->Draw();
0066   TLine * line = new TLine(489,0.,2132,0);
0067   line->SetLineStyle(2);
0068   line->Draw("same");
0069
0070   TCanvas *cDijetMassPulls=new TCanvas();
0071   TH1D* h_DijetMassPulls = new TH1D(*h_CorDijetMass);
0072   h_DijetMassPulls->SetTitle("(Data- Fit)/Error;dijet mass (GeV)");
0073
0074   //Fill the histogram of the data minus the fit's values
0075   for (int bin=1; bin<=h_CorDijetMass->GetNbinsX(); bin++){
0076     double data_val = h_CorDijetMass->GetBinContent(bin);
0077     double fit_val = massfit->Eval(h_CorDijetMass->GetBinCenter(bin));
0078     double err_val  = h_CorDijetMass->GetBinError(bin);
0079     std::cout<<h_CorDijetMass->GetBinCenter(bin)<<": "<<data_val<<"\n";
0080     // Skip bins with no data value
0081     if (data_val != 0.0) {
0082       h_DijetMassPulls->SetBinContent(bin, (data_val - fit_val)/err_val );
0083       h_DijetMassPulls->SetBinError(bin, err_val /err_val );
0084     }
0085   }
0086   h_DijetMassPulls->Draw();
0087   h_DijetMassPulls->GetYaxis()->SetRangeUser(-2.5,2.5);
0088   line->Draw("same");
0089
0090   //------------------Dijet Centrality Ratio---------------------------
0091
0092
0093   //Get pointers to some histograms we will want
0094   TH1D* h_InnerDijetMass=(TH1D*)gROOT->FindObject("hInnerDijetMass");
0095   TH1D* h_OuterDijetMass=(TH1D*)gROOT->FindObject("hOuterDijetMass");
0096   //Error bars, please.
0097   h_InnerDijetMass->Sumw2();
0098   h_OuterDijetMass->Sumw2();
0099
0100
0101   //Make a home for the pretty plots
0102   TCanvas *cInnerOuterCounts=new TCanvas("cInnerOuterCounts","cInnerOuterCounts");
0103   // Set the current TPad to "Logy()". Not an English word.
0104   gPad->SetLogy();
0105   h_InnerDijetMass->SetLineColor(2);
0106   h_InnerDijetMass->SetMarkerStyle(25);
0107   h_InnerDijetMass->SetMarkerSize(0.7);
0108   h_InnerDijetMass->SetMarkerColor(2);
0109   h_OuterDijetMass->SetMarkerStyle(20);
0110   h_OuterDijetMass->SetMarkerSize(0.3);
0111   h_InnerDijetMass->Draw();
0112   h_OuterDijetMass->Draw("same");
0113
0114
0115   TCanvas *cDijetDeltaEtaRatio=new TCanvas("cDijetDeltaEtaRatio","cDijetDeltaEtaRatio");
0116   // Make the dijet delta eta ratio in a histogram
0117   TH1D* h_DijetDeltaEtaRatio = (TH1D*)h_InnerDijetMass->Clone();
0118   h_DijetDeltaEtaRatio->Divide(h_OuterDijetMass);
0119   h_DijetDeltaEtaRatio->SetMarkerStyle(20);
0120   h_DijetDeltaEtaRatio->SetTitle("Dijet |#Delta#eta| Ratio; dijet mass (GeV)");
0121   h_DijetDeltaEtaRatio->GetYaxis()->SetRangeUser(0,2.5);
0122   h_DijetDeltaEtaRatio->Draw();
0123
0124   TF1* flatline = new TF1("flatline", "[0]",489,2132);
0125   h_DijetDeltaEtaRatio->Fit(flatline, "R");
0126   flatline->Draw("same");
0127
0128   return;
0129 }