StFlowMephi/CalcNewXY.C

#include <iostream>
#include <TFile.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TF1.h>
#include <TGraphErrors.h>

const int Niter = 20;
const int NptBins = 15;
const Double_t ptBins[NptBins + 1] = {0.2, 0.4, 0.6, 0.8, 1., 1.2, 1.4, 1.6, 1.8, 2., 2.2, 2.4, 2.6, 2.8, 3., 3.2};

void CalcNewXY(const Char_t *inFileName, const Char_t *outFileName)
{
  TFile *fi = new TFile(inFileName, "read");
  fi->cd();

  TH2D *hNewXYCut[NptBins];
  TH1D *hNewX[NptBins], *hNewYPion[NptBins], *hNewYKaon[NptBins];

  TF1 *gausPionX[NptBins], *gausPionY[NptBins];
  TF1 *gausKaonX[NptBins], *gausKaonY[NptBins];
  TF1 *gausPiBgX[NptBins], *gausPiBgY[NptBins];
  TF1 *gausKaBgX[NptBins], *gausKaBgY[NptBins];

  TF1 *gausDoubleX[NptBins], *gausDoublePionY[NptBins], *gausDoubleKaonY[NptBins];

  TF1 *polPionX, *polKaonX, *polPionY, *polKaonY, *polMeanX;

  TGraphErrors *grSigmPionX, *grSigmPionY, *grSigmKaonX, *grSigmKaonY;
  TGraphErrors *grMeanKaonX;
  std::vector<Double_t> vXX, vXYPion, vXeX, vXeYPion, vXYKaon, vXeYKaon;
  std::vector<Double_t> vYX, vYYPion, vYeX, vYeYPion, vYYKaon, vYeYKaon;
  std::vector<Double_t> vMeanKaonX, vMeanKaoneX;

  for (int i = 0; i < NptBins; i++)
  {
    hNewXYCut[i] = (TH2D *)fi->Get(Form("hPionNewXYCut%i", i));

    gausPionX[i] = new TF1(Form("gausPionX%i", i), "gaus", -0.1, 0.1);
    gausKaonX[i] = new TF1(Form("gausKaonX%i", i), "gaus", 0.15, 0.3);
    gausPionY[i] = new TF1(Form("gausPionY%i", i), "gaus", -0.1, 0.1);
    gausKaonY[i] = new TF1(Form("gausKaonY%i", i), "gaus", -0.1, 0.1);

    gausPiBgX[i] = new TF1(Form("gausPiBgX%i", i), "gaus", -0.2, 0.2);
    gausKaBgX[i] = new TF1(Form("gausKaBgX%i", i), "gaus", 0.1, 0.4);
    gausPiBgY[i] = new TF1(Form("gausPiBgY%i", i), "gaus", -0.2, 0.2);
    gausKaBgY[i] = new TF1(Form("gausKaBgY%i", i), "gaus", -0.2, 0.2);

    gausDoubleX[i] = new TF1(Form("gausDoubleX%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)", -0.2, 0.8);
    gausDoublePionY[i] = new TF1(Form("gausDoublePionY%i", i), "gaus(0)+gaus(3)", -0.4, 0.4);
    gausDoubleKaonY[i] = new TF1(Form("gausDoubleKaonY%i", i), "gaus(0)+gaus(3)", -0.4, 0.4);
  }

  for (int i = 0; i < NptBins; i++)
  {
    hNewX[i] = (TH1D *)hNewXYCut[i]->ProjectionX(Form("hNewX%i", i));
    hNewYKaon[i] = (TH1D *)hNewXYCut[i]->ProjectionY(Form("hNewYKaon%i", i),(Int_t)hNewXYCut[i]->GetXaxis()->FindBin(0.16),(Int_t)hNewXYCut[i]->GetXaxis()->FindBin(0.5));
    hNewYPion[i] = (TH1D *)hNewXYCut[i]->ProjectionY(Form("hNewYPion%i", i),(Int_t)hNewXYCut[i]->GetXaxis()->FindBin(-0.2),(Int_t)hNewXYCut[i]->GetXaxis()->FindBin(0.2));
  }

  polPionX = new TF1("polPionX","pol3",0.2,3.2);
  polKaonX = new TF1("polKaonX","pol3",0.2,3.2);
  polPionY = new TF1("polPionY","pol3",0.2,3.2);
  polKaonY = new TF1("polKaonY","pol3",0.2,3.2);
  polMeanX = new TF1("polMeanX","pol3",0.2,2.);

  Double_t parX[12], parYPion[6], parYKaon[6];

  std::cout << "Fitting procedure begins" << std::endl;

  for (int i = 0; i < NptBins; i++)
  {
    std::cout << "pT-range " << i << std::endl;
    for (int j = 0; j < Niter; j++)
    {
      hNewX[i]->Fit(gausPionX[i], "RQ0");
      hNewX[i]->Fit(gausKaonX[i], "RQ0");
      hNewX[i]->Fit(gausPiBgX[i], "RQ0");
      hNewX[i]->Fit(gausKaBgX[i], "RQ0");

      hNewYPion[i]->Fit(gausPionY[i], "RQ0");
      hNewYKaon[i]->Fit(gausKaonY[i], "RQ0");
      hNewYPion[i]->Fit(gausPiBgY[i], "RQ0");
      hNewYKaon[i]->Fit(gausKaBgY[i], "RQ0");
    }
    gausPionX[i]->GetParameters(&parX[0]);
    gausKaonX[i]->GetParameters(&parX[3]);
    gausPiBgX[i]->GetParameters(&parX[6]);
    gausKaBgX[i]->GetParameters(&parX[9]);

    gausPionY[i]->GetParameters(&parYPion[0]);
    gausKaonY[i]->GetParameters(&parYKaon[0]);
    gausPiBgY[i]->GetParameters(&parYPion[3]);
    gausKaBgY[i]->GetParameters(&parYKaon[3]);

    gausDoubleX[i]->SetParameters(parX);
    gausDoublePionY[i]->SetParameters(parYPion);
    gausDoubleKaonY[i]->SetParameters(parYKaon);
    for (int j = 0; j < Niter; j++)
    {
      if (j != Niter - 1)
      {
        hNewX[i]->Fit(gausDoubleX[i], "RQ0");
        hNewYPion[i]->Fit(gausDoublePionY[i], "RQ0");
        hNewYKaon[i]->Fit(gausDoubleKaonY[i], "RQ0");
      }
      if (j == Niter - 1)
      {
        hNewX[i]->Fit(gausDoubleX[i], "R");
        hNewYPion[i]->Fit(gausDoublePionY[i], "R");
        hNewYKaon[i]->Fit(gausDoubleKaonY[i], "R");
      }
    }

    vXX.push_back(0.5 * (ptBins[i] + ptBins[i + 1]));
    vYX.push_back(0.5 * (ptBins[i] + ptBins[i + 1]));
    vXeX.push_back(0.);
    vYeX.push_back(0.);

    //Pion sigma (NewX)
    if ((gausDoubleX[i]->GetParameter(8) < 0 && gausDoubleX[i]->GetParameter(2) >= 0) || (gausDoubleX[i]->GetParameter(2) >= 0 && gausDoubleX[i]->GetParameter(2) < gausDoubleX[i]->GetParameter(8)))
    {
      vXYPion.push_back(gausDoubleX[i]->GetParameter(2));
      vXeYPion.push_back(gausDoubleX[i]->GetParError(2));
    }
    else
    {
      vXYPion.push_back(gausDoubleX[i]->GetParameter(8));
      vXeYPion.push_back(gausDoubleX[i]->GetParError(8));
    }
    //Kaon sigma (NewX)
    if (gausDoubleX[i]->GetParameter(11) < 0 || (gausDoubleX[i]->GetParameter(5) >= 0 && gausDoubleX[i]->GetParameter(5) < gausDoubleX[i]->GetParameter(11)))
    {
      vXYKaon.push_back(gausDoubleX[i]->GetParameter(5));
      vXeYKaon.push_back(gausDoubleX[i]->GetParError(5));
      vMeanKaonX.push_back(gausDoubleX[i]->GetParameter(4));
      vMeanKaoneX.push_back(gausDoubleX[i]->GetParError(4));
    }
    else
    {
      vXYKaon.push_back(gausDoubleX[i]->GetParameter(11));
      vXeYKaon.push_back(gausDoubleX[i]->GetParError(11));
      vMeanKaonX.push_back(gausDoubleX[i]->GetParameter(10));
      vMeanKaoneX.push_back(gausDoubleX[i]->GetParError(10));
    }

    //Pion sigma (NewY)
    if (gausDoublePionY[i]->GetParameter(5) < 0 || (gausDoublePionY[i]->GetParameter(2) >= 0 && gausDoublePionY[i]->GetParameter(2) < gausDoublePionY[i]->GetParameter(5)))
    {
      vYYPion.push_back(gausDoublePionY[i]->GetParameter(2));
      vYeYPion.push_back(gausDoublePionY[i]->GetParError(2));
    }
    else
    {
      vYYPion.push_back(gausDoublePionY[i]->GetParameter(5));
      vYeYPion.push_back(gausDoublePionY[i]->GetParError(5));
    }
    //Kaon sigma (NewY)
    if (gausDoubleKaonY[i]->GetParameter(5) < 0 || (gausDoubleKaonY[i]->GetParameter(2) >= 0 && gausDoubleKaonY[i]->GetParameter(2) < gausDoubleKaonY[i]->GetParameter(5)))
    {
      vYYKaon.push_back(gausDoubleKaonY[i]->GetParameter(2));
      vYeYKaon.push_back(gausDoubleKaonY[i]->GetParError(2));
    }
    else
    {
      vYYKaon.push_back(gausDoubleKaonY[i]->GetParameter(5));
      vYeYKaon.push_back(gausDoubleKaonY[i]->GetParError(5));
    }
  }
  grSigmPionX = new TGraphErrors(vXX.size(), &vXX[0], &vXYPion[0], &vXeX[0], &vXeYPion[0]);
  grSigmKaonX = new TGraphErrors(vXX.size(), &vXX[0], &vXYKaon[0], &vXeX[0], &vXeYKaon[0]);
  grSigmPionY = new TGraphErrors(vYX.size(), &vYX[0], &vYYPion[0], &vYeX[0], &vYeYPion[0]);
  grSigmKaonY = new TGraphErrors(vYX.size(), &vYX[0], &vYYKaon[0], &vYeX[0], &vYeYKaon[0]);

  grMeanKaonX = new TGraphErrors(vXX.size(), &vXX[0], &vMeanKaonX[0], &vXeX[0], &vMeanKaoneX[0]);
  
  grSigmPionX->SetName("grSigmPionX");
  grSigmKaonX->SetName("grSigmKaonX");
  grSigmPionY->SetName("grSigmPionY");
  grSigmKaonY->SetName("grSigmKaonY");

  grMeanKaonX->SetName("grMeanKaonX");

  // Removing negative points
  for (int i=0;i<vXX.size();i++)
  {
    if (grSigmPionX->GetY()[i] < 0 || (i!=0 && grSigmPionX->GetY()[i] > 5*grSigmPionX->GetY()[i-1]) ) {grSigmPionX->RemovePoint(i); i--;}
    if (grSigmKaonX->GetY()[i] < 0 || (i!=0 && grSigmKaonX->GetY()[i] > 5*grSigmKaonX->GetY()[i-1]) ) {grSigmKaonX->RemovePoint(i); i--;}
    if (grSigmPionY->GetY()[i] < 0 || (i!=0 && grSigmPionY->GetY()[i] > 5*grSigmPionY->GetY()[i-1]) ) {grSigmPionY->RemovePoint(i); i--;}
    if (grSigmKaonY->GetY()[i] < 0 || (i!=0 && grSigmKaonY->GetY()[i] > 5*grSigmKaonY->GetY()[i-1]) ) {grSigmKaonY->RemovePoint(i); i--;}
  }

  for (int j = 0; j < Niter; j++)
  {
    if (j != Niter - 1)
    {
      grSigmPionX->Fit(polPionX,"RQ0");
      grSigmKaonX->Fit(polKaonX,"RQ0");
      grSigmPionY->Fit(polPionY,"RQ0");
      grSigmKaonY->Fit(polKaonY,"RQ0");

      grMeanKaonX->Fit(polMeanX,"RQ0");
    }
    if (j == Niter - 1)
    {
      grSigmPionX->Fit(polPionX,"R");
      grSigmKaonX->Fit(polKaonX,"R");
      grSigmPionY->Fit(polPionY,"R");
      grSigmKaonY->Fit(polKaonY,"R");

      grMeanKaonX->Fit(polMeanX,"R");
    }
  }

  TFile *fo = new TFile(outFileName, "recreate");
  fo->cd();
  for (int i = 0; i < NptBins; i++)
  {
    hNewX[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    gausDoubleX[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    hNewYPion[i]->Write();
    hNewYKaon[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    gausDoublePionY[i]->Write();
    gausDoubleKaonY[i]->Write();
  }
  grSigmPionX->Write();
  grSigmKaonX->Write();
  grSigmPionY->Write();
  grSigmKaonY->Write();
  grMeanKaonX->Write();

  polPionX->Write();
  polKaonX->Write();
  polPionY->Write();
  polKaonY->Write();
  polMeanX->Write();
}