#include <RooDataHist.h>
#include <RooHistPdf.h>
#include <RooMsgService.h>
+#include <TMath.h>
#include "AliLnPt.h"
#include "B2.h"
-#ifdef HAVE_ROOUNFOLD
- #include "/opt/RooUnfold/src/RooUnfoldResponse.h"
- #include "/opt/RooUnfold/src/RooUnfoldBayes.h"
-#endif
-
ClassImp(AliLnPt)
AliLnPt::AliLnPt(const TString& particle, Double_t trigEff, const TString& inputFilename, const TString& outputFilename, const TString& otag, const TString& corrFilename, const TString& corrtag)
, fOutputTag(otag)
, fCorrFilename(corrFilename)
, fCorrTag(corrtag)
-, fLowPtBin(1)
-, fHiPtBin(19)
-, fPidM2(0)
-, fUnfolding(0)
-, fNIter(4)
+, fPtMin(0.5)
+, fPtMax(3.0)
+, fPid(0)
, fSecondaries(1)
, fEfficiency(1)
, fIsOnlyGen(0)
-, fYMin(-0.5)
-, fYMax(0.5)
-, fLowM2Bin(9)
-, fHiM2Bin(17)
-, fMinM2Bkg(2.2)
-, fMaxM2Bkg(5.)
-, fMinM2tpc(2.)
-, fMaxM2tpc(6.)
+, fPtPid(1.)
+, fBkgMin(2.2)
+, fBkgMax(5.)
+, fIntMin(2.)
+, fIntMax(6.)
, fMakeStats(1)
, fMCtoINEL(0)
, fVtxFactor(1)
, fFitFrac(1)
, fFdwnCorr(1)
, fSameFdwn(0)
+, fPidProc(kMassSquared)
+, fPidEff(1)
+, fDebugLevel(0)
{
//
// constructor
//
TH1::SetDefaultSumw2(); // switch on histogram errors
- // disable verbose in RooFit
- RooMsgService::instance().setSilentMode(1);
- RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
+ if(fDebugLevel < 3)
+ {
+ // disable verbose in RooFit
+ RooMsgService::instance().setSilentMode(1);
+ RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
+ }
}
AliLnPt::~AliLnPt()
//
// Get the true pt distribution
//
- using namespace std;
-
TFile* finput = new TFile(fInputFilename.Data(), "read");
if (finput->IsZombie()) exit(1);
TString species = fParticle;
species.ReplaceAll("Anti","");
- TH1D* hStats = (TH1D*)FindObj(finput, species + "_Stats");
+ TH1D* hStats = FindObj<TH1D>(finput, species + "_Stats");
- TH1D* hCorPt = 0;
+ TH1D* hCorrPt = 0;
if(fIsOnlyGen)
{
- hCorPt = (TH1D*)((TH1D*)FindObj(finput, fParticle + "_Gen_PhS_Prim_Pt"))->Clone( Form("%s_Cor_Pt", fParticle.Data()));
+ hCorrPt = dynamic_cast<TH1D*>( (FindObj<TH1D>(finput, fParticle + "_Gen_PhS_Prim_Pt"))->Clone( Form("%s_Cor_Pt", fParticle.Data())) );
}
- else // correct the reconstructed pt
+ else
{
fcorr = new TFile(fCorrFilename.Data(),"read");
if (fcorr->IsZombie()) exit(1);
fdebug->mkdir(fOutputTag.Data());
- // -------------- get pointers to the required histograms ---------
+ // pointers to the required histograms
- // pt distributions
- TH1D* hOrigPt = (TH1D*)FindObj(finput, fParticle + "_PID_Pt");
- TH2D* hM2Pt = (TH2D*)FindObj(finput, fParticle + "_PID_M2_Pt");
+ TH1D* hOrigPt = FindObj<TH1D>(finput, fParticle + "_PID_Pt");
- // corrections
-#ifdef HAVE_ROOUNFOLD
- TH2D* hResponseMtx = (TH2D*)FindObj(fcorr, fCorrTag, fParticle + "_Response_Matrix");
- TH1D* hMeasuredPt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_Measured_Pt");
- TH1D* hTruePt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_True_Pt");
-#endif
+ TH2D* hPidPt = 0;
+ if(fPidProc == kTimeOfFlight) hPidPt = FindObj<TH2D>(finput, fParticle + "_PID_DTime_Pt");
+ else if(fPidProc == kMassSquared) hPidPt = FindObj<TH2D>(finput, fParticle + "_PID_M2_Pt");
+ else if(fPidProc == kMassDifference) hPidPt = FindObj<TH2D>(finput, fParticle + "_PID_DM2_Pt");
+ else fPid = 0;
- TH1D* hFracMatPt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_Frac_Mat_Pt");
- TH1D* hFracFdwnPt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_Frac_Fdwn_Pt");
- TF1* fncFracMatPt = (TF1*)FindObj(fcorr, fCorrTag, fParticle + "_Frac_Mat_Fit_Pt");
- TF1* fncFracFdwnPt = (TF1*)FindObj(fcorr, fCorrTag, fParticle + "_Frac_Fdwn_Fit_Pt");
+ // corrections
+ TH1D* hFracMatPt = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Frac_Mat_Pt");
+ TH1D* hFracFdwnPt = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Frac_Fdwn_Pt");
+ TF1* fncFracMatPt = FindObj<TF1>(fcorr, fCorrTag, fParticle + "_Frac_Mat_Fit_Pt");
+ TF1* fncFracFdwnPt = FindObj<TF1>(fcorr, fCorrTag, fParticle + "_Frac_Fdwn_Fit_Pt");
if(fSameFdwn)
{
- hFracFdwnPt = (TH1D*)FindObj(fcorr, fCorrTag, species + "_Frac_Fdwn_Pt");
- fncFracFdwnPt = (TF1*)FindObj(fcorr, fCorrTag, species + "_Frac_Fdwn_Fit_Pt");
+ hFracFdwnPt = FindObj<TH1D>(fcorr, fCorrTag, species + "_Frac_Fdwn_Pt");
+ fncFracFdwnPt = FindObj<TF1>(fcorr, fCorrTag, species + "_Frac_Fdwn_Fit_Pt");
}
- TH1D* hEffVtxPt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_Eff_Vtx_Pt");
- TH1D* hEffAccTrkPt = (TH1D*)FindObj(fcorr, fCorrTag, fParticle + "_Eff_AccTrk_Pt");
+ TH1D* hEffVtxPt = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Eff_Vtx_Pt");
+ TH1D* hEffAccTrkPt = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Eff_AccTrk_Pt");
- // -------------- apply corrections -----------
+ // apply corrections
fdebug->cd(fOutputTag.Data());
- TH1D* hPidPt = 0;
- if(fPidM2 && (fHiPtBin > fLowM2Bin))
+ TH1D* hPidCorrPt = 0;
+ if(fPid)
{
- Int_t hibin = (fHiPtBin > fHiM2Bin) ? fHiM2Bin : fHiPtBin;
- hPidPt = this->PID(hOrigPt, hM2Pt, fLowM2Bin, hibin, fParticle + "_M2Corr_Pt");
+ if(fDebugLevel>1) std::cout << "\t\tpid correction" << std::endl;
+
+ hPidCorrPt = this->PID(hOrigPt, hPidPt, fPtPid, fPtMax, fParticle + "_PidCorr_Pt");
+ hPidCorrPt->Scale(1./fPidEff);
}
else
{
- hPidPt = (TH1D*)hOrigPt->Clone(Form("%s_M2Corr_Pt", fParticle.Data()));
+ hPidCorrPt = dynamic_cast<TH1D*>(hOrigPt->Clone(Form("%s_PidCorr_Pt", fParticle.Data())));
}
- TH1D* hSecCorPt = 0;
+ TH1D* hSecCorrPt = 0;
if(fSecondaries)
{
+ if(fDebugLevel>1) std::cout << "\t\tremoval of secondaries" << std::endl;
+
if(fFitFrac)
{
- hSecCorPt = this->Secondaries(hPidPt, fncFracMatPt, fncFracFdwnPt, fParticle + "_SecCor_Pt");
+ hSecCorrPt = this->Secondaries(hPidCorrPt, fncFracMatPt, fncFracFdwnPt, fParticle + "_SecCorr_Pt");
}
else
{
- hSecCorPt = this->Secondaries(hPidPt, hFracMatPt, hFracFdwnPt, fParticle + "_SecCor_Pt");
+ hSecCorrPt = this->Secondaries(hPidCorrPt, hFracMatPt, hFracFdwnPt, fParticle + "_SecCorr_Pt");
}
}
else
{
- hSecCorPt = (TH1D*)hPidPt->Clone(Form("%s_SecCor_Pt", fParticle.Data()));
+ hSecCorrPt = dynamic_cast<TH1D*>(hPidCorrPt->Clone(Form("%s_SecCorr_Pt", fParticle.Data())));
}
- TH1D* hUnfoldedPt = 0;
-#ifdef HAVE_ROOUNFOLD
- if(fUnfolding)
- {
- hUnfoldedPt = this->Unfolding(hSecCorPt, hMeasuredPt, hTruePt, hResponseMtx, fParticle + "_Unfolded_Pt", fNIter);
- }
- else
-#endif
- {
- hUnfoldedPt = (TH1D*)hSecCorPt->Clone(Form("%s_Unfolded_Pt", fParticle.Data()));
- }
-
- TH1D* hEffCorPt = 0;
+ TH1D* hEffCorrPt = 0;
if(fEfficiency)
{
+ if(fDebugLevel>1) std::cout << "\t\tefficiency correction" << std::endl;
+
if(fMCtoINEL)
{
- TH1D* hStatsMC = (TH1D*)FindObj(fcorr, species + "_Stats");
+ TH1D* hStatsMC = FindObj<TH1D>(fcorr, fCorrTag, species + "_Stats");
fVtxFactor = this->GetVertexCorrection(hStats, hStatsMC);
}
- hEffCorPt = this->Efficiency(hUnfoldedPt, hEffVtxPt, hEffAccTrkPt, fParticle + "_EffCor_Pt");
+ hEffCorrPt = this->Efficiency(hSecCorrPt, hEffVtxPt, hEffAccTrkPt, fParticle + "_EffCorr_Pt");
}
else
{
- hEffCorPt = (TH1D*)hUnfoldedPt->Clone(Form("%s_EffCor_Pt", fParticle.Data()));
+ hEffCorrPt = dynamic_cast<TH1D*>(hSecCorrPt->Clone(Form("%s_EffCorr_Pt", fParticle.Data())));
}
- // --------- corrected pt --------
+ // corrected pt
- hCorPt = (TH1D*) hEffCorPt->Clone(Form("%s_Cor_Pt", fParticle.Data()));
+ hCorrPt = dynamic_cast<TH1D*>(hEffCorrPt->Clone(Form("%s_Corr_Pt", fParticle.Data())));
- // ---------- debug -----
+ // debug histograms
fdebug->cd(fOutputTag.Data());
hOrigPt->Write(); // original distribution
- hPidPt->Write();
- hUnfoldedPt->Write();
- hSecCorPt->Write();
- hEffCorPt->Write();
- hCorPt->Write();
-
- // ------------- clean --------------
-
- delete hPidPt;
- delete hUnfoldedPt;
- delete hSecCorPt;
- delete hEffCorPt;
+ hPidCorrPt->Write();
+ hSecCorrPt->Write();
+ hEffCorrPt->Write();
+ hCorrPt->Write();
+
+ // free memory
+ delete hPidCorrPt;
+ delete hSecCorrPt;
+ delete hEffCorrPt;
}
- // ---------- save the final pt ----------
+ // save final pt
- TH1D* hPt = (TH1D*) hCorPt->Clone(Form("%s_Pt", fParticle.Data()));
+ TH1D* hPt = dynamic_cast<TH1D*>(hCorrPt->Clone(Form("%s_Pt", fParticle.Data())));
hPt->SetTitle(fParticle.Data());
hPt->Reset();
hPt->Sumw2();
- for(Int_t i=fLowPtBin; i<fHiPtBin; ++i)
+ Int_t lowbin = hPt->GetXaxis()->FindFixBin(fPtMin);
+ Int_t hibin = hPt->GetXaxis()->FindFixBin(fPtMax);
+
+ for(Int_t i=lowbin; i<hibin; ++i)
{
- hPt->SetBinContent(i, hCorPt->GetBinContent(i));
- hPt->SetBinError(i, hCorPt->GetBinError(i));
+ hPt->SetBinContent(i, hCorrPt->GetBinContent(i));
+ hPt->SetBinError(i, hCorrPt->GetBinError(i));
}
foutput->cd();
hPt->Write();
delete hPt;
- delete hCorPt;
+ delete hCorrPt;
delete fdebug;
delete fcorr;
- // --------- stats -------------------
+ // stats
if(fMakeStats)
{
Double_t AliLnPt::GetM2Width(Double_t pt, Double_t m) const
{
//
-// parameterization of the expected M2 width
-// as a function of pt from LHC12a5b
+// expected M2 width parameterization from LHC12a5b
//
Double_t c0 = 1.90212e-03;
Double_t c1 = 1.29814e-02;
return c0/(pt*pt) + c1*(pt*pt/(m*m) + 1.);
}
-RooWorkspace* AliLnPt::GetM2Model(Double_t pt, Double_t m, const TString& name, Double_t max) const
+Double_t AliLnPt::GetExpectedDT(Double_t pt, Double_t m) const
{
//
-// deuteron m2 model for TPC pid
-// gaussian + exponential background of pi, K and p
+// estimate of the expected time of flight (ns) for mass hypothesis m
+//
+ Double_t l = 370.; // cm
+ Double_t c = 2.99792458e+1; // cm/ns
+
+ return (l/c)*TMath::Sqrt(1.+(m/pt)*(m/pt));
+}
+
+RooWorkspace* AliLnPt::GetToFModel(Double_t pt, const TString& name) const
+{
+//
+// model for the time of flight distribution
//
using namespace RooFit;
RooWorkspace* w = new RooWorkspace(name.Data());
- RooRealVar m2("m2", "m^{2} (GeV^{2}/c^{4})", fMinM2Bkg, fMaxM2Bkg);
- w->import(m2);
+ // signal
+ w->factory(Form("AliLnGaussianExpTail::Sd(x[%f,%f], mu[0,-0.06,0.06], sigma[0.120,0.05,0.3], tau[0.1,0.08,0.4])", fBkgMin, fBkgMax));
- Double_t sigma = GetM2Width(pt, m);
+ // background
+ Double_t expTd = this->GetExpectedDT(pt, GetMass("deuteron"));
+ Double_t expDTp = this->GetExpectedDT(pt, GetMass("proton")) - expTd;
+ Double_t expDTk = this->GetExpectedDT(pt, GetMass("kaon")) - expTd;
+ Double_t expDTpi = this->GetExpectedDT(pt, GetMass("pion")) - expTd;
+
+ w->factory(Form("AliLnGaussianExpTail::ProtonBkg(x, muP[%f,%f-0.2,%f+0.2], widthP[0.120,0.05,0.3], tauP[0.08,0.01,0.2])",expDTp, expDTp, expDTp));
+ w->factory(Form("AliLnGaussianExpTail::KaonBkg(x, muK[%f,%f-0.2,%f+0.2], widthK[0.120,0.05,0.3], tauP[0.08,0.01,0.2])",expDTk, expDTk, expDTk));
+ w->factory(Form("AliLnGaussianExpTail::PionBkg(x, muPi[%f,%f-0.2,%f+0.2], widthPi[0.120,0.05,0.3], tauPi[0.08,0.01,0.2])",expDTpi, expDTpi, expDTpi));
+
+ w->factory("SUM::Bkg( Npi[0.3,0,1.e+9]*PionBkg, Nk[0.24,0,1.e+9]*KaonBkg, Np[0.47,0,1.e+9]*ProtonBkg, Np[0.47,0,1.e+9]*ProtonBkg)");
+
+ w->factory("SUM::model( Nd[0,1.e+9]*Sd, Nbkg[0,1.e+9]*Bkg )");
+
+ w->var("x")->SetTitle("t - <t> (ns)");
+
+ return w;
+}
+
+RooWorkspace* AliLnPt::GetM2Model(Double_t pt, Double_t m, const TString& name, Double_t max) const
+{
+//
+// mass squared model for deuterons
+//
+ using namespace RooFit;
- // signal (m is shifted to the left)
- w->factory(Form("Gaussian::Sd(m2, mean[%f,%f,%f], width[%f,%f,%f])",m*m, m*m-0.3, m*m+0.05, sigma, sigma-0.05, sigma+0.05));
+ RooWorkspace* w = new RooWorkspace(name.Data());
// background
- w->factory("Exponential::PionBkg(m2,lambda1[-0.5,-100.,0.])");
- w->factory("Exponential::KaonBkg(m2,lambda2[-0.5,-100.,0.])");
- w->factory("Exponential::ProtonBkg(m2,lambda3[-0.5,-100.,0.])");
- w->factory(Form("SUM::Bkg( Npi[0.3,0,%f]*PionBkg, Nk[0.24,0,%f]*KaonBkg, Np[0.47,0,%f]*ProtonBkg )",max,max,max));
+ w->factory(Form("Exponential::PionBkg(x[%f,%f],lambdaPi[-0.3,-10.,0.])", fBkgMin, fBkgMax));
+ w->factory("Exponential::KaonBkg(x,lambdaK[-0.4,-10.,0.])");
+ w->factory("Exponential::ProtonBkg(x,lambdaP[-0.5,-10.,0.])");
+
+ w->factory("SUM::Bkg( Npi[0.3,0.,1.]*PionBkg, Nk[0.24,0.,1.]*KaonBkg, Np[0.47,0.,1.]*ProtonBkg )");
- w->factory(Form("SUM::model( Nd[0,%f]*Sd, Nbkg[0,%f]*Bkg )",max,max));
+ // signal
+
+ Double_t sigma = GetM2Width(pt, m);
+ Double_t mm = (fPidProc == kMassSquared) ? m*m : 0;
+
+ if(pt<2.0)
+ {
+ w->factory(Form("AliLnM2::Sd(x,mu[%f,%f,%f], sigma[%f,%f,%f],tau[0.1,0.05,0.4],p[%f,%f,%f])",mm, mm-0.15, mm+0.1, sigma, 0.1, sigma+0.1, 1.16*pt,pt,1.33*pt));
+ }
+ else // neglect non-gaussian tails
+ {
+ w->factory(Form("Gaussian::Sd(x, mu[%f,%f,%f], sigma[%f,%f,%f])",mm, mm-0.15, mm+0.1, sigma, 0.1, sigma+0.1));
+ }
+
+ // model
+ w->factory(Form("SUM::model( Nd[0.,%f]*Sd, Nbkg[0.,%f]*Bkg )", max, max));
+
+ w->var("x")->SetTitle("#it{m}^{2} (GeV^{2}/c^{4})");
return w;
}
-void AliLnPt::GetPtFromM2(TH1D* hPt, Int_t lowbin, Int_t hibin, const TH2D* hM2Pt, Double_t m2min, Double_t m2max) const
+void AliLnPt::GetPtFromPid(TH1D* hPt, Double_t ptmin, Double_t ptmax, const TH2D* hPidPt, Double_t pidMin, Double_t pidMax) const
{
//
-// Fill the pt from lowbin to hibin integrating
-// the m2-pt distribution from m2min to m2max
+// Integrate pid distribution
//
+ Int_t lowbin = hPt->GetXaxis()->FindFixBin(ptmin);
+ Int_t hibin = hPt->GetXaxis()->FindFixBin(ptmax);
+
for(Int_t i=lowbin; i<hibin; ++i)
{
- TH1D* hM2 = hM2Pt->ProjectionY(Form("_M2_%02d",i),i,i);
+ TH1D* hPid = hPidPt->ProjectionY(Form("%s_PID_%02d",fParticle.Data(),i),i,i);
- Int_t m2minbin = hM2->GetXaxis()->FindFixBin(m2min);
- Int_t m2maxbin = hM2->GetXaxis()->FindFixBin(m2max); // m2 has a non-gaussian tail on the right
+ Int_t pidLowBin = hPid->GetXaxis()->FindFixBin(pidMin);
+ Int_t pidHiBin = hPid->GetXaxis()->FindFixBin(pidMax);
Double_t err = 0;
- Double_t n = hM2->IntegralAndError(m2minbin, m2maxbin, err);
+ Double_t n = hPid->IntegralAndError(pidLowBin, pidHiBin, err);
hPt->SetBinContent(i, n);
hPt->SetBinError(i, err);
- delete hM2;
+ delete hPid;
}
}
-TH1D* AliLnPt::PID(const TH1D* hPt, const TH2D* hM2Pt, Int_t lowbin, Int_t hibin, const TString& name)
+TH1D* AliLnPt::PID(const TH1D* hPt, const TH2D* hPidPt2, Double_t ptmin, Double_t ptmax, const TString& name)
{
//
-// Remove contamination in the m^2 distribution
+// Remove contamination in the pt distribution
//
+ TH1D* hPidCorrPt = dynamic_cast<TH1D*>(hPt->Clone(name.Data()));
+ if(hPidCorrPt == 0) return 0;
+
+ hPidCorrPt->Reset();
+ hPidCorrPt->Sumw2();
+
+ // set binning
const Int_t kNBin = 2;
- Double_t m = GetMass(fParticle);
- TH1D* hPidPt = (TH1D*) hPt->Clone(name.Data());
- hPidPt->Reset();
- hPidPt->Sumw2();
+ TH2D* hPidPt = dynamic_cast<TH2D*>(hPidPt2->Clone("_pid_pt_"));
+ hPidPt->Rebin2D(1,kNBin);
- // first bins are not contaminated since the identification is clear
- // integrate around the m2 value
+ // integrate around the mean value
- this->GetPtFromM2(hPidPt, 0, lowbin, hM2Pt, fMinM2tpc, fMaxM2tpc);
+ this->GetPtFromPid(hPidCorrPt, 0, ptmin, hPidPt, fIntMin, fIntMax);
- // the remaining bins are contaminated
- // slice the m2-pt and fit each slice with a gaussian + exponentials bkg
+ // slice the m2-pt
using namespace RooFit;
+ Int_t lowbin = hPt->GetXaxis()->FindFixBin(ptmin);
+ Int_t hibin = hPt->GetXaxis()->FindFixBin(ptmax);
+
+ Double_t m = GetMass(fParticle);
+
for(Int_t i=lowbin; i<hibin; ++i)
{
- TH1D* hM2 = hM2Pt->ProjectionY(Form("%s_PID_M2_%02d",fParticle.Data(),i),i,i);
- hM2->Rebin(kNBin);
+ TH1D* hPid = hPidPt->ProjectionY(Form("%s_PID_%02d",fParticle.Data(),i),i,i);
- RooWorkspace* w = GetM2Model(hPt->GetBinCenter(i), m, Form("%s_M2_%02d",fParticle.Data(),i),hM2->Integral());
+ RooWorkspace* w = (fPidProc==kTimeOfFlight) ? GetToFModel( hPt->GetBinCenter(i), Form("%s_M2_%02d",fParticle.Data(),i)) : GetM2Model(hPt->GetBinCenter(i), m, Form("%s_M2_%02d",fParticle.Data(),i),hPid->Integral());
- RooRealVar* m2 = w->var("m2");
- RooDataHist data("data", "data to fit", *m2, hM2);
+ RooRealVar* x = w->var("x");
+ RooDataHist data("data", "data to fit", *x, hPid);
w->import(data);
w->pdf("model")->fitTo(data, Minos(kTRUE), PrintLevel(-1), Verbose(kFALSE), Warnings(kFALSE), PrintEvalErrors(-1) );
- hPidPt->SetBinContent(i, w->var("Nd")->getVal());
- hPidPt->SetBinError(i, w->var("Nd")->getError());
+ hPidCorrPt->SetBinContent(i, w->var("Nd")->getVal());
+ hPidCorrPt->SetBinError(i, w->var("Nd")->getError());
w->Write();
delete w;
- delete hM2;
+ delete hPid;
}
- return hPidPt;
+ delete hPidPt;
+ return hPidCorrPt;
}
-#ifdef HAVE_ROOUNFOLD
-
-TH1D* AliLnPt::Unfolding(const TH1D* hPt, const TH1D* hMeasuredPt, const TH1D* hTruePt, const TH2D* hResponseMtx, const TString& name, Int_t iterations) const
-{
-//
-// Bayesian unfolding
-//
- RooUnfoldResponse* matrix = new RooUnfoldResponse(hMeasuredPt, hTruePt, hResponseMtx);
-
- RooUnfoldBayes* unfold = new RooUnfoldBayes(matrix, hPt, iterations);
-
- TH1D* hUnfoldedPt = (TH1D*)unfold->Hreco();
- hUnfoldedPt->SetName(name.Data());
- hUnfoldedPt->SetXTitle("p_{T} (GeV/c)");
-
- delete unfold;
- delete matrix;
-
- return hUnfoldedPt;
-}
-
-#endif
-
TH1D* AliLnPt::Secondaries(const TH1D* hPt, const TH1D* hFracMatPt, const TH1D* hFracFdwnPt, const TString& name) const
{
//
TH1D* AliLnPt::Secondaries(const TH1D* hPt, TF1* fncMatPt, TF1* fncFdwnPt, const TString& name) const
{
//
-// remove contamination of secondaries using fitted fractions
+// remove contamination of secondaries using a function fitted to the fractions
// N_prim = N/(1+fmat+fdwn)
//
TF1* fnc = new TF1("_secondaries_","1 + [0]*exp(-[1]*x)+[2] + [3]*exp(-[4]*x)+[5]", 0, 10);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff