Belle II Software development
PXDBackgroundModule.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#include <pxd/modules/pxdBackground/PXDBackgroundModule.h>
9
10#include <framework/datastore/DataStore.h>
11#include <framework/datastore/StoreObjPtr.h>
12#include <framework/datastore/StoreArray.h>
13#include <framework/datastore/RelationArray.h>
14
15#include <framework/dataobjects/FileMetaData.h>
16#include <framework/dataobjects/BackgroundMetaData.h>
17#include <framework/gearbox/Const.h>
18#include <mdst/dataobjects/MCParticle.h>
19#include <pxd/dataobjects/PXDSimHit.h>
20#include <pxd/dataobjects/PXDTrueHit.h>
21#include <pxd/dataobjects/PXDDigit.h>
22#include <pxd/dataobjects/PXDCluster.h>
23#include <pxd/dataobjects/PXDEnergyDepositionEvent.h>
24#include <pxd/dataobjects/PXDNeutronFluxEvent.h>
25#include <pxd/dataobjects/PXDOccupancyEvent.h>
26#include <cmath>
27#include <fstream>
28#include <boost/format.hpp>
29
30using namespace std;
31using boost::format;
32using namespace Belle2;
33using namespace Belle2::PXD;
34
35
36
37//-----------------------------------------------------------------
38// Register the Module
39//-----------------------------------------------------------------
40REG_MODULE(PXDBackground);
41
42
43//-----------------------------------------------------------------
44// Implementation
45//-----------------------------------------------------------------
46
57{
58 //Set module properties
59 setDescription("PXD background module");
60 setPropertyFlags(c_ParallelProcessingCertified); // specify this flag if you need parallel processing
61 // FIXME: This information can in principle be extracted from bg files, though not trivially.
62 addParam("componentName", m_componentName, "Background component name to process", m_componentName);
63 addParam("componentTime", m_componentTime, "Background component time", m_componentTime);
64 addParam("integrationTime", m_integrationTime, "PXD integration time", m_integrationTime);
65 addParam("doseReportingLevel", m_doseReportingLevel, "0 - no data, 1 - summary only, 2 - summary + ntuple", m_doseReportingLevel);
66 addParam("nfluxReportingLevel", m_nfluxReportingLevel, "0 - no data, 1 - summary only, 2 - summary + ntuple",
68 addParam("occupancyReportingLevel", m_occupancyReportingLevel, "0 - no data, 1 - summary only, 2 - summary + ntuple",
70 addParam("outputDirectory", m_outputDirectoryName, "Name of output directory", m_outputDirectoryName);
71}
72
73const ROOT::Math::XYZVector& PXDBackgroundModule::pointToGlobal(VxdID sensorID, const ROOT::Math::XYZVector& local)
74{
75 static ROOT::Math::XYZVector result(0, 0, 0);
76
77 const PXD::SensorInfo& info = getInfo(sensorID);
78 result = info.pointToGlobal(local);
79 return result;
80}
81
82const ROOT::Math::XYZVector& PXDBackgroundModule::vectorToGlobal(VxdID sensorID, const ROOT::Math::XYZVector& local)
83{
84 static ROOT::Math::XYZVector result(0, 0, 0);
85
86 const PXD::SensorInfo& info = getInfo(sensorID);
87 result = info.vectorToGlobal(local);
88 return result;
89}
90
94
96{
97 //Register collections
105
106 RelationArray relDigitsMCParticles(storeDigits, storeMCParticles);
107 RelationArray relDigitsTrueHits(storeDigits, storeTrueHits);
108 RelationArray relMCParticlesTrueHits(storeMCParticles, storeTrueHits);
109 RelationArray relTrueHitsSimHits(storeTrueHits, storeSimHits);
110
111 // Add two new StoreArrays
113 storeEnergyDeposits.registerInDataStore();
115 storeNeutronFluxes.registerInDataStore();
117 storeOccupancyEvents.registerInDataStore();
118
119 //Store names to speed up creation later
120 m_storeFileMetaDataName = storeFileMetaData.getName();
121 m_storeBgMetaDataName = storeBgMetaData.getName();
122 m_storeMCParticlesName = storeMCParticles.getName();
123 m_storeSimHitsName = storeSimHits.getName();
124 m_storeTrueHitsName = storeTrueHits.getName();
125 m_storeDigitsName = storeDigits.getName();
126 m_relDigitsMCParticlesName = relDigitsMCParticles.getName();
127 m_relDigitsTrueHitsName = relDigitsTrueHits.getName();
128 m_relParticlesTrueHitsName = relMCParticlesTrueHits.getName();
129 m_relTrueHitsSimHitsName = relTrueHitsSimHits.getName();
130 m_storeEnergyDepositsName = storeEnergyDeposits.getName();
131 m_storeNeutronFluxesName = storeNeutronFluxes.getName();
132
135}
136
140
142{
143 //Register collections
146 const StoreArray<MCParticle> storeMCParticles(m_storeMCParticlesName);
147 const StoreArray<PXDSimHit> storeSimHits(m_storeSimHitsName);
148 const StoreArray<PXDTrueHit> storeTrueHits(m_storeTrueHitsName);
149 const StoreArray<PXDDigit> storeDigits(m_storeDigitsName);
150 const StoreArray<PXDCluster> storeClsuters(m_storeClustersName);
151
152 // Add two new StoreArrays
156
157 // Relations
158 RelationArray relDigitsMCParticles(storeDigits, storeMCParticles, m_relDigitsMCParticlesName);
159 RelationArray relDigitsTrueHits(storeDigits, storeTrueHits, m_relDigitsTrueHitsName);
160 RelationArray relTrueHitsSimHits(storeTrueHits, storeSimHits, m_relTrueHitsSimHitsName);
161 RelationArray relTrueHitsMCParticles(storeMCParticles, storeTrueHits, m_relParticlesTrueHitsName);
162
163 // unsigned long numberOfEvents = storeFileMetaData->getNEvents();
164 double currentComponentTime = storeBgMetaData->getRealTime();
165 if (currentComponentTime != m_componentTime)
166 B2FATAL("Mismatch in component times:\n"
167 << "Steering file: " << m_componentTime << "\n"
168 << "Background file: " << currentComponentTime);
169
170 VxdID currentSensorID(0);
171 double currentSensorThickness(0);
172 double currentSensorArea(0);
173
174 // Exposition and dose
176 B2DEBUG(100, "Expo and dose");
177 currentSensorID.setID(0);
178 double currentSensorMass(0);
179 for (const PXDSimHit& hit : storeSimHits) {
180 // Update if we have a new sensor
181 VxdID sensorID = hit.getSensorID();
182 if (sensorID != currentSensorID) {
183 currentSensorID = sensorID;
184 currentSensorThickness = getSensorThickness(currentSensorID);
185 currentSensorMass = getSensorMass(currentSensorID);
186 currentSensorArea = getSensorArea(currentSensorID);
187 }
188 double hitEnergy = hit.getElectrons() * Const::ehEnergy;
189 // Dose in Gy/smy, normalize by sensor mass
190 m_sensorData[currentSensorID].m_dose +=
191 (hitEnergy / Unit::J) / (currentSensorMass / 1000) * (c_smy / currentComponentTime);
192 // Exposition in GeV/cm2/s
193 m_sensorData[currentSensorID].m_expo += hitEnergy / currentSensorArea / (currentComponentTime / Unit::s);
195 const ROOT::Math::XYZVector localPos = hit.getPosIn();
196 const ROOT::Math::XYZVector globalPos = pointToGlobal(currentSensorID, localPos);
197 float globalPosXYZ[3];
198 globalPos.GetCoordinates(globalPosXYZ);
199 storeEnergyDeposits.appendNew(
200 sensorID.getLayerNumber(), sensorID.getLadderNumber(), sensorID.getSensorNumber(),
201 hit.getPDGcode(), hit.getGlobalTime(),
202 localPos.X(), localPos.Y(), globalPosXYZ, hitEnergy,
203 (hitEnergy / Unit::J) / (currentSensorMass / 1000) / (currentComponentTime / Unit::s),
204 (hitEnergy / Unit::J) / currentSensorArea / (currentComponentTime / Unit::s)
205 );
206 }
207 }
208 }
209
210 // Neutron flux
212 B2DEBUG(100, "Neutron flux");
213 currentSensorID.setID(0);
214 for (const PXDTrueHit& hit : storeTrueHits) {
215 VxdID sensorID = hit.getSensorID();
216 // Update if we are on a new sensor
217 if (sensorID != currentSensorID) {
218 currentSensorID = sensorID;
219 currentSensorThickness = getSensorThickness(currentSensorID);
220 currentSensorArea = getSensorArea(currentSensorID);
221 }
222 // J(TrueHit) = abs(step)/thickness * correctionFactor;
223 ROOT::Math::XYZVector entryPos(hit.getEntryU(), hit.getEntryV(), hit.getEntryW());
224 ROOT::Math::XYZVector exitPos(hit.getExitU(), hit.getExitV(), hit.getExitW());
225 double stepLength = (exitPos - entryPos).R();
226 // Identify what particle we've got. We need type and kinetic energy.
227 // TODO: TrueHit must carry pdg or SimHit must carry energy.
228 // NOTE: MCParticles may get remapped, then SimHits still carry correct pdg.
229 const PXDSimHit* simhit = hit.getRelatedTo<PXDSimHit>();
230 if (!simhit) { //either something is very wrong, or we just don't have the relation. Try to find an appropriate SimHit manually.
231 double minDistance = 1.0e10;
232 for (const PXDSimHit& related : storeSimHits) {
233 double distance = (entryPos - related.getPosIn()).R();
234 if (distance < minDistance) {
235 minDistance = distance;
236 simhit = &related;
237 }
238 }
239 }
240 if (!simhit) {
241 B2WARNING("No related PXDSimHit found");
242 continue; //skip this true hit if the simhit is still null after setting it manually
243 }
244 // FIXME: Is there a difference between positrons and electrons wrt. NIEL?
245 // We fill neutronFluxBars with summary NIEL deposit for all kinds of particles by layer and component.
246 // Fluency plots are by component and are deposition histograms for a particular type of particle and compoonent.
247 // Special treatment of corrupt p's in TrueHits:
248 ROOT::Math::XYZVector hitMomentum(hit.getMomentum());
249 hitMomentum.SetX(std::isfinite(hitMomentum.X()) ? hitMomentum.X() : 0.0);
250 hitMomentum.SetY(std::isfinite(hitMomentum.Y()) ? hitMomentum.Y() : 0.0);
251 hitMomentum.SetZ(std::isfinite(hitMomentum.Z()) ? hitMomentum.Z() : 0.0);
252 int pdg = abs(simhit->getPDGcode());
253 double kineticEnergy(0.0);
254 double nielWeight(0.0);
255 if (pdg == Const::neutron.getPDGCode()) {
256 double m0 = Const::neutronMass;
257 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
258 nielWeight = m_nielNeutrons->getNielFactor(kineticEnergy / Unit::MeV);
259 }
260 if (pdg == Const::proton.getPDGCode()) {
261 double m0 = Const::protonMass;
262 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
263 nielWeight = m_nielProtons->getNielFactor(kineticEnergy / Unit::MeV);
264 }
265 if (pdg == Const::pi0.getPDGCode() || pdg == Const::pion.getPDGCode()) {
266 double m0 = Const::pi0Mass;
267 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
268 nielWeight = m_nielPions->getNielFactor(kineticEnergy / Unit::MeV);
269 }
270 if (pdg == Const::electron.getPDGCode()) {
271 double m0 = Const::electronMass;
272 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
273 nielWeight = m_nielElectrons->getNielFactor(kineticEnergy / Unit::MeV);
274 }
275 if (pdg == Const::photon.getPDGCode()) {
276 double m0 = 0.0;
277 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
278 }
279
280 // Only set weight for supported particles
281 nielWeight = std::isfinite(nielWeight) ? nielWeight : 0.0;
282 m_sensorData[currentSensorID].m_neutronFlux += nielWeight * stepLength / currentSensorThickness / currentSensorArea /
283 currentComponentTime * c_smy;
284
285 // Store data in a PXDNeutronFluxEvent object
287 ROOT::Math::XYZVector localPos(hit.getU(), hit.getV(), hit.getW());
288 const ROOT::Math::XYZVector globalPos = pointToGlobal(currentSensorID, localPos);
289 float globalPosXYZ[3];
290 globalPos.GetCoordinates(globalPosXYZ);
291 ROOT::Math::XYZVector localMom = hit.getMomentum();
292 const ROOT::Math::XYZVector globalMom = vectorToGlobal(currentSensorID, localMom);
293 float globalMomXYZ[3];
294 globalMom.GetCoordinates(globalMomXYZ);
295 storeNeutronFluxes.appendNew(
296 sensorID.getLayerNumber(), sensorID.getLadderNumber(), sensorID.getSensorNumber(),
297 simhit->getPDGcode(), simhit->getGlobalTime(),
298 hit.getU(), hit.getV(), globalPosXYZ, globalMomXYZ, kineticEnergy,
299 stepLength, nielWeight,
300 stepLength / currentSensorThickness / currentSensorArea / (currentComponentTime / Unit::s),
301 nielWeight * stepLength / currentSensorThickness / currentSensorArea / (currentComponentTime / Unit::s)
302 );
303 }
304 }
305 }
306
307 // Occupancy
309 B2DEBUG(100, "Fired pixels");
310 currentSensorID.setID(0);
311 double currentSensorCut = 0;
312 // Store fired pixels: count number of digits over threshold
313 std::map<VxdID, std::vector<float> > firedPixels;
314 for (const PXDDigit& storeDigit : storeDigits) {
315 // Filter out digits with signals below zero-suppression threshold
316 // ARE THERE SUCH DIGITS?
317 VxdID sensorID = storeDigit.getSensorID();
318 if (sensorID != currentSensorID) {
319 currentSensorID = sensorID;
320 auto info = getInfo(sensorID);
321 currentSensorCut = info.getChargeThreshold();
322 }
323 B2DEBUG(30, "Digit charge: " << storeDigit.getCharge() << " threshold: " << currentSensorCut);
324 if (storeDigit.getCharge() < currentSensorCut) continue;
325 B2DEBUG(30, "Passed.");
326 firedPixels[sensorID].push_back(storeDigit.getCharge());
327 }
328 // Process the map
329 for (auto idAndSet : firedPixels) {
330 VxdID sensorID = idAndSet.first;
331 double sensorArea = getSensorArea(sensorID);
332 int nFired = idAndSet.second.size();
333 double fired = nFired / (currentComponentTime / Unit::s) / sensorArea;
334 m_sensorData[sensorID].m_fired += fired;
335 }
336
337 B2DEBUG(100, "Occupancy");
338 currentSensorID.setID(0);
339 int nPixels = 0;
340 for (auto cluster : storeClsuters) {
341 VxdID sensorID = cluster.getSensorID();
342 if (currentSensorID != sensorID) {
343 currentSensorID = sensorID;
344 auto info = getInfo(sensorID);
345 nPixels = info.getUCells() * info.getVCells();
346 }
347
348 double w_acceptance = m_integrationTime / currentComponentTime;
349 double occupancy = 1.0 / nPixels * cluster.getSize();
350 m_sensorData[sensorID].m_occupancy += w_acceptance * occupancy;
351
353 storeOccupancyEvents.appendNew(
354 sensorID.getLayerNumber(), sensorID.getLadderNumber(),
355 sensorID.getSensorNumber(),
356 cluster.getU(), cluster.getV(), cluster.getSize(),
357 cluster.getCharge(), occupancy
358 );
359 }
360 }
361 }
362}
363
367
368
370{
371 // Write out m_data
372 ofstream outfile;
373 string outfileName(m_outputDirectoryName + m_componentName + "_summary.txt");
374 outfile.open(outfileName.c_str(), ios::out | ios::trunc);
375 outfile << "component_name\t"
376 << "component_time\t"
377 << "layer\t"
378 << "ladder\t"
379 << "sensor\t"
380 << "dose\t"
381 << "expo\t"
382 << "neutronFlux\t"
383 << "fired\t"
384 << "occupancy"
385 << endl;
386 double componentTime = m_componentTime / Unit::us;
387 for (auto vxdSensor : m_sensorData) {
388 outfile << m_componentName.c_str() << "\t"
389 << componentTime << "\t"
390 << vxdSensor.first.getLayerNumber() << "\t"
391 << vxdSensor.first.getLadderNumber() << "\t"
392 << vxdSensor.first.getSensorNumber() << "\t"
393 << vxdSensor.second.m_dose << "\t"
394 << vxdSensor.second.m_expo << "\t"
395 << vxdSensor.second.m_neutronFlux << "\t"
396 << vxdSensor.second.m_fired << "\t"
397 << vxdSensor.second.m_occupancy
398 << endl;
399 }
400 outfile << endl;
401}
double R
typedef autogenerated by FFTW
static const ParticleType neutron
neutron particle
Definition Const.h:675
static const ParticleType pi0
neutral pion particle
Definition Const.h:674
static const ChargedStable pion
charged pion particle
Definition Const.h:661
static const double electronMass
electron mass
Definition Const.h:685
static const double neutronMass
neutron mass
Definition Const.h:692
static const ChargedStable proton
proton particle
Definition Const.h:663
static const double ehEnergy
Energy needed to create an electron-hole pair in Si at std.
Definition Const.h:697
static const double protonMass
proton mass
Definition Const.h:689
static const ParticleType photon
photon particle
Definition Const.h:673
static const double pi0Mass
neutral pion mass
Definition Const.h:691
static const ChargedStable electron
electron particle
Definition Const.h:659
@ c_Persistent
Object is available during entire execution time.
Definition DataStore.h:60
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition Module.cc:208
Module()
Constructor.
Definition Module.cc:30
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition Module.h:80
The PXD digit class.
Definition PXDDigit.h:27
Class PXDSimHit - Geant4 simulated hit for the PXD.
Definition PXDSimHit.h:24
Class PXDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition PXDTrueHit.h:31
static const unsigned short c_reportNTuple
Summary and NTuple.
std::string m_storeFileMetaDataName
Name of the persistent FileMetaData object.
unsigned short m_nfluxReportingLevel
0 - no data, 1 - summary only, 2 - ntuple
const ROOT::Math::XYZVector & vectorToGlobal(VxdID sensorID, const ROOT::Math::XYZVector &local)
Convert local vector coordinates to global.
const ROOT::Math::XYZVector & pointToGlobal(VxdID sensorID, const ROOT::Math::XYZVector &local)
Convert local sensor coordinates to global.
std::string m_storeEnergyDepositsName
PXDEnergyDepositEvents StoreArray name.
const double c_smy
Seconds in snowmass year.
unsigned short m_doseReportingLevel
0 - no data, 1 - summary only, 2 - ntuple
std::string m_storeOccupancyEventsName
PXDOccupancyEvents StoreArray name.
std::map< VxdID, SensorData > m_sensorData
Struct to hold sensor-wise background data.
virtual void initialize() override
Initialize module.
std::string m_relTrueHitsSimHitsName
PXDTrueHitsToPXDSimHits RelationArray name.
double m_integrationTime
Integration time of PXD.
std::string m_relDigitsMCParticlesName
StoreArray name of PXDDigits to MCParticles relation.
const std::string c_niel_neutronFile
NIEL-correction file for neutrons.
virtual void event() override
Event processing.
double getSensorMass(VxdID sensorID) const
Return mass of the sensor with the given sensor ID.
static const unsigned short c_reportNone
No reporting.
virtual void endRun() override
End-of-run tasks.
std::string m_storeNeutronFluxesName
PXDNeutronFluxEvents StoreArray name.
std::string m_storeTrueHitsName
PXDTrueHits StoreArray name.
virtual void terminate() override
Final summary and cleanup.
std::unique_ptr< TNiel > m_nielNeutrons
Pointer to Niel table for neutrons.
std::string m_storeMCParticlesName
MCParticles StoreArray name.
std::unique_ptr< TNiel > m_nielProtons
Pointer to Niel table for protons.
double getSensorThickness(VxdID sensorID) const
Return thickness of the sensor with the given sensor ID.
const PXD::SensorInfo & getInfo(VxdID sensorID) const
This is a shortcut to getting PXD::SensorInfo from the GeoCache.
std::unique_ptr< TNiel > m_nielPions
Pointer to Niel table for pions.
virtual void beginRun() override
Start-of-run initializations.
const std::string c_niel_electronFile
NIEL-correction file for electrons.
std::string m_storeBgMetaDataName
Name of the persistent BackgroundMetaDta object.
double getSensorArea(VxdID sensorID) const
Return area of the sensor with the given sensor ID.
std::string m_outputDirectoryName
Path to directory where output data will be stored.
std::string m_storeDigitsName
PXDDigits StoreArray name.
std::string m_storeClustersName
PXDClusters StoreArray name.
std::unique_ptr< TNiel > m_nielElectrons
Pointer to Niel table for electrons.
const std::string c_niel_protonFile
NIEL-correction file for protons.
const std::string c_niel_pionFile
NIEL-correction file for pions.
std::string m_relParticlesTrueHitsName
MCParticlesToPXDTrueHits RelationArray name.
std::string m_componentName
Name of the current bg component.
double m_componentTime
Time of current component.
unsigned short m_occupancyReportingLevel
0 - no data, 1 - summary only, 2 - ntuple
std::string m_storeSimHitsName
PXDSimHits StoreArray name.
std::string m_relDigitsTrueHitsName
StoreArray name of PXDDigits to PXDTrueHits relation.
Specific implementation of SensorInfo for PXD Sensors which provides additional pixel specific inform...
Definition SensorInfo.h:23
Low-level class to create/modify relations between StoreArrays.
const std::string & getName() const
Return name under which the object is saved in the DataStore.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Accessor to arrays stored in the data store.
Definition StoreArray.h:113
T * appendNew()
Construct a new T object at the end of the array.
Definition StoreArray.h:246
Type-safe access to single objects in the data store.
Definition StoreObjPtr.h:96
static const double us
[microsecond]
Definition Unit.h:97
static const double J
[joule]
Definition Unit.h:116
static const double MeV
[megaelectronvolt]
Definition Unit.h:114
static const double s
[second]
Definition Unit.h:95
float getGlobalTime() const override
Return the time of the electron deposition.
Definition VXDSimHit.h:78
int getPDGcode() const
Return the PDG code of the particle causing the electron deposition.
Definition VXDSimHit.h:68
Class to uniquely identify a any structure of the PXD and SVD.
Definition VxdID.h:33
baseType getSensorNumber() const
Get the sensor id.
Definition VxdID.h:100
void setID(baseType id)
Set the unique id.
Definition VxdID.h:105
baseType getLadderNumber() const
Get the ladder id.
Definition VxdID.h:98
baseType getLayerNumber() const
Get the layer id.
Definition VxdID.h:96
TNiel - the class providing values for NIEL factors.
Definition niel_fun.h:17
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Namespace to encapsulate code needed for simulation and reconstrucion of the PXD.
Abstract base class for different kinds of events.
STL namespace.