Belle II Software development
HitBasedT0Extractor.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 <tracking/eventTimeExtraction/findlets/HitBasedT0Extractor.h>
9#include <tracking/eventTimeExtraction/findlets/BaseEventTimeExtractor.icc.h>
10
11#include <tracking/trackFindingCDC/utilities/StringManipulation.h>
12
13#include <framework/core/ModuleParamList.h>
14
15#include <framework/logging/Logger.h>
16#include <framework/utilities/ScopeGuard.h>
17
18#include <TH1D.h>
19#include <TF1.h>
20#include <TFitResult.h>
21#include <TCanvas.h>
22
23#include <boost/lexical_cast.hpp>
24
25#include <memory>
26
27using namespace Belle2;
28using namespace TrackFindingCDC;
29
31{
32 return "Extracts the T0 time of an event only using CDC Hits";
33}
34
40
42 const std::string& prefix)
43{
44 moduleParamList->addParameter(prefixed(prefix, "searchWindow"),
46 "Size of the time distance (in ns) to search for t0 in both direction from the current best t0",
48
49 moduleParamList->addParameter(prefixed(prefix, "fitWindow"),
51 "Size of the time distance (in ns) to used for fitting",
53
54 moduleParamList->addParameter(prefixed(prefix, "refitWindow"),
56 "Size of the time distance (in ns) to used for re-fitting",
58
59 moduleParamList->addParameter(prefixed(prefix, "binCountTimeHistogram"),
61 "Number of bins in the timing histogram used for fitting",
63
64 moduleParamList->addParameter(prefixed(prefix, "rejectByBackgroundFlag"),
66 "Don't consider hits if they have the background flag set",
68
69 moduleParamList->addParameter(prefixed(prefix, "rejectIfNotTakenFlag"),
71 "Don't consider hits which have not been assigned during track finding. The CDC track finding has "
72 "to be run before for this flag to be useful.",
74
75 moduleParamList->addParameter(prefixed(prefix, "rejectIfChiSquareLargerThan"),
77 "consider the t0 fit failed which have larger chi2 than this number",
79
80 moduleParamList->addParameter(prefixed(prefix, "rejectIfUncertaintyLargerThan"),
82 "consider the t0 fit if the uncertainty on t0 is larger than this value",
84
85 moduleParamList->addParameter(prefixed(prefix, "storeAllFits"),
87 "store images for all fits",
89
90 moduleParamList->addParameter(prefixed(prefix, "minHitCount"),
92 "Minimum amount of hits which is required to try the extraction",
94
95 Super::exposeParameters(moduleParamList, prefix);
96}
97
98void HitBasedT0Extractor::apply(std::vector<CDCWireHit>& inputWireHits)
99{
100 const auto timeHistogramName = "HitBasedT0Extractor_time_hist";
101 const std::string debugImageName = "HitBasedT0Extractor_debug_" + boost::lexical_cast<std::string>
102 (m_eventMetaData->getEvent()) + ".png";
103
104 auto timingHistgram = TH1D(timeHistogramName, timeHistogramName,
107
108 // Enable batch mode - we do not want to show the canvas etc.
109 auto batchGuard = ScopeGuard::guardBatchMode();
110
111 TCanvas canvas(debugImageName.c_str(), debugImageName.c_str(), 800, 600);
112
113 if (inputWireHits.size() == 0) {
114 B2WARNING("No input CDC hits available for the CDC hit based t0 extraction.");
115 return;
116 }
117
118 for (auto const& wireHit : inputWireHits) {
119
121 && wireHit.getAutomatonCell().hasBackgroundFlag())
122 continue;
123
125 && !wireHit.getAutomatonCell().hasTakenFlag())
126 continue;
127
128 // the taken flag is also set for background hits. Therefore we must
129 // also hits to be also classified as background to not get flooded
130 // with background hits
132 && wireHit.getAutomatonCell().hasTakenFlag()
133 && wireHit.getAutomatonCell().hasBackgroundFlag())
134 continue;
135
136 // Attention: at this stage we use the drift time that was set when the WireHits were created.
137 // We will just *assume* that it was zero at this time!
138 timingHistgram.Fill(wireHit.getDriftTime());
139 }
140
141 if (timingHistgram.GetEntries() < m_param_minHitCount) {
142 B2DEBUG(25, "Only " << timingHistgram.GetEntries() << " hits satisfied the requirements for t0 extraction, " << m_param_minHitCount
143 << " are required.");
144 return;
145 }
146
147 // add an overall offset of 1 to not have to care about empty bins in
148 // the cumulated plot
149 timingHistgram.SetBinContent(1, timingHistgram.GetBinContent(1) + 1);
150
151 B2DEBUG(25,
152 "Filled histogram with " << timingHistgram.GetEntries() << " Entries");
153
154 std::unique_ptr<TH1> cumTimingHistogram(timingHistgram.GetCumulative());
155 // detach histogram from any directory so we can delete it ourselves
156 cumTimingHistogram->SetDirectory(0);
157
158 // set the error for each bin taking into account that the bin content is the
159 // cumulated content of previous bins
160 double errPrevBin = 0.0f;
161 double errThisBin = 0.0f;
162 for (int i = 0; i < cumTimingHistogram->GetNbinsX(); i++) {
163 if (i > 0) {
164 const double prevEntries = cumTimingHistogram->GetBinContent(i - 1);
165 const double addedEntries = cumTimingHistogram->GetBinContent(i) - prevEntries;
166 if (addedEntries > 0.0) {
167 // combine the error of the previous bin with the new entries of this bin
168 const double errNewEntries = 1.0 / std::sqrt(addedEntries);
169 errThisBin = errPrevBin + errNewEntries;
170 }
171
172 } else {
173 if (cumTimingHistogram->GetBinContent(i) > 0) {
174 errThisBin = 1.0 / std::sqrt(cumTimingHistogram->GetBinContent(i));
175 }
176 }
177
178 cumTimingHistogram->SetBinError(i, errThisBin);
179 errPrevBin = errThisBin;
180 }
181
183 cumTimingHistogram->Draw();
184 }
185
186 auto rangeBkg = std::make_pair(-m_param_fitWindow, -m_param_searchWindow);
187 auto rangeSig = std::make_pair(m_param_searchWindow, m_param_fitWindow);
188 // fit the background and signal side of the time distribution
189 TF1 fitfuncBkg = TF1("HitBasedT0Extractor_fit_bkg", "[0]*x + [1]",
190 rangeBkg.first, rangeBkg.second);
191 TF1 fitfuncSig = TF1("HitBasedT0Extractor_fit_sig", "[0]*x + [1]",
192 rangeSig.first, rangeSig.second);
193
194 auto fitresBkg = cumTimingHistogram->Fit(&fitfuncBkg, "LQS", "",
195 rangeBkg.first, rangeBkg.second);
196 auto fitresSig = cumTimingHistogram->Fit(&fitfuncSig, "LQS", "",
197 rangeSig.first, rangeSig.second);
198
200 fitfuncBkg.Draw("SAME");
201 fitfuncSig.Draw("SAME");
202 }
203
204 TF1 fitfuncSegmented = TF1("HitBasedT0Extractor_fit_seg",
205 "[1]*(x+TMath::Abs(x+[3])) + [2]*(x-TMath::Abs(x+[3])) + [0]",
207
208 if (fitresSig->IsValid() && fitresBkg->IsValid()) {
209 double t0_estimate = (fitresBkg->Parameter(1) - fitresSig->Parameter(1))
210 / (fitresSig->Parameter(0) - fitresBkg->Parameter(0));
211
212 // apply segmented fit
213 std::array<double, 4> fit_params = { 0, // = overall background offset
214 fitresSig->Parameter(0), // signal hits slope
215 fitresBkg->Parameter(0), // background hits slope
216 0.0 // breaking point shift
217 };
218
219 // use t0 estimate if it is something useful
220 if (std::abs(t0_estimate) < m_param_searchWindow) {
221 fit_params[3] = t0_estimate;
222 }
223
224 fitfuncSegmented.SetParameters(fit_params.data());
225
226 auto fitresFull = cumTimingHistogram->Fit(&fitfuncSegmented, "QS", "", -m_param_fitWindow, m_param_fitWindow);
228 fitfuncSegmented.Draw("SAME");
229 }
230
231 // refit with a fixed window around the extracted t0 to remove a possible bias
232 // because we fitted a large part of the signal side in the first iteration
233 const double fitted_t0_first = -fitresFull->Parameter(3);
234 bool refitSuccess = false;
235 if (std::abs(fitted_t0_first) < m_param_searchWindow) {
236 fitresFull = cumTimingHistogram->Fit(&fitfuncSegmented, "QS", "",
237 fitted_t0_first - m_param_refitWindow, fitted_t0_first + m_param_refitWindow);
238 refitSuccess = true;
239 } else {
240 B2DEBUG(25, "First t0 estimate not in proper range" << fitted_t0_first);
241 }
242
244 fitfuncSegmented.Draw("SAME");
245 }
246
247 // actually store this result ?
248 if (refitSuccess && fitresFull->IsValid()) {
249 const double fitted_t0 = -fitresFull->Parameter(3);
250 const double fitted_t0_error = fitresFull->Error(3);
251
252 const double norm_chi2 = fitresFull->Chi2() / double(fitresFull->Ndf());
253
254 B2DEBUG(25, "T0 fit with t0 " << fitted_t0 << " +- " << fitted_t0_error << " and normalized chi2 " << norm_chi2 << " and " <<
255 timingHistgram.GetEntries() << " hits");
256
257 // check if all the criteria required for a "good fit" have been met
258 if (norm_chi2 > m_param_rejectIfChiSquareLargerThan) {
259 B2DEBUG(25,
260 "T0 fit has too large Chi2 " << fitresFull->Chi2());
261 } else if (std::abs(fitted_t0_error) > m_param_rejectIfUncertaintyLargerThan) {
262 B2DEBUG(25,
263 "T0 fit has too large error " << fitted_t0_error);
264 } else {
265
266 // Since drift times are corrected for EventT0 (in RealisticTDCCountTranslator), if any other T0 modules were executed before, add the EventT0 offset back.
267 // This leads to "absolute" event T0 determination which should be consistent with other T0 modules.
268 double lastEventT0 = m_eventT0->hasEventT0() ? m_eventT0->getEventT0() : 0;
269 EventT0::EventT0Component eventT0Component(fitted_t0 + lastEventT0, fitted_t0_error, Const::CDC, "hit based", norm_chi2);
270 m_eventT0->addTemporaryEventT0(eventT0Component);
271 m_eventT0->setEventT0(eventT0Component);
272 m_wasSuccessful = true;
273 B2DEBUG(25,
274 "Successful t0 extraction with CDC hits: " << fitted_t0 << " +- " << fitted_t0_error);
275 }
276 } else {
277 B2DEBUG(25,
278 "Cannot fit t0 from CDC hits only. Won't set EventT0 for now.");
279 }
280 } else {
281 B2DEBUG(25,
282 "Cannot extract background or signal segment because fit failed. Won't set EventT0 for now.");
283 }
284
286 canvas.Draw();
287 canvas.SaveAs(debugImageName.c_str());
288 }
289}
void initialize() override final
Initialize the event meta data.
std::string getDescription() override final
Short description of the findlet.
float m_param_rejectIfChiSquareLargerThan
largest allowable chi2 value
int m_param_binCountTimeHistogram
number of bins in the timing histogram
bool m_param_rejectByBackgroundFlag
don't use hits for the fit which have been flagged as background
void apply(std::vector< TrackFindingCDC::CDCWireHit > &inputWireHits) override final
Collects all Wire Hits and executes the t0 fit.
StoreObjPtr< EventMetaData > m_eventMetaData
access to event nr for debugging purposes
float m_param_searchWindow
the window (+-m_param_searchWindow) in ns where to search for the best t0
bool m_param_storeAllFits
store an image of the histogram and fit result, for debug purpose only
bool m_param_rejectIfNotTakenFlag
don't use hits for the fit which have not been assigned to any track
void exposeParameters(ModuleParamList *moduleParamList, const std::string &prefix) override final
Expose the parameters to a module.
float m_param_rejectIfUncertaintyLargerThan
largest allowable uncertainty value
float m_param_fitWindow
the space (+-m_param_fitWindow) in ns used to fit the t0
unsigned int m_param_minHitCount
minimum number of hits
float m_param_refitWindow
the width of the window in ns (+- m_param_refitWindow) used to refit the final t0 value
The Module parameter list class.
static ScopeGuard guardBatchMode(bool batchMode=true)
Create a ScopeGuard to turn ROOT into batch mode and restore the initial batch mode status after the ...
Definition ScopeGuard.h:352
void addParameter(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module list.
virtual void exposeParameters(ModuleParamList *moduleParamList, const std::string &prefix) override
Abstract base class for different kinds of events.
Structure for storing the extracted event t0s together with its detector and its uncertainty.
Definition EventT0.h:33