StFlowMephi/CalcPID.C

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

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};

const double pionMassSqr = TMath::Power(0.13957061, 2);
const double kaonMassSqr = TMath::Power(0.493677, 2);
const double protMassSqr = TMath::Power(0.938272081, 2);

const int Niter = 1; //20

Double_t student(Double_t *x, Double_t *par)
{
  return par[0] * TMath::Gamma((par[1] + 1) / 2) / (TMath::Gamma(par[1] / 2) * par[2] * TMath::Sqrt(TMath::Pi() * par[1])) *
         TMath::Power(1 + 1 / par[1] * TMath::Power((x[0] - par[3]) / par[2], 2), -(par[1] + 1) / 2);
}

Double_t student3(Double_t *x, Double_t *par)
{
  return par[0] * TMath::Gamma((par[1] + 1) / 2) / (TMath::Gamma(par[1] / 2) * par[2] * TMath::Sqrt(TMath::Pi() * par[1])) *
             TMath::Power(1 + 1 / par[1] * TMath::Power((x[0] - par[3]) / par[2], 2), -(par[1] + 1) / 2) +
         par[4] * TMath::Gamma((par[5] + 1) / 2) / (TMath::Gamma(par[5] / 2) * par[6] * TMath::Sqrt(TMath::Pi() * par[5])) *
             TMath::Power(1 + 1 / par[1] * TMath::Power((x[0] - par[7]) / par[6], 2), -(par[5] + 1) / 2) +
         par[8] * TMath::Gamma((par[9] + 1) / 2) / (TMath::Gamma(par[9] / 2) * par[10] * TMath::Sqrt(TMath::Pi() * par[9])) *
             TMath::Power(1 + 1 / par[9] * TMath::Power((x[0] - par[11]) / par[10], 2), -(par[9] + 1) / 2);
}

