9#include <svd/modules/svdTimeGrouping/SVDTimeGroupingModule.h>
12#include <framework/logging/Logger.h>
13#include <framework/utilities/FileSystem.h>
16#include <svd/geometry/SensorInfo.h>
17#include <svd/dataobjects/SVDEventInfo.h>
37 "SVDEventInfo collection name.", std::string(
"SVDEventInfo"));
41 "use SVDRecoConfiguration from DB",
bool(
true));
43 "if true, module is enabled for 6-sample DAQ mode",
bool(
false));
45 "if true, module is enabled for 3-sample DAQ mode",
bool(
false));
47 "use SVDTimeGroupingConfiguration from DB",
bool(
true));
50 addParam(
"tRangeLow",
m_usedPars.tRange[0],
"This sets the x- range of histogram [ns].",
52 addParam(
"tRangeHigh",
m_usedPars.tRange[1],
"This sets the x+ range of histogram [ns].",
55 "Time bin width is 1/rebinningFactor ns. Disables the module if set zero",
58 "Number of Gaussian sigmas (= hardcoded resolutions) used to fill the time histogram for each cluster.",
63 "Lower limit of cluster time sigma for the fit for the peak-search [ns].",
66 "Upper limit of cluster time sigma for the fit for the peak-search [ns].",
69 "half width of the range in which the fit for the peak-search is performed [ns].",
72 "Evaluate and remove gauss upto N sigma.",
75 "Minimum fraction of candidates in a peak (wrt to the highest peak) considered for fitting in the peak-search.",
78 "Maximum number of groups to be accepted.",
83 "Expected time of the signal [ns].",
86 "Expected low range of signal hits [ns].",
89 "Expected high range of signal hits [ns].",
92 "Group prominence is weighted with exponential weight with a lifetime defined by this parameter [ns].",
97 "Accept clusters upto N sigma.",
100 "Write group info into SVDClusters.",
105 "Assign groups to under and overflow.",
109 m_usedPars.clsSigma[0][0] = {2.0417, 2.3606, 2.1915, 1.9810, 1.8042, 1.6205};
110 m_usedPars.clsSigma[0][1] = {3.5880, 3.4526, 2.9363, 2.6833, 2.5342, 2.2895};
111 m_usedPars.clsSigma[1][0] = {2.1069, 2.0530, 1.9895, 1.8720, 1.6453, 1.5905};
112 m_usedPars.clsSigma[1][1] = {3.3919, 2.2280, 2.1177, 2.0852, 1.9968, 1.9914};
113 m_usedPars.clsSigma[2][0] = {1.6863, 1.9920, 1.8498, 1.7737, 1.6320, 1.5629};
114 m_usedPars.clsSigma[2][1] = {3.2798, 3.2243, 2.9404, 2.7911, 2.6331, 2.5666};
124 B2FATAL(
"no valid configuration found for SVD reconstruction");
126 B2DEBUG(20,
"SVDRecoConfiguration: from now on we are using " <<
m_recoConfig->get_uniqueID());
133 B2INFO(
"SVDTimeGrouping : SVDCluster groupId is assigned for 6-sample DAQ mode.");
135 B2INFO(
"SVDTimeGrouping : SVDCluster groupId is not assigned for 6-sample DAQ mode.");
138 B2INFO(
"SVDTimeGrouping : SVDCluster groupId is assigned for 3-sample DAQ mode.");
140 B2INFO(
"SVDTimeGrouping : SVDCluster groupId is not assigned for 3-sample DAQ mode.");
146 B2FATAL(
"no valid configuration found for SVD reconstruction");
148 B2DEBUG(20,
"SVDRecoConfiguration: from now on we are using " <<
m_recoConfig->get_uniqueID());
150 TString timeRecoWith6SamplesAlgorithm =
m_recoConfig->getTimeRecoWith6Samples();
151 TString timeRecoWith3SamplesAlgorithm =
m_recoConfig->getTimeRecoWith3Samples();
154 B2FATAL(
"no valid configuration found for SVDTimeGrouping");
156 B2DEBUG(20,
"SVDTimeGroupingConfiguration: from now on we are using " <<
m_groupingConfig->get_uniqueID());
170 B2DEBUG(20,
"SVDTimeGroupingModule \nsvdClusters: " <<
m_svdClusters.getName());
186 if (!eventinfo) B2ERROR(
"No SVDEventInfo!");
187 int numberOfAcquiredSamples = eventinfo->getNSamples();
190 if (numberOfAcquiredSamples == 6) {
193 }
else if (numberOfAcquiredSamples == 3) {
207 std::vector<GroupInfo> groupInfoVector;
235 double tmpRange[2] = {std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()};
236 for (
int ij = 0; ij < totClusters; ij++) {
238 if (std::isnan(tmpRange[0]) || clsTime > tmpRange[0]) tmpRange[0] = clsTime;
239 if (std::isnan(tmpRange[1]) || clsTime < tmpRange[1]) tmpRange[1] = clsTime;
243 if (tRangeHigh > tmpRange[0]) tRangeHigh = tmpRange[0];
244 if (tRangeLow < tmpRange[1]) tRangeLow = tmpRange[1];
246 int nBin = tRangeHigh - tRangeLow;
247 if (nBin < 1) nBin = 1;
249 if (nBin < 2) nBin = 2;
250 B2DEBUG(21,
"tRange: [" << tRangeLow <<
"," << tRangeHigh <<
"], nBin: " << nBin);
252 hist = TH1D(
"h_clsTime",
"h_clsTime", nBin, tRangeLow, tRangeHigh);
253 hist.GetXaxis()->SetLimits(tRangeLow, tRangeHigh);
255 for (
int ij = 0; ij < totClusters; ij++) {
259 double gSigma = (clsSize >= int(
m_usedPars.clsSigma[sType][isUcls].size()) ?
261 m_usedPars.clsSigma[sType][isUcls][clsSize - 1]);
275 double maxIntegral = 0.;
278 int roughCleaningCounter = 0;
282 int maxBin = hist.GetMaximumBin();
283 double maxBinCenter = hist.GetBinCenter(maxBin);
284 double maxBinContent = hist.GetBinContent(maxBin);
288 maxBinCenter >
m_usedPars.expectedSignalTime[0] && maxBinCenter <
m_usedPars.expectedSignalTime[2])
289 maxPeak = maxBinContent;
291 if (maxPeak != 0 && maxBinContent < maxPeak *
m_usedPars.fracThreshold) { amDone =
true;
continue;}
296 TF1 ngaus(
"ngaus",
myGaus,
297 hist.GetXaxis()->GetXmin(), hist.GetXaxis()->GetXmax(), 3);
300 double maxPar0 = maxBinContent * 2.50662827 *
m_usedPars.fitRangeHalfWidth;
301 ngaus.SetParameter(0, maxBinContent);
302 ngaus.SetParLimits(0,
305 ngaus.SetParameter(1, maxBinCenter);
306 ngaus.SetParLimits(1,
307 maxBinCenter -
m_usedPars.fitRangeHalfWidth * 0.2,
308 maxBinCenter +
m_usedPars.fitRangeHalfWidth * 0.2);
309 ngaus.SetParameter(2,
m_usedPars.fitRangeHalfWidth);
310 ngaus.SetParLimits(2,
316 int status = hist.Fit(
"ngaus",
"NQ0",
"",
324 ngaus.GetParameter(0),
325 ngaus.GetParameter(1),
326 std::fabs(ngaus.GetParameter(2))
334 if (roughCleaningCounter++ >
m_usedPars.maxGroups) amDone =
true;
339 if (maxPeak != 0 && maxIntegral == 0) maxIntegral = pars[0];
341 if (maxIntegral != 0 && pars[0] < maxIntegral *
m_usedPars.fracThreshold) { amDone =
true;
continue;}
348 groupInfoVector.push_back(
GroupInfo(pars[0], pars[1], pars[2]));
349 B2DEBUG(21,
" group " <<
int(groupInfoVector.size())
350 <<
" pars[0] " << pars[0] <<
" pars[1] " << pars[1] <<
" pars[2] " << pars[2]);
352 if (
int(groupInfoVector.size()) >=
m_usedPars.maxGroups) { amDone =
true;
continue;}
357 if (roughCleaningCounter++ >
m_usedPars.maxGroups) amDone =
true;
369 for (
int ij =
int(groupInfoVector.size()) - 2; ij >= 0; ij--) {
370 keyGroup = groupInfoVector[ij];
371 double keyGroupIntegral = std::get<0>(keyGroup);
372 double keyGroupCenter = std::get<1>(keyGroup);
373 bool isKeyGroupSignal =
true;
374 if (keyGroupIntegral != 0. &&
375 (keyGroupCenter <
m_usedPars.expectedSignalTime[0] || keyGroupCenter >
m_usedPars.expectedSignalTime[2]))
376 isKeyGroupSignal =
false;
377 if (isKeyGroupSignal)
continue;
380 while (kj <
int(groupInfoVector.size())) {
381 double otherGroupIntegral = std::get<0>(groupInfoVector[kj]);
382 double otherGroupCenter = std::get<1>(groupInfoVector[kj]);
383 bool isOtherGroupSignal =
true;
384 if (otherGroupIntegral != 0. &&
385 (otherGroupCenter <
m_usedPars.expectedSignalTime[0] || otherGroupCenter >
m_usedPars.expectedSignalTime[2]))
386 isOtherGroupSignal =
false;
387 if (!isOtherGroupSignal && (otherGroupIntegral > keyGroupIntegral))
break;
388 groupInfoVector[kj - 1] = groupInfoVector[kj];
391 groupInfoVector[kj - 1] = keyGroup;
400 for (
int ij = 1; ij < int(groupInfoVector.size()); ij++) {
401 keyGroup = groupInfoVector[ij];
402 double keyGroupIntegral = std::get<0>(keyGroup);
403 if (keyGroupIntegral <= 0)
break;
404 double keyGroupCenter = std::get<1>(keyGroup);
405 bool isKeyGroupSignal =
true;
406 if (keyGroupIntegral > 0 &&
407 (keyGroupCenter <
m_usedPars.expectedSignalTime[0] || keyGroupCenter >
m_usedPars.expectedSignalTime[2]))
408 isKeyGroupSignal =
false;
409 if (!isKeyGroupSignal)
break;
411 double keyWt = keyGroupIntegral * TMath::Exp(-std::fabs(keyGroupCenter -
m_usedPars.expectedSignalTime[1]) /
415 double otherGroupIntegral = std::get<0>(groupInfoVector[kj]);
416 double otherGroupCenter = std::get<1>(groupInfoVector[kj]);
417 double grWt = otherGroupIntegral * TMath::Exp(-std::fabs(otherGroupCenter -
m_usedPars.expectedSignalTime[1]) /
419 if (grWt > keyWt)
break;
420 groupInfoVector[kj + 1] = groupInfoVector[kj];
423 groupInfoVector[kj + 1] = keyGroup;
432 double tRangeLow = hist.GetXaxis()->GetXmin();
433 double tRangeHigh = hist.GetXaxis()->GetXmax();
436 if (
int(groupInfoVector.size()) == 0)
437 for (
int jk = 0; jk < totClusters; jk++)
442 for (
int ij = 0; ij < int(groupInfoVector.size()); ij++) {
445 std::get<0>(groupInfoVector[ij]),
446 std::get<1>(groupInfoVector[ij]),
447 std::get<2>(groupInfoVector[ij])
450 if (pars[2] == 0 && ij !=
int(groupInfoVector.size()) - 1)
continue;
455 double lowestAcceptedTime = pars[1] -
m_usedPars.acceptSigmaN * pars[2];
456 double highestAcceptedTime = pars[1] +
m_usedPars.acceptSigmaN * pars[2];
457 if (lowestAcceptedTime < tRangeLow) lowestAcceptedTime = tRangeLow;
458 if (highestAcceptedTime > tRangeHigh) highestAcceptedTime = tRangeHigh;
459 B2DEBUG(21,
" group " << ij
460 <<
" lowestAcceptedTime " << lowestAcceptedTime
461 <<
" highestAcceptedTime " << highestAcceptedTime);
464 for (
int jk = 0; jk < totClusters; jk++) {
468 clsTime >= lowestAcceptedTime && clsTime <= highestAcceptedTime) {
478 B2DEBUG(29,
" accepted cluster " << jk
479 <<
" clsTime " << clsTime
480 <<
" GroupId " <<
m_svdClusters[jk]->getTimeGroupId().back());
484 B2DEBUG(29,
" rejected cluster " << jk
485 <<
" clsTime " << clsTime);
487 if (ij ==
int(groupInfoVector.size()) - 1 &&
490 if (
m_usedPars.includeOutOfRangeClusters && clsTime < tRangeLow)
492 else if (
m_usedPars.includeOutOfRangeClusters && clsTime > tRangeHigh)
498 B2DEBUG(29,
" leftover cluster " << jk
499 <<
" GroupId " <<
m_svdClusters[jk]->getTimeGroupId().back());
void setDescription(const std::string &description)
Sets the description of the module.
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
bool m_useParamFromDB
if true use the configuration from SVDTimeGroupingConfiguration DB.
std::string m_svdClustersName
SVDCluster collection name.
SVDTimeGroupingParameters m_usedParsIn3Samples
module parameter values for 3-sample DAQ taken from SVDTimeGroupingConfiguration dbobject.
SVDTimeGroupingParameters m_usedPars
module parameter values used.
virtual void initialize() override
Init the module.
StoreArray< SVDCluster > m_svdClusters
the storeArray for svdClusters as member, is faster than recreating link for each event
virtual void event() override
EventWise jobs.
SVDTimeGroupingParameters m_usedParsIn6Samples
module parameter values for 6-sample DAQ taken from SVDTimeGroupingConfiguration dbobject.
void resizeToMaxSize(std::vector< GroupInfo > &groupInfoVector)
increase the size of vector to max, this helps in sorting
void sortBackgroundGroups(std::vector< GroupInfo > &groupInfoVector)
Sort Background Groups.
bool m_forceGroupingFromDB
if true use configuration from the SVDRecConfiguration DB.
void beginRun() override
configure
void assignGroupIdsToClusters(TH1D &hist, std::vector< GroupInfo > &groupInfoVector)
Assign groupId to the clusters.
std::string m_svdEventInfoName
Name of the collection to use for the SVDEventInfo.
void createAndFillHistorgram(TH1D &hist)
Create Histogram and Fill cluster time in it.
void sortSignalGroups(std::vector< GroupInfo > &groupInfoVector)
Sort Signals.
SVDTimeGroupingModule()
Constructor.
void searchGausPeaksInHistogram(TH1D &hist, std::vector< GroupInfo > &groupInfoVector)
Find Gaussian components in a Histogram.
DBObjPtr< SVDRecoConfiguration > m_recoConfig
SVD Reconstruction Configuration payload.
bool m_isEnabledIn6Samples
Enables the module if true for 6-sample DAQ mode.
DBObjPtr< SVDTimeGroupingConfiguration > m_groupingConfig
SVDTimeGrouping Configuration payload.
bool m_isEnabledIn3Samples
Enables the module if true for 3-sample DAQ mode.
Type-safe access to single objects in the data store.
bool isValid() const
Check whether the object was created.
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
std::tuple< double, double, double > GroupInfo
typedef to be used to store Gauss parameters (integral, center, sigma)
int getSensorType(const VxdID &sensorID)
Get Sensor Type of SVD sensors.
double myGaus(const double *x, const double *par)
Gauss function to be used in the fit.
void subtractGausFromHistogram(TH1D &hist, const double &integral, const double ¢er, const double &sigma, const double &sigmaN)
Subtract a Gaussian from a histogram.
void addGausToHistogram(TH1D &hist, const double &integral, const double ¢er, const double &sigma, const double &sigmaN, const bool &isAddition=true)
Add (or Subtract) a Gaussian to (or from) a histogram.
Abstract base class for different kinds of events.