Belle II Software  release-05-01-25
eclBhabhaTAlgorithm Class Reference

Calibrate ecl crystals using bhabha events. More...

#include <eclBhabhaTAlgorithm.h>

Inheritance diagram for eclBhabhaTAlgorithm:
Collaboration diagram for eclBhabhaTAlgorithm:

Public Types

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

Public Member Functions

 eclBhabhaTAlgorithm ()
 ..Constructor
 
virtual ~eclBhabhaTAlgorithm ()
 ..Destructor
 
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. More...
 
Calibration::ExpRun convertPyExpRun (PyObject *pyObj)
 Performs the conversion of PyObject to ExpRun. More...
 
std::string getCollectorName () const
 Alias for prefix. More...
 
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. More...
 
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. More...
 
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 algoithm (set by developers in constructor)
 
bool loadInputJson (const std::string &jsonString)
 Load the m_inputJson variable from a string (useful from Python interface). The rturn 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>.
 

Public Attributes

int cellIDLo
 Fit crystals with cellID0 in the inclusive range [cellIDLo,cellIDHi].
 
int cellIDHi
 Fit crystals with cellID0 in the inclusive range [cellIDLo,cellIDHi].
 
double meanCleanRebinFactor
 Rebinning factor for mean calculation.
 
double meanCleanCutMinFactor
 After rebinning, create a mask for bins that have values less than meanCleanCutMinFactor times the maximum bin value. More...
 
int crateIDLo
 Fit crates with crateID0 in the inclusive range [crateIDLo,crateIDHi].
 
int crateIDHi
 Fit crates with crateID0 in the inclusive range [crateIDLo,crateIDHi].
 
bool debugOutput
 Save every histogram and fitted function to debugFilename.
 
std::string debugFilenameBase
 Name of file with debug output, eclBhabhaTAlgorithm.root by default.
 
std::string collectorName
 Name of the collector.
 

Protected Member Functions

virtual EResult calibrate ()
 ..Run algorithm on events More...
 
void setInputFileNames (std::vector< std::string > inputFileNames)
 Set the input file names used for this algorithm. More...
 
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 interally 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
 

Private Member Functions

std::string getExpRunString (Calibration::ExpRun &expRun) const
 Gets the "exp.run" string repr. of (exp,run)
 
std::string getFullObjectPath (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

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)
 

Detailed Description

Calibrate ecl crystals using bhabha events.

Definition at line 39 of file eclBhabhaTAlgorithm.h.

Member Enumeration Documentation

◆ EResult

enum EResult
inherited

The result of calibration.

Enumerator
c_OK 

Finished successfuly =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 50 of file CalibrationAlgorithm.h.

Member Function Documentation

◆ calibrate()

CalibrationAlgorithm::EResult calibrate ( )
protectedvirtual

..Run algorithm on events

Put root into batch mode so that we don't try to open a graphics window

Write out job parameters


Implements CalibrationAlgorithm.

Definition at line 36 of file eclBhabhaTAlgorithm.cc.

