Belle II Software development
SVDdEdxCalibrationAlgorithm Class Reference

Class implementing the SVD dEdx calibration algorithm. More...

#include <SVDdEdxCalibrationAlgorithm.h>

Inheritance diagram for SVDdEdxCalibrationAlgorithm:
CalibrationAlgorithm

Public Types

enum  EResult {
  c_OK ,
  c_Iterate ,
  c_NotEnoughData ,
  c_Failure ,
  c_Undefined
}
 The result of calibration. More...
 

Public Member Functions

 SVDdEdxCalibrationAlgorithm ()
 Constructor.
 
virtual ~SVDdEdxCalibrationAlgorithm ()
 Destructor.
 
void setMonitoringPlots (bool value=false)
 function to enable plotting
 
void setNumDEdxBins (const int &value)
 set the number of dEdx bins for the payloads
 
void setNumPBins (const int &value)
 set the number of momentum bins for the payloads
 
void setNumBGBins (const int &value)
 set the number of beta*gamma bins for the fits
 
void setDEdxCutoff (const double &value)
 set the upper edge of the dEdx binning for the payloads
 
void setMinEvtsPerTree (const double &value)
 set the upper edge of the dEdx binning for the payloads
 
void setNEventsToGenerate (const int &value)
 set the number of events to generate, per momentum bin, for the payloads
 
void setCustomProfile (bool value=true)
 reimplement the profile histogram calculation
 
void setFixUnstableFitParameter (bool value=true)
 In the dEdx:betagamma fit, there is one free parameter that makes fit convergence poor.
 
void setUsePionBGFunctionForEverything (bool value=false)
 use the pion beta*gamma function for other hadrons
 
void setUseProtonBGFunctionForEverything (bool value=false)
 use the proton beta*gamma function for other hadrons
 
std::string getPrefix () const
 Get the prefix used for getting calibration data.
 
bool checkPyExpRun (PyObject *pyObj)
 Checks that a PyObject can be successfully converted to an ExpRun type.
 
Calibration::ExpRun convertPyExpRun (PyObject *pyObj)
 Performs the conversion of PyObject to ExpRun.
 
std::string getCollectorName () const
 Alias for prefix.
 
void setPrefix (const std::string &prefix)
 Set the prefix used to identify datastore objects.
 
void setInputFileNames (PyObject *inputFileNames)
 Set the input file names used for this algorithm from a Python list.
 
PyObject * getInputFileNames ()
 Get the input file names used for this algorithm and pass them out as a Python list of unicode strings.
 
std::vector< Calibration::ExpRun > getRunListFromAllData () const
 Get the complete list of runs from inspection of collected data.
 
RunRange getRunRangeFromAllData () const
 Get the complete RunRange from inspection of collected data.
 
IntervalOfValidity getIovFromAllData () const
 Get the complete IoV from inspection of collected data.
 
void fillRunToInputFilesMap ()
 Fill the mapping of ExpRun -> Files.
 
std::string getGranularity () const
 Get the granularity of collected data.
 
EResult execute (std::vector< Calibration::ExpRun > runs={}, int iteration=0, IntervalOfValidity iov=IntervalOfValidity())
 Runs calibration over vector of runs for a given iteration.
 
EResult execute (PyObject *runs, int iteration=0, IntervalOfValidity iov=IntervalOfValidity())
 Runs calibration over Python list of runs. Converts to C++ and then calls the other execute() function.
 
std::list< Database::DBImportQuery > & getPayloads ()
 Get constants (in TObjects) for database update from last execution.
 
std::list< Database::DBImportQuerygetPayloadValues ()
 Get constants (in TObjects) for database update from last execution but passed by VALUE.
 
bool commit ()
 Submit constants from last calibration into database.
 
bool commit (std::list< Database::DBImportQuery > payloads)
 Submit constants from a (potentially previous) set of payloads.
 
const std::string & getDescription () const
 Get the description of the algorithm (set by developers in constructor)
 
bool loadInputJson (const std::string &jsonString)
 Load the m_inputJson variable from a string (useful from Python interface). The return bool indicates success or failure.
 
const std::string dumpOutputJson () const
 Dump the JSON string of the output JSON object.
 
const std::vector< Calibration::ExpRun > findPayloadBoundaries (std::vector< Calibration::ExpRun > runs, int iteration=0)
 Used to discover the ExpRun boundaries that you want the Python CAF to execute on. This is optional and only used in some.
 
