StFlowMephi/CalcRes.C

#include <iostream>

#include <TMath.h>
#include <TH1D.h>
#include <TProfile.h>
#include <TFile.h>
#include <Math/SpecFuncMathMore.h>

// Add libMathMore library for GetResZDC() & GetChiZDC()
R__LOAD_LIBRARY(libMathMore.so);

const int fNharmonics = 2;

Double_t GetRes(Double_t _chi)
{
  Double_t con = TMath::Sqrt(TMath::Pi() / 2) / 2;
  Double_t arg = _chi * _chi / 4.;
  Double_t res = con * _chi * exp(-arg) * (TMath::BesselI0(arg) + TMath::BesselI1(arg));
  return res;
}

Double_t GetResZDC(Double_t _chi, Int_t _harm)
{
  Double_t con = TMath::Sqrt(TMath::Pi() / 2) / 2;
  Double_t arg = _chi * _chi / 4.;
  if (_harm == 1)
  {
    Double_t res = con * _chi * exp(-arg) * (TMath::BesselI0(arg) + TMath::BesselI1(arg));
    return res;
  }
  else 
  {
    if (_harm == 2)
    {
      Double_t res = con * _chi * exp(-arg) * (ROOT::Math::cyl_bessel_i(0.5, arg) + ROOT::Math::cyl_bessel_i(1.5, arg));
      return res;
    }
    else
    {
      std::cout << "Wrong harmonic in resolution extrapolation! " << std::endl;
      return -999.;
    }
  }
}

Double_t GetChi(Double_t _res, Int_t accuracy = 50)
{
  Double_t chi = 2.0;
  Double_t delta = 1.0;
  for (int i = 0; i < accuracy; i++)
  {
    if (GetRes(chi) < _res)
    {
      chi = chi + delta;
    }
    else
    {
      chi = chi - delta;
    }
    delta = delta / 2.;
  }
  return chi;
}

Double_t GetChiZDC(Double_t _res, Int_t _harm, Int_t accuracy = 50)
{
  Double_t chi = 2.0;
  Double_t delta = 1.0;
  for (int i = 0; i < accuracy; i++)
  {
    if (GetResZDC(chi, _harm) < _res)
    {
      chi = chi + delta;
    }
    else
    {
      chi = chi - delta;
    }
    delta = delta / 2.;
  }
  return chi;
}

