Belle II Software  release-08-01-10
SVDCoGTimeCalibrationAlgorithm.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 #include <svd/calibration/SVDCoGTimeCalibrationAlgorithm.h>
9 
10 #include <svd/dbobjects/SVDCoGCalibrationFunction.h>
11 #include <svd/calibration/SVDCoGTimeCalibrations.h>
12 
13 #include <TF1.h>
14 #include <TProfile.h>
15 #include <TH1F.h>
16 #include <TH2F.h>
17 #include <TH3F.h>
18 #include <framework/logging/Logger.h>
19 #include <iostream>
20 #include <TString.h>
21 #include <TFitResult.h>
22 
23 using namespace std;
24 using namespace Belle2;
25 
26 SVDCoGTimeCalibrationAlgorithm::SVDCoGTimeCalibrationAlgorithm(const std::string& str) :
27  CalibrationAlgorithm("SVDTimeCalibrationCollector")
28 {
29  setDescription("SVDCoGTimeCalibration calibration algorithm");
30  m_id = str;
31 }
32 
34 {
35 
36  gROOT->SetBatch(true);
37 
38  auto timeCal = new Belle2::SVDCoGCalibrationFunction();
39  auto payload = new Belle2::SVDCoGTimeCalibrations::t_payload(*timeCal, m_id);
40 
41  std::unique_ptr<TF1> pol1(new TF1("pol1", "[0] + [1]*x", -10, 80));
42  pol1->SetParameters(-40, 0.9);
43  std::unique_ptr<TF1> pol3(new TF1("pol3", "[0] + [1]*x + [2]*x*x + [3]*x*x*x", -10, 80));
44  pol3->SetParLimits(0, -200, 0);
45  pol3->SetParLimits(1, 0, 10);
46  pol3->SetParLimits(2, -1, 0);
47  pol3->SetParLimits(3, 0, 1);
48  std::unique_ptr<TF1> pol5(new TF1("pol5", "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + [5]*x*x*x*x*x", -100, 100));
49  pol5->SetParameters(-50, 1.5, 0.01, 0.0001, 0.00001, 0.000001);
50 
51  FileStat_t info;
52  int cal_rev = 1;
53  while (gSystem->GetPathInfo(Form("algorithm_6SampleCoG_output_rev_%d.root", cal_rev), info) == 0)
54  cal_rev++;
55  std::unique_ptr<TFile> f(new TFile(Form("algorithm_6SampleCoG_output_rev_%d.root", cal_rev), "RECREATE"));
56 
57  auto m_tree = new TTree(Form("rev_%d", cal_rev), "RECREATE");
58  int layer_num, ladder_num, sensor_num, view, ndf;
59  float a, b, c, d, a_err, b_err, c_err, d_err, chi2, p;
60  m_tree->Branch("layer", &layer_num, "layer/I");
61  m_tree->Branch("ladder", &ladder_num, "ladder/I");
62  m_tree->Branch("sensor", &sensor_num, "sensor/I");
63  m_tree->Branch("isU", &view, "isU/I");
64  m_tree->Branch("a", &a, "a/F");
65  m_tree->Branch("b", &b, "b/F");
66  m_tree->Branch("c", &c, "c/F");
67  m_tree->Branch("d", &d, "d/F");
68  m_tree->Branch("a_err", &a_err, "a_err/F");
69  m_tree->Branch("b_err", &b_err, "b_err/F");
70  m_tree->Branch("c_err", &c_err, "c_err/F");
71  m_tree->Branch("d_err", &d_err, "d_err/F");
72  m_tree->Branch("chi2", &chi2, "chi2/F");
73  m_tree->Branch("ndf", &ndf, "ndf/I");
74  m_tree->Branch("p", &p, "p/F");
75 
76  auto __hEventT0vsCoG__ = getObjectPtr<TH3F>("__hEventT0vsCoG__");
77  auto __hEventT0__ = getObjectPtr<TH2F>("__hEventT0__");
78  auto __hEventT0NoSync__ = getObjectPtr<TH2F>("__hEventT0NoSync__");
79  auto __hBinToSensorMap__ = getObjectPtr<TH1F>("__hBinToSensorMap__");
80 
81  for (int ij = 0; ij < (__hBinToSensorMap__->GetNbinsX()); ij++) {
82  {
83  {
84  {
85 
86  auto binLabel = __hBinToSensorMap__->GetXaxis()->GetBinLabel(ij + 1);
87  char side;
88  std::sscanf(binLabel, "L%dL%dS%d%c", &layer_num, &ladder_num, &sensor_num, &side);
89  view = 0;
90  if (side == 'U')
91  view = 1;
92 
93  B2INFO("Projecting for Sensor: " << binLabel << " with Bin Number: " << ij + 1);
94 
95  __hEventT0vsCoG__->GetZaxis()->SetRange(ij + 1, ij + 1);
96  auto hEventT0vsCoG = (TH2D*)__hEventT0vsCoG__->Project3D("yxe");
97  auto hEventT0 = (TH1D*)__hEventT0__->ProjectionX("hEventT0_tmp", ij + 1, ij + 1);
98  auto hEventT0nosync = (TH1D*)__hEventT0NoSync__->ProjectionX("hEventT0NoSync_tmp", ij + 1, ij + 1);
99 
100  hEventT0vsCoG->SetName(Form("eventT0vsCoG__L%dL%dS%d%c", layer_num, ladder_num, sensor_num, side));
101  hEventT0->SetName(Form("eventT0__L%dL%dS%d%c", layer_num, ladder_num, sensor_num, side));
102  hEventT0nosync->SetName(Form("eventT0nosync__L%dL%dS%d%c", layer_num, ladder_num, sensor_num, side));
103 
104  char sidePN = (side == 'U' ? 'P' : 'N');
105  hEventT0vsCoG->SetTitle(Form("EventT0Sync vs rawTime in %d.%d.%d %c/%c", layer_num, ladder_num, sensor_num, side, sidePN));
106  hEventT0->SetTitle(Form("EventT0Sync in %d.%d.%d %c/%c", layer_num, ladder_num, sensor_num, side, sidePN));
107  hEventT0nosync->SetTitle(Form("EventT0NoSync in %d.%d.%d %c/%c", layer_num, ladder_num, sensor_num, side, sidePN));
108 
109  hEventT0vsCoG->SetDirectory(0);
110  hEventT0->SetDirectory(0);
111  hEventT0nosync->SetDirectory(0);
112 
113  B2INFO("Histogram: " << hEventT0vsCoG->GetName() <<
114  " Entries (n. clusters): " << hEventT0vsCoG->GetEntries());
115  if (layer_num == 3 && hEventT0vsCoG->GetEntries() < m_minEntries) {
116  B2INFO("Histogram: " << hEventT0vsCoG->GetName() <<
117  " Entries (n. clusters): " << hEventT0vsCoG->GetEntries() <<
118  " Entries required: " << m_minEntries);
119  B2WARNING("Not enough data, adding one run to the collector");
120  f->Close();
121  gSystem->Unlink(Form("algorithm_6SampleCoG_output_rev_%d.root", cal_rev));
122  return c_NotEnoughData;
123  }
124  if (layer_num != 3 && hEventT0vsCoG->GetEntries() < m_minEntries / 10) {
125  B2INFO("Histogram: " << hEventT0vsCoG->GetName() <<
126  " Entries (n. clusters): " << hEventT0vsCoG->GetEntries() <<
127  " Entries required: " << m_minEntries / 10);
128  B2WARNING("Not enough data, adding one run to the collector");
129  f->Close();
130  gSystem->Unlink(Form("algorithm_6SampleCoG_output_rev_%d.root", cal_rev));
131  return c_NotEnoughData;
132  }
133  for (int i = 1; i <= hEventT0vsCoG->GetNbinsX(); i++) {
134  for (int j = 1; j <= hEventT0vsCoG->GetNbinsY(); j++) {
135  if (hEventT0vsCoG->GetBinContent(i, j) < max(2, int(hEventT0vsCoG->GetEntries() * 0.001))) {
136  hEventT0vsCoG->SetBinContent(i, j, 0);
137  }
138  }
139  }
140  TProfile* pfx = hEventT0vsCoG->ProfileX();
141  std::string name = "pfx_" + std::string(hEventT0vsCoG->GetName());
142  pfx->SetName(name.c_str());
143  TFitResultPtr tfr = pfx->Fit("pol3", "RQSM");
144  double par[4];
145  pol3->GetParameters(par);
147  /*
148  pfx->Fit("pol1", "RQ");
149  double par[4];
150  pol1->GetParameters(par);
151  par[2] = 0;
152  par[3] = 0;
153  */
154  // double meanT0 = hEventT0->GetMean();
155  // double meanT0NoSync = hEventT0nosync->GetMean();
156  timeCal->set_current(1);
157  // timeCal->set_current(2);
158  timeCal->set_pol3parameters(par[0], par[1], par[2], par[3]);
159  payload->set(layer_num, ladder_num, sensor_num, bool(view), 1, *timeCal);
160  f->cd();
161  hEventT0->Write();
162  hEventT0vsCoG->Write();
163  hEventT0nosync->Write();
164  pfx->Write();
165 
166  delete pfx;
167  delete hEventT0vsCoG;
168  delete hEventT0;
169  delete hEventT0nosync;
170 
171  if (tfr.Get() == nullptr || (tfr->Status() != 0 && tfr->Status() != 4 && tfr->Status() != 4000)) {
172  f->Close();
173  B2FATAL("Fit to the histogram failed in SVDCoGTimeCalibrationAlgorithm. "
174  << "Check the 2-D histogram to clarify the reason.");
175  } else {
176  a = par[0]; b = par[1]; c = par[2]; d = par[3];
177  a_err = tfr->ParError(0); b_err = tfr->ParError(1); c_err = tfr->ParError(2); d_err = tfr->ParError(3);
178  chi2 = tfr->Chi2();
179  ndf = tfr->Ndf();
180  p = tfr->Prob();
181  m_tree->Fill();
182  }
183  }
184  }
185  }
186  }
187  m_tree->Write();
188  f->Close();
189  saveCalibration(payload, "SVDCoGTimeCalibrations");
190 
191  //delete f;
192 
193  // probably not needed - would trigger re-doing the collection
194  // if ( ... too large corrections ... ) return c_Iterate;
195  return c_OK;
196 }
197 
198 bool SVDCoGTimeCalibrationAlgorithm::isBoundaryRequired(const Calibration::ExpRun& currentRun)
199 {
200  float meanRawTimeL3V = 0;
201  // auto eventT0Hist = getObjectPtr<TH1F>("hEventT0FromCDC");
202  auto rawTimeL3V = getObjectPtr<TH1F>("hRawTimeL3V");
203  // float meanEventT0 = eventT0Hist->GetMean();
204  if (!rawTimeL3V) {
206  meanRawTimeL3V = m_previousRawTimeMeanL3V.value();
207  } else {
208  if (rawTimeL3V->GetEntries() > m_minEntries)
209  meanRawTimeL3V = rawTimeL3V->GetMean();
210  else {
212  meanRawTimeL3V = m_previousRawTimeMeanL3V.value();
213  }
214  }
216  B2INFO("Setting start payload boundary to be the first run ("
217  << currentRun.first << "," << currentRun.second << ")");
218  m_previousRawTimeMeanL3V.emplace(meanRawTimeL3V);
219 
220  return true;
221  } else if (abs(meanRawTimeL3V - m_previousRawTimeMeanL3V.value()) > m_allowedTimeShift) {
222  B2INFO("Histogram mean has shifted from " << m_previousRawTimeMeanL3V.value()
223  << " to " << meanRawTimeL3V << ". We are requesting a new payload boundary for ("
224  << currentRun.first << "," << currentRun.second << ")");
225  m_previousRawTimeMeanL3V.emplace(meanRawTimeL3V);
226  return true;
227  } else {
228  return false;
229  }
230 }
Base class for calibration algorithms.
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
EResult
The result of calibration.
@ c_OK
Finished successfuly =0 in Python.
@ c_NotEnoughData
Needs more data =2 in Python.
class to contain the CoG Time calibrations
std::string m_id
Parameter given to set the UniqueID of the payload.
std::optional< float > m_previousRawTimeMeanL3V
CoG time mean of the previous run for V side of layer 3.
float m_minEntries
Set the minimun number of entries required in the histograms of layer 3.
virtual EResult calibrate() override
Run algo on data.
virtual bool isBoundaryRequired(const Calibration::ExpRun &currentRun) override
If the event T0 changes significantly return true.
SVDCalibrationsBase< SVDCalibrationsScalar< SVDCoGCalibrationFunction > > t_payload
typedef for the SVDCoGCalibrationFunction payload of all SVD sensors
Abstract base class for different kinds of events.