Belle II Software development
SVDdEdxCalibrationAlgorithm.h
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#pragma once
10
11#include <calibration/CalibrationAlgorithm.h>
12
13#include <TH2F.h>
14#include <TF1.h>
15#include <TList.h>
16#include <TTree.h>
17#include <TProfile.h>
18#include <TDatabasePDG.h>
19
20namespace Belle2 {
29 public:
33
35
39
41
45 void setMonitoringPlots(bool value = false) { m_isMakePlots = value; }
46
50 void setNumDEdxBins(const int& value) { m_numDEdxBins = value; }
51
55 void setNumPBins(const int& value) { m_numPBins = value; }
56
60 void setNumBGBins(const int& value) { m_numBGBins = value; }
61
65 void setDEdxCutoff(const double& value) { m_dedxCutoff = value; }
66
70 void setMinEvtsPerTree(const double& value) { m_MinEvtsPerTree = value; }
71
75 void setNEventsToGenerate(const int& value) { m_NToGenerate = value; }
76
80 void setCustomProfile(bool value = true) { m_CustomProfile = value; }
81
85 void setFixUnstableFitParameter(bool value = true) { m_FixUnstableFitParameter = value; }
86
91
96
97 protected:
101 virtual EResult calibrate() override;
102
103 private:
105 TTree* LambdaMassFit(std::shared_ptr<TTree> preselTree);
106 std::unique_ptr<TList> LambdaHistogramming(TTree* inputTree);
107 TTree* DstarMassFit(std::shared_ptr<TTree> preselTree);
108 std::unique_ptr<TList> DstarHistogramming(TTree* inputTree);
109 std::unique_ptr<TList> GammaHistogramming(std::shared_ptr<TTree> preselTree);
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>
112 ttreeGeneric);
113 int m_numDEdxBins = 100;
114 int m_numPBins = 69;
115 int m_numBGBins = 69;
116 double m_dedxCutoff = 5.e6;
117 double m_dedxMaxPossible = 7.e6;
119 100;
121 500000;
124 0;
126 0;
128 1;
129
130 const double m_ElectronPDGMass = TDatabasePDG::Instance()->GetParticle(11)->Mass();
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();
134 const double m_ProtonPDGMass = TDatabasePDG::Instance()->GetParticle(2212)->Mass();
135 const double m_DeuteronPDGMass = TDatabasePDG::Instance()->GetParticle(1000010020)->Mass();
136
140 std::vector<double> CreatePBinningScheme()
141 {
142 std::vector<double> pbins;
143 pbins.reserve(m_numPBins + 1);
144 pbins.push_back(0.0);
145 pbins.push_back(0.05);
146
147 for (int iBin = 2; iBin <= m_numPBins; iBin++) {
148 if (iBin <= 19)
149 pbins.push_back(0.025 + 0.025 * iBin);
150 else if (iBin <= 59)
151 pbins.push_back(pbins.at(19) + 0.05 * (iBin - 19));
152 else
153 pbins.push_back(pbins.at(59) + 0.3 * (iBin - 59));
154 }
155
156 return pbins;
157 }
158
163 TH2F* Normalise2DHisto(TH2F* HistoToNormalise)
164 {
165 for (int pbin = 0; pbin <= m_numPBins + 1; pbin++) {
166 for (int dedxbin = 0; dedxbin <= m_numDEdxBins + 1; dedxbin++) {
167 // get rid of the bins with negative weights
168 if (HistoToNormalise->GetBinContent(pbin, dedxbin) <= 1) {
169 HistoToNormalise->SetBinContent(pbin, dedxbin, 0);
170 };
171 }
172 // create a projection (1D histogram) in a given momentum bin
173
174 TH1D* MomentumSlice = (TH1D*)HistoToNormalise->ProjectionY("slice_tr", pbin, pbin);
175 // normalise, but ignore the cases with empty histograms
176 if (MomentumSlice->Integral(0, m_numDEdxBins + 1) > 0) {
177 MomentumSlice->Scale(1. / MomentumSlice->Integral(0, m_numDEdxBins + 1));
178 }
179 // fill back the 2D histogram with the result
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));
183 }
184 }
185 return HistoToNormalise;
186 }
187
191 TH2F* PrepareNewHistogram(TH2F* DataHistogram, TString NewName, TF1* betagamma_function, TF1* ResolutionFunctionOriginal,
192 double bias_correction)
193 {
194 TF1* ResolutionFunction = (TF1*) ResolutionFunctionOriginal->Clone(Form("%sClone",
195 ResolutionFunctionOriginal->GetName())); // to avoid modifying the resolution function
196 ResolutionFunction->SetRange(0, m_dedxMaxPossible); // allow the function to take values outside the histogram range
197 TH2F* DataHistogramNew = (TH2F*) DataHistogram->Clone(NewName);
198
199 DataHistogramNew->Reset();
200
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);
204
205 // create a projection (1D histogram) in a given momentum bin
206 TH1D* MomentumSlice = (TH1D*)DataHistogramNew->ProjectionY("slice", pbin, pbin);
207
208 // fill manually (instead of FillRandom) to also preserve events in the overflow bin
209 // this is needed for the correct normalisation
210 for (int iEvent = 0; iEvent < m_NToGenerate; iEvent++) {
211 MomentumSlice->Fill(ResolutionFunction->GetRandom());
212 }
213
214 // normalise each momentum slice to unity, but ignore the cases with empty histograms
215 if (MomentumSlice->Integral(0, m_numDEdxBins + 1) > 0) {
216 MomentumSlice->Scale(1. / MomentumSlice->Integral(0, m_numDEdxBins + 1));
217 }
218 // fill back the 2D histo with the result
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));
222 }
223 }
224
225 return DataHistogramNew;
226
227 }
228
232 TH1F* PrepareProfile(TH2F* DataHistogram, TString NewName)
233 {
234// define our resolution function: Crystal Ball. Parameter [1] is the mean and [2] is the relative width.
235 TF1* ResolutionFunction = new TF1("ResolutionFunction", "[0]*ROOT::Math::crystalball_function(x,[4],[3],[2]*[1],[1])", 100e3,
236 7000e3);
237
238
239 ResolutionFunction->SetNpx(1000);
240
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);
247
248 ResolutionFunction->SetRange(0, m_dedxMaxPossible); // allow the function to take values outside the histogram range
249
250 TH1F* DataHistogramNew = (TH1F*)DataHistogram->ProfileX()->ProjectionX();
251 TH1F* DataHistogramClone = (TH1F*)DataHistogramNew->Clone(Form("%sClone",
252 DataHistogramNew->GetName())); // preserve the original profile for uncertainty cross-checks
253 DataHistogramNew->SetName(NewName);
254 DataHistogramNew->SetTitle(NewName);
255 DataHistogramNew->Reset();
256
257 for (int pbin = 1; pbin <= m_numBGBins; pbin++) {
258 // create a projection (1D histogram) in a given momentum bin
259 TH1F* MomentumSlice = (TH1F*)DataHistogram->ProjectionY("slice", pbin, pbin);
260
261 if (MomentumSlice->Integral() < 1) continue;
262// guesstimate the starting fit values
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);
266// fit each slice to extract the mean
267 MomentumSlice->Fit(ResolutionFunction, "RQI");
268
269 double stat_error = DataHistogramClone->GetBinError(pbin);
270// fill back the 1D histo with the result
271 double bincontent = ResolutionFunction->GetParameter(1);
272 double binerror = ResolutionFunction->GetParError(1);
273
274 binerror = std::max(binerror, stat_error);
275
276 DataHistogramNew->SetBinContent(pbin, bincontent);
277 DataHistogramNew->SetBinError(pbin, binerror);
278 }
279
280 return DataHistogramNew;
281
282 }
283
284 };
285
286} // namespace Belle2
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
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
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.