Belle II Software  release-05-02-19
SVDHitRateCounter.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Marko Staric, Hikaru Tanigawa, Ludovico Massaccesi *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 // Own include
12 #include <background/modules/BeamBkgHitRateMonitor/SVDHitRateCounter.h>
13 
14 // framework aux
15 #include <framework/logging/Logger.h>
16 #include <hlt/softwaretrigger/core/FinalTriggerDecisionCalculator.h>
17 #include <framework/gearbox/Const.h>
18 #include <framework/gearbox/Unit.h>
19 #include <vxd/geometry/GeoCache.h>
20 
21 using namespace std;
22 
23 #define LOGRATIO(x,y) (x) << " / " << (y) << " = " << ((x) * 100 / (y)) << "%"
24 
25 namespace Belle2 {
30  namespace Background {
31 
32  void SVDHitRateCounter::initialize(TTree* tree)
33  {
34  // register collection(s) as optional, your detector might be excluded in DAQ
35  m_digits.isOptional(m_svdShaperDigitsName);
36  m_clusters.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  int nSamp = m_eventInfo.isValid() ? m_eventInfo->getNSamples() : 6; // Assume 6 samples if real value unknown
245  if (nSamp < 2) nSamp = 2; // Avoid division by zero if nSamp = 1. Not sure it's correct, but I don't expect nSamp = 1 to happen.
246  double integrationTimeSeconds = (nSamp - 1) / c_SVDSamplingClockFrequency / Unit::s;
247  double chargePerUnitTime = cluster.getCharge() / integrationTimeSeconds;
248  if (cluster.isUCluster()) {
249  rates_energyU.layerAverageRates[layer] += chargePerUnitTime;
250  rates_energyU.layerLadderAverageRates[layer][ladder] += chargePerUnitTime;
251  rates_energyU.layerSensorAverageRates[layer][sensor] += chargePerUnitTime;
252  rates_energyU.averageRate += chargePerUnitTime;
253  if (layer == 0)
254  rates_energyU.l3LadderSensorAverageRates[ladder][sensor] += chargePerUnitTime;
255  clustersU.layerAverageRates[layer]++;
256  clustersU.layerLadderAverageRates[layer][ladder]++;
257  clustersU.layerSensorAverageRates[layer][sensor]++;
258  clustersU.averageRate++;
259  if (layer == 0)
260  clustersU.l3LadderSensorAverageRates[ladder][sensor]++;
261  } else {
262  rates_energyV.layerAverageRates[layer] += chargePerUnitTime;
263  rates_energyV.layerLadderAverageRates[layer][ladder] += chargePerUnitTime;
264  rates_energyV.layerSensorAverageRates[layer][sensor] += chargePerUnitTime;
265  rates_energyV.averageRate += chargePerUnitTime;
266  if (layer == 0)
267  rates_energyV.l3LadderSensorAverageRates[ladder][sensor] += chargePerUnitTime;
268  clustersV.layerAverageRates[layer]++;
269  clustersV.layerLadderAverageRates[layer][ladder]++;
270  clustersV.layerSensorAverageRates[layer][sensor]++;
271  clustersV.averageRate++;
272  if (layer == 0)
273  clustersV.l3LadderSensorAverageRates[ladder][sensor]++;
274  }
275  }
276 
277  // set flag to true to indicate the rates are valid
278  rates_highE.valid = true;
279  rates_lowE.valid = true;
280  clustersU.valid = true;
281  clustersV.valid = true;
282  rates_energyU.valid = true;
283  rates_energyV.valid = true;
284  }
285 
286  }
287 
288  void SVDHitRateCounter::normalize(unsigned timeStamp)
289  {
290  B2DEBUG(10, "SVDHitRateCounter: normalize()");
291  // copy buffer element
292  m_rates = m_buffer[timeStamp];
293  m_ratesU = m_bufferU[timeStamp];
294  m_ratesV = m_bufferV[timeStamp];
295  m_rates_highE = m_buffer_highE[timeStamp];
296  m_rates_lowE = m_buffer_lowE[timeStamp];
297  m_clustersU = m_buffer_clustersU[timeStamp];
298  m_clustersV = m_buffer_clustersV[timeStamp];
299  m_rates_energyU = m_buffer_energyU[timeStamp];
300  m_rates_energyV = m_buffer_energyV[timeStamp];
301 
302  SVDHitRateCounter::normalizeRates(m_rates);
303  SVDHitRateCounter::normalizeRates(m_ratesU, true);
304  SVDHitRateCounter::normalizeRates(m_ratesV, false, true);
305  SVDHitRateCounter::normalizeRates(m_rates_highE);
306  SVDHitRateCounter::normalizeRates(m_rates_lowE);
307  SVDHitRateCounter::normalizeRates(m_clustersU, true);
308  SVDHitRateCounter::normalizeRates(m_clustersV, false, true);
309  SVDHitRateCounter::normalizeEnergyRates(m_rates_energyU);
310  SVDHitRateCounter::normalizeEnergyRates(m_rates_energyV);
311  }
312 
313  void SVDHitRateCounter::normalizeRates(TreeStruct& rates, bool isU, bool isV)
314  {
315  if (not rates.valid) return;
316 
317  // normalize -> nHits on each segment in single event
318  rates.normalize();
319 
320  // Take the correct active strips counter
321  const auto& activeStrips = isU ? m_activeStripsU : (isV ? m_activeStripsV : m_activeStrips);
322  const auto& layerActiveStrips = isU ? m_layerActiveStripsU : (isV ? m_layerActiveStripsV : m_layerActiveStrips);
323  const auto& layerLadderActiveStrips = isU ? m_layerLadderActiveStripsU
324  : (isV ? m_layerLadderActiveStripsV : m_layerLadderActiveStrips);
325  const auto& layerSensorActiveStrips = isU ? m_layerSensorActiveStripsU
326  : (isV ? m_layerSensorActiveStripsV : m_layerSensorActiveStrips);
327  const auto& l3LadderSensorActiveStrips = isU ? m_l3LadderSensorActiveStripsU
328  : (isV ? m_l3LadderSensorActiveStripsV : m_l3LadderSensorActiveStrips);
329 
330  // convert to occupancy [%]
331  rates.averageRate /= activeStrips / 100.0;
332  for (int layer = 0; layer < m_nLayers; layer++) {
333  rates.layerAverageRates[layer] /= layerActiveStrips[layer] / 100.0;
334  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++) {
335  rates.layerLadderAverageRates[layer][ladder] /= layerLadderActiveStrips[layer][ladder] / 100.0;
336  }
337  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++) {
338  rates.layerSensorAverageRates[layer][sensor] /= layerSensorActiveStrips[layer][sensor] / 100.0;
339  }
340  }
341  int layer = 0;
342  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++) {
343  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++) {
344  rates.l3LadderSensorAverageRates[ladder][sensor] /= l3LadderSensorActiveStrips[ladder][sensor] / 100.0;
345  }
346  }
347  }
348 
349  void SVDHitRateCounter::normalizeEnergyRates(TreeStruct& rates)
350  {
351  static const double ehEnergyJoules = Const::ehEnergy / Unit::J;
352  // Convert charge/s to mrad/s by multiplying by energy/pair and dividing by the mass
353  static const double conv = ehEnergyJoules * 100e3;
354 
355  if (!rates.valid) return;
356 
357  rates.normalize(); // Divide by nEvents
358  // Convert to dose rate [mrad/s]
359  rates.averageRate *= conv / m_massKg;
360  for (int layer = 0; layer < m_nLayers; layer++) {
361  rates.layerAverageRates[layer] *= conv / m_layerMassKg[layer];
362  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++)
363  rates.layerLadderAverageRates[layer][ladder] *= conv / m_layerLadderMassKg[layer][ladder];
364  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++)
365  rates.layerSensorAverageRates[layer][sensor] *= conv / m_layerSensorMassKg[layer][sensor];
366  }
367  int layer = 0;
368  for (int ladder = 0; ladder < m_nLadders[layer]; ladder++)
369  for (int sensor = 0; sensor < m_nSensors[layer]; sensor++)
370  rates.l3LadderSensorAverageRates[ladder][sensor] *= conv / massOfSensor(layer, ladder, sensor);
371  }
372 
373  double SVDHitRateCounter::massOfSensor(int layer, int ladder, int sensor)
374  {
375  static const double rho_Si = 2.329 * Unit::g_cm3;
376  auto& sensorInfo = VXD::GeoCache::getInstance().getSensorInfo(VxdID(layer + 3, ladder + 1, sensor + 1));
377  double length = sensorInfo.getLength();
378  double width = (sensorInfo.getForwardWidth() + sensorInfo.getBackwardWidth()) / 2.0;
379  double thickness = sensorInfo.getWSize();
380  return length * width * thickness * rho_Si / 1e3;
381  }
382 
383  bool SVDHitRateCounter::isStripActive(const VxdID& sensorID, const bool& isU,
384  const unsigned short& strip)
385  {
386  return ((m_ignoreHotStripsPayload || !m_HotStripsCalib.isHot(sensorID, isU, strip))
387  && (m_ignoreMaskedStripsPayload || !m_FADCMaskedStrips.isMasked(sensorID, isU, strip)));
388  }
389 
390  } // background namespace
392 } // Belle2 namespace
393 
Belle2::PXD::rho_Si
const double rho_Si
Silicon density in g cm^-3.
Definition: PXDUtilities.h:42
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::Background::SVDHitRateCounter::TreeStruct
tree structure
Definition: SVDHitRateCounter.h:51
Belle2::Background::SVDHitRateCounter::TreeStruct::l3LadderSensorAverageRates
float l3LadderSensorAverageRates[7][2]
Layer 3 sensors [#ladder][#sensor].
Definition: SVDHitRateCounter.h:56
Belle2::VxdID::getLadderNumber
baseType getLadderNumber() const
Get the ladder id.
Definition: VxdID.h:108
Belle2::Background::SVDHitRateCounter::TreeStruct::normalize
void normalize()
normalize accumulated hits to single event
Definition: SVDHitRateCounter.h:63
Belle2::Background::SVDHitRateCounter::TreeStruct::averageRate
float averageRate
total SVD average occupancy
Definition: SVDHitRateCounter.h:55
Belle2::Background::SVDHitRateCounter::TreeStruct::layerLadderAverageRates
float layerLadderAverageRates[4][16]
[#layer][#ladder]
Definition: SVDHitRateCounter.h:53
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::Background::SVDHitRateCounter::TreeStruct::layerSensorAverageRates
float layerSensorAverageRates[4][5]
[#layer][#sensor]
Definition: SVDHitRateCounter.h:54
Belle2::VxdID::getSensorNumber
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:110
Belle2::Background::SVDHitRateCounter::TreeStruct::valid
bool valid
status: true = rates valid
Definition: SVDHitRateCounter.h:58
Belle2::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:106
Belle2::Background::SVDHitRateCounter::TreeStruct::layerAverageRates
float layerAverageRates[4]
layer average occupancy
Definition: SVDHitRateCounter.h:52