HadronBgPrep Class Reference

Class to prepare sample for fitting in beta gamma bins. More...

#include <HadronBgPrep.h>

Public Member Functions
	HadronBgPrep ()
	Constructor: Sets the description, the properties and the parameters of the algorithm.
virtual	~HadronBgPrep ()
	Destructor.
	HadronBgPrep (int bgbins, double upperbg, double lowerbg, int cosbins, double uppercos, double lowercos, int injbins, double lowerinj, double upperinj, int nhitbins, double lowernhit, double uppernhit, double cut)
	Constructor: set the input variables.
void	prepareSample (std::shared_ptr< TTree > hadron, TFile *&outfile, const std::string &suffix, const std::string &bgcurvefile, const std::string &bgsigmafile, const std::string &pdg, bool ismakePlots)
	function to prepare sample for monitoring plots, bg curve fitting and sigma vs ionz fitting
void	defineHisto (std::vector< TH1F * > &htemp, const std::string &svar, const std::string &stype, const std::string &pdg)
	function to define histograms
void	plotDist (std::map< int, std::vector< TH1F * > > &hist, const std::string &suffix, int bins)
	function to plot the map of histograms
void	plotDist (std::vector< TH1F * > &hist, const std::string &suffix, int nbins)
	function to plot the histograms
void	setPars (TFile *&outfile, std::string pdg, std::vector< TH1F * > &hdedx_bg, std::vector< TH1F * > &hchi_bg, std::vector< TH1F * > &hionzsigma_bg, std::map< int, std::vector< TH1F * > > &hchi_inj)
	function to fill the parameters like mean and reso in the tree
void	fitGaussianWRange (TH1F *&temphist, gstatus &status, double sigmaR)
	function to perform gauss fit for input histogram
void	fit (TH1F *&hist, const std::string &pdg, gstatus &status)
	function to fit the histograms
void	printCanvasCos (std::map< int, std::vector< TH1F * > > &hchicos_allbg, std::map< int, std::vector< TH1F * > > &hchicos_1by3bg, std::map< int, std::vector< TH1F * > > &hchicos_2by3bg, std::map< int, std::vector< TH1F * > > &hchicos_3by3bg, const std::string &particle, const std::string &suffix)
	function to draw the dedx vs costh histograms
void	FormatGraph (TGraphErrors &gr, int flag, const std::string &name="")
	function to set graph cosmetics
void	clearVars ()
	function to clear the variables
void	deleteHistos (std::vector< TH1F * > &htemp)
	function to delete the histograms
double	getParticleMass (const std::string &particle)
	function to get the particle mass

Private Attributes
std::vector< double >	m_sumcos
	variables to add cos values
std::vector< double >	m_sumbg
	variables to add bg values
std::vector< double >	m_sumres_square
	variables to add square of resolution
std::vector< double >	m_suminj
	variables to add injection time
std::vector< double >	m_means
	mean variable
std::vector< double >	m_errors
	error variable
std::vector< int >	m_sumsize
	size of the bg bins
std::vector< int >	m_injsize
	size of the injection bins
int	m_bgBins
	bins for bg
double	m_bgMin
	min range of bg
double	m_bgMax
	max range of bg
int	m_injBins
	bins for injection time
double	m_injMin
	min range of injection time
double	m_injMax
	max range of injection time
int	m_cosBins
	bins for cosine
double	m_cosMin
	min range of cosine
double	m_cosMax
	max range of cosine
int	m_nhitBins
	bins for nhits
double	m_nhitMin
	min range of nhits
double	m_nhitMax
	max range of nhits
double	m_cut
	cut to clean protons

Detailed Description

Class to prepare sample for fitting in beta gamma bins.

Definition at line 43 of file HadronBgPrep.h.

Constructor & Destructor Documentation

HadronBgPrep() [1/2]

HadronBgPrep

(

)

Constructor: Sets the description, the properties and the parameters of the algorithm.

Definition at line 12 of file HadronBgPrep.cc.

{
  m_bgBins = 15;
  m_bgMin = 0.1;
  m_bgMax = 40;
  m_cosBins = 18;
  m_cosMin = -1.0;
  m_cosMax = 1.0;
  m_injBins = 20;
  m_injMin = 0;
  m_injMax = 80000;
  m_nhitBins = 10;
  m_nhitMin = 7;
  m_nhitMax = 39;
  m_cut = 0.5;
}

~HadronBgPrep()

virtual ~HadronBgPrep ( )

inlinevirtual

Destructor.

Definition at line 55 of file HadronBgPrep.h.

{};

HadronBgPrep() [2/2]

HadronBgPrep	(	int	bgbins,
		double	upperbg,
		double	lowerbg,
		int	cosbins,
		double	uppercos,
		double	lowercos,
		int	injbins,
		double	lowerinj,
		double	upperinj,
		int	nhitbins,
		double	lowernhit,
		double	uppernhit,
		double	cut )

Constructor: set the input variables.

Definition at line 28 of file HadronBgPrep.cc.

{
  m_bgBins = bgbins;
  m_bgMin = lowerbg;
  m_bgMax = upperbg;
  m_cosBins = cosbins;
  m_cosMin = lowercos;
  m_cosMax = uppercos;
  m_injBins = injbins;
  m_injMin = lowerinj;
  m_injMax = upperinj;
  m_nhitBins = nhitbins;
  m_nhitMin = lowernhit;
  m_nhitMax = uppernhit;
  m_cut = cut;
}

Member Function Documentation

clearVars()

void clearVars ( )

inline

function to clear the variables

Definition at line 138 of file HadronBgPrep.h.

    {
      m_sumcos.clear();
      m_sumbg.clear();
      m_sumres_square.clear();
      m_sumsize.clear();
      m_means.clear();
      m_errors.clear();
      m_suminj.clear();
    }

defineHisto()

void defineHisto	(	std::vector< TH1F * > &	htemp,
		const std::string &	svar,
		const std::string &	stype,
		const std::string &	pdg )

function to define histograms

Definition at line 249 of file HadronBgPrep.cc.

