Belle II Software  release-05-01-25
PXDBgTupleProducerModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Benjamin Schwenker *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 
12 
13 #include <pxd/modules/pxdBackground/PXDBgTupleProducerModule.h>
14 #include <framework/datastore/StoreObjPtr.h>
15 #include <framework/datastore/StoreArray.h>
16 #include <framework/dataobjects/EventMetaData.h>
17 #include <pxd/reconstruction/PXDGainCalibrator.h>
18 #include <pxd/reconstruction/PXDPixelMasker.h>
19 
20 #include <pxd/dataobjects/PXDDigit.h>
21 #include <pxd/dataobjects/PXDCluster.h>
22 #include <boost/format.hpp>
23 
24 #include <TFile.h>
25 #include <TTree.h>
26 
27 using namespace std;
28 using boost::format;
29 using namespace Belle2;
30 using namespace Belle2::PXD;
31 
32 
33 
34 //-----------------------------------------------------------------
35 // Register the Module
36 //-----------------------------------------------------------------
37 REG_MODULE(PXDBgTupleProducer)
38 
39 
40 //-----------------------------------------------------------------
41 // Implementation
42 //-----------------------------------------------------------------
43 
45  , m_nPXDSensors(0), m_hasPXDData(false)
46 {
47  //Set module properties
48  setDescription("PXD background tuple producer module");
49  addParam("integrationTime", m_integrationTime, "PXD integration time in micro seconds", double(20));
50  addParam("timePeriod", m_timePeriod, "Period for background time series in seconds.", double(1));
51  addParam("outputFileName", m_outputFileName, "Output file name", string("beast_tuple.root"));
52  addParam("maskDeadPixels", m_maskDeadPixels, "Correct bg rates by known dead pixels", bool(true));
53  addParam("nBinsU", m_nBinsU, "Number of regions per sensor along u side", int(1));
54  addParam("nBinsV", m_nBinsV, "Number of regions per sensor along v side", int(6));
55 }
56 
57 
58 void PXDBgTupleProducerModule::initialize()
59 {
60  //Register collections
61  StoreArray<PXDCluster> storeClusters(m_storeClustersName);
62  StoreArray<PXDDigit> storeDigits(m_storeDigitsName);
63 
64  //Make sure the EventMetaData already exists.
65  StoreObjPtr<EventMetaData>().isRequired();
66 
67  //Store names to speed up creation later
68  m_storeDigitsName = storeDigits.getName();
69 
70  // PXD integration time
71  m_integrationTime *= Unit::us;
72 
73  // Period for time series
74  m_timePeriod *= Unit::s;
75 
76  // So far, we did not see PXD data
77  m_hasPXDData = false;
78 
79  //Pointer to GeoTools instance
80  auto gTools = VXD::GeoCache::getInstance().getGeoTools();
81  if (gTools->getNumberOfPXDLayers() == 0) {
82  B2WARNING("Missing geometry for PXD, PXD-masking is skiped.");
83  }
84  m_nPXDSensors = gTools->getNumberOfPXDSensors();
85 
86  // Initialize m_sensorData with empty sensorData for all sensors
87  for (int i = 0; i < m_nPXDSensors; i++) {
88  VxdID sensorID = gTools->getSensorIDFromPXDIndex(i);
89  m_sensorData[sensorID] = SensorData();
90  // Start value for minOccupancy should be one not zero
91  m_sensorData[sensorID].m_minOccupancy = 1.0;
92  // Initialize counters for subdivisions per sensor
93  m_sensorData[sensorID].m_regionExpoMap = vector<double>(m_nBinsU * m_nBinsV, 0);
94  m_sensorData[sensorID].m_regionDoseMap = vector<double>(m_nBinsU * m_nBinsV, 0);
95  m_sensorData[sensorID].m_regionSoftPhotonFluxMap = vector<double>(m_nBinsU * m_nBinsV, 0);
96  m_sensorData[sensorID].m_regionChargedParticleFluxMap = vector<double>(m_nBinsU * m_nBinsV, 0);
97  m_sensorData[sensorID].m_regionHardPhotonFluxMap = vector<double>(m_nBinsU * m_nBinsV, 0);
98  }
99 }
100 
101 void PXDBgTupleProducerModule::beginRun()
102 {
103  // Compute the sensitive area for all PXD sensors
104  for (auto const& pair2 : m_sensorData) {
105  auto const& sensorID = pair2.first;
106  auto info = getInfo(sensorID);
107 
108  // Compute nominal number of pixel per sensor
109  m_sensitivePixelMap[sensorID] = info.getUCells() * info.getVCells();
110  // Compute nominal area per sensor
111  m_sensitiveAreaMap[sensorID] = getSensorArea(sensorID);
112 
113  for (int uBin = 0; uBin < m_nBinsU; ++uBin) {
114  for (int vBin = 0; vBin < m_nBinsV; ++vBin) {
115  std::pair<VxdID, int> key(sensorID, getRegionID(uBin, vBin));
116  // Compute nominal number of pixel per sensor subregion
117  m_regionSensitivePixelMap[key] = info.getUCells() * info.getVCells() / m_nBinsU / m_nBinsV;
118  // Compute nominal area per sensor subregion
119  m_regionSensitiveAreaMap[key] = getRegionArea(sensorID, vBin);
120  }
121  }
122 
123  if (m_maskDeadPixels) {
124  for (int ui = 0; ui < info.getUCells(); ++ui) {
125  for (int vi = 0; vi < info.getVCells(); ++vi) {
126  if (PXD::PXDPixelMasker::getInstance().pixelDead(sensorID, ui, vi)
127  || !PXD::PXDPixelMasker::getInstance().pixelOK(sensorID, ui, vi)) {
128  m_sensitivePixelMap[sensorID] -= 1;
129  m_sensitiveAreaMap[sensorID] -= info.getVPitch(info.getVCellPosition(vi)) * info.getUPitch();
130  int uBin = PXDGainCalibrator::getInstance().getBinU(sensorID, ui, vi, m_nBinsU);
131  int vBin = PXDGainCalibrator::getInstance().getBinV(sensorID, vi, m_nBinsV);
132  std::pair<VxdID, int> key(sensorID, getRegionID(uBin, vBin));
133  m_regionSensitivePixelMap[key] -= 1;
134  m_regionSensitiveAreaMap[key] -= info.getVPitch(info.getVCellPosition(vi)) * info.getUPitch();
135  }
136  }
137  }
138  }
139 
140  if (m_sensitivePixelMap[sensorID] == 0) {
141  B2WARNING("All pixels from Sensor=" << sensorID << " are masked.");
142  }
143 
144  for (int uBin = 0; uBin < m_nBinsU; ++uBin) {
145  for (int vBin = 0; vBin < m_nBinsV; ++vBin) {
146  std::pair<VxdID, int> key(sensorID, getRegionID(uBin, vBin));
147  if (m_regionSensitivePixelMap[key] == 0) {
148  B2WARNING("All pixels from subregion uBin=" << uBin << " vBin=" << vBin << " on Sensor=" << sensorID << " are masked.");
149  }
150  }
151  }
152 
153  }
154 }
155 
156 void PXDBgTupleProducerModule::event()
157 {
158  //Register collections
159  const StoreArray<PXDCluster> storeClusters(m_storeClustersName);
160  const StoreArray<PXDDigit> storeDigits(m_storeDigitsName);
161 
162  //Get the event meta data
163  StoreObjPtr<EventMetaData> eventMetaDataPtr;
164 
165  // Compute the curent one second timestamp
166  unsigned long long int ts = eventMetaDataPtr->getTime() / m_timePeriod;
167 
168  // If needed, add a new one second block to buffer
169  auto iter = m_buffer.find(ts);
170  if (iter == m_buffer.end()) {
171  m_buffer[ts] = m_sensorData;
172  }
173 
174  // Empty map for computing event wise occupancy
175  std::map<VxdID, double> occupancyMap;
176  auto gTools = VXD::GeoCache::getInstance().getGeoTools();
177  for (int i = 0; i < m_nPXDSensors; i++) {
178  VxdID sensorID = gTools->getSensorIDFromPXDIndex(i);
179  occupancyMap[sensorID] = 0.0;
180  }
181 
182  // Check if there is PXD data
183  if (storeDigits.getEntries() > 0) {
184  m_hasPXDData = true;
185  }
186 
187  for (const PXDDigit& storeDigit : storeDigits) {
188  VxdID sensorID = storeDigit.getSensorID();
189  double ADUToEnergy = PXDGainCalibrator::getInstance().getADUToEnergy(sensorID, storeDigit.getUCellID(), storeDigit.getVCellID());
190  double hitEnergy = storeDigit.getCharge() * ADUToEnergy;
191 
192  if (m_sensitivePixelMap[sensorID] != 0) {
193  occupancyMap[sensorID] += 1.0 / m_sensitivePixelMap[sensorID];
194  }
195  m_buffer[ts][sensorID].m_dose += (hitEnergy / Unit::J);
196  m_buffer[ts][sensorID].m_expo += hitEnergy;
197 
198  int uBin = PXDGainCalibrator::getInstance().getBinU(sensorID, storeDigit.getUCellID(), storeDigit.getVCellID(), m_nBinsU);
199  int vBin = PXDGainCalibrator::getInstance().getBinV(sensorID, storeDigit.getVCellID(), m_nBinsV);
200  int regionID = getRegionID(uBin, vBin);
201  m_buffer[ts][sensorID].m_regionDoseMap[regionID] += (hitEnergy / Unit::J);
202  m_buffer[ts][sensorID].m_regionExpoMap[regionID] += hitEnergy;
203  }
204 
205  for (auto& pair : m_buffer[ts]) {
206  auto& sensorID = pair.first;
207  auto& bgdata = pair.second;
208  bgdata.m_run = eventMetaDataPtr->getRun();
209  // Check if there is actually data for this sensor
210  if (occupancyMap.find(sensorID) != occupancyMap.end()) {
211  bgdata.m_nEvents += 1;
212  bgdata.m_meanOccupancy += occupancyMap[sensorID];
213  if (occupancyMap[sensorID] > bgdata.m_maxOccupancy) {
214  bgdata.m_maxOccupancy = occupancyMap[sensorID];
215  }
216  if (occupancyMap[sensorID] < bgdata.m_minOccupancy) {
217  bgdata.m_minOccupancy = occupancyMap[sensorID];
218  }
219  }
220  }
221 
222  for (const PXDCluster& cluster : storeClusters) {
223  // Update if we have a new sensor
224  VxdID sensorID = cluster.getSensorID();
225  auto info = getInfo(sensorID);
226 
227  auto cluster_uID = info.getUCellID(cluster.getU());
228  auto cluster_vID = info.getVCellID(cluster.getV());
229  int uBin = PXDGainCalibrator::getInstance().getBinU(sensorID, cluster_uID, cluster_vID, m_nBinsU);
230  int vBin = PXDGainCalibrator::getInstance().getBinV(sensorID, cluster_vID, m_nBinsV);
231  int regionID = getRegionID(uBin, vBin);
232  double ADUToEnergy = PXDGainCalibrator::getInstance().getADUToEnergy(sensorID, cluster_uID, cluster_vID);
233  double clusterEnergy = cluster.getCharge() * ADUToEnergy;
234 
235  if (cluster.getSize() == 1 && clusterEnergy < 10000 * Unit::eV && clusterEnergy > 6000 * Unit::eV) {
236  m_buffer[ts][sensorID].m_softPhotonFlux += 1.0;
237  m_buffer[ts][sensorID].m_regionSoftPhotonFluxMap[regionID] += 1.0;
238  } else if (cluster.getSize() == 1 && clusterEnergy > 10000 * Unit::eV) {
239  m_buffer[ts][sensorID].m_hardPhotonFlux += 1.0;
240  m_buffer[ts][sensorID].m_regionHardPhotonFluxMap[regionID] += 1.0;
241  } else if (cluster.getSize() > 1 && clusterEnergy > 10000 * Unit::eV) {
242  m_buffer[ts][sensorID].m_chargedParticleFlux += 1.0;
243  m_buffer[ts][sensorID].m_regionChargedParticleFluxMap[regionID] += 1.0;
244  }
245  }
246 }
247 
248 void PXDBgTupleProducerModule::terminate()
249 {
250  // Create beast tuple
251  if (m_hasPXDData) {
252  TFile* rfile = new TFile(m_outputFileName.c_str(), "RECREATE");
253  TTree* treeBEAST = new TTree("tout", "BEAST data tree");
254 
255  unsigned int ts = 0;
256  treeBEAST->Branch("ts", &(ts));
257 
258  for (auto& pair : m_sensorData) {
259  auto& sensorID = pair.first;
260  auto& bgdata = pair.second;
261  string sensorDescr = sensorID;
262  treeBEAST->Branch(str(format("pxd_%1%_run") % sensorDescr).c_str(), &(bgdata.m_run));
263  treeBEAST->Branch(str(format("pxd_%1%_nEvents") % sensorDescr).c_str(), &(bgdata.m_nEvents));
264  treeBEAST->Branch(str(format("pxd_%1%_minOccupancy") % sensorDescr).c_str(), &(bgdata.m_minOccupancy));
265  treeBEAST->Branch(str(format("pxd_%1%_maxOccupancy") % sensorDescr).c_str(), &(bgdata.m_maxOccupancy));
266  treeBEAST->Branch(str(format("pxd_%1%_meanOccupancy") % sensorDescr).c_str(), &(bgdata.m_meanOccupancy));
267  treeBEAST->Branch(str(format("pxd_%1%_exposition") % sensorDescr).c_str(), &(bgdata.m_expo));
268  treeBEAST->Branch(str(format("pxd_%1%_dose") % sensorDescr).c_str(), &(bgdata.m_dose));
269  treeBEAST->Branch(str(format("pxd_%1%_softPhotonFlux") % sensorDescr).c_str(), &(bgdata.m_softPhotonFlux));
270  treeBEAST->Branch(str(format("pxd_%1%_hardPhotonFlux") % sensorDescr).c_str(), &(bgdata.m_hardPhotonFlux));
271  treeBEAST->Branch(str(format("pxd_%1%_chargedParticleFlux") % sensorDescr).c_str(),
272  &(bgdata.m_chargedParticleFlux));
273 
274  for (int uBin = 0; uBin < m_nBinsU; ++uBin) {
275  for (int vBin = 0; vBin < m_nBinsV; ++vBin) {
276  int regionID = getRegionID(uBin, vBin);
277  treeBEAST->Branch(str(format("pxd_%1%_region_%2%_%3%_exposition") % sensorDescr % uBin % vBin).c_str(),
278  &(bgdata.m_regionExpoMap[regionID]));
279  treeBEAST->Branch(str(format("pxd_%1%_region_%2%_%3%_dose") % sensorDescr % uBin % vBin).c_str(),
280  &(bgdata.m_regionDoseMap[regionID]));
281  treeBEAST->Branch(str(format("pxd_%1%_region_%2%_%3%_softPhotonFlux") % sensorDescr % uBin % vBin).c_str(),
282  &(bgdata.m_regionSoftPhotonFluxMap[regionID]));
283  treeBEAST->Branch(str(format("pxd_%1%_region_%2%_%3%_hardPhotonFlux") % sensorDescr % uBin % vBin).c_str(),
284  &(bgdata.m_regionHardPhotonFluxMap[regionID]));
285  treeBEAST->Branch(str(format("pxd_%1%_region_%2%_%3%_chargedParticleFlux") % sensorDescr % uBin % vBin).c_str(),
286  &(bgdata.m_regionChargedParticleFluxMap[regionID]));
287  }
288  }
289  }
290 
291  // Write timestamp and background rates into TTree
292  for (auto const& pair1 : m_buffer) {
293  auto const& timestamp = pair1.first;
294  auto const& sensors = pair1.second;
295 
296  // Set variables for dumping into tree
297  ts = timestamp;
298  for (auto const& pair2 : sensors) {
299  auto const& sensorID = pair2.first;
300  auto const& bgdata = pair2.second;
301  double currentComponentTime = bgdata.m_nEvents * m_integrationTime;
302  const PXD::SensorInfo& info = getInfo(sensorID);
303  double currentSensorMass = m_sensitiveAreaMap[sensorID] * info.getThickness() * c_densitySi;
304  double currentSensorArea = m_sensitiveAreaMap[sensorID];
305  m_sensorData[sensorID] = bgdata;
306  // Some bg rates are still in wrong units. We have to fix this now.
307  m_sensorData[sensorID].m_meanOccupancy = bgdata.m_meanOccupancy / bgdata.m_nEvents;
308 
309  if (currentSensorArea > 0) {
310  m_sensorData[sensorID].m_dose *= (1.0 / (currentComponentTime / Unit::s)) * (1000 / currentSensorMass);
311  m_sensorData[sensorID].m_expo *= (1.0 / currentSensorArea) * (1.0 / (currentComponentTime / Unit::s));
312  m_sensorData[sensorID].m_softPhotonFlux *= (1.0 / currentSensorArea) * (1.0 / (currentComponentTime / Unit::s));
313  m_sensorData[sensorID].m_hardPhotonFlux *= (1.0 / currentSensorArea) * (1.0 / (currentComponentTime / Unit::s));
314  m_sensorData[sensorID].m_chargedParticleFlux *= (1.0 / currentSensorArea) * (1.0 / (currentComponentTime / Unit::s));
315 
316  for (int regionID = 0; regionID < m_nBinsU * m_nBinsV; ++regionID) {
317  std::pair<VxdID, int> key(sensorID, regionID);
318  double currentRegionMass = m_regionSensitiveAreaMap[key] * info.getThickness() * c_densitySi;
319  double currentRegionArea = m_regionSensitiveAreaMap[key];
320  if (currentRegionArea > 0) {
321  m_sensorData[sensorID].m_regionDoseMap[regionID] *= (1.0 / currentComponentTime) * (1000 / currentRegionMass);
322  m_sensorData[sensorID].m_regionExpoMap[regionID] *= (1.0 / currentRegionArea) * (1.0 / (currentComponentTime / Unit::s));
323  m_sensorData[sensorID].m_regionSoftPhotonFluxMap[regionID] *= (1.0 / currentRegionArea) * (1.0 / (currentComponentTime / Unit::s));
324  m_sensorData[sensorID].m_regionHardPhotonFluxMap[regionID] *= (1.0 / currentRegionArea) * (1.0 / (currentComponentTime / Unit::s));
325  m_sensorData[sensorID].m_regionChargedParticleFluxMap[regionID] *= (1.0 / currentRegionArea) * (1.0 /
326  (currentComponentTime / Unit::s));
327  }
328  }
329  }
330  }
331  // Dump variables into tree
332  treeBEAST->Fill();
333  }
334 
335  // Write output tuple
336  rfile->cd();
337  treeBEAST->Write();
338  rfile->Close();
339  }
340 }
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::PXDDigit
The PXD digit class.
Definition: PXDDigit.h:38
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::PXD::SensorInfo
Specific implementation of SensorInfo for PXD Sensors which provides additional pixel specific inform...
Definition: SensorInfo.h:34
Belle2::PXD::PXDBgTupleProducerModule
PXD Background Tuple Producer.
Definition: PXDBgTupleProducerModule.h:50
Belle2::PXD
Namespace to encapsulate code needed for simulation and reconstrucion of the PXD.
Definition: PXDCalibrationUtilities.h:28
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::PXD::PXDBgTupleProducerModule::SensorData
Struct to hold data of an PXD sensor.
Definition: PXDBgTupleProducerModule.h:55
Belle2::PXDCluster
The PXD Cluster class This class stores all information about reconstructed PXD clusters The position...
Definition: PXDCluster.h:41
Belle2::StoreArray< PXDCluster >
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226