test/Macros/Plot_mass.C

0001 #include "TH1F.h"
0002 #include "TProfile.h"
0003 #include "TProfile.h"
0004 #include "TCanvas.h"
0005 #include "TFile.h"
0006 #include "TStyle.h"
0007 #include "TROOT.h"
0008
0009 #include <iostream>
0010
0011
0012 /// Helper class containing the histograms
0013 class histos
0014 {
0015 public:
0016   histos(TFile * inputFile)
0017   {
0018     mass         = dynamic_cast<TH1F*>(inputFile->Get("hRecBestResAllEvents_Mass"));
0019     TDirectory * MassPdir = dynamic_cast<TDirectory*>(inputFile->Get("Mass_P"));
0020     massProb     = dynamic_cast<TProfile*>(MassPdir->Get("Mass_PProf"));
0021     // TDirectory * MassPdir = dynamic_cast<TDirectory*>(inputFile->Get("Mass_Probability"));
0022     // massProb     = dynamic_cast<TProfile*>(MassPdir->Get("Mass_Probability"));
0023     TDirectory * massFinePdir = dynamic_cast<TDirectory*>(inputFile->Get("Mass_fine_P"));
0024     massFineProb = dynamic_cast<TProfile*>(massFinePdir->Get("Mass_fine_PProf"));
0025     // TDirectory * massFinePdir = dynamic_cast<TDirectory*>(inputFile->Get("Mass_fine_Probability"));
0026     // massFineProb = dynamic_cast<TProfile*>(massFinePdir->Get("Mass_fine_Probability"));
0027     likePt       = dynamic_cast<TProfile*>(inputFile->Get("hLikeVSMu_LikelihoodVSPt_prof"));
0028     likePhi      = dynamic_cast<TProfile*>(inputFile->Get("hLikeVSMu_LikelihoodVSPhi_prof"));
0029     likeEta      = dynamic_cast<TProfile*>(inputFile->Get("hLikeVSMu_LikelihoodVSEta_prof"));
0030   }
0031   TH1F * mass;
0032   TProfile * massProb;
0033   TProfile * massFineProb;
0034   TProfile * likePt;
0035   TProfile * likePhi;
0036   TProfile * likeEta;
0037 };
0038
0039 /// Helper function to compute the bins corresponding to an interval in the x axis
0040 int getXbins(const TH1 * h, const double & xMin, const double & xMax)
0041 {
0042   // To get the correct integral for rescaling determine the borders
0043   const TAxis * xAxis = h->GetXaxis();
0044   double xAxisMin = xAxis->GetXmin();
0045   // The bins have all the same width, therefore we can use this
0046   double binWidth = xAxis->GetBinWidth(1);
0047   int xMinBin = int((xMin - xAxisMin)/binWidth) + 1;
0048   int xMaxBin = int((xMax - xAxisMin)/binWidth) + 1;
0049   // std::cout << "xMinBin = " << xMinBin << std::endl;
0050   // std::cout << "xMaxBin = " << xMaxBin << std::endl;
0051   // std::cout << "binWidth = " << binWidth << std::endl;
0052   return( xMaxBin - xMinBin );
0053 }
0054
0055 /// Helper function to draw mass and mass probability histograms
0056 void drawMasses(const double ResMass, const double ResHalfWidth, histos & h, const int ires, const int rebin = 1)
0057 {
0058   TH1F * mass = (TH1F*)h.mass->Clone();
0059   TProfile * massProb = 0;
0060   mass->Rebin(rebin);
0061   // Use massProb for the Z and fineMass for the other resonances
0062   if( ires == 0 ) massProb = (TProfile*)h.massProb->Clone();
0063   else massProb = (TProfile*)h.massFineProb->Clone();
0064   mass->SetAxisRange(ResMass - ResHalfWidth, ResMass + ResHalfWidth);
0065   mass->SetLineColor(kRed);
0066
0067   // Set a consistent binning
0068   int massXbins = getXbins( mass, (ResMass - ResHalfWidth), (ResMass + ResHalfWidth) );
0069   int massProbXbins = getXbins( massProb, (ResMass - ResHalfWidth), (ResMass + ResHalfWidth) );
0070
0071   if( massProbXbins > massXbins && massXbins != 0 ) {
0072     std::cout << "massProbXbins("<<massProbXbins<<") > " << "massXbins("<<massXbins<<")" << std::endl;
0073     std::cout << "massProb = " << massProb << std::endl;
0074     std::cout << "mass = " << mass << std::endl;
0075     massProb->Rebin(massProbXbins/massXbins);
0076   }
0077   else if( massXbins > massProbXbins && massProbXbins != 0 ) {
0078     std::cout << "massXbins("<<massXbins<<") > " << "massProbXbins("<<massProbXbins<<")" << std::endl;
0079     std::cout << "massProb = " << massProb << std::endl;
0080     std::cout << "mass = " << mass << std::endl;
0081     mass->Rebin(massXbins/massProbXbins);
0082   }
0083
0084   mass->SetAxisRange(ResMass - ResHalfWidth, ResMass + ResHalfWidth);
0085   massProb->SetAxisRange(ResMass - ResHalfWidth, ResMass + ResHalfWidth);
0086
0087   double massProbIntegral = massProb->Integral("width");
0088   // double massProbIntegral = massProb->Integral();
0089   double normFactor = 1.;
0090   if( massProbIntegral != 0 ) normFactor = (mass->Integral("width"))/massProbIntegral;
0091   // if( massProbIntegral != 0 ) normFactor = (mass->Integral())/massProbIntegral;
0092
0093   massProb->SetLineColor(kBlue);
0094   massProb->Scale(normFactor);
0095   if( ires == 3 ) {
0096     std::cout << "massProbIntegralWidth = " << massProb->Integral("width") << std::endl;
0097     std::cout << "massIntegralWidth = " << mass->Integral("width") << std::endl;
0098     std::cout << "massProbIntegral = " << massProb->Integral() << std::endl;
0099     std::cout << "massIntegral = " << mass->Integral() << std::endl;
0100   }
0101
0102   mass->SetMarkerColor(kRed);
0103   // ATTENTION: put so that the maximum is correct
0104   if( (massProb->GetMaximum()) > mass->GetMaximum() ) {
0105     massProb->Draw("PE");
0106     mass->Draw("PESAMEHISTO");
0107   }
0108   else {
0109     mass->Draw("PE");
0110     massProb->Draw("PESAMEHISTO");
0111   }
0112 }
0113
0114 /// Helper function to draw likelihood histograms
0115 void drawLike(TProfile * likeHisto1, TProfile * likeHisto2)
0116 {
0117   likeHisto1->SetMinimum(-8);
0118   likeHisto1->SetMaximum(-4);
0119   likeHisto1->Draw("");
0120   likeHisto2->SetLineColor(kRed);
0121   likeHisto2->Draw("SAME");
0122 }
0123
0124 /// Function drawing the histograms
0125 void Plot_mass(const TString & fileNameBefore = "0", const TString & fileNameAfter = "1", const int rebin = 1) {
0126
0127   gROOT->SetBatch(true);
0128   gStyle->SetOptStat ("111111");
0129   gStyle->SetOptFit (1);
0130
0131   double ResHalfWidth[6] = {20., 0.5, 0.5, 0.5, 0.2, 0.2};
0132   // double ResHalfWidth[6] = {20., 0.25, 0.25, 0.25, 0.2, 0.2};
0133   // double ResHalfWidth[6] = {20., 0.35, 0.35, 0.35, 0.2, 0.2};
0134   double ResMass[6] = {90.986, 10.3552, 10.0233, 9.4603, 3.68609, 3.0969};
0135
0136   TFile * inputFile1 = new TFile(fileNameBefore+"_MuScleFit.root", "READ");
0137   TFile * inputFile2 = new TFile(fileNameAfter+"_MuScleFit.root", "READ");
0138   histos histos1(inputFile1);
0139   histos histos2(inputFile2);
0140
0141   TFile * outputFile = new TFile("plotMassOutput.root", "RECREATE");
0142
0143   TCanvas * Allres = new TCanvas ("Allres", "All resonances", 600, 600);
0144   Allres->Divide (2,3);
0145   for (int ires=0; ires<6; ires++) {
0146     Allres->cd(ires+1);
0147     drawMasses(ResMass[ires], ResHalfWidth[ires], histos1, ires, rebin);
0148   }
0149   Allres->Write();
0150   TCanvas * AllresTogether = new TCanvas ("AllresTogether", "All resonances together", 600, 600);
0151   AllresTogether->Divide(2,2);
0152   // The J/Psi + Psi2S
0153   AllresTogether->cd(1);
0154   drawMasses(ResMass[5], 2, histos1, 3, rebin);
0155  // Draw also the Upsilons
0156   AllresTogether->cd(2);
0157   drawMasses(ResMass[2], 2, histos1, 3, rebin);
0158   // All low Pt resonances
0159   AllresTogether->cd(3);
0160   drawMasses(6., 6, histos1, 3, rebin);
0161   // All resonances
0162   AllresTogether->cd(4);
0163   drawMasses(50., 50, histos1, 0, rebin);
0164   AllresTogether->Write();
0165
0166   TCanvas * Allres2 = new TCanvas ("Allres2", "All resonances", 600, 600);
0167   Allres2->Divide (2,3);
0168   for (int ires=0; ires<6; ires++) {
0169     Allres2->cd(ires+1);
0170     drawMasses(ResMass[ires], ResHalfWidth[ires], histos2, ires, rebin);
0171   }
0172   Allres2->Write();
0173   TCanvas * AllresTogether2 = new TCanvas ("AllresTogether2", "All resonances together", 600, 600);
0174   AllresTogether2->Divide(2,2);
0175   // The J/Psi + Psi2S
0176   AllresTogether2->cd(1);
0177   drawMasses(ResMass[5], 2, histos2, 3, rebin);
0178   // Draw also the Upsilons
0179   AllresTogether2->cd(2);
0180   drawMasses(ResMass[2], 2, histos2, 3, rebin);
0181   // All low Pt resonances
0182   AllresTogether2->cd(3);
0183   drawMasses(6., 6, histos2, 3, rebin);
0184   // All resonances
0185   AllresTogether2->cd(4);
0186   drawMasses(50., 50, histos2, 0, rebin);
0187   AllresTogether2->Write();
0188
0189   TCanvas * LR = new TCanvas ("LR", "Likelihood and resolution before and after corrections", 600, 600 );
0190   LR->Divide (2,3);
0191
0192   LR->cd(1);
0193   drawLike(histos1.likePt, histos2.likePt);
0194   LR->cd(3);
0195   drawLike(histos1.likeEta, histos2.likeEta);
0196   LR->cd(5);
0197   drawLike(histos1.likePhi, histos2.likePhi);
0198
0199   LR->Write();
0200
0201   TCanvas * G = new TCanvas ("G", "Mass before and after corrections", 600, 600 );
0202   G->Divide (3,3);
0203
0204   TH1F * MuZ = histos1.mass;
0205   TH1F * McZ = histos2.mass;
0206
0207   G->cd(1);
0208   MuZ->SetAxisRange (2.8, 3.4);
0209   MuZ->SetLineColor (kBlue);
0210   McZ->SetLineColor (kRed);
0211   MuZ->DrawCopy ("HISTO");
0212   McZ->DrawCopy ("SAME");
0213
0214   G->cd(2);
0215   MuZ->SetAxisRange (9., 10.);
0216   McZ->SetAxisRange (9., 10.);
0217   MuZ->DrawCopy ("HISTO");
0218   McZ->DrawCopy ("SAME");
0219
0220   G->cd(3);
0221   MuZ->SetAxisRange (70., 110.);
0222   McZ->SetAxisRange (70., 110.);
0223   MuZ->DrawCopy ("HISTO");
0224   McZ->DrawCopy ("SAME");
0225
0226   G->cd(4);
0227   MuZ->SetAxisRange (2.8, 3.4);
0228   MuZ->Fit ("gaus", "", "", 3., 3.2);
0229   MuZ->DrawCopy();
0230   G->cd(7);
0231   McZ->SetAxisRange (2.8, 3.4);
0232   McZ->Fit ("gaus", "", "", 3., 3.2);
0233   McZ->DrawCopy();
0234   G->cd(5);
0235   MuZ->SetAxisRange (9., 10.);
0236   MuZ->Fit ("gaus", "", "", 9., 10.);
0237   MuZ->DrawCopy();
0238   G->cd(8);
0239   McZ->SetAxisRange (9., 10.);
0240   McZ->Fit ("gaus", "", "", 9., 10.);
0241   McZ->DrawCopy();
0242   G->cd(6);
0243   MuZ->SetAxisRange (70., 110.);
0244   MuZ->Fit ("gaus", "", "", 70., 110.);
0245   MuZ->DrawCopy();
0246   G->cd(9);
0247   McZ->SetAxisRange (70., 110.);
0248   McZ->Fit ("gaus", "", "", 70., 110.);
0249   McZ->DrawCopy();
0250
0251   G->Write();
0252
0253   TCanvas * Spectrum = new TCanvas ("Spectrum", "Mass before and after corrections", 600, 600 );
0254   Spectrum->cd(1);
0255   MuZ->SetAxisRange(1.,15.);
0256   McZ->SetAxisRange(1.,15.);
0257   McZ->SetLineColor(kRed);
0258   McZ->DrawCopy("HISTO");
0259   Spectrum->cd(2);
0260   MuZ->SetLineColor(kBlack);
0261   MuZ->DrawCopy("SAMEHISTO");
0262
0263   Spectrum->Write();
0264
0265   outputFile->Close();
0266 }
0267