{
  int nbdEdx = 200, nbchi = 200;
  double dedxMax = 20.0;
 
  if (pdg == "pion") {
    nbdEdx = 400, dedxMax = 4.0;
  } else if (pdg == "kaon") {
    nbdEdx = 500, dedxMax = 5.0;
  } else if (pdg == "proton") {
    nbdEdx = 1500, dedxMax = 30.0;
    nbchi = 350;
  } else if (pdg == "muon") {
    nbdEdx = 300, dedxMax = 3.0;
  } else if (pdg == "electron") {
    nbdEdx = 200, dedxMax = 2.0;
  }
 
  int bins;
  double min, max;
 
  if (stype == "bg") bins = m_bgBins, min = m_bgMin, max = m_bgMax ;
  else if (stype == "inj_0" || stype == "inj_1") bins = m_injBins, min = m_injMin, max = m_injMax ;
  else if (stype == "nhit") bins = m_nhitBins, min = m_nhitMin, max = m_nhitMax ;
  else bins = m_cosBins, min = m_cosMin, max = m_cosMax ;
 
  double step = (max - min) / bins;
  if (stype == "nhit") step = (max - min + 1) / bins;
 
  for (int j = 0; j < bins; ++j) {
 
    double start = min + j * step;
    double end = start + step;
    if (stype == "nhit") end = int((start + step) * 0.99999);
 
    // initialize the histograms
    std::string histname = Form("%s_%s_%s_%d", pdg.data(), svar.data(), stype.data(), j);
    std::string title = Form("%s_%s_%s (%.02f, %.02f)", pdg.data(), svar.data(), stype.data(), start, end);
 
    if (svar == "dedx")
      htemp.push_back(new TH1F(histname.data(), title.data(), nbdEdx, 0, dedxMax));
    else if (svar == "chi")
      htemp.push_back(new TH1F(histname.data(), title.data(), nbchi, -10.0, 10.0));
    else
      htemp.push_back(new TH1F(histname.data(), title.data(), 300, -3.0, 3.0));
 
    htemp[j]->GetXaxis()->SetTitle(Form("%s", svar.data()));
    htemp[j]->GetYaxis()->SetTitle("Entries");
  }
}

deleteHistos()

void deleteHistos ( std::vector< TH1F * > & htemp )

inline

function to delete the histograms

Definition at line 152 of file HadronBgPrep.h.

    {
      for (auto& hist : htemp) hist->Delete();
    }

fit()

void fit	(	TH1F *&	hist,
		const std::string &	pdg,
		gstatus &	status )

function to fit the histograms

Definition at line 501 of file HadronBgPrep.cc.

{
 
  if (pdg == "pion") fitGaussianWRange(hist, status, 1.0);
  else fitGaussianWRange(hist, status, 2.0);
 
  hist->SetFillColorAlpha(kAzure + 1, 0.30);
 
  if (status == OK) {
    double mean = hist->GetFunction("gaus")->GetParameter(1);
    double meanerr = hist->GetFunction("gaus")->GetParError(1);
    double width = hist->GetFunction("gaus")->GetParameter(2);
 
    std::string title = Form("#mu_{fit}: %0.03f #pm %0.03f, #sigma_{fit}: %0.03f", mean, meanerr, width);
    hist->SetTitle(Form("%s, %s", hist->GetTitle(), title.data()));
  }
 
}

fitGaussianWRange()

void fitGaussianWRange	(	TH1F *&	temphist,
		gstatus &	status,
		double	sigmaR )

function to perform gauss fit for input histogram

Definition at line 521 of file HadronBgPrep.cc.

