Belle II Software  release-05-01-25
PXDMCBgTupleProducerModule.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/PXDMCBgTupleProducerModule.h>
14 #include <framework/datastore/DataStore.h>
15 #include <framework/datastore/StoreObjPtr.h>
16 #include <framework/datastore/StoreArray.h>
17 #include <framework/dataobjects/BackgroundMetaData.h>
18 #include <pxd/reconstruction/PXDGainCalibrator.h>
19 #include <pxd/reconstruction/PXDPixelMasker.h>
20 
21 #include <pxd/dataobjects/PXDDigit.h>
22 #include <pxd/dataobjects/PXDCluster.h>
23 #include <boost/format.hpp>
24 
25 #include <TFile.h>
26 #include <TTree.h>
27 
28 using namespace std;
29 using boost::format;
30 using namespace Belle2;
31 using namespace Belle2::PXD;
32 
33 
34 
35 //-----------------------------------------------------------------
36 // Register the Module
37 //-----------------------------------------------------------------
38 REG_MODULE(PXDMCBgTupleProducer)
39 
40 
41 //-----------------------------------------------------------------
42 // Implementation
43 //-----------------------------------------------------------------
44 
45 PXDMCBgTupleProducerModule::PXDMCBgTupleProducerModule() : Module()
46  , m_hasPXDData(false), m_componentTime(0)
47 {
48  //Set module properties
49  setDescription("PXD background tuple producer module");
50  addParam("integrationTime", m_integrationTime, "PXD integration time in micro seconds", double(20));
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  addParam("overrideComponentTime", m_overrideComponentTime, "User specified component time in micro seconds", double(0.0));
56 }
57 
58 
59 
60 void PXDMCBgTupleProducerModule::initialize()
61 {
62  //Register collections
63  StoreArray<PXDCluster> storeClusters(m_storeClustersName);
64  StoreArray<PXDDigit> storeDigits(m_storeDigitsName);
65  StoreObjPtr<BackgroundMetaData> storeBgMetaData(m_storeBgMetaDataName, DataStore::c_Persistent);
66 
67  //Store names to speed up creation later
68  m_storeDigitsName = storeDigits.getName();
69  m_storeBgMetaDataName = storeBgMetaData.getName();
70 
71  // PXD integration time
72  m_integrationTime *= Unit::us;
73 
74  m_overrideComponentTime *= Unit::us;
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  int nPXDSensors = gTools->getNumberOfPXDSensors();
85 
86  // Initialize m_sensorData with empty sensorData for all sensors
87  for (int i = 0; i < nPXDSensors; i++) {
88  VxdID sensorID = gTools->getSensorIDFromPXDIndex(i);
89  m_sensorData[sensorID] = SensorData();
90  // Initialize counters for subdivisions per sensor
91  m_sensorData[sensorID].m_regionExpoMap = vector<double>(m_nBinsU * m_nBinsV, 0);
92  m_sensorData[sensorID].m_regionDoseMap = vector<double>(m_nBinsU * m_nBinsV, 0);
93  m_sensorData[sensorID].m_regionSoftPhotonFluxMap = vector<double>(m_nBinsU * m_nBinsV, 0);
94  m_sensorData[sensorID].m_regionChargedParticleFluxMap = vector<double>(m_nBinsU * m_nBinsV, 0);
95  m_sensorData[sensorID].m_regionHardPhotonFluxMap = vector<double>(m_nBinsU * m_nBinsV, 0);
96  }
97 }
98 
99 void PXDMCBgTupleProducerModule::beginRun()
100 {
101  // Compute the sensitive area for all PXD sensors
102  for (auto const& pair2 : m_sensorData) {
103  auto const& sensorID = pair2.first;
104  auto info = getInfo(sensorID);
105 
106  // Compute nominal number of pixel per sensor
107  m_sensitivePixelMap[sensorID] = info.getUCells() * info.getVCells();
108  // Compute nominal area per sensor
109  m_sensitiveAreaMap[sensorID] = getSensorArea(sensorID);
110 
111  for (int uBin = 0; uBin < m_nBinsU; ++uBin) {
112  for (int vBin = 0; vBin < m_nBinsV; ++vBin) {
113  std::pair<VxdID, int> key(sensorID, getRegionID(uBin, vBin));
114  // Compute nominal number of pixel per sensor subregion
115  m_regionSensitivePixelMap[key] = info.getUCells() * info.getVCells() / m_nBinsU / m_nBinsV;
116  // Compute nominal area per sensor subregion
117  m_regionSensitiveAreaMap[key] = getRegionArea(sensorID, vBin);
118  }
119  }
120 
121  if (m_maskDeadPixels) {
122  for (int ui = 0; ui < info.getUCells(); ++ui) {
123  for (int vi = 0; vi < info.getVCells(); ++vi) {
124  if (PXD::PXDPixelMasker::getInstance().pixelDead(sensorID, ui, vi)
125  || !PXD::PXDPixelMasker::getInstance().pixelOK(sensorID, ui, vi)) {
126  m_sensitivePixelMap[sensorID] -= 1;
127  m_sensitiveAreaMap[sensorID] -= info.getVPitch(info.getVCellPosition(vi)) * info.getUPitch();
128  int uBin = PXDGainCalibrator::getInstance().getBinU(sensorID, ui, vi, m_nBinsU);
129  int vBin = PXDGainCalibrator::getInstance().getBinV(sensorID, vi, m_nBinsV);
130  std::pair<VxdID, int> key(sensorID, getRegionID(uBin, vBin));
131  m_regionSensitivePixelMap[key] -= 1;
132  m_regionSensitiveAreaMap[key] -= info.getVPitch(info.getVCellPosition(vi)) * info.getUPitch();
133  }
134  }
135  }
136  }
137 
138  if (m_sensitivePixelMap[sensorID] == 0) {
139  B2WARNING("All pixels from Sensor=" << sensorID << " are masked.");
140  }
141 
142  for (int uBin = 0; uBin < m_nBinsU; ++uBin) {
143  for (int vBin = 0; vBin < m_nBinsV; ++vBin) {
144  std::pair<VxdID, int> key(sensorID, getRegionID(uBin, vBin));
145  if (m_regionSensitivePixelMap[key] == 0) {
146  B2WARNING("All pixels from subregion uBin=" << uBin << " vBin=" << vBin << " on Sensor=" << sensorID << " are masked.");
147  }
148  }
149  }
150 
151  }
152 }
153 
154 void PXDMCBgTupleProducerModule::event()
155 {
156  //Register collections
157  const StoreArray<PXDCluster> storeClusters(m_storeClustersName);
158  const StoreArray<PXDDigit> storeDigits(m_storeDigitsName);
159  const StoreObjPtr<BackgroundMetaData> storeBgMetaData(m_storeBgMetaDataName, DataStore::c_Persistent);
160 
161  // Set the real time
162  m_componentTime = storeBgMetaData->getRealTime();
163 
164  // Empty map for computing event wise occupancy
165  std::map<VxdID, double> occupancyMap;
166 
167  // Check if there is PXD data
168  if (storeDigits.getEntries() > 0) {
169  m_hasPXDData = true;
170  }
171 
172  for (const PXDDigit& storeDigit : storeDigits) {
173  VxdID sensorID = storeDigit.getSensorID();
174  double ADUToEnergy = PXDGainCalibrator::getInstance().getADUToEnergy(sensorID, storeDigit.getUCellID(), storeDigit.getVCellID());
175  double hitEnergy = storeDigit.getCharge() * ADUToEnergy;
176 
177  if (m_sensitivePixelMap[sensorID] != 0) {
178  occupancyMap[sensorID] += 1.0 / m_sensitivePixelMap[sensorID];
179  }
180  m_sensorData[sensorID].m_dose += (hitEnergy / Unit::J);
181  m_sensorData[sensorID].m_expo += hitEnergy;
182 
183  int uBin = PXDGainCalibrator::getInstance().getBinU(sensorID, storeDigit.getUCellID(), storeDigit.getVCellID(), m_nBinsU);
184  int vBin = PXDGainCalibrator::getInstance().getBinV(sensorID, storeDigit.getVCellID(), m_nBinsV);
185  int regionID = getRegionID(uBin, vBin);
186  m_sensorData[sensorID].m_regionDoseMap[regionID] += (hitEnergy / Unit::J);
187  m_sensorData[sensorID].m_regionExpoMap[regionID] += hitEnergy;
188  }
189 
190  for (auto& pair : m_sensorData) {
191  auto& sensorID = pair.first;
192  auto& bgdata = pair.second;
193 
194  // Check if there is actually data for this sensor
195  if (occupancyMap.find(sensorID) != occupancyMap.end()) {
196  bgdata.m_meanOccupancy += occupancyMap[sensorID];
197  }
198  }
199 
200  for (const PXDCluster& cluster : storeClusters) {
201  // Update if we have a new sensor
202  VxdID sensorID = cluster.getSensorID();
203  auto info = getInfo(sensorID);
204 
205  auto cluster_uID = info.getUCellID(cluster.getU());
206  auto cluster_vID = info.getVCellID(cluster.getV());
207  int uBin = PXDGainCalibrator::getInstance().getBinU(sensorID, cluster_uID, cluster_vID, m_nBinsU);
208  int vBin = PXDGainCalibrator::getInstance().getBinV(sensorID, cluster_vID, m_nBinsV);
209  int regionID = getRegionID(uBin, vBin);
210  double ADUToEnergy = PXDGainCalibrator::getInstance().getADUToEnergy(sensorID, cluster_uID, cluster_vID);
211  double clusterEnergy = cluster.getCharge() * ADUToEnergy;
212 
213  if (cluster.getSize() == 1 && clusterEnergy < 10000 * Unit::eV && clusterEnergy > 6000 * Unit::eV) {
214  m_sensorData[sensorID].m_softPhotonFlux += 1.0;
215  m_sensorData[sensorID].m_regionSoftPhotonFluxMap[regionID] += 1.0;
216  } else if (cluster.getSize() == 1 && clusterEnergy > 10000 * Unit::eV) {
217  m_sensorData[sensorID].m_hardPhotonFlux += 1.0;
218  m_sensorData[sensorID].m_regionHardPhotonFluxMap[regionID] += 1.0;
219  } else if (cluster.getSize() > 1 && clusterEnergy > 10000 * Unit::eV) {
220  m_sensorData[sensorID].m_chargedParticleFlux += 1.0;
221  m_sensorData[sensorID].m_regionChargedParticleFluxMap[regionID] += 1.0;
222  }
223  }
224 }
225 
226 void PXDMCBgTupleProducerModule::terminate()
227 {
228  // Create beast tuple
229  if (m_hasPXDData) {
230  TFile* rfile = new TFile(m_outputFileName.c_str(), "RECREATE");
231  TTree* treeBEAST = new TTree("tout", "BEAST data tree");
232 
233  double currentComponentTime = m_componentTime;
234  if (m_overrideComponentTime > 0.0) currentComponentTime = m_overrideComponentTime;
235 
236  B2RESULT("Total real time is " << currentComponentTime / Unit::us << " microseconds.");
237  B2RESULT("This is equivalent to " << currentComponentTime / m_integrationTime << " random triggered events.");
238 
239  for (auto& pair : m_sensorData) {
240  auto& sensorID = pair.first;
241  auto& bgdata = pair.second;
242  const PXD::SensorInfo& info = getInfo(sensorID);
243  double currentSensorMass = m_sensitiveAreaMap[sensorID] * info.getThickness() * c_densitySi;
244  double currentSensorArea = m_sensitiveAreaMap[sensorID];
245 
246  // Finalize computation of rates
247  m_sensorData[sensorID].m_meanOccupancy = bgdata.m_meanOccupancy * (m_integrationTime / currentComponentTime);
248 
249  if (currentSensorArea > 0) {
250  m_sensorData[sensorID].m_dose *= (1.0 / (currentComponentTime / Unit::s)) * (1000 / currentSensorMass);
251  m_sensorData[sensorID].m_expo *= (1.0 / currentSensorArea) * (1.0 / (currentComponentTime / Unit::s));
252  m_sensorData[sensorID].m_softPhotonFlux *= (1.0 / currentSensorArea) * (1.0 / (currentComponentTime / Unit::s));
253  m_sensorData[sensorID].m_hardPhotonFlux *= (1.0 / currentSensorArea) * (1.0 / (currentComponentTime / Unit::s));
254  m_sensorData[sensorID].m_chargedParticleFlux *= (1.0 / currentSensorArea) * (1.0 / (currentComponentTime / Unit::s));
255 
256  for (int regionID = 0; regionID < m_nBinsU * m_nBinsV; ++regionID) {
257  std::pair<VxdID, int> key(sensorID, regionID);
258  double currentRegionMass = m_regionSensitiveAreaMap[key] * info.getThickness() * c_densitySi;
259  double currentRegionArea = m_regionSensitiveAreaMap[key];
260  if (currentRegionArea > 0) {
261  m_sensorData[sensorID].m_regionDoseMap[regionID] *= (1.0 / currentComponentTime) * (1000 / currentRegionMass);
262  m_sensorData[sensorID].m_regionExpoMap[regionID] *= (1.0 / currentRegionArea) * (1.0 / (currentComponentTime / Unit::s));
263  m_sensorData[sensorID].m_regionSoftPhotonFluxMap[regionID] *= (1.0 / currentRegionArea) * (1.0 / (currentComponentTime / Unit::s));
264  m_sensorData[sensorID].m_regionHardPhotonFluxMap[regionID] *= (1.0 / currentRegionArea) * (1.0 / (currentComponentTime / Unit::s));
265  m_sensorData[sensorID].m_regionChargedParticleFluxMap[regionID] *= (1.0 / currentRegionArea) * (1.0 /
266  (currentComponentTime / Unit::s));
267  }
268  }
269  }
270 
271  // Prepare output tree
272  string sensorDescr = sensorID;
273  treeBEAST->Branch(str(format("pxd_%1%_meanOccupancy") % sensorDescr).c_str(), &(bgdata.m_meanOccupancy));
274  treeBEAST->Branch(str(format("pxd_%1%_exposition") % sensorDescr).c_str(), &(bgdata.m_expo));
275  treeBEAST->Branch(str(format("pxd_%1%_dose") % sensorDescr).c_str(), &(bgdata.m_dose));
276  treeBEAST->Branch(str(format("pxd_%1%_softPhotonFlux") % sensorDescr).c_str(), &(bgdata.m_softPhotonFlux));
277  treeBEAST->Branch(str(format("pxd_%1%_hardPhotonFlux") % sensorDescr).c_str(), &(bgdata.m_hardPhotonFlux));
278  treeBEAST->Branch(str(format("pxd_%1%_chargedParticleFlux") % sensorDescr).c_str(),
279  &(bgdata.m_chargedParticleFlux));
280 
281  for (int uBin = 0; uBin < m_nBinsU; ++uBin) {
282  for (int vBin = 0; vBin < m_nBinsV; ++vBin) {
283  int regionID = getRegionID(uBin, vBin);
284  treeBEAST->Branch(str(format("pxd_%1%_region_%2%_%3%_exposition") % sensorDescr % uBin % vBin).c_str(),
285  &(bgdata.m_regionExpoMap[regionID]));
286  treeBEAST->Branch(str(format("pxd_%1%_region_%2%_%3%_dose") % sensorDescr % uBin % vBin).c_str(),
287  &(bgdata.m_regionDoseMap[regionID]));
288  treeBEAST->Branch(str(format("pxd_%1%_region_%2%_%3%_softPhotonFlux") % sensorDescr % uBin % vBin).c_str(),
289  &(bgdata.m_regionSoftPhotonFluxMap[regionID]));
290  treeBEAST->Branch(str(format("pxd_%1%_region_%2%_%3%_hardPhotonFlux") % sensorDescr % uBin % vBin).c_str(),
291  &(bgdata.m_regionHardPhotonFluxMap[regionID]));
292  treeBEAST->Branch(str(format("pxd_%1%_region_%2%_%3%_chargedParticleFlux") % sensorDescr % uBin % vBin).c_str(),
293  &(bgdata.m_regionChargedParticleFluxMap[regionID]));
294  }
295  }
296  }
297 
298  // Dump variables into tree
299  treeBEAST->Fill();
300  // Write output tuple
301  rfile->cd();
302  treeBEAST->Write();
303  rfile->Close();
304  }
305 }
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::PXD::PXDMCBgTupleProducerModule
PXD MC Background Tuple Producer.
Definition: PXDMCBgTupleProducerModule.h:50
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
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::PXDCluster
The PXD Cluster class This class stores all information about reconstructed PXD clusters The position...
Definition: PXDCluster.h:41
Belle2::StoreArray< PXDCluster >
Belle2::PXD::PXDMCBgTupleProducerModule::SensorData
Struct to hold data of an PXD sensor.
Definition: PXDMCBgTupleProducerModule.h:55
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226