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>
20 #include <pxd/dataobjects/PXDDigit.h>
21 #include <pxd/dataobjects/PXDCluster.h>
22 #include <boost/format.hpp>
45 , m_nPXDSensors(0), m_hasPXDData(false)
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));
58 void PXDBgTupleProducerModule::initialize()
68 m_storeDigitsName = storeDigits.getName();
71 m_integrationTime *= Unit::us;
74 m_timePeriod *= Unit::s;
80 auto gTools = VXD::GeoCache::getInstance().getGeoTools();
81 if (gTools->getNumberOfPXDLayers() == 0) {
82 B2WARNING(
"Missing geometry for PXD, PXD-masking is skiped.");
84 m_nPXDSensors = gTools->getNumberOfPXDSensors();
87 for (
int i = 0; i < m_nPXDSensors; i++) {
88 VxdID sensorID = gTools->getSensorIDFromPXDIndex(i);
91 m_sensorData[sensorID].m_minOccupancy = 1.0;
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);
101 void PXDBgTupleProducerModule::beginRun()
104 for (
auto const& pair2 : m_sensorData) {
105 auto const& sensorID = pair2.first;
106 auto info = getInfo(sensorID);
109 m_sensitivePixelMap[sensorID] = info.getUCells() * info.getVCells();
111 m_sensitiveAreaMap[sensorID] = getSensorArea(sensorID);
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));
117 m_regionSensitivePixelMap[key] = info.getUCells() * info.getVCells() / m_nBinsU / m_nBinsV;
119 m_regionSensitiveAreaMap[key] = getRegionArea(sensorID, vBin);
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();
140 if (m_sensitivePixelMap[sensorID] == 0) {
141 B2WARNING(
"All pixels from Sensor=" << sensorID <<
" are masked.");
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.");
156 void PXDBgTupleProducerModule::event()
166 unsigned long long int ts = eventMetaDataPtr->getTime() / m_timePeriod;
169 auto iter = m_buffer.find(ts);
170 if (iter == m_buffer.end()) {
171 m_buffer[ts] = m_sensorData;
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;
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;
192 if (m_sensitivePixelMap[sensorID] != 0) {
193 occupancyMap[sensorID] += 1.0 / m_sensitivePixelMap[sensorID];
195 m_buffer[ts][sensorID].m_dose += (hitEnergy / Unit::J);
196 m_buffer[ts][sensorID].m_expo += hitEnergy;
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;
205 for (
auto& pair : m_buffer[ts]) {
206 auto& sensorID = pair.first;
207 auto& bgdata = pair.second;
208 bgdata.m_run = eventMetaDataPtr->getRun();
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];
216 if (occupancyMap[sensorID] < bgdata.m_minOccupancy) {
217 bgdata.m_minOccupancy = occupancyMap[sensorID];
222 for (
const PXDCluster& cluster : storeClusters) {
224 VxdID sensorID = cluster.getSensorID();
225 auto info = getInfo(sensorID);
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;
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;
248 void PXDBgTupleProducerModule::terminate()
252 TFile* rfile =
new TFile(m_outputFileName.c_str(),
"RECREATE");
253 TTree* treeBEAST =
new TTree(
"tout",
"BEAST data tree");
256 treeBEAST->Branch(
"ts", &(ts));
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));
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]));
292 for (
auto const& pair1 : m_buffer) {
293 auto const& timestamp = pair1.first;
294 auto const& sensors = pair1.second;
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;
303 double currentSensorMass = m_sensitiveAreaMap[sensorID] * info.getThickness() * c_densitySi;
304 double currentSensorArea = m_sensitiveAreaMap[sensorID];
305 m_sensorData[sensorID] = bgdata;
307 m_sensorData[sensorID].m_meanOccupancy = bgdata.m_meanOccupancy / bgdata.m_nEvents;
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));
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));