Belle II Software  release-05-02-19
TOPXTalkChargeShareSetterModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2018 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Rikuya Okuto *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 //own include
12 #include <top/modules/TOPXTalkChargeShareSetter/TOPXTalkChargeShareSetterModule.h>
13 
14 // framework - DataStore
15 #include <framework/datastore/StoreArray.h>
16 
17 // Dataobject classes
18 #include <top/dataobjects/TOPRawDigit.h>
19 #include <top/dataobjects/TOPRawWaveform.h>
20 #include <top/dataobjects/TOPDigit.h>
21 
22 #include <TMath.h>
23 
24 using namespace std;
25 
26 using namespace Belle2;
27 
28 //-----------------------------------------------------------------
29 // Register the Module
30 //-----------------------------------------------------------------
31 REG_MODULE(TOPXTalkChargeShareSetter)
32 
33 //-----------------------------------------------------------------
34 // Implementation
35 //-----------------------------------------------------------------
36 
38 {
39  // Set module properties
40  setDescription("Crosstalk & chargeshare flag setter");
41  setPropertyFlags(c_ParallelProcessingCertified);
42 
43  // Add parameter
44  addParam("sumChargeShare", m_sumChargeShare,
45  "sum up charge of PrimaryChargeShare and SecondaryChargeShare", (bool)false);
46  addParam("timeCut", m_timeCut,
47  "cut range of hittiming for chargeshare flag [ns]", (float) 1);
48  addParam("preValleyDepthLoose", m_preValleyDepthLoose,
49  "loose threshold for depth of pre valley [ADC counts], for corss talk identification", 20);
50  addParam("preValleyDepthTight", m_preValleyDepthTight,
51  "tight threshold for depth of pre valley [ADC counts], identified as cross talk with loose threshold for the second peak amplitude",
52  50);
53  addParam("2ndPeakAmplitudeLoose", m_2ndPeakAmplitudeLoose,
54  "loose threshold for amplitude of the second peak [ADC counts] for cross talk identification. Defined as ADC count difference between the valley just after the main peak and the second peak, Used when the \"preValleyDepthTight\" was satisfied",
55  30);
56  addParam("2ndPeakAmplitudeRatioTight", m_2ndPeakAmplitudeRatioTight,
57  "tight threshold for amplitude ratio of the second peak to the main peak height [ADC counts]", 0.2);
58  addParam("nSampleBefore", m_nSampleBefore,
59  "the number of samples by which the pre-valley should exist from the CFD timing, used for cross talk identification."
60  "Set negative value not to check existence of the pre-valley.", 5);
61  addParam("nSampleAfter", m_nSampleAfter,
62  "the number of samples by which the second peak should exist from the CFD timing, used for cross talk identification."
63  "Set negative value not to check existence of the second peak.", 10);
64 
65 }
66 
67 TOPXTalkChargeShareSetterModule::~TOPXTalkChargeShareSetterModule() {}
68 
69 void TOPXTalkChargeShareSetterModule::initialize()
70 {
71  StoreArray<TOPDigit> digits;
72  digits.isRequired();
73 
74  if (m_preValleyDepthLoose < 0 || m_preValleyDepthTight < 0 || m_nSampleBefore < 0) m_checkPreValleyForXTalkId = false;
75  if (m_2ndPeakAmplitudeLoose < 0 || m_2ndPeakAmplitudeRatioTight < 0 || m_nSampleAfter < 0) m_checkPostValleyForXTalkId = false;
76 }
77 
78 void TOPXTalkChargeShareSetterModule::beginRun()
79 {
80 }
81 
82 void TOPXTalkChargeShareSetterModule::event()
83 {
84  StoreArray<TOPDigit> digits;
85 
86  //Set Cross talk events and charge Share events
87  //
88  std::map<int, TOPDigit*> hitInfoMap;
89  for (auto& digit : digits) {
90 
91  if (digit.getHitQuality() == TOPDigit::c_Junk) continue;
92 
93  const auto* rawDigit = digit.getRelated<TOPRawDigit>();
94  if (!rawDigit) continue;
95 
96  const auto* waveform = rawDigit->getRelated<TOPRawWaveform>();
97  if (waveform and isCrossTalk(waveform->getWaveform(), TMath::FloorNint(digit.getRawTime()), digit.getPulseHeight()))
98  digit.setHitQuality(TOPDigit::c_CrossTalk);
99  else {
100  //identify a ringing which follows a cross talk hit in the same channle as cross talk
101  int slotId = digit.getModuleID();
102  unsigned int channelId = digit.getChannel();
103  double rawTime = digit.getRawTime();
104  for (auto& digit2 : digits) {
105  if (digit2.getChannel() != channelId || digit2.getModuleID() != slotId) continue;
106  if (digit2.getHitQuality() == TOPDigit::c_CrossTalk
107  and TMath::Abs(rawTime - digit2.getRawTime() - m_nCrossTalkRingingSamples / 2.) < m_nCrossTalkRingingSamples / 2.)
108  digit.setHitQuality(TOPDigit::c_CrossTalk);
109  }//for(digit2)
110  }//if(waveform and isCrossTalk) else
111 
112  if (digit.getHitQuality() == TOPDigit::c_CrossTalk) continue;
113 
114  int globalPixelId = (digit.getModuleID() - 1) * 512 + digit.getPixelID() - 1;
115  while (hitInfoMap.count(globalPixelId) > 0)globalPixelId += 10000;
116  hitInfoMap[globalPixelId] = &digit;
117  }//for(digit)
118 
119 
120 
121 
122  for (auto& digit : digits) {
123 
124  if (digit.getHitQuality() == TOPDigit::c_Junk
125  || digit.getHitQuality() == TOPDigit::c_CrossTalk) continue;
126 
127  short pixelId = digit.getPixelID();
128  short slotId = digit.getModuleID();
129  short pmtId = digit.getPMTNumber();
130  double hitTime = digit.getTime();
131  double charge = digit.getIntegral();
132  bool isPrimaryChargeShare = true;
133  std::vector<TOPDigit*> vSecondaryCandidates;
134 
135  int adjacentPixelIds[] = { pixelId - 1 - c_NPixelsPerRow, pixelId - c_NPixelsPerRow, pixelId + 1 - c_NPixelsPerRow, pixelId + 1,
136  pixelId + 1 + c_NPixelsPerRow, pixelId + c_NPixelsPerRow, pixelId - 1 + c_NPixelsPerRow, pixelId - 1
137  };
138 
139  for (const auto& adjacentPixelId : adjacentPixelIds) {
140  if (adjacentPixelId > 0 && adjacentPixelId <= 512) {
141  int globalPixelId = (slotId - 1) * 512 + adjacentPixelId - 1;
142  while (hitInfoMap.count(globalPixelId) > 0) {
143 
144  if (pmtId != hitInfoMap[globalPixelId]->getPMTNumber()
145  or TMath::Abs(hitTime - hitInfoMap[globalPixelId]->getTime()) > m_timeCut) {
146  globalPixelId += 10000; continue;
147  }
148 
149  float adjacentIntegral = hitInfoMap[globalPixelId]->getIntegral();
150 
151  if (charge > adjacentIntegral
152  or (charge == adjacentIntegral && pixelId > hitInfoMap[globalPixelId]->getPixelID())) {
153  vSecondaryCandidates.push_back(hitInfoMap[globalPixelId]);
154  hitInfoMap[globalPixelId]->setSecondaryChargeShare();
155  } else {
156  isPrimaryChargeShare = false;
157  break;
158  }
159  globalPixelId += 10000;
160  }//while(hitInfoMap.count(globalPixelId)>0)
161  if (!isPrimaryChargeShare)break;
162  }//if(adjacentPixelId exist)
163  }//for(adjacentPixelIds)
164 
165  if (isPrimaryChargeShare && vSecondaryCandidates.size() > 0) {
166  digit.setPrimaryChargeShare();
167 
168  if (m_sumChargeShare) {
169  for (auto& vSecondaryCandidate : vSecondaryCandidates) {
170  digit.setPulseHeight(digit.getPulseHeight() + vSecondaryCandidate->getPulseHeight());
171  vSecondaryCandidate->setPulseHeight(0.0);
172 
173  digit.setIntegral(digit.getIntegral() + vSecondaryCandidate->getIntegral());
174  vSecondaryCandidate->setIntegral(0.0);
175  }//for(vSecondaryCandidates)
176  }
177  }
178  }//for(digits)
179 }//event()
180 
181 void TOPXTalkChargeShareSetterModule::endRun()
182 {
183 }
184 
185 void TOPXTalkChargeShareSetterModule::terminate()
186 {
187 }
188 
189 
190 bool TOPXTalkChargeShareSetterModule::isCrossTalk(std::vector<short> wfm, int iRawTime, short height)
191 {
192 
193  int nWfmSampling = wfm.size();
194  for (int iWin = 0 ; iWin < 16 ; iWin++) {
195  //int jRawTime = iRawTime + (iWin - 20) * (TOPRawDigit::c_WindowSize)/4;//scan offset by 16-sample step
196  int jRawTime = iRawTime - TMath::FloorNint(iRawTime / (TOPRawDigit::c_WindowSize)) * (TOPRawDigit::c_WindowSize) +
197  (TOPRawDigit::c_WindowSize) / 4 * iWin;
198  if (jRawTime > 0 && jRawTime < nWfmSampling - 1 && wfm[jRawTime] < height / 2. && wfm[jRawTime + 1] > height / 2.) {
199  bool preValleyExist = false;
200  short preValleyDepth = -1;
201  if (!m_checkPreValleyForXTalkId) preValleyExist = true;
202  else {
203  for (int iSample = jRawTime ; iSample - 1 > 0 ; iSample--) {
204  if (jRawTime - iSample > m_nSampleBefore) return false;
205  else if (wfm[iSample] - wfm[iSample - 1] >= 0) continue;
206  else {
207  preValleyDepth = (-1) * wfm[iSample];
208  if (preValleyDepth < m_preValleyDepthLoose) return false;
209  else {
210  if (!m_checkPostValleyForXTalkId) return true;
211  preValleyExist = true;
212  break;
213  }
214  }
215  }//for( iSample )
216  }//if( m_checkPreValleyForXTalkId )
217  if (!preValleyExist) return false;
218 
219  //check ringing (oscillation pattern) in trailing edge
220  short sign = 1;
221  short valley_adc = 9999;
222  for (int jSample = jRawTime ; jSample < nWfmSampling - 1 ; jSample++) {
223  if (jSample - jRawTime > m_nSampleAfter) return false;
224  if ((wfm[jSample + 1] - wfm[jSample])*sign > 0) continue;
225  else { //when peak or valley is found
226  if (sign < 0 && valley_adc > height + 1) //in case a valley is found
227  valley_adc = wfm[jSample];
228  if (sign > 0 && valley_adc < height) {//in case of second peak
229  if (wfm[jSample] - valley_adc > (height * m_2ndPeakAmplitudeRatioTight)
230  || (preValleyDepth > m_preValleyDepthTight
231  && (wfm[jSample] - valley_adc) > m_2ndPeakAmplitudeLoose))
232  return true;
233  else return false;
234  }
235  sign = (sign > 0 ? -1 : 1);
236  }
237  }//for( jSample0 )
238  return false;
239  }//if( jRawTime )
240  }//for( iWin )
241 
242  return false;
243 }
244 
245 
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::RelationsInterface::getRelated
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
Definition: RelationsObject.h:280
Belle2::TOPRawWaveform
Class to store raw data waveforms.
Definition: TOPRawWaveform.h:34
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::TOPXTalkChargeShareSetterModule
Crosstalk & chargeshare flag setter.
Definition: TOPXTalkChargeShareSetterModule.h:39
Belle2::TOPRawDigit
Class to store unpacked raw data (hits in feature-extraction format) It provides also calculation of ...
Definition: TOPRawDigit.h:34
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33