{
  double histmean = temphist->GetMean();
  double histrms = temphist->GetRMS();
  temphist->GetXaxis()->SetRangeUser(histmean - 5.0 * histrms, histmean + 5.0 * histrms);
 
  int fs = temphist->Fit("gaus", "Q");
  if (fs != 0) {
    B2INFO(Form("\tFit (round 1) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
    status = Failed;
    return;
  } else {
    double mean = temphist->GetFunction("gaus")->GetParameter(1);
    double width = temphist->GetFunction("gaus")->GetParameter(2);
    temphist->GetXaxis()->SetRangeUser(mean - 5.0 * width, mean + 5.0 * width);
    fs = temphist->Fit("gaus", "QR", "", mean - sigmaR * width, mean + sigmaR * width);
    if (fs != 0) {
      B2INFO(Form("\tFit (round 2) for hist (%s) failed (status = %d)", temphist->GetName(), fs));
      status = Failed;
      return;
    } else {
      temphist->GetXaxis()->SetRangeUser(mean - 5.0 * width, mean + 5.0 * width);
      B2INFO(Form("\tFit for hist (%s) successful (status = %d)", temphist->GetName(), fs));
      status = OK;
    }
  }
}

FormatGraph()

void FormatGraph ( TGraphErrors & gr, int flag, const std::string & name = "" )

inline

function to set graph cosmetics

Definition at line 110 of file HadronBgPrep.h.

    {
      if (flag == 0) {
        gr.SetTitle(Form("%s;cos(#theta);#chi_{mean}", name.data()));
        gr.SetMarkerStyle(24);
        gr.SetMarkerColor(2);
        gr.SetMarkerSize(.7);
        gr.SetMaximum(1.0);
        gr.SetMinimum(-1.0);
        gr.GetXaxis()->SetRangeUser(-1, 1);
      } else if (flag == 1) {
        gr.SetMarkerStyle(22);
        gr.SetMarkerColor(9);
        gr.SetMarkerSize(.7);
      } else if (flag == 2) {
        gr.SetTitle(Form("%s;cos(#theta);#sigma", name.data()));
        gr.SetMarkerStyle(24);
        gr.SetMarkerColor(2);
        gr.SetMarkerSize(.7);
        gr.SetMaximum(2);
        gr.SetMinimum(0);
        gr.GetXaxis()->SetRangeUser(-1, 1);
      }
    }

getParticleMass()

double getParticleMass ( const std::string & particle )

inline

function to get the particle mass

Definition at line 160 of file HadronBgPrep.h.

    {
      double mass;
      if (particle == "pion") mass = Const::pion.getMass();
      else if (particle == "kaon") mass = Const::kaon.getMass();
      else if (particle == "proton") mass = Const::proton.getMass();
      else if (particle == "muon") mass = Const::muon.getMass();
      else if (particle == "electron") mass = Const::electron.getMass();
      else mass = 0.0;
      return mass;
    }

plotDist() [1/2]

void plotDist	(	std::map< int, std::vector< TH1F * > > &	hist,
		const std::string &	suffix,
		int	bins )

function to plot the map of histograms

Definition at line 301 of file HadronBgPrep.cc.

{
 
  TCanvas ctmp(Form("cdcdedx_%s", sname.data()), "", 1200, 600);
  ctmp.Divide(2, 1);
  ctmp.SetBatch(kTRUE);
 
  std::stringstream psname;
  psname << Form("plots/HadronPrep/%s.pdf[", sname.data());
  ctmp.Print(psname.str().c_str());
  psname.str("");
  psname << Form("plots/HadronPrep/%s.pdf", sname.data());
  for (int j = 0; j < bins; ++j) {
 
    for (int i = 0 ; i < 2; ++i) {
      ctmp.cd(i + 1);
      hist[i][j]->SetFillColorAlpha(i + 5, 0.25);
      hist[i][j]->Draw();
    }
    ctmp.Print(psname.str().c_str());
  }
  psname.str("");
  psname << Form("plots/HadronPrep/%s.pdf]", sname.data());
  ctmp.Print(psname.str().c_str());
}

plotDist() [2/2]

void plotDist	(	std::vector< TH1F * > &	hist,
		const std::string &	suffix,
		int	nbins )

function to plot the histograms

Definition at line 328 of file HadronBgPrep.cc.

{
 
  TCanvas ctmp(Form("cdcdedx_%s", sname.data()), "", 1200, 1200);
  ctmp.Divide(2, 2);
  ctmp.SetBatch(kTRUE);
 
  std::stringstream psname;
  psname << Form("plots/HadronPrep/%s.pdf[", sname.data());
  ctmp.Print(psname.str().c_str());
  psname.str("");
  psname << Form("plots/HadronPrep/%s.pdf", sname.data());
 
  for (int i = 0 ; i < nbins; ++i) {
    ctmp.cd(i % 4 + 1);
    hist[i]->SetFillColor(kYellow - 9);
    hist[i]->Draw();
 
    if ((i + 1) % 4 == 0 || (i + 1) == nbins) {
      ctmp.Print(psname.str().c_str());
      ctmp.Clear("D");
    }
  }
  psname.str("");
  psname << Form("plots/HadronPrep/%s.pdf]", sname.data());
  ctmp.Print(psname.str().c_str());
}

prepareSample()

void prepareSample	(	std::shared_ptr< TTree >	hadron,
		TFile *&	outfile,
		const std::string &	suffix,
		const std::string &	bgcurvefile,
		const std::string &	bgsigmafile,
		const std::string &	pdg,
		bool	ismakePlots )

function to prepare sample for monitoring plots, bg curve fitting and sigma vs ionz fitting

Definition at line 46 of file HadronBgPrep.cc.

{
 
  std::vector<TH1F*>  hdedx_bg, hchi_bg, hionzsigma_bg;
  std::map<int, std::vector<TH1F*>>  hchi_inj, hchicos_allbg, hchicos_1by3bg, hchicos_2by3bg, hchicos_3by3bg;
 
  //define histograms
  defineHisto(hdedx_bg, "dedx", "bg", pdg.data());
  defineHisto(hchi_bg, "chi", "bg", pdg.data());
  defineHisto(hionzsigma_bg, "ionzsigma", "bg", pdg.data());
  for (int i = 0; i < 2; ++i) {
 
    defineHisto(hchi_inj[i], "chi", Form("inj_%d", i), pdg.data());
 
    std::string charge = "pos";
    if (i == 1) charge = "neg";
 
    defineHisto(hchicos_allbg[i], "chi", Form("%s_allbg_cos", charge.data()), pdg.data());
    defineHisto(hchicos_1by3bg[i], "chi", Form("%s_1b3bg_cos", charge.data()), pdg.data());
    defineHisto(hchicos_2by3bg[i], "chi", Form("%s_2b3bg_cos", charge.data()), pdg.data());
    defineHisto(hchicos_3by3bg[i], "chi", Form("%s_3b3bg_cos", charge.data()), pdg.data());
  }
 
 
  // Create some containers to calculate averages
  for (int i = 0; i < m_bgBins; ++i) {
    m_sumcos.push_back(0.);
    m_sumbg.push_back(0.);
    m_sumres_square.push_back(0.);
    m_sumsize.push_back(0.);
    m_means.push_back(0.);
    m_errors.push_back(0.);
  }
 
  for (int i = 0; i < m_injBins; ++i) {
    m_injsize.push_back(0.);
    m_suminj.push_back(0.);
  }
 
  // --------------------------------------------------
  // LOOP OVER EVENTS AND FILL CONTAINERS
  // --------------------------------------------------
  // Fill the histograms to be fitted
 
  double dedxnosat;    // dE/dx w/o HS correction (use to get new parameters)
  double p;       // track momentum
  double costh;   // cosine of track polar angle
  double timereso;
  int nhit;       // number of hits on this track
  double injtime; //injection time
  double isher;
  int charge;
 
  hadron->SetBranchAddress("dedxnosat", &dedxnosat); //must be without any HS correction
  hadron->SetBranchAddress("p", &p);
  hadron->SetBranchAddress("costh", &costh);
  hadron->SetBranchAddress("timereso", &timereso);
  hadron->SetBranchAddress("nhits", &nhit);
  hadron->SetBranchAddress("injtime", &injtime);
  hadron->SetBranchAddress("injring", &isher);
  hadron->SetBranchAddress("charge", &charge);
 
  double mass = getParticleMass(pdg);
  if (mass == 0.0) B2FATAL("Mass of particle " << pdg.data() << " is zero");
 
  const double cosstep = (m_cosMax - m_cosMin) / m_cosBins;
  const double tstep = (m_injMax - m_injMin) / m_injBins;
 
  CDCDedxMeanPred mbg;
  mbg.setParameters(bgcurvefile);
  CDCDedxSigmaPred sbg;
  sbg.setParameters(bgsigmafile);
 
  CDCDedxHadSat had;
  had.setParameters("sat-pars.fit.txt");
 
  unsigned int entries = hadron->GetEntries();
 
  for (unsigned int index = 0; index < entries; ++index) {
 
    hadron->GetEntry(index);
 
    int chg = (charge < 0) ? 1 : 0;
    double bg = fabs(p) / mass;
 
    // clean up bad events and restrict the momentum range
    if (nhit <  0 || nhit  > 100)continue;
    if (injtime <= 0 || isher < 0)continue;
 
    if (dedxnosat <= 0)continue;
    if (costh != costh)continue;
    if (bg < m_bgMin || bg > m_bgMax)continue;
 
    if (fabs(p) > 7.0)continue;
 
    if (pdg == "proton")  if ((dedxnosat - 0.45) * abs(p) * abs(p) < m_cut)continue;
 
    int bgBin = static_cast<int>((bg - m_bgMin) / (m_bgMax - m_bgMin) * m_bgBins);
    bgBin = std::min(m_bgBins - 1, std::max(0, bgBin));
 
    if (bgBin >= m_bgBins) {
      B2WARNING("bgBin out of range: " << bgBin
                << " (valid range: 0-" << (m_bgBins - 1) << ")");
    }
 
    double dedx_new = had.D2I(costh, had.I2D(costh, 1.0) * dedxnosat);
 
    hdedx_bg[bgBin]->Fill(dedx_new);
 
    //get predicted value of dEdx mean
    double dedx_cur = mbg.getMean(bg);
    double dedx_res = sbg.getSigma(dedx_new, nhit, costh, timereso);
 
    if (dedx_res == 0) {
      B2INFO("RESOLUTION IS ZERO!!!");
      continue;
    }
 
    //Modified chi_new values
    double chi_new  = (dedx_new - dedx_cur) / dedx_res;
 
    hchi_bg[bgBin]->Fill(chi_new);
 
    double ionz_res = sbg.cosPrediction(costh) * sbg.nhitPrediction(nhit) * timereso;
 
    if (ionz_res == 0) {
      B2INFO("RESOLUTION costh * nhit * timereso IS ZERO!!!");
      continue;
    }
 
    hionzsigma_bg[bgBin]->Fill((dedx_new - dedx_cur) / ionz_res);
 
    m_sumcos[bgBin] += costh;
    m_sumbg[bgBin] += bg;
    m_sumres_square[bgBin] += pow(dedx_res, 2);
    m_sumsize[bgBin] += 1;
 
    // make histograms of dE/dx vs. cos(theta) for validation
    int icos = static_cast<int>((costh + 1) / cosstep);
    if (icos >= m_cosBins) {
      B2WARNING("cosBin (icos) out of range: " << icos
                << " (valid range: 0-" << (m_cosBins - 1) << ")");
    }
    icos = std::min(m_cosBins - 1, icos);
 
    hchicos_allbg[chg][icos]->Fill(chi_new);
 
    if (bgBin <= int(m_bgBins / 3)) hchicos_1by3bg[chg][icos]->Fill(chi_new);
    else if (bgBin <= int(2 * m_bgBins / 3)) hchicos_2by3bg[chg][icos]->Fill(chi_new);
    else hchicos_3by3bg[chg][icos]->Fill(chi_new);
 
    if (injtime > m_injMax) injtime = m_injMax - 10.0;
    int wr = 0;
    if (isher > 0.5) wr = 1;
    int injBin = static_cast<int>((injtime - m_injMin) / tstep);
    if (injBin >= m_injBins) {
      B2WARNING("injBin out of range: " << injBin
                << " (valid range: 0-" << (m_injBins - 1) << ")");
    }
    injBin = std::min(m_injBins - 1, std::max(0, injBin));
    hchi_inj[wr][injBin]->Fill(chi_new);
 
    m_suminj[injBin] += injtime;
    m_injsize[injBin] += 1;
 
  } // end of event loop
 
  // --------------------------------------------------
  // FIT IN BINS OF BETA-GAMMA AND COS(THETA)
  // --------------------------------------------------
  // fit the histograms with Gaussian functions
  // and extract the means and errors
 
  setPars(outfile, pdg, hdedx_bg, hchi_bg, hionzsigma_bg, hchi_inj);
 
  // Plot the histograms
  if (ismakePlots) {
 
    plotDist(hdedx_bg, Form("fits_dedx_inbg_%s_%s", suffix.data(), pdg.data()), m_bgBins);
    plotDist(hchi_bg, Form("fits_chi_inbg_%s_%s", suffix.data(), pdg.data()), m_bgBins);
    plotDist(hionzsigma_bg, Form("fits_ionzreso_inbg_%s_%s", suffix.data(), pdg.data()), m_bgBins);
    plotDist(hchi_inj, Form("fits_chi_ininj_%s_%s", suffix.data(), pdg.data()), m_injBins);
    printCanvasCos(hchicos_allbg, hchicos_1by3bg, hchicos_2by3bg, hchicos_3by3bg, pdg, suffix);
  }
 
  // delete histograms
  clearVars();
  deleteHistos(hdedx_bg);
  deleteHistos(hchi_bg);
  deleteHistos(hionzsigma_bg);
  for (int i = 0; i < 2; ++i) {
    deleteHistos(hchi_inj[i]);
    deleteHistos(hchicos_allbg[i]);
    deleteHistos(hchicos_1by3bg[i]);
    deleteHistos(hchicos_2by3bg[i]);
    deleteHistos(hchicos_3by3bg[i]);
  }
 
}

printCanvasCos()

void printCanvasCos	(	std::map< int, std::vector< TH1F * > > &	hchicos_allbg,
		std::map< int, std::vector< TH1F * > > &	hchicos_1by3bg,
		std::map< int, std::vector< TH1F * > > &	hchicos_2by3bg,
		std::map< int, std::vector< TH1F * > > &	hchicos_3by3bg,
		const std::string &	particle,
		const std::string &	suffix )

function to draw the dedx vs costh histograms

Definition at line 553 of file HadronBgPrep.cc.

{
 
  std::vector<std::vector<double>> chicos(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> sigmacos(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> chicoserr(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> sigmacoserr(2, std::vector<double>(m_cosBins));
 
  std::vector<std::vector<double>> chicos_1b3bg(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> sigmacos_1b3bg(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> chicos_1b3bgerr(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> sigmacos_1b3bgerr(2, std::vector<double>(m_cosBins));
 
  std::vector<std::vector<double>> chicos_2b3bg(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> sigmacos_2b3bg(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> chicos_2b3bgerr(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> sigmacos_2b3bgerr(2, std::vector<double>(m_cosBins));
 
  std::vector<std::vector<double>> chicos2(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> sigmacos2(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> chicos2err(2, std::vector<double>(m_cosBins));
  std::vector<std::vector<double>> sigmacos2err(2, std::vector<double>(m_cosBins));
 
  for (int c = 0; c < 2; ++c) {
    for (int i = 0; i < m_cosBins; ++i) {
      if (hchicos_allbg[c][i]->Integral() > 100) {
        gstatus allbgstat;
        fit(hchicos_allbg[c][i],  pdg.data(), allbgstat);
        if (allbgstat == OK) {
          chicos[c][i] = hchicos_allbg[c][i]->GetFunction("gaus")->GetParameter(1);
          chicoserr[c][i] = hchicos_allbg[c][i]->GetFunction("gaus")->GetParError(1);
          sigmacos[c][i] = hchicos_allbg[c][i]->GetFunction("gaus")->GetParameter(2);
          sigmacoserr[c][i] = hchicos_allbg[c][i]->GetFunction("gaus")->GetParError(2);
        } else { chicos[c][i] = 0.0; chicoserr[c][i] = 0.0; sigmacos[c][i] = 0.0; sigmacoserr[c][i] = 0.0;}
      }
 
      if (hchicos_1by3bg[c][i]->Integral() > 100) {
        gstatus all1bgstat;
 
        fit(hchicos_1by3bg[c][i],  pdg.data(), all1bgstat);
        if (all1bgstat == OK) {
          chicos_1b3bg[c][i] = hchicos_1by3bg[c][i]->GetFunction("gaus")->GetParameter(1);
          chicos_1b3bgerr[c][i] = hchicos_1by3bg[c][i]->GetFunction("gaus")->GetParError(1);
          sigmacos_1b3bg[c][i] = hchicos_1by3bg[c][i]->GetFunction("gaus")->GetParameter(2);
          sigmacos_1b3bgerr[c][i] = hchicos_1by3bg[c][i]->GetFunction("gaus")->GetParError(2);
        } else { chicos_1b3bg[c][i] = 0.0; chicos_1b3bgerr[c][i] = 0.0; sigmacos_1b3bg[c][i] = 0.0; sigmacos_1b3bgerr[c][i] = 0.0;}
 
      }
 
      if (hchicos_2by3bg[c][i]->Integral() > 100) {
        gstatus all2bgstat;
        fit(hchicos_2by3bg[c][i],  pdg.data(), all2bgstat);
        if (all2bgstat == OK) {
          chicos_2b3bg[c][i] = hchicos_2by3bg[c][i]->GetFunction("gaus")->GetParameter(1);
          chicos_2b3bgerr[c][i] = hchicos_2by3bg[c][i]->GetFunction("gaus")->GetParError(1);
          sigmacos_2b3bg[c][i] = hchicos_2by3bg[c][i]->GetFunction("gaus")->GetParameter(2);
          sigmacos_2b3bgerr[c][i] = hchicos_2by3bg[c][i]->GetFunction("gaus")->GetParError(2);
        } else { chicos_2b3bg[c][i] = 0.0; chicos_2b3bgerr[c][i] = 0.0; sigmacos_2b3bg[c][i] = 0.0; sigmacos_2b3bgerr[c][i] = 0.0;}
 
      }
 
      if (hchicos_3by3bg[c][i]->Integral() > 100) {
        gstatus all3bgstat;
        fit(hchicos_3by3bg[c][i],  pdg.data(), all3bgstat);
        if (all3bgstat == OK) {
          chicos2[c][i] = hchicos_3by3bg[c][i]->GetFunction("gaus")->GetParameter(1);
          chicos2err[c][i] = hchicos_3by3bg[c][i]->GetFunction("gaus")->GetParError(1);
          sigmacos2[c][i] = hchicos_3by3bg[c][i]->GetFunction("gaus")->GetParameter(2);
          sigmacos2err[c][i] = hchicos_3by3bg[c][i]->GetFunction("gaus")->GetParError(2);
        } else { chicos2[c][i] = 0.0; chicos2err[c][i] = 0.0; sigmacos2[c][i] = 0.0; sigmacos2err[c][i] = 0.0;}
 
      }
    }
  }
 
  plotDist(hchicos_allbg,  Form("fits_chi_vscos_allbg_%s_%s", pdg.data(), suffix.data()), m_cosBins);
  plotDist(hchicos_1by3bg, Form("fits_chi_vscos_1by3bg_%s_%s", pdg.data(), suffix.data()), m_cosBins);
  plotDist(hchicos_2by3bg,  Form("fits_chi_vscos_2by3bg_%s_%s", pdg.data(), suffix.data()), m_cosBins);
  plotDist(hchicos_3by3bg, Form("fits_chi_vscos_3by3bg_%s_%s", pdg.data(), suffix.data()), m_cosBins);
 
  const double cosstep = 2.0 / m_cosBins;
  std::vector<double> cosArray(m_cosBins), cosArrayErr(m_cosBins);
  for (int i = 0; i < m_cosBins; ++i) {
    cosArray[i] = -1.0 + (i * cosstep + cosstep / 2.0); //finding bin centre
    cosArrayErr[i] = 0.0;
  }
 
  TGraphErrors grchicos(m_cosBins, cosArray.data(), chicos[0].data(), cosArrayErr.data(), chicoserr[0].data());
  TGraphErrors grchicos_1b3bg(m_cosBins, cosArray.data(), chicos_1b3bg[0].data(), cosArrayErr.data(), chicos_1b3bgerr[0].data());
  TGraphErrors grchicos_2b3bg(m_cosBins, cosArray.data(), chicos_2b3bg[0].data(), cosArrayErr.data(), chicos_2b3bgerr[0].data());
  TGraphErrors grchicos2(m_cosBins, cosArray.data(), chicos2[0].data(), cosArrayErr.data(), chicos2err[0].data());
 
  TGraphErrors grchicosn(m_cosBins, cosArray.data(), chicos[1].data(), cosArrayErr.data(), chicoserr[1].data());
  TGraphErrors grchicos_1b3bgn(m_cosBins, cosArray.data(), chicos_1b3bg[1].data(), cosArrayErr.data(), chicos_1b3bgerr[1].data());
  TGraphErrors grchicos_2b3bgn(m_cosBins, cosArray.data(), chicos_2b3bg[1].data(), cosArrayErr.data(), chicos_2b3bgerr[1].data());
  TGraphErrors grchicos2n(m_cosBins, cosArray.data(), chicos2[1].data(), cosArrayErr.data(), chicos2err[1].data());
 
  TGraphErrors grsigmacos(m_cosBins, cosArray.data(), sigmacos[0].data(), cosArrayErr.data(), sigmacoserr[0].data());
  TGraphErrors grsigmacos_1b3bg(m_cosBins, cosArray.data(), sigmacos_1b3bg[0].data(), cosArrayErr.data(),
                                sigmacos_1b3bgerr[0].data());
  TGraphErrors grsigmacos_2b3bg(m_cosBins, cosArray.data(), sigmacos_2b3bg[0].data(), cosArrayErr.data(),
                                sigmacos_2b3bgerr[0].data());
  TGraphErrors grsigmacos2(m_cosBins, cosArray.data(), sigmacos2[0].data(), cosArrayErr.data(), sigmacos2err[0].data());
 
  TGraphErrors grsigmacosn(m_cosBins, cosArray.data(), sigmacos[1].data(), cosArrayErr.data(), sigmacoserr[1].data());
  TGraphErrors grsigmacos_1b3bgn(m_cosBins, cosArray.data(), sigmacos_1b3bg[1].data(), cosArrayErr.data(),
                                 sigmacos_1b3bgerr[1].data());
  TGraphErrors grsigmacos_2b3bgn(m_cosBins, cosArray.data(), sigmacos_2b3bg[1].data(), cosArrayErr.data(),
                                 sigmacos_2b3bgerr[1].data());
  TGraphErrors grsigmacos2n(m_cosBins, cosArray.data(), sigmacos2[1].data(), cosArrayErr.data(), sigmacos2err[1].data());
 
  TLine line0(-1, 0, 1, 0);
  line0.SetLineStyle(kDashed);
  line0.SetLineColor(kRed);
 
  TLine line1(-1, 1, 1, 1);
  line1.SetLineStyle(kDashed);
  line1.SetLineColor(kRed);
 
  TString ptype("");
  double mass = 0.0;
  if (pdg == "pion") {ptype += "#pi"; mass = Const::pion.getMass();}
  else if (pdg == "kaon") {ptype += "K"; mass = Const::kaon.getMass();}
  else if (pdg == "proton") { ptype += "p"; mass = Const::proton.getMass();}
  else if (pdg == "muon") {ptype += "#mu"; mass = Const::muon.getMass();}
  else if (pdg == "electron") {ptype += "e"; mass = Const::electron.getMass();}
  if (mass == 0.0) B2FATAL("Mass of particle " << pdg.data() << " is zero");
 
  int bglow = 0, bghigh = 0;
  double bgstep = (m_bgMax - m_bgMin) / m_bgBins;
 
  TLegend lchi(0.7, 0.75, 0.8, 0.85);
  lchi.SetBorderSize(0);
  lchi.SetFillColor(kYellow);
  lchi.AddEntry(&grchicos, ptype + "^{+}", "pc");
  lchi.AddEntry(&grchicosn, ptype + "^{-}", "pc");
 
  TCanvas cchi(Form("cchi_%s", suffix.data()), "cchi", 700, 600);
  cchi.Divide(2, 2);
  cchi.cd(1);
  FormatGraph(grchicos, 0, Form("all (%s) bg bins, p =(%0.02f, %0.02f)", pdg.data(), m_bgMin * mass, m_bgMax * mass));
  FormatGraph(grchicosn, 1);
  grchicos.Draw("AP");
  grchicosn.Draw("P,same");
  line0.Draw("same");
  lchi.Draw("same");
 
  cchi.cd(2);
  bghigh = int(m_bgBins / 3);
  FormatGraph(grchicos_1b3bg, 0, Form("first 1/3 bg bins, p =(%0.02f, %0.02f)", m_bgMin * mass, bghigh * bgstep * mass));
  FormatGraph(grchicos_1b3bgn, 1);
  grchicos_1b3bg.Draw("AP");
  grchicos_1b3bgn.Draw("P,same");
  line0.Draw("same");
  lchi.Draw("same");
 
  cchi.cd(3);
  bglow = int(m_bgBins / 3);
  bghigh = int(2 * m_bgBins / 3);
  FormatGraph(grchicos_2b3bg, 0, Form("second 1/3 bg bins, p =(%0.02f, %0.02f)", bglow * bgstep * mass, bghigh * bgstep * mass));
  FormatGraph(grchicos_2b3bgn, 1);
  grchicos_2b3bg.Draw("AP");
  grchicos_2b3bgn.Draw("P,same");
  line0.Draw("same");
  lchi.Draw("same");
 
  cchi.cd(4);
  bglow = int(2 * m_bgBins / 3);
  FormatGraph(grchicos2, 0, Form("third 1/3 bg bins, p =(%0.02f, %0.02f)", bglow * bgstep * mass, m_bgMax * mass));
  FormatGraph(grchicos2n, 1);
  grchicos2.Draw("AP");
  grchicos2n.Draw("P,same");
  line0.Draw("same");
  lchi.Draw("same");
  cchi.SaveAs(Form("plots/HadronCal/Monitoring/mean_chi_vscos_inbg_%s_%s.pdf", pdg.data(), suffix.data()));
 
  cchi.cd(1);
  FormatGraph(grsigmacos, 2, Form("all (%s) bg bins, p =(%0.02f, %0.02f)", pdg.data(), m_bgMin * mass, m_bgMax * mass));
  FormatGraph(grsigmacosn, 1);
  grsigmacos.Draw("AP");
  grsigmacosn.Draw("P,same");
  line1.Draw("same");
  lchi.Draw("same");
 
  cchi.cd(2);
  bghigh = int(m_bgBins / 3);
  FormatGraph(grsigmacos_1b3bg, 2, Form("first 1/3 bg bins, p =(%0.02f, %0.02f)", m_bgMin * mass, bghigh * bgstep * mass));
  FormatGraph(grsigmacos_1b3bgn, 1);
  grsigmacos_1b3bg.Draw("AP");
  grsigmacos_1b3bgn.Draw("P,same");
  line1.Draw("same");
  lchi.Draw("same");
 
  cchi.cd(3);
  bglow = int(m_bgBins / 3);
  bghigh = int(2 * m_bgBins / 3);
  FormatGraph(grsigmacos_2b3bg, 2, Form("second 1/3 bg bins, p =(%0.02f, %0.02f)", bglow * bgstep * mass, bghigh * bgstep * mass));
  FormatGraph(grsigmacos_2b3bgn, 1);
  grsigmacos_2b3bg.Draw("AP");
  grsigmacos_2b3bgn.Draw("P,same");
  line1.Draw("same");
  lchi.Draw("same");
 
  cchi.cd(4);
  bglow = int(2 * m_bgBins / 3);
  FormatGraph(grsigmacos2, 2, Form("third 1/3 bg bins, p =(%0.02f, %0.02f)", bglow * bgstep * mass, m_bgMax * mass));
  FormatGraph(grsigmacos2n, 1);
  grsigmacos2.Draw("AP");
  grsigmacos2n.Draw("P,same");
  line1.Draw("same");
  lchi.Draw("same");
  cchi.SaveAs(Form("plots/HadronCal/Monitoring/sigma_chi_vscos_inbg_%s_%s.pdf", pdg.data(), suffix.data()));
 
}

setPars()

void setPars	(	TFile *&	outfile,
		std::string	pdg,
		std::vector< TH1F * > &	hdedx_bg,
		std::vector< TH1F * > &	hchi_bg,
		std::vector< TH1F * > &	hionzsigma_bg,
		std::map< int, std::vector< TH1F * > > &	hchi_inj )

function to fill the parameters like mean and reso in the tree

Definition at line 357 of file HadronBgPrep.cc.

{
  outfile->cd();
 
  TTree* satTree = new TTree(Form("%s", pdg.data()), "dE/dx m_means and m_errors");
  double satbg;          // beta-gamma value for this bin
  double satcosth;       // cos(theta) value for this bin
  double satdedx;        // mean dE/dx value for this bin
  double satdedxerr;     // error on ^
  double satdedxwidth;// width of ^ distribution
 
  double satchi;         // mean chi value for this bin
  double satchierr;      // error on ^
  double satchiwidth;    // width of ^ distribution
  double satchiwidth_err;
  double sationzres;   // width of dedx reso
 
  double satbg_avg;      // average beta-gamma value for this sample
  double satcosth_avg;   // average cos(theta) value for this sample
  double satdedxres_avg; // average dE/dx error squared for this sample
 
  satTree->Branch("bg", &satbg, "bg/D");
  satTree->Branch("costh", &satcosth, "costh/D");
 
  //dEdx related variables
  satTree->Branch("dedx", &satdedx, "dedx/D");
  satTree->Branch("dedxerr", &satdedxerr, "dedxerr/D");
  satTree->Branch("dedxwidth", &satdedxwidth, "dedxwidth/D");
 
  // Chi related variables
  satTree->Branch("chimean", &satchi, "chimean/D");
  satTree->Branch("chimean_err", &satchierr, "chimean_err/D");
  satTree->Branch("chisigma", &satchiwidth, "chisigma/D");
  satTree->Branch("chisigma_err", &satchiwidth_err, "chisigma_err/D");
 
  satTree->Branch("ionzres", &sationzres, "ionzres/D");
 
  // Other variables
  satTree->Branch("bg_avg", &satbg_avg, "bg_avg/D");
  satTree->Branch("costh_avg", &satcosth_avg, "costh_avg/D");
  satTree->Branch("dedxres_avg", &satdedxres_avg, "dedxres_avg/D");
 
  double bgstep = (m_bgMax - m_bgMin) / m_bgBins;
 
  for (int i = 0; i < m_bgBins; ++i) {
 
    // fill some details for this bin
    satbg = m_bgMin + 0.5 * bgstep + i * bgstep;
    satcosth = 0.0;
 
    satbg_avg = (m_sumsize[i] > 0) ? m_sumbg[i] / m_sumsize[i] : 0.0;
    satcosth_avg = (m_sumsize[i] > 0) ? m_sumcos[i] / m_sumsize[i] : 0.0;
    satdedxres_avg = (m_sumsize[i] > 0) ? m_sumres_square[i] / m_sumsize[i] : 0.0;
 
    //1. -------------------------
    // fit the dE/dx distribution in bins of beta-gamma
    gstatus bgstat;
    fit(hdedx_bg[i],  pdg.data(), bgstat);
    TF1* f = hdedx_bg[i]->GetFunction("gaus");
    if (bgstat == OK && f) {
      const auto mean = f->GetParameter(1);
      satdedx = mean;
      m_means[i] = mean;
      const auto err = f->GetParError(1);
      satdedxerr = err;
      m_errors[i] = err;
      satdedxwidth = f->GetParameter(2);
    } else { satdedx = 0.0; satdedxerr = 0.0; satdedxwidth = 0.0;}
 
    //2. -------------------------
    // fit the chi distribution  in bins of beta-gamma
    gstatus chistat;
    fit(hchi_bg[i], pdg.data(), chistat);
    f = hchi_bg[i]->GetFunction("gaus");
    if (chistat == OK && f) {
      satchi = f->GetParameter(1);
      satchierr  = f->GetParError(1);
      satchiwidth = f->GetParameter(2);
      satchiwidth_err = f->GetParError(2);
    } else { satchi = 0.0; satchierr = 0.0; satchiwidth = 0.0; satchiwidth_err = 0.0;}
 
 
    //3. -------------------------
    // fit the chi distribution  in bins of beta-gamma
    gstatus ionstat;
    fit(hionzsigma_bg[i], pdg.data(), ionstat);
    if (ionstat == OK) {
      sationzres = hionzsigma_bg[i]->GetFunction("gaus")->GetParameter(2);
    } else sationzres = 0.0;
 
    // fill the tree for this bin
    satTree->Fill();
  }
  // write out the data to file
  satTree->Write();
 
  // --------------------------------------------------
  // FIT IN BINS OF INJECTION TIME FOR VALIDATION
  // --------------------------------------------------
 
  std::string svar = "ler";
  for (int ir = 0; ir < 2; ir++) {
 
    if (ir == 1) svar = "her";
 
    TTree* injTree = new TTree(Form("%s_%s", pdg.data(), svar.data()), "chi m_means and m_errors");
    double inj_avg;
    double mean;
    double mean_err;
    double sigma;
    double sigma_err;
 
    injTree->Branch("inj_avg", &inj_avg, "inj_avg/D");
    injTree->Branch("chimean", &mean, "chimean/D");
    injTree->Branch("chimean_err", &mean_err, "chimean_err/D");
    injTree->Branch("chisigma", &sigma, "chisigma/D");
    injTree->Branch("chisigma_err", &sigma_err, "chisigma_err/D");
 
    for (int i = 0; i < m_injBins; ++i) {
 
      inj_avg = (m_injsize[i] > 0) ? m_suminj[i] / m_injsize[i] : 0.0;
 
 
      // fit the dE/dx distribution in bins of injection time'
      gstatus injstat;
      fit(hchi_inj[ir][i],  pdg.data(), injstat);
      if (injstat == OK) {
        mean = hchi_inj[ir][i]->GetFunction("gaus")->GetParameter(1);
        mean_err = hchi_inj[ir][i]->GetFunction("gaus")->GetParError(1);
        sigma = hchi_inj[ir][i]->GetFunction("gaus")->GetParameter(2);
        sigma_err = hchi_inj[ir][i]->GetFunction("gaus")->GetParError(2);
      } else {
        mean = 0.0; mean_err = 0.0; sigma = 0.0; sigma_err = 0.0;
      }
      injTree->Fill();
    }
 
    injTree->Write();
  }
  outfile->Write();
}

Member Data Documentation

m_bgBins

int m_bgBins

private

bins for bg

Definition at line 184 of file HadronBgPrep.h.

m_bgMax

double m_bgMax

private

max range of bg

Definition at line 186 of file HadronBgPrep.h.

m_bgMin

double m_bgMin

private

min range of bg

Definition at line 185 of file HadronBgPrep.h.

m_cosBins

int m_cosBins

private

bins for cosine

Definition at line 192 of file HadronBgPrep.h.

m_cosMax

double m_cosMax

private

max range of cosine

Definition at line 194 of file HadronBgPrep.h.

m_cosMin

double m_cosMin

private

min range of cosine

Definition at line 193 of file HadronBgPrep.h.

m_cut

double m_cut

private

cut to clean protons

Definition at line 200 of file HadronBgPrep.h.

m_errors

std::vector<double> m_errors

private

error variable

Definition at line 179 of file HadronBgPrep.h.

m_injBins

int m_injBins

private

bins for injection time

Definition at line 188 of file HadronBgPrep.h.

m_injMax

double m_injMax

private

max range of injection time

Definition at line 190 of file HadronBgPrep.h.

m_injMin

double m_injMin

private

min range of injection time

Definition at line 189 of file HadronBgPrep.h.

m_injsize

std::vector<int> m_injsize

private

size of the injection bins

Definition at line 182 of file HadronBgPrep.h.

m_means

std::vector<double> m_means

private

mean variable

Definition at line 178 of file HadronBgPrep.h.

m_nhitBins

int m_nhitBins

private

bins for nhits

Definition at line 196 of file HadronBgPrep.h.

m_nhitMax

double m_nhitMax

private

max range of nhits

Definition at line 198 of file HadronBgPrep.h.

m_nhitMin

double m_nhitMin

private

min range of nhits

Definition at line 197 of file HadronBgPrep.h.

m_sumbg

std::vector<double> m_sumbg

private

variables to add bg values

Definition at line 175 of file HadronBgPrep.h.

m_sumcos

std::vector<double> m_sumcos

private

variables to add cos values

Definition at line 174 of file HadronBgPrep.h.

m_suminj

std::vector<double> m_suminj

private

variables to add injection time

Definition at line 177 of file HadronBgPrep.h.

m_sumres_square

std::vector<double> m_sumres_square

private

variables to add square of resolution

Definition at line 176 of file HadronBgPrep.h.

m_sumsize

std::vector<int> m_sumsize

private

size of the bg bins

Definition at line 181 of file HadronBgPrep.h.

---

The documentation for this class was generated from the following files:

• cdc/calibration/CDCdEdx/include/HadronBgPrep.h
• cdc/calibration/CDCdEdx/src/HadronBgPrep.cc