Belle II Software  release-08-01-10
SVDHitRateCounter.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 
9 // Own header.
10 #include <background/modules/BeamBkgHitRateMonitor/SVDHitRateCounter.h>
11 
12 // framework aux
13 #include <framework/logging/Logger.h>
14 #include <hlt/softwaretrigger/core/FinalTriggerDecisionCalculator.h>
15 #include <framework/gearbox/Const.h>
16 #include <framework/gearbox/Unit.h>
17 #include <vxd/geometry/GeoCache.h>
18 
19 using namespace std;
20 
21 #define LOGRATIO(x,y) (x) << " / " << (y) << " = " << ((x) * 100 / (y)) << "%"
22 
23 namespace Belle2 {
28  namespace Background {
29 
30  void SVDHitRateCounter::initialize(TTree* tree)
31  {
32  // register collection(s) as optional, your detector might be excluded in DAQ
33  m_digits.isOptional(m_svdShaperDigitsName);
34  m_clusters.isOptional();
35  m_resultStoreObjectPointer.isOptional();
36  m_eventInfo.isOptional();
37 
38  B2DEBUG(10, "SVDHitRateCounter: initialize()");
39  // set branch address
40  tree->Branch("svd", &m_rates,
41  "layerAverageRates[4]/F:layerLadderAverageRates[4][16]/F:layerSensorAverageRates[4][5]:averageRate/F:l3LadderSensorAverageRates[7][2]/F:numEvents/I:valid/O");
42  tree->Branch("svdU", &m_ratesU,
43  "layerAverageRates[4]/F:layerLadderAverageRates[4][16]/F:layerSensorAverageRates[4][5]:averageRate/F:l3LadderSensorAverageRates[7][2]/F:numEvents/I:valid/O");
44  tree->Branch("svdV", &m_ratesV,
45  "layerAverageRates[4]/F:layerLadderAverageRates[4][16]/F:layerSensorAverageRates[4][5]:averageRate/F:l3LadderSensorAverageRates[7][2]/F:numEvents/I:valid/O");
46  tree->Branch("svd_highE", &m_rates_highE,
47  "layerAverageRates[4]/F:layerLadderAverageRates[4][16]/F:layerSensorAverageRates[4][5]:averageRate/F:l3LadderSensorAverageRates[7][2]/F:numEvents/I:valid/O");
48  tree->Branch("svd_lowE", &m_rates_lowE,
49  "layerAverageRates[4]/F:layerLadderAverageRates[4][16]/F:layerSensorAverageRates[4][5]:averageRate/F:l3LadderSensorAverageRates[7][2]/F:numEvents/I:valid/O");
50  tree->Branch("svd_clustersU", &m_clustersU,
51  "layerAverageRates[4]/F:layerLadderAverageRates[4][16]/F:layerSensorAverageRates[4][5]:averageRate/F:l3LadderSensorAverageRates[7][2]/F:numEvents/I:valid/O");
52  tree->Branch("svd_clustersV", &m_clustersV,
53  "layerAverageRates[4]/F:layerLadderAverageRates[4][16]/F:layerSensorAverageRates[4][5]:averageRate/F:l3LadderSensorAverageRates[7][2]/F:numEvents/I:valid/O");
54  tree->Branch("svd_energyU", &m_rates_energyU,
55  "layerAverageRates[4]/F:layerLadderAverageRates[4][16]/F:layerSensorAverageRates[4][5]:averageRate/F:l3LadderSensorAverageRates[7][2]/F:numEvents/I:valid/O");
56  tree->Branch("svd_energyV", &m_rates_energyV,
57  "layerAverageRates[4]/F:layerLadderAverageRates[4][16]/F:layerSensorAverageRates[4][5]:averageRate/F:l3LadderSensorAverageRates[7][2]/F:numEvents/I:valid/O");
58 
59  // count active strips
60  for (int layer = 0; layer < m_nLayers; layer++) {
61  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++) {
62  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++) {
63  VxdID sensorID(layer + 3, ladder + 1, sensor + 1);
64  for (bool isU : {true, false}) {
65  int nStrips = nStripsOnLayerSide(layer, isU);
66  for (int strip = 0; strip < nStrips; strip++) {
67  if (isStripActive(sensorID, isU, strip)) {
68  m_activeStrips ++;
69  m_layerActiveStrips[layer] ++;
70  m_layerLadderActiveStrips[layer][ladder] ++;
71  m_layerSensorActiveStrips[layer][sensor] ++;
72  if (isU) {
73  m_activeStripsU++;
74  m_layerActiveStripsU[layer]++;
75  m_layerLadderActiveStripsU[layer][ladder]++;
76  m_layerSensorActiveStripsU[layer][sensor]++;
77  } else {
78  m_activeStripsV++;
79  m_layerActiveStripsV[layer]++;
80  m_layerLadderActiveStripsV[layer][ladder]++;
81  m_layerSensorActiveStripsV[layer][sensor]++;
82  }
83  }
84  }
85  }
86  }
87  }
88  }
89  int layer = 0;
90  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++) {
91  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++) {
92  VxdID sensorID(layer + 3, ladder + 1, sensor + 1);
93  for (bool isU : {true, false}) {
94  int nStrips = nStripsOnLayerSide(layer, isU);
95  for (int strip = 0; strip < nStrips; strip++) {
96  if (isStripActive(sensorID, isU, strip)) {
97  m_l3LadderSensorActiveStrips[ladder][sensor] ++;
98  if (isU)
99  m_l3LadderSensorActiveStripsU[ladder][sensor]++;
100  else
101  m_l3LadderSensorActiveStripsV[ladder][sensor]++;
102  }
103  }
104  }
105  }
106  }
107  B2INFO("SVD active strips = " << LOGRATIO(m_activeStrips, 223744));
108  for (layer = 0; layer < m_nLayers; layer++)
109  B2INFO(" Active strips L" << layer + 3 << ".X.X = "
110  << LOGRATIO(m_layerActiveStrips[layer], m_nLadders[layer] * m_nSensors[layer] * (nStripsOnLayerSide(layer, false) + 768)));
111  for (layer = 0; layer < m_nLayers; layer++)
112  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++)
113  B2INFO(" Active strips L" << layer + 3 << "." << ladder + 1 << ".X = "
114  << LOGRATIO(m_layerLadderActiveStrips[layer][ladder], m_nSensors[layer] * (nStripsOnLayerSide(layer, false) + 768)));
115  for (layer = 0; layer < m_nLayers; layer++)
116  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++)
117  B2INFO(" Active strips L" << layer + 3 << ".X." << sensor + 1 << " = "
118  << LOGRATIO(m_layerSensorActiveStrips[layer][sensor], m_nLadders[layer] * (nStripsOnLayerSide(layer, false) + 768)));
119  layer = 0;
120  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++)
121  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++)
122  B2INFO(" Active strips L3." << ladder + 1 << "." << sensor + 1 << " = "
123  << LOGRATIO(m_l3LadderSensorActiveStrips[ladder][sensor], 2 * 768));
124 
125  // Compute active mass
126  for (layer = 0; layer < m_nLayers; layer++) {
127  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++) {
128  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++) {
129  double mass = massOfSensor(layer, ladder, sensor);
130  m_massKg += mass;
131  m_layerMassKg[layer] += mass;
132  m_layerLadderMassKg[layer][ladder] += mass;
133  m_layerSensorMassKg[layer][sensor] += mass;
134  }
135  }
136  }
137  }
138 
139  void SVDHitRateCounter::clear()
140  {
141  m_buffer.clear();
142  }
143 
144  void SVDHitRateCounter::accumulate(unsigned timeStamp)
145  {
146  B2DEBUG(10, "SVDHitRateCounter: accumulate()");
147 
148  // check if the event has passed HLT filter
149  if (m_resultStoreObjectPointer.isValid()) {
150  const bool eventAccepted = SoftwareTrigger::FinalTriggerDecisionCalculator::getFinalTriggerDecision(*m_resultStoreObjectPointer);
151  if (!eventAccepted) return;
152  }
153 
154  // check if data are available
155  if (m_digits.isValid()) {
156 
157  // get buffer element
158  auto& rates = m_buffer[timeStamp];
159  auto& ratesU = m_bufferU[timeStamp];
160  auto& ratesV = m_bufferV[timeStamp];
161 
162  // increment event counter
163  rates.numEvents++;
164  ratesU.numEvents++;
165  ratesV.numEvents++;
166 
167  // accumulate hits
168  for (const auto& digit : m_digits) {
169  // select digits to count (usualy only good ones)
170  VxdID sensorID = digit.getSensorID();
171  int layer = sensorID.getLayerNumber() - 3;
172  int ladder = sensorID.getLadderNumber() - 1;
173  int sensor = sensorID.getSensorNumber() - 1;
174  rates.layerAverageRates[layer] ++;
175  rates.layerLadderAverageRates[layer][ladder] ++;
176  rates.layerSensorAverageRates[layer][sensor] ++;
177  rates.averageRate ++;
178  if (layer == 0)
179  rates.l3LadderSensorAverageRates[ladder][sensor] ++;
180 
181  if (digit.isUStrip()) {
182  ratesU.layerAverageRates[layer]++;
183  ratesU.layerLadderAverageRates[layer][ladder]++;
184  ratesU.layerSensorAverageRates[layer][sensor]++;
185  ratesU.averageRate++;
186  if (layer == 0)
187  ratesU.l3LadderSensorAverageRates[ladder][sensor]++;
188  } else {
189  ratesV.layerAverageRates[layer]++;
190  ratesV.layerLadderAverageRates[layer][ladder]++;
191  ratesV.layerSensorAverageRates[layer][sensor]++;
192  ratesV.averageRate++;
193  if (layer == 0)
194  ratesV.l3LadderSensorAverageRates[ladder][sensor]++;
195  }
196  }
197 
198  // set flag to true to indicate the rates are valid
199  rates.valid = true;
200  ratesU.valid = true;
201  ratesV.valid = true;
202  }
203 
204  // check if data are available
205  if (m_clusters.isValid()) {
206 
207  // get buffer element
208  auto& rates_highE = m_buffer_highE[timeStamp];
209  auto& rates_lowE = m_buffer_lowE[timeStamp];
210  auto& clustersU = m_buffer_clustersU[timeStamp];
211  auto& clustersV = m_buffer_clustersV[timeStamp];
212  auto& rates_energyU = m_buffer_energyU[timeStamp];
213  auto& rates_energyV = m_buffer_energyV[timeStamp];
214 
215  // increment event counter
216  rates_highE.numEvents++;
217  rates_lowE.numEvents++;
218  clustersU.numEvents++;
219  clustersV.numEvents++;
220  rates_energyU.numEvents++;
221  rates_energyV.numEvents++;
222 
223  // accumulate clusters
224  for (const auto& cluster : m_clusters) {
225  VxdID sensorID = cluster.getSensorID();
226  int layer = sensorID.getLayerNumber() - 3;
227  int ladder = sensorID.getLadderNumber() - 1;
228  int sensor = sensorID.getSensorNumber() - 1;
229  if (cluster.getCharge() > m_thrCharge) {
230  rates_highE.layerAverageRates[layer] ++;
231  rates_highE.layerLadderAverageRates[layer][ladder] ++;
232  rates_highE.layerSensorAverageRates[layer][sensor] ++;
233  rates_highE.averageRate ++;
234  if (layer == 0)
235  rates_highE.l3LadderSensorAverageRates[ladder][sensor] ++;
236  } else {
237  rates_lowE.layerAverageRates[layer] ++;
238  rates_lowE.layerLadderAverageRates[layer][ladder] ++;
239  rates_lowE.layerSensorAverageRates[layer][sensor] ++;
240  rates_lowE.averageRate ++;
241  if (layer == 0)
242  rates_lowE.l3LadderSensorAverageRates[ladder][sensor] ++;
243  }
244 
245  // Compute integration time for obtaining the dose rate
246  int nSamp = 6; // Fallback value
247  double svdSamplingClock = c_SVDSamplingClockFrequency; // Fallback value
248  // Replace number of samples with real value if available
249  if (m_eventInfo.isValid())
250  nSamp = m_eventInfo->getNSamples();
251  else
252  B2WARNING("SVDEventInfo not available: assuming 6 samples were taken.");
253  // Avoid division by zero if nSamp = 1. Valid values are only 6 and 3,
254  // so this line is only meant to prevent a crash in case of corrupted data.
255  if (nSamp < 2) nSamp = 2;
256  // Replace sampling frequency with real value if available
257  if (m_clockSettings.isValid())
258  svdSamplingClock = m_clockSettings->getClockFrequency(Const::SVD, "sampling");
259  else
260  B2WARNING("HardwareClockSettings not available: using approximated SVD sampling clock frequency.");
261  double integrationTimeSeconds = (nSamp - 1) / svdSamplingClock / Unit::s;
262  double chargePerUnitTime = cluster.getCharge() / integrationTimeSeconds;
263  if (cluster.isUCluster()) {
264  rates_energyU.layerAverageRates[layer] += chargePerUnitTime;
265  rates_energyU.layerLadderAverageRates[layer][ladder] += chargePerUnitTime;
266  rates_energyU.layerSensorAverageRates[layer][sensor] += chargePerUnitTime;
267  rates_energyU.averageRate += chargePerUnitTime;
268  if (layer == 0)
269  rates_energyU.l3LadderSensorAverageRates[ladder][sensor] += chargePerUnitTime;
270  clustersU.layerAverageRates[layer]++;
271  clustersU.layerLadderAverageRates[layer][ladder]++;
272  clustersU.layerSensorAverageRates[layer][sensor]++;
273  clustersU.averageRate++;
274  if (layer == 0)
275  clustersU.l3LadderSensorAverageRates[ladder][sensor]++;
276  } else {
277  rates_energyV.layerAverageRates[layer] += chargePerUnitTime;
278  rates_energyV.layerLadderAverageRates[layer][ladder] += chargePerUnitTime;
279  rates_energyV.layerSensorAverageRates[layer][sensor] += chargePerUnitTime;
280  rates_energyV.averageRate += chargePerUnitTime;
281  if (layer == 0)
282  rates_energyV.l3LadderSensorAverageRates[ladder][sensor] += chargePerUnitTime;
283  clustersV.layerAverageRates[layer]++;
284  clustersV.layerLadderAverageRates[layer][ladder]++;
285  clustersV.layerSensorAverageRates[layer][sensor]++;
286  clustersV.averageRate++;
287  if (layer == 0)
288  clustersV.l3LadderSensorAverageRates[ladder][sensor]++;
289  }
290  }
291 
292  // set flag to true to indicate the rates are valid
293  rates_highE.valid = true;
294  rates_lowE.valid = true;
295  clustersU.valid = true;
296  clustersV.valid = true;
297  rates_energyU.valid = true;
298  rates_energyV.valid = true;
299  }
300 
301  }
302 
303  void SVDHitRateCounter::normalize(unsigned timeStamp)
304  {
305  B2DEBUG(10, "SVDHitRateCounter: normalize()");
306  // copy buffer element
307  m_rates = m_buffer[timeStamp];
308  m_ratesU = m_bufferU[timeStamp];
309  m_ratesV = m_bufferV[timeStamp];
310  m_rates_highE = m_buffer_highE[timeStamp];
311  m_rates_lowE = m_buffer_lowE[timeStamp];
312  m_clustersU = m_buffer_clustersU[timeStamp];
313  m_clustersV = m_buffer_clustersV[timeStamp];
314  m_rates_energyU = m_buffer_energyU[timeStamp];
315  m_rates_energyV = m_buffer_energyV[timeStamp];
316 
317  SVDHitRateCounter::normalizeRates(m_rates);
318  SVDHitRateCounter::normalizeRates(m_ratesU, true);
319  SVDHitRateCounter::normalizeRates(m_ratesV, false, true);
320  SVDHitRateCounter::normalizeRates(m_rates_highE);
321  SVDHitRateCounter::normalizeRates(m_rates_lowE);
322  SVDHitRateCounter::normalizeRates(m_clustersU, true);
323  SVDHitRateCounter::normalizeRates(m_clustersV, false, true);
324  SVDHitRateCounter::normalizeEnergyRates(m_rates_energyU);
325  SVDHitRateCounter::normalizeEnergyRates(m_rates_energyV);
326  }
327 
328  void SVDHitRateCounter::normalizeRates(TreeStruct& rates, bool isU, bool isV)
329  {
330  if (not rates.valid) return;
331 
332  // normalize -> nHits on each segment in single event
333  rates.normalize();
334 
335  // Take the correct active strips counter
336  const auto& activeStrips = isU ? m_activeStripsU : (isV ? m_activeStripsV : m_activeStrips);
337  const auto& layerActiveStrips = isU ? m_layerActiveStripsU : (isV ? m_layerActiveStripsV : m_layerActiveStrips);
338  const auto& layerLadderActiveStrips = isU ? m_layerLadderActiveStripsU
339  : (isV ? m_layerLadderActiveStripsV : m_layerLadderActiveStrips);
340  const auto& layerSensorActiveStrips = isU ? m_layerSensorActiveStripsU
341  : (isV ? m_layerSensorActiveStripsV : m_layerSensorActiveStrips);
342  const auto& l3LadderSensorActiveStrips = isU ? m_l3LadderSensorActiveStripsU
343  : (isV ? m_l3LadderSensorActiveStripsV : m_l3LadderSensorActiveStrips);
344 
345  // convert to occupancy [%]
346  rates.averageRate /= activeStrips / 100.0;
347  for (int layer = 0; layer < m_nLayers; layer++) {
348  rates.layerAverageRates[layer] /= layerActiveStrips[layer] / 100.0;
349  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++) {
350  rates.layerLadderAverageRates[layer][ladder] /= layerLadderActiveStrips[layer][ladder] / 100.0;
351  }
352  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++) {
353  rates.layerSensorAverageRates[layer][sensor] /= layerSensorActiveStrips[layer][sensor] / 100.0;
354  }
355  }
356  int layer = 0;
357  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++) {
358  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++) {
359  rates.l3LadderSensorAverageRates[ladder][sensor] /= l3LadderSensorActiveStrips[ladder][sensor] / 100.0;
360  }
361  }
362  }
363 
364  void SVDHitRateCounter::normalizeEnergyRates(TreeStruct& rates)
365  {
366  static const double ehEnergyJoules = Const::ehEnergy / Unit::J;
367  // Convert charge/s to mrad/s by multiplying by energy/pair and dividing by the mass
368  static const double conv = ehEnergyJoules * 100e3;
369 
370  if (!rates.valid) return;
371 
372  rates.normalize(); // Divide by nEvents
373  // Convert to dose rate [mrad/s]
374  rates.averageRate *= conv / m_massKg;
375  for (int layer = 0; layer < m_nLayers; layer++) {
376  rates.layerAverageRates[layer] *= conv / m_layerMassKg[layer];
377  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++)
378  rates.layerLadderAverageRates[layer][ladder] *= conv / m_layerLadderMassKg[layer][ladder];
379  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++)
380  rates.layerSensorAverageRates[layer][sensor] *= conv / m_layerSensorMassKg[layer][sensor];
381  }
382  int layer = 0;
383  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++)
384  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++)
385  rates.l3LadderSensorAverageRates[ladder][sensor] *= conv / massOfSensor(layer, ladder, sensor);
386  }
387 
388  double SVDHitRateCounter::massOfSensor(int layer, int ladder, int sensor)
389  {
390  static const double rho_Si = 2.329 * Unit::g_cm3;
391  auto& sensorInfo = VXD::GeoCache::getInstance().getSensorInfo(VxdID(layer + 3, ladder + 1, sensor + 1));
392  double length = sensorInfo.getLength();
393  double width = (sensorInfo.getForwardWidth() + sensorInfo.getBackwardWidth()) / 2.0;
394  double thickness = sensorInfo.getWSize();
395  return length * width * thickness * rho_Si / 1e3;
396  }
397 
398  bool SVDHitRateCounter::isStripActive(const VxdID& sensorID, const bool& isU,
399  const unsigned short& strip)
400  {
401  return ((m_ignoreHotStripsPayload || !m_HotStripsCalib.isHot(sensorID, isU, strip))
402  && (m_ignoreMaskedStripsPayload || !m_FADCMaskedStrips.isMasked(sensorID, isU, strip)));
403  }
404 
405  } // background namespace
407 } // Belle2 namespace
408 
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:33
Belle2::VxdID::getSensorNumber
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:100
Belle2::VxdID::getLadderNumber
baseType getLadderNumber() const
Get the ladder id.
Definition: VxdID.h:98
Belle2::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:96
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17
Belle2::Background::SVDHitRateCounter::TreeStruct::layerAverageRates
float layerAverageRates[4]
layer average occupancy
Definition: SVDHitRateCounter.h:44
Belle2::Background::SVDHitRateCounter::TreeStruct::layerLadderAverageRates
float layerLadderAverageRates[4][16]
[#layer][#ladder]
Definition: SVDHitRateCounter.h:45
Belle2::Background::SVDHitRateCounter::TreeStruct::layerSensorAverageRates
float layerSensorAverageRates[4][5]
[#layer][#sensor]
Definition: SVDHitRateCounter.h:46
Belle2::Background::SVDHitRateCounter::TreeStruct::averageRate
float averageRate
total SVD average occupancy
Definition: SVDHitRateCounter.h:47
Belle2::Background::SVDHitRateCounter::TreeStruct::normalize
void normalize()
normalize accumulated hits to single event
Definition: SVDHitRateCounter.h:55
Belle2::Background::SVDHitRateCounter::TreeStruct::l3LadderSensorAverageRates
float l3LadderSensorAverageRates[7][2]
Layer 3 sensors [#ladder][#sensor].
Definition: SVDHitRateCounter.h:48