Belle II Software  release-05-01-25
SVDCoGTimeEstimatorModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2017 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Michael De Nuccio, Giulia Casarosa *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <svd/modules/svdReconstruction/SVDCoGTimeEstimatorModule.h>
12 
13 using namespace Belle2;
14 using namespace std;
15 
16 //-----------------------------------------------------------------
17 // Register the Module
18 //-----------------------------------------------------------------
19 REG_MODULE(SVDCoGTimeEstimator)
20 
21 //-----------------------------------------------------------------
22 // Implementation
23 //-----------------------------------------------------------------
24 
26 {
27  setDescription("From SVDShaperDigit to SVDRecoDigit. Strip charge is evaluated as the max of the 6 samples; hit time is evaluated as a corrected Centre of Gravity (CoG) time.");
28  setPropertyFlags(c_ParallelProcessingCertified);
29 
30  addParam("SVDEventInfo", m_svdEventInfoName,
31  "SVDEventInfo name", string(""));
32  addParam("ShaperDigits", m_storeShaperDigitsName,
33  "ShaperDigits collection name", string(""));
34  addParam("RecoDigits", m_storeRecoDigitsName,
35  "RecoDigits collection name", string(""));
36  addParam("StripPeakTimeCorrection", m_corrPeakTime,
37  "Correct for the different peaking times of the strips, obtained from local run calibration", true);
38  addParam("CalibrationWithEventT0", m_calEventT0,
39  "Use the timing informations of the EventT0 in order to calibrate the CoG.",
40  true);
41 
42 }
43 
44 
46 {
47 }
48 
49 
51 {
52  m_storeTrueHits.isOptional();
53  m_storeMCParticles.isOptional();
54 
55  //Inizialization of needed store array
56  m_storeShaper.isRequired(m_storeShaperDigitsName);
57 
58  if (!m_storeSVDEvtInfo.isOptional(m_svdEventInfoName)) m_svdEventInfoName = "SVDEventInfoSim";
59  m_storeSVDEvtInfo.isRequired(m_svdEventInfoName);
60 
61  //Initialize the new RecoDigit
62  m_storeReco.registerInDataStore(m_storeRecoDigitsName, DataStore::c_ErrorIfAlreadyRegistered);
63 
64  RelationArray relRecoDigitShaperDigits(m_storeReco, m_storeShaper);
65  relRecoDigitShaperDigits.registerInDataStore();
66  RelationArray relRecoDigitTrueHits(m_storeReco, m_storeTrueHits);
67  RelationArray relRecoDigitMCParticles(m_storeReco, m_storeMCParticles);
68  RelationArray relShaperDigitTrueHits(m_storeShaper, m_storeTrueHits);
69  RelationArray relShaperDigitMCParticles(m_storeShaper, m_storeMCParticles);
70 
71  //Relations to simulation objects only if the ancestor relations exist
72  if (relShaperDigitTrueHits.isOptional())
73  relRecoDigitTrueHits.registerInDataStore();
74  if (relShaperDigitMCParticles.isOptional())
75  relRecoDigitMCParticles.registerInDataStore();
76 
77  //Store names to speed up creation later
78  m_storeTrueHitsName = m_storeTrueHits.getName();
79  m_storeMCParticlesName = m_storeMCParticles.getName();
80 
81  m_relRecoDigitShaperDigitName = relRecoDigitShaperDigits.getName();
82  m_relRecoDigitTrueHitName = relRecoDigitTrueHits.getName();
83  m_relRecoDigitMCParticleName = relRecoDigitMCParticles.getName();
84  m_relShaperDigitTrueHitName = relShaperDigitTrueHits.getName();
85  m_relShaperDigitMCParticleName = relShaperDigitMCParticles.getName();
86 
87  B2DEBUG(1, " 1. COLLECTIONS:");
88  B2DEBUG(1, " --> MCParticles: " << m_storeMCParticlesName);
89  B2DEBUG(1, " --> Digits: " << m_storeShaperDigitsName);
90  B2DEBUG(1, " --> RecoDigits: " << m_storeRecoDigitsName);
91  B2DEBUG(1, " --> TrueHits: " << m_storeTrueHitsName);
92  B2DEBUG(1, " --> DigitMCRel: " << m_relShaperDigitMCParticleName);
93  B2DEBUG(1, " --> RecoDigitMCRel: " << m_relRecoDigitMCParticleName);
94  B2DEBUG(1, " --> RecoDigitDigitRel: " << m_relRecoDigitShaperDigitName);
95  B2DEBUG(1, " --> DigitTrueRel: " << m_relShaperDigitTrueHitName);
96  B2DEBUG(1, " --> RecoDigitTrueRel: " << m_relRecoDigitTrueHitName);
97 
98 }
100 {
101 }
102 
103 
105 {
106 
108  std::vector<float> probabilities = {0.5};
109 
110  // If no digits or no SVDEventInfo, nothing to do
111  if (!m_storeShaper || !m_storeShaper.getEntries() || !m_storeSVDEvtInfo.isValid()) return;
112 
113  SVDModeByte modeByte = m_storeSVDEvtInfo->getModeByte();
114  size_t nDigits = m_storeShaper.getEntries();
115 
116  RelationArray relRecoDigitShaperDigit(m_storeReco, m_storeShaper,
117  m_relRecoDigitShaperDigitName);
118  if (relRecoDigitShaperDigit) relRecoDigitShaperDigit.clear();
119 
120 
121  RelationArray relShaperDigitMCParticle(m_storeShaper, m_storeMCParticles, m_relShaperDigitMCParticleName);
122  RelationArray relShaperDigitTrueHit(m_storeShaper, m_storeTrueHits, m_relShaperDigitTrueHitName);
123 
124  RelationArray relRecoDigitMCParticle(m_storeReco, m_storeMCParticles,
125  m_relRecoDigitMCParticleName);
126  if (relRecoDigitMCParticle) relRecoDigitMCParticle.clear();
127 
128 
129  RelationArray relRecoDigitTrueHit(m_storeReco, m_storeTrueHits,
130  m_relRecoDigitTrueHitName);
131  if (relRecoDigitTrueHit) relRecoDigitTrueHit.clear();
132 
133  //Build lookup tables for relations
134  createRelationLookup(relShaperDigitMCParticle, m_mcRelation, nDigits);
135  createRelationLookup(relShaperDigitTrueHit, m_trueRelation, nDigits);
136 
137  //start loop on SVDSHaperDigits
139 
140  m_NumberOfAPVSamples = m_storeSVDEvtInfo->getNSamples();
141 
142  B2DEBUG(1, "number of APV samples = " << m_NumberOfAPVSamples);
143 
144  for (const SVDShaperDigit& shaper : m_storeShaper) {
145 
146  m_StopCreationReco = false;
147 
148 
149  if (m_StopCreationReco)
150  continue;
151 
152  samples_vec = shaper.getSamples();
153 
154  //retrieve the VxdID, sensor and cellID of the current RecoDigit
155  VxdID thisSensorID = shaper.getSensorID();
156  bool thisSide = shaper.isUStrip();
157  int thisCellID = shaper.getCellID();
158 
159  //call of the functions doomed to calculate the required quantities
160  m_weightedMeanTime = CalculateWeightedMeanPeakTime(samples_vec);
161  if (m_StopCreationReco)
162  continue;
163  m_amplitude = CalculateAmplitude(samples_vec);
164  m_amplitudeError = CalculateAmplitudeError(thisSensorID, thisSide, thisCellID);
165 
166  //need the amplitudeError in ADC as the noise of the strip to computer the error on time
167  m_weightedMeanTimeError = CalculateWeightedMeanPeakTimeError(samples_vec);
168 
169  m_chi2 = CalculateChi2();
170 
171  //check too high ADC
172  if (m_amplitude > 255)
173  B2DEBUG(10, "ERROR: m_amplitude = " << m_amplitude << ", should be <= 255");
174 
175  //CALIBRATION
176  //convert ADC into #e- and apply offset to shift estimated peak time to hit time (to be completed)
177  m_amplitude = m_PulseShapeCal.getChargeFromADC(thisSensorID, thisSide, thisCellID, m_amplitude);
178  m_amplitudeError = m_PulseShapeCal.getChargeFromADC(thisSensorID, thisSide, thisCellID, m_amplitudeError);
179 
180  if (m_corrPeakTime)
181  m_weightedMeanTime -= m_PulseShapeCal.getPeakTime(thisSensorID, thisSide, thisCellID);
182  SVDModeByte::baseType triggerBin = modeByte.getTriggerBin();
183 
184  if (m_calEventT0) {
185  m_weightedMeanTime = m_TimeCal.getCorrectedTime(thisSensorID, thisSide, thisCellID, m_weightedMeanTime, triggerBin);
186  m_weightedMeanTimeError = m_TimeCal.getCorrectedTimeError(thisSensorID, thisSide, thisCellID, m_weightedMeanTime,
187  m_weightedMeanTimeError, triggerBin);
188  }
189 
190  //check high charges and too high ADC
191  if (m_amplitude > 100000) {
192  B2DEBUG(42, "Charge = " << m_amplitude);
193  B2DEBUG(42, "corresponding ADC = " << m_PulseShapeCal.getADCFromCharge(thisSensorID, thisSide, thisCellID, m_amplitude));
194  B2DEBUG(42, "thisLayerNumber = " << thisSensorID.getLayerNumber());
195  B2DEBUG(42, "thisLadderNumber = " << thisSensorID.getLadderNumber());
196  B2DEBUG(42, "thisSensorNumber = " << thisSensorID.getSensorNumber());
197  B2DEBUG(42, "thisSide = " << thisSide);
198  B2DEBUG(42, "thisCellID = " << thisCellID);
199  B2DEBUG(42, "-----");
200  }
201 
202  //recording of the RecoDigit
203  m_storeReco.appendNew(SVDRecoDigit(shaper.getSensorID(), shaper.isUStrip(), shaper.getCellID(), m_amplitude, m_amplitudeError,
204  m_weightedMeanTime, m_weightedMeanTimeError, probabilities, m_chi2, modeByte));
205 
206  //Add digit to the RecoDigit->ShaperDigit relation list
207  int recoDigitIndex = m_storeReco.getEntries() - 1;
208  vector<pair<unsigned int, float> > digit_weights;
209  digit_weights.reserve(1);
210  digit_weights.emplace_back(shaper.getArrayIndex(), 1.0);
211  relRecoDigitShaperDigit.add(recoDigitIndex, digit_weights.begin(), digit_weights.end());
212 
213  // Finally, we save the RecoDigit and its relations.
214  map<unsigned int, float> mc_relations;
215  map<unsigned int, float> truehit_relations;
216 
217  // Store relations to MCParticles and SVDTrueHits
218  fillRelationMap(m_mcRelation, mc_relations, shaper.getArrayIndex());
219  fillRelationMap(m_trueRelation, truehit_relations, shaper.getArrayIndex());
220 
221  //Create relations to the cluster
222  if (!mc_relations.empty()) {
223  relRecoDigitMCParticle.add(recoDigitIndex, mc_relations.begin(), mc_relations.end());
224  }
225  if (!truehit_relations.empty()) {
226  relRecoDigitTrueHit.add(recoDigitIndex, truehit_relations.begin(), truehit_relations.end());
227  }
228 
229  }
230 }
231 
232 
234 {
235 }
236 
237 
239 {
240 }
241 
242 
244 {
245  float averagetime = 0;
246  float sumAmplitudes = 0;
247  //calculate weighted average time
248  for (int k = 0; k < m_NumberOfAPVSamples; k ++) {
249  averagetime += k * samples[k];
250  sumAmplitudes += samples[k];
251  }
252  if (sumAmplitudes != 0) {
253  averagetime /= (sumAmplitudes);
254  averagetime *= m_DeltaT;
255  } else {
256  averagetime = -1;
257  m_StopCreationReco = true;
258  B2WARNING("Trying to divide by 0 (ZERO)! Sum of amplitudes is NULL! Skipping this SVDShaperDigit!");
259  }
260 
261  return averagetime;
262 }
263 
265 {
266  float amplitude = 0;
267  //calculate amplitude
268  for (int k = 0; k < m_NumberOfAPVSamples; k ++) {
269  if (samples[k] > amplitude)
270  amplitude = samples[k];
271  }
272 
273  return amplitude;
274 }
275 
277 {
278 
279  //assuming that noise of the samples are totally UNcorrelated
280  //in MC this hypothesis is correct
281 
282  //sum of samples amplitudes
283  float Atot = 0;
284  //sum of time residuals squared
285  float tmpResSq = 0;
286 
287  for (int k = 0; k < m_NumberOfAPVSamples; k ++) {
288  Atot += samples[k];
289  tmpResSq += TMath::Power(k * m_DeltaT - m_weightedMeanTime, 2);
290  }
291 
292  return m_amplitudeError / Atot * TMath::Sqrt(tmpResSq);
293 
294 }
295 
296 float SVDCoGTimeEstimatorModule::CalculateAmplitudeError(VxdID ThisSensorID, bool ThisSide, int ThisCellID)
297 {
298  float stripnoise;
299  stripnoise = m_NoiseCal.getNoise(ThisSensorID, ThisSide, ThisCellID);
300 
301  return stripnoise;
302 }
303 
305 {
306  return 0.01;
307 }
308 
309 
310 
312  RelationLookup& lookup, size_t digits)
313 {
314  lookup.clear();
315  //If we don't have a relation we don't build a lookuptable
316  if (!relation) return;
317  //Resize to number of digits and set all values
318  lookup.resize(digits);
319  for (const auto& element : relation) {
320  lookup[element.getFromIndex()] = &element;
321  }
322 }
323 
325  std::map<unsigned int, float>& relation, unsigned int index)
326 {
327  //If the lookup table is not empty and the element is set
328  if (!lookup.empty() && lookup[index]) {
329  const RelationElement& element = *lookup[index];
330  const unsigned int size = element.getSize();
331  //Add all Relations to the map
332  for (unsigned int i = 0; i < size; ++i) {
333  //negative weights are from ignored particles, we don't like them and
334  //thus ignore them :D
335  if (element.getWeight(i) < 0) continue;
336  relation[element.getToIndex(i)] += element.getWeight(i);
337  }
338  }
339 }
340 
341 
342 
343 
344 
345 
346 
347 
348 
349 
350 
351 
352 
353 
354 
355 
356 
357 
358 
359 
360 
361 
362 
363 
364 
365 
366 
367 
368 
369 
370 
371 
Belle2::SVDCoGTimeEstimatorModule::~SVDCoGTimeEstimatorModule
virtual ~SVDCoGTimeEstimatorModule()
default destructor
Definition: SVDCoGTimeEstimatorModule.cc:45
Belle2::SVDCoGTimeEstimatorModule::endRun
virtual void endRun() override
This method is called if the current run ends.
Definition: SVDCoGTimeEstimatorModule.cc:233
Belle2::RelationArray::clear
void clear() override
Clear all elements from the relation.
Definition: RelationArray.h:279
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:72
Belle2::SVDCoGTimeEstimatorModule::CalculateAmplitudeError
float CalculateAmplitudeError(VxdID ThisSensorID, bool ThisSide, int ThisCellID)
Function to calculate the amplitude error as the noise of the strip.
Definition: SVDCoGTimeEstimatorModule.cc:296
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::SVDCoGTimeEstimatorModule::RelationLookup
std::vector< const RelationElement * > RelationLookup
Container for a RelationArray Lookup table.
Definition: SVDCoGTimeEstimatorModule.h:58
Belle2::SVDModeByte::getTriggerBin
baseType getTriggerBin() const
Get the triggerBin id.
Definition: SVDModeByte.h:150
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::SVDCoGTimeEstimatorModule::initialize
virtual void initialize() override
Initialize the SVDCoGTimeEstimator.
Definition: SVDCoGTimeEstimatorModule.cc:50
Belle2::VxdID::getLadderNumber
baseType getLadderNumber() const
Get the ladder id.
Definition: VxdID.h:108
Belle2::StoreAccessorBase::isOptional
bool isOptional(const std::string &name="")
Tell the DataStore about an optional input.
Definition: StoreAccessorBase.h:95
Belle2::StoreAccessorBase::getName
const std::string & getName() const
Return name under which the object is saved in the DataStore.
Definition: StoreAccessorBase.h:130
Belle2::SVDModeByte
Class to store SVD mode information.
Definition: SVDModeByte.h:79
Belle2::StoreAccessorBase::registerInDataStore
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Definition: StoreAccessorBase.h:54
Belle2::SVDCoGTimeEstimatorModule::terminate
virtual void terminate() override
This method is called at the end of the event processing.
Definition: SVDCoGTimeEstimatorModule.cc:238
Belle2::SVDCoGTimeEstimatorModule::event
virtual void event() override
This method is the core of the SVDCoGTimeEstimator.
Definition: SVDCoGTimeEstimatorModule.cc:104
Belle2::SVDCoGTimeEstimatorModule::CalculateWeightedMeanPeakTimeError
float CalculateWeightedMeanPeakTimeError(Belle2::SVDShaperDigit::APVFloatSamples samples)
Function to calculate the peak time error.
Definition: SVDCoGTimeEstimatorModule.cc:276
Belle2::SVDShaperDigit
The SVD ShaperDigit class.
Definition: SVDShaperDigit.h:46
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::RelationElement
Class to store a single element of a relation.
Definition: RelationElement.h:33
Belle2::SVDShaperDigit::APVFloatSamples
std::array< APVFloatSampleType, c_nAPVSamples > APVFloatSamples
array of APVFloatSampleType objects
Definition: SVDShaperDigit.h:63
Belle2::SVDCoGTimeEstimatorModule::fillRelationMap
void fillRelationMap(const RelationLookup &lookup, std::map< unsigned int, float > &relation, unsigned int index)
Add the relation from a given SVDShaperDigit index to a map.
Definition: SVDCoGTimeEstimatorModule.cc:324
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::SVDModeByte::baseType
uint8_t baseType
The base integer type for SVDModeByte.
Definition: SVDModeByte.h:82
Belle2::SVDCoGTimeEstimatorModule::CalculateChi2
float CalculateChi2()
Function to calculate chi2, that is not used here, so just set at 0.01.
Definition: SVDCoGTimeEstimatorModule.cc:304
Belle2::RelationArray::add
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
Definition: RelationArray.h:291
Belle2::SVDCoGTimeEstimatorModule::beginRun
virtual void beginRun() override
Called when entering a new run.
Definition: SVDCoGTimeEstimatorModule.cc:99
Belle2::VxdID::getSensorNumber
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:110
Belle2::SVDCoGTimeEstimatorModule::CalculateAmplitude
float CalculateAmplitude(Belle2::SVDShaperDigit::APVFloatSamples samples)
Function to calculate the amplitude of the shaper, obtained as the largest of the 6 samples.
Definition: SVDCoGTimeEstimatorModule.cc:264
Belle2::DataStore::c_ErrorIfAlreadyRegistered
@ c_ErrorIfAlreadyRegistered
If the object/array was already registered, produce an error (aborting initialisation).
Definition: DataStore.h:74
Belle2::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:106
Belle2::SVDRecoDigit
The SVD RecoDigit class.
Definition: SVDRecoDigit.h:54
Belle2::SVDCoGTimeEstimatorModule::CalculateWeightedMeanPeakTime
float CalculateWeightedMeanPeakTime(Belle2::SVDShaperDigit::APVFloatSamples samples)
Function to calculate the peak time, obtained as the weighted mean of the time of the samples,...
Definition: SVDCoGTimeEstimatorModule.cc:243
Belle2::SVDCoGTimeEstimatorModule::createRelationLookup
void createRelationLookup(const RelationArray &relation, RelationLookup &lookup, size_t digits)
Create lookup maps for relations FIXME: This has to be significantly simplified here,...
Definition: SVDCoGTimeEstimatorModule.cc:311
Belle2::SVDCoGTimeEstimatorModule
This module builds the SVDRecoDigits (calibrated and fitted strips) from the SVDShaperDigits.
Definition: SVDCoGTimeEstimatorModule.h:54