void Fit6gausMassSqr(const Char_t *inFileName, const Char_t *outFileName)
{
  TFile *fi = new TFile(inFileName, "read");

  TH2D *hNsigMsqrPion[NptBins];
  TH2D *hNsigMsqrKaon[NptBins];
  TH2D *hNsigMsqrProton[NptBins];

  TH1D *hMsqrPion[NptBins];
  TH1D *hMsqrKaon[NptBins];
  TH1D *hMsqrProton[NptBins];

  TF1 *gausPion[NptBins];
  TF1 *gausKaon[NptBins];
  TF1 *gausProton[NptBins];
  TF1 *gausPiBg[NptBins];
  TF1 *gausKaBg[NptBins];
  TF1 *gausProtBg[NptBins];

  TF1 *gausTriplePion[NptBins];
  TF1 *gausTripleKaon[NptBins];
  TF1 *gausTripleProton[NptBins];

  TH1D *hMsqrSigmaPion = new TH1D("hMsqrSigmaPion", "hMsqrSigmaPion", NptBins, 0.2, 3.2);
  TH1D *hMsqrSigmaKaon = new TH1D("hMsqrSigmaKaon", "hMsqrSigmaKaon", NptBins, 0.2, 3.2);
  TH1D *hMsqrSigmaProton = new TH1D("hMsqrSigmaProton", "hMsqrSigmaProton", NptBins, 0.2, 3.2);

  for (int i = 0; i < NptBins; i++)
  {
    hNsigMsqrPion[i] = (TH2D *)fi->Get(Form("hNsigPionMSqr%i", i));
    hNsigMsqrKaon[i] = (TH2D *)fi->Get(Form("hNsigKaonMSqr%i", i));
    hNsigMsqrProton[i] = (TH2D *)fi->Get(Form("hNsigProtonMSqr%i", i));

    gausPion[i] = new TF1(Form("gausPion%i", i), "gaus", pionMassSqr * 0.5, pionMassSqr * 1.5);
    gausKaon[i] = new TF1(Form("gausKaon%i", i), "gaus", kaonMassSqr * 0.5, kaonMassSqr * 1.5);
    gausProton[i] = new TF1(Form("gausProton%i", i), "gaus", protMassSqr * 0.5, protMassSqr * 1.5);
    gausPiBg[i] = new TF1(Form("gausPiBg%i", i), "gaus", pionMassSqr * 0.2, pionMassSqr * 1.8);
    gausKaBg[i] = new TF1(Form("gausKaBg%i", i), "gaus", kaonMassSqr * 0.2, kaonMassSqr * 1.8);
    gausProtBg[i] = new TF1(Form("gausProtBg%i", i), "gaus", protMassSqr * 0.2, protMassSqr * 1.8);

    gausTriplePion[i] = new TF1(Form("gausTriplePion%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)+gaus(15)", -0.5, 1.5);
    gausTripleKaon[i] = new TF1(Form("gausTripleKaon%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)+gaus(15)", -0.5, 1.5);
    gausTripleProton[i] = new TF1(Form("gausTripleProton%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)+gaus(15)", -0.5, 1.5);
  }

  int bin1 = 0, bin2 = 0;
  Double_t parPion[18], parKaon[18], parProton[18];

  for (int i = 0; i < NptBins; i++)
  {
    bin1 = hNsigMsqrPion[i]->GetXaxis()->FindBin(-3.); // 651;
    bin2 = hNsigMsqrPion[i]->GetXaxis()->FindBin(3.);  // 750;
    hMsqrPion[i] = (TH1D *)hNsigMsqrPion[i]->ProjectionY(Form("hMsqrPion%i", i), bin1, bin2);
    hMsqrKaon[i] = (TH1D *)hNsigMsqrKaon[i]->ProjectionY(Form("hMsqrKaon%i", i), bin1, bin2);
    hMsqrProton[i] = (TH1D *)hNsigMsqrProton[i]->ProjectionY(Form("hMsqrProton%i", i), bin1, bin2);
  }

  //Pion
  for (int i = 0; i < NptBins; i++)
  {
    std::cout << "Pions, Pt range: " << i << std::endl;
    for (int j=0;j<Niter;j++)
    {
      hMsqrPion[i]->Fit(gausPion[i], "RQ0");
      hMsqrPion[i]->Fit(gausKaon[i], "RQ0");
      hMsqrPion[i]->Fit(gausProton[i], "RQ0");
      hMsqrPion[i]->Fit(gausPiBg[i], "RQ0");
      hMsqrPion[i]->Fit(gausKaBg[i], "RQ0");
      hMsqrPion[i]->Fit(gausProtBg[i], "RQ0");
    }

    gausPion[i]->GetParameters(&parPion[0]);
    gausKaon[i]->GetParameters(&parPion[3]);
    gausProton[i]->GetParameters(&parPion[6]);
    gausPiBg[i]->GetParameters(&parPion[9]);
    gausKaBg[i]->GetParameters(&parPion[12]);
    gausProtBg[i]->GetParameters(&parPion[15]);

    gausTriplePion[i]->SetParameters(parPion);

    for (int j=0;j<Niter;j++)
    {
      if (j != Niter-1) hMsqrPion[i]->Fit(gausTriplePion[i], "RQ0");
      if (j == Niter-1) hMsqrPion[i]->Fit(gausTriplePion[i], "R");
    }
    if (gausTriplePion[i]->GetParameter(2) > 0 && 
        (gausTriplePion[i]->GetParameter(2) < gausTriplePion[i]->GetParameter(11) ||
        gausTriplePion[i]->GetParameter(11) < 0) )
    {
      hMsqrSigmaPion->SetBinContent(i+1, gausTriplePion[i]->GetParameter(2));
      hMsqrSigmaPion->SetBinError(i+1, gausTriplePion[i]->GetParError(2));
    }
    else
    {
      hMsqrSigmaPion->SetBinContent(i+1, gausTriplePion[i]->GetParameter(11));
      hMsqrSigmaPion->SetBinError(i+1, gausTriplePion[i]->GetParError(11));
    }
    
  }

  //Kaon
  for (int i = 0; i < NptBins; i++)
  {
    std::cout << "Kaons, Pt range: " << i << std::endl;
    for (int j=0;j<Niter;j++)
    {
      hMsqrKaon[i]->Fit(gausPion[i], "RQ0");
      hMsqrKaon[i]->Fit(gausKaon[i], "RQ0");
      hMsqrKaon[i]->Fit(gausProton[i], "RQ0");
      hMsqrKaon[i]->Fit(gausPiBg[i], "RQ0");
      hMsqrKaon[i]->Fit(gausKaBg[i], "RQ0");
      hMsqrKaon[i]->Fit(gausProtBg[i], "RQ0");
    }

    gausPion[i]->GetParameters(&parKaon[0]);
    gausKaon[i]->GetParameters(&parKaon[3]);
    gausProton[i]->GetParameters(&parKaon[6]);
    gausPiBg[i]->GetParameters(&parKaon[9]);
    gausKaBg[i]->GetParameters(&parKaon[12]);
    gausProtBg[i]->GetParameters(&parKaon[15]);

    gausTripleKaon[i]->SetParameters(parKaon);

    for (int j=0;j<Niter;j++)
    {
      if (j != Niter-1) hMsqrKaon[i]->Fit(gausTripleKaon[i], "RQ0");
      if (j == Niter-1) hMsqrKaon[i]->Fit(gausTripleKaon[i], "R");
    }
    if (gausTripleKaon[i]->GetParameter(5) > 0 && 
        (gausTripleKaon[i]->GetParameter(5) < gausTripleKaon[i]->GetParameter(14) ||
        gausTripleKaon[i]->GetParameter(14) < 0) )
    {
      hMsqrSigmaKaon->SetBinContent(i+1, gausTripleKaon[i]->GetParameter(5));
      hMsqrSigmaKaon->SetBinError(i+1, gausTripleKaon[i]->GetParError(5));
    }
    else
    {
      hMsqrSigmaKaon->SetBinContent(i+1, gausTripleKaon[i]->GetParameter(14));
      hMsqrSigmaKaon->SetBinError(i+1, gausTripleKaon[i]->GetParError(14));
    }
  }

  //Proton
  for (int i = 0; i < NptBins; i++)
  {
    std::cout << "Protons, Pt range: " << i << std::endl;
    for (int j=0;j<Niter;j++)
    {
      hMsqrProton[i]->Fit(gausPion[i], "RQ0");
      hMsqrProton[i]->Fit(gausKaon[i], "RQ0");
      hMsqrProton[i]->Fit(gausProton[i], "RQ0");
      hMsqrProton[i]->Fit(gausPiBg[i], "RQ0");
      hMsqrProton[i]->Fit(gausKaBg[i], "RQ0");
      hMsqrProton[i]->Fit(gausProtBg[i], "RQ0");
    }

    gausPion[i]->GetParameters(&parProton[0]);
    gausKaon[i]->GetParameters(&parProton[3]);
    gausProton[i]->GetParameters(&parProton[6]);
    gausPiBg[i]->GetParameters(&parProton[9]);
    gausKaBg[i]->GetParameters(&parProton[12]);
    gausProtBg[i]->GetParameters(&parProton[15]);

    gausTripleProton[i]->SetParameters(parProton);

    for (int j=0;j<Niter;j++)
    {
      if (j != Niter-1) hMsqrProton[i]->Fit(gausTripleProton[i], "RQ0");
      if (j == Niter-1) hMsqrProton[i]->Fit(gausTripleProton[i], "R");
    }
    if (gausTripleProton[i]->GetParameter(8) > 0 && 
        (gausTripleProton[i]->GetParameter(8) < gausTripleProton[i]->GetParameter(17) ||
         gausTripleProton[i]->GetParameter(17) < 0) )
    {
      hMsqrSigmaProton->SetBinContent(i+1, gausTripleProton[i]->GetParameter(8));
      hMsqrSigmaProton->SetBinError(i+1, gausTripleProton[i]->GetParError(8));
    }
    else
    {
      hMsqrSigmaProton->SetBinContent(i+1, gausTripleProton[i]->GetParameter(17));
      hMsqrSigmaProton->SetBinError(i+1, gausTripleProton[i]->GetParError(17));
    }
  }

  TF1 *polPion = new TF1("polPion", "pol3", 0.4, 3.2);
  TF1 *polKaon = new TF1("polKaon", "pol3", 0.2, 3.2);
  TF1 *polProton = new TF1("polProton", "pol3", 0.2, 3.2);

  for (int j=0;j<Niter;j++)
  {
    if (j == Niter-1)
    {
      hMsqrSigmaPion->Fit(polPion, "R");
      hMsqrSigmaKaon->Fit(polKaon, "R");
      hMsqrSigmaProton->Fit(polProton, "R");
    }
    else
    {
      hMsqrSigmaPion->Fit(polPion, "R0Q");
      hMsqrSigmaKaon->Fit(polKaon, "R0Q");
      hMsqrSigmaProton->Fit(polProton, "R0Q");
    }
    
  }

  TFile *fo = new TFile(outFileName, "recreate");
  fo->cd();
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Write();
    hMsqrKaon[i]->Write();
    hMsqrProton[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    gausTriplePion[i]->Write();
    gausTripleKaon[i]->Write();
    gausTripleProton[i]->Write();
  }
  hMsqrSigmaPion->Write();
  hMsqrSigmaKaon->Write();
  hMsqrSigmaProton->Write();

  polPion->Write();
  polKaon->Write();
  polProton->Write();
}

void Fit5gausMassSqr(const Char_t *inFileName, const Char_t *outFileName)
{
  TFile *fi = new TFile(inFileName, "read");

  TH2D *hNsigMsqrPion[NptBins];
  TH2D *hNsigMsqrKaon[NptBins];
  TH2D *hNsigMsqrProton[NptBins];

  TH1D *hMsqrPion[NptBins];
  TH1D *hMsqrKaon[NptBins];
  TH1D *hMsqrProton[NptBins];

  TF1 *gausPion[NptBins];
  TF1 *gausKaon[NptBins];
  TF1 *gausProton[NptBins];
  TF1 *gausMesonBg[NptBins];
  TF1 *gausBosonBg[NptBins];

  TF1 *gausTriplePion[NptBins];
  TF1 *gausTripleKaon[NptBins];
  TF1 *gausTripleProton[NptBins];

  TH1D *hMsqrSigmaPion = new TH1D("hMsqrSigmaPion", "hMsqrSigmaPion", NptBins, 0.15, 3.75);
  TH1D *hMsqrSigmaKaon = new TH1D("hMsqrSigmaKaon", "hMsqrSigmaKaon", NptBins, 0.15, 3.75);
  TH1D *hMsqrSigmaProton = new TH1D("hMsqrSigmaProton", "hMsqrSigmaProton", NptBins, 0.15, 3.75);

  for (int i = 0; i < NptBins; i++)
  {
    hNsigMsqrPion[i] = (TH2D *)fi->Get(Form("hNsigPionMSqr%i", i));
    hNsigMsqrKaon[i] = (TH2D *)fi->Get(Form("hNsigKaonMSqr%i", i));
    hNsigMsqrProton[i] = (TH2D *)fi->Get(Form("hNsigProtonMSqr%i", i));

    gausPion[i] = new TF1(Form("gausPion%i", i), "gaus", pionMassSqr * 0., pionMassSqr * 2.);
    gausKaon[i] = new TF1(Form("gausKaon%i", i), "gaus", kaonMassSqr * 0.5, kaonMassSqr * 1.5);
    gausProton[i] = new TF1(Form("gausProton%i", i), "gaus", protMassSqr * 0.5, protMassSqr * 1.5);
    gausMesonBg[i] = new TF1(Form("gausMesonBg%i", i), "gaus", pionMassSqr * 0., kaonMassSqr * 1.5);
    gausBosonBg[i] = new TF1(Form("gausBosonBg%i", i), "gaus", kaonMassSqr * 1.5, protMassSqr * 1.5);

    gausTriplePion[i] = new TF1(Form("gausTriplePion%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)", -0.5, 1.5);
    gausTripleKaon[i] = new TF1(Form("gausTripleKaon%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)", -0.5, 1.5);
    gausTripleProton[i] = new TF1(Form("gausTripleProton%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)", -0.5, 1.5);
  }

  int bin1 = 0, bin2 = 0;
  Double_t parPion[15], parKaon[15], parProton[15];

  for (int i = 0; i < NptBins; i++)
  {
    bin1 = hNsigMsqrPion[i]->GetXaxis()->FindBin(-3.); // 651;
    bin2 = hNsigMsqrPion[i]->GetXaxis()->FindBin(3.);  // 750;
    hMsqrPion[i] = (TH1D *)hNsigMsqrPion[i]->ProjectionY(Form("hMsqrPion%i", i), bin1, bin2);
    hMsqrKaon[i] = (TH1D *)hNsigMsqrKaon[i]->ProjectionY(Form("hMsqrKaon%i", i), bin1, bin2);
    hMsqrProton[i] = (TH1D *)hNsigMsqrProton[i]->ProjectionY(Form("hMsqrProton%i", i), bin1, bin2);
  }

  //Pion
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Fit(gausPion[i], "RQ0");
    hMsqrPion[i]->Fit(gausKaon[i], "RQ0");
    hMsqrPion[i]->Fit(gausProton[i], "RQ0");
    hMsqrPion[i]->Fit(gausMesonBg[i], "RQ0");
    hMsqrPion[i]->Fit(gausBosonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parPion[0]);
    gausKaon[i]->GetParameters(&parPion[3]);
    gausProton[i]->GetParameters(&parPion[6]);
    gausMesonBg[i]->GetParameters(&parPion[9]);
    gausBosonBg[i]->GetParameters(&parPion[12]);

    gausTriplePion[i]->SetParameters(parPion);

    hMsqrPion[i]->Fit(gausTriplePion[i], "R");
    hMsqrSigmaPion->SetBinContent(i, gausTriplePion[i]->GetParameter(2));
    hMsqrSigmaPion->SetBinError(i, gausTriplePion[i]->GetParError(2));
  }

  //Kaon
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrKaon[i]->Fit(gausPion[i], "RQ0");
    hMsqrKaon[i]->Fit(gausKaon[i], "RQ0");
    hMsqrKaon[i]->Fit(gausProton[i], "RQ0");
    hMsqrKaon[i]->Fit(gausMesonBg[i], "RQ0");
    hMsqrKaon[i]->Fit(gausBosonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parKaon[0]);
    gausKaon[i]->GetParameters(&parKaon[3]);
    gausProton[i]->GetParameters(&parKaon[6]);
    gausMesonBg[i]->GetParameters(&parKaon[9]);
    gausBosonBg[i]->GetParameters(&parKaon[12]);

    gausTripleKaon[i]->SetParameters(parKaon);

    hMsqrKaon[i]->Fit(gausTripleKaon[i], "R");
    hMsqrSigmaKaon->SetBinContent(i, gausTripleKaon[i]->GetParameter(5));
    hMsqrSigmaKaon->SetBinError(i, gausTripleKaon[i]->GetParError(5));
  }

  //Proton
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrProton[i]->Fit(gausPion[i], "RQ0");
    hMsqrProton[i]->Fit(gausKaon[i], "RQ0");
    hMsqrProton[i]->Fit(gausProton[i], "RQ0");
    hMsqrProton[i]->Fit(gausMesonBg[i], "RQ0");
    hMsqrProton[i]->Fit(gausBosonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parProton[0]);
    gausKaon[i]->GetParameters(&parProton[3]);
    gausProton[i]->GetParameters(&parProton[6]);
    gausMesonBg[i]->GetParameters(&parProton[9]);
    gausBosonBg[i]->GetParameters(&parProton[12]);

    gausTripleProton[i]->SetParameters(parProton);

    hMsqrProton[i]->Fit(gausTripleProton[i], "R");
    hMsqrSigmaProton->SetBinContent(i, gausTripleProton[i]->GetParameter(8));
    hMsqrSigmaProton->SetBinError(i, gausTripleProton[i]->GetParError(8));
  }

  TF1 *polPion = new TF1("polPion", "pol8", 0.15, 3.7);
  TF1 *polKaon = new TF1("polKaon", "pol8", 0.15, 3.);
  TF1 *polProton = new TF1("polProton", "pol8", 0.3, 3.7);

  hMsqrSigmaPion->Fit(polPion, "R");
  hMsqrSigmaKaon->Fit(polKaon, "R");
  hMsqrSigmaProton->Fit(polProton, "R");

  TFile *fo = new TFile(outFileName, "recreate");
  fo->cd();
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Write();
    hMsqrKaon[i]->Write();
    hMsqrProton[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    gausTriplePion[i]->Write();
    gausTripleKaon[i]->Write();
    gausTripleProton[i]->Write();
  }
  hMsqrSigmaPion->Write();
  hMsqrSigmaKaon->Write();
  hMsqrSigmaProton->Write();

  polPion->Write();
  polKaon->Write();
  polProton->Write();
}

void Fit4gausMassSqr(const Char_t *inFileName, const Char_t *outFileName)
{
  TFile *fi = new TFile(inFileName, "read");

  TH2D *hNsigMsqrPion[NptBins];
  TH2D *hNsigMsqrKaon[NptBins];
  TH2D *hNsigMsqrProton[NptBins];

  TH1D *hMsqrPion[NptBins];
  TH1D *hMsqrKaon[NptBins];
  TH1D *hMsqrProton[NptBins];

  TF1 *gausPion[NptBins];
  TF1 *gausKaon[NptBins];
  TF1 *gausProton[NptBins];
  TF1 *gausMesonBg[NptBins];

  TF1 *gausTriplePion[NptBins];
  TF1 *gausTripleKaon[NptBins];
  TF1 *gausTripleProton[NptBins];

<<<<<<< HEAD
  TH1D *hMsqrSigmaPion = new TH1D("hMsqrSigmaPion", "hMsqrSigmaPion", NptBins, 0.15, 3.75);
  TH1D *hMsqrSigmaKaon = new TH1D("hMsqrSigmaKaon", "hMsqrSigmaKaon", NptBins, 0.15, 3.75);
  TH1D *hMsqrSigmaProton = new TH1D("hMsqrSigmaProton", "hMsqrSigmaProton", NptBins, 0.15, 3.75);
=======
  TH1D *hMsqrSigmaPion = new TH1D("hMsqrSigmaPion","hMsqrSigmaPion",NptBins,0.15,6.15);
  TH1D *hMsqrSigmaKaon = new TH1D("hMsqrSigmaKaon","hMsqrSigmaKaon",NptBins,0.15,6.15);
  TH1D *hMsqrSigmaProton = new TH1D("hMsqrSigmaProton","hMsqrSigmaProton",NptBins,0.15,6.15);
>>>>>>> 55931e7eab0e390e397eb65905a5af66aea18769

  for (int i = 0; i < NptBins; i++)
  {
    hNsigMsqrPion[i] = (TH2D *)fi->Get(Form("hNsigPionMSqr%i", i));
    hNsigMsqrKaon[i] = (TH2D *)fi->Get(Form("hNsigKaonMSqr%i", i));
    hNsigMsqrProton[i] = (TH2D *)fi->Get(Form("hNsigProtonMSqr%i", i));

    gausPion[i] = new TF1(Form("gausPion%i", i), "gaus", pionMassSqr * 0., pionMassSqr * 2.);
    gausKaon[i] = new TF1(Form("gausKaon%i", i), "gaus", kaonMassSqr * 0.5, kaonMassSqr * 1.5);
    gausProton[i] = new TF1(Form("gausProton%i", i), "gaus", protMassSqr * 0.5, protMassSqr * 1.5);
    gausMesonBg[i] = new TF1(Form("gausMesonBg%i", i), "gaus", pionMassSqr * 0., kaonMassSqr * 1.5);

    gausTriplePion[i] = new TF1(Form("gausTriplePion%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)", -0.5, 1.5);
    gausTripleKaon[i] = new TF1(Form("gausTripleKaon%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)", -0.5, 1.5);
    gausTripleProton[i] = new TF1(Form("gausTripleProton%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)", -0.5, 1.5);
  }

  int bin1 = 0, bin2 = 0;
  Double_t parPion[15], parKaon[15], parProton[15];

  for (int i = 0; i < NptBins; i++)
  {
<<<<<<< HEAD
    bin1 = hNsigMsqrPion[i]->GetXaxis()->FindBin(-3.); // 651;
    bin2 = hNsigMsqrPion[i]->GetXaxis()->FindBin(3.);  // 750;
    hMsqrPion[i] = (TH1D *)hNsigMsqrPion[i]->ProjectionY(Form("hMsqrPion%i", i), bin1, bin2);
    hMsqrKaon[i] = (TH1D *)hNsigMsqrKaon[i]->ProjectionY(Form("hMsqrKaon%i", i), bin1, bin2);
    hMsqrProton[i] = (TH1D *)hNsigMsqrProton[i]->ProjectionY(Form("hMsqrProton%i", i), bin1, bin2);
=======
    bin1 = 651; //hNsigMsqrPion[i]->FindBin(-2.5);
    bin2 = 750; //hNsigMsqrPion[i]->FindBin(2.5);
    hMsqrPion[i] = (TH1D*) hNsigMsqrPion[i]->ProjectionY(Form("hMsqrPion%i",i),bin1,bin2);
    hMsqrKaon[i] = (TH1D*) hNsigMsqrKaon[i]->ProjectionY(Form("hMsqrKaon%i",i),bin1,bin2);
    hMsqrProton[i] = (TH1D*) hNsigMsqrProton[i]->ProjectionY(Form("hMsqrProton%i",i),bin1,bin2);
>>>>>>> 55931e7eab0e390e397eb65905a5af66aea18769
  }

  //Pion
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Fit(gausPion[i], "RQ0");
    hMsqrPion[i]->Fit(gausKaon[i], "RQ0");
    hMsqrPion[i]->Fit(gausProton[i], "RQ0");
    hMsqrPion[i]->Fit(gausMesonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parPion[0]);
    gausKaon[i]->GetParameters(&parPion[3]);
    gausProton[i]->GetParameters(&parPion[6]);
    gausMesonBg[i]->GetParameters(&parPion[9]);

    gausTriplePion[i]->SetParameters(parPion);

    hMsqrPion[i]->Fit(gausTriplePion[i], "R");
    hMsqrSigmaPion->SetBinContent(i, gausTriplePion[i]->GetParameter(2));
    hMsqrSigmaPion->SetBinError(i, gausTriplePion[i]->GetParError(2));
  }

  //Kaon
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrKaon[i]->Fit(gausPion[i], "RQ0");
    hMsqrKaon[i]->Fit(gausKaon[i], "RQ0");
    hMsqrKaon[i]->Fit(gausProton[i], "RQ0");
    hMsqrKaon[i]->Fit(gausMesonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parKaon[0]);
    gausKaon[i]->GetParameters(&parKaon[3]);
    gausProton[i]->GetParameters(&parKaon[6]);
    gausMesonBg[i]->GetParameters(&parKaon[9]);

    gausTripleKaon[i]->SetParameters(parKaon);

    hMsqrKaon[i]->Fit(gausTripleKaon[i], "R");
    hMsqrSigmaKaon->SetBinContent(i, gausTripleKaon[i]->GetParameter(5));
    hMsqrSigmaKaon->SetBinError(i, gausTripleKaon[i]->GetParError(5));
  }

  //Proton
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrProton[i]->Fit(gausPion[i], "RQ0");
    hMsqrProton[i]->Fit(gausKaon[i], "RQ0");
    hMsqrProton[i]->Fit(gausProton[i], "RQ0");
    hMsqrProton[i]->Fit(gausMesonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parProton[0]);
    gausKaon[i]->GetParameters(&parProton[3]);
    gausProton[i]->GetParameters(&parProton[6]);
    gausMesonBg[i]->GetParameters(&parProton[9]);

    gausTripleProton[i]->SetParameters(parProton);

    hMsqrProton[i]->Fit(gausTripleProton[i], "R");
    hMsqrSigmaProton->SetBinContent(i, gausTripleProton[i]->GetParameter(8));
    hMsqrSigmaProton->SetBinError(i, gausTripleProton[i]->GetParError(8));
  }

  TF1 *polPion = new TF1("polPion", "pol8", 0.15, 3.7);
  TF1 *polKaon = new TF1("polKaon", "pol8", 0.15, 3.);
  TF1 *polProton = new TF1("polProton", "pol8", 0.3, 3.7);

  hMsqrSigmaPion->Fit(polPion, "R");
  hMsqrSigmaKaon->Fit(polKaon, "R");
  hMsqrSigmaProton->Fit(polProton, "R");

  TFile *fo = new TFile(outFileName, "recreate");
  fo->cd();
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Write();
    hMsqrKaon[i]->Write();
    hMsqrProton[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    gausTriplePion[i]->Write();
    gausTripleKaon[i]->Write();
    gausTripleProton[i]->Write();
  }
  hMsqrSigmaPion->Write();
  hMsqrSigmaKaon->Write();
  hMsqrSigmaProton->Write();

  polPion->Write();
  polKaon->Write();
  polProton->Write();
}

void Fit3gausMassSqr(const Char_t *inFileName, const Char_t *outFileName)
{
  TFile *fi = new TFile(inFileName, "read");

  TH2D *hNsigMsqrPion[NptBins];
  TH2D *hNsigMsqrKaon[NptBins];
  TH2D *hNsigMsqrProton[NptBins];

  TH1D *hMsqrPion[NptBins];
  TH1D *hMsqrKaon[NptBins];
  TH1D *hMsqrProton[NptBins];

  TF1 *gausPion[NptBins];
  TF1 *gausKaon[NptBins];
  TF1 *gausProton[NptBins];

  TF1 *gausTriplePion[NptBins];
  TF1 *gausTripleKaon[NptBins];
  TF1 *gausTripleProton[NptBins];

  TH1D *hMsqrSigmaPion = new TH1D("hMsqrSigmaPion", "hMsqrSigmaPion", NptBins, 0.2, 3.2);
  TH1D *hMsqrSigmaKaon = new TH1D("hMsqrSigmaKaon", "hMsqrSigmaKaon", NptBins, 0.2, 3.2);
  TH1D *hMsqrSigmaProton = new TH1D("hMsqrSigmaProton", "hMsqrSigmaProton", NptBins, 0.2, 3.2);

  for (int i = 0; i < NptBins; i++)
  {
    hNsigMsqrPion[i] = (TH2D *)fi->Get(Form("hNsigPionMSqr%i", i));
    hNsigMsqrKaon[i] = (TH2D *)fi->Get(Form("hNsigKaonMSqr%i", i));
    hNsigMsqrProton[i] = (TH2D *)fi->Get(Form("hNsigProtonMSqr%i", i));

    gausPion[i] = new TF1(Form("gausPion%i", i), "gaus", pionMassSqr * 0., pionMassSqr * 2.);
    gausKaon[i] = new TF1(Form("gausKaon%i", i), "gaus", kaonMassSqr * 0.5, kaonMassSqr * 1.5);
    gausProton[i] = new TF1(Form("gausProton%i", i), "gaus", protMassSqr * 0.5, protMassSqr * 1.5);

    gausTriplePion[i] = new TF1(Form("gausTriplePion%i", i), "gaus(0)+gaus(3)+gaus(6)", -0.5, 1.5);
    gausTripleKaon[i] = new TF1(Form("gausTripleKaon%i", i), "gaus(0)+gaus(3)+gaus(6)", -0.5, 1.5);
    gausTripleProton[i] = new TF1(Form("gausTripleProton%i", i), "gaus(0)+gaus(3)+gaus(6)", -0.5, 1.5);
  }

  int bin1 = 0, bin2 = 0;
  Double_t parPion[9], parKaon[9], parProton[9];

  for (int i = 0; i < NptBins; i++)
  {
    bin1 = hNsigMsqrPion[i]->GetXaxis()->FindBin(-3.); // 651;
    bin2 = hNsigMsqrPion[i]->GetXaxis()->FindBin(3.);  // 750;
    hMsqrPion[i] = (TH1D *)hNsigMsqrPion[i]->ProjectionY(Form("hMsqrPion%i", i), bin1, bin2);
    hMsqrKaon[i] = (TH1D *)hNsigMsqrKaon[i]->ProjectionY(Form("hMsqrKaon%i", i), bin1, bin2);
    hMsqrProton[i] = (TH1D *)hNsigMsqrProton[i]->ProjectionY(Form("hMsqrProton%i", i), bin1, bin2);
  }

  //Pion
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Fit(gausPion[i], "RQ0");
    hMsqrPion[i]->Fit(gausKaon[i], "RQ0");
    hMsqrPion[i]->Fit(gausProton[i], "RQ0");

    gausPion[i]->GetParameters(&parPion[0]);
    gausKaon[i]->GetParameters(&parPion[3]);
    gausProton[i]->GetParameters(&parPion[6]);

    gausTriplePion[i]->SetParameters(parPion);

    hMsqrPion[i]->Fit(gausTriplePion[i], "R");
    hMsqrSigmaPion->SetBinContent(i, gausTriplePion[i]->GetParameter(2));
    hMsqrSigmaPion->SetBinError(i, gausTriplePion[i]->GetParError(2));
  }

  //Kaon
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrKaon[i]->Fit(gausPion[i], "RQ0");
    hMsqrKaon[i]->Fit(gausKaon[i], "RQ0");
    hMsqrKaon[i]->Fit(gausProton[i], "RQ0");

    gausPion[i]->GetParameters(&parKaon[0]);
    gausKaon[i]->GetParameters(&parKaon[3]);
    gausProton[i]->GetParameters(&parKaon[6]);

    gausTripleKaon[i]->SetParameters(parKaon);

    hMsqrKaon[i]->Fit(gausTripleKaon[i], "R");
    hMsqrSigmaKaon->SetBinContent(i, gausTripleKaon[i]->GetParameter(5));
    hMsqrSigmaKaon->SetBinError(i, gausTripleKaon[i]->GetParError(5));
  }

  //Proton
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrProton[i]->Fit(gausPion[i], "RQ0");
    hMsqrProton[i]->Fit(gausKaon[i], "RQ0");
    hMsqrProton[i]->Fit(gausProton[i], "RQ0");

    gausPion[i]->GetParameters(&parProton[0]);
    gausKaon[i]->GetParameters(&parProton[3]);
    gausProton[i]->GetParameters(&parProton[6]);

    gausTripleProton[i]->SetParameters(parProton);

    hMsqrProton[i]->Fit(gausTripleProton[i], "R");
    hMsqrSigmaProton->SetBinContent(i, gausTripleProton[i]->GetParameter(8));
    hMsqrSigmaProton->SetBinError(i, gausTripleProton[i]->GetParError(8));
  }

  TF1 *polPion = new TF1("polPion", "pol4", 0.15, 3.7);
  TF1 *polKaon = new TF1("polKaon", "pol8", 0.15, 3.);
  TF1 *polProton = new TF1("polProton", "pol8", 0.15, 3.7);

  hMsqrSigmaPion->Fit(polPion, "R");
  hMsqrSigmaKaon->Fit(polKaon, "R");
  hMsqrSigmaProton->Fit(polProton, "R");

  TFile *fo = new TFile(outFileName, "recreate");
  fo->cd();
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Write();
    hMsqrKaon[i]->Write();
    hMsqrProton[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    gausTriplePion[i]->Write();
    gausTripleKaon[i]->Write();
    gausTripleProton[i]->Write();
  }
  hMsqrSigmaPion->Write();
  hMsqrSigmaKaon->Write();
  hMsqrSigmaProton->Write();

  polPion->Write();
  polKaon->Write();
  polProton->Write();
}

