PWGJE/macros/PID/extractMCmuonToElRatio.C

   1 #include "TCanvas.h"
   2 #include "TF1.h"
   3 #include "TFile.h"
   4 #include "TH1D.h"
   5 #include "THnSparse.h"
   6 #include "TObjArray.h"
   7
   8 #include <iostream>
   9
  10 #include "THnSparseDefinitions.h"
  11
  12 Int_t extractMCmuonToElRatio(const TString pathNameData, const TString listName, const Double_t lowerJetPt, const Double_t upperJetPt,
  13                              const Double_t lowerCentrality, const Double_t upperCentrality, const Bool_t jetHandling)
  14 {
  15   TFile* fileData = TFile::Open(pathNameData.Data());
  16   if (!fileData) {
  17     printf("Failed to open data file \"%s\"\n", pathNameData.Data());
  18     return -1;
  19   }
  20
  21   TObjArray* histList = (TObjArray*)(fileData->Get(listName.Data()));
  22   if (!histList) {
  23     std::cout << "Failed to load list \"" << listName.Data() << "\"!" << std::endl;
  24     return -1;
  25   }
  26
  27   // Extract the data histogram
  28   THnSparse* hPIDdata = dynamic_cast<THnSparse*>(histList->FindObject("hPIDdataAll"));
  29   if (!hPIDdata) {
  30     std::cout << "Failed to load data histo!" << std::endl;
  31     return -1;
  32   }
  33
  34   // Set proper errors, if not yet calculated
  35   if (!hPIDdata->GetCalculateErrors()) {
  36     std::cout << "Re-calculating errors of " << hPIDdata->GetName() << "..." << std::endl;
  37     hPIDdata->Sumw2();
  38     Long64_t nBinsTHnSparse = hPIDdata->GetNbins();
  39     Double_t binContent = 0;
  40
  41     for (Long64_t bin = 0; bin < nBinsTHnSparse; bin++) {
  42       binContent = hPIDdata->GetBinContent(bin);
  43       hPIDdata->SetBinError(bin, TMath::Sqrt(binContent));
  44     }
  45   }
  46
  47
  48   Int_t lowerJetPtBinLimit = -1;
  49   Int_t upperJetPtBinLimit = -1;
  50   Bool_t restrictJetPtAxis = kFALSE;
  51   Double_t actualLowerJetPt = -1.;
  52   Double_t actualUpperJetPt = -1.;
  53
  54   if (lowerJetPt >= 0 && upperJetPt >= 0) {
  55     // Add subtract a very small number to avoid problems with values right on the border between to bins
  56     lowerJetPtBinLimit = hPIDdata->GetAxis(kPidJetPt)->FindBin(lowerJetPt + 0.001);
  57     upperJetPtBinLimit = hPIDdata->GetAxis(kPidJetPt)->FindBin(upperJetPt - 0.001);
  58
  59     // Check if the values look reasonable
  60     if (lowerJetPtBinLimit <= upperJetPtBinLimit && lowerJetPtBinLimit >= 1 &&
  61         upperJetPtBinLimit <= hPIDdata->GetAxis(kPidJetPt)->GetNbins()) {
  62       actualLowerJetPt = hPIDdata->GetAxis(kPidJetPt)->GetBinLowEdge(lowerJetPtBinLimit);
  63       actualUpperJetPt = hPIDdata->GetAxis(kPidJetPt)->GetBinUpEdge(upperJetPtBinLimit);
  64
  65       restrictJetPtAxis = kTRUE;
  66     }
  67     else {
  68       std::cout << std::endl;
  69       std::cout << "Requested jet pT range out of limits or upper and lower limit are switched!" << std::endl;
  70       return -1;
  71     }
  72   }
  73
  74   std::cout << "jet pT: ";
  75   if (restrictJetPtAxis) {
  76     std::cout << actualLowerJetPt << " - " << actualUpperJetPt << std::endl;
  77     hPIDdata->GetAxis(kPidJetPt)->SetRangeUser(actualLowerJetPt, actualUpperJetPt);
  78   }
  79   else {
  80     std::cout << "All" << std::endl;
  81   }
  82
  83
  84   // If desired, restrict centrality axis
  85   Int_t lowerCentralityBinLimit = -1;
  86   Int_t upperCentralityBinLimit = -1;
  87   Bool_t restrictCentralityAxis = kFALSE;
  88   Double_t actualLowerCentrality = -1.;
  89   Double_t actualUpperCentrality = -1.;
  90
  91   if (lowerCentrality >= -1 && upperCentrality >= -1) {
  92     // Add subtract a very small number to avoid problems with values right on the border between to bins
  93     lowerCentralityBinLimit = hPIDdata->GetAxis(kPidCentrality)->FindBin(lowerCentrality + 0.001);
  94     upperCentralityBinLimit = hPIDdata->GetAxis(kPidCentrality)->FindBin(upperCentrality - 0.001);
  95
  96     // Check if the values look reasonable
  97     if (lowerCentralityBinLimit <= upperCentralityBinLimit && lowerCentralityBinLimit >= 1
  98         && upperCentralityBinLimit <= hPIDdata->GetAxis(kPidCentrality)->GetNbins()) {
  99       actualLowerCentrality = hPIDdata->GetAxis(kPidCentrality)->GetBinLowEdge(lowerCentralityBinLimit);
 100       actualUpperCentrality = hPIDdata->GetAxis(kPidCentrality)->GetBinUpEdge(upperCentralityBinLimit);
 101
 102       restrictCentralityAxis = kTRUE;
 103     }
 104     else {
 105       std::cout << std::endl;
 106       std::cout << "Requested centrality range out of limits or upper and lower limit are switched!" << std::endl;
 107       return -1;
 108     }
 109   }
 110
 111   std::cout << "centrality: ";
 112   if (restrictCentralityAxis) {
 113     std::cout << actualLowerCentrality << " - " << actualUpperCentrality << std::endl;
 114      hPIDdata->GetAxis(kPidCentrality)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
 115   }
 116   else {
 117     std::cout << "All" << std::endl;
 118   }
 119
 120   // Obtain MC information about particle yields
 121   hPIDdata->GetAxis(kSelectSpecies)->SetRange(1, 1); // Do not count each particle more than once
 122
 123   hPIDdata->GetAxis(kPidMCpid)->SetRange(3, 3);
 124   TH1D* hMCmuonYield = (TH1D*)hPIDdata->Projection(kPidPt, "e");
 125   hMCmuonYield->SetName("hMCmuonYield");
 126
 127   TCanvas* cYields = new TCanvas;
 128   cYields->SetLogx(kTRUE);
 129   cYields->SetLogy(kTRUE);
 130   hMCmuonYield->SetLineColor(kOrange);
 131   hMCmuonYield->GetXaxis()->SetRangeUser(0.15, 50.);
 132   hMCmuonYield->Draw();
 133
 134   hPIDdata->GetAxis(kPidMCpid)->SetRange(1, 1);
 135   TH1D* hMCelectronYield = (TH1D*)hPIDdata->Projection(kPidPt, "e");
 136   hMCelectronYield->SetName("hMCelectronYield");
 137
 138   hMCelectronYield->SetLineColor(kMagenta);
 139   hMCelectronYield->Draw("same");
 140
 141   TH1D* hMCmuonToElRatio = new TH1D(*hMCmuonYield);
 142   hMCmuonToElRatio->Divide(hMCelectronYield);
 143
 144   TCanvas* cFit = new TCanvas;
 145   cFit->SetLogx(kTRUE);
 146   hMCmuonToElRatio->GetYaxis()->SetRangeUser(0., 1.);
 147   hMCmuonToElRatio->Draw();
 148
 149   TF1 f("f", "[0]+[1]/TMath::Min(x, [4])+[2]*TMath::Min(x, [4])+[3]*TMath::Min(x, [4])*TMath::Min(x, [4])+[5]*TMath::Min(x, [4])*TMath::Min(x, [4])*TMath::Min(x, [4])+[6]*(x>[7])*TMath::Min(x-[7], [8]-[7])",
 150         0.01, 50.);
 151
 152   f.SetParameters(-0.688, 0.042, 4.52, -6.2, 0.53, 2.89, 0.035, 2.3, 6.);
 153   restrictJetPtAxis = kTRUE;
 154   if (jetHandling)
 155     f.FixParameter(8, 6.);
 156
 157   f.SetParLimits(4, 0.2, 2.0);
 158   f.SetParLimits(7, 1.0, 6.0);
 159
 160   hMCmuonToElRatio->Fit(&f, jetHandling ? "+" : "+W", "", 0.15, restrictJetPtAxis ? actualUpperJetPt : 15.0);
 161
 162   return 0;
 163 }