37 {
39  gROOT->SetBatch();
40 
41 
43  B2INFO("eclBhabhaTAlgorithm parameters:");
44  B2INFO("cellIDLo = " << cellIDLo);
45  B2INFO("cellIDHi = " << cellIDHi);
46  B2INFO("meanCleanRebinFactor = " << meanCleanRebinFactor);
47  B2INFO("meanCleanCutMinFactor = " << meanCleanCutMinFactor);
48  B2INFO("crateIDLo = " << crateIDLo);
49  B2INFO("crateIDHi = " << crateIDHi);
50 
51 
52  /* Histogram with the data collected by eclBhabhaTCollectorModule */
53 
54  auto TimevsCrysPrevCrateCalibPrevCrystCalib = getObjectPtr<TH2F>("TimevsCrysPrevCrateCalibPrevCrystCalib");
55  auto TimevsCratePrevCrateCalibPrevCrystCalib = getObjectPtr<TH2F>("TimevsCratePrevCrateCalibPrevCrystCalib");
56  auto TimevsCrysNoCalibrations = getObjectPtr<TH2F>("TimevsCrysNoCalibrations");
57  auto TimevsCrysPrevCrateCalibNoCrystCalib = getObjectPtr<TH2F>("TimevsCrysPrevCrateCalibNoCrystCalib");
58  auto TimevsCrateNoCrateCalibPrevCrystCalib = getObjectPtr<TH2F>("TimevsCrateNoCrateCalibPrevCrystCalib");
59 
60  // Collect other plots just for reference - combines all the runs for these plots.
61  auto cutflow = getObjectPtr<TH1F>("cutflow");
62 
63 
64  if (!TimevsCrysNoCalibrations) return c_Failure;
65 
68  // Make tool for mapping ecl crystal to other ecl objects
69  // e.g. crates, shapers, etc.
70  // Need to load information about the event/run/experiment to get the right database information
71  int expNumberForCrates = -1;
72  int runNumberForCrates = -1;
73  int eventNumberForCrates = 1;
74  for (auto expRun : getRunList()) {
75  expNumberForCrates = expRun.first;
76  runNumberForCrates = expRun.second;
77  }
78  updateDBObjPtrs(eventNumberForCrates, runNumberForCrates, expNumberForCrates);
79  unique_ptr<ECLChannelMapper> crystalMapper(new ECL::ECLChannelMapper());
80  crystalMapper->initFromDB();
81 
82 
83  TFile* histfile = 0;
84  unique_ptr<TTree> tree_crystal(new TTree("tree_crystal", "Debug data from bhabha time calibration algorithm for crystals"));
85 
86  unique_ptr<TTree> tree_crate(new TTree("tree_crate", "Debug data from bhabha time calibration algorithm for crates"));
87  int tree_cid;
88 
89  // Vector of time offsets to be saved in the database.
90  vector<float> t_offsets;
91  // Vector of time offset uncertainties to be saved in the database.
92  vector<float> t_offsets_unc;
93 
94  vector<float> t_offsets_prev; // previous time offsets
95 
96  // temp copies of offsets
97  vector<float> t_offsets_temp;
98  vector<float> t_offsets_unc_temp;
99 
100 
101  int minNumEntries = 40;
102 
103 
104  double mean = 0;
105  double sigma = -1;
106  double mean_unc = 0;
107  int crystalCalibSaved = 0;
108  double tsPrev = 0;
109 
110 
111  bool minRunNumBool = false;
112  bool maxRunNumBool = false;
113  int minRunNum = -1;
114  int maxRunNum = -1;
115  int minExpNum = -1;
116  int maxExpNum = -1;
117  for (auto expRun : getRunList()) {
118  int expNumber = expRun.first;
119  int runNumber = expRun.second;
120  if (!minRunNumBool) {
121  minExpNum = expNumber;
122  minRunNum = runNumber;
123  minRunNumBool = true;
124  }
125  if (!maxRunNumBool) {
126  maxExpNum = expNumber;
127  maxRunNum = runNumber;
128  maxRunNumBool = true;
129  }
130  if (((minRunNum > runNumber) && (minExpNum >= expNumber)) ||
131  (minExpNum > expNumber)) {
132  minExpNum = expNumber;
133  minRunNum = runNumber;
134  }
135  if (((maxRunNum < runNumber) && (maxExpNum <= expNumber)) ||
136  (maxExpNum < expNumber))
137 
138  {
139  maxExpNum = expNumber;
140  maxRunNum = runNumber;
141  }
142  }
143 
144  B2INFO("debugFilenameBase = " << debugFilenameBase);
145  string runNumsString = string("_") + to_string(minExpNum) + "_" + to_string(minRunNum) + string("-") + to_string(
146  maxExpNum) + "_" + to_string(maxRunNum);
147  string debugFilename = debugFilenameBase + runNumsString + string(".root");
148 
149 
150  B2INFO("Debug output rootfile: " << debugFilename);
151  histfile = new TFile(debugFilename.c_str(), "recreate");
152 
153 
154  TimevsCrysPrevCrateCalibPrevCrystCalib ->Write();
155  TimevsCratePrevCrateCalibPrevCrystCalib->Write();
156  TimevsCrysNoCalibrations ->Write();
157  TimevsCrysPrevCrateCalibNoCrystCalib ->Write();
158  TimevsCrateNoCrateCalibPrevCrystCalib ->Write();
159 
160  cutflow->Write();
161 
162 
163  if (debugOutput) {
164  tree_crystal->Branch("cid", &tree_cid)->SetTitle("Cell ID, 1..8736");
165  tree_crystal->Branch("ts", &mean)->SetTitle("Time offset mean, ts, ns");
166  tree_crystal->Branch("tsUnc", &mean_unc)->SetTitle("Error of time ts mean, ns.");
167  tree_crystal->Branch("tsSigma", &sigma)->SetTitle("Sigma of time ts distribution, ns");
168  tree_crystal->Branch("crystalCalibSaved",
169  &crystalCalibSaved)->SetTitle("0=crystal skipped, 1=crystal calib saved (num entries based)");
170  tree_crystal->Branch("tsPrev", &tsPrev)->SetTitle("Previous crystal time offset, ts, ns");
171  tree_crystal->SetAutoSave(10);
172  }
173 
174 
175  double hist_tmin = TimevsCrysNoCalibrations->GetYaxis()->GetXmin();
176  double hist_tmax = TimevsCrysNoCalibrations->GetYaxis()->GetXmax();
177 
178  double time_fit_min = hist_tmax; // Set min value to largest possible value so that it gets reset
179  double time_fit_max = hist_tmin; // Set max value to smallest possible value so that it gets reset
180 
181  B2INFO("hist_tmin = " << hist_tmin);
182  B2INFO("hist_tmax = " << hist_tmax);
183 
184 
185  const double TICKS_TO_NS = 1.0 / (4.0 * EclConfiguration::m_rf) *
186  1e3; // 1/(4fRF) = 0.4913 ns/clock tick, where fRF is the accelerator RF frequency, fRF=508.889 MHz. Same for all crystals. Proper accurate value
187 
188  // The ts and tcrate database values are filled once per tcol instance so count the number of times that the database values
189  // were summed together by the histogram merging process and extract out the original values again.
190  auto databaseCounter = getObjectPtr<TH1F>("databaseCounter");
191  float numTimesFilled = databaseCounter->GetBinContent(1);
192  B2INFO("Number of times database histograms were merged = " << numTimesFilled);
193 
194 
195  auto TsDatabase = getObjectPtr<TH1F>("TsDatabase");
196  auto TsDatabaseUnc = getObjectPtr<TH1F>("TsDatabaseUnc");
197  for (int i = 1; i <= 8736; i++) {
198  t_offsets.push_back(TsDatabase->GetBinContent(i - 1) / numTimesFilled);
199  t_offsets_prev.push_back(TsDatabase->GetBinContent(i) / numTimesFilled);
200 
201  B2INFO(" TsDatabase->GetBinContent(i) = " << TsDatabase->GetBinContent(i));
202  B2INFO("t_offsets_prev = " << t_offsets_prev[i - 1]);
203 
204 
205  t_offsets_unc.push_back(TsDatabaseUnc->GetBinContent(i - 1) / numTimesFilled);
206  }
207 
208 
209  for (int crys_id = cellIDLo; crys_id <= cellIDHi; crys_id++) {
210  mean = 0;
211  mean_unc = -1;
212  sigma = -1;
213  crystalCalibSaved = 0;
214 
215  double database_mean = 0;
216  double database_mean_unc = 0;
217 
218  B2INFO("Crystal id = " << crys_id);
219 
220 
221  vector<bool> ts_new_was_set(8736, false);
222 
223 
224  /* Determining which bins to mask out for mean calculation
225  */
226 
227  TH1D* h_time = TimevsCrysPrevCrateCalibNoCrystCalib->ProjectionY("h_time_psi", crys_id, crys_id);
228  TH1D* h_timeMask = (TH1D*)h_time->Clone();
229  TH1D* h_timeMasked = (TH1D*)h_time->Clone();
230  TH1D* h_timeRebin = (TH1D*)h_time->Clone();
231 
232  // Do rebinning and cleaning of some bins but only if the user selection values call for it since it slows the code down
233  if (meanCleanRebinFactor != 1 || meanCleanCutMinFactor != 1) {
234 
235  h_timeRebin->Rebin(meanCleanRebinFactor);
236 
237  h_timeMask->Scale(0.0); // set all bins to being masked off
238 
239  time_fit_min = hist_tmax; // Set min value to largest possible value so that it gets reset
240  time_fit_max = hist_tmin; // Set max value to smallest possible value so that it gets reset
241 
242  // Find value of bin with max value
243  double histRebin_max = h_timeRebin->GetMaximum();
244 
245  bool maskedOutNonZeroBin = false;
246  // Loop over all bins to find those with content less than a certain threshold. Mask the non-rebinned histogram for the corresponding bins
247  for (int bin = 1; bin <= h_timeRebin->GetNbinsX(); bin++) {
248  for (int rebinCounter = 1; rebinCounter <= meanCleanRebinFactor; rebinCounter++) {
249  int nonRebinnedBinNumber = (bin - 1) * meanCleanRebinFactor + rebinCounter;
250  if (nonRebinnedBinNumber < h_time->GetNbinsX()) {
251  if (h_timeRebin->GetBinContent(bin) >= histRebin_max * meanCleanCutMinFactor) {
252  h_timeMask->SetBinContent(nonRebinnedBinNumber, 1);
253 
254  // Save the lower and upper edges of the rebin histogram time range for fitting purposes
255  double x_lower = h_timeRebin->GetXaxis()->GetBinLowEdge(bin);
256  double x_upper = h_timeRebin->GetXaxis()->GetBinUpEdge(bin);
257  if (x_lower < time_fit_min) {
258  time_fit_min = x_lower;
259  }
260  if (x_upper > time_fit_max) {
261  time_fit_max = x_upper;
262  }
263 
264  } else {
265  if (h_time->GetBinContent(nonRebinnedBinNumber) > 0) {
266  B2INFO("Setting bin " << nonRebinnedBinNumber << " from " << h_timeMasked->GetBinContent(nonRebinnedBinNumber) << " to 0");
267  maskedOutNonZeroBin = true;
268  }
269  h_timeMasked->SetBinContent(nonRebinnedBinNumber, 0);
270  }
271  }
272  }
273  }
274  B2INFO("Bins with non-zero values have been masked out: " << maskedOutNonZeroBin);
275  h_timeMasked->ResetStats();
276  h_timeMask->ResetStats();
277 
278  }
279 
280  // Calculate mean from masked histogram
281  double default_meanMasked = h_timeMasked->GetMean();
282  //double default_meanMasked_unc = h_timeMasked->GetMeanError();
283  B2INFO("default_meanMasked = " << default_meanMasked);
284 
285 
286  // Get the overall mean and standard deviation of the distribution within the plot. This doesn't require a fit.
287  double default_mean = h_time->GetMean();
288  double default_mean_unc = h_time->GetMeanError();
289  double default_sigma = h_time->GetStdDev();
290 
291  B2INFO("Fitting crystal between " << time_fit_min << " and " << time_fit_max);
292 
293  // gaus(0) is a substitute for [0]*exp(-0.5*((x-[1])/[2])**2)
294  TF1* gaus = new TF1("func", "gaus(0)", time_fit_min, time_fit_max);
295  gaus->SetParNames("numCrystalHitsNormalization", "mean", "sigma");
296  /*
297  gaus->ReleaseParameter(0); // number of crystals
298  gaus->ReleaseParameter(1); // mean
299  gaus->ReleaseParameter(2); // standard deviation
300  */
301 
302  double hist_max = h_time->GetMaximum();
303 
304  //=== Estimate initial value of sigma as std dev.
305  double stddev = h_time->GetStdDev();
306  sigma = stddev;
307  mean = default_mean;
308 
309  //=== Setting parameters for initial iteration
310  gaus->SetParameter(0, hist_max / 2.);
311  gaus->SetParameter(1, mean);
312  gaus->SetParameter(2, sigma);
313  // L -- Use log likelihood method
314  // I -- Use integral of function in bin instead of value at bin center // not using
315  // R -- Use the range specified in the function range
316  // B -- Fix one or more parameters with predefined function // not using
317  // Q -- Quiet mode
318 
319  h_timeMasked->Fit(gaus, "LQR"); // L for likelihood, R for x-range, Q for fit quiet mode
320 
321  double fit_mean = gaus->GetParameter(1);
322  double fit_mean_unc = gaus->GetParError(1);
323  double fit_sigma = gaus->GetParameter(2);
324 
325  double meanDiff = fit_mean - default_mean;
326  double meanUncDiff = fit_mean_unc - default_mean_unc;
327  double sigmaDiff = fit_sigma - default_sigma;
328 
329  bool good_fit = false;
330 
331  if ((fabs(meanDiff) > 10) ||
332  (fabs(meanUncDiff) > 10) ||
333  (fabs(sigmaDiff) > 10) ||
334  (fit_mean_unc > 3) ||
335  (fit_sigma < 0.1) ||
336  (fit_mean < time_fit_min) ||
337  (fit_mean > time_fit_max)) {
338  B2INFO("Crystal id = " << crys_id);
339  B2INFO("fit mean, default mean = " << fit_mean << ", " << default_mean);
340  B2INFO("fit mean unc, default mean unc = " << fit_mean_unc << ", " << default_mean_unc);
341  B2INFO("fit sigma, default sigma = " << fit_sigma << ", " << default_sigma);
342 
343  B2INFO("crystal fit mean - hist mean = " << meanDiff);
344  B2INFO("fit mean unc. - hist mean unc. = " << meanUncDiff);
345  B2INFO("fit sigma - hist sigma = " << sigmaDiff);
346 
347  B2INFO("fit_mean = " << fit_mean);
348  B2INFO("time_fit_min = " << time_fit_min);
349  B2INFO("time_fit_max = " << time_fit_max);
350 
351  if (fabs(meanDiff) > 10) B2INFO("fit mean diff too large");
352  if (fabs(meanUncDiff) > 10) B2INFO("fit mean unc diff too large");
353  if (fabs(sigmaDiff) > 10) B2INFO("fit mean sigma diff too large");
354  if (fit_mean_unc > 3) B2INFO("fit mean unc too large");
355  if (fit_sigma < 0.1) B2INFO("fit sigma too small");
356 
357  } else {
358  good_fit = true;
359  }
360 
361 
362 
363  // Set the tree_crystal values - ignore fit values !!!!!!!!!!!!!!!!
364  sigma = default_sigma;
365 
366 
367  int numEntries = h_time->GetEntries();
368  // If number of entries in histogram is greater than X then use the statistical information from the data otherwise leave crystal uncalibrated. Histograms are still shown though.
369  // ALSO require the that fits are good.
370  if ((numEntries >= minNumEntries) && good_fit) {
371  crystalCalibSaved = 1;
372  database_mean = fit_mean;
373  database_mean_unc = fit_mean_unc;
374  } else {
375  database_mean = default_mean;
376  database_mean_unc = default_mean_unc;
377  }
378 
379  if (numEntries < minNumEntries) B2INFO("Number of entries less than minimum");
380  if (numEntries == 0) B2INFO("Number of entries == 0");
381 
382 
383  tree_cid = crys_id;
384 
385  // For the database, convert back from ns to ADC ticks.
386  t_offsets[crys_id - 1] = database_mean / TICKS_TO_NS;
387  t_offsets_unc[crys_id - 1] = database_mean_unc / TICKS_TO_NS;
388 
389 
390  histfile->WriteTObject(h_time, (std::string("h_time_psi") + std::to_string(crys_id)).c_str());
391  histfile->WriteTObject(h_timeMasked, (std::string("h_time_psi_masked") + std::to_string(crys_id)).c_str());
392 
393  mean = database_mean;
394  mean_unc = database_mean_unc;
395 
396  tsPrev = t_offsets_prev[crys_id - 1] * TICKS_TO_NS;
397 
398  delete gaus;
399  tree_crystal->Fill();
400  }
401 
402  ECLCrystalCalib* BhabhaTCalib = new ECLCrystalCalib();
403  BhabhaTCalib->setCalibVector(t_offsets, t_offsets_unc);
404 
405 
406  // Save the information to the payload if there is at least one crystal
407  // begin calibrated.
408  if (cellIDLo <= cellIDHi) {
409  saveCalibration(BhabhaTCalib, "ECLCrystalTimeOffset");
410  B2DEBUG(30, "crystal payload made");
411  }
412 
413 
414  B2DEBUG(30, "end of crystal start of crate corrections .....");
415 
416 
417  /* CRATE CORRECTIONS */
418 
419  hist_tmin = TimevsCrateNoCrateCalibPrevCrystCalib->GetYaxis()->GetXmin();
420  hist_tmax = TimevsCrateNoCrateCalibPrevCrystCalib->GetYaxis()->GetXmax();
421 
422  B2DEBUG(30, "Found min/max of X axis of TimevsCrateNoCrateCalibPrevCrystCalib");
423 
424  // Vector of time offsets to be saved in the database.
425 
426  auto TcrateDatabase = getObjectPtr<TH1F>("TcrateDatabase");
427 
428 
429  B2DEBUG(30, "Retrieved Ts and Tcrate histograms from tcol root file");
430 
431  double database_mean_crate = 0;
432  double database_mean_crate_unc = 0;
433 
434  double mean_crate = 0;
435  double sigma_crate = -1;
436 
437  vector<float> tcrate_mean_new(52, 0.0);
438  vector<float> tcrate_mean_unc_new(52, 0.0);
439  vector<float> tcrate_sigma_new(52, 0.0);
440  vector<float> tcrate_mean_prev(52, 0.0);
441  vector<float> tcrate_sigma_prev(52, 0.0);
442  vector<bool> tcrate_new_was_set(52, false);
443 
444  B2DEBUG(30, "crate vectors set");
445 
446 
447 
448  // Crate time calibration constants are saved per crystal so read them per crystal
449  // and save as one entry per crate in the array
450  for (int crys_id = 1; crys_id <= 8736; crys_id++) {
451  int crate_id_from_crystal = crystalMapper->getCrateID(crys_id);
452  tcrate_mean_prev[crate_id_from_crystal - 1] = TcrateDatabase->GetBinContent(crys_id) / numTimesFilled;
453  }
454 
455  for (int crate_id = crateIDLo; crate_id <= crateIDHi; crate_id++) {
456 
457  B2DEBUG(30, "Start of crate id = " << crate_id);
458 
459  TH1D* h_time_crate = TimevsCrateNoCrateCalibPrevCrystCalib->ProjectionY("h_time_psi_crate", crate_id, crate_id);
460  TH1D* h_time_crate_mask = (TH1D*)h_time_crate->Clone();
461  TH1D* h_time_crate_masked = (TH1D*)h_time_crate->Clone();
462  TH1D* h_time_crate_rebin = (TH1D*)h_time_crate->Clone();
463 
464 
465 
466  // Do rebinning and cleaning of some bins but only if the user selection values call for it since it slows the code down
467  if (meanCleanRebinFactor != 1 || meanCleanCutMinFactor != 1) {
468 
469  h_time_crate_rebin->Rebin(meanCleanRebinFactor);
470  h_time_crate_mask->Scale(0.0); // set all bins to being masked off
471 
472  time_fit_min = hist_tmax; // Set min value to largest possible value so that it gets reset
473  time_fit_max = hist_tmin; // Set max value to smallest possible value so that it gets reset
474 
475  // Find value of bin with max value
476  double histRebin_max = h_time_crate_rebin->GetMaximum();
477 
478  bool maskedOutNonZeroBin = false;
479  // Loop over all bins to find those with content less than a certain threshold. Mask the non-rebinned histogram for the corresponding bins
480  for (int bin = 1; bin <= h_time_crate_rebin->GetNbinsX(); bin++) {
481  for (int rebinCounter = 1; rebinCounter <= meanCleanRebinFactor; rebinCounter++) {
482  int nonRebinnedBinNumber = (bin - 1) * meanCleanRebinFactor + rebinCounter;
483  if (nonRebinnedBinNumber < h_time_crate->GetNbinsX()) {
484  if (h_time_crate_rebin->GetBinContent(bin) >= histRebin_max * meanCleanCutMinFactor) {
485  h_time_crate_mask->SetBinContent(nonRebinnedBinNumber, 1);
486 
487  // Save the lower and upper edges of the rebin histogram time range for fitting purposes
488  double x_lower = h_time_crate_rebin->GetXaxis()->GetBinLowEdge(bin);
489  double x_upper = h_time_crate_rebin->GetXaxis()->GetBinUpEdge(bin);
490  if (x_lower < time_fit_min) {
491  time_fit_min = x_lower;
492  }
493  if (x_upper > time_fit_max) {
494  time_fit_max = x_upper;
495  }
496  } else {
497  if (h_time_crate->GetBinContent(nonRebinnedBinNumber) > 0) {
498  B2INFO("Setting bin " << nonRebinnedBinNumber << " from " << h_time_crate_masked->GetBinContent(nonRebinnedBinNumber) << " to 0");
499  maskedOutNonZeroBin = true;
500  }
501  h_time_crate_masked->SetBinContent(nonRebinnedBinNumber, 0);
502  }
503  }
504  }
505  }
506  B2INFO("Bins with non-zero values have been masked out: " << maskedOutNonZeroBin);
507  h_time_crate_masked->ResetStats();
508  h_time_crate_mask->ResetStats();
509 
510  }
511 
512 
513 
514  B2DEBUG(30, "crate loop - projected h_time_psi_crate");
515 
516 
517  double default_mean_crate = h_time_crate_masked->GetMean();
518  double default_mean_crate_unc = h_time_crate_masked->GetMeanError();
519  double default_sigma_crate = h_time_crate_masked->GetStdDev();
520  B2INFO("Fitting crate between " << time_fit_min << " and " << time_fit_max);
521  TF1* gaus = new TF1("func", "gaus(0)", time_fit_min, time_fit_max);
522  gaus->SetParNames("numCrateHisNormalization", "mean", "sigma");
523  double hist_max = h_time_crate->GetMaximum();
524  double stddev = h_time_crate->GetStdDev();
525  sigma_crate = stddev;
526  mean_crate = default_mean_crate;
527  gaus->SetParameter(0, hist_max / 2.);
528  gaus->SetParameter(1, mean_crate);
529  gaus->SetParameter(2, sigma_crate);
530 
531  h_time_crate_masked->Fit(gaus, "LQR"); // L for likelihood, R for x-range, Q for fit quiet mode
532 
533  double fit_mean_crate = gaus->GetParameter(1);
534  double fit_mean_crate_unc = gaus->GetParError(1);
535  double fit_sigma_crate = gaus->GetParameter(2);
536 
537  double meanDiff = fit_mean_crate - default_mean_crate;
538  double meanUncDiff = fit_mean_crate_unc - default_mean_crate_unc;
539  double sigmaDiff = fit_sigma_crate - default_sigma_crate;
540 
541  bool good_fit = false;
542 
543  B2DEBUG(30, "Crate id = " << crate_id << " with crate mean = " << default_mean_crate << " +- " << fit_mean_crate_unc);
544 
545  if ((fabs(meanDiff) > 7) ||
546  (fabs(meanUncDiff) > 7) ||
547  (fabs(sigmaDiff) > 7) ||
548  (fit_mean_crate_unc > 3) ||
549  (fit_sigma_crate < 0.1) ||
550  (fit_mean_crate < time_fit_min) ||
551  (fit_mean_crate > time_fit_max)) {
552  B2INFO("Crate id = " << crate_id);
553  B2INFO("fit mean, default mean = " << fit_mean_crate << ", " << default_mean_crate);
554  B2INFO("fit sigma, default sigma = " << fit_sigma_crate << ", " << default_sigma_crate);
555 
556  B2INFO("crate fit mean - hist mean = " << meanDiff);
557  B2INFO("fit mean unc. - hist mean unc. = " << meanUncDiff);
558  B2INFO("fit sigma - hist sigma = " << sigmaDiff);
559  B2INFO("fit_mean_crate = " << fit_mean_crate);
560  B2INFO("time_fit_min = " << time_fit_min);
561  B2INFO("time_fit_max = " << time_fit_max);
562  } else {
563  good_fit = true;
564  }
565 
566  int numEntries = h_time_crate->GetEntries();
567  B2INFO("Number entries = " << numEntries);
568  if ((numEntries >= minNumEntries) && good_fit) {
569 
570  database_mean_crate = fit_mean_crate;
571  database_mean_crate_unc = fit_mean_crate_unc;
572 
573  } else {
574  database_mean_crate = default_mean_crate;
575  database_mean_crate_unc = default_mean_crate_unc;
576  }
577  if (numEntries == 0) {
578  database_mean_crate = 0;
579  database_mean_crate_unc = 0;
580  }
581 
582  tcrate_mean_new[crate_id - 1] = database_mean_crate;
583  tcrate_mean_unc_new[crate_id - 1] = database_mean_crate_unc;
584  tcrate_sigma_new[crate_id - 1] = fit_sigma_crate;
585  tcrate_new_was_set[crate_id - 1] = true;
586 
587 
588  histfile->WriteTObject(h_time_crate, (std::string("h_time_psi_crate") + std::to_string(crate_id)).c_str());
589  histfile->WriteTObject(h_time_crate_masked, (std::string("h_time_psi_crate_masked") + std::to_string(crate_id)).c_str());
590  histfile->WriteTObject(h_time_crate_rebin, (std::string("h_time_psi_crate_rebinned") + std::to_string(crate_id)).c_str());
591 
592  delete gaus;
593  }
594 
595  B2DEBUG(30, "crate histograms made");
596 
597 
598  // Save database for crates
599  // Vector of time offsets to be saved in the database.
600  vector<float> t_offsets_crate;
601  // Vector of time offset uncertainties to be saved in the database.
602  vector<float> t_offsets_crate_unc;
603  for (int i = 1; i <= 8736; i++) {
604  t_offsets_crate.push_back(0);
605  t_offsets_crate_unc.push_back(0);
606  }
607 
608 
609  for (int crys_id = 1; crys_id <= 8736; crys_id++) {
610 
611  int crate_id_from_crystal = crystalMapper->getCrateID(crys_id);
612  if (tcrate_new_was_set[crate_id_from_crystal - 1]) {
613  t_offsets_crate[crys_id - 1] = tcrate_mean_new[crate_id_from_crystal - 1] / TICKS_TO_NS;
614  t_offsets_crate_unc[crys_id - 1] = tcrate_mean_unc_new[crate_id_from_crystal - 1] / TICKS_TO_NS;
615 
616  } else {
617  t_offsets_crate[crys_id - 1] = tcrate_mean_prev[crate_id_from_crystal - 1];
618  B2INFO("used old crate mean but zeroed uncertainty since not saved in root file");
619  }
620 
621  }
622 
623  ECLCrystalCalib* BhabhaTCrateCalib = new ECLCrystalCalib();
624  BhabhaTCrateCalib->setCalibVector(t_offsets_crate, t_offsets_crate_unc);
625 
626 
627  // Save the information to the payload if there is at least one crate
628  // begin calibrated.
629  if (crateIDLo <= crateIDHi) {
630  saveCalibration(BhabhaTCrateCalib, "ECLCrateTimeOffset");
631  B2DEBUG(30, "crate payload made");
632  }
633 
634 
635  int tree_crateid;
636  int tree_runNum;
637  double tree_tcrate_mean;
638  double tree_tcrate_mean_unc;
639  double tree_tcrate_sigma;
640  double tree_tcrate_meanPrev;
641 
642  tree_crate->Branch("runNum", &tree_runNum)->SetTitle("Run number, 0..infinity and beyond!");
643  tree_crate->Branch("crateid", &tree_crateid)->SetTitle("Crate id, 1..52");
644  tree_crate->Branch("tcrate", &tree_tcrate_mean)->SetTitle("Crate time offset mean, tcrate, ns");
645  tree_crate->Branch("tcratePrev", &tree_tcrate_meanPrev)->SetTitle("Previous crate time offset mean, tcrate, ns");
646  tree_crate->Branch("tcrate_unc", &tree_tcrate_mean_unc)->SetTitle("Error of time tcrate mean, ns.");
647  tree_crate->Branch("tcrate_sigma", &tree_tcrate_sigma)->SetTitle("Sigma of time tcrate distribution, ns");
648  tree_crate->SetAutoSave(10);
649 
650 
651  for (auto expRun : getRunList()) {
652  // Key command to make sure your DBObjPtrs are correct
653  B2INFO("run num, exp num: " << expRun.second << ", " << expRun.first);
654  int runNumber = expRun.second;
655 
656  for (int crate_id = 1; crate_id <= 52; crate_id++) {
657  if (tcrate_new_was_set[crate_id - 1]) {
658  tree_runNum = runNumber;
659  tree_crateid = crate_id;
660  tree_tcrate_mean = tcrate_mean_new[crate_id - 1];
661  tree_tcrate_mean_unc = tcrate_mean_unc_new[crate_id - 1];
662  tree_tcrate_sigma = tcrate_sigma_new[crate_id - 1];
663  tree_tcrate_meanPrev = tcrate_mean_prev[crate_id - 1] * TICKS_TO_NS;
664  tree_crate->Fill();
665  }
666  }
667  }
668 
669  B2DEBUG(30, "end of crate corrections .....");
670 
671  tree_crystal->Write();
672  tree_crate->Write();
673 
674  histfile->Close();
675 
676  tree_crystal = 0;
677  tree_crate = 0;
678 
679  B2INFO("Finished talgorithm");
680  return c_OK;
681 }

