ROOT logo
#include "TCanvas.h"
#include "TF1.h"
#include "TFile.h"
#include "TH1D.h"
#include "THnSparse.h"
#include "TObjArray.h"

#include <iostream>

#include "THnSparseDefinitions.h"

Int_t extractMCmuonToElRatio(const TString pathNameData, const TString listName, const Double_t lowerJetPt, const Double_t upperJetPt,
                             const Double_t lowerCentrality, const Double_t upperCentrality, const Bool_t jetHandling)
{
  TFile* fileData = TFile::Open(pathNameData.Data());
  if (!fileData) {
    printf("Failed to open data file \"%s\"\n", pathNameData.Data());
    return -1;
  }
  
  TObjArray* histList = (TObjArray*)(fileData->Get(listName.Data()));
  if (!histList) {
    std::cout << "Failed to load list \"" << listName.Data() << "\"!" << std::endl;
    return -1;
  }
  
  // Extract the data histogram
  THnSparse* hPIDdata = dynamic_cast<THnSparse*>(histList->FindObject("hPIDdataAll"));
  if (!hPIDdata) {
    std::cout << "Failed to load data histo!" << std::endl;
    return -1;
  }
  
  // Set proper errors, if not yet calculated
  if (!hPIDdata->GetCalculateErrors()) {
    std::cout << "Re-calculating errors of " << hPIDdata->GetName() << "..." << std::endl;
    hPIDdata->Sumw2();
    Long64_t nBinsTHnSparse = hPIDdata->GetNbins();
    Double_t binContent = 0;
    
    for (Long64_t bin = 0; bin < nBinsTHnSparse; bin++) {
      binContent = hPIDdata->GetBinContent(bin);
      hPIDdata->SetBinError(bin, TMath::Sqrt(binContent));
    }
  }
  
  
  Int_t lowerJetPtBinLimit = -1;
  Int_t upperJetPtBinLimit = -1;
  Bool_t restrictJetPtAxis = kFALSE;
  Double_t actualLowerJetPt = -1.;
  Double_t actualUpperJetPt = -1.;
  
  if (lowerJetPt >= 0 && upperJetPt >= 0) {
    // Add subtract a very small number to avoid problems with values right on the border between to bins
    lowerJetPtBinLimit = hPIDdata->GetAxis(kPidJetPt)->FindBin(lowerJetPt + 0.001);
    upperJetPtBinLimit = hPIDdata->GetAxis(kPidJetPt)->FindBin(upperJetPt - 0.001);
    
    // Check if the values look reasonable
    if (lowerJetPtBinLimit <= upperJetPtBinLimit && lowerJetPtBinLimit >= 1 &&
        upperJetPtBinLimit <= hPIDdata->GetAxis(kPidJetPt)->GetNbins()) {
      actualLowerJetPt = hPIDdata->GetAxis(kPidJetPt)->GetBinLowEdge(lowerJetPtBinLimit);
      actualUpperJetPt = hPIDdata->GetAxis(kPidJetPt)->GetBinUpEdge(upperJetPtBinLimit);

      restrictJetPtAxis = kTRUE;
    }
    else {
      std::cout << std::endl;
      std::cout << "Requested jet pT range out of limits or upper and lower limit are switched!" << std::endl;
      return -1;
    }
  }
  
  std::cout << "jet pT: ";
  if (restrictJetPtAxis) {
    std::cout << actualLowerJetPt << " - " << actualUpperJetPt << std::endl;
    hPIDdata->GetAxis(kPidJetPt)->SetRangeUser(actualLowerJetPt, actualUpperJetPt); 
  }
  else {
    std::cout << "All" << std::endl;
  }
  
  
  // If desired, restrict centrality axis
  Int_t lowerCentralityBinLimit = -1;
  Int_t upperCentralityBinLimit = -1;
  Bool_t restrictCentralityAxis = kFALSE;
  Double_t actualLowerCentrality = -1.;
  Double_t actualUpperCentrality = -1.;
  
  if (lowerCentrality >= -1 && upperCentrality >= -1) {
    // Add subtract a very small number to avoid problems with values right on the border between to bins
    lowerCentralityBinLimit = hPIDdata->GetAxis(kPidCentrality)->FindBin(lowerCentrality + 0.001);
    upperCentralityBinLimit = hPIDdata->GetAxis(kPidCentrality)->FindBin(upperCentrality - 0.001);
    
    // Check if the values look reasonable
    if (lowerCentralityBinLimit <= upperCentralityBinLimit && lowerCentralityBinLimit >= 1
        && upperCentralityBinLimit <= hPIDdata->GetAxis(kPidCentrality)->GetNbins()) {
      actualLowerCentrality = hPIDdata->GetAxis(kPidCentrality)->GetBinLowEdge(lowerCentralityBinLimit);
      actualUpperCentrality = hPIDdata->GetAxis(kPidCentrality)->GetBinUpEdge(upperCentralityBinLimit);

      restrictCentralityAxis = kTRUE;
    }
    else {
      std::cout << std::endl;
      std::cout << "Requested centrality range out of limits or upper and lower limit are switched!" << std::endl;
      return -1;
    }
  }
  
  std::cout << "centrality: ";
  if (restrictCentralityAxis) {
    std::cout << actualLowerCentrality << " - " << actualUpperCentrality << std::endl;
     hPIDdata->GetAxis(kPidCentrality)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
  }
  else {
    std::cout << "All" << std::endl;
  }
  
  // Obtain MC information about particle yields
  hPIDdata->GetAxis(kSelectSpecies)->SetRange(1, 1); // Do not count each particle more than once
  
  hPIDdata->GetAxis(kPidMCpid)->SetRange(3, 3);
  TH1D* hMCmuonYield = (TH1D*)hPIDdata->Projection(kPidPt, "e");
  hMCmuonYield->SetName("hMCmuonYield");
  
  TCanvas* cYields = new TCanvas;
  cYields->SetLogx(kTRUE);
  cYields->SetLogy(kTRUE);
  hMCmuonYield->SetLineColor(kOrange);
  hMCmuonYield->GetXaxis()->SetRangeUser(0.15, 50.);
  hMCmuonYield->Draw();
  
  hPIDdata->GetAxis(kPidMCpid)->SetRange(1, 1);
  TH1D* hMCelectronYield = (TH1D*)hPIDdata->Projection(kPidPt, "e");
  hMCelectronYield->SetName("hMCelectronYield");
  
  hMCelectronYield->SetLineColor(kMagenta);
  hMCelectronYield->Draw("same");
  
  TH1D* hMCmuonToElRatio = new TH1D(*hMCmuonYield);
  hMCmuonToElRatio->Divide(hMCelectronYield);
  
  TCanvas* cFit = new TCanvas;
  cFit->SetLogx(kTRUE);
  hMCmuonToElRatio->GetYaxis()->SetRangeUser(0., 1.);
  hMCmuonToElRatio->Draw();
  
  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])",
        0.01, 50.);
  
  f.SetParameters(-0.688, 0.042, 4.52, -6.2, 0.53, 2.89, 0.035, 2.3, 6.);
  restrictJetPtAxis = kTRUE;
  if (jetHandling)
    f.FixParameter(8, 6.);
  
  f.SetParLimits(4, 0.2, 2.0);
  f.SetParLimits(7, 1.0, 6.0);
  
  hMCmuonToElRatio->Fit(&f, jetHandling ? "+" : "+W", "", 0.15, restrictJetPtAxis ? actualUpperJetPt : 15.0);
  
  return 0;
 extractMCmuonToElRatio.C:1
 extractMCmuonToElRatio.C:2
 extractMCmuonToElRatio.C:3
 extractMCmuonToElRatio.C:4
 extractMCmuonToElRatio.C:5
 extractMCmuonToElRatio.C:6
 extractMCmuonToElRatio.C:7
 extractMCmuonToElRatio.C:8
 extractMCmuonToElRatio.C:9
 extractMCmuonToElRatio.C:10
 extractMCmuonToElRatio.C:11
 extractMCmuonToElRatio.C:12
 extractMCmuonToElRatio.C:13
 extractMCmuonToElRatio.C:14
 extractMCmuonToElRatio.C:15
 extractMCmuonToElRatio.C:16
 extractMCmuonToElRatio.C:17
 extractMCmuonToElRatio.C:18
 extractMCmuonToElRatio.C:19
 extractMCmuonToElRatio.C:20
 extractMCmuonToElRatio.C:21
 extractMCmuonToElRatio.C:22
 extractMCmuonToElRatio.C:23
 extractMCmuonToElRatio.C:24
 extractMCmuonToElRatio.C:25
 extractMCmuonToElRatio.C:26
 extractMCmuonToElRatio.C:27
 extractMCmuonToElRatio.C:28
 extractMCmuonToElRatio.C:29
 extractMCmuonToElRatio.C:30
 extractMCmuonToElRatio.C:31
 extractMCmuonToElRatio.C:32
 extractMCmuonToElRatio.C:33
 extractMCmuonToElRatio.C:34
 extractMCmuonToElRatio.C:35
 extractMCmuonToElRatio.C:36
 extractMCmuonToElRatio.C:37
 extractMCmuonToElRatio.C:38
 extractMCmuonToElRatio.C:39
 extractMCmuonToElRatio.C:40
 extractMCmuonToElRatio.C:41
 extractMCmuonToElRatio.C:42
 extractMCmuonToElRatio.C:43
 extractMCmuonToElRatio.C:44
 extractMCmuonToElRatio.C:45
 extractMCmuonToElRatio.C:46
 extractMCmuonToElRatio.C:47
 extractMCmuonToElRatio.C:48
 extractMCmuonToElRatio.C:49
 extractMCmuonToElRatio.C:50
 extractMCmuonToElRatio.C:51
 extractMCmuonToElRatio.C:52
 extractMCmuonToElRatio.C:53
 extractMCmuonToElRatio.C:54
 extractMCmuonToElRatio.C:55
 extractMCmuonToElRatio.C:56
 extractMCmuonToElRatio.C:57
 extractMCmuonToElRatio.C:58
 extractMCmuonToElRatio.C:59
 extractMCmuonToElRatio.C:60
 extractMCmuonToElRatio.C:61
 extractMCmuonToElRatio.C:62
 extractMCmuonToElRatio.C:63
 extractMCmuonToElRatio.C:64
 extractMCmuonToElRatio.C:65
 extractMCmuonToElRatio.C:66
 extractMCmuonToElRatio.C:67
 extractMCmuonToElRatio.C:68
 extractMCmuonToElRatio.C:69
 extractMCmuonToElRatio.C:70
 extractMCmuonToElRatio.C:71
 extractMCmuonToElRatio.C:72
 extractMCmuonToElRatio.C:73
 extractMCmuonToElRatio.C:74
 extractMCmuonToElRatio.C:75
 extractMCmuonToElRatio.C:76
 extractMCmuonToElRatio.C:77
 extractMCmuonToElRatio.C:78
 extractMCmuonToElRatio.C:79
 extractMCmuonToElRatio.C:80
 extractMCmuonToElRatio.C:81
 extractMCmuonToElRatio.C:82
 extractMCmuonToElRatio.C:83
 extractMCmuonToElRatio.C:84
 extractMCmuonToElRatio.C:85
 extractMCmuonToElRatio.C:86
 extractMCmuonToElRatio.C:87
 extractMCmuonToElRatio.C:88
 extractMCmuonToElRatio.C:89
 extractMCmuonToElRatio.C:90
 extractMCmuonToElRatio.C:91
 extractMCmuonToElRatio.C:92
 extractMCmuonToElRatio.C:93
 extractMCmuonToElRatio.C:94
 extractMCmuonToElRatio.C:95
 extractMCmuonToElRatio.C:96
 extractMCmuonToElRatio.C:97
 extractMCmuonToElRatio.C:98
 extractMCmuonToElRatio.C:99
 extractMCmuonToElRatio.C:100
 extractMCmuonToElRatio.C:101
 extractMCmuonToElRatio.C:102
 extractMCmuonToElRatio.C:103
 extractMCmuonToElRatio.C:104
 extractMCmuonToElRatio.C:105
 extractMCmuonToElRatio.C:106
 extractMCmuonToElRatio.C:107
 extractMCmuonToElRatio.C:108
 extractMCmuonToElRatio.C:109
 extractMCmuonToElRatio.C:110
 extractMCmuonToElRatio.C:111
 extractMCmuonToElRatio.C:112
 extractMCmuonToElRatio.C:113
 extractMCmuonToElRatio.C:114
 extractMCmuonToElRatio.C:115
 extractMCmuonToElRatio.C:116
 extractMCmuonToElRatio.C:117
 extractMCmuonToElRatio.C:118
 extractMCmuonToElRatio.C:119
 extractMCmuonToElRatio.C:120
 extractMCmuonToElRatio.C:121
 extractMCmuonToElRatio.C:122
 extractMCmuonToElRatio.C:123
 extractMCmuonToElRatio.C:124
 extractMCmuonToElRatio.C:125
 extractMCmuonToElRatio.C:126
 extractMCmuonToElRatio.C:127
 extractMCmuonToElRatio.C:128
 extractMCmuonToElRatio.C:129
 extractMCmuonToElRatio.C:130
 extractMCmuonToElRatio.C:131
 extractMCmuonToElRatio.C:132
 extractMCmuonToElRatio.C:133
 extractMCmuonToElRatio.C:134
 extractMCmuonToElRatio.C:135
 extractMCmuonToElRatio.C:136
 extractMCmuonToElRatio.C:137
 extractMCmuonToElRatio.C:138
 extractMCmuonToElRatio.C:139
 extractMCmuonToElRatio.C:140
 extractMCmuonToElRatio.C:141
 extractMCmuonToElRatio.C:142
 extractMCmuonToElRatio.C:143
 extractMCmuonToElRatio.C:144
 extractMCmuonToElRatio.C:145
 extractMCmuonToElRatio.C:146
 extractMCmuonToElRatio.C:147
 extractMCmuonToElRatio.C:148
 extractMCmuonToElRatio.C:149
 extractMCmuonToElRatio.C:150
 extractMCmuonToElRatio.C:151
 extractMCmuonToElRatio.C:152
 extractMCmuonToElRatio.C:153
 extractMCmuonToElRatio.C:154
 extractMCmuonToElRatio.C:155
 extractMCmuonToElRatio.C:156
 extractMCmuonToElRatio.C:157
 extractMCmuonToElRatio.C:158
 extractMCmuonToElRatio.C:159
 extractMCmuonToElRatio.C:160
 extractMCmuonToElRatio.C:161
 extractMCmuonToElRatio.C:162
 extractMCmuonToElRatio.C:163