void CalcRes(const Char_t *resFileName = "./oResolutionTest.root", const Char_t *resFitName = "./oResFitTest.root")
{

  //const std::map<int,int> CentBins [] = {{1,75},{2,65},{3,55},{4,45},{5,35},{6,25},{7,15},{8,7.5},{9,2.5}};
  const int NGaps = 1;
  TProfile *pRes2Sqr[NGaps];
  TProfile *pRes3Sqr[NGaps];
  TProfile *pResSP2Sqr[NGaps];
  TProfile *pResSP3Sqr[NGaps];
  TProfile *pNormSP2Sqr[NGaps];
  TProfile *pNormSP3Sqr[NGaps];
  TProfile *pResBBC[fNharmonics + 1];
  TProfile *pResZDC[fNharmonics + 1];
  TProfile *p_res2_EP_3sub_AB;
  TProfile *p_res2_EP_3sub_AC;
  TProfile *p_res2_EP_3sub_BC;
  TProfile *p_res2_EP_3sub1_AB;
  TProfile *p_res2_EP_3sub1_AC;
  TProfile *p_res2_EP_3sub1_BC;

  TH1D *hResBBCHalf[fNharmonics + 1];
  TH1D *hResBBCFULL[fNharmonics + 1];

  TH1D *hResZDCHalf[fNharmonics + 1];
  TH1D *hResZDCFULL[fNharmonics + 1];

  TH1D *hRes2Half[NGaps];
  TH1D *hRes3Half[NGaps];
  TH1D *hRes2FULL[NGaps];
  TH1D *hRes3FULL[NGaps];
  TH1D *hRes2Chi[NGaps];
  TH1D *hRes3Chi[NGaps];

  TH1D *hResSP2Half[NGaps];
  TH1D *hResSP3Half[NGaps];
  TH1D *hResSP2FULL[NGaps];
  TH1D *hResSP3FULL[NGaps];
  TH1D *hResSP2Chi[NGaps];
  TH1D *hResSP3Chi[NGaps];

  TH1D *hResEP3subTPC;
  TH1D *hResEP3subBBCEast;
  TH1D *hResEP3subBBCWest;
  TH1D *hResEP3sub1BBC;
  TH1D *hResEP3sub1TPCEast;
  TH1D *hResEP3sub1TPCWest;

  TFile *fi = new TFile(resFileName, "read");
  fi->cd();

  for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  {
    pResBBC[iHarm] = (TProfile *)fi->Get(Form("p_resBBC_EP%i", iHarm));

    hResBBCHalf[iHarm] = new TH1D(Form("hResBBCHalf%i", iHarm), Form("hResBBCHalf%i", iHarm), 16, -0.5, 15.5);
    hResBBCFULL[iHarm] = new TH1D(Form("hResBBCFULL%i", iHarm), Form("hResBBCFULL%i", iHarm), 16, -0.5, 15.5);

    pResZDC[iHarm] = (TProfile *)fi->Get(Form("p_resZDC_EP%i", iHarm));

    hResZDCHalf[iHarm] = new TH1D(Form("hResZDCHalf%i", iHarm), Form("hResZDCHalf%i", iHarm), 16, -0.5, 15.5);
    hResZDCFULL[iHarm] = new TH1D(Form("hResZDCFULL%i", iHarm), Form("hResZDCFULL%i", iHarm), 16, -0.5, 15.5);
  }

  p_res2_EP_3sub_AB = (TProfile *)fi->Get("p_res2_EP_3sub_AB");
  p_res2_EP_3sub_AC = (TProfile *)fi->Get("p_res2_EP_3sub_AC");
  p_res2_EP_3sub_BC = (TProfile *)fi->Get("p_res2_EP_3sub_BC");

  p_res2_EP_3sub1_AB = (TProfile *)fi->Get("p_res2_EP_3sub1_AB");
  p_res2_EP_3sub1_AC = (TProfile *)fi->Get("p_res2_EP_3sub1_AC");
  p_res2_EP_3sub1_BC = (TProfile *)fi->Get("p_res2_EP_3sub1_BC");

  hResEP3subTPC = new TH1D("hResEP3subTPC", "hResEP3subTPC", 16, -0.5, 15.5);
  hResEP3subBBCEast = new TH1D("hResEP3subBBCEast", "hResEP3subBBCEast", 16, -0.5, 15.5);
  hResEP3subBBCWest = new TH1D("hResEP3subBBCWest", "hResEP3subBBCWest", 16, -0.5, 15.5);

  hResEP3sub1BBC = new TH1D("hResEP3sub1BBC", "hResEP3sub1BBC", 16, -0.5, 15.5);
  hResEP3sub1TPCEast = new TH1D("hResEP3sub1TPCEast", "hResEP3sub1TPCEast", 16, -0.5, 15.5);
  hResEP3sub1TPCWest = new TH1D("hResEP3sub1TPCWest", "hResEP3sub1TPCWest", 16, -0.5, 15.5);

  for (int i = 0; i < NGaps; i++)
  {
    pRes2Sqr[i] = (TProfile *)fi->Get(Form("p_res2_EP_gap%i", i));
    pRes3Sqr[i] = (TProfile *)fi->Get(Form("p_res3_EP_gap%i", i));

    pResSP2Sqr[i] = (TProfile *)fi->Get(Form("p_res2_SP_gap%i", i));
    pResSP3Sqr[i] = (TProfile *)fi->Get(Form("p_res3_SP_gap%i", i));
    pNormSP2Sqr[i] = (TProfile *)fi->Get(Form("p_res2_SP_Norm_gap%i", i));
    pNormSP3Sqr[i] = (TProfile *)fi->Get(Form("p_res3_SP_Norm_gap%i", i));

    hRes2Half[i] = new TH1D(Form("hRes2Half%i", i), Form("hRes2Half%i", i), 16, -0.5, 15.5);
    hRes3Half[i] = new TH1D(Form("hRes3Half%i", i), Form("hRes3Half%i", i), 16, -0.5, 15.5);
    hRes2FULL[i] = new TH1D(Form("hRes2FULL%i", i), Form("hRes2FULL%i", i), 16, -0.5, 15.5);
    hRes3FULL[i] = new TH1D(Form("hRes3FULL%i", i), Form("hRes3FULL%i", i), 16, -0.5, 15.5);
    hRes2Chi[i] = new TH1D(Form("hRes2Chi%i", i), Form("hRes2Chi%i", i), 16, -0.5, 15.5);
    hRes3Chi[i] = new TH1D(Form("hRes3Chi%i", i), Form("hRes3Chi%i", i), 16, -0.5, 15.5);

    hResSP2Half[i] = new TH1D(Form("hResSP2Half%i", i), Form("hResSP2Half%i", i), 16, -0.5, 15.5);
    hResSP3Half[i] = new TH1D(Form("hResSP3Half%i", i), Form("hResSP3Half%i", i), 16, -0.5, 15.5);
    hResSP2FULL[i] = new TH1D(Form("hResSP2FULL%i", i), Form("hResSP2FULL%i", i), 16, -0.5, 15.5);
    hResSP3FULL[i] = new TH1D(Form("hResSP3FULL%i", i), Form("hResSP3FULL%i", i), 16, -0.5, 15.5);
    hResSP2Chi[i] = new TH1D(Form("hResSP2Chi%i", i), Form("hResSP2Chi%i", i), 16, -0.5, 15.5);
    hResSP3Chi[i] = new TH1D(Form("hResSP3Chi%i", i), Form("hResSP3Chi%i", i), 16, -0.5, 15.5);
  }

  int bin = 0;
  Double_t AB, AC, BC, eAB, eAC, eBC;

  for (int ibin = 1; ibin <= p_res2_EP_3sub_AB->GetNbinsX(); ibin++)
  {
    AB = p_res2_EP_3sub_AB->GetBinContent(ibin);
    AC = p_res2_EP_3sub_AC->GetBinContent(ibin);
    BC = p_res2_EP_3sub_BC->GetBinContent(ibin);
    eAB = p_res2_EP_3sub_AB->GetBinError(ibin);
    eAC = p_res2_EP_3sub_AC->GetBinError(ibin);
    eBC = p_res2_EP_3sub_BC->GetBinError(ibin);

    // std::cout << "bin " << ibin << "AB: " << AB << " AC: " << AC << " BC: " << BC << "\t| ResSqr: " << (AB*AC/BC) << "\t| Res: " << TMath::Sqrt(AB*AC/BC) << std::endl;

    if (AB * AC / BC >= 0)
    {
      hResEP3subTPC->SetBinContent(ibin, TMath::Sqrt(AB * AC / BC));
      hResEP3subTPC->SetBinError(ibin, 0.5 * TMath::Sqrt(AC * eAB * eAB / (BC * AB) +
                                                         AB * eAC * eAC / (BC * AC) +
                                                         AB * AC * eBC * eBC / (BC * BC * BC)));
    }
    if (AB * BC / AC >= 0)
    {
      hResEP3subBBCEast->SetBinContent(ibin, TMath::Sqrt(AB * BC / AC));
      hResEP3subBBCEast->SetBinError(ibin, 0.5 * TMath::Sqrt(BC * eAB * eAB / (AC * AB) +
                                                             AB * eBC * eBC / (AC * BC) +
                                                             AB * BC * eAC * eAC / (AC * AC * AC)));
    }
    if (AC * BC / AB >= 0)
    {
      hResEP3subBBCWest->SetBinContent(ibin, TMath::Sqrt(AC * BC / AB));
      hResEP3subBBCWest->SetBinError(ibin, 0.5 * TMath::Sqrt(BC * eAC * eAC / (AB * AC) +
                                                             AC * eBC * eBC / (AB * BC) +
                                                             AC * BC * eAB * eAB / (AB * AB * AB)));
    }
  }

  for (int ibin = 1; ibin <= p_res2_EP_3sub1_AB->GetNbinsX(); ibin++)
  {
    AB = p_res2_EP_3sub1_AB->GetBinContent(ibin);
    AC = p_res2_EP_3sub1_AC->GetBinContent(ibin);
    BC = p_res2_EP_3sub1_BC->GetBinContent(ibin);
    eAB = p_res2_EP_3sub1_AB->GetBinError(ibin);
    eAC = p_res2_EP_3sub1_AC->GetBinError(ibin);
    eBC = p_res2_EP_3sub1_BC->GetBinError(ibin);

    // std::cout << "bin " << ibin << "AB: " << AB << " AC: " << AC << " BC: " << BC << "\t| ResSqr: " << (AB*AC/BC) << "\t| Res: " << TMath::Sqrt(AB*AC/BC) << std::endl;

    if (AB * AC / BC >= 0)
    {
      hResEP3sub1BBC->SetBinContent(ibin, TMath::Sqrt(AB * AC / BC));
      hResEP3sub1BBC->SetBinError(ibin, 0.5 * TMath::Sqrt(AC * eAB * eAB / (BC * AB) +
                                                          AB * eAC * eAC / (BC * AC) +
                                                          AB * AC * eBC * eBC / (BC * BC * BC)));
    }
    if (AB * BC / AC >= 0)
    {
      hResEP3sub1TPCEast->SetBinContent(ibin, TMath::Sqrt(AB * BC / AC));
      hResEP3sub1TPCEast->SetBinError(ibin, 0.5 * TMath::Sqrt(BC * eAB * eAB / (AC * AB) +
                                                              AB * eBC * eBC / (AC * BC) +
                                                              AB * BC * eAC * eAC / (AC * AC * AC)));
    }
    if (AC * BC / AB >= 0)
    {
      hResEP3sub1TPCWest->SetBinContent(ibin, TMath::Sqrt(AC * BC / AB));
      hResEP3sub1TPCWest->SetBinError(ibin, 0.5 * TMath::Sqrt(BC * eAC * eAC / (AB * AC) +
                                                              AC * eBC * eBC / (AB * BC) +
                                                              AC * BC * eAB * eAB / (AB * AB * AB)));
    }
  }

  //Fill ResHalf
  for (int i = 0; i < NGaps; i++)
  {
    for (int ibin = 1; ibin <= pRes2Sqr[i]->GetNbinsX(); ibin++)
    {
      hRes2Half[i]->SetBinContent(ibin, TMath::Sqrt(pRes2Sqr[i]->GetBinContent(ibin)));
      hRes2Half[i]->SetBinError(ibin, (Double_t)(0.5 * pRes2Sqr[i]->GetBinError(ibin) / TMath::Sqrt(pRes2Sqr[i]->GetBinContent(ibin))));

      hRes3Half[i]->SetBinContent(ibin, TMath::Sqrt(pRes3Sqr[i]->GetBinContent(ibin)));
      hRes3Half[i]->SetBinError(ibin, (Double_t)(0.5 * pRes3Sqr[i]->GetBinError(ibin) / TMath::Sqrt(pRes3Sqr[i]->GetBinContent(ibin))));
    }
    for (int ibin = 1; ibin <= pResSP2Sqr[i]->GetNbinsX(); ibin++)
    {
      if (pNormSP2Sqr[i]->GetBinContent(ibin) != 0 && pResSP2Sqr[i]->GetBinContent(ibin) > 0)
      {
        hResSP2Half[i]->SetBinContent(ibin, TMath::Sqrt(pResSP2Sqr[i]->GetBinContent(ibin) / pNormSP2Sqr[i]->GetBinContent(ibin)));
        hResSP2Half[i]->SetBinError(ibin, (Double_t)(0.5 * TMath::Sqrt(
                                                               TMath::Power(pResSP2Sqr[i]->GetBinError(ibin), 2) / (pResSP2Sqr[i]->GetBinContent(ibin) * pNormSP2Sqr[i]->GetBinContent(ibin)) +
                                                               pResSP2Sqr[i]->GetBinContent(ibin) * TMath::Power(pNormSP2Sqr[i]->GetBinError(ibin), 2) / TMath::Power(pNormSP2Sqr[i]->GetBinContent(ibin), 3))));
      }

      if (pNormSP3Sqr[i]->GetBinContent(ibin) != 0 && pResSP3Sqr[i]->GetBinContent(ibin) > 0)
      {
        hResSP3Half[i]->SetBinContent(ibin, TMath::Sqrt(pResSP3Sqr[i]->GetBinContent(ibin) / pNormSP3Sqr[i]->GetBinContent(ibin)));
        hResSP3Half[i]->SetBinError(ibin, (Double_t)(0.5 * TMath::Sqrt(
                                                               TMath::Power(pResSP3Sqr[i]->GetBinError(ibin), 2) / (pResSP3Sqr[i]->GetBinContent(ibin) * pNormSP3Sqr[i]->GetBinContent(ibin)) +
                                                               pResSP3Sqr[i]->GetBinContent(ibin) * TMath::Power(pNormSP3Sqr[i]->GetBinError(ibin), 2) / TMath::Power(pNormSP3Sqr[i]->GetBinContent(ibin), 3))));
      }
    }
  }
  for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  {
    for (int ibin = 1; ibin <= pResBBC[iHarm]->GetNbinsX(); ibin++)
    {
      hResBBCHalf[iHarm]->SetBinContent(ibin, TMath::Sqrt(TMath::Abs(pResBBC[iHarm]->GetBinContent(ibin))));
      hResBBCHalf[iHarm]->SetBinError(ibin, 0.5 * pResBBC[iHarm]->GetBinError(ibin) / TMath::Sqrt(TMath::Abs(pResBBC[iHarm]->GetBinContent(ibin))));
    }
    for (int ibin = 1; ibin <= pResZDC[iHarm]->GetNbinsX(); ibin++)
    {
      if (iHarm == 0)
      {
        hResZDCHalf[iHarm]->SetBinContent(ibin, TMath::Sqrt(TMath::Abs(pResZDC[iHarm]->GetBinContent(ibin))));
        hResZDCHalf[iHarm]->SetBinError(ibin, 0.5 * pResZDC[iHarm]->GetBinError(ibin) / TMath::Sqrt(TMath::Abs(pResZDC[iHarm]->GetBinContent(ibin))));
      }
    }
  }

  Double_t chi = 0.;

  //Fill Chi & FULL
  for (int i = 0; i < NGaps; i++)
  {
    for (int ibin = 1; ibin <= pRes2Sqr[i]->GetNbinsX(); ibin++)
    {
      chi = GetChi(hRes2Half[i]->GetBinContent(ibin));
      hRes2Chi[i]->SetBinContent(ibin, GetRes(chi));
      hRes2Chi[i]->SetBinError(ibin, hRes2Half[i]->GetBinError(ibin));
      hRes2FULL[i]->SetBinContent(ibin, GetRes(TMath::Sqrt(2) * chi));
      hRes2FULL[i]->SetBinError(ibin, TMath::Sqrt(2) * hRes2Half[i]->GetBinError(ibin));

      chi = GetChi(hRes3Half[i]->GetBinContent(ibin));
      hRes3Chi[i]->SetBinContent(ibin, GetRes(chi));
      hRes3Chi[i]->SetBinError(ibin, hRes3Half[i]->GetBinError(ibin));
      hRes3FULL[i]->SetBinContent(ibin, GetRes(TMath::Sqrt(2) * chi));
      hRes3FULL[i]->SetBinError(ibin, TMath::Sqrt(2) * hRes3Half[i]->GetBinError(ibin));
    }
    for (int ibin = 1; ibin <= pResSP2Sqr[i]->GetNbinsX(); ibin++)
    {
      chi = GetChi(hResSP2Half[i]->GetBinContent(ibin));
      hResSP2Chi[i]->SetBinContent(ibin, GetRes(chi));
      hResSP2Chi[i]->SetBinError(ibin, hResSP2Half[i]->GetBinError(ibin));
      hResSP2FULL[i]->SetBinContent(ibin, GetRes(TMath::Sqrt(2) * chi));
      hResSP2FULL[i]->SetBinError(ibin, TMath::Sqrt(2) * hResSP2Half[i]->GetBinError(ibin));

      chi = GetChi(hResSP3Half[i]->GetBinContent(ibin));
      hResSP3Chi[i]->SetBinContent(ibin, GetRes(chi));
      hResSP3Chi[i]->SetBinError(ibin, hResSP3Half[i]->GetBinError(ibin));
      hResSP3FULL[i]->SetBinContent(ibin, GetRes(TMath::Sqrt(2) * chi));
      hResSP3FULL[i]->SetBinError(ibin, TMath::Sqrt(2) * hResSP3Half[i]->GetBinError(ibin));
    }
  }
  for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  {
    for (int ibin = 1; ibin <= pResBBC[iHarm]->GetNbinsX(); ibin++)
    {
      chi = GetChi(hResBBCHalf[iHarm]->GetBinContent(ibin));
      hResBBCFULL[iHarm]->SetBinContent(ibin, GetRes(TMath::Sqrt(2) * chi));
      hResBBCFULL[iHarm]->SetBinError(ibin, TMath::Sqrt(2) * hResBBCHalf[iHarm]->GetBinError(ibin));
    }
  }

  for (int ibin = 1; ibin <= hResZDCHalf[0]->GetNbinsX(); ibin++)
  {
    chi = GetChiZDC(hResZDCHalf[0]->GetBinContent(ibin), 1);
    std::cout << "chi_1[" << ibin << "] = " << chi;
    hResZDCFULL[0]->SetBinContent(ibin, GetResZDC(TMath::Sqrt(2) * chi, 1));
    hResZDCFULL[0]->SetBinError(ibin, TMath::Sqrt(2) * hResZDCHalf[0]->GetBinError(ibin));
    std::cout << "\t Res_1[" << ibin << "] = " << GetResZDC(TMath::Sqrt(2) * chi, 1) << std::endl;
  }
  for (int ibin = 1; ibin <= hResZDCHalf[0]->GetNbinsX(); ibin++)
  {
    chi = GetChiZDC(hResZDCHalf[0]->GetBinContent(ibin), 1);
    std::cout << "chi_2[" << ibin << "] = " << chi;
    hResZDCFULL[1]->SetBinContent(ibin, GetResZDC(TMath::Sqrt(2) * chi, 2));
    hResZDCFULL[1]->SetBinError(ibin, GetResZDC(TMath::Sqrt(2) * chi, 2)/ hResZDCHalf[0]->GetBinContent(ibin) * hResZDCHalf[0]->GetBinError(ibin));
    std::cout << "\t Res_2[" << ibin << "] = " << GetResZDC(TMath::Sqrt(2) * chi, 2) << std::endl;
    hResZDCHalf[1]->SetBinContent(ibin, GetResZDC(chi, 2));
    hResZDCHalf[1]->SetBinError(ibin, GetResZDC(chi, 2)/ hResZDCHalf[0]->GetBinContent(ibin) * hResZDCHalf[0]->GetBinError(ibin));
  }

  std::cout << "chi_1[] = chi_2[]. If it isn't the case, check the code." << std::endl;

  TFile *fo = new TFile(resFitName, "recreate");
  fo->cd();

  for (int i = 0; i < NGaps; i++)
  {
    hRes2Half[i]->Write();
    hRes2Chi[i]->Write();
    hRes2FULL[i]->Write();
  }

  for (int i = 0; i < NGaps; i++)
  {
    hRes3Half[i]->Write();
    hRes3Chi[i]->Write();
    hRes3FULL[i]->Write();
  }

  for (int i = 0; i < NGaps; i++)
  {
    hResSP2Half[i]->Write();
    hResSP2Chi[i]->Write();
    hResSP2FULL[i]->Write();
  }

  for (int i = 0; i < NGaps; i++)
  {
    hResSP3Half[i]->Write();
    hResSP3Chi[i]->Write();
    hResSP3FULL[i]->Write();
  }
  for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  {
    hResBBCHalf[iHarm]->Write();
    hResBBCFULL[iHarm]->Write();
  }
  for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  {
    hResZDCHalf[iHarm]->Write();
    hResZDCFULL[iHarm]->Write();
  }
  hResEP3subTPC->Write();
  hResEP3subBBCEast->Write();
  hResEP3subBBCWest->Write();
  hResEP3sub1BBC->Write();
  hResEP3sub1TPCEast->Write();
  hResEP3sub1TPCWest->Write();
}