fitmultigamma_1p_cut.C

#include "Math/SpecFuncMathCore.h"
#include "TAxis.h"
#include "TDirectory.h"
#include "TF1.h"
#include "TGraph.h"
#include "TH1.h"
#include "TMath.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TTimeStamp.h"
Int_t idebug = 0;

// Total fitting function declaration
const int max_n = 126;
int n_steer[max_n];
//	0	1	2	3	4	5
enum
{
  kFlag,
  kMuPed,
  kMu,
  kMuN,
  kSigN,
  kMuDum,
  kSigDum,
  kADum,
  kAN
};

Int_t idump = 0;
Double_t binw = 1;

TF1 *fcor = 0;

Int_t threeint = 0;

std::map<int, TGraph> graphs;
int fill_gr = 0;

int npeak = 1; // Number of neutron peaks in the fit
TGraph *gr = 0;
Double_t fitfun(Double_t *x, Double_t *part)
{
  double xx = x[0] - part[kMuPed];
  if (xx <= 0)
    return 0;
  // Parameters are AMPLITUDES of the various peaks
  Double_t fitval = 0.;
  // If we consider pileup we need to update the amplitudes
  // otherwise they are ok
  // TF1 *fc=new
  // TF1("fc","[0]*exp(x*log([2])/[1]*[2])/ROOT::Math::tgamma(x/[1]*[2]+1)/(exp([2]*log([2]))/ROOT::Math::tgamma([2]+1))",0,100)
  double *par = part;
  double beta = part[kMuN] / part[kSigN] / part[kSigN];
  for (Int_t i = 0; i < npeak; i++)
  {
    if (par[kAN + i] > 0)
    {
      double ln_N = TMath::Log(part[kAN + i]);
      double mu = part[kMuN] * double(i + 1);
      double bm = beta * mu;
      double mypeak = 0;
      if (bm < 150)
      {
        double fval = ln_N - beta * xx + (bm - 1) * TMath::Log(xx) + bm * TMath::Log(beta);
        mypeak = TMath::Exp(fval) / TMath::Gamma(bm);
      }
      else
      {
        double fval = ln_N - beta * xx + (bm - 1) * TMath::Log(xx) + bm * TMath::Log(beta) - (bm * TMath::Log(bm) - bm + 0.5 * TMath::Log(2. * TMath::Pi() / bm));
        mypeak = TMath::Exp(fval);
      }
      fitval += mypeak;
      if (fill_gr)
      {
        graphs[i + 1].SetPoint(graphs[i + 1].GetN(), x[0], mypeak);
      }
    }
  }
  if (part[kMuDum] > 0)
  {
    double mybckg = part[kADum] * TMath::Exp(-0.5 * (x[0] - part[kMuDum]) * (x[0] - part[kMuDum]) / part[kSigDum] / part[kSigDum]);
    fitval += mybckg;
    if (fill_gr)
    {
      graphs[0].SetPoint(graphs[0].GetN(), x[0], mybckg);
    }
  }
  // printf("%g %g %g %g %g %g  %g->
  // ",par[0],par[1],par[2],par[3],par[4],par[5],par[6]); printf("x=%g
  // f=%g\n",x[0], fitval);
  if (gr)
    gr->SetPoint(gr->GetN(), x[0], fitval);
  return fitval;
}

TH1 *fit_h = nullptr;
int fit_status = 0;

void fitmultigamma_1p_cut(const char *hname, unsigned int iflag, Double_t mped = 0, Double_t m1 = 9, Double_t s1 = 2, Int_t npk = 10, double xfmin = 0, double xdummy = 0, int irun = 0, TString side = "", TString channel = "",  double mu = 0)
{
  npeak = npk;
  printf("fitmultigamma(hname=\"%s\", mped=%g, m1=%g, s1=%g, npeak=%d, mu=%g, RunNumber=%d\n", hname, mped, m1, s1, npeak, mu, irun);
  if (npeak > max_n)
  {
    printf("Please increase max_n=%d above npeak=%d\n", max_n, npeak);
    return;
  }
  TH1 *hn = new TH1F("h1", "h1", 200, 0, 200);
  TCanvas *c1 = new TCanvas("c1", "c1", 1);
  fit_h = nullptr;
  TH1F *h = (TH1F *)gDirectory->Get(hname);
  if (!h)
  {
    printf("Cannot find %s\n", hname);
    gDirectory->pwd();
    gDirectory->ls();
    return;
  }
  fit_h = h;
  if (gr)
    gr->Set(0);
  else
    gr = new TGraph();
    h->SetStats(0);
  // if(npeak>5)gStyle->SetOptFit(0);
  // else gStyle->SetOptFit(1);
  // gStyle->SetOptStat(11111111);
  ////gStyle->SetOptStat("ne");
  ////gStyle->SetOptFit(1);

  Int_t np = npeak + kAN; // Number of fit parameters

  TF1 *f = new TF1("f", fitfun, xfmin > 0 ? xfmin : h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax(), np);
  f->SetLineWidth(2);
  f->SetLineColor(2);
  f->FixParameter(kFlag, 0); // Funzione usata nel fit
  f->SetNpx(1000);
  f->SetParName(kMuPed, "#mu_{ped}");
  f->SetParName(kMuN, "#mu_{1n}");
  if (iflag & 0x1)
  {
    f->SetParameter(kMuPed, mped);
    f->SetParLimits(kMuPed, -5, 5);
  }
  else
  {
    f->FixParameter(kMuPed, mped);
    printf("muped = %g\n", mped);
  }
  f->SetParameter(kMuN, m1);
  printf("m1 = %g\n", m1);
  f->SetParameter(kSigN, s1);
  printf("s1 = %g\n", s1);
  f->SetParLimits(kSigN, 0.25 * s1, 4 * s1);
  f->SetParName(kSigN, "#sigma_{1n}");
  f->SetParName(kMu, "#mu");
  f->SetParLimits(kMu, 0, 1);
  if (mu <= 0)
    f->FixParameter(kMu, -mu);
  else
    f->FixParameter(kMu, mu);
  for (int ipeak = 0; ipeak < npeak; ipeak++)
  {

    f->SetParameter(kAN + ipeak, h->GetBinContent(h->GetXaxis()->FindBin(double(ipeak + 1) * m1 + mped)));
    f->SetParLimits(kAN + ipeak, 0, 100 * h->GetMaximum());
    //    f->SetParameter(kAN+ipeak, 1);
    f->SetParName(kAN + ipeak, TString::Format("A_{%d}", ipeak + 1));
    printf("%d %g ", ipeak, double(ipeak + 1) * m1 + mped);
    printf("%g\n", f->GetParameter(kAN + ipeak));
  }

  double xmin = 0, xmax = 90;
  for (int i = 0; i < 100; i++)
  {
    gr->SetPoint(gr->GetN(), i, f->Eval(i));
    //    hn->Fill(f->GetRandom());
  }
  gPad->SetLogy(1);
  h->Draw();
  //  f->Draw("same");
  if (xdummy != 0)
  {
    f->SetParName(kMuDum, "#mu bck");
    f->SetParName(kSigDum, "#sigma bck");
    f->SetParName(kADum, "A_{bck}");
    if (xdummy > 0)
    {
      f->SetParameter(kMuDum, xdummy);
      printf("m1 = %g\n", m1);
      f->SetParLimits(kMuDum, 0, m1 * 0.7);
    }
    else
    {
      f->FixParameter(kMuDum, -xdummy);
    }
    f->SetParameter(kSigDum, s1);
    f->SetParLimits(kSigDum, 10, 200);
    f->SetParameter(kADum, h->GetBinContent(h->GetXaxis()->FindBin(TMath::Abs(xdummy) * 1.5)));
    f->SetParLimits(kADum, 0, h->GetBinContent(h->GetXaxis()->FindBin(TMath::Abs(xdummy)) * 2));
  }
  else
  {
    f->FixParameter(kMuDum, 0);
    f->FixParameter(kSigDum, 0);
    f->FixParameter(kADum, 0);
  }
  h->Fit("f", "R");
  // double sigma = TMath::Sqrt(f->GetParameter(kMuN)/f->GetParameter(kSigN));
  // printf("sigma = %f\n",sigma);
  double sigma = f->GetParameter(kSigN);
  f->SetRange(TMath::Min(xfmin, f->GetParameter(kMuPed) + f->GetParameter(kMuN) - 2.5 * sigma), f->GetParameter(kMuPed) + f->GetParameter(kMuN) * npk);
  h->GetXaxis()->SetRangeUser(0, f->GetParameter(kMuPed) + f->GetParameter(kMuN) * npk);
  h->Fit("f", "R,L");
  fill_gr = 1;
  for (Int_t i = 0; i <= npeak; i++)
  {
    graphs[i].Set(0);
  }
  for (int ibin = 1; ibin < h->GetNbinsX(); ibin++)
  {
    double x = h->GetXaxis()->GetBinCenter(ibin);
    if (x >= f->GetXmin() && x <= f->GetXmax())
    {
      fitfun(&x, f->GetParameters());
    }
  }
  fill_gr = 0;
  for (Int_t i = 0; i < npeak; i++)
  {
    auto &gr = graphs[i];
    gr.SetName(TString::Format("%dn", i));
    gr.SetLineColor(i % 8 + 2);
    gr.SetLineWidth(1);
    if (gr.GetN() > 0)
    {
      gr.DrawClone("l");
      h->GetListOfFunctions()->Add(&gr);
    }
    else
    {
      printf("No data for gr %d\n", i);
    }
  }
  //    gr->Draw("al");







  TLine *line = new TLine(f->GetParameter(3), h->GetMinimum(), f->GetParameter(3), h->GetMaximum());
  line->SetLineColor(kBlack); // Set line color to red
  line->SetLineWidth(3);      // Set line width
  line->Draw("same");

  std::ostringstream RunNumber;
  RunNumber << "Run number: " << irun;
 
  TText *text0 = new TText(f->GetParameter(3) + 20, h->GetMaximum(), RunNumber.str().c_str());
  text0->SetTextColor(kBlack); // Set text color
  text0->SetTextSize(0.04);    // Set the text size
  text0->Draw();

  std::ostringstream position;
  position << "position  : " << f->GetParameter(3);

  std::ostringstream sigma2;
  sigma2 << "sigma     : " << f->GetParameter(4);

  std::ostringstream resolution;
  resolution << "resolution: " << (f->GetParameter(4) / f->GetParameter(3)) * 100 << " %";

  TText *text1 = new TText(f->GetParameter(3) + 20, h->GetMaximum() / 2, position.str().c_str());
  text1->SetTextColor(kBlack); // Set text color
  text1->SetTextSize(0.04);    // Set the text size
  text1->Draw();

  TText *text2 = new TText(f->GetParameter(3) + 20, h->GetMaximum() / 4, sigma2.str().c_str());
  text2->SetTextColor(kBlack); // Set text color
  text2->SetTextSize(0.04);    // Set the text size
  text2->Draw();

  TText *text3 = new TText(f->GetParameter(3) + 20, h->GetMaximum() / 8, resolution.str().c_str());
  text3->SetTextColor(kBlack); // Set text color
  text3->SetTextSize(0.04);    // Set the text size
  text3->Draw();

  h->GetListOfFunctions()->Add(line);
  h->GetListOfFunctions()->Add(text0);
  h->GetListOfFunctions()->Add(text1);
  h->GetListOfFunctions()->Add(text2);
  h->GetListOfFunctions()->Add(text3);
  TFile *f2 = new TFile(TString::Format("results/%d/1pcut/h_TDC_ZP%s_%s_A_H_cut_fit.root", irun, side.Data(), channel.Data()), "RECREATE");
  h->Write();














  c1->Update();
  c1->Show();
  printf("mu_1n===%f\n", f->GetParameter(kMuN));
  printf("rms_1n===%f\n", f->GetParameter(kSigN));
  printf("cnt_1n===%f\n", f->GetParameter(kAN));
  printf("cnt_2n===%f\n", f->GetParameter(kAN + 1));
  printf("cnt_3n===%f\n", f->GetParameter(kAN + 2));
  ////c1->Print("c1.pdf");
  c1->Print(TString::Format("results/%d/1pcut/Pietro_Peak1n_ZN%s_%s_%d.pdf",irun,side.Data(), channel.Data(),irun));
  //c1->Print(TString::Format("results/%d/1pcut/Peak1p_cut_ZN_%d.pdf",irun,irun));
  ////gSystem->Sleep(500);
}