template<>
std::shared_ptr< TTree > getObjectPtr (const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
 Specialization of getObjectPtr<TTree>.
 

Protected Member Functions

virtual EResult calibrate () override
 run algorithm on data
 
void setInputFileNames (std::vector< std::string > inputFileNames)
 Set the input file names used for this algorithm.
 
virtual bool isBoundaryRequired (const Calibration::ExpRun &)
 Given the current collector data, make a decision about whether or not this run should be the start of a payload boundary.
 
virtual void boundaryFindingSetup (std::vector< Calibration::ExpRun >, int)
 If you need to make some changes to your algorithm class before 'findPayloadBoundaries' is run, make them in this function.
 
virtual void boundaryFindingTearDown ()
 Put your algorithm back into a state ready for normal execution if you need to.
 
const std::vector< Calibration::ExpRun > & getRunList () const
 Get the list of runs for which calibration is called.
 
int getIteration () const
 Get current iteration.
 
std::vector< std::string > getVecInputFileNames () const
 Get the input file names used for this algorithm as a STL vector.
 
template<class T>
std::shared_ptr< T > getObjectPtr (const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
 Get calibration data object by name and list of runs, the Merge function will be called to generate the overall object.
 
template<class T>
std::shared_ptr< T > getObjectPtr (std::string name)
 Get calibration data object (for all runs the calibration is requested for) This function will only work during or after execute() has been called once.
 
template<>
shared_ptr< TTree > getObjectPtr (const string &name, const vector< ExpRun > &requestedRuns)
 We cheekily cast the TChain to TTree for the returned pointer so that the user never knows Hopefully this doesn't cause issues if people do low level stuff to the tree...
 
std::string getGranularityFromData () const
 Get the granularity of collected data.
 
void saveCalibration (TClonesArray *data, const std::string &name)
 Store DBArray payload with given name with default IOV.
 
void saveCalibration (TClonesArray *data, const std::string &name, const IntervalOfValidity &iov)
 Store DBArray with given name and custom IOV.
 
void saveCalibration (TObject *data)
 Store DB payload with default name and default IOV.
 
void saveCalibration (TObject *data, const IntervalOfValidity &iov)
 Store DB payload with default name and custom IOV.
 
void saveCalibration (TObject *data, const std::string &name)
 Store DB payload with given name with default IOV.
 
void saveCalibration (TObject *data, const std::string &name, const IntervalOfValidity &iov)
 Store DB payload with given name and custom IOV.
 
void updateDBObjPtrs (const unsigned int event, const int run, const int experiment)
 Updates any DBObjPtrs by calling update(event) for DBStore.
 
void setDescription (const std::string &description)
 Set algorithm description (in constructor)
 
void clearCalibrationData ()
 Clear calibration data.
 
Calibration::ExpRun getAllGranularityExpRun () const
 Returns the Exp,Run pair that means 'Everything'. Currently unused.
 
void resetInputJson ()
 Clears the m_inputJson member variable.
 
void resetOutputJson ()
 Clears the m_outputJson member variable.
 
template<class T>
void setOutputJsonValue (const std::string &key, const T &value)
 Set a key:value pair for the outputJson object, expected to used internally during calibrate()
 
template<class T>
const T getOutputJsonValue (const std::string &key) const
 Get a value using a key from the JSON output object, not sure why you would want to do this.
 
template<class T>
const T getInputJsonValue (const std::string &key) const
 Get an input JSON value using a key. The normal exceptions are raised when the key doesn't exist.
 
const nlohmann::json & getInputJsonObject () const
 Get the entire top level JSON object. We explicitly say this must be of object type so that we might pick.
 
bool inputJsonKeyExists (const std::string &key) const
 Test for a key in the input JSON object.
 

Protected Attributes

std::vector< Calibration::ExpRun > m_boundaries
 When using the boundaries functionality from isBoundaryRequired, this is used to store the boundaries. It is cleared when.
 

Private Member Functions

TTree * LambdaMassFit (std::shared_ptr< TTree > preselTree)
 Mass fit for Lambda->ppi.
 
std::unique_ptr< TList > LambdaHistogramming (TTree *inputTree)
 produce histograms for protons
 
TTree * DstarMassFit (std::shared_ptr< TTree > preselTree)
 Mass fit for D*->Dpi.
 
std::unique_ptr< TList > DstarHistogramming (TTree *inputTree)
 produce histograms for K/pi
 
std::unique_ptr< TList > GammaHistogramming (std::shared_ptr< TTree > preselTree)
 produce histograms for e
 
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
 
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 momentum bin is 1.
 
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 function that encodes dEdx resolution.
 
TH1F * PrepareProfile (TH2F *DataHistogram, TString NewName)
 Reimplement the Profile histogram calculation for a 2D histogram.
 
std::string getExpRunString (Calibration::ExpRun &expRun) const
 Gets the "exp.run" string repr. of (exp,run)
 
std::string getFullObjectPath (const std::string &name, Calibration::ExpRun expRun) const
 constructs the full TDirectory + Key name of an object in a TFile based on its name and exprun
 

Private Attributes

bool m_isMakePlots
 produce plots for monitoring
 
int m_numDEdxBins = 100
 the number of dEdx bins for the payloads
 
int m_numPBins = 69
 the number of momentum bins for the payloads
 
int m_numBGBins = 69
 the number of beta*gamma bins for the profile and fitting
 
double m_dedxCutoff = 5.e6
 the upper edge of the dEdx binning for the payloads
 
double m_dedxMaxPossible = 7.e6
 the approximate max possible value of dEdx
 
int m_MinEvtsPerTree
 number of events in TTree below which we don't try to fit
 
int m_NToGenerate
 the number of events to be generated in each momentum bin in the new payloads
 
bool m_CustomProfile = 1
 reimplement profile histogram calculation instead of the ROOT implementation?
 
bool m_UsePionBGFunctionForEverything
 Assume that the dEdx:betagamma trend is the same for all hadrons; use the pion trend as representative.
 
bool m_UseProtonBGFunctionForEverything
 Assume that the dEdx:betagamma trend is the same for all hadrons; use the proton trend as representative.
 
bool m_FixUnstableFitParameter
 In the dEdx:betagamma fit, there is one free parameter that makes fit convergence poor.
 
const double m_ElectronPDGMass = TDatabasePDG::Instance()->GetParticle(11)->Mass()
 PDG mass for the electron.
 
const double m_MuonPDGMass = TDatabasePDG::Instance()->GetParticle(13)->Mass()
 PDG mass for the muon.
 
const double m_PionPDGMass = TDatabasePDG::Instance()->GetParticle(211)->Mass()
 PDG mass for the charged pion.
 
const double m_KaonPDGMass = TDatabasePDG::Instance()->GetParticle(321)->Mass()
 PDG mass for the charged kaon.
 
const double m_ProtonPDGMass = TDatabasePDG::Instance()->GetParticle(2212)->Mass()
 PDG mass for the proton.
 
const double m_DeuteronPDGMass = TDatabasePDG::Instance()->GetParticle(1000010020)->Mass()
 PDG mass for the deuteron.
 
std::vector< std::string > m_inputFileNames
 List of input files to the Algorithm, will initially be user defined but then gets the wildcards expanded during execute()
 
std::map< Calibration::ExpRun, std::vector< std::string > > m_runsToInputFiles
 Map of Runs to input files. Gets filled when you call getRunRangeFromAllData, gets cleared when setting input files again.
 
std::string m_granularityOfData
 Granularity of input data. This only changes when the input files change so it isn't specific to an execution.
 
ExecutionData m_data
 Data specific to a SINGLE execution of the algorithm. Gets reset at the beginning of execution.
 
std::string m_description {""}
 Description of the algorithm.
 
std::string m_prefix {""}
 The name of the TDirectory the collector objects are contained within.
 
nlohmann::json m_jsonExecutionInput = nlohmann::json::object()
 Optional input JSON object used to make decisions about how to execute the algorithm code.
 
nlohmann::json m_jsonExecutionOutput = nlohmann::json::object()
 Optional output JSON object that can be set during the execution by the underlying algorithm code.
 

Static Private Attributes

static const Calibration::ExpRun m_allExpRun = make_pair(-1, -1)
 allExpRun
 

Detailed Description

Class implementing the SVD dEdx calibration algorithm.

Definition at line 28 of file SVDdEdxCalibrationAlgorithm.h.

Member Enumeration Documentation

◆ EResult

enum EResult
inherited

The result of calibration.

Enumerator
c_OK 

Finished successfully =0 in Python.

c_Iterate 

Needs iteration =1 in Python.

c_NotEnoughData 

Needs more data =2 in Python.

c_Failure 

Failed =3 in Python.

c_Undefined 

Not yet known (before execution) =4 in Python.

Definition at line 40 of file CalibrationAlgorithm.h.

40 {
41 c_OK,
42 c_Iterate,
43 c_NotEnoughData,
44 c_Failure,
45 c_Undefined
46 };

Constructor & Destructor Documentation

◆ SVDdEdxCalibrationAlgorithm()

Constructor.

Definition at line 41 of file SVDdEdxCalibrationAlgorithm.cc.

41 : CalibrationAlgorithm("SVDdEdxCollector"),
42 m_isMakePlots(true)
43{
44 setDescription("SVD dE/dx calibration algorithm");
45}
void setDescription(const std::string &description)
Set algorithm description (in constructor)
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
bool m_isMakePlots
produce plots for monitoring

◆ ~SVDdEdxCalibrationAlgorithm()

virtual ~SVDdEdxCalibrationAlgorithm ( )
inlinevirtual

Destructor.

Definition at line 40 of file SVDdEdxCalibrationAlgorithm.h.

40{}

Member Function Documentation

◆ boundaryFindingSetup()

virtual void boundaryFindingSetup ( std::vector< Calibration::ExpRun > ,
int  )
inlineprotectedvirtualinherited

If you need to make some changes to your algorithm class before 'findPayloadBoundaries' is run, make them in this function.

Reimplemented in PXDAnalyticGainCalibrationAlgorithm, PXDValidationAlgorithm, SVD3SampleCoGTimeCalibrationAlgorithm, SVD3SampleELSTimeCalibrationAlgorithm, SVDCoGTimeCalibrationAlgorithm, TestBoundarySettingAlgorithm, and TestCalibrationAlgorithm.

Definition at line 252 of file CalibrationAlgorithm.h.

252{};

◆ boundaryFindingTearDown()

virtual void boundaryFindingTearDown ( )
inlineprotectedvirtualinherited

Put your algorithm back into a state ready for normal execution if you need to.

Definition at line 257 of file CalibrationAlgorithm.h.

257{};

◆ calibrate()

CalibrationAlgorithm::EResult calibrate ( )
overrideprotectedvirtual

run algorithm on data

Implements CalibrationAlgorithm.

Definition at line 48 of file SVDdEdxCalibrationAlgorithm.cc.

49{
50 gROOT->SetBatch(true);
51
52 const auto exprun = getRunList()[0];
53 B2INFO("ExpRun used for calibration: " << exprun.first << " " << exprun.second);
54
55 auto payload = new Belle2::SVDdEdxPDFs();
56
57 // Get data objects
58 auto ttreeLambda = getObjectPtr<TTree>("Lambda");
59 auto ttreeDstar = getObjectPtr<TTree>("Dstar");
60 auto ttreeGamma = getObjectPtr<TTree>("Gamma");
61 auto ttreeGeneric = getObjectPtr<TTree>("Generic");
62
63 if ((ttreeLambda->GetEntries() < m_MinEvtsPerTree) || (ttreeDstar->GetEntries() < m_MinEvtsPerTree)
64 || (ttreeGamma->GetEntries() < m_MinEvtsPerTree)) {
65 B2WARNING("Not enough data for calibration.");
66 return c_NotEnoughData;
67 }
68
69 // call the calibration function
70 std::unique_ptr<TList> GeneratedList = GenerateNewHistograms(ttreeLambda, ttreeDstar, ttreeGamma, ttreeGeneric);
71
72 TH2F* histoE = (TH2F*) GeneratedList->FindObject("Electron2DHistogramNew");
73 TH2F* histoMu = (TH2F*) GeneratedList->FindObject("Muon2DHistogramNew");
74 TH2F* histoPi = (TH2F*) GeneratedList->FindObject("Pion2DHistogramNew");
75 TH2F* histoK = (TH2F*) GeneratedList->FindObject("Kaon2DHistogramNew");
76 TH2F* histoP = (TH2F*) GeneratedList->FindObject("Proton2DHistogramNew");
77 TH2F* histoDeut = (TH2F*) GeneratedList->FindObject("Deuteron2DHistogramNew");
78
79 std::vector<double> pbins = CreatePBinningScheme();
80 TH2F hEmpty("hEmpty", "A histogram returned if we cannot calibrate", m_numPBins, pbins.data(), m_numDEdxBins, 0, m_dedxCutoff);
81 for (int pbin = 0; pbin <= m_numPBins + 1; pbin++) {
82 for (int dedxbin = 0; dedxbin <= m_numDEdxBins + 1; dedxbin++) {
83 hEmpty.SetBinContent(pbin, dedxbin, 0.01);
84 };
85 }
86
87 B2INFO("Histograms are ready, proceed to creating the payload object...");
88 std::vector<TH2F*> hDedxPDFs(6);
89
90 std::array<std::string, 6> part = {"Electron", "Muon", "Pion", "Kaon", "Proton", "Deuteron"};
91
92 std::unique_ptr<TCanvas> candEdx(new TCanvas("candEdx", "SVD dEdx payloads", 1200, 700));
93 candEdx->Divide(3, 2);
94 gStyle->SetOptStat(11);
95
96 for (bool trunmean : {false, true}) {
97 for (int iPart = 0; iPart < 6; iPart++) {
98 if (iPart == 0 && trunmean) {
99 hDedxPDFs[iPart] = histoE;
100 hDedxPDFs[iPart]->SetName("hist_d1_11_trunc");
101 } else if (iPart == 1 && trunmean) {
102 hDedxPDFs[iPart] = histoMu;
103 hDedxPDFs[iPart]->SetName("hist_d1_13_trunc");
104 } else if (iPart == 2 && trunmean) {
105 hDedxPDFs[iPart] = histoPi;
106 hDedxPDFs[iPart]->SetName("hist_d1_211_trunc");
107 } else if (iPart == 3 && trunmean) {
108 hDedxPDFs[iPart] = histoK;
109 hDedxPDFs[iPart]->SetName("hist_d1_321_trunc");
110 } else if (iPart == 4 && trunmean) {
111 hDedxPDFs[iPart] = histoP;
112 hDedxPDFs[iPart]->SetName("hist_d1_2212_trunc");
113 } else if (iPart == 5 && trunmean) {
114 hDedxPDFs[iPart] = histoDeut;
115 hDedxPDFs[iPart]->SetName("hist_d1_1000010020_trunc");
116 } else if (iPart == 0 && !trunmean) {
117 hDedxPDFs[iPart] = &hEmpty;
118 hDedxPDFs[iPart]->SetName("hist_d1_11");
119 } else if (iPart == 1 && !trunmean) {
120 hDedxPDFs[iPart] = &hEmpty;
121 hDedxPDFs[iPart]->SetName("hist_d1_13");
122 } else if (iPart == 2 && !trunmean) {
123 hDedxPDFs[iPart] = &hEmpty;
124 hDedxPDFs[iPart]->SetName("hist_d1_211");
125 } else if (iPart == 3 && !trunmean) {
126 hDedxPDFs[iPart] = &hEmpty;
127 hDedxPDFs[iPart]->SetName("hist_d1_321");
128 } else if (iPart == 4 && !trunmean) {
129 hDedxPDFs[iPart] = &hEmpty;
130 hDedxPDFs[iPart]->SetName("hist_d1_2212");
131 } else if (iPart == 5 && !trunmean) {
132 hDedxPDFs[iPart] = &hEmpty;
133 hDedxPDFs[iPart]->SetName("hist_d1_1000010020");
134 } else
135 hDedxPDFs[iPart] = &hEmpty;
136 payload->setPDF(*hDedxPDFs[iPart], iPart, trunmean);
137
138 candEdx->cd(iPart + 1);
139 hDedxPDFs[iPart]->SetTitle(Form("%s; p(GeV/c) of %s; dE/dx", hDedxPDFs[iPart]->GetTitle(), part[iPart].data()));
140 hDedxPDFs[iPart]->DrawCopy("colz");
141 }
142
143 if (m_isMakePlots) {
144 candEdx->SaveAs("PlotsSVDdEdxPDFs_wTruncMean.pdf");
145 std::unique_ptr<TList> l(new TList());
146 for (int iPart = 0; iPart < 6; iPart++) {
147 l->Add(hDedxPDFs[iPart]);
148 }
149
150 TFile SVDdEdxPDFsPlotFile("PlotsSVDdEdxPDFs_wTruncMean.root", "RECREATE");
151 l->Write("histlist", TObject::kSingleKey);
152 SVDdEdxPDFsPlotFile.Close();
153 }
154
155 // candEdx->SetTitle(Form("Likehood dist. of charged particles from %s, trunmean = %s", idet.data(), check.str().data()));
156 }
157
158 saveCalibration(payload, "SVDdEdxPDFs");
159 B2INFO("SVD dE/dx calibration done!");
160
161 return c_OK;
162}
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
@ c_OK
Finished successfully =0 in Python.
@ c_NotEnoughData
Needs more data =2 in Python.
int m_numPBins
the number of momentum bins for the payloads
std::vector< double > CreatePBinningScheme()
build the binning scheme for the momentum
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
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
double m_dedxCutoff
the upper edge of the dEdx binning for the payloads
std::shared_ptr< T > getObjectPtr(const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
Get calibration data object by name and list of runs, the Merge function will be called to generate t...

◆ checkPyExpRun()

bool checkPyExpRun ( PyObject * pyObj)
inherited

Checks that a PyObject can be successfully converted to an ExpRun type.

Checks if the PyObject can be converted to ExpRun.

Definition at line 28 of file CalibrationAlgorithm.cc.

29{
30 // Is it a sequence?
31 if (PySequence_Check(pyObj)) {
32 Py_ssize_t nObj = PySequence_Length(pyObj);
33 // Does it have 2 objects in it?
34 if (nObj != 2) {
35 B2DEBUG(29, "ExpRun was a Python sequence which didn't have exactly 2 entries!");
36 return false;
37 }
38 PyObject* item1, *item2;
39 item1 = PySequence_GetItem(pyObj, 0);
40 item2 = PySequence_GetItem(pyObj, 1);
41 // Did the GetItem work?
42 if ((item1 == NULL) || (item2 == NULL)) {
43 B2DEBUG(29, "A PyObject pointer was NULL in the sequence");
44 return false;
45 }
46 // Are they longs?
47 if (PyLong_Check(item1) && PyLong_Check(item2)) {
48 long value1, value2;
49 value1 = PyLong_AsLong(item1);
50 value2 = PyLong_AsLong(item2);
51 if (((value1 == -1) || (value2 == -1)) && PyErr_Occurred()) {
52 B2DEBUG(29, "An error occurred while converting the PyLong to long");
53 return false;
54 }
55 } else {
56 B2DEBUG(29, "One or more of the PyObjects in the ExpRun wasn't a long");
57 return false;
58 }
59 // Make sure to kill off the reference GetItem gave us responsibility for
60 Py_DECREF(item1);
61 Py_DECREF(item2);
62 } else {
63 B2DEBUG(29, "ExpRun was not a Python sequence.");
64 return false;
65 }
66 return true;
67}

◆ clearCalibrationData()

void clearCalibrationData ( )
inlineprotectedinherited

Clear calibration data.

Definition at line 324 of file CalibrationAlgorithm.h.

324{m_data.clearCalibrationData();}

◆ commit() [1/2]

bool commit ( )
inherited

Submit constants from last calibration into database.

Definition at line 302 of file CalibrationAlgorithm.cc.

303{
304 if (getPayloads().empty())
305 return false;
306 list<Database::DBImportQuery> payloads = getPayloads();
307 B2INFO("Committing " << payloads.size() << " payloads to database.");
308 return Database::Instance().storeData(payloads);
309}
std::list< Database::DBImportQuery > & getPayloads()
Get constants (in TObjects) for database update from last execution.
static Database & Instance()
Instance of a singleton Database.
Definition Database.cc:41
bool storeData(const std::string &name, TObject *object, const IntervalOfValidity &iov)
Store an object in the database.
Definition Database.cc:140

◆ commit() [2/2]

bool commit ( std::list< Database::DBImportQuery > payloads)
inherited

Submit constants from a (potentially previous) set of payloads.

Definition at line 311 of file CalibrationAlgorithm.cc.

312{
313 if (payloads.empty())
314 return false;
315 return Database::Instance().storeData(payloads);
316}

◆ convertPyExpRun()

ExpRun convertPyExpRun ( PyObject * pyObj)
inherited

Performs the conversion of PyObject to ExpRun.

Converts the PyObject to an ExpRun. We've preoviously checked the object so this assumes a lot about the PyObject.

Definition at line 70 of file CalibrationAlgorithm.cc.

71{
72 ExpRun expRun;
73 PyObject* itemExp, *itemRun;
74 itemExp = PySequence_GetItem(pyObj, 0);
75 itemRun = PySequence_GetItem(pyObj, 1);
76 expRun.first = PyLong_AsLong(itemExp);
77 Py_DECREF(itemExp);
78 expRun.second = PyLong_AsLong(itemRun);
79 Py_DECREF(itemRun);
80 return expRun;
81}

◆ CreatePBinningScheme()

std::vector< double > CreatePBinningScheme ( )
inlineprivate

build the binning scheme for the momentum

Definition at line 140 of file SVDdEdxCalibrationAlgorithm.h.

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 }

◆ DstarHistogramming()

std::unique_ptr< TList > DstarHistogramming ( TTree * inputTree)
private

produce histograms for K/pi

Definition at line 523 of file SVDdEdxCalibrationAlgorithm.cc.

524{
525 gROOT->SetBatch(true);
526 inputTree->SetEstimate(-1);
527 std::vector<double> pbins = CreatePBinningScheme();
528
529 TH2F* hDstarKMomentum = new TH2F("hist_d1_321_truncMomentum", "hist_d1_321_trunc;Momentum [GeV/c];dEdx [arb. units]", m_numPBins,
530 pbins.data(),
532 // the pion payload
533 TH2F* hDstarPiMomentum = new TH2F("hist_d1_211_truncMomentum", "hist_d1_211_trunc;Momentum [GeV/c];dEdx [arb. units]", m_numPBins,
534 pbins.data(),
536
537 inputTree->Draw("KaonSVDdEdx:KaonSVDdEdxTrackMomentum>>hist_d1_321_truncMomentum", "nSignalDstar_sw * (KaonSVDdEdx>0)", "goff");
538 // the pion one will be built from both pions in the Dstar decay tree
539 TH2F* hDstarPiPart1Momentum = (TH2F*)hDstarPiMomentum->Clone("hist_d1_211_truncPart1Momentum");
540 TH2F* hDstarPiPart2Momentum = (TH2F*)hDstarPiMomentum->Clone("hist_d1_211_truncPart2Momentum");
541
542 inputTree->Draw("PionDSVDdEdx:PionDSVDdEdxTrackMomentum>>hist_d1_211_truncPart1Momentum", "nSignalDstar_sw * (PionDSVDdEdx>0)",
543 "goff");
544 inputTree->Draw("SlowPionSVDdEdx:SlowPionSVDdEdxTrackMomentum>>hist_d1_211_truncPart2Momentum",
545 "nSignalDstar_sw * (SlowPionSVDdEdx>0)",
546 "goff");
547 hDstarPiMomentum->Add(hDstarPiPart1Momentum);
548 hDstarPiMomentum->Add(hDstarPiPart2Momentum);
549
550 inputTree->Draw(Form("KaonSVDdEdxTrackMomentum/%f", m_KaonPDGMass), "", "goff",
551 ((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins);
552 double* KaonMomentumDataset = inputTree->GetV1();
553 TKDTreeBinning* kdBinsK = new TKDTreeBinning(((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, KaonMomentumDataset,
555 const double* binsMinEdgesKOriginal = kdBinsK->SortOneDimBinEdges();
556 double* binsMinEdgesK = const_cast<double*>(binsMinEdgesKOriginal);
557 binsMinEdgesK[0] = 0.1;
558 binsMinEdgesK[m_numBGBins + 1] = 50.;
559
560// get a distribution that contains both pions to get a typical kinematics for the binning scheme
561 inputTree->Draw(Form("SlowPionSVDdEdxTrackMomentum/%f * (event %% 2 == 0) + PionDSVDdEdxTrackMomentum/%f * (event %% 2 ==1)",
562 m_PionPDGMass, m_PionPDGMass), "", "goff",
563 ((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins);
564 double* PionMomentumDataset = inputTree->GetV1();
565
566 TKDTreeBinning* kdBinsPi = new TKDTreeBinning(((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, PionMomentumDataset,
568 const double* binsMinEdgesPiOriginal = kdBinsPi->SortOneDimBinEdges();
569 double* binsMinEdgesPi = const_cast<double*>(binsMinEdgesPiOriginal);
570 binsMinEdgesPi[0] = 0.1;
571 binsMinEdgesPi[m_numBGBins + 1] = 50.;
572
573 TH2F* hDstarKBetaGamma = new TH2F("hist_d1_321_truncBetaGamma", "hist_d1_321_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
575 binsMinEdgesK,
577 // the pion payload
578 TH2F* hDstarPiBetaGamma = new TH2F("hist_d1_211_truncBetaGamma", "hist_d1_211_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
580 binsMinEdgesPi,
582
583 inputTree->Draw(Form("KaonSVDdEdx:KaonSVDdEdxTrackMomentum/%f>>hist_d1_321_truncBetaGamma", m_KaonPDGMass),
584 "nSignalDstar_sw * (KaonSVDdEdx>0)", "goff");
585 // the pion one will be built from both pions in the Dstar decay tree
586 TH2F* hDstarPiPart1BetaGamma = (TH2F*)hDstarPiBetaGamma->Clone("hist_d1_211_truncPart1BetaGamma");
587 TH2F* hDstarPiPart2BetaGamma = (TH2F*)hDstarPiBetaGamma->Clone("hist_d1_211_truncPart2BetaGamma");
588
589 inputTree->Draw(Form("PionDSVDdEdx:PionDSVDdEdxTrackMomentum/%f>>hist_d1_211_truncPart1BetaGamma", m_PionPDGMass),
590 "nSignalDstar_sw * (PionDSVDdEdx>0)",
591 "goff");
592 inputTree->Draw(Form("SlowPionSVDdEdx:SlowPionSVDdEdxTrackMomentum/%f>>hist_d1_211_truncPart2BetaGamma", m_PionPDGMass),
593 "nSignalDstar_sw * (SlowPionSVDdEdx>0)", "goff");
594 hDstarPiBetaGamma->Add(hDstarPiPart1BetaGamma);
595 hDstarPiBetaGamma->Add(hDstarPiPart2BetaGamma);
596
597
598
599 // produce the 1D profiles
600
601
602 TH1F* PionProfileMomentum = (TH1F*)hDstarPiMomentum->ProfileX("PionProfileMomentum");
603 PionProfileMomentum->SetTitle("PionProfile");
604 PionProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
605 PionProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
606 PionProfileMomentum->GetYaxis()->SetTitle("dE/dx");
607 PionProfileMomentum->SetLineColor(kRed);
608
609 TH1F* PionProfileBetaGamma = (TH1F*)hDstarPiBetaGamma->ProfileX("PionProfileBetaGamma");
610 if (m_CustomProfile) {
611 PionProfileBetaGamma = PrepareProfile(hDstarPiBetaGamma, "PionProfileBetaGamma");
612 }
613 PionProfileBetaGamma->SetTitle("PionProfile");
614 PionProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
615 PionProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
616 PionProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
617 PionProfileBetaGamma->SetLineColor(kRed);
618
619
620 TH1F* KaonProfileMomentum = (TH1F*)hDstarKMomentum->ProfileX("KaonProfileMomentum");
621 KaonProfileMomentum->SetTitle("KaonProfile");
622 KaonProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
623 KaonProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
624 KaonProfileMomentum->GetYaxis()->SetTitle("dE/dx");
625 KaonProfileMomentum->SetLineColor(kRed);
626
627
628 TH1F* KaonProfileBetaGamma = (TH1F*)hDstarKBetaGamma->ProfileX("KaonProfileBetaGamma");
629 if (m_CustomProfile) {
630 KaonProfileBetaGamma = PrepareProfile(hDstarKBetaGamma, "KaonProfileBetaGamma");
631 }
632 KaonProfileBetaGamma->SetTitle("KaonProfile");
633
634 KaonProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
635 KaonProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
636 KaonProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
637 KaonProfileBetaGamma->SetLineColor(kRed);
638
639
640
641 // normalisation
642 hDstarKMomentum->Sumw2();
643 hDstarKMomentum = Normalise2DHisto(hDstarKMomentum);
644
645 hDstarPiMomentum->Sumw2();
646 hDstarPiMomentum = Normalise2DHisto(hDstarPiMomentum);
647
648 std::unique_ptr<TList> histList(new TList);
649 histList->Add(KaonProfileMomentum);
650 histList->Add(KaonProfileBetaGamma);
651 histList->Add(hDstarKMomentum);
652
653 histList->Add(PionProfileMomentum);
654 histList->Add(PionProfileBetaGamma);
655 histList->Add(hDstarPiMomentum);
656
657 if (m_isMakePlots) {
658 TFile DstarHistogrammingPlotFile("SVDdEdxCalibrationDstarHistogramming.root", "RECREATE");
659 histList->Write();
660 DstarHistogrammingPlotFile.Close();
661 }
662
663 return histList;
664}
TH1F * PrepareProfile(TH2F *DataHistogram, TString NewName)
Reimplement the Profile histogram calculation for a 2D histogram.
TH2F * Normalise2DHisto(TH2F *HistoToNormalise)
Normalise a given dEdx:momentum histogram in each momentum bin, so that sum of entries in each moment...
int m_numBGBins
the number of beta*gamma bins for the profile and fitting
double m_dedxMaxPossible
the approximate max possible value of dEdx
const double m_KaonPDGMass
PDG mass for the charged kaon.
bool m_CustomProfile
reimplement profile histogram calculation instead of the ROOT implementation?
const double m_PionPDGMass
PDG mass for the charged pion.

◆ DstarMassFit()

TTree * DstarMassFit ( std::shared_ptr< TTree > preselTree)
private

Mass fit for D*->Dpi.

Definition at line 387 of file SVDdEdxCalibrationAlgorithm.cc.

388{
389 B2INFO("Configuring the Dstar fit...");
390 gROOT->SetBatch(true);
391 RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
392
393 RooRealVar deltaM("deltaM", "m(D*)-m(D^{0})", 0.139545, 0.151, "GeV/c^{2}");
394
395 RooRealVar KaonMomentum("KaonMomentum", "momentum for Kaon (GeV)", -1.e8, 1.e8);
396 RooRealVar KaonSVDdEdxTrackMomentum("KaonSVDdEdxTrackMomentum", "momentum for Kaon (GeV), from the track", -1.e8, 1.e8);
397 RooRealVar KaonSVDdEdx("KaonSVDdEdx", "", -1.e8, 1.e8);
398 RooRealVar PionDMomentum("PionDMomentum", "momentum for pion (GeV)", -1.e8, 1.e8);
399 RooRealVar PionDSVDdEdxTrackMomentum("PionDSVDdEdxTrackMomentum", "momentum for pion (GeV), from the track", -1.e8, 1.e8);
400 RooRealVar PionDSVDdEdx("PionDSVDdEdx", "", -1.e8, 1.e8);
401 RooRealVar SlowPionMomentum("SlowPionMomentum", "momentum for slow pion (GeV)", -1.e8, 1.e8);
402 RooRealVar SlowPionSVDdEdxTrackMomentum("SlowPionSVDdEdxTrackMomentum", "momentum for slow pion (GeV), from the track", -1.e8,
403 1.e8);
404 RooRealVar SlowPionSVDdEdx("SlowPionSVDdEdx", "", -1.e8, 1.e8);
405
406 RooRealVar exp("exp", "experiment number", 0, 1.e5);
407 RooRealVar run("run", "run number", 0, 1.e8);
408 RooRealVar event("event", "event number", 0, 1.e10);
409
410 auto variables = new RooArgSet();
411 variables->add(deltaM);
412 variables->add(KaonMomentum);
413 variables->add(KaonSVDdEdxTrackMomentum);
414 variables->add(KaonSVDdEdx);
415 variables->add(PionDMomentum);
416 variables->add(PionDSVDdEdxTrackMomentum);
417 variables->add(PionDSVDdEdx);
418 variables->add(SlowPionMomentum);
419 variables->add(SlowPionSVDdEdxTrackMomentum);
420 variables->add(SlowPionSVDdEdx);
421 variables->add(exp);
422 variables->add(run);
423 variables->add(event);
424
425 RooDataSet* DstarDataset = new RooDataSet("DstarDataset", "DstarDataset", *variables, Import(*preselTree));
426
427 if (DstarDataset->sumEntries() == 0) {
428 B2FATAL("The Dstar dataset is empty, stopping here");
429 }
430
431 RooPlot* DstarFitFrame = DstarDataset->plotOn(deltaM.frame());
432
433 RooRealVar GaussMean("GaussMean", "GaussMean", 0.145, 0.140, 0.150);
434 RooRealVar GaussSigma1("GaussSigma1", "GaussSigma1", 0.01, 1.e-4, 1.0);
435 RooGaussian DstarGauss1("DstarGauss1", "DstarGauss1", deltaM, GaussMean, GaussSigma1);
436 RooRealVar GaussSigma2("GaussSigma2", "GaussSigma2", 0.001, 1.e-4, 1.0);
437 RooGaussian DstarGauss2("DstarGauss2", "DstarGauss2", deltaM, GaussMean, GaussSigma2);
438 RooRealVar fracGaussYield("fracGaussYield", "Fraction of two Gaussians", 0.75, 0.0, 1.0);
439 RooAddPdf DstarSignalPDF("DstarSignalPDF", "DstarGauss1+DstarGauss2", RooArgList(DstarGauss1, DstarGauss2), fracGaussYield);
440
441 RooRealVar dm0Bkg("dm0Bkg", "dm0", 0.13957018, 0.130, 0.140);
442 RooRealVar aBkg("aBkg", "a", -0.0784, -0.08, 3.0);
443 RooRealVar bBkg("bBkg", "b", -0.444713, -0.5, 0.4);
444 RooRealVar cBkg("cBkg", "c", 0.3);
445 RooDstD0BG DstarBkgPDF("DstarBkgPDF", "DstarBkgPDF", deltaM, dm0Bkg, cBkg, aBkg, bBkg);
446 RooRealVar nSignalDstar("nSignalDstar", "signal yield", 0.5 * preselTree->GetEntries(), 0, preselTree->GetEntries());
447 RooRealVar nBkgDstar("nBkgDstar", "background yield", 0.5 * preselTree->GetEntries(), 0, preselTree->GetEntries());
448 RooAddPdf totalPDFDstar("totalPDFDstar", "totalPDFDstar pdf", RooArgList(DstarSignalPDF, DstarBkgPDF),
449 RooArgList(nSignalDstar, nBkgDstar));
450
451 B2INFO("Dstar: Start fitting...");
452 RooFitResult* DstarFitResult = totalPDFDstar.fitTo(*DstarDataset, Save(kTRUE), PrintLevel(-1));
453
454 int status = DstarFitResult->status();
455 int covqual = DstarFitResult->covQual();
456 double diff = nSignalDstar.getValV() + nBkgDstar.getValV() - DstarDataset->sumEntries();
457
458 B2INFO("Dstar: Fit status: " << status << "; covariance quality: " << covqual);
459 // if the fit is not healthy, try again once before giving up, with a slightly different setup:
460 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalDstar.getError() < sqrt(nSignalDstar.getValV()))
461 || (nSignalDstar.getError() > (nSignalDstar.getValV()))) {
462
463 DstarFitResult = totalPDFDstar.fitTo(*DstarDataset, Save(), Strategy(2), Offset(1));
464 status = DstarFitResult->status();
465 covqual = DstarFitResult->covQual();
466 diff = nSignalDstar.getValV() + nBkgDstar.getValV() - DstarDataset->sumEntries();
467 B2INFO("Dstar: Updated fit status: " << status << "; covariance quality: " << covqual);
468 }
469
470 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalDstar.getError() < sqrt(nSignalDstar.getValV()))
471 || (nSignalDstar.getError() > (nSignalDstar.getValV()))) {
472 B2WARNING("Dstar: Fit problem: fit status " << status << "; sum of component yields minus the dataset yield is " << diff <<
473 "; signal yield is " << nSignalDstar.getValV() << ", while its uncertainty is " << nSignalDstar.getError());
474 }
475 if (covqual < 2) {
476 B2INFO("Dstar: Fit warning: covariance quality " << covqual);
477 }
478
479 totalPDFDstar.plotOn(DstarFitFrame, LineColor(TColor::GetColor("#4575b4")));
480
481 double chisquare = DstarFitFrame->chiSquare();
482 B2INFO("Dstar: Fit chi2 = " << chisquare);
483 totalPDFDstar.paramOn(DstarFitFrame, Layout(0.63, 0.96, 0.93), Format("NEU", AutoPrecision(2)));
484 DstarFitFrame->getAttText()->SetTextSize(0.03);
485
486 totalPDFDstar.plotOn(DstarFitFrame, Components("DstarSignalPDF"), LineColor(TColor::GetColor("#d73027")));
487 totalPDFDstar.plotOn(DstarFitFrame, Components("DstarBkgPDF"), LineColor(TColor::GetColor("#fc8d59")));
488 totalPDFDstar.plotOn(DstarFitFrame, LineColor(TColor::GetColor("#4575b4")));
489
490 DstarFitFrame->GetXaxis()->SetTitle("#Deltam [GeV/c^{2}]");
491 if (m_isMakePlots) {
492 std::unique_ptr<TCanvas> canvDstar(new TCanvas("canvDstar", "canvDstar"));
493 canvDstar->cd();
494
495 DstarFitFrame->Draw();
496
497 canvDstar->Print("SVDdEdxCalibrationFitDstar.pdf");
498 TFile DstarFitPlotFile("SVDdEdxCalibrationDstarFitPlotFile.root", "RECREATE");
499 canvDstar->Write();
500 DstarFitPlotFile.Close();
501 }
502
504
505 RooStats::SPlot* sPlotDatasetDstar = new RooStats::SPlot("sData", "An SPlot", *DstarDataset, &totalPDFDstar,
506 RooArgList(nSignalDstar, nBkgDstar));
507
508 for (int iEvt = 0; iEvt < 5; iEvt++) {
509 if (TMath::Abs(sPlotDatasetDstar->GetSWeight(iEvt, "nSignalDstar") + sPlotDatasetDstar->GetSWeight(iEvt, "nBkgDstar") - 1) > 5.e-3)
510 B2FATAL("Dstar: sPlot error: sum of weights not equal to 1");
511 }
512
513 RooDataSet* DstarDatasetSWeighted = new RooDataSet(DstarDataset->GetName(), DstarDataset->GetTitle(), DstarDataset,
514 *DstarDataset->get());
515
516 TTree* treeDstarSWeighted = DstarDatasetSWeighted->GetClonedTree();
517 treeDstarSWeighted->SetName("treeDstarSWeighted");
518
519 B2INFO("Dstar: sPlot done. Proceed to histogramming");
520 return treeDstarSWeighted;
521}
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28

◆ dumpOutputJson()

const std::string dumpOutputJson ( ) const
inlineinherited

Dump the JSON string of the output JSON object.

Definition at line 223 of file CalibrationAlgorithm.h.

223{return m_jsonExecutionOutput.dump();}

◆ execute() [1/2]

CalibrationAlgorithm::EResult execute ( PyObject * runs,
int iteration = 0,
IntervalOfValidity iov = IntervalOfValidity() )
inherited

Runs calibration over Python list of runs. Converts to C++ and then calls the other execute() function.

Definition at line 83 of file CalibrationAlgorithm.cc.

84{
85 B2DEBUG(29, "Running execute() using Python Object as input argument");
86 // Reset the execution specific data in case the algorithm was previously called
87 m_data.reset();
88 m_data.setIteration(iteration);
89 vector<ExpRun> vecRuns;
90 // Is it a list?
91 if (PySequence_Check(runs)) {
92 boost::python::handle<> handle(boost::python::borrowed(runs));
93 boost::python::list listRuns(handle);
94
95 int nList = boost::python::len(listRuns);
96 for (int iList = 0; iList < nList; ++iList) {
97 boost::python::object pyExpRun(listRuns[iList]);
98 if (!checkPyExpRun(pyExpRun.ptr())) {
99 B2ERROR("Received Python ExpRuns couldn't be converted to C++");
100 m_data.setResult(c_Failure);
101 return c_Failure;
102 } else {
103 vecRuns.push_back(convertPyExpRun(pyExpRun.ptr()));
104 }
105 }
106 } else {
107 B2ERROR("Tried to set the input runs but we didn't receive a Python sequence object (list,tuple).");
108 m_data.setResult(c_Failure);
109 return c_Failure;
110 }
111 return execute(vecRuns, iteration, iov);
112}
bool checkPyExpRun(PyObject *pyObj)
Checks that a PyObject can be successfully converted to an ExpRun type.
EResult execute(std::vector< Calibration::ExpRun > runs={}, int iteration=0, IntervalOfValidity iov=IntervalOfValidity())
Runs calibration over vector of runs for a given iteration.
Calibration::ExpRun convertPyExpRun(PyObject *pyObj)
Performs the conversion of PyObject to ExpRun.
ExecutionData m_data
Data specific to a SINGLE execution of the algorithm. Gets reset at the beginning of execution.

◆ execute() [2/2]

CalibrationAlgorithm::EResult execute ( std::vector< Calibration::ExpRun > runs = {},
int iteration = 0,
IntervalOfValidity iov = IntervalOfValidity() )
inherited

Runs calibration over vector of runs for a given iteration.

You can also specify the IoV to save the database payload as. By default the Algorithm will create an IoV from your requested ExpRuns, or from the overall ExpRuns of the input data if you haven't specified ExpRuns in this function.

No checks are performed to make sure that a IoV you specify matches the data you ran over, it simply labels the IoV to commit to the database later.

Definition at line 114 of file CalibrationAlgorithm.cc.

115{
116 // Check if we are calling this function directly and need to reset, or through Python where it was already done.
117 if (m_data.getResult() != c_Undefined) {
118 m_data.reset();
119 m_data.setIteration(iteration);
120 }
121
122 if (m_inputFileNames.empty()) {
123 B2ERROR("There aren't any input files set. Please use CalibrationAlgorithm::setInputFiles()");
124 m_data.setResult(c_Failure);
125 return c_Failure;
126 }
127
128 // Did we receive runs to execute over explicitly?
129 if (!(runs.empty())) {
130 for (auto expRun : runs) {
131 B2DEBUG(29, "ExpRun requested = (" << expRun.first << ", " << expRun.second << ")");
132 }
133 // We've asked explicitly for certain runs, but we should check if the data granularity is 'run'
134 if (strcmp(getGranularity().c_str(), "all") == 0) {
135 B2ERROR(("The data is collected with granularity=all (exp=-1,run=-1), but you seem to request calibration for specific runs."
136 " We'll continue but using ALL the input data given instead of the specific runs requested."));
137 }
138 } else {
139 // If no runs are provided, infer the runs from all collected data
140 runs = getRunListFromAllData();
141 // Let's check that we have some now
142 if (runs.empty()) {
143 B2ERROR("No collected data in input files.");
144 m_data.setResult(c_Failure);
145 return c_Failure;
146 }
147 for (auto expRun : runs) {
148 B2DEBUG(29, "ExpRun requested = (" << expRun.first << ", " << expRun.second << ")");
149 }
150 }
151
152 m_data.setRequestedRuns(runs);
153 if (iov.empty()) {
154 // If no user specified IoV we use the IoV from the executed run list
155 iov = IntervalOfValidity(runs[0].first, runs[0].second, runs[runs.size() - 1].first, runs[runs.size() - 1].second);
156 }
157 m_data.setRequestedIov(iov);
158 // After here, the getObject<...>(...) helpers start to work
159
161 m_data.setResult(result);
162 return result;
163}
std::vector< Calibration::ExpRun > getRunListFromAllData() const
Get the complete list of runs from inspection of collected data.
std::vector< std::string > m_inputFileNames
List of input files to the Algorithm, will initially be user defined but then gets the wildcards expa...
EResult
The result of calibration.
@ c_Undefined
Not yet known (before execution) =4 in Python.
virtual EResult calibrate()=0
Run algo on data - pure virtual: needs to be implemented.
std::string getGranularity() const
Get the granularity of collected data.

◆ fillRunToInputFilesMap()

void fillRunToInputFilesMap ( )
inherited

Fill the mapping of ExpRun -> Files.

Definition at line 330 of file CalibrationAlgorithm.cc.

331{
332 m_runsToInputFiles.clear();
333 // Save TDirectory to change back at the end
334 TDirectory* dir = gDirectory;
335 RunRange* runRange;
336 // Construct the TDirectory name where we expect our objects to be
337 string runRangeObjName(getPrefix() + "/" + RUN_RANGE_OBJ_NAME);
338 for (const auto& fileName : m_inputFileNames) {
339 //Open TFile to get the objects
340 unique_ptr<TFile> f;
341 f.reset(TFile::Open(fileName.c_str(), "READ"));
342 runRange = dynamic_cast<RunRange*>(f->Get(runRangeObjName.c_str()));
343 if (runRange) {
344 // Insert or extend the run -> file mapping for this ExpRun
345 auto expRuns = runRange->getExpRunSet();
346 for (const auto& expRun : expRuns) {
347 auto runFiles = m_runsToInputFiles.find(expRun);
348 if (runFiles != m_runsToInputFiles.end()) {
349 (runFiles->second).push_back(fileName);
350 } else {
351 m_runsToInputFiles.insert(std::make_pair(expRun, std::vector<std::string> {fileName}));
352 }
353 }
354 } else {
355 B2WARNING("Missing a RunRange object for file: " << fileName);
356 }
357 }
358 dir->cd();
359}
std::string getPrefix() const
Get the prefix used for getting calibration data.
std::map< Calibration::ExpRun, std::vector< std::string > > m_runsToInputFiles
Map of Runs to input files. Gets filled when you call getRunRangeFromAllData, gets cleared when setti...
const std::set< Calibration::ExpRun > & getExpRunSet()
Get access to the stored set.
Definition RunRange.h:64

◆ findPayloadBoundaries()

const std::vector< ExpRun > findPayloadBoundaries ( std::vector< Calibration::ExpRun > runs,
int iteration = 0 )
inherited

Used to discover the ExpRun boundaries that you want the Python CAF to execute on. This is optional and only used in some.

Definition at line 520 of file CalibrationAlgorithm.cc.

521{
522 m_boundaries.clear();
523 if (m_inputFileNames.empty()) {
524 B2ERROR("There aren't any input files set. Please use CalibrationAlgorithm::setInputFiles()");
525 return m_boundaries;
526 }
527 // Reset the internal execution data just in case something is hanging around
528 m_data.reset();
529 if (runs.empty()) {
530 // Want to loop over all runs we could possibly know about
531 runs = getRunListFromAllData();
532 }
533 // Let's check that we have some now
534 if (runs.empty()) {
535 B2ERROR("No collected data in input files.");
536 return m_boundaries;
537 }
538 // In order to find run boundaries we must have collected with data granularity == 'run'
539 if (strcmp(getGranularity().c_str(), "all") == 0) {
540 B2ERROR("The data is collected with granularity='all' (exp=-1,run=-1), and we can't use that to find run boundaries.");
541 return m_boundaries;
542 }
543 m_data.setIteration(iteration);
544 // User defined setup function
545 boundaryFindingSetup(runs, iteration);
546 std::vector<ExpRun> runList;
547 // Loop over run list and call derived class "isBoundaryRequired" member function
548 for (auto currentRun : runs) {
549 runList.push_back(currentRun);
550 m_data.setRequestedRuns(runList);
551 // After here, the getObject<...>(...) helpers start to work
552 if (isBoundaryRequired(currentRun)) {
553 m_boundaries.push_back(currentRun);
554 }
555 // Only want run-by-run
556 runList.clear();
557 // Don't want memory hanging around
558 m_data.clearCalibrationData();
559 }
560 m_data.reset();
562 return m_boundaries;
563}
std::vector< Calibration::ExpRun > m_boundaries
When using the boundaries functionality from isBoundaryRequired, this is used to store the boundaries...
virtual void boundaryFindingTearDown()
Put your algorithm back into a state ready for normal execution if you need to.
virtual void boundaryFindingSetup(std::vector< Calibration::ExpRun >, int)
If you need to make some changes to your algorithm class before 'findPayloadBoundaries' is run,...
virtual bool isBoundaryRequired(const Calibration::ExpRun &)
Given the current collector data, make a decision about whether or not this run should be the start o...

◆ GammaHistogramming()

std::unique_ptr< TList > GammaHistogramming ( std::shared_ptr< TTree > preselTree)
private

produce histograms for e

Definition at line 667 of file SVDdEdxCalibrationAlgorithm.cc.

668{
669 B2INFO("Histogramming the converted photon selection...");
670 gROOT->SetBatch(true);
671
672
673 if (preselTree->GetEntries() == 0) {
674 B2FATAL("The Gamma tree is empty, stopping here");
675 }
676 preselTree->SetEstimate(-1);
677 std::vector<double> pbins = CreatePBinningScheme();
678
679
680 TH2F* hGammaEMomentum = new TH2F("hist_d1_11_truncMomentum", "hist_d1_11_trunc;Momentum [GeV/c];dEdx [arb. units]", m_numPBins,
681 pbins.data(), m_numDEdxBins, 0, m_dedxCutoff);
682
683 TH2F* hGammaEPart1Momentum = (TH2F*)hGammaEMomentum->Clone("hist_d1_11_truncPart1Momentum");
684 TH2F* hGammaEPart2Momentum = (TH2F*)hGammaEMomentum->Clone("hist_d1_11_truncPart2Momentum");
685
686 preselTree->Draw("FirstElectronSVDdEdx:FirstElectronSVDdEdxTrackMomentum>>hist_d1_11_truncPart1Momentum",
687 "FirstElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2", "goff");
688 preselTree->Draw("SecondElectronSVDdEdx:SecondElectronSVDdEdxTrackMomentum>>hist_d1_11_truncPart2Momentum",
689 "SecondElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2", "goff");
690 hGammaEMomentum->Add(hGammaEPart1Momentum);
691 hGammaEMomentum->Add(hGammaEPart2Momentum);
692
693// get a distribution that contains both pions to get a typical kinematics for the binning scheme
694 preselTree->Draw(
695 Form("FirstElectronSVDdEdxTrackMomentum/%f* (event %% 2==0) + SecondElectronSVDdEdxTrackMomentum/%f* (event %% 2==1)",
697 "", "goff", ((preselTree->GetEntries()) / m_numBGBins)*m_numBGBins);
698 double* ElectronMomentumDataset = preselTree->GetV1();
699
700 TKDTreeBinning* kdBinsE = new TKDTreeBinning(((preselTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, ElectronMomentumDataset,
702 const double* binsMinEdgesEOriginal = kdBinsE->SortOneDimBinEdges();
703 double* binsMinEdgesE = const_cast<double*>(binsMinEdgesEOriginal);
704 binsMinEdgesE[0] = 0.;
705 binsMinEdgesE[m_numBGBins + 1] = 10000.;
706
707
708 TH2F* hGammaEBetaGamma = new TH2F("hist_d1_11_truncBetaGamma", "hist_d1_11_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
710 binsMinEdgesE,
712 TH2F* hGammaEPart1BetaGamma = (TH2F*)hGammaEBetaGamma->Clone("hist_d1_11_truncPart1BetaGamma");
713 TH2F* hGammaEPart2BetaGamma = (TH2F*)hGammaEBetaGamma->Clone("hist_d1_11_truncPart2BetaGamma");
714
715 preselTree->Draw(Form("FirstElectronSVDdEdx:FirstElectronSVDdEdxTrackMomentum/%f>>hist_d1_11_truncPart1BetaGamma",
717 "FirstElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2 && FirstElectronSVDdEdx<1.8e6", "goff");
718 preselTree->Draw(Form("SecondElectronSVDdEdx:SecondElectronSVDdEdxTrackMomentum/%f>>hist_d1_11_truncPart2BetaGamma",
720 "SecondElectronSVDdEdx>0 && DIRA>0.995 && dr>1.2 && SecondElectronSVDdEdx<1.8e6", "goff");
721 hGammaEBetaGamma->Add(hGammaEPart1BetaGamma);
722 hGammaEBetaGamma->Add(hGammaEPart2BetaGamma);
723
724
725 // produce the 1D profile (for data-MC comparisons)
726 TH1F* ElectronProfileMomentum = (TH1F*)hGammaEMomentum->ProfileX("ElectronProfileMomentum");
727 ElectronProfileMomentum->SetTitle("ElectronProfile");
728 ElectronProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
729 ElectronProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
730 ElectronProfileMomentum->GetYaxis()->SetTitle("dE/dx");
731 ElectronProfileMomentum->SetLineColor(kRed);
732
733// beta*gamma profile for the fit
734 TH1F* ElectronProfileBetaGamma = (TH1F*)hGammaEBetaGamma->ProfileX("ElectronProfileBetaGamma");
735 if (m_CustomProfile) {
736 ElectronProfileBetaGamma = PrepareProfile(hGammaEBetaGamma, "ElectronProfileBetaGamma");
737 }
738 ElectronProfileBetaGamma->SetTitle("ElectronProfile");
739
740 ElectronProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
741 ElectronProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
742 ElectronProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
743 ElectronProfileBetaGamma->SetLineColor(kRed);
744
745
746 hGammaEMomentum = Normalise2DHisto(hGammaEMomentum);
747
748 std::unique_ptr<TList> histList(new TList);
749 histList->Add(ElectronProfileMomentum);
750 histList->Add(ElectronProfileBetaGamma);
751 histList->Add(hGammaEMomentum);
752
753 if (m_isMakePlots) {
754 TFile GammaHistogrammingPlotFile("SVDdEdxCalibrationGammaHistogramming.root", "RECREATE");
755 histList->Write();
756 GammaHistogrammingPlotFile.Close();
757 }
758
759 return histList;
760
761}
const double m_ElectronPDGMass
PDG mass for the electron.

◆ GenerateNewHistograms()

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 )
private

generate high-statistics histograms

Definition at line 763 of file SVDdEdxCalibrationAlgorithm.cc.

766{
767 gROOT->SetBatch(true);
768 gStyle->SetOptStat(0);
769// run the background subtraction and histogramming parts
770 TTree* treeLambda = LambdaMassFit(ttreeLambda);
771 std::unique_ptr<TList> HistListLambda = LambdaHistogramming(treeLambda);
772 TH1F* ProtonProfileBetaGamma = (TH1F*) HistListLambda->FindObject("ProtonProfileBetaGamma");
773 TH2F* Proton2DHistogram = (TH2F*) HistListLambda->FindObject("hist_d1_2212_truncMomentum");
774
775 TTree* treeDstar = DstarMassFit(ttreeDstar);
776 std::unique_ptr<TList> HistListDstar = DstarHistogramming(treeDstar);
777 TH1F* PionProfileBetaGamma = (TH1F*) HistListDstar->FindObject("PionProfileBetaGamma");
778 TH2F* Pion2DHistogram = (TH2F*) HistListDstar->FindObject("hist_d1_211_truncMomentum");
779 TH1F* KaonProfileBetaGamma = (TH1F*) HistListDstar->FindObject("KaonProfileBetaGamma");
780 TH2F* Kaon2DHistogram = (TH2F*) HistListDstar->FindObject("hist_d1_321_truncMomentum");
781
782 std::unique_ptr<TList> HistListGamma = GammaHistogramming(ttreeGamma);
783 TH1F* ElectronProfileBetaGamma = (TH1F*) HistListGamma->FindObject("ElectronProfileBetaGamma");
784 TH2F* Electron2DHistogram = (TH2F*) HistListGamma->FindObject("hist_d1_11_truncMomentum");
785
786 int cred = TColor::GetColor("#e31a1c");
787 PionProfileBetaGamma->SetMarkerSize(4);
788 PionProfileBetaGamma->SetLineWidth(2);
789 PionProfileBetaGamma->SetMarkerColor(cred);
790 PionProfileBetaGamma->SetLineColor(cred);
791
792 int cpink = TColor::GetColor("#807dba");
793 KaonProfileBetaGamma->SetMarkerSize(4);
794 KaonProfileBetaGamma->SetLineWidth(2);
795 KaonProfileBetaGamma->SetMarkerColor(cpink);
796 KaonProfileBetaGamma->SetLineColor(cpink);
797
798 int cblue = TColor::GetColor("#084594");
799 ProtonProfileBetaGamma->SetMarkerSize(4);
800 ProtonProfileBetaGamma->SetLineWidth(2);
801 ProtonProfileBetaGamma->SetMarkerColor(cblue);
802 ProtonProfileBetaGamma->SetLineColor(cblue);
803
804 int cgreen = TColor::GetColor("#238b45");
805 ElectronProfileBetaGamma->SetMarkerSize(4);
806 ElectronProfileBetaGamma->SetLineWidth(2);
807 ElectronProfileBetaGamma->SetMarkerColor(cgreen);
808 ElectronProfileBetaGamma->SetLineColor(cgreen);
809
810// prepare the fitting
811
812 PionProfileBetaGamma->GetYaxis()->SetRangeUser(5.e5, 5.5e6);
813 KaonProfileBetaGamma->GetYaxis()->SetRangeUser(5.e5, 5.5e6);
814 ProtonProfileBetaGamma->GetYaxis()->SetRangeUser(5.e5, 5.5e6);
815
816// enhance the proton histogram (which ends at beta*gamma around 3) by adding pion data above this value – otherwise the fit is very unstable
817 auto PionEdges = PionProfileBetaGamma->GetXaxis()->GetXbins()->GetArray();
818 auto ProtonEdges = ProtonProfileBetaGamma->GetXaxis()->GetXbins()->GetArray();
819
820 std::vector<float> CombinedEdgesVector;
821
822 double borderline = 3.;
823
824 for (int i = 0; i < ProtonProfileBetaGamma->GetNbinsX() + 1; i++)
825 if (ProtonEdges[i] < borderline) CombinedEdgesVector.push_back(ProtonEdges[i]);
826
827
828 for (int i = 0; i < PionProfileBetaGamma->GetNbinsX() + 1; i++)
829 if (PionEdges[i] > borderline) CombinedEdgesVector.push_back(PionEdges[i]);
830
831
832 TH1F* CombinedHistogramPAndPi = new TH1F("CombinedHistogramPAndPi", "histo_for_fit", CombinedEdgesVector.size() - 1,
833 CombinedEdgesVector.data());
834
835 int iterator = 1;
836 for (int i = 1; i < ProtonProfileBetaGamma->GetNbinsX() + 1; i++)
837 if (ProtonEdges[i - 1] < borderline) {
838 CombinedHistogramPAndPi->SetBinContent(i, ProtonProfileBetaGamma->GetBinContent(i));
839 CombinedHistogramPAndPi->SetBinError(i, ProtonProfileBetaGamma->GetBinError(i));
840 iterator++;
841 }
842
843 for (int i = 1; i < PionProfileBetaGamma->GetNbinsX() + 1; i++)
844 if (PionEdges[i - 1] > borderline) {
845
846 CombinedHistogramPAndPi->SetBinContent(iterator, PionProfileBetaGamma->GetBinContent(i));
847 CombinedHistogramPAndPi->SetBinError(iterator, PionProfileBetaGamma->GetBinError(i));
848 iterator++;
849 }
850
851// define the beta*gamma vs momentum function
852 TF1* BetaGammaFunctionPion = new TF1("BetaGammaFunctionPion", "[0] + [1] * x/[2] + [5]/(x^2/[2]^2 + [3])**[4] + [6]* (x/[2])**0.5",
853 0.01, 25.);
854
855 BetaGammaFunctionPion->SetNpx(1000);
856
857 BetaGammaFunctionPion->SetParameters(5.e5, 2.e3, 1, 0.15, 1.2, 6.e5, 3.e5);
858
859 BetaGammaFunctionPion->SetParLimits(0, 3.e5, 7.e5);
860 BetaGammaFunctionPion->SetParLimits(1, -3.e4, 1.e4);
861 BetaGammaFunctionPion->SetParLimits(3, 0.1, 0.2);
862 BetaGammaFunctionPion->SetParLimits(4, 0.9, 1.6);
863 BetaGammaFunctionPion->SetParLimits(5, 3.e5, 7.e5);
864 BetaGammaFunctionPion->SetParLimits(6, 0., 1.e6);
865 BetaGammaFunctionPion->FixParameter(2, 1);
866 if (m_FixUnstableFitParameter) BetaGammaFunctionPion->FixParameter(3, 0.15);
867
868// fit it to the pion data
869 ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
870 auto FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0SI", "", 0.4, 25);
871
872 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
873 BetaGammaFunctionPion->FixParameter(3, 0.15);
874 FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0SI", "", 0.4, 25);
875 }
876 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
877 BetaGammaFunctionPion->FixParameter(3, 0.15);
878 FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0SI", "", 0.45, 25);
879 }
880 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
881 BetaGammaFunctionPion->FixParameter(3, 0.15);
882 FitResultBetaGammaPion = PionProfileBetaGamma->Fit("BetaGammaFunctionPion", "0S", "", 0.5, 25);
883 }
884
885 B2INFO("BetaGamma fit for pions done. Fit status: " << FitResultBetaGammaPion->Status());
886 // B2INFO(FitResultBetaGammaPion->Print(Belle2::LogConfig::c_Info));
887 B2INFO("Fit parameters:");
888 B2INFO("p0: " << BetaGammaFunctionPion->GetParameter(0) << " +- " << BetaGammaFunctionPion->GetParError(0));
889 B2INFO("p1: " << BetaGammaFunctionPion->GetParameter(1) << " +- " << BetaGammaFunctionPion->GetParError(1));
890 B2INFO("p2: " << BetaGammaFunctionPion->GetParameter(2) << " +- " << BetaGammaFunctionPion->GetParError(2));
891 B2INFO("p3: " << BetaGammaFunctionPion->GetParameter(3) << " +- " << BetaGammaFunctionPion->GetParError(3));
892 B2INFO("p4: " << BetaGammaFunctionPion->GetParameter(4) << " +- " << BetaGammaFunctionPion->GetParError(4));
893 B2INFO("p5: " << BetaGammaFunctionPion->GetParameter(5) << " +- " << BetaGammaFunctionPion->GetParError(5));
894 B2INFO("p6: " << BetaGammaFunctionPion->GetParameter(6) << " +- " << BetaGammaFunctionPion->GetParError(6));
895
896 // repeat the same for kaons
897 TF1* BetaGammaFunctionKaon = new TF1("BetaGammaFunctionKaon", "[0] + [1] * x/[2] + [5]/(x^2/[2]^2 + [3])**[4]+ [6]* (x/[2])**0.5",
898 0.01, 25.);
899
900 BetaGammaFunctionKaon->SetNpx(1000);
901 BetaGammaFunctionKaon->SetParameters(5.e5, 2.e3, 1, 0.15, 1.2, 6.e5, 3.e5);
902
903 BetaGammaFunctionKaon->SetParLimits(0, 3.e5, 7.e5);
904 BetaGammaFunctionKaon->SetParLimits(1, -3.e4, 1.e4);
905 BetaGammaFunctionKaon->SetParLimits(3, 0.1, 0.2);
906 BetaGammaFunctionKaon->SetParLimits(4, 0.9, 1.6);
907 BetaGammaFunctionKaon->SetParLimits(5, 3.e5, 7.e5);
908 BetaGammaFunctionKaon->SetParLimits(6, 0., 1.e6);
909
910 BetaGammaFunctionKaon->FixParameter(2, 1);
911 if (m_FixUnstableFitParameter) BetaGammaFunctionKaon->FixParameter(3, 0.15);
912
913 BetaGammaFunctionKaon->SetLineColor(KaonProfileBetaGamma->GetMarkerColor());
914
915 auto FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0SI", "", 0.4, 8.5);
916
917 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
918 BetaGammaFunctionKaon->FixParameter(3, 0.15);
919 FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0SI", "", 0.4, 8.5);
920 }
921 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
922 BetaGammaFunctionKaon->FixParameter(3, 0.15);
923 FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0SI", "", 0.45, 8);
924 }
925 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
926 BetaGammaFunctionKaon->FixParameter(3, 0.15);
927 FitResultBetaGammaKaon = KaonProfileBetaGamma->Fit("BetaGammaFunctionKaon", "0S", "", 0.5, 8);
928 }
929
930 B2INFO("BetaGamma fit for kaons done. Fit status: " << FitResultBetaGammaKaon->Status());
931 B2INFO("Fit parameters:");
932 B2INFO("p0: " << BetaGammaFunctionKaon->GetParameter(0) << " +- " << BetaGammaFunctionKaon->GetParError(0));
933 B2INFO("p1: " << BetaGammaFunctionKaon->GetParameter(1) << " +- " << BetaGammaFunctionKaon->GetParError(1));
934 B2INFO("p2: " << BetaGammaFunctionKaon->GetParameter(2) << " +- " << BetaGammaFunctionKaon->GetParError(2));
935 B2INFO("p3: " << BetaGammaFunctionKaon->GetParameter(3) << " +- " << BetaGammaFunctionKaon->GetParError(3));
936 B2INFO("p4: " << BetaGammaFunctionKaon->GetParameter(4) << " +- " << BetaGammaFunctionKaon->GetParError(4));
937 B2INFO("p5: " << BetaGammaFunctionKaon->GetParameter(5) << " +- " << BetaGammaFunctionKaon->GetParError(5));
938 B2INFO("p6: " << BetaGammaFunctionKaon->GetParameter(6) << " +- " << BetaGammaFunctionKaon->GetParError(6));
939
940 // repeat the same for protons
941 TF1* BetaGammaFunctionProton = new TF1("BetaGammaFunctionProton",
942 "[0] + [1] * x/[2] + [5]/(x^2/[2]^2 + [3])**[4]+ [6]* (x/[2])**0.5", 0.01, 25.);
943
944 BetaGammaFunctionProton->SetNpx(1000);
945
946 BetaGammaFunctionProton->SetParameters(5.e5, 2.e3, 1, 0.15, 1.2, 6.e5, 3.e5);
947
948 BetaGammaFunctionProton->SetParLimits(0, 3.e5, 7.e5);
949 BetaGammaFunctionProton->SetParLimits(1, -3.e4, 1.e4);
950 BetaGammaFunctionProton->SetParLimits(3, 0.1, 0.2);
951 BetaGammaFunctionProton->SetParLimits(4, 0.9, 1.6);
952 BetaGammaFunctionProton->SetParLimits(5, 3.e5, 7.e5);
953 BetaGammaFunctionProton->SetParLimits(6, 0., 1.e6);
954
955 BetaGammaFunctionProton->FixParameter(2, 1);
956 if (m_FixUnstableFitParameter) BetaGammaFunctionProton->FixParameter(3, 0.15);
957
958 BetaGammaFunctionProton->SetLineColor(ProtonProfileBetaGamma->GetMarkerColor());
959
960 auto FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0SI", "", 0.45, 15);
961
962 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
963 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
964 BetaGammaFunctionProton->FixParameter(3, 0.15);
965 FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0SI", "", 0.45, 15);
966 }
967 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
968 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
969 BetaGammaFunctionProton->FixParameter(3, 0.15);
970 FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0SI", "", 0.45, 10);
971 }
972 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
973 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
974 BetaGammaFunctionProton->FixParameter(3, 0.15);
975 FitResultBetaGammaProton = CombinedHistogramPAndPi->Fit("BetaGammaFunctionProton", "0S", "", 0.5, 10);
976 }
977
978 B2INFO("BetaGamma fit for protons done. Fit status: " << FitResultBetaGammaProton->Status());
979 B2INFO("Fit parameters:");
980 B2INFO("p0: " << BetaGammaFunctionProton->GetParameter(0) << " +- " << BetaGammaFunctionProton->GetParError(0));
981 B2INFO("p1: " << BetaGammaFunctionProton->GetParameter(1) << " +- " << BetaGammaFunctionProton->GetParError(1));
982 B2INFO("p2: " << BetaGammaFunctionProton->GetParameter(2) << " +- " << BetaGammaFunctionProton->GetParError(2));
983 B2INFO("p3: " << BetaGammaFunctionProton->GetParameter(3) << " +- " << BetaGammaFunctionProton->GetParError(3));
984 B2INFO("p4: " << BetaGammaFunctionProton->GetParameter(4) << " +- " << BetaGammaFunctionProton->GetParError(4));
985 B2INFO("p5: " << BetaGammaFunctionProton->GetParameter(5) << " +- " << BetaGammaFunctionProton->GetParError(5));
986 B2INFO("p6: " << BetaGammaFunctionProton->GetParameter(6) << " +- " << BetaGammaFunctionProton->GetParError(6));
987
988 if (m_isMakePlots) {
989// plot a summary of all beta*gamma fits for hadrons
990 std::unique_ptr<TCanvas> CombinedCanvasHadrons(new TCanvas("CombinedCanvasHadrons", "Hadron beta*gamma fits", 10, 10, 1000, 700));
991 gStyle->SetOptFit(1111);
992
993 PionProfileBetaGamma->Draw();
994 PionProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionPion);
995 KaonProfileBetaGamma->Draw("SAME");
996 KaonProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionKaon);
997 ProtonProfileBetaGamma->Draw("SAME");
998 ProtonProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionProton);
999// BetaGammaFunctionPion->Draw("SAME");
1000// BetaGammaFunctionKaon->Draw("SAME");
1001// BetaGammaFunctionProton->Draw("SAME");
1002 auto legend = new TLegend(0.4, 0.7, 0.65, 0.9);
1003 legend->AddEntry(PionProfileBetaGamma, "Pions", "lep");
1004 legend->AddEntry(KaonProfileBetaGamma, "Kaons", "lep");
1005 legend->AddEntry(ProtonProfileBetaGamma, "Protons", "lep");
1006 legend->Draw();
1007
1008 gPad->SetLogx();
1009 gPad->SetLogy();
1010
1011 CombinedCanvasHadrons->Print("HadronBetaGammaFits.pdf");
1012 TFile HadronFitPlotFile("SVDdEdxCalibrationHadronFitPlotFile.root", "RECREATE");
1013 PionProfileBetaGamma->Write();
1014 KaonProfileBetaGamma->Write();
1015 ProtonProfileBetaGamma->Write();
1016 CombinedHistogramPAndPi->Write();
1017 BetaGammaFunctionPion->Write();
1018 BetaGammaFunctionKaon->Write();
1019 BetaGammaFunctionProton->Write();
1020 CombinedCanvasHadrons->Write();
1021 HadronFitPlotFile.Close();
1022 }
1023
1024
1025 // in case we assume that all hadrons are equal
1027 BetaGammaFunctionKaon = (TF1*) BetaGammaFunctionPion->Clone("BetaGammaFunctionKaon");
1028 BetaGammaFunctionProton = (TF1*) BetaGammaFunctionPion->Clone("BetaGammaFunctionProton");
1029 }
1031 BetaGammaFunctionKaon = (TF1*) BetaGammaFunctionProton->Clone("BetaGammaFunctionKaon");
1032 BetaGammaFunctionPion = (TF1*) BetaGammaFunctionProton->Clone("BetaGammaFunctionPion");
1033 }
1034
1035 // sanity checks: are all fits ok?
1036 if ((FitResultBetaGammaProton->Status() > 1) || (BetaGammaFunctionProton->Eval(1) < 5.e5)
1037 || (BetaGammaFunctionProton->Eval(1) > 5.e6)) {
1038 if (FitResultBetaGammaPion->Status() == 0) {
1039 BetaGammaFunctionProton = (TF1*) BetaGammaFunctionPion->Clone("BetaGammaFunctionProton");
1040 } else if (FitResultBetaGammaKaon->Status() == 0) {
1041 BetaGammaFunctionProton = (TF1*) BetaGammaFunctionKaon->Clone("BetaGammaFunctionProton");
1042 } else {
1043 B2WARNING("Problem with the beta*gamma fit for protons, reverting to the default values");
1044 BetaGammaFunctionProton->SetParameters(450258, -10900.8, 1, 0.126797, 1.155, 641907, 86304.5);
1045 }
1046 }
1047
1048 if ((FitResultBetaGammaKaon->Status() > 1) || (BetaGammaFunctionKaon->Eval(1) < 5.e5) || (BetaGammaFunctionKaon->Eval(1) > 5.e6)) {
1049 if (FitResultBetaGammaProton->Status() == 0) {
1050 BetaGammaFunctionKaon = (TF1*) BetaGammaFunctionProton->Clone("BetaGammaFunctionKaon");
1051 } else if (FitResultBetaGammaPion->Status() == 0) {
1052 BetaGammaFunctionKaon = (TF1*) BetaGammaFunctionPion->Clone("BetaGammaFunctionKaon");
1053 } else {
1054 B2WARNING("Problem with the beta*gamma fit for kaons, reverting to the default values");
1055 BetaGammaFunctionKaon->SetParameters(543386, 3013.81, 1, 0.135517, 1.19742, 619509, 15484.4);
1056 }
1057 }
1058
1059 if ((FitResultBetaGammaPion->Status() > 1) || (BetaGammaFunctionPion->Eval(1) < 5.e5) || (BetaGammaFunctionPion->Eval(1) > 5.e6)) {
1060 if (FitResultBetaGammaKaon->Status() == 0) {
1061 BetaGammaFunctionPion = (TF1*) BetaGammaFunctionKaon->Clone("BetaGammaFunctionPion");
1062 } else if (FitResultBetaGammaProton->Status() == 0) {
1063 BetaGammaFunctionPion = (TF1*) BetaGammaFunctionProton->Clone("BetaGammaFunctionPion");
1064 } else {
1065 B2WARNING("Problem with the beta*gamma fit for pions, reverting to the default values");
1066 BetaGammaFunctionPion->SetParameters(537623, -1937.62, 1, 0.15292, 1.23803, 623678, 30400.9);
1067 }
1068 }
1069
1070
1071// electrons
1072
1073 TF1* BetaGammaFunctionElectron = new TF1("BetaGammaFunctionElectron", "[0] + [1]* x", 1, 10000.);
1074 BetaGammaFunctionElectron->SetParameters(6.e5, -1);
1075 BetaGammaFunctionElectron->SetParLimits(0, 3.e5, 8.e5);
1076 BetaGammaFunctionElectron->SetParLimits(1, -1.e5, 1.e5);
1077 auto FitResultBetaGammaElectron = ElectronProfileBetaGamma->Fit("BetaGammaFunctionElectron", "0SI", "", 100, 8000);
1078
1079
1080 if ((FitResultBetaGammaElectron->Status() > 1) || (BetaGammaFunctionElectron->Eval(1) < 3.e5)
1081 || (BetaGammaFunctionElectron->Eval(1) > 5.e6)) {
1082 FitResultBetaGammaElectron = ElectronProfileBetaGamma->Fit("BetaGammaFunctionElectron", "0S", "", 100, 10000);
1083 }
1084 B2INFO("BetaGamma fit for electrons done. Fit status: " << FitResultBetaGammaElectron->Status());
1085 B2INFO("Fit parameters:");
1086 B2INFO("p0: " << BetaGammaFunctionElectron->GetParameter(0) << " +- " << BetaGammaFunctionElectron->GetParError(0));
1087 B2INFO("p1: " << BetaGammaFunctionElectron->GetParameter(1) << " +- " << BetaGammaFunctionElectron->GetParError(1));
1088
1089
1090 ElectronProfileBetaGamma->SetMarkerSize(4);
1091 ElectronProfileBetaGamma->SetLineWidth(2);
1092 ElectronProfileBetaGamma->GetYaxis()->SetRangeUser(5e5, 1e6);
1093 ElectronProfileBetaGamma->GetListOfFunctions()->Add(BetaGammaFunctionElectron);
1094
1095 if (m_isMakePlots) {
1096 std::unique_ptr<TCanvas> ElectronCanvas(new TCanvas("ElectronCanvas", "Electron histogram", 10, 10, 1000, 700));
1097 ElectronProfileBetaGamma->Draw();
1098
1099 gPad->SetLogx();
1100
1101 ElectronCanvas->Print("ElectronBetaGammaFits.pdf");
1102 TFile ElectronFitPlotFile("SVDdEdxCalibrationElectronFitPlotFile.root", "RECREATE");
1103 ElectronProfileBetaGamma->Write();
1104 BetaGammaFunctionElectron->Write();
1105 ElectronCanvas->Write();
1106 ElectronFitPlotFile.Close();
1107 }
1108
1109 TF1* MomentumFunctionElectron = (TF1*) BetaGammaFunctionElectron->Clone("MomentumFunctionElectron");
1110 MomentumFunctionElectron->SetParameter(2, m_ElectronPDGMass);
1111 MomentumFunctionElectron->SetRange(0.01, 5.5);
1112 MomentumFunctionElectron->SetLineColor(kRed);
1113 MomentumFunctionElectron->SetLineWidth(4);
1114
1115 TF1* MomentumFunctionPion = (TF1*) BetaGammaFunctionPion->Clone("MomentumFunctionPion");
1116 MomentumFunctionPion->SetParameter(2, m_PionPDGMass);
1117 MomentumFunctionPion->SetRange(0.01, 5.5);
1118 MomentumFunctionPion->SetLineColor(kRed);
1119 MomentumFunctionPion->SetLineWidth(4);
1120
1121 TF1* MomentumFunctionProton = (TF1*) BetaGammaFunctionProton->Clone("MomentumFunctionProton");
1122 MomentumFunctionProton->SetParameter(2, m_ProtonPDGMass);
1123 MomentumFunctionProton->SetRange(0.01, 5.5);
1124 MomentumFunctionProton->SetLineColor(kRed);
1125 MomentumFunctionProton->SetLineWidth(4);
1126
1127 TF1* MomentumFunctionKaon = (TF1*) BetaGammaFunctionKaon->Clone("MomentumFunctionKaon");
1128 MomentumFunctionKaon->SetParameter(2, m_KaonPDGMass);
1129 MomentumFunctionKaon->SetRange(0.01, 5.5);
1130 MomentumFunctionKaon->SetLineColor(kRed);
1131 MomentumFunctionKaon->SetLineWidth(4);
1132
1133 gStyle->SetOptFit(1111);
1134 std::unique_ptr<TCanvas> CanvasOverlays(new TCanvas("CanvasOverlays", "overlays", 1300, 1000));
1135 CanvasOverlays->Divide(2, 2);
1136 CanvasOverlays->cd(1); Electron2DHistogram->Draw(); MomentumFunctionElectron->Draw("SAME");
1137 CanvasOverlays->cd(2); Pion2DHistogram->Draw(); MomentumFunctionPion->Draw("SAME");
1138 CanvasOverlays->cd(3); Kaon2DHistogram->Draw(); MomentumFunctionKaon->Draw("SAME");
1139 CanvasOverlays->cd(4); Proton2DHistogram->Draw(); MomentumFunctionProton->Draw("SAME");
1140 CanvasOverlays->Print("SVDdEdxOverlaysFitsHistos.pdf");
1141
1142 TF1* MomentumFunctionDeuteron = (TF1*) BetaGammaFunctionProton->Clone("MomentumFunctionDeuteron");
1143 MomentumFunctionDeuteron->SetParameter(2, m_DeuteronPDGMass);
1144 MomentumFunctionDeuteron->SetRange(0.01, 5.5);
1145 MomentumFunctionDeuteron->SetLineColor(kRed);
1146
1147 TF1* MomentumFunctionMuon = (TF1*) BetaGammaFunctionPion->Clone("MomentumFunctionMuon");
1148 MomentumFunctionMuon->SetParameter(2, m_MuonPDGMass);
1149 MomentumFunctionMuon->SetRange(0.01, 5.5);
1150 MomentumFunctionMuon->SetLineColor(kRed);
1151
1152// overlay all fits in one plot
1153
1154 std::unique_ptr<TCanvas> OverlayAllTracksCanvas(new TCanvas("OverlayAllTracksCanvas", "The Ultimate Plot", 10, 10, 1000, 700));
1155
1156 TH2F* AllTracksHistogram = new TH2F("AllTracksHistogram", "AllTracksHistogram;Momentum [GeV/c];dEdx [arb. units]", 1000, 0.05, 5,
1157 1000, 2.e5, 6.e6);
1158
1159 ttreeGeneric->Draw("TrackSVDdEdx:TrackSVDdEdxTrackMomentum>>AllTracksHistogram", "TracknSVDHits>7", "goff");
1160 AllTracksHistogram->Draw("COLZ");
1161 AllTracksHistogram->GetXaxis()->SetTitle("Momentum [GeV/c]");
1162 AllTracksHistogram->GetYaxis()->SetTitle("dE/dx [arbitrary units]");
1163 MomentumFunctionElectron->Draw("SAME");
1164 MomentumFunctionMuon->Draw("SAME");
1165 MomentumFunctionPion->Draw("SAME");
1166 MomentumFunctionKaon->Draw("SAME");
1167 MomentumFunctionProton->Draw("SAME");
1168 MomentumFunctionDeuteron->Draw("SAME");
1169 OverlayAllTracksCanvas->SetLogx();
1170 OverlayAllTracksCanvas->SetLogz();
1171
1172 OverlayAllTracksCanvas->Print("SVDdEdxAllTracksWithFits.pdf");
1173 TFile OverlayAllTracksPlotFile("SVDdEdxCalibrationOverlayAllTracks.root", "RECREATE");
1174 AllTracksHistogram->Write();
1175 MomentumFunctionElectron->Write();
1176 MomentumFunctionMuon->Write();
1177 MomentumFunctionPion->Write();
1178 MomentumFunctionKaon->Write();
1179 MomentumFunctionProton->Write();
1180 MomentumFunctionDeuteron->Write();
1181 OverlayAllTracksCanvas->Write();
1182 OverlayAllTracksPlotFile.Close();
1183
1184
1185// resolution studies //
1186
1187// For resolution measurement, we need to take a ProjectionY of the data histograms in the momentum range where the dEdx is flat vs momentum. We use our educated guess of the flat range (e.g. 0.6-1 GeV for pions) and FindBin to figure out which bin numbers those momentum values correspond to.
1188 double PionRangeMin = 0.6;
1189 double PionRangeMax = 1.;
1190 double KaonRangeMin = 1.9;
1191 double KaonRangeMax = 3;
1192 double ElectronRangeMin = 1.;
1193 double ElectronRangeMax = 1.4;
1194
1195 auto PionResolutionHistogram = Pion2DHistogram->ProjectionY("PionResolutionHistogram",
1196 Pion2DHistogram->GetXaxis()->FindBin(PionRangeMin),
1197 Pion2DHistogram->GetXaxis()->FindBin(PionRangeMax));
1198 auto ElectronResolutionHistogram = Electron2DHistogram->ProjectionY("ElectronResolutionHistogram",
1199 Electron2DHistogram->GetXaxis()->FindBin(ElectronRangeMin), Electron2DHistogram->GetXaxis()->FindBin(ElectronRangeMax));
1200 auto KaonResolutionHistogram = Kaon2DHistogram->ProjectionY("KaonResolutionHistogram",
1201 Kaon2DHistogram->GetXaxis()->FindBin(KaonRangeMin),
1202 Kaon2DHistogram->GetXaxis()->FindBin(KaonRangeMax));
1203// for protons, there is not enough data in the flat range.
1204
1205
1206 TF1* PionResolutionFunction = new TF1("PionResolutionFunction",
1207 "[0]*TMath::Landau(x, [1], [1]*[2])*TMath::Gaus(x, [1], [1]*[2]*[4]) + [3]*TMath::Gaus(x, [1], [1]*[2]*[5])", 100e3, 1500e3);
1208// parameter [1] is the mean of the Landau
1209// parameter [2] is the relative resolution (w.r.t. the mean) of the Landau
1210// parameters [4]-[5] are the relative resolution of Gauss contributions w.r.t. that of Landau
1211// parameters [0] and [3] are fractions of the two components
1212 PionResolutionFunction->SetParameters(1, 6.e5, 0.1, 0.5, 2, 1);
1213 PionResolutionFunction->SetParLimits(0, 0, 500);
1214 PionResolutionFunction->SetParLimits(1, 3.e5, 8.e5);
1215 PionResolutionFunction->SetParLimits(2, 0, 1);
1216 PionResolutionFunction->SetParLimits(3, 0, 500);
1217 PionResolutionFunction->SetParLimits(4, 0, 7);
1218 PionResolutionFunction->SetParLimits(5, 1, 7);
1219 PionResolutionFunction->SetNpx(1000);
1220 auto FitResultResolutionPion = PionResolutionHistogram->Fit(PionResolutionFunction, "RSI");
1221
1222 B2INFO("relative resolution for pions: " << PionResolutionFunction->GetParameter(2));
1223 B2INFO("resolution for pions: fit status" << FitResultResolutionPion->Status());
1224
1225 TF1* KaonResolutionFunction = new TF1("KaonResolutionFunction",
1226 "[0]*TMath::Landau(x, [1], [1]*[2])*TMath::Gaus(x, [1], [1]*[2]*[4]) + [3]*TMath::Gaus(x, [1], [1]*[2]*[5])", 100e3, 1500e3);
1227
1228
1229 KaonResolutionFunction->SetParameters(1, 6.e5, 0.1, 0.5, 2, 1);
1230 KaonResolutionFunction->SetParLimits(0, 0, 500);
1231 KaonResolutionFunction->SetParLimits(1, 3.e5, 8.e5);
1232 KaonResolutionFunction->SetParLimits(2, 0, 1);
1233 KaonResolutionFunction->SetParLimits(3, 0, 500);
1234 KaonResolutionFunction->SetParLimits(4, 0, 7);
1235 KaonResolutionFunction->SetParLimits(5, 1, 7);
1236 KaonResolutionFunction->SetNpx(1000);
1237 auto FitResultResolutionKaon = KaonResolutionHistogram->Fit(KaonResolutionFunction, "RSI");
1238
1239 B2INFO("relative resolution for kaons: " << KaonResolutionFunction->GetParameter(2));
1240 B2INFO("resolution for kaons: fit status" << FitResultResolutionKaon->Status());
1241
1242 if ((FitResultResolutionKaon->Status() > 1)
1243 && (FitResultResolutionPion->Status() <= 1)) KaonResolutionFunction = (TF1*)
1244 PionResolutionFunction->Clone("KaonResolutionFunction");
1245
1246
1247
1248
1249 TF1* ElectronResolutionFunction = new TF1("ElectronResolutionFunction",
1250 "[0]*TMath::Landau(x, [1], [1]*[2])*TMath::Gaus(x, [1], [1]*[2]*[4]) + [3]*TMath::Gaus(x, [1], [1]*[2]*[5])", 50e3, 1500e3);
1251
1252
1253 ElectronResolutionFunction->SetParameters(1, 6.e5, 0.1, 0.5, 2, 1);
1254 ElectronResolutionFunction->SetParLimits(0, 0, 500);
1255 ElectronResolutionFunction->SetParLimits(1, 3.e5, 8.e5);
1256 ElectronResolutionFunction->SetParLimits(2, 0, 1);
1257 ElectronResolutionFunction->SetParLimits(3, 0, 500);
1258 ElectronResolutionFunction->SetParLimits(4, 0, 7);
1259 ElectronResolutionFunction->SetParLimits(5, 1, 7);
1260 ElectronResolutionFunction->SetNpx(1000);
1261 auto FitResultResolutionElectron = ElectronResolutionHistogram->Fit(ElectronResolutionFunction, "RSI");
1262
1263 B2INFO("relative resolution for electrons: " << ElectronResolutionFunction->GetParameter(2));
1264 B2INFO("resolution for electrons: fit status" << FitResultResolutionElectron->Status());
1265
1266 // plot all the resolution fits
1267 if (m_isMakePlots) {
1268 TCanvas* CanvasResolutions = new TCanvas("CanvasResolutions", "Resolutions", 1200, 650);
1269 CanvasResolutions->Divide(3, 1);
1270 CanvasResolutions->cd(1); PionResolutionHistogram->Draw();
1271 CanvasResolutions->cd(2); KaonResolutionHistogram->Draw();
1272 CanvasResolutions->cd(3); ElectronResolutionHistogram->Draw();
1273
1274 CanvasResolutions->Print("SVDdEdxResolutions.pdf");
1275 TFile OverlayResolutionsPlotFile("SVDdEdxCalibrationResolutions.root", "RECREATE");
1276 PionResolutionHistogram->Write();
1277 KaonResolutionHistogram->Write();
1278 ElectronResolutionHistogram->Write();
1279 CanvasResolutions->Write();
1280 OverlayResolutionsPlotFile.Close();
1281 }
1282
1283// evaluate the bias correction:
1284// difference between the MomentumFunction prediction and the mean of the resolution function in the flat part
1285// it should be of the order -1e4, i.e. about -1% of the absolute dEdx value
1286 double BiasCorrectionPion = PionResolutionFunction->GetParameter(1) - MomentumFunctionPion->Eval((
1287 PionRangeMax + PionRangeMin) / 2.);
1288 B2INFO("BiasCorrectionPion = " << BiasCorrectionPion);
1289
1290// generate a new pion payload using the MomentumFunctionPion, PionResolutionFunction and the bias correction
1291 TH2F* Pion2DHistogramNew = PrepareNewHistogram(Pion2DHistogram, Form("%sNew", Pion2DHistogram->GetName()), MomentumFunctionPion,
1292 PionResolutionFunction, BiasCorrectionPion);
1293
1294// sanity check: residual between the generated distribution and the data one
1295 TH2F* Pion2DHistogramResidual = (TH2F*) Pion2DHistogram->Clone("Pion2DHistogramResidual");
1296 Pion2DHistogramResidual->Add(Pion2DHistogramNew, Pion2DHistogram, 1, -1);
1297 Pion2DHistogramResidual->SetMinimum(-0.15);
1298 Pion2DHistogramResidual->SetMaximum(0.15);
1299
1300 // repeat, for kaons
1301 double BiasCorrectionKaon = KaonResolutionFunction->GetParameter(1) - MomentumFunctionKaon->Eval((
1302 KaonRangeMax + KaonRangeMin) / 2.);
1303 B2INFO("BiasCorrectionKaon = " << BiasCorrectionKaon);
1304
1305 // for protons, we compare the flat part of the MomentumFunction (~3 GeV) with the mean of the kaon resolution function
1306 // as there's not enough stats in the flat part to extract proton resolution from data
1307 double BiasCorrectionProton = KaonResolutionFunction->GetParameter(1) - MomentumFunctionProton->Eval(3.);
1308 B2INFO("BiasCorrectionProton = " << BiasCorrectionProton);
1309
1310 if ((BiasCorrectionProton / BiasCorrectionKaon) > 1.5) BiasCorrectionProton =
1311 BiasCorrectionKaon; // probably something went wrong due to low statistics
1312
1313 // back to kaons: generate a new payload
1314 TH2F* Kaon2DHistogramNew = PrepareNewHistogram(Kaon2DHistogram, Form("%sNew", Kaon2DHistogram->GetName()), MomentumFunctionKaon,
1315 KaonResolutionFunction, BiasCorrectionKaon);
1316// residual generated - data for kaons
1317 TH2F* Kaon2DHistogramResidual = (TH2F*) Kaon2DHistogram->Clone("Kaon2DHistogramResidual");
1318 Kaon2DHistogramResidual->Add(Kaon2DHistogramNew, Kaon2DHistogram, 1, -1);
1319 Kaon2DHistogramResidual->SetMinimum(-0.15);
1320 Kaon2DHistogramResidual->SetMaximum(0.15);
1321
1322 // same for protons (we use the kaon resolution function as explained above)
1323 TH2F* Proton2DHistogramNew = PrepareNewHistogram(Proton2DHistogram, Form("%sNew", Proton2DHistogram->GetName()),
1324 MomentumFunctionProton,
1325 KaonResolutionFunction, BiasCorrectionProton);
1326
1327// residual for protons
1328 TH2F* Proton2DHistogramResidual = (TH2F*) Proton2DHistogram->Clone("Proton2DHistogramResidual");
1329 Proton2DHistogramResidual->Add(Proton2DHistogramNew, Proton2DHistogram, 1, -1);
1330 Proton2DHistogramResidual->SetMinimum(-0.15);
1331 Proton2DHistogramResidual->SetMaximum(0.15);
1332
1333// deuterons: same as protons, but use the MomentumFunctionDeuteron
1334 TH2F* Deuteron2DHistogramNew = PrepareNewHistogram(Proton2DHistogram, "Deuteron2DHistogramNew", MomentumFunctionDeuteron,
1335 KaonResolutionFunction,
1336 BiasCorrectionKaon);
1337 Deuteron2DHistogramNew->SetTitle("hist_d1_1000010020_trunc");
1338
1339// muons: same as pions, but use the MomentumFunctionMuon
1340 TH2F* Muon2DHistogramNew = PrepareNewHistogram(Pion2DHistogram, "Muon2DHistogramNew", MomentumFunctionMuon, PionResolutionFunction,
1341 BiasCorrectionPion);
1342 Muon2DHistogramNew->SetTitle("hist_d1_13_trunc");
1343
1344// same for electrons
1345 double BiasCorrectionElectron = ElectronResolutionFunction->GetParameter(1) - MomentumFunctionElectron->Eval((
1346 ElectronRangeMax + ElectronRangeMin) / 2.);
1347 B2INFO("BiasCorrectionElectron = " << BiasCorrectionElectron);
1348 TH2F* Electron2DHistogramNew = PrepareNewHistogram(Electron2DHistogram, Form("%sNew", Electron2DHistogram->GetName()),
1349 MomentumFunctionElectron,
1350 ElectronResolutionFunction, BiasCorrectionElectron);
1351
1352 TH2F* Electron2DHistogramResidual = (TH2F*) Electron2DHistogram->Clone("Electron2DHistogramResidual");
1353 Electron2DHistogramResidual->Add(Electron2DHistogramNew, Electron2DHistogram, 1, -1);
1354 Electron2DHistogramResidual->SetMinimum(-0.15);
1355 Electron2DHistogramResidual->SetMaximum(0.15);
1356
1357 Electron2DHistogramNew->SetName("Electron2DHistogramNew");
1358 Muon2DHistogramNew->SetName("Muon2DHistogramNew");
1359 Pion2DHistogramNew->SetName("Pion2DHistogramNew");
1360 Kaon2DHistogramNew->SetName("Kaon2DHistogramNew");
1361 Proton2DHistogramNew->SetName("Proton2DHistogramNew");
1362 Deuteron2DHistogramNew->SetName("Deuteron2DHistogramNew");
1363
1364// plot the summary of all the distributions
1365 if (m_isMakePlots) {
1366 TCanvas* CanvasSummaryGenerated = new TCanvas("CanvasSummaryGenerated", "Generated payloads", 1700, 850);
1367 CanvasSummaryGenerated->Divide(3, 2);
1368 CanvasSummaryGenerated->cd(1); Electron2DHistogramNew->Draw("COLZ");
1369 CanvasSummaryGenerated->cd(2); Muon2DHistogramNew->Draw("COLZ");
1370 CanvasSummaryGenerated->cd(3); Pion2DHistogramNew->Draw("COLZ");
1371 CanvasSummaryGenerated->cd(4); Kaon2DHistogramNew->Draw("COLZ");
1372 CanvasSummaryGenerated->cd(5); Proton2DHistogramNew->Draw("COLZ");
1373 CanvasSummaryGenerated->cd(6); Deuteron2DHistogramNew->Draw("COLZ");
1374
1375 CanvasSummaryGenerated->Print("SVDdEdxGeneratedPayloads.pdf");
1376 TFile SummaryGeneratedPlotFile("SVDdEdxCalibrationSummaryGenerated.root", "RECREATE");
1377 Electron2DHistogramNew->Write();
1378 Muon2DHistogramNew->Write();
1379 Pion2DHistogramNew->Write();
1380 Kaon2DHistogramNew->Write();
1381 Proton2DHistogramNew->Write();
1382 Deuteron2DHistogramNew->Write();
1383 SummaryGeneratedPlotFile.Close();
1384
1385
1386 TCanvas* CanvasSummaryData = new TCanvas("CanvasSummaryData", "Data distributions", 1700, 850);
1387 CanvasSummaryData->Divide(3, 2);
1388 CanvasSummaryData->cd(1); Electron2DHistogram->Draw("COLZ");
1389 CanvasSummaryData->cd(3); Pion2DHistogram->Draw("COLZ");
1390 CanvasSummaryData->cd(4); Kaon2DHistogram->Draw("COLZ");
1391 CanvasSummaryData->cd(5); Proton2DHistogram->Draw("COLZ");
1392
1393 CanvasSummaryData->Print("SVDdEdxDataDistributions.pdf");
1394 TFile SummaryDataPlotFile("SVDdEdxCalibrationSummaryData.root", "RECREATE");
1395 Electron2DHistogram->Write();
1396 Pion2DHistogram->Write();
1397 Kaon2DHistogram->Write();
1398 Proton2DHistogram->Write();
1399 SummaryDataPlotFile.Close();
1400
1401
1402 TCanvas* CanvasSummaryResiduals = new TCanvas("CanvasSummaryResiduals", "Residuals", 1700, 850);
1403 CanvasSummaryResiduals->Divide(3, 2);
1404 CanvasSummaryResiduals->cd(1); Electron2DHistogramResidual->Draw("COLZ");
1405 CanvasSummaryResiduals->cd(3); Pion2DHistogramResidual->Draw("COLZ");
1406 CanvasSummaryResiduals->cd(4); Kaon2DHistogramResidual->Draw("COLZ");
1407 CanvasSummaryResiduals->cd(5); Proton2DHistogramResidual->Draw("COLZ");
1408
1409
1410 CanvasSummaryResiduals->Print("SVDdEdxResiduals.pdf");
1411 TFile SummaryResidualsPlotFile("SVDdEdxCalibrationSummaryResiduals.root", "RECREATE");
1412 Electron2DHistogramResidual->Write();
1413 Pion2DHistogramResidual->Write();
1414 Kaon2DHistogramResidual->Write();
1415 Proton2DHistogramResidual->Write();
1416 SummaryResidualsPlotFile.Close();
1417 }
1418
1419
1420 // return all the generated payloads
1421 std::unique_ptr<TList> histList(new TList);
1422 histList->Add(Electron2DHistogramNew);
1423 histList->Add(Muon2DHistogramNew);
1424 histList->Add(Pion2DHistogramNew);
1425 histList->Add(Kaon2DHistogramNew);
1426 histList->Add(Proton2DHistogramNew);
1427 histList->Add(Deuteron2DHistogramNew);
1428
1429 return histList;
1430
1431}
const double m_MuonPDGMass
PDG mass for the muon.
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.
bool m_UseProtonBGFunctionForEverything
Assume that the dEdx:betagamma trend is the same for all hadrons; use the proton trend as representat...
std::unique_ptr< TList > DstarHistogramming(TTree *inputTree)
produce histograms for K/pi
TTree * LambdaMassFit(std::shared_ptr< TTree > preselTree)
Mass fit for Lambda->ppi.
std::unique_ptr< TList > LambdaHistogramming(TTree *inputTree)
produce histograms for protons
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
TTree * DstarMassFit(std::shared_ptr< TTree > preselTree)
Mass fit for D*->Dpi.
bool m_UsePionBGFunctionForEverything
Assume that the dEdx:betagamma trend is the same for all hadrons; use the pion trend as representativ...
const double m_ProtonPDGMass
PDG mass for the proton.

◆ getAllGranularityExpRun()

Calibration::ExpRun getAllGranularityExpRun ( ) const
inlineprotectedinherited

Returns the Exp,Run pair that means 'Everything'. Currently unused.

Definition at line 327 of file CalibrationAlgorithm.h.

327{return m_allExpRun;}

◆ getCollectorName()

std::string getCollectorName ( ) const
inlineinherited

Alias for prefix.

For convenience and less writing, we say developers to set this to default collector module name in constructor of base class. One can however use the dublets of collector+algorithm multiple times with different settings. To bind these together correctly, the prefix has to be set the same for algo and collector. So we call the setter setPrefix rather than setModuleName or whatever. This getter will work out of the box for default cases -> return the name of module you have to add to your path to collect data for this algorithm.

Definition at line 164 of file CalibrationAlgorithm.h.

164{return getPrefix();}

◆ getDescription()

const std::string & getDescription ( ) const
inlineinherited

Get the description of the algorithm (set by developers in constructor)

Definition at line 216 of file CalibrationAlgorithm.h.

216{return m_description;}

◆ getExpRunString()

string getExpRunString ( Calibration::ExpRun & expRun) const
privateinherited

Gets the "exp.run" string repr. of (exp,run)

Definition at line 254 of file CalibrationAlgorithm.cc.

255{
256 string expRunString;
257 expRunString += to_string(expRun.first);
258 expRunString += ".";
259 expRunString += to_string(expRun.second);
260 return expRunString;
261}

◆ getFullObjectPath()

string getFullObjectPath ( const std::string & name,
Calibration::ExpRun expRun ) const
privateinherited

constructs the full TDirectory + Key name of an object in a TFile based on its name and exprun

Definition at line 263 of file CalibrationAlgorithm.cc.

264{
265 string dirName = getPrefix() + "/" + name;
266 string objName = name + "_" + getExpRunString(expRun);
267 return dirName + "/" + objName;
268}
std::string getExpRunString(Calibration::ExpRun &expRun) const
Gets the "exp.run" string repr. of (exp,run)

◆ getGranularity()

std::string getGranularity ( ) const
inlineinherited

Get the granularity of collected data.

Definition at line 188 of file CalibrationAlgorithm.h.

188{return m_granularityOfData;};

◆ getGranularityFromData()

string getGranularityFromData ( ) const
protectedinherited

Get the granularity of collected data.

Definition at line 383 of file CalibrationAlgorithm.cc.

384{
385 // Save TDirectory to change back at the end
386 TDirectory* dir = gDirectory;
387 RunRange* runRange;
388 string runRangeObjName(getPrefix() + "/" + RUN_RANGE_OBJ_NAME);
389 // We only check the first file
390 string fileName = m_inputFileNames[0];
391 unique_ptr<TFile> f;
392 f.reset(TFile::Open(fileName.c_str(), "READ"));
393 runRange = dynamic_cast<RunRange*>(f->Get(runRangeObjName.c_str()));
394 if (!runRange) {
395 B2FATAL("The input file " << fileName << " does not contain a RunRange object at "
396 << runRangeObjName << ". Please set your input files to exclude it.");
397 return "";
398 }
399 string granularity = runRange->getGranularity();
400 dir->cd();
401 return granularity;
402}
std::string getGranularity() const
Gets the m_granularity.
Definition RunRange.h:110

◆ getInputFileNames()

PyObject * getInputFileNames ( )
inherited

Get the input file names used for this algorithm and pass them out as a Python list of unicode strings.

Definition at line 245 of file CalibrationAlgorithm.cc.

246{
247 PyObject* objInputFileNames = PyList_New(m_inputFileNames.size());
248 for (size_t i = 0; i < m_inputFileNames.size(); ++i) {
249 PyList_SetItem(objInputFileNames, i, Py_BuildValue("s", m_inputFileNames[i].c_str()));
250 }
251 return objInputFileNames;
252}

◆ getInputJsonObject()

const nlohmann::json & getInputJsonObject ( ) const
inlineprotectedinherited

Get the entire top level JSON object. We explicitly say this must be of object type so that we might pick.

Definition at line 357 of file CalibrationAlgorithm.h.

357{return m_jsonExecutionInput;}

◆ getInputJsonValue()

template<class T>
const T getInputJsonValue ( const std::string & key) const
inlineprotectedinherited

Get an input JSON value using a key. The normal exceptions are raised when the key doesn't exist.

Definition at line 350 of file CalibrationAlgorithm.h.

351 {
352 return m_jsonExecutionInput.at(key);
353 }

◆ getIovFromAllData()

IntervalOfValidity getIovFromAllData ( ) const
inherited

Get the complete IoV from inspection of collected data.

Definition at line 325 of file CalibrationAlgorithm.cc.

326{
328}
RunRange getRunRangeFromAllData() const
Get the complete RunRange from inspection of collected data.
IntervalOfValidity getIntervalOfValidity()
Make IntervalOfValidity from the set, spanning all runs. Works because sets are sorted by default.
Definition RunRange.h:70

◆ getIteration()

int getIteration ( ) const
inlineprotectedinherited

Get current iteration.

Definition at line 269 of file CalibrationAlgorithm.h.

269{ return m_data.getIteration(); }

◆ getObjectPtr()

template<class T>
std::shared_ptr< T > getObjectPtr ( std::string name)
inlineprotectedinherited

Get calibration data object (for all runs the calibration is requested for) This function will only work during or after execute() has been called once.

Definition at line 285 of file CalibrationAlgorithm.h.

286 {
287 if (m_runsToInputFiles.size() == 0)
288 fillRunToInputFilesMap();
289 return getObjectPtr<T>(name, m_data.getRequestedRuns());
290 }

◆ getOutputJsonValue()

template<class T>
const T getOutputJsonValue ( const std::string & key) const
inlineprotectedinherited

Get a value using a key from the JSON output object, not sure why you would want to do this.

Definition at line 342 of file CalibrationAlgorithm.h.

343 {
344 return m_jsonExecutionOutput.at(key);
345 }

◆ getPayloads()

std::list< Database::DBImportQuery > & getPayloads ( )
inlineinherited

Get constants (in TObjects) for database update from last execution.

Definition at line 204 of file CalibrationAlgorithm.h.

204{return m_data.getPayloads();}

◆ getPayloadValues()

std::list< Database::DBImportQuery > getPayloadValues ( )
inlineinherited

Get constants (in TObjects) for database update from last execution but passed by VALUE.

Definition at line 207 of file CalibrationAlgorithm.h.

207{return m_data.getPayloadValues();}

◆ getPrefix()

std::string getPrefix ( ) const
inlineinherited

Get the prefix used for getting calibration data.

Definition at line 146 of file CalibrationAlgorithm.h.

146{return m_prefix;}

◆ getRunList()

const std::vector< Calibration::ExpRun > & getRunList ( ) const
inlineprotectedinherited

Get the list of runs for which calibration is called.

Definition at line 266 of file CalibrationAlgorithm.h.

266{return m_data.getRequestedRuns();}

◆ getRunListFromAllData()

vector< ExpRun > getRunListFromAllData ( ) const
inherited

Get the complete list of runs from inspection of collected data.

Definition at line 318 of file CalibrationAlgorithm.cc.

319{
320 RunRange runRange = getRunRangeFromAllData();
321 set<ExpRun> expRunSet = runRange.getExpRunSet();
322 return vector<ExpRun>(expRunSet.begin(), expRunSet.end());
323}

◆ getRunRangeFromAllData()

RunRange getRunRangeFromAllData ( ) const
inherited

Get the complete RunRange from inspection of collected data.

Definition at line 361 of file CalibrationAlgorithm.cc.

362{
363 // Save TDirectory to change back at the end
364 TDirectory* dir = gDirectory;
365 RunRange runRange;
366 // Construct the TDirectory name where we expect our objects to be
367 string runRangeObjName(getPrefix() + "/" + RUN_RANGE_OBJ_NAME);
368 for (const auto& fileName : m_inputFileNames) {
369 //Open TFile to get the objects
370 unique_ptr<TFile> f;
371 f.reset(TFile::Open(fileName.c_str(), "READ"));
372 RunRange* runRangeOther = dynamic_cast<RunRange*>(f->Get(runRangeObjName.c_str()));
373 if (runRangeOther) {
374 runRange.merge(runRangeOther);
375 } else {
376 B2WARNING("Missing a RunRange object for file: " << fileName);
377 }
378 }
379 dir->cd();
380 return runRange;
381}
virtual void merge(const RunRange *other)
Implementation of merging - other is added to the set (union)
Definition RunRange.h:52

◆ getVecInputFileNames()

std::vector< std::string > getVecInputFileNames ( ) const
inlineprotectedinherited

Get the input file names used for this algorithm as a STL vector.

Definition at line 275 of file CalibrationAlgorithm.h.

275{return m_inputFileNames;}

◆ inputJsonKeyExists()

bool inputJsonKeyExists ( const std::string & key) const
inlineprotectedinherited

Test for a key in the input JSON object.

Definition at line 360 of file CalibrationAlgorithm.h.

360{return m_jsonExecutionInput.count(key);}

◆ isBoundaryRequired()

virtual bool isBoundaryRequired ( const Calibration::ExpRun & )
inlineprotectedvirtualinherited

Given the current collector data, make a decision about whether or not this run should be the start of a payload boundary.

Reimplemented in PXDAnalyticGainCalibrationAlgorithm, PXDValidationAlgorithm, SVD3SampleCoGTimeCalibrationAlgorithm, SVD3SampleELSTimeCalibrationAlgorithm, SVDCoGTimeCalibrationAlgorithm, TestBoundarySettingAlgorithm, and TestCalibrationAlgorithm.

Definition at line 243 of file CalibrationAlgorithm.h.

244 {
245 B2ERROR("You didn't implement a isBoundaryRequired() member function in your CalibrationAlgorithm but you are calling it!");
246 return false;
247 }

◆ LambdaHistogramming()

std::unique_ptr< TList > LambdaHistogramming ( TTree * inputTree)
private

produce histograms for protons

Definition at line 311 of file SVDdEdxCalibrationAlgorithm.cc.

312{
313 gROOT->SetBatch(true);
314 inputTree->SetEstimate(-1);
315 std::vector<double> pbins = CreatePBinningScheme();
316
317 TH2F* hLambdaPMomentum = new TH2F("hist_d1_2212_truncMomentum", "hist_d1_2212_trunc;Momentum [GeV/c];dEdx [arb. units]",
318 m_numPBins, pbins.data(), m_numDEdxBins, 0,
320
321 inputTree->Draw("ProtonSVDdEdx:ProtonSVDdEdxTrackMomentum>>hist_d1_2212_truncMomentum",
322 "nSignalLambda_sw * (ProtonSVDdEdx>0) * (ProtonSVDdEdxTrackMomentum>0.13)", "goff");
323
324// create isopopulated beta*gamma binning
325 inputTree->Draw(Form("ProtonSVDdEdxTrackMomentum/%f", m_ProtonPDGMass), "", "goff",
326 ((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins);
327 double* ProtonMomentumDataset = inputTree->GetV1();
328
329 TKDTreeBinning* kdBinsP = new TKDTreeBinning(((inputTree->GetEntries()) / m_numBGBins)*m_numBGBins, 1, ProtonMomentumDataset,
331 const double* binsMinEdgesP_pointer = kdBinsP->SortOneDimBinEdges();
332 double* binsMinEdgesP = const_cast<double*>(binsMinEdgesP_pointer);
333
334
335 binsMinEdgesP[0] = 0.1;
336 binsMinEdgesP[m_numBGBins + 1] = 50.;
337
338
339 TH2F* hLambdaPBetaGamma = new TH2F("hist_d1_2212_truncBetaGamma", "hist_d1_2212_truncBetaGamma;#beta*#gamma;dEdx [arb. units]",
340 m_numBGBins, binsMinEdgesP, m_numDEdxBins,
341 0,
343
344 inputTree->Draw(Form("ProtonSVDdEdx:ProtonSVDdEdxTrackMomentum/%f>>hist_d1_2212_truncBetaGamma", m_ProtonPDGMass),
345 "nSignalLambda_sw * (ProtonSVDdEdx>0) * (ProtonSVDdEdxTrackMomentum>0.13) * (ProtonSVDdEdx>1.2e6 - 1.e6*ProtonSVDdEdxTrackMomentum)",
346 "goff");
347
348 // produce the 1D profile
349 // momentum: for data-MC comparisons
350
351 TH1F* ProtonProfileMomentum = (TH1F*)hLambdaPMomentum->ProfileX("ProtonProfileMomentum");
352 ProtonProfileMomentum->SetTitle("ProtonProfile");
353 ProtonProfileMomentum->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
354 ProtonProfileMomentum->GetXaxis()->SetTitle("Momentum, GeV/c");
355 ProtonProfileMomentum->GetYaxis()->SetTitle("dE/dx");
356 ProtonProfileMomentum->SetLineColor(kRed);
357
358 //beta*gamma: for the fit
359 TH1F* ProtonProfileBetaGamma = (TH1F*)hLambdaPBetaGamma->ProfileX("ProtonProfileBetaGamma");
360 if (m_CustomProfile) {
361 ProtonProfileBetaGamma = PrepareProfile(hLambdaPBetaGamma, "ProtonProfileBetaGamma");
362 }
363 ProtonProfileBetaGamma->SetTitle("ProtonProfile");
364 ProtonProfileBetaGamma->GetYaxis()->SetRangeUser(0, m_dedxCutoff);
365 ProtonProfileBetaGamma->GetXaxis()->SetTitle("#beta*#gamma");
366 ProtonProfileBetaGamma->GetYaxis()->SetTitle("dE/dx");
367 ProtonProfileBetaGamma->SetLineColor(kRed);
368
369
370 // for each momentum bin, normalize the pdf
371 hLambdaPMomentum = Normalise2DHisto(hLambdaPMomentum);
372
373 std::unique_ptr<TList> histList(new TList);
374 histList->Add(ProtonProfileMomentum);
375 histList->Add(ProtonProfileBetaGamma);
376 histList->Add(hLambdaPMomentum);
377
378 if (m_isMakePlots) {
379 TFile LambdaHistogrammingPlotFile("SVDdEdxCalibrationLambdaHistogramming.root", "RECREATE");
380 histList->Write();
381 LambdaHistogrammingPlotFile.Close();
382 }
383
384 return histList;
385}

◆ LambdaMassFit()

TTree * LambdaMassFit ( std::shared_ptr< TTree > preselTree)
private

Mass fit for Lambda->ppi.

Definition at line 164 of file SVDdEdxCalibrationAlgorithm.cc.

165{
166 B2INFO("Configuring the Lambda fit...");
167 gROOT->SetBatch(true);
168 RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
169
170 RooRealVar InvM("InvM", "m(p^{+}#pi^{-})", 1.1, 1.13, "GeV/c^{2}");
171
172 RooRealVar ProtonMomentum("ProtonMomentum", "momentum for p", -1.e8, 1.e8);
173 RooRealVar ProtonSVDdEdxTrackMomentum("ProtonSVDdEdxTrackMomentum", "momentum for p", -1.e8, 1.e8);
174 RooRealVar ProtonSVDdEdx("ProtonSVDdEdx", "", -1.e8, 1.e8);
175 RooRealVar ProtonSVDdEdxTrackCosTheta("ProtonSVDdEdxTrackCosTheta", "", -10., 10.);
176
177 RooRealVar exp("exp", "experiment number", 0, 1.e5);
178 RooRealVar run("run", "run number", 0, 1.e7);
179
180 auto variables = new RooArgSet();
181
182 variables->add(InvM);
183
184 variables->add(ProtonMomentum);
185 variables->add(ProtonSVDdEdxTrackMomentum);
186 variables->add(ProtonSVDdEdx);
187 variables->add(ProtonSVDdEdxTrackCosTheta);
188 variables->add(exp);
189 variables->add(run);
190
191 RooDataSet* LambdaDataset = new RooDataSet("LambdaDataset", "LambdaDataset", *variables, Import(*preselTree));
192
193 if (LambdaDataset->sumEntries() == 0) {
194 B2FATAL("The Lambda dataset is empty, stopping here");
195 }
196
197 // the signal PDF; might be revisited at a later point
198
199 RooRealVar GaussMean("GaussMean", " GaussMean", 1.116, 1.111, 1.12);
200 RooRealVar GaussSigma("GaussSigma", "#sigma_{1}", 3.e-3, 3.e-5, 10.e-3);
201 RooGaussian LambdaGauss("LambdaGauss", "LambdaGauss", InvM, GaussMean, GaussSigma);
202
203 /* temporary RooRealVar sigmaBifurGaussL1 and sigmaBifurGaussR1 to replace
204 * RooRealVar resolutionParamL("resolutionParamL", "resolutionParamL", 0.4, 5.e-4, 1.0);
205 * RooRealVar resolutionParamR("resolutionParamR", "resolutionParamR", 0.4, 5.e-4, 1.0);
206 * RooFormulaVar sigmaBifurGaussL1("sigmaBifurGaussL1", "resolutionParamL*GaussSigma", RooArgSet(resolutionParamL, GaussSigma));
207 * RooFormulaVar sigmaBifurGaussR1("sigmaBifurGaussR1", "resolutionParamR*GaussSigma", RooArgSet(resolutionParamR, GaussSigma));
208 */
209 RooRealVar sigmaBifurGaussL1("sigmaBifurGaussL1", "sigma left", 0.4 * 3.e-3, 3.e-5, 10.e-3);
210 RooRealVar sigmaBifurGaussR1("sigmaBifurGaussR1", "sigma right", 0.4 * 3.e-3, 3.e-5, 10.e-3);
211 RooBifurGauss LambdaBifurGauss("LambdaBifurGauss", "LambdaBifurGauss", InvM, GaussMean, sigmaBifurGaussL1, sigmaBifurGaussR1);
212
213 /* temporary RooRealVar sigmaBifurGaussL2 to replace
214 * RooRealVar resolutionParam2("resolutionParam2", "resolutionParam2", 0.2, 5.e-4, 1.0);
215 * sigmaBifurGaussL2("sigmaBifurGaussL2", "resolutionParam2*GaussSigma", RooArgSet(resolutionParam2, GaussSigma));
216 */
217 RooRealVar sigmaBifurGaussL2("sigmaBifurGaussL2", "sigmaBifurGaussL2", 0.2 * 3.e-3, 3.e-5, 10.e-3);
218 RooGaussian LambdaBifurGauss2("LambdaBifurGauss2", "LambdaBifurGauss2", InvM, GaussMean, sigmaBifurGaussL2);
219
220 RooRealVar fracBifurGaussYield("fracBifurGaussYield", "fracBifurGaussYield", 0.3, 5.e-4, 1.0);
221 RooRealVar fracGaussYield("fracGaussYield", "fracGaussYield", 0.8, 5.e-4, 1.0);
222
223 RooAddPdf LambdaCombinedBifurGauss("LambdaCombinedBifurGauss", "LambdaBifurGauss + LambdaBifurGauss2 ", RooArgList(LambdaBifurGauss,
224 LambdaBifurGauss2), RooArgList(fracBifurGaussYield));
225
226 RooAddPdf LambdaSignalPDF("LambdaSignalPDF", "LambdaCombinedBifurGauss + LambdaGauss", RooArgList(LambdaCombinedBifurGauss,
227 LambdaGauss), RooArgList(fracGaussYield));
228
229 // Background PDF
230 RooRealVar BkgPolyCoef0("BkgPolyCoef0", "BkgPolyCoef0", 0.1, 0., 1.5);
231 RooRealVar BkgPolyCoef1("BkgPolyCoef1", "BkgPolyCoef1", -0.5, -1.5, -1.e-3);
232 RooChebychev BkgPolyPDF("BkgPolyPDF", "BkgPolyPDF", InvM, RooArgList(BkgPolyCoef0, BkgPolyCoef1));
233
234 RooRealVar nSignalLambda("nSignalLambda", "nSignalLambda", 0.6 * preselTree->GetEntries(), 0., 0.99 * preselTree->GetEntries());
235 RooRealVar nBkgLambda("nBkgLambda", "nBkgLambda", 0.4 * preselTree->GetEntries(), 0., 0.99 * preselTree->GetEntries());
236 RooAddPdf totalPDFLambda("totalPDFLambda", "totalPDFLambda pdf", RooArgList(LambdaSignalPDF, BkgPolyPDF),
237 RooArgList(nSignalLambda, nBkgLambda));
238
239 B2INFO("Lambda: Start fitting...");
240 RooFitResult* LambdaFitResult = totalPDFLambda.fitTo(*LambdaDataset, Save(kTRUE), PrintLevel(-1));
241
242 int status = LambdaFitResult->status();
243 int covqual = LambdaFitResult->covQual();
244 double diff = nSignalLambda.getValV() + nBkgLambda.getValV() - LambdaDataset->sumEntries();
245
246 B2INFO("Lambda: Fit status: " << status << "; covariance quality: " << covqual);
247 // if the fit is not healthy, try again once before giving up, with a slightly different setup:
248 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalLambda.getError() < sqrt(nSignalLambda.getValV()))
249 || (nSignalLambda.getError() > (nSignalLambda.getValV()))) {
250
251 LambdaFitResult = totalPDFLambda.fitTo(*LambdaDataset, Save(), Strategy(2), Offset(1));
252 status = LambdaFitResult->status();
253 covqual = LambdaFitResult->covQual();
254 diff = nSignalLambda.getValV() + nBkgLambda.getValV() - LambdaDataset->sumEntries();
255 B2INFO("Lambda: updated fit status: " << status << "; covariance quality: " << covqual);
256 }
257
258 if ((status > 0) || (TMath::Abs(diff) > 1.) || (nSignalLambda.getError() < sqrt(nSignalLambda.getValV()))
259 || (nSignalLambda.getError() > (nSignalLambda.getValV()))) {
260 B2WARNING("Lambda: Fit problem: fit status " << status << "; sum of component yields minus the dataset yield is " << diff <<
261 "; signal yield is " << nSignalLambda.getValV() << ", while its uncertainty is " << nSignalLambda.getError());
262 }
263 if (covqual < 2) {
264 B2INFO("Lambda: Fit warning: covariance quality " << covqual);
265 }
266
267 if (m_isMakePlots) {
268 std::unique_ptr<TCanvas> canvLambda(new TCanvas("canvLambda", "canvLambda"));
269 canvLambda->cd();
270 RooPlot* LambdaFitFrame = LambdaDataset->plotOn(InvM.frame(130));
271 totalPDFLambda.plotOn(LambdaFitFrame, LineColor(TColor::GetColor("#4575b4")));
272
273 double chisquare = LambdaFitFrame->chiSquare();
274 B2INFO("Lambda: Fit chi2 = " << chisquare);
275 totalPDFLambda.paramOn(LambdaFitFrame, Layout(0.6, 0.96, 0.93), Format("NEU", AutoPrecision(2)));
276 LambdaFitFrame->getAttText()->SetTextSize(0.03);
277
278 totalPDFLambda.plotOn(LambdaFitFrame, Components("LambdaSignalPDF"), LineColor(TColor::GetColor("#d73027")));
279 totalPDFLambda.plotOn(LambdaFitFrame, Components("BkgPolyPDF"), LineColor(TColor::GetColor("#fc8d59")));
280 totalPDFLambda.plotOn(LambdaFitFrame, LineColor(TColor::GetColor("#4575b4")));
281
282 LambdaFitFrame->GetXaxis()->SetTitle("m(p#pi^{-}) (GeV/c^{2})");
283
284 LambdaFitFrame->Draw();
285
286
287 canvLambda->Print("SVDdEdxCalibrationFitLambda.pdf");
288 TFile LambdaFitPlotFile("SVDdEdxCalibrationLambdaFitPlotFile.root", "RECREATE");
289 canvLambda->Write();
290 LambdaFitPlotFile.Close();
291 }
292 RooStats::SPlot* sPlotDatasetLambda = new RooStats::SPlot("sData", "An SPlot", *LambdaDataset, &totalPDFLambda,
293 RooArgList(nSignalLambda, nBkgLambda));
294
295 for (int iEvt = 0; iEvt < 5; iEvt++) {
296 if (TMath::Abs(sPlotDatasetLambda->GetSWeight(iEvt, "nSignalLambda") + sPlotDatasetLambda->GetSWeight(iEvt,
297 "nBkgLambda") - 1) > 5.e-3)
298 B2FATAL("Lambda: sPlot error: sum of weights not equal to 1");
299 }
300
301 RooDataSet* LambdaDatasetSWeighted = new RooDataSet(LambdaDataset->GetName(), LambdaDataset->GetTitle(), LambdaDataset,
302 *LambdaDataset->get());
303
304 TTree* treeLambdaSWeighted = LambdaDatasetSWeighted->GetClonedTree();
305 treeLambdaSWeighted->SetName("treeLambdaSWeighted");
306
307 B2INFO("Lambda: sPlot done. Proceed to histogramming");
308 return treeLambdaSWeighted;
309}

◆ loadInputJson()

bool loadInputJson ( const std::string & jsonString)
inherited

Load the m_inputJson variable from a string (useful from Python interface). The return bool indicates success or failure.

Definition at line 502 of file CalibrationAlgorithm.cc.

503{
504 try {
505 auto jsonInput = nlohmann::json::parse(jsonString);
506 // Input string has an object (dict) as the top level object?
507 if (jsonInput.is_object()) {
508 m_jsonExecutionInput = jsonInput;
509 return true;
510 } else {
511 B2ERROR("JSON input string isn't an object type i.e. not a '{}' at the top level.");
512 return false;
513 }
514 } catch (nlohmann::json::parse_error&) {
515 B2ERROR("Parsing of JSON input string failed");
516 return false;
517 }
518}
nlohmann::json m_jsonExecutionInput
Optional input JSON object used to make decisions about how to execute the algorithm code.

◆ Normalise2DHisto()

TH2F * Normalise2DHisto ( TH2F * HistoToNormalise)
inlineprivate

Normalise a given dEdx:momentum histogram in each momentum bin, so that sum of entries in each momentum bin is 1.

Note that this accounts for entries in the underflow/overflow bins.

Definition at line 163 of file SVDdEdxCalibrationAlgorithm.h.

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 }

◆ PrepareNewHistogram()

TH2F * PrepareNewHistogram ( TH2F * DataHistogram,
TString NewName,
TF1 * betagamma_function,
TF1 * ResolutionFunctionOriginal,
double bias_correction )
inlineprivate

Generate a new dEdx:momentum histogram from a function that encodes dEdx:momentum trend and a function that encodes dEdx resolution.

Definition at line 191 of file SVDdEdxCalibrationAlgorithm.h.

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 }

◆ PrepareProfile()

TH1F * PrepareProfile ( TH2F * DataHistogram,
TString NewName )
inlineprivate

Reimplement the Profile histogram calculation for a 2D histogram.

The standard ROOT implementation takes the mean Y at a given X, which produces biased results in case of rapidly-rising distributions with non-Gaussian resolutions. We fit in slices to extract the mean more precisely.

Definition at line 232 of file SVDdEdxCalibrationAlgorithm.h.

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 }

◆ resetInputJson()

void resetInputJson ( )
inlineprotectedinherited

Clears the m_inputJson member variable.

Definition at line 330 of file CalibrationAlgorithm.h.

330{m_jsonExecutionInput.clear();}

◆ resetOutputJson()

void resetOutputJson ( )
inlineprotectedinherited

Clears the m_outputJson member variable.

Definition at line 333 of file CalibrationAlgorithm.h.

333{m_jsonExecutionOutput.clear();}

◆ saveCalibration() [1/6]

void saveCalibration ( TClonesArray * data,
const std::string & name )
protectedinherited

Store DBArray payload with given name with default IOV.

Definition at line 297 of file CalibrationAlgorithm.cc.

298{
299 saveCalibration(data, name, m_data.getRequestedIov());
300}

◆ saveCalibration() [2/6]

void saveCalibration ( TClonesArray * data,
const std::string & name,
const IntervalOfValidity & iov )
protectedinherited

Store DBArray with given name and custom IOV.

Definition at line 276 of file CalibrationAlgorithm.cc.

277{
278 B2DEBUG(29, "Saving calibration TClonesArray '" << name << "' to payloads list.");
279 getPayloads().emplace_back(name, data, iov);
280}

◆ saveCalibration() [3/6]

void saveCalibration ( TObject * data)
protectedinherited

Store DB payload with default name and default IOV.

Definition at line 287 of file CalibrationAlgorithm.cc.

288{
289 saveCalibration(data, DataStore::objectName(data->IsA(), ""));
290}
static std::string objectName(const TClass *t, const std::string &name)
Return the storage name for an object of the given TClass and name.
Definition DataStore.cc:150

◆ saveCalibration() [4/6]

void saveCalibration ( TObject * data,
const IntervalOfValidity & iov )
protectedinherited

Store DB payload with default name and custom IOV.

Definition at line 282 of file CalibrationAlgorithm.cc.

283{
284 saveCalibration(data, DataStore::objectName(data->IsA(), ""), iov);
285}

◆ saveCalibration() [5/6]

void saveCalibration ( TObject * data,
const std::string & name )
protectedinherited

Store DB payload with given name with default IOV.

Definition at line 292 of file CalibrationAlgorithm.cc.

293{
294 saveCalibration(data, name, m_data.getRequestedIov());
295}

◆ saveCalibration() [6/6]

void saveCalibration ( TObject * data,
const std::string & name,
const IntervalOfValidity & iov )
protectedinherited

Store DB payload with given name and custom IOV.

Definition at line 270 of file CalibrationAlgorithm.cc.

271{
272 B2DEBUG(29, "Saving calibration TObject = '" << name << "' to payloads list.");
273 getPayloads().emplace_back(name, data, iov);
274}

◆ setCustomProfile()

void setCustomProfile ( bool value = true)
inline

reimplement the profile histogram calculation

Definition at line 80 of file SVDdEdxCalibrationAlgorithm.h.

80{ m_CustomProfile = value; }

◆ setDEdxCutoff()

void setDEdxCutoff ( const double & value)
inline

set the upper edge of the dEdx binning for the payloads

Definition at line 65 of file SVDdEdxCalibrationAlgorithm.h.

65{ m_dedxCutoff = value; }

◆ setDescription()

void setDescription ( const std::string & description)
inlineprotectedinherited

Set algorithm description (in constructor)

Definition at line 321 of file CalibrationAlgorithm.h.

321{m_description = description;}

◆ setFixUnstableFitParameter()

void setFixUnstableFitParameter ( bool value = true)
inline

In the dEdx:betagamma fit, there is one free parameter that makes fit convergence poor.

It is ok to fix it, unless the dEdx behavior changes radically with time.

Definition at line 85 of file SVDdEdxCalibrationAlgorithm.h.

85{ m_FixUnstableFitParameter = value; }

◆ setInputFileNames() [1/2]

void setInputFileNames ( PyObject * inputFileNames)
inherited

Set the input file names used for this algorithm from a Python list.

Set the input file names used for this algorithm and resolve the wildcards.

Definition at line 166 of file CalibrationAlgorithm.cc.

167{
168 // The reasoning for this very 'manual' approach to extending the Python interface
169 // (instead of using boost::python) is down to my fear of putting off final users with
170 // complexity on their side.
171 //
172 // I didn't want users that inherit from this class to be forced to use boost and
173 // to have to define a new python module just to use the CAF. A derived class from
174 // from a boost exposed class would need to have its own boost python module definition
175 // to allow access from a steering file and to the base class functions (I think).
176 // I also couldn't be bothered to write a full framework to get around the issue in a similar
177 // way to Module()...maybe there's an easy way.
178 //
179 // But this way we can allow people to continue using their ROOT implemented classes and inherit
180 // easily from this one. But add in a few helper functions that work with Python objects
181 // created in their steering file i.e. instead of being forced to use STL objects as input
182 // to the algorithm.
183 if (PyList_Check(inputFileNames)) {
184 boost::python::handle<> handle(boost::python::borrowed(inputFileNames));
185 boost::python::list listInputFileNames(handle);
186 auto vecInputFileNames = PyObjConvUtils::convertPythonObject(listInputFileNames, vector<string>());
187 setInputFileNames(vecInputFileNames);
188 } else {
189 B2ERROR("Tried to set the input files but we didn't receive a Python list.");
190 }
191}
void setInputFileNames(PyObject *inputFileNames)
Set the input file names used for this algorithm from a Python list.
Scalar convertPythonObject(const boost::python::object &pyObject, Scalar)
Convert from Python to given type.

◆ setInputFileNames() [2/2]

void setInputFileNames ( std::vector< std::string > inputFileNames)
protectedinherited

Set the input file names used for this algorithm.

Set the input file names used for this algorithm and resolve the wildcards.

Definition at line 194 of file CalibrationAlgorithm.cc.

195{
196 // A lot of code below is tweaked from RootInputModule::initialize,
197 // since we're basically copying the functionality anyway.
198 if (inputFileNames.empty()) {
199 B2WARNING("You have called setInputFileNames() with an empty list. Did you mean to do that?");
200 return;
201 }
202 auto tmpInputFileNames = RootIOUtilities::expandWordExpansions(inputFileNames);
203
204 // We'll use a set to enforce sorted unique file paths as we check them
205 set<string> setInputFileNames;
206 // Check that files exist and convert to absolute paths
207 for (auto path : tmpInputFileNames) {
208 string fullPath = fs::absolute(path).string();
209 if (fs::exists(fullPath)) {
210 setInputFileNames.insert(fs::canonical(fullPath).string());
211 } else {
212 B2WARNING("Couldn't find the file " << path);
213 }
214 }
215
216 if (setInputFileNames.empty()) {
217 B2WARNING("No valid files specified!");
218 return;
219 } else {
220 // Reset the run -> files map as our files are likely different
221 m_runsToInputFiles.clear();
222 }
223
224 // Open TFile to check they can be accessed by ROOT
225 TDirectory* dir = gDirectory;
226 for (const string& fileName : setInputFileNames) {
227 unique_ptr<TFile> f;
228 try {
229 f.reset(TFile::Open(fileName.c_str(), "READ"));
230 } catch (logic_error&) {
231 //this might happen for ~invaliduser/foo.root
232 //actually undefined behaviour per standard, reported as ROOT-8490 in JIRA
233 }
234 if (!f || !f->IsOpen()) {
235 B2FATAL("Couldn't open input file " + fileName);
236 }
237 }
238 dir->cd();
239
240 // Copy the entries of the set to a vector
241 m_inputFileNames = vector<string>(setInputFileNames.begin(), setInputFileNames.end());
243}
std::string m_granularityOfData
Granularity of input data. This only changes when the input files change so it isn't specific to an e...
std::string getGranularityFromData() const
Get the granularity of collected data.
std::vector< std::string > expandWordExpansions(const std::vector< std::string > &filenames)
Performs wildcard expansion using wordexp(), returns matches.

◆ setMinEvtsPerTree()

void setMinEvtsPerTree ( const double & value)
inline

set the upper edge of the dEdx binning for the payloads

Definition at line 70 of file SVDdEdxCalibrationAlgorithm.h.

70{ m_MinEvtsPerTree = value; }

◆ setMonitoringPlots()

void setMonitoringPlots ( bool value = false)
inline

function to enable plotting

Definition at line 45 of file SVDdEdxCalibrationAlgorithm.h.

45{ m_isMakePlots = value; }

◆ setNEventsToGenerate()

void setNEventsToGenerate ( const int & value)
inline

set the number of events to generate, per momentum bin, for the payloads

Definition at line 75 of file SVDdEdxCalibrationAlgorithm.h.

75{ m_NToGenerate = value; }

◆ setNumBGBins()

void setNumBGBins ( const int & value)
inline

set the number of beta*gamma bins for the fits

Definition at line 60 of file SVDdEdxCalibrationAlgorithm.h.

60{ m_numBGBins = value; }

◆ setNumDEdxBins()

void setNumDEdxBins ( const int & value)
inline

set the number of dEdx bins for the payloads

Definition at line 50 of file SVDdEdxCalibrationAlgorithm.h.

50{ m_numDEdxBins = value; }

◆ setNumPBins()

void setNumPBins ( const int & value)
inline

set the number of momentum bins for the payloads

Definition at line 55 of file SVDdEdxCalibrationAlgorithm.h.

55{ m_numPBins = value; }

◆ setOutputJsonValue()

template<class T>
void setOutputJsonValue ( const std::string & key,
const T & value )
inlineprotectedinherited

Set a key:value pair for the outputJson object, expected to used internally during calibrate()

Definition at line 337 of file CalibrationAlgorithm.h.

337{m_jsonExecutionOutput[key] = value;}

◆ setPrefix()

void setPrefix ( const std::string & prefix)
inlineinherited

Set the prefix used to identify datastore objects.

Definition at line 167 of file CalibrationAlgorithm.h.

167{m_prefix = prefix;}

◆ setUsePionBGFunctionForEverything()

void setUsePionBGFunctionForEverything ( bool value = false)
inline

use the pion beta*gamma function for other hadrons

Definition at line 90 of file SVDdEdxCalibrationAlgorithm.h.

90{ m_UsePionBGFunctionForEverything = value; }

◆ setUseProtonBGFunctionForEverything()

void setUseProtonBGFunctionForEverything ( bool value = false)
inline

use the proton beta*gamma function for other hadrons

Definition at line 95 of file SVDdEdxCalibrationAlgorithm.h.

95{ m_UseProtonBGFunctionForEverything = value; }

◆ updateDBObjPtrs()

void updateDBObjPtrs ( const unsigned int event = 1,
const int run = 0,
const int experiment = 0 )
protectedinherited

Updates any DBObjPtrs by calling update(event) for DBStore.

Definition at line 404 of file CalibrationAlgorithm.cc.

405{
406 // Construct an EventMetaData object but NOT in the Datastore
407 EventMetaData emd(event, run, experiment);
408 // Explicitly update while avoiding registering a Datastore object
410 // Also update the intra-run objects to the event at the same time (maybe unnecessary...)
412}
static DBStore & Instance()
Instance of a singleton DBStore.
Definition DBStore.cc:26
void updateEvent()
Updates all intra-run dependent objects.
Definition DBStore.cc:140
void update()
Updates all objects that are outside their interval of validity.
Definition DBStore.cc:77

Member Data Documentation

◆ m_allExpRun

const ExpRun m_allExpRun = make_pair(-1, -1)
staticprivateinherited

allExpRun

Definition at line 364 of file CalibrationAlgorithm.h.

◆ m_boundaries

std::vector<Calibration::ExpRun> m_boundaries
protectedinherited

When using the boundaries functionality from isBoundaryRequired, this is used to store the boundaries. It is cleared when.

Definition at line 261 of file CalibrationAlgorithm.h.

◆ m_CustomProfile

bool m_CustomProfile = 1
private

reimplement profile histogram calculation instead of the ROOT implementation?

Definition at line 122 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_data

ExecutionData m_data
privateinherited

Data specific to a SINGLE execution of the algorithm. Gets reset at the beginning of execution.

Definition at line 382 of file CalibrationAlgorithm.h.

◆ m_dedxCutoff

double m_dedxCutoff = 5.e6
private

the upper edge of the dEdx binning for the payloads

Definition at line 116 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_dedxMaxPossible

double m_dedxMaxPossible = 7.e6
private

the approximate max possible value of dEdx

Definition at line 117 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_description

std::string m_description {""}
privateinherited

Description of the algorithm.

Definition at line 385 of file CalibrationAlgorithm.h.

385{""};

◆ m_DeuteronPDGMass

const double m_DeuteronPDGMass = TDatabasePDG::Instance()->GetParticle(1000010020)->Mass()
private

PDG mass for the deuteron.

Definition at line 135 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_ElectronPDGMass

const double m_ElectronPDGMass = TDatabasePDG::Instance()->GetParticle(11)->Mass()
private

PDG mass for the electron.

Definition at line 130 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_FixUnstableFitParameter

bool m_FixUnstableFitParameter
private
Initial value:
=
1

In the dEdx:betagamma fit, there is one free parameter that makes fit convergence poor.

It is ok to fix it, unless the dEdx behavior changes radically with time.

Definition at line 127 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_granularityOfData

std::string m_granularityOfData
privateinherited

Granularity of input data. This only changes when the input files change so it isn't specific to an execution.

Definition at line 379 of file CalibrationAlgorithm.h.

◆ m_inputFileNames

std::vector<std::string> m_inputFileNames
privateinherited

List of input files to the Algorithm, will initially be user defined but then gets the wildcards expanded during execute()

Definition at line 373 of file CalibrationAlgorithm.h.

◆ m_isMakePlots

bool m_isMakePlots
private

produce plots for monitoring

Definition at line 104 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_jsonExecutionInput

nlohmann::json m_jsonExecutionInput = nlohmann::json::object()
privateinherited

Optional input JSON object used to make decisions about how to execute the algorithm code.

Definition at line 397 of file CalibrationAlgorithm.h.

◆ m_jsonExecutionOutput

nlohmann::json m_jsonExecutionOutput = nlohmann::json::object()
privateinherited

Optional output JSON object that can be set during the execution by the underlying algorithm code.

Definition at line 403 of file CalibrationAlgorithm.h.

◆ m_KaonPDGMass

const double m_KaonPDGMass = TDatabasePDG::Instance()->GetParticle(321)->Mass()
private

PDG mass for the charged kaon.

Definition at line 133 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_MinEvtsPerTree

int m_MinEvtsPerTree
private
Initial value:
=
100

number of events in TTree below which we don't try to fit

Definition at line 118 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_MuonPDGMass

const double m_MuonPDGMass = TDatabasePDG::Instance()->GetParticle(13)->Mass()
private

PDG mass for the muon.

Definition at line 131 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_NToGenerate

int m_NToGenerate
private
Initial value:
=
500000

the number of events to be generated in each momentum bin in the new payloads

Definition at line 120 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_numBGBins

int m_numBGBins = 69
private

the number of beta*gamma bins for the profile and fitting

Definition at line 115 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_numDEdxBins

int m_numDEdxBins = 100
private

the number of dEdx bins for the payloads

Definition at line 113 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_numPBins

int m_numPBins = 69
private

the number of momentum bins for the payloads

Definition at line 114 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_PionPDGMass

const double m_PionPDGMass = TDatabasePDG::Instance()->GetParticle(211)->Mass()
private

PDG mass for the charged pion.

Definition at line 132 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_prefix

std::string m_prefix {""}
privateinherited

The name of the TDirectory the collector objects are contained within.

Definition at line 388 of file CalibrationAlgorithm.h.

388{""};

◆ m_ProtonPDGMass

const double m_ProtonPDGMass = TDatabasePDG::Instance()->GetParticle(2212)->Mass()
private

PDG mass for the proton.

Definition at line 134 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_runsToInputFiles

std::map<Calibration::ExpRun, std::vector<std::string> > m_runsToInputFiles
privateinherited

Map of Runs to input files. Gets filled when you call getRunRangeFromAllData, gets cleared when setting input files again.

Definition at line 376 of file CalibrationAlgorithm.h.

◆ m_UsePionBGFunctionForEverything

bool m_UsePionBGFunctionForEverything
private
Initial value:
=
0

Assume that the dEdx:betagamma trend is the same for all hadrons; use the pion trend as representative.

Definition at line 123 of file SVDdEdxCalibrationAlgorithm.h.

◆ m_UseProtonBGFunctionForEverything

bool m_UseProtonBGFunctionForEverything
private
Initial value:
=
0

Assume that the dEdx:betagamma trend is the same for all hadrons; use the proton trend as representative.

Definition at line 125 of file SVDdEdxCalibrationAlgorithm.h.


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