11#include <calibration/CalibrationAlgorithm.h>
18#include <TDatabasePDG.h>
110 std::unique_ptr<TList>
GenerateNewHistograms(std::shared_ptr<TTree> ttreeLambda, std::shared_ptr<TTree> ttreeDstar,
111 std::shared_ptr<TTree> ttreeGamma, std::shared_ptr<TTree>
131 const double m_MuonPDGMass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
132 const double m_PionPDGMass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
133 const double m_KaonPDGMass = TDatabasePDG::Instance()->GetParticle(321)->Mass();
142 std::vector<double> pbins;
144 pbins.push_back(0.0);
145 pbins.push_back(0.05);
147 for (
int iBin = 2; iBin <=
m_numPBins; iBin++) {
149 pbins.push_back(0.025 + 0.025 * iBin);
151 pbins.push_back(pbins.at(19) + 0.05 * (iBin - 19));
153 pbins.push_back(pbins.at(59) + 0.3 * (iBin - 59));
165 for (
int pbin = 0; pbin <=
m_numPBins + 1; pbin++) {
166 for (
int dedxbin = 0; dedxbin <=
m_numDEdxBins + 1; dedxbin++) {
168 if (HistoToNormalise->GetBinContent(pbin, dedxbin) <= 1) {
169 HistoToNormalise->SetBinContent(pbin, dedxbin, 0);
174 TH1D* MomentumSlice = (TH1D*)HistoToNormalise->ProjectionY(
"slice_tr", pbin, pbin);
177 MomentumSlice->Scale(1. / MomentumSlice->Integral(0,
m_numDEdxBins + 1));
180 for (
int dedxbin = 0; dedxbin <=
m_numDEdxBins + 1; dedxbin++) {
181 HistoToNormalise->SetBinContent(pbin, dedxbin, MomentumSlice->GetBinContent(dedxbin));
182 HistoToNormalise->SetBinError(pbin, dedxbin, MomentumSlice->GetBinError(dedxbin));
185 return HistoToNormalise;
191 TH2F*
PrepareNewHistogram(TH2F* DataHistogram, TString NewName, TF1* betagamma_function, TF1* ResolutionFunctionOriginal,
192 double bias_correction)
194 TF1* ResolutionFunction = (TF1*) ResolutionFunctionOriginal->Clone(Form(
"%sClone",
195 ResolutionFunctionOriginal->GetName()));
197 TH2F* DataHistogramNew = (TH2F*) DataHistogram->Clone(NewName);
199 DataHistogramNew->Reset();
201 for (
int pbin = 1; pbin <=
m_numPBins + 1; pbin++) {
202 double mean_dEdx_value = betagamma_function->Eval(DataHistogramNew->GetXaxis()->GetBinCenter(pbin));
203 ResolutionFunction->FixParameter(1, mean_dEdx_value + bias_correction);
206 TH1D* MomentumSlice = (TH1D*)DataHistogramNew->ProjectionY(
"slice", pbin, pbin);
211 MomentumSlice->Fill(ResolutionFunction->GetRandom());
216 MomentumSlice->Scale(1. / MomentumSlice->Integral(0,
m_numDEdxBins + 1));
219 for (
int dedxbin = 0; dedxbin <=
m_numDEdxBins + 1; dedxbin++) {
220 DataHistogramNew->SetBinContent(pbin, dedxbin, MomentumSlice->GetBinContent(dedxbin));
221 DataHistogramNew->SetBinError(pbin, dedxbin, MomentumSlice->GetBinError(dedxbin));
225 return DataHistogramNew;
235 TF1* ResolutionFunction =
new TF1(
"ResolutionFunction",
"[0]*ROOT::Math::crystalball_function(x,[4],[3],[2]*[1],[1])", 100e3,
239 ResolutionFunction->SetNpx(1000);
241 ResolutionFunction->SetParameters(1000, 6.e5, 0.1, 1, 1);
242 ResolutionFunction->SetParLimits(0, 0, 1.e6);
243 ResolutionFunction->SetParLimits(1, 3.e5, 7.e6);
244 ResolutionFunction->SetParLimits(2, 0, 10);
245 ResolutionFunction->SetParLimits(3, 0.01, 100);
246 ResolutionFunction->SetParLimits(4, 0.01, 100);
250 TH1F* DataHistogramNew = (TH1F*)DataHistogram->ProfileX()->ProjectionX();
251 TH1F* DataHistogramClone = (TH1F*)DataHistogramNew->Clone(Form(
"%sClone",
252 DataHistogramNew->GetName()));
253 DataHistogramNew->SetName(NewName);
254 DataHistogramNew->SetTitle(NewName);
255 DataHistogramNew->Reset();
259 TH1F* MomentumSlice = (TH1F*)DataHistogram->ProjectionY(
"slice", pbin, pbin);
261 if (MomentumSlice->Integral() < 1)
continue;
263 ResolutionFunction->SetParameter(1, MomentumSlice->GetMean());
264 ResolutionFunction->SetParameter(2, MomentumSlice->GetStdDev() / MomentumSlice->GetMean());
265 ResolutionFunction->SetRange(MomentumSlice->GetMean() * 0.2, MomentumSlice->GetMean() * 1.75);
267 MomentumSlice->Fit(ResolutionFunction,
"RQI");
269 double stat_error = DataHistogramClone->GetBinError(pbin);
271 double bincontent = ResolutionFunction->GetParameter(1);
272 double binerror = ResolutionFunction->GetParError(1);
274 binerror = std::max(binerror, stat_error);
276 DataHistogramNew->SetBinContent(pbin, bincontent);
277 DataHistogramNew->SetBinError(pbin, binerror);
280 return DataHistogramNew;
EResult
The result of calibration.
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
void setMinEvtsPerTree(const double &value)
set the upper edge of the dEdx binning for the payloads
TH1F * PrepareProfile(TH2F *DataHistogram, TString NewName)
Reimplement the Profile histogram calculation for a 2D histogram.
int m_numPBins
the number of momentum bins for the payloads
const double m_MuonPDGMass
PDG mass for the muon.
const double m_ElectronPDGMass
PDG mass for the electron.
void setNumDEdxBins(const int &value)
set the number of dEdx bins for the payloads
virtual ~SVDdEdxCalibrationAlgorithm()
Destructor.
bool m_FixUnstableFitParameter
In the dEdx:betagamma fit, there is one free parameter that makes fit convergence poor.
const double m_DeuteronPDGMass
PDG mass for the deuteron.
void setFixUnstableFitParameter(bool value=true)
In the dEdx:betagamma fit, there is one free parameter that makes fit convergence poor.
std::vector< double > CreatePBinningScheme()
build the binning scheme for the momentum
SVDdEdxCalibrationAlgorithm()
Constructor.
TH2F * Normalise2DHisto(TH2F *HistoToNormalise)
Normalise a given dEdx:momentum histogram in each momentum bin, so that sum of entries in each moment...
bool m_UseProtonBGFunctionForEverything
Assume that the dEdx:betagamma trend is the same for all hadrons; use the proton trend as representat...
int m_numBGBins
the number of beta*gamma bins for the profile and fitting
void setDEdxCutoff(const double &value)
set the upper edge of the dEdx binning for the payloads
std::unique_ptr< TList > DstarHistogramming(TTree *inputTree)
produce histograms for K/pi
double m_dedxMaxPossible
the approximate max possible value of dEdx
bool m_isMakePlots
produce plots for monitoring
int m_MinEvtsPerTree
number of events in TTree below which we don't try to fit
int m_numDEdxBins
the number of dEdx bins for the payloads
TTree * LambdaMassFit(std::shared_ptr< TTree > preselTree)
Mass fit for Lambda->ppi.
std::unique_ptr< TList > LambdaHistogramming(TTree *inputTree)
produce histograms for protons
const double m_KaonPDGMass
PDG mass for the charged kaon.
int m_NToGenerate
the number of events to be generated in each momentum bin in the new payloads
TH2F * PrepareNewHistogram(TH2F *DataHistogram, TString NewName, TF1 *betagamma_function, TF1 *ResolutionFunctionOriginal, double bias_correction)
Generate a new dEdx:momentum histogram from a function that encodes dEdx:momentum trend and a functio...
std::unique_ptr< TList > GammaHistogramming(std::shared_ptr< TTree > preselTree)
produce histograms for e
void setNumBGBins(const int &value)
set the number of beta*gamma bins for the fits
TTree * DstarMassFit(std::shared_ptr< TTree > preselTree)
Mass fit for D*->Dpi.
void setCustomProfile(bool value=true)
reimplement the profile histogram calculation
bool m_UsePionBGFunctionForEverything
Assume that the dEdx:betagamma trend is the same for all hadrons; use the pion trend as representativ...
void setNumPBins(const int &value)
set the number of momentum bins for the payloads
std::unique_ptr< TList > GenerateNewHistograms(std::shared_ptr< TTree > ttreeLambda, std::shared_ptr< TTree > ttreeDstar, std::shared_ptr< TTree > ttreeGamma, std::shared_ptr< TTree > ttreeGeneric)
generate high-statistics histograms
virtual EResult calibrate() override
run algorithm on data
bool m_CustomProfile
reimplement profile histogram calculation instead of the ROOT implementation?
void setUsePionBGFunctionForEverything(bool value=false)
use the pion beta*gamma function for other hadrons
double m_dedxCutoff
the upper edge of the dEdx binning for the payloads
void setMonitoringPlots(bool value=false)
function to enable plotting
const double m_ProtonPDGMass
PDG mass for the proton.
void setUseProtonBGFunctionForEverything(bool value=false)
use the proton beta*gamma function for other hadrons
void setNEventsToGenerate(const int &value)
set the number of events to generate, per momentum bin, for the payloads
const double m_PionPDGMass
PDG mass for the charged pion.
Abstract base class for different kinds of events.