◆ 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 21 of file CalibrationAlgorithm.cc.

◆ 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 63 of file CalibrationAlgorithm.cc.

◆ execute()

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.

◆ 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 algorihtm.

Definition at line 174 of file CalibrationAlgorithm.h.

◆ 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 159 of file CalibrationAlgorithm.cc.

◆ 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 187 of file CalibrationAlgorithm.cc.

Member Data Documentation

◆ meanCleanCutMinFactor

double meanCleanCutMinFactor

After rebinning, create a mask for bins that have values less than meanCleanCutMinFactor times the maximum bin value.

Expand mask and apply to non-rebinned histogram.

Definition at line 53 of file eclBhabhaTAlgorithm.h.


The documentation for this class was generated from the following files:
Belle2::ECL::eclBhabhaTAlgorithm::cellIDHi
int cellIDHi
Fit crystals with cellID0 in the inclusive range [cellIDLo,cellIDHi].
Definition: eclBhabhaTAlgorithm.h:51
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition: CalibrationAlgorithm.cc:290
Belle2::ECL::eclBhabhaTAlgorithm::debugOutput
bool debugOutput
Save every histogram and fitted function to debugFilename.
Definition: eclBhabhaTAlgorithm.h:58
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::ECLCrystalCalib
General DB object to store one calibration number per ECL crystal.
Definition: ECLCrystalCalib.h:34
Belle2::ECL::eclBhabhaTAlgorithm::crateIDLo
int crateIDLo
Fit crates with crateID0 in the inclusive range [crateIDLo,crateIDHi].
Definition: eclBhabhaTAlgorithm.h:56
Belle2::CalibrationAlgorithm::updateDBObjPtrs
void updateDBObjPtrs(const unsigned int event, const int run, const int experiment)
Updates any DBObjPtrs by calling update(event) for DBStore.
Definition: CalibrationAlgorithm.cc:398
Belle2::ECL::eclBhabhaTAlgorithm::crateIDHi
int crateIDHi
Fit crates with crateID0 in the inclusive range [crateIDLo,crateIDHi].
Definition: eclBhabhaTAlgorithm.h:57
Belle2::ECL::ECLChannelMapper
This class provides access to ECL channel map that is either a) Loaded from the database (see ecl/dbo...
Definition: ECLChannelMapper.h:36
Belle2::CalibrationAlgorithm::c_Failure
@ c_Failure
Failed =3 in Python.
Definition: CalibrationAlgorithm.h:54
Belle2::ECL::EclConfiguration::m_rf
static constexpr double m_rf
accelerating RF, http://ptep.oxfordjournals.org/content/2013/3/03A006.full.pdf
Definition: EclConfiguration.h:45
Belle2::CalibrationAlgorithm::getRunList
const std::vector< Calibration::ExpRun > & getRunList() const
Get the list of runs for which calibration is called.
Definition: CalibrationAlgorithm.h:276
Belle2::ECL::eclBhabhaTAlgorithm::cellIDLo
int cellIDLo
Fit crystals with cellID0 in the inclusive range [cellIDLo,cellIDHi].
Definition: eclBhabhaTAlgorithm.h:50
Belle2::ECL::eclBhabhaTAlgorithm::meanCleanRebinFactor
double meanCleanRebinFactor
Rebinning factor for mean calculation.
Definition: eclBhabhaTAlgorithm.h:52
Belle2::ECLCrystalCalib::setCalibVector
void setCalibVector(const std::vector< float > &CalibConst, const std::vector< float > &CalibConstUnc)
Set vector of constants with uncertainties.
Definition: ECLCrystalCalib.h:48
Belle2::ECL::eclBhabhaTAlgorithm::debugFilenameBase
std::string debugFilenameBase
Name of file with debug output, eclBhabhaTAlgorithm.root by default.
Definition: eclBhabhaTAlgorithm.h:60
Belle2::ECL::eclBhabhaTAlgorithm::meanCleanCutMinFactor
double meanCleanCutMinFactor
After rebinning, create a mask for bins that have values less than meanCleanCutMinFactor times the ma...
Definition: eclBhabhaTAlgorithm.h:53