void FitCombgausMassSqr(const Char_t *inFileName, const Char_t *outFileName)
{
  TFile *fi = new TFile(inFileName, "read");

  TH2D *hNsigMsqrPion[NptBins];
  TH2D *hNsigMsqrKaon[NptBins];
  TH2D *hNsigMsqrProton[NptBins];

  TH1D *hMsqrPion[NptBins];
  TH1D *hMsqrKaon[NptBins];
  TH1D *hMsqrProton[NptBins];

  TF1 *gausPion[NptBins];
  TF1 *gausKaon[NptBins];
  TF1 *gausProton[NptBins];
  TF1 *gausMesonBg[NptBins];
  TF1 *gausBosonBg[NptBins];

  TF1 *gausTriplePion[NptBins];
  TF1 *gausTripleKaon[NptBins];
  TF1 *gausTripleProton[NptBins];

  TH1D *hMsqrSigmaPion = new TH1D("hMsqrSigmaPion", "hMsqrSigmaPion", NptBins, 0.15, 3.75);
  TH1D *hMsqrSigmaKaon = new TH1D("hMsqrSigmaKaon", "hMsqrSigmaKaon", NptBins, 0.15, 3.75);
  TH1D *hMsqrSigmaProton = new TH1D("hMsqrSigmaProton", "hMsqrSigmaProton", NptBins, 0.15, 3.75);

  for (int i = 0; i < NptBins; i++)
  {
    hNsigMsqrPion[i] = (TH2D *)fi->Get(Form("hNsigPionMSqr%i", i));
    hNsigMsqrKaon[i] = (TH2D *)fi->Get(Form("hNsigKaonMSqr%i", i));
    hNsigMsqrProton[i] = (TH2D *)fi->Get(Form("hNsigProtonMSqr%i", i));

    gausPion[i] = new TF1(Form("gausPion%i", i), "gaus", pionMassSqr * 0., pionMassSqr * 2.);
    gausKaon[i] = new TF1(Form("gausKaon%i", i), "gaus", kaonMassSqr * 0.5, kaonMassSqr * 1.5);
    gausProton[i] = new TF1(Form("gausProton%i", i), "gaus", protMassSqr * 0.5, protMassSqr * 1.5);
    gausMesonBg[i] = new TF1(Form("gausMesonBg%i", i), "gaus", pionMassSqr * 0., kaonMassSqr * 1.5);
    gausBosonBg[i] = new TF1(Form("gausBosonBg%i", i), "gaus", kaonMassSqr * 1.5, protMassSqr * 1.5);

    gausTriplePion[i] = new TF1(Form("gausTriplePion%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)", -0.5, 1.5);
    gausTripleKaon[i] = new TF1(Form("gausTripleKaon%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)", -0.5, 1.5);
    gausTripleProton[i] = new TF1(Form("gausTripleProton%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)", -0.5, 1.5);
  }

  int bin1 = 0, bin2 = 0;
  Double_t parPion[15], parKaon[15], parProton[12];

  for (int i = 0; i < NptBins; i++)
  {
    bin1 = hNsigMsqrPion[i]->GetXaxis()->FindBin(-3.); // 651;
    bin2 = hNsigMsqrPion[i]->GetXaxis()->FindBin(3.);  // 750;
    hMsqrPion[i] = (TH1D *)hNsigMsqrPion[i]->ProjectionY(Form("hMsqrPion%i", i), bin1, bin2);
    hMsqrKaon[i] = (TH1D *)hNsigMsqrKaon[i]->ProjectionY(Form("hMsqrKaon%i", i), bin1, bin2);
    hMsqrProton[i] = (TH1D *)hNsigMsqrProton[i]->ProjectionY(Form("hMsqrProton%i", i), bin1, bin2);
  }

  //Pion
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Fit(gausPion[i], "RQ0");
    hMsqrPion[i]->Fit(gausKaon[i], "RQ0");
    hMsqrPion[i]->Fit(gausProton[i], "RQ0");
    hMsqrPion[i]->Fit(gausMesonBg[i], "RQ0");
    hMsqrPion[i]->Fit(gausBosonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parPion[0]);
    gausKaon[i]->GetParameters(&parPion[3]);
    gausProton[i]->GetParameters(&parPion[6]);
    gausMesonBg[i]->GetParameters(&parPion[9]);
    gausBosonBg[i]->GetParameters(&parPion[12]);

    gausTriplePion[i]->SetParameters(parPion);

    hMsqrPion[i]->Fit(gausTriplePion[i], "R");
    hMsqrSigmaPion->SetBinContent(i, gausTriplePion[i]->GetParameter(2));
    hMsqrSigmaPion->SetBinError(i, gausTriplePion[i]->GetParError(2));
  }

  //Kaon
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrKaon[i]->Fit(gausPion[i], "RQ0");
    hMsqrKaon[i]->Fit(gausKaon[i], "RQ0");
    hMsqrKaon[i]->Fit(gausProton[i], "RQ0");
    hMsqrKaon[i]->Fit(gausMesonBg[i], "RQ0");
    hMsqrKaon[i]->Fit(gausBosonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parKaon[0]);
    gausKaon[i]->GetParameters(&parKaon[3]);
    gausProton[i]->GetParameters(&parKaon[6]);
    gausMesonBg[i]->GetParameters(&parKaon[9]);
    gausBosonBg[i]->GetParameters(&parKaon[12]);

    gausTripleKaon[i]->SetParameters(parKaon);

    hMsqrKaon[i]->Fit(gausTripleKaon[i], "R");
    hMsqrSigmaKaon->SetBinContent(i, gausTripleKaon[i]->GetParameter(5));
    hMsqrSigmaKaon->SetBinError(i, gausTripleKaon[i]->GetParError(5));
  }

  //Proton
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrProton[i]->Fit(gausPion[i], "RQ0");
    hMsqrProton[i]->Fit(gausKaon[i], "RQ0");
    hMsqrProton[i]->Fit(gausProton[i], "RQ0");
    hMsqrProton[i]->Fit(gausMesonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parProton[0]);
    gausKaon[i]->GetParameters(&parProton[3]);
    gausProton[i]->GetParameters(&parProton[6]);
    gausMesonBg[i]->GetParameters(&parProton[9]);

    gausTripleProton[i]->SetParameters(parProton);

    hMsqrProton[i]->Fit(gausTripleProton[i], "R");
    hMsqrSigmaProton->SetBinContent(i, gausTripleProton[i]->GetParameter(8));
    hMsqrSigmaProton->SetBinError(i, gausTripleProton[i]->GetParError(8));
  }

  TF1 *polPion = new TF1("polPion", "pol8", 0.15, 3.7);
  TF1 *polKaon = new TF1("polKaon", "pol8", 0.15, 3.);
  TF1 *polProton = new TF1("polProton", "pol8", 0.3, 3.7);

  hMsqrSigmaPion->Fit(polPion, "R");
  hMsqrSigmaKaon->Fit(polKaon, "R");
  hMsqrSigmaProton->Fit(polProton, "R");

  TFile *fo = new TFile(outFileName, "recreate");
  fo->cd();
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Write();
    hMsqrKaon[i]->Write();
    hMsqrProton[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    gausTriplePion[i]->Write();
    gausTripleKaon[i]->Write();
    gausTripleProton[i]->Write();
  }
  hMsqrSigmaPion->Write();
  hMsqrSigmaKaon->Write();
  hMsqrSigmaProton->Write();

  polPion->Write();
  polKaon->Write();
  polProton->Write();
}

