test/Macros/CompareProbs.cc

0001 #include "TFile.h"
0002 #include "TCanvas.h"
0003 #include "TH2D.h"
0004 #include "TProfile2D.h"
0005 #include <iostream>
0006 #include <map>
0007 #include <sstream>
0008 #include <cmath>
0009
0010 /**
0011  * This macro compares the pdf used in the likelihood fit with the mass distributions. <br>
0012  * It can be used to check that the fit converged well in different bins of pt, eta, phi...
0013  */
0014
0015 void projectHisto(const TH2 * histo, std::map<unsigned int, TH1D*> & profileMap)
0016 {
0017   unsigned int binsX = histo->GetNbinsX();
0018   unsigned int binsY = histo->GetNbinsY();
0019   std::cout << "binsX = " << binsX << " binsY = " << binsY << std::endl;
0020   for( unsigned int i=1; i<binsX; ++i ) {
0021     std::stringstream ss;
0022     ss << i;
0023     TH1D * projected = histo->ProjectionY(TString(histo->GetName())+"_"+ss.str(), i, i);
0024     profileMap.insert(std::make_pair(i, projected));
0025     // if( i == 20 ) {
0026     //   projected->Draw();
0027     // }
0028   }
0029 }
0030
0031 void CompareProbs()
0032 {
0033   TFile * inputFile = new TFile("2_MuScleFit.root", "READ");
0034
0035   // Take the pdf histogram
0036   TProfile2D * profile = (TProfile2D*)inputFile->FindObjectAny("hMassProbVsMu_fine_MassVSEta");
0037   TH2D * probHisto = profile->ProjectionXY();
0038   probHisto->SetName(TString(profile->GetName())+"_TH2D");
0039
0040   std::map<unsigned int, TH1D*> projectProb;
0041   projectHisto(probHisto, projectProb);
0042
0043   // Take the mass histogram
0044   TH2F * massHisto = (TH2F*)inputFile->FindObjectAny("hRecBestResVSMu_MassVSEta");
0045   std::map<unsigned int, TH1D*> projectMass;
0046   projectHisto(massHisto, projectMass);
0047
0048   // Draw the histograms
0049   if( projectProb.size() != projectMass.size() ) {
0050     std::cout << "Error: bins for pdf = " << projectProb.size() << " != bins for mass = " << projectMass.size() << std::endl;
0051     return;
0052   }
0053
0054   std::map<unsigned int, TH1D*>::const_iterator pIt = projectProb.begin();
0055   unsigned int totPads = 0;
0056   for( ; pIt != projectProb.end(); ++pIt ) {
0057     if( pIt->second->GetEntries() != 0 ) ++totPads;
0058   }
0059   TCanvas * newCanvas = new TCanvas("newCanvas", "newCanvas", 1000, 800);
0060   std::cout << "totPads = " << totPads << ", totPads%2 = " << totPads%2 << std::endl;
0061   unsigned int psq(sqrt(totPads));
0062   if( psq*psq == totPads ) newCanvas->Divide(psq, psq);
0063   else if( (psq+1)*psq >= totPads ) newCanvas->Divide(psq+1, psq);
0064   else newCanvas->Divide(psq+1, psq+1);
0065   newCanvas->Draw();
0066   pIt = projectProb.begin();
0067   unsigned int pad = 1;
0068   for( ; pIt != projectProb.end(); ++pIt ) {
0069     if( pIt->second->GetEntries() == 0 ) continue;
0070     // if( pIt->first != 20 ) continue;
0071     newCanvas->cd(pad);
0072
0073     // double xMin = 3.0969-0.2;
0074     // double xMax = 3.0969+0.2;
0075     double xMin = 2.85;
0076     double xMax = 3.3;
0077
0078     TH1D * mP = pIt->second;
0079     mP->Scale(1/mP->Integral(mP->FindBin(xMin), mP->FindBin(xMax)));
0080     mP->SetLineColor(2);
0081     mP->SetMarkerColor(2);
0082     mP->SetMarkerStyle(2);
0083     mP->SetMarkerSize(0.2);
0084
0085     TH1D * mM = projectMass[pIt->first];
0086     mM->Scale(1/mM->Integral(mM->FindBin(xMin), mM->FindBin(xMax)));
0087
0088     mP->Draw();
0089     std::cout << "pad = " << pad << std::endl;
0090     mP->GetXaxis()->SetRangeUser(xMin, xMax);
0091     mM->Draw("same");
0092     ++pad;
0093   }
0094 }