Belle II Software  release-05-01-25
TOPCommonT0BFAlgorithm.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * *
6  * Author: The Belle II Collaboration *
7  * Contributors: Marko Staric *
8  * *
9  * This software is provided "as is" without any warranty. *
10  **************************************************************************/
11 
12 #include <top/calibration/TOPCommonT0BFAlgorithm.h>
13 #include <top/dbobjects/TOPCalCommonT0.h>
14 #include <math.h>
15 #include <TROOT.h>
16 #include <TF1.h>
17 #include <TFile.h>
18 #include <TTree.h>
19 #include <string>
20 
21 using namespace std;
22 
23 namespace Belle2 {
28  namespace TOP {
29 
30  TOPCommonT0BFAlgorithm::TOPCommonT0BFAlgorithm():
31  CalibrationAlgorithm("TOPCommonT0BFCollector")
32  {
33  setDescription("Calibration algorithm for common T0 calibration "
34  "with a fit of bunch finder residuals (method BF)");
35 
36  }
37 
39  {
40  gROOT->SetBatch();
41 
42  // get histogram to fit and check if statistics is sufficient
43 
44  auto h = getHistogram();
45  if (not h) {
46  B2ERROR("TOPCommonT0BFAlgorithm: no histogram to fit");
47  return c_NotEnoughData;
48  }
49 
50  int numEntries = h->GetSumOfWeights();
51  if (numEntries < m_minEntries) {
52  B2ERROR("TOPCommonT0BFAlgorithm: too few entries to fit the histogram");
53  return c_NotEnoughData;
54  }
55 
56  m_bunchTimeSep = h->GetXaxis()->GetXmax() - h->GetXaxis()->GetXmin();
57 
58  // construct file name, open output root file and book output tree
59 
60  const auto& expRun = getRunList();
61  string expNo = to_string(expRun[0].first);
62  while (expNo.length() < 4) expNo.insert(0, "0");
63  string runNo = to_string(expRun[0].second);
64  while (runNo.length() < 5) runNo.insert(0, "0");
65  string outputFileName = "commonT0-e" + expNo + "-r" + runNo + ".root";
66  auto* file = TFile::Open(outputFileName.c_str(), "recreate");
67 
68  auto* tree = new TTree("tree", "common T0 calibration results");
69  tree->Branch<int>("expNum", &m_expNo);
70  tree->Branch<int>("runNum", &m_runNo);
71  tree->Branch<int>("runLast", &m_runLast);
72  tree->Branch<float>("fitted_offset", &m_fittedOffset);
73  tree->Branch<float>("offset", &m_offset);
74  tree->Branch<float>("offsetErr", &m_offsetError);
75  tree->Branch<float>("sigma", &m_sigmaCore);
76  tree->Branch<float>("sigmaTail", &m_sigmaTail);
77  tree->Branch<float>("tailFract", &m_tailFract);
78  tree->Branch<float>("chi2", &m_chi2);
79  tree->Branch<float>("integral", &m_integral);
80  tree->Branch<int>("nEvt", &m_numEvents);
81  tree->Branch<int>("fitStatus", &m_fitStatus);
82 
83  m_expNo = expRun[0].first;
84  m_runNo = expRun[0].second;
85  m_runLast = expRun.back().second;
86 
87  // fit histogram
88 
89  if (numEntries > m_cutoffEntries and m_tailFractInit > 0) {
90  int status = fitDoubleGaus(h);
91  if (status != 0) fitSingleGaus(h);
92  } else {
93  fitSingleGaus(h);
94  }
95  tree->Fill();
96 
97  // write the results and close the file
98 
99  h->Write();
100  file->Write();
101  file->Close();
102 
103  // check the results and return if not good enough
104 
105  if (m_fitStatus != 0 or m_offsetError > m_minError) return c_NotEnoughData;
106 
107  // otherwise create and import payload to DB
108 
109  auto* commonT0 = new TOPCalCommonT0(m_offset, m_offsetError);
110  saveCalibration(commonT0);
111 
112  return c_OK;
113  }
114 
115  std::shared_ptr<TH1F> TOPCommonT0BFAlgorithm::getHistogram()
116  {
117  auto h1a = getObjectPtr<TH1F>("offset_a");
118  auto h1b = getObjectPtr<TH1F>("offset_b");
119  if (not h1a) {
120  B2ERROR("TOPCommonT0BFAlgorithm: histogram 'offset_a' not found");
121  return 0;
122  }
123  if (not h1b) {
124  B2ERROR("TOPCommonT0BFAlgorithm: histogram 'offset_b' not found");
125  return 0;
126  }
127 
128  int halfbins = h1a->GetNbinsX() / 2;
129  if (abs(h1a->GetMaximumBin() - halfbins) < abs(h1b->GetMaximumBin() - halfbins)) {
130  return h1a;
131  } else {
132  return h1b;
133  }
134  }
135 
136  int TOPCommonT0BFAlgorithm::fitSingleGaus(std::shared_ptr<TH1F> h)
137  {
138  double sum = h->GetSumOfWeights();
139  if (sum < 5) return 5;
140  double maxPosition = h->GetBinCenter(h->GetMaximumBin());
141  double binWidth = h->GetBinWidth(1);
142  double xmin = h->GetXaxis()->GetXmin();
143  double xmax = h->GetXaxis()->GetXmax();
144 
145  auto* func = new TF1("func1g",
146  "[0] + [1]/sqrt(2*pi)/[3]*exp(-0.5*((x-[2])/[3])**2)",
147  xmin, xmax);
148  func->SetParameter(0, 0);
149  func->SetParameter(1, sum * binWidth);
150  func->SetParameter(2, maxPosition);
151  func->SetParameter(3, m_sigmaCoreInit);
152 
153  int status = h->Fit(func, "LRSQ");
154 
155  m_fittedOffset = func->GetParameter(2);
157  m_offsetError = func->GetParError(2);
158  m_sigmaCore = func->GetParameter(3);
159  m_sigmaTail = 0;
160  m_tailFract = 0;
161  m_chi2 = func->GetChisquare() / func->GetNDF();
162  m_integral = func->GetParameter(1) / binWidth;
163  m_numEvents = sum;
164  m_fitStatus = status;
165 
166  return status;
167  }
168 
169 
170  int TOPCommonT0BFAlgorithm::fitDoubleGaus(std::shared_ptr<TH1F> h)
171  {
172  double sum = h->GetSumOfWeights();
173  if (sum < 7) return 7;
174  double maxPosition = h->GetBinCenter(h->GetMaximumBin());
175  double binWidth = h->GetBinWidth(1);
176  double xmin = h->GetXaxis()->GetXmin();
177  double xmax = h->GetXaxis()->GetXmax();
178 
179  auto* func = new TF1("func2g",
180  "[0] + [1]*((1-[4])/sqrt(2*pi)/[3]*exp(-0.5*((x-[2])/[3])**2)"
181  "+ [4]/sqrt(2*pi)/[5]*exp(-0.5*((x-[2])/[5])**2))",
182  xmin, xmax);
183  func->SetParameter(0, 0);
184  func->SetParameter(1, sum * binWidth);
185  func->SetParameter(2, maxPosition);
186  func->SetParameter(3, m_sigmaCoreInit);
187  func->SetParameter(4, m_tailFractInit);
188  func->SetParameter(5, m_sigmaTailInit);
189 
190  int status = h->Fit(func, "LRSQ");
191 
192  m_fittedOffset = func->GetParameter(2);
194  m_offsetError = func->GetParError(2);
195  m_sigmaCore = func->GetParameter(3);
196  m_sigmaTail = func->GetParameter(5);
197  m_tailFract = func->GetParameter(4);
198  m_chi2 = func->GetChisquare() / func->GetNDF();
199  m_integral = func->GetParameter(1) / binWidth;
200  m_numEvents = sum;
201  m_fitStatus = status;
202 
203  return status;
204  }
205 
206 
207  } // end namespace TOP
209 } // end namespace Belle2
Belle2::TOP::TOPCommonT0BFAlgorithm::fitDoubleGaus
int fitDoubleGaus(std::shared_ptr< TH1F > h)
Fit double gaus w/ same mean + constant.
Definition: TOPCommonT0BFAlgorithm.cc:170
Belle2::TOP::TOPCommonT0BFAlgorithm::m_tailFract
float m_tailFract
tail fraction (set to 0 for single gauss fit)
Definition: TOPCommonT0BFAlgorithm.h:119
Belle2::TOP::TOPCommonT0BFAlgorithm::m_offsetError
float m_offsetError
error on fitted offset
Definition: TOPCommonT0BFAlgorithm.h:116
Belle2::TOP::TOPCommonT0BFAlgorithm::m_chi2
float m_chi2
normalized chi^2
Definition: TOPCommonT0BFAlgorithm.h:120
Belle2::TOP::TOPCommonT0BFAlgorithm::calibrate
virtual EResult calibrate() final
algorithm implementation
Definition: TOPCommonT0BFAlgorithm.cc:38
Belle2::TOP::TOPCommonT0BFAlgorithm::m_numEvents
int m_numEvents
number of all events in the histogram
Definition: TOPCommonT0BFAlgorithm.h:122
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::TOP::TOPCommonT0BFAlgorithm::getHistogram
std::shared_ptr< TH1F > getHistogram()
Returns histogram to be fitted.
Definition: TOPCommonT0BFAlgorithm.cc:115
Belle2::TOP::TOPCommonT0BFAlgorithm::m_sigmaCore
float m_sigmaCore
core gaussian sigma
Definition: TOPCommonT0BFAlgorithm.h:117
Belle2::TOP::TOPCommonT0BFAlgorithm::m_fitStatus
int m_fitStatus
fit status
Definition: TOPCommonT0BFAlgorithm.h:123
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:331
Belle2::TOP::TOPCommonT0BFAlgorithm::m_sigmaTail
float m_sigmaTail
tail gaussian sigma (set to 0 for single gauss fit)
Definition: TOPCommonT0BFAlgorithm.h:118
Belle2::TOP::TOPCommonT0BFAlgorithm::m_fittedOffset
float m_fittedOffset
fitted offset
Definition: TOPCommonT0BFAlgorithm.h:114
Belle2::TOP::TOPCommonT0BFAlgorithm::m_expNo
int m_expNo
experiment number
Definition: TOPCommonT0BFAlgorithm.h:111
Belle2::TOP::TOPCommonT0BFAlgorithm::m_runLast
int m_runLast
last run number
Definition: TOPCommonT0BFAlgorithm.h:113
Belle2::TOP::TOPCommonT0BFAlgorithm::m_minError
double m_minError
minimal commonT0 uncertainty [ns] to declare c_OK
Definition: TOPCommonT0BFAlgorithm.h:108
Belle2::TOP::TOPCommonT0BFAlgorithm::m_integral
float m_integral
number of fitted signal events
Definition: TOPCommonT0BFAlgorithm.h:121
Belle2::TOP::TOPCommonT0BFAlgorithm::m_runNo
int m_runNo
first run number
Definition: TOPCommonT0BFAlgorithm.h:112
Belle2::TOP::TOPCommonT0BFAlgorithm::m_offset
float m_offset
wrap-around of fitted offset (= common T0)
Definition: TOPCommonT0BFAlgorithm.h:115
Belle2::TOP::TOPCommonT0BFAlgorithm::m_cutoffEntries
int m_cutoffEntries
cutoff entries for single/double gaussian fit
Definition: TOPCommonT0BFAlgorithm.h:107
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TOP::TOPCommonT0BFAlgorithm::m_sigmaTailInit
double m_sigmaTailInit
tail gaussian sigma [ns]
Definition: TOPCommonT0BFAlgorithm.h:105
Belle2::TOP::TOPCommonT0BFAlgorithm::m_tailFractInit
double m_tailFractInit
fraction of tail gaussian
Definition: TOPCommonT0BFAlgorithm.h:106
Belle2::TOP::TOPCommonT0BFAlgorithm::m_minEntries
int m_minEntries
minimal number of histogram entries to perform fit
Definition: TOPCommonT0BFAlgorithm.h:103
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
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::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:53
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47
Belle2::TOPCalCommonT0
Common T0 calibration constant.
Definition: TOPCalCommonT0.h:33
Belle2::TOP::TOPCommonT0BFAlgorithm::fitSingleGaus
int fitSingleGaus(std::shared_ptr< TH1F > h)
Fit single gaus + constant.
Definition: TOPCommonT0BFAlgorithm.cc:136
Belle2::TOP::TOPCommonT0BFAlgorithm::m_bunchTimeSep
double m_bunchTimeSep
bunch separation in time [ns]
Definition: TOPCommonT0BFAlgorithm.h:126
Belle2::TOP::TOPCommonT0BFAlgorithm::m_sigmaCoreInit
double m_sigmaCoreInit
core gaussian sigma [ns]
Definition: TOPCommonT0BFAlgorithm.h:104