void FitComb1gausMassSqr(const Char_t *inFileName, const Char_t *outFileName)
{
  TFile *fi = new TFile(inFileName, "read");

  TH2D *hNsigMsqrPion[NptBins];
  TH2D *hNsigMsqrKaon[NptBins];
  TH2D *hNsigMsqrProton[NptBins];

  TH1D *hMsqrPion[NptBins];
  TH1D *hMsqrKaon[NptBins];
  TH1D *hMsqrProton[NptBins];

  TF1 *gausPion[NptBins];
  TF1 *gausKaon[NptBins];
  TF1 *gausProton[NptBins];
  TF1 *gausMesonBg[NptBins];
  TF1 *gausBosonBg[NptBins];

  TF1 *gausTriplePion[NptBins];
  TF1 *gausTripleKaon[NptBins];
  TF1 *gausTripleProton[NptBins];

  TH1D *hMsqrSigmaPion = new TH1D("hMsqrSigmaPion", "hMsqrSigmaPion", NptBins, 0.15, 3.75);
  TH1D *hMsqrSigmaKaon = new TH1D("hMsqrSigmaKaon", "hMsqrSigmaKaon", NptBins, 0.15, 3.75);
  TH1D *hMsqrSigmaProton = new TH1D("hMsqrSigmaProton", "hMsqrSigmaProton", NptBins, 0.15, 3.75);

  for (int i = 0; i < NptBins; i++)
  {
    hNsigMsqrPion[i] = (TH2D *)fi->Get(Form("hNsigPionMSqr%i", i));
    hNsigMsqrKaon[i] = (TH2D *)fi->Get(Form("hNsigKaonMSqr%i", i));
    hNsigMsqrProton[i] = (TH2D *)fi->Get(Form("hNsigProtonMSqr%i", i));

    gausPion[i] = new TF1(Form("gausPion%i", i), "gaus", pionMassSqr * 0., pionMassSqr * 2.);
    gausKaon[i] = new TF1(Form("gausKaon%i", i), "gaus", kaonMassSqr * 0.5, kaonMassSqr * 1.5);
    gausProton[i] = new TF1(Form("gausProton%i", i), "gaus", protMassSqr * 0.5, protMassSqr * 1.5);
    gausMesonBg[i] = new TF1(Form("gausMesonBg%i", i), "gaus", pionMassSqr * 0., kaonMassSqr * 1.5);
    gausBosonBg[i] = new TF1(Form("gausBosonBg%i", i), "gaus", kaonMassSqr * 1.5, protMassSqr * 1.5);

    gausTriplePion[i] = new TF1(Form("gausTriplePion%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)", -0.5, 1.5);
    gausTripleKaon[i] = new TF1(Form("gausTripleKaon%i", i), "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)", -0.5, 1.5);
    gausTripleProton[i] = new TF1(Form("gausTripleProton%i", i), "gaus(0)+gaus(3)+gaus(6)", -0.5, 1.5);
  }

  int bin1 = 0, bin2 = 0;
  Double_t parPion[15], parKaon[15], parProton[9];

  for (int i = 0; i < NptBins; i++)
  {
    bin1 = hNsigMsqrPion[i]->GetXaxis()->FindBin(-3.); // 651;
    bin2 = hNsigMsqrPion[i]->GetXaxis()->FindBin(3.);  // 750;
    hMsqrPion[i] = (TH1D *)hNsigMsqrPion[i]->ProjectionY(Form("hMsqrPion%i", i), bin1, bin2);
    hMsqrKaon[i] = (TH1D *)hNsigMsqrKaon[i]->ProjectionY(Form("hMsqrKaon%i", i), bin1, bin2);
    hMsqrProton[i] = (TH1D *)hNsigMsqrProton[i]->ProjectionY(Form("hMsqrProton%i", i), bin1, bin2);
  }

  //Pion
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Fit(gausPion[i], "RQ0");
    hMsqrPion[i]->Fit(gausKaon[i], "RQ0");
    hMsqrPion[i]->Fit(gausProton[i], "RQ0");
    hMsqrPion[i]->Fit(gausMesonBg[i], "RQ0");
    hMsqrPion[i]->Fit(gausBosonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parPion[0]);
    gausKaon[i]->GetParameters(&parPion[3]);
    gausProton[i]->GetParameters(&parPion[6]);
    gausMesonBg[i]->GetParameters(&parPion[9]);
    gausBosonBg[i]->GetParameters(&parPion[12]);

    gausTriplePion[i]->SetParameters(parPion);

    hMsqrPion[i]->Fit(gausTriplePion[i], "R");
    hMsqrSigmaPion->SetBinContent(i, gausTriplePion[i]->GetParameter(2));
    hMsqrSigmaPion->SetBinError(i, gausTriplePion[i]->GetParError(2));
  }

  //Kaon
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrKaon[i]->Fit(gausPion[i], "RQ0");
    hMsqrKaon[i]->Fit(gausKaon[i], "RQ0");
    hMsqrKaon[i]->Fit(gausProton[i], "RQ0");
    hMsqrKaon[i]->Fit(gausMesonBg[i], "RQ0");
    hMsqrKaon[i]->Fit(gausBosonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parKaon[0]);
    gausKaon[i]->GetParameters(&parKaon[3]);
    gausProton[i]->GetParameters(&parKaon[6]);
    gausMesonBg[i]->GetParameters(&parKaon[9]);
    gausBosonBg[i]->GetParameters(&parKaon[12]);

    gausTripleKaon[i]->SetParameters(parKaon);

    hMsqrKaon[i]->Fit(gausTripleKaon[i], "R");
    hMsqrSigmaKaon->SetBinContent(i, gausTripleKaon[i]->GetParameter(5));
    hMsqrSigmaKaon->SetBinError(i, gausTripleKaon[i]->GetParError(5));
  }

  //Proton
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrProton[i]->Fit(gausPion[i], "RQ0");
    hMsqrProton[i]->Fit(gausKaon[i], "RQ0");
    hMsqrProton[i]->Fit(gausProton[i], "RQ0");
    hMsqrProton[i]->Fit(gausMesonBg[i], "RQ0");

    gausPion[i]->GetParameters(&parProton[0]);
    gausKaon[i]->GetParameters(&parProton[3]);
    gausProton[i]->GetParameters(&parProton[6]);

    gausTripleProton[i]->SetParameters(parProton);

    hMsqrProton[i]->Fit(gausTripleProton[i], "R");
    hMsqrSigmaProton->SetBinContent(i, gausTripleProton[i]->GetParameter(8));
    hMsqrSigmaProton->SetBinError(i, gausTripleProton[i]->GetParError(8));
  }

  TF1 *polPion = new TF1("polPion", "pol8", 0.15, 3.7);
  TF1 *polKaon = new TF1("polKaon", "pol8", 0.15, 3.);
  TF1 *polProton = new TF1("polProton", "pol8", 0.3, 3.7);

  hMsqrSigmaPion->Fit(polPion, "R");
  hMsqrSigmaKaon->Fit(polKaon, "R");
  hMsqrSigmaProton->Fit(polProton, "R");

  TFile *fo = new TFile(outFileName, "recreate");
  fo->cd();
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Write();
    hMsqrKaon[i]->Write();
    hMsqrProton[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    gausTriplePion[i]->Write();
    gausTripleKaon[i]->Write();
    gausTripleProton[i]->Write();
  }
  hMsqrSigmaPion->Write();
  hMsqrSigmaKaon->Write();
  hMsqrSigmaProton->Write();

  polPion->Write();
  polKaon->Write();
  polProton->Write();
}

void QA1DgausMassSqr(const Char_t *inFileName, const Char_t *outFileName)
{
  TFile *fi = new TFile(inFileName, "read");

  TH2D *hNsigMsqrPion[NptBins];
  TH2D *hNsigMsqrKaon[NptBins];
  TH2D *hNsigMsqrProton[NptBins];

  TH1D *hMsqrPion[NptBins];
  TH1D *hMsqrKaon[NptBins];
  TH1D *hMsqrProton[NptBins];

  TH1D *hMsqrPionAfter[NptBins];
  TH1D *hMsqrKaonAfter[NptBins];
  TH1D *hMsqrProtonAfter[NptBins];

  TF1 *gausPion[NptBins];
  TF1 *gausKaon[NptBins];
  TF1 *gausProton[NptBins];

  TF1 *gausTriplePion[NptBins];
  TF1 *gausTripleKaon[NptBins];
  TF1 *gausTripleProton[NptBins];

  TH1D *hMsqrSigmaPion = new TH1D("hMsqrSigmaPion", "hMsqrSigmaPion", NptBins, 0.15, 3.75);
  TH1D *hMsqrSigmaKaon = new TH1D("hMsqrSigmaKaon", "hMsqrSigmaKaon", NptBins, 0.15, 3.75);
  TH1D *hMsqrSigmaProton = new TH1D("hMsqrSigmaProton", "hMsqrSigmaProton", NptBins, 0.15, 3.75);

  for (int i = 0; i < NptBins; i++)
  {
    hNsigMsqrPion[i] = (TH2D *)fi->Get(Form("hNsigPionMSqr%i", i));
    hNsigMsqrKaon[i] = (TH2D *)fi->Get(Form("hNsigKaonMSqr%i", i));
    hNsigMsqrProton[i] = (TH2D *)fi->Get(Form("hNsigProtonMSqr%i", i));

    hMsqrPionAfter[i] = (TH1D *)fi->Get(Form("hNsigPionMSqrAfter%i", i));
    hMsqrKaonAfter[i] = (TH1D *)fi->Get(Form("hNsigKaonMSqrAfter%i", i));
    hMsqrProtonAfter[i] = (TH1D *)fi->Get(Form("hNsigProtonMSqrAfter%i", i));

    gausPion[i] = new TF1(Form("gausPion%i", i), "gaus", pionMassSqr * 0., pionMassSqr * 2.);
    gausKaon[i] = new TF1(Form("gausKaon%i", i), "gaus", kaonMassSqr * 0.5, kaonMassSqr * 1.5);
    gausProton[i] = new TF1(Form("gausProton%i", i), "gaus", protMassSqr * 0.5, protMassSqr * 1.5);

    gausTriplePion[i] = new TF1(Form("gausTriplePion%i", i), "gaus(0)+gaus(3)+gaus(6)", -0.5, 1.5);
    gausTripleKaon[i] = new TF1(Form("gausTripleKaon%i", i), "gaus(0)+gaus(3)+gaus(6)", -0.5, 1.5);
    gausTripleProton[i] = new TF1(Form("gausTripleProton%i", i), "gaus(0)+gaus(3)+gaus(6)", -0.5, 1.5);
  }

  int bin1 = 0, bin2 = 0;
  Double_t parPion[9], parKaon[9], parProton[9];

  for (int i = 0; i < NptBins; i++)
  {
    bin1 = 651; //hNsigMsqrPion[i]->FindBin(-2.5);
    bin2 = 750; //hNsigMsqrPion[i]->FindBin(2.5);
    hMsqrPion[i] = (TH1D *)hNsigMsqrPion[i]->ProjectionY(Form("hMsqrPion%i", i), bin1, bin2);
    hMsqrKaon[i] = (TH1D *)hNsigMsqrKaon[i]->ProjectionY(Form("hMsqrKaon%i", i), bin1, bin2);
    hMsqrProton[i] = (TH1D *)hNsigMsqrProton[i]->ProjectionY(Form("hMsqrProton%i", i), bin1, bin2);
  }

  //Pion
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Fit(gausPion[i], "RQ0");
    hMsqrPion[i]->Fit(gausKaon[i], "RQ0");
    hMsqrPion[i]->Fit(gausProton[i], "RQ0");

    gausPion[i]->GetParameters(&parPion[0]);
    gausKaon[i]->GetParameters(&parPion[3]);
    gausProton[i]->GetParameters(&parPion[6]);

    gausTriplePion[i]->SetParameters(parPion);

    hMsqrPion[i]->Fit(gausTriplePion[i], "R");
    hMsqrSigmaPion->SetBinContent(i, gausTriplePion[i]->GetParameter(2));
    hMsqrSigmaPion->SetBinError(i, gausTriplePion[i]->GetParError(2));
  }

  //Kaon
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrKaon[i]->Fit(gausPion[i], "RQ0");
    hMsqrKaon[i]->Fit(gausKaon[i], "RQ0");
    hMsqrKaon[i]->Fit(gausProton[i], "RQ0");

    gausPion[i]->GetParameters(&parKaon[0]);
    gausKaon[i]->GetParameters(&parKaon[3]);
    gausProton[i]->GetParameters(&parKaon[6]);

    gausTripleKaon[i]->SetParameters(parKaon);

    hMsqrKaon[i]->Fit(gausTripleKaon[i], "R");
    hMsqrSigmaKaon->SetBinContent(i, gausTripleKaon[i]->GetParameter(5));
    hMsqrSigmaKaon->SetBinError(i, gausTripleKaon[i]->GetParError(5));
  }

  //Proton
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrProton[i]->Fit(gausPion[i], "RQ0");
    hMsqrProton[i]->Fit(gausKaon[i], "RQ0");
    hMsqrProton[i]->Fit(gausProton[i], "RQ0");

    gausPion[i]->GetParameters(&parProton[0]);
    gausKaon[i]->GetParameters(&parProton[3]);
    gausProton[i]->GetParameters(&parProton[6]);

    gausTripleProton[i]->SetParameters(parProton);

    hMsqrProton[i]->Fit(gausTripleProton[i], "R");
    hMsqrSigmaProton->SetBinContent(i, gausTripleProton[i]->GetParameter(8));
    hMsqrSigmaProton->SetBinError(i, gausTripleProton[i]->GetParError(8));
  }

  TF1 *polPion = new TF1("polPion", "pol8", 0.15, 3.7);
  TF1 *polKaon = new TF1("polKaon", "pol8", 0.15, 3.);
  TF1 *polProton = new TF1("polProton", "pol8", 0.3, 3.7);

  hMsqrSigmaPion->Fit(polPion, "R");
  hMsqrSigmaKaon->Fit(polKaon, "R");
  hMsqrSigmaProton->Fit(polProton, "R");

  TFile *fo = new TFile(outFileName, "recreate");
  fo->cd();
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Write();
    hMsqrKaon[i]->Write();
    hMsqrProton[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPionAfter[i]->Write();
    hMsqrKaonAfter[i]->Write();
    hMsqrProtonAfter[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    gausTriplePion[i]->Write();
    gausTripleKaon[i]->Write();
    gausTripleProton[i]->Write();
  }
  hMsqrSigmaPion->Write();
  hMsqrSigmaKaon->Write();
  hMsqrSigmaProton->Write();

  polPion->Write();
  polKaon->Write();
  polProton->Write();
}

void Fit3studentMassSqr(const Char_t *inFileName, const Char_t *outFileName)
{
  const std::pair<Double_t, Double_t> PionRange[NptBins] = {{0.01, 0.03}, {0.00, 0.04}, {0.00, 0.04}, {-0.02, 0.06}, {-0.05, 0.08}, {-0.06, 0.1}, {-0.08, 0.1}, {-0.08, 0.1}, {-0.08, 0.1}, {-0.08, 0.1}, {-0.08, 0.1}, {-0.08, 0.1}, {-0.08, 0.1}, {-0.08, 0.1}, {-0.08, 0.1}};
  const std::pair<Double_t, Double_t> KaonRange[NptBins] = {{0.22, 0.27}, {0.22, 0.27}, {0.21, 0.28}, {0.20, 0.30}, {0.18, 0.30}, {0.17, 0.32}, {0.15, 0.35}, {0.15, 0.35}, {0.15, 0.35}, {0.15, 0.35}, {0.15, 0.35}, {0.15, 0.35}, {0.15, 0.35}, {0.15, 0.35}, {0.15, 0.35}};
  const std::pair<Double_t, Double_t> ProtonRange[NptBins] = {{0.80, 0.98}, {0.80, 0.98}, {0.80, 0.96}, {0.80, 0.96}, {0.80, 0.97}, {0.75, 1.05}, {0.75, 1.05}, {0.75, 1.05}, {0.75, 1.05}, {0.75, 1.05}, {0.75, 1.05}, {0.75, 1.05}, {0.75, 1.05}, {0.75, 1.05}, {0.75, 1.05}};

  TFile *fi = new TFile(inFileName, "read");

  TH2D *hNsigMsqrPion[NptBins];
  TH2D *hNsigMsqrKaon[NptBins];
  TH2D *hNsigMsqrProton[NptBins];

  TH1D *hMsqrPion[NptBins];
  TH1D *hMsqrKaon[NptBins];
  TH1D *hMsqrProton[NptBins];

  TF1 *gausPion[NptBins];
  TF1 *gausKaon[NptBins];
  TF1 *gausProton[NptBins];
  TF1 *funcPion[NptBins];
  TF1 *funcKaon[NptBins];
  TF1 *funcProton[NptBins];

  TF1 *gausTriplePion[NptBins];
  TF1 *gausTripleKaon[NptBins];
  TF1 *gausTripleProton[NptBins];
  TF1 *funcTriplePion[NptBins];
  TF1 *funcTripleKaon[NptBins];
  TF1 *funcTripleProton[NptBins];

  // TH1D *hMsqrSigmaPion = new TH1D("hMsqrSigmaPion", "hMsqrSigmaPion", NptBins, 0.15, 3.75);
  // TH1D *hMsqrSigmaKaon = new TH1D("hMsqrSigmaKaon", "hMsqrSigmaKaon", NptBins, 0.15, 3.75);
  // TH1D *hMsqrSigmaProton = new TH1D("hMsqrSigmaProton", "hMsqrSigmaProton", NptBins, 0.15, 3.75);

  std::vector<Double_t> vPionX, vPionY, vPionEX, vPionEY;
  std::vector<Double_t> vKaonX, vKaonY, vKaonEX, vKaonEY;
  std::vector<Double_t> vProtonX, vProtonY, vProtonEX, vProtonEY;

  for (int i = 0; i<NptBins;i++)
  {
    vPionX.push_back(ptBins[i]+0.5*(ptBins[i+1]-ptBins[i]));
    vKaonX.push_back(ptBins[i]+0.5*(ptBins[i+1]-ptBins[i]));
    vProtonX.push_back(ptBins[i]+0.5*(ptBins[i+1]-ptBins[i]));

    vPionEX.push_back(0.);
    vKaonEX.push_back(0.);
    vProtonEX.push_back(0.);
  }

  for (int i = 0; i < NptBins; i++)
  {
    hNsigMsqrPion[i] = (TH2D *)fi->Get(Form("hNsigPionMSqr%i", i));
    hNsigMsqrKaon[i] = (TH2D *)fi->Get(Form("hNsigKaonMSqr%i", i));
    hNsigMsqrProton[i] = (TH2D *)fi->Get(Form("hNsigProtonMSqr%i", i));

    gausPion[i] = new TF1(Form("gausPion%i", i), "gaus", PionRange[i].first, PionRange[i].second);
    gausKaon[i] = new TF1(Form("gausKaon%i", i), "gaus", KaonRange[i].first, KaonRange[i].second);
    gausProton[i] = new TF1(Form("gausProton%i", i), "gaus", ProtonRange[i].first, ProtonRange[i].second);

    funcPion[i] = new TF1(Form("funcPion%i", i), student, PionRange[i].first, PionRange[i].second, 4);
    funcKaon[i] = new TF1(Form("funcKaon%i", i), student, KaonRange[i].first, KaonRange[i].second, 4);
    funcProton[i] = new TF1(Form("funcProton%i", i), student, ProtonRange[i].first, ProtonRange[i].second, 4);

    gausTriplePion[i] = new TF1(Form("gausTriplePion%i", i), "gaus(0)+gaus(3)+gaus(6)", PionRange[i].first, ProtonRange[i].second);
    gausTripleKaon[i] = new TF1(Form("gausTripleKaon%i", i), "gaus(0)+gaus(3)+gaus(6)", PionRange[i].first, ProtonRange[i].second);
    gausTripleProton[i] = new TF1(Form("gausTripleProton%i", i), "gaus(0)+gaus(3)+gaus(6)", PionRange[i].first, ProtonRange[i].second);

    funcTriplePion[i] = new TF1(Form("funcTriplePion%i", i), student3, PionRange[i].first, ProtonRange[i].second, 12);
    funcTripleKaon[i] = new TF1(Form("funcTripleKaon%i", i), student3, PionRange[i].first, ProtonRange[i].second, 12);
    funcTripleProton[i] = new TF1(Form("funcTripleProton%i", i), student3, PionRange[i].first, ProtonRange[i].second, 12);
  }

  int bin1 = 0, bin2 = 0;
  Double_t parPion[12], parKaon[12], parProton[12];
  Double_t parGausPion[9], parGausKaon[9], parGausProton[9];

  for (int i = 0; i < NptBins; i++)
  {
    bin1 = hNsigMsqrPion[i]->GetXaxis()->FindBin(-3.); // 651;
    bin2 = hNsigMsqrPion[i]->GetXaxis()->FindBin(3.);  // 750;
    hMsqrPion[i] = (TH1D *)hNsigMsqrPion[i]->ProjectionY(Form("hMsqrPion%i", i), bin1, bin2);
    hMsqrKaon[i] = (TH1D *)hNsigMsqrKaon[i]->ProjectionY(Form("hMsqrKaon%i", i), bin1, bin2);
    hMsqrProton[i] = (TH1D *)hNsigMsqrProton[i]->ProjectionY(Form("hMsqrProton%i", i), bin1, bin2);
  }

  //Pion
  for (int i = 0; i < NptBins; i++)
  {
    if (i < 3)
    {
      for (int j=0;j<Niter;j++)
      {
        hMsqrPion[i]->Fit(gausPion[i], "RQ0");
        hMsqrPion[i]->Fit(gausKaon[i], "RQ0");
        hMsqrPion[i]->Fit(gausProton[i], "RQ0");
      }

      gausPion[i]->GetParameters(&parGausPion[0]);
      gausKaon[i]->GetParameters(&parGausPion[3]);
      gausProton[i]->GetParameters(&parGausPion[6]);

      gausTriplePion[i]->SetParameters(parGausPion);

      for (int j=0;j<Niter;j++)
      {
        if (j!= Niter-1) hMsqrPion[i]->Fit(gausTriplePion[i], "R0Q");
        if (j== Niter-1) hMsqrPion[i]->Fit(gausTriplePion[i], "R");
      }
      vPionY.push_back(gausTriplePion[i]->GetParameter(2));
      vPionEY.push_back(gausTriplePion[i]->GetParError(2));
      // hMsqrSigmaPion->SetBinContent(i + 1, gausTriplePion[i]->GetParameter(2));
      // hMsqrSigmaPion->SetBinError(i + 1, gausTriplePion[i]->GetParError(2));
    }
    if (i >= 3)
    {
      hMsqrPion[i]->Fit(gausPion[i], "RQ0");
      hMsqrPion[i]->Fit(gausKaon[i], "RQ0");
      hMsqrPion[i]->Fit(gausProton[i], "RQ0");

      funcPion[i]->SetParameter(0, gausPion[i]->GetParameter(0));
      funcPion[i]->SetParameter(1, 100.);
      funcPion[i]->SetParameter(2, gausPion[i]->GetParameter(2));
      funcPion[i]->SetParameter(3, gausPion[i]->GetParameter(1));
      funcKaon[i]->SetParameter(0, gausKaon[i]->GetParameter(0));
      funcKaon[i]->SetParameter(1, 100.);
      funcKaon[i]->SetParameter(2, gausKaon[i]->GetParameter(2));
      funcKaon[i]->SetParameter(3, gausKaon[i]->GetParameter(1));
      funcProton[i]->SetParameter(0, gausProton[i]->GetParameter(0));
      funcProton[i]->SetParameter(1, 100.);
      funcProton[i]->SetParameter(2, gausProton[i]->GetParameter(2));
      funcProton[i]->SetParameter(3, gausProton[i]->GetParameter(1));

      for (int j=0;j<Niter;j++)
      {
        hMsqrPion[i]->Fit(funcPion[i], "RQ0");
        hMsqrPion[i]->Fit(funcKaon[i], "RQ0");
        hMsqrPion[i]->Fit(funcProton[i], "RQ0");
      }

      funcPion[i]->GetParameters(&parPion[0]);
      funcKaon[i]->GetParameters(&parPion[4]);
      funcProton[i]->GetParameters(&parPion[8]);

      funcTriplePion[i]->SetParameters(parPion);

      for (int j=0;j<Niter;j++)
      {
        if (j!=Niter-1) hMsqrPion[i]->Fit(funcTriplePion[i], "RQ0");
        if (j==Niter-1) hMsqrPion[i]->Fit(funcTriplePion[i], "R");
      }
      vPionY.push_back(funcTriplePion[i]->GetParameter(2));
      vPionEY.push_back(funcTriplePion[i]->GetParError(2));
      // hMsqrSigmaPion->SetBinContent(i + 1, funcTriplePion[i]->GetParameter(2));
      // hMsqrSigmaPion->SetBinError(i + 1, funcTriplePion[i]->GetParError(2));
    }
  }

  //Kaon
  for (int i = 0; i < NptBins; i++)
  {
    if (i < 3)
    {
      for (int j=0;j<Niter;j++)
      {
        hMsqrKaon[i]->Fit(gausPion[i], "RQ0");
        hMsqrKaon[i]->Fit(gausKaon[i], "RQ0");
        hMsqrKaon[i]->Fit(gausProton[i], "RQ0");
      }

      gausPion[i]->GetParameters(&parGausKaon[0]);
      gausKaon[i]->GetParameters(&parGausKaon[3]);
      gausProton[i]->GetParameters(&parGausKaon[6]);

      gausTripleKaon[i]->SetParameters(parGausKaon);

      for (int j=0;j<Niter;j++)
      {
        hMsqrKaon[i]->Fit(gausTripleKaon[i], "R0Q");
      }
      vKaonY.push_back(gausTripleKaon[i]->GetParameter(5));
      vKaonEY.push_back(gausTripleKaon[i]->GetParError(5));
      // hMsqrSigmaKaon->SetBinContent(i + 1, gausTripleKaon[i]->GetParameter(5));
      // hMsqrSigmaKaon->SetBinError(i + 1, gausTripleKaon[i]->GetParError(5));
    }
    if (i >= 3)
    {
      hMsqrKaon[i]->Fit(gausPion[i], "RQ0");
      hMsqrKaon[i]->Fit(gausKaon[i], "RQ0");
      hMsqrKaon[i]->Fit(gausProton[i], "RQ0");

      funcPion[i]->SetParameter(0, gausPion[i]->GetParameter(0));
      funcPion[i]->SetParameter(1, 100.);
      funcPion[i]->SetParameter(2, gausPion[i]->GetParameter(2));
      funcPion[i]->SetParameter(3, gausPion[i]->GetParameter(1));
      funcKaon[i]->SetParameter(0, gausKaon[i]->GetParameter(0));
      funcKaon[i]->SetParameter(1, 100.);
      funcKaon[i]->SetParameter(2, gausKaon[i]->GetParameter(2));
      funcKaon[i]->SetParameter(3, gausKaon[i]->GetParameter(1));
      funcProton[i]->SetParameter(0, gausProton[i]->GetParameter(0));
      funcProton[i]->SetParameter(1, 100.);
      funcProton[i]->SetParameter(2, gausProton[i]->GetParameter(2));
      funcProton[i]->SetParameter(3, gausProton[i]->GetParameter(1));

      for (int j=0;j<Niter;j++)
      {
        hMsqrKaon[i]->Fit(funcPion[i], "RQ0");
        hMsqrKaon[i]->Fit(funcKaon[i], "RQ0");
        hMsqrKaon[i]->Fit(funcProton[i], "RQ0");
      }

      funcPion[i]->GetParameters(&parKaon[0]);
      funcKaon[i]->GetParameters(&parKaon[4]);
      funcProton[i]->GetParameters(&parKaon[8]);

      funcTripleKaon[i]->SetParameters(parKaon);

      for (int j=0;j<Niter;j++)
      {
        hMsqrKaon[i]->Fit(funcTripleKaon[i], "RQ0");
      }
      vKaonY.push_back(funcTripleKaon[i]->GetParameter(6));
      vKaonEY.push_back(funcTripleKaon[i]->GetParError(6));
      // hMsqrSigmaKaon->SetBinContent(i + 1, funcTripleKaon[i]->GetParameter(6));
      // hMsqrSigmaKaon->SetBinError(i + 1, funcTripleKaon[i]->GetParError(6));
    }
  }

  //Proton
  for (int i = 0; i < NptBins; i++)
  {
    if (i < 3)
    {
      for (int j=0;j<Niter;j++)
      {
        hMsqrProton[i]->Fit(gausPion[i], "RQ0");
        hMsqrProton[i]->Fit(gausKaon[i], "RQ0");
        hMsqrProton[i]->Fit(gausProton[i], "RQ0");
      }

      gausPion[i]->GetParameters(&parGausProton[0]);
      gausKaon[i]->GetParameters(&parGausProton[3]);
      gausProton[i]->GetParameters(&parGausProton[6]);

      gausTripleProton[i]->SetParameters(parGausProton);

      for (int j=0;j<Niter;j++)
      {
        hMsqrProton[i]->Fit(gausTripleProton[i], "RQ0");
      }
      vKaonY.push_back(gausTripleKaon[i]->GetParameter(8));
      vKaonEY.push_back(gausTripleKaon[i]->GetParError(8));
      // hMsqrSigmaProton->SetBinContent(i + 1, gausTripleProton[i]->GetParameter(8));
      // hMsqrSigmaProton->SetBinError(i + 1, gausTripleProton[i]->GetParError(8));
    }
    if (i >= 3)
    {
      hMsqrProton[i]->Fit(gausPion[i], "RQ0");
      hMsqrProton[i]->Fit(gausKaon[i], "RQ0");
      hMsqrProton[i]->Fit(gausProton[i], "RQ0");

      funcPion[i]->SetParameter(0, gausPion[i]->GetParameter(0));
      funcPion[i]->SetParameter(1, 100.);
      funcPion[i]->SetParameter(2, gausPion[i]->GetParameter(2));
      funcPion[i]->SetParameter(3, gausPion[i]->GetParameter(1));
      funcKaon[i]->SetParameter(0, gausKaon[i]->GetParameter(0));
      funcKaon[i]->SetParameter(1, 100.);
      funcKaon[i]->SetParameter(2, gausKaon[i]->GetParameter(2));
      funcKaon[i]->SetParameter(3, gausKaon[i]->GetParameter(1));
      funcProton[i]->SetParameter(0, gausProton[i]->GetParameter(0));
      funcProton[i]->SetParameter(1, 100.);
      funcProton[i]->SetParameter(2, gausProton[i]->GetParameter(2));
      funcProton[i]->SetParameter(3, gausProton[i]->GetParameter(1));

      for (int j=0;j<Niter;j++)
      {
        hMsqrProton[i]->Fit(funcPion[i], "RQ0");
        hMsqrProton[i]->Fit(funcKaon[i], "RQ0");
        hMsqrProton[i]->Fit(funcProton[i], "RQ0");
      }

      funcPion[i]->GetParameters(&parProton[0]);
      funcKaon[i]->GetParameters(&parProton[4]);
      funcProton[i]->GetParameters(&parProton[8]);

      funcTripleProton[i]->SetParameters(parProton);

      for (int j=0;j<Niter;j++)
      {
        hMsqrProton[i]->Fit(funcTripleProton[i], "RQ0");
      }
      vProtonY.push_back(funcTripleProton[i]->GetParameter(10));
      vProtonEY.push_back(funcTripleProton[i]->GetParError(10));
      // hMsqrSigmaProton->SetBinContent(i + 1, funcTripleProton[i]->GetParameter(10));
      // hMsqrSigmaProton->SetBinError(i + 1, funcTripleProton[i]->GetParError(10));
    }
  }

  TGraphErrors *hMsqrSigmaPion = new TGraphErrors(vPionX.size(),&vPionX[0],&vPionY[0],&vPionEX[0],&vPionEY[0]);
  TGraphErrors *hMsqrSigmaKaon = new TGraphErrors(vKaonX.size(),&vKaonX[0],&vKaonY[0],&vKaonEX[0],&vKaonEY[0]);
  TGraphErrors *hMsqrSigmaProton = new TGraphErrors(vProtonX.size(),&vProtonX[0],&vProtonY[0],&vProtonEX[0],&vProtonEY[0]);

  TF1 *polPion = new TF1("polPion", "pol3", 0.2, 3.2);
  TF1 *polKaon = new TF1("polKaon", "pol3", 0.2, 3.2);
  TF1 *polProton = new TF1("polProton", "pol3", 0.2, 3.2);

  // hMsqrSigmaPion->Fit(polPion, "R");
  // hMsqrSigmaKaon->Fit(polKaon, "R");
  // hMsqrSigmaProton->Fit(polProton, "R");

  TFile *fo = new TFile(outFileName, "recreate");
  fo->cd();
  for (int i = 0; i < NptBins; i++)
  {
    hMsqrPion[i]->Write();
    hMsqrKaon[i]->Write();
    hMsqrProton[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    if (i<3)
    {
      gausTriplePion[i]->Write();
      gausTripleKaon[i]->Write();
      gausTripleProton[i]->Write();
    }
    if (i>=3)
    {
      funcTriplePion[i]->Write();
      funcTripleKaon[i]->Write();
      funcTripleProton[i]->Write();
    }
  }
  hMsqrSigmaPion->Write();
  hMsqrSigmaKaon->Write();
  hMsqrSigmaProton->Write();

  // polPion->Write();
  // polKaon->Write();
  // polProton->Write();
}

void Fit2Dgaus(const Char_t *inFileName, const Char_t *outFileName)
{

  const double pionNsigMean[NptBins] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
  const double kaonNsigMean[NptBins] = {12., 6., 3., 1., -0.3, -1., -1.2, -1.5, -1.5, -1.5, -1.7, -2., -2., -2., -2.};
  const double protNsigMean[NptBins] = {22., 15., 9., 5., 3., 1.5, 0.2, -0.5, -1., -1.5, -1.8, -2., -2.2, -2.2, -2.5};

  const double pionNsigWidth[NptBins] = {4., 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5};
  const double kaonNsigWidth[NptBins] = {8., 5., 4., 4., 3.7, 3., 3., 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5};
  const double protNsigWidth[NptBins] = {8., 10., 7., 5., 4., 3.5, 3., 2.6, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5};

  const std::pair<double,double> pionMsqrRange [NptBins] = {{-0.039,0.078},{-0.039,0.078},{-0.039,0.078},{-0.039,0.078},
                                                            {-0.039,0.078},{-0.039,0.078},{-0.039,0.078},{-0.1,0.14},
                                                            {-0.039,0.078},{-0.039,0.078},{-0.039,0.078},{-0.036,0.019},
                                                            {-0.039,0.078},{-0.039,0.078},{-0.039,0.078}};

  const std::pair<double,double> kaonMsqrRange [NptBins] = {{0.098,0.38},{0.098,0.38},{0.098,0.38},{0.098,0.38},
                                                            {0.098,0.38},{0.098,0.38},{0.098,0.38},{0.13,0.4},
                                                            {0.098,0.38},{0.098,0.38},{0.098,0.38},{0.098,0.38},
                                                            {0.098,0.38},{0.098,0.38},{0.098,0.38}};

  const std::pair<double,double> protMsqrRange [NptBins] = {{0.35,1.4},{0.35,1.4},{0.35,1.4},{0.35,1.4},
                                                            {0.35,1.4},{0.35,1.4},{0.35,1.4},{0.35,1.4},
                                                            {0.35,1.4},{0.35,1.4},{0.35,1.4},{0.35,1.4},
                                                            {0.35,1.4},{0.35,1.4},{0.35,1.4}};

  const std::pair<double,double> tripleNsigRange [NptBins] = {{-6., 28.},{-6., 28.},{-6., 28.},{-6., 28.},
                                                              {-6., 28.},{-6., 28.},{-6., 28.},{-5., 4.},
                                                              {-5., 4.},{-6., 28.},{-6., 28.},{-5., 3.},
                                                              {-5., 3.},{-5., 3.},{-5., 3.}};

  const std::pair<double,double> tripleMsqrRange [NptBins] = {{-0.5,1.5},{-0.5,1.5},{-0.5,1.5},{-0.5,1.5},
                                                              {-0.5,1.5},{-0.5,1.5},{-0.5,1.5},{-0.4,1.1},
                                                              {-0.3,1.3},{-0.5,1.5},{-0.5,1.5},{-0.5,1.5},
                                                              {-0.5,1.5},{-0.5,1.5},{-0.5,1.5}};

  TFile *fi = new TFile(inFileName, "read");

  TH2F *hNsigMsqr[NptBins];

  for (int i = 0; i < NptBins; i++)
  {
    hNsigMsqr[i] = (TH2F *)fi->Get(Form("hNsigPionMSqr%i", i));
  }

  TF2 *gausPion[NptBins];
  TF2 *gausKaon[NptBins];
  TF2 *gausProton[NptBins];
  TF2 *gausPiBg[NptBins];
  TF2 *gausKaBg[NptBins];
  TF2 *gausProtBg[NptBins];
  TF2 *gausTriple[NptBins];

  TStopwatch timer;
  timer.Start();

  for (int i = 0; i < NptBins; i++)
  {
    gausPion[i] = new TF2(Form("gausPion%i", i), "xygaus", pionNsigMean[i] - pionNsigWidth[i], pionNsigMean[i] + pionNsigWidth[i], pionMsqrRange[i].first, pionMsqrRange[i].second);
    gausKaon[i] = new TF2(Form("gausKaon%i", i), "xygaus", kaonNsigMean[i] - kaonNsigWidth[i], kaonNsigMean[i] + kaonNsigWidth[i], kaonMsqrRange[i].first, kaonMsqrRange[i].second);
    gausProton[i] = new TF2(Form("gausProton%i", i), "xygaus", protNsigMean[i] - protNsigWidth[i], protNsigMean[i] + protNsigWidth[i], protMsqrRange[i].first, protMsqrRange[i].second);
    gausPiBg[i] = new TF2(Form("gausPibg%i", i), "xygaus", pionNsigMean[i] - 2*pionNsigWidth[i], pionNsigMean[i] + 2*pionNsigWidth[i], 0.5*pionMsqrRange[i].first, 2*pionMsqrRange[i].second);
    gausKaBg[i] = new TF2(Form("gausKabg%i", i), "xygaus", kaonNsigMean[i] - 2*kaonNsigWidth[i], kaonNsigMean[i] + 2*kaonNsigWidth[i], 0.5*kaonMsqrRange[i].first, 2*kaonMsqrRange[i].second);
    gausProtBg[i] = new TF2(Form("gausProtBg%i", i), "xygaus", protNsigMean[i] - 2*protNsigWidth[i], protNsigMean[i] + 2*protNsigWidth[i], 0.5*protMsqrRange[i].first, 2*protMsqrRange[i].second);
    gausTriple[i] = new TF2(Form("gausTriple%i", i), "xygaus(0)+xygaus(5)+xygaus(10)+xygaus(15)+xygaus(20)+xygaus(25)", tripleNsigRange[i].first, tripleNsigRange[i].second, tripleMsqrRange[i].first, tripleMsqrRange[i].second);
  }

  // gausPion[0] = new TF2(Form("gausPion%i",0),"xygaus",-3.,3.,0.,0.4);
  // gausKaon[0] = new TF2(Form("gausKaon%i",0),"xygaus",4.,20.,0.2,0.28);
  // gausProton[0] = new TF2(Form("gausProton%i",0),"xygaus",6.8,32.85,0.54,1.25);
  // gausProton[0]->SetParameters(7.97,20.85,0.888,0.04366);
  // gausTriple[0] = new TF2(Form("gausTriple%i",0),"xygaus(0)+xygaus(5)+xygaus(10)",-5.,32.85,-0.5,2.);

  int TestStop = 4;

  for (int i = 0; i < NptBins; i++)
  {
    std::cout << "Pions, Pt bin " << i << std::endl;
    for (int j=0;j<Niter;j++)
    {
      std::cout << "2+2D fit. Iter " << j << ": ";
      timer.Print();
      timer.Continue();
      hNsigMsqr[i]->Fit(gausPion[i], "RQ0");
      hNsigMsqr[i]->Fit(gausKaon[i], "RQ0");
      hNsigMsqr[i]->Fit(gausProton[i], "RQ0");
      hNsigMsqr[i]->Fit(gausPiBg[i], "RQ0");
      hNsigMsqr[i]->Fit(gausKaBg[i], "RQ0");
      hNsigMsqr[i]->Fit(gausProtBg[i], "RQ0");
    }

    Double_t par[30];
    gausPion[i]->GetParameters(&par[0]);
    gausKaon[i]->GetParameters(&par[5]);
    gausProton[i]->GetParameters(&par[10]);
    gausPiBg[i]->GetParameters(&par[15]);
    gausKaBg[i]->GetParameters(&par[20]);
    gausProtBg[i]->GetParameters(&par[25]);
    gausTriple[i]->SetParameters(par);

    for (int j=0;j<Niter;j++)
    {
      std::cout << "Triple 2+2D fit. Iter " << j << ": ";
      timer.Print();
      timer.Continue();
      if (j!=Niter-1) hNsigMsqr[i]->Fit(gausTriple[i], "RQ0");
      if (j==Niter-1) hNsigMsqr[i]->Fit(gausTriple[i], "R");
    }
  }

  TFile *fo = new TFile(outFileName, "recreate");
  fo->cd();
  for (int i = 0; i < NptBins; i++)
  {
    hNsigMsqr[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    gausPion[i]->Write();
    gausKaon[i]->Write();
    gausProton[i]->Write();
  }
  for (int i = 0; i < NptBins; i++)
  {
    gausTriple[i]->Write();
  }
}

void CalcPID(const Char_t *inFileName, const Char_t *outFileName, const Int_t mode)
{
  if (mode == 0)
  {
    FitCombgausMassSqr(inFileName, outFileName);
  }
  if (mode == 1)
  {
    FitComb1gausMassSqr(inFileName, outFileName);
  }
  if (mode == 3)
  {
    Fit3gausMassSqr(inFileName, outFileName);
  }
  if (mode == 4)
  {
    Fit4gausMassSqr(inFileName, outFileName);
  }
  if (mode == 5)
  {
    Fit5gausMassSqr(inFileName, outFileName);
  }
  if (mode == 6)
  {
    Fit6gausMassSqr(inFileName, outFileName);
  }
  if (mode == 13)
  {
    Fit3studentMassSqr(inFileName, outFileName);
  }
  if (mode == 26)
  {
    Fit2Dgaus(inFileName, outFileName);
  }
  if (mode == 100)
  {
    QA1DgausMassSqr(inFileName, outFileName);
  }
}