9#include <pxd/modules/pxdSimulation/PXDDigitizerModule.h>
10#include <vxd/geometry/GeoCache.h>
11#include <vxd/geometry/GeoTools.h>
12#include <pxd/reconstruction/PXDGainCalibrator.h>
13#include <pxd/reconstruction/PXDPixelMasker.h>
15#include <framework/logging/Logger.h>
16#include <framework/gearbox/Unit.h>
17#include <framework/gearbox/Const.h>
18#include <framework/datastore/DataStore.h>
19#include <framework/datastore/StoreArray.h>
20#include <framework/datastore/RelationArray.h>
21#include <framework/datastore/RelationIndex.h>
23#include <mdst/dataobjects/MCParticle.h>
24#include <pxd/dataobjects/PXDTrueHit.h>
25#include <pxd/dataobjects/PXDDigit.h>
55 "Noise added by the electronics, set in ENC", 150.0);
58 "MCParticle collection name",
string(
""));
66 "Apply Poisson smearing of electrons collected on pixels?",
true);
73 "Split Signalpoints in smaller groups of N electrons (in e)", 100);
75 "Time step for tracking electron groups in readout plane (in ns)",
78 "Maximum number of steps when propagating electrons", 200);
79 addParam(
"Gq",
m_gq,
"Gq of a pixel in nA/electron", 0.6);
80 addParam(
"ADCUnit",
m_ADCUnit,
"Slope of the linear ADC transfer curve in nA/ADU", 130.0);
84 "Time window during which the PXD is not collecting charge in nano seconds", 400.0);
86 "Digits from gated rows not sent to DHH ",
true);
88 "Time window during which digits from gated rows are not sent to DHH in nano seconds", 1400.0);
90 "Constant time delay between bunch crossing and switching on triggergate in nano seconds", 0.0);
128 m_nGates = gTools->getNumberOfPXDReadoutGates();
135 B2DEBUG(20,
"PXDDigitizer Parameters (in system units, *=calculated +=set in xml):");
136 B2DEBUG(20,
" --> ElectronicEffects: " << (
m_applyNoise ?
"true" :
"false"));
137 B2DEBUG(20,
" --> ElectronicNoise: " <<
m_elNoise <<
" e-");
138 B2DEBUG(20,
" --> +ChargeThreshold: " <<
"set in xml by sensor, nominal 4 ADU");
139 B2DEBUG(20,
" --> *NoiseFraction: " <<
"set in xml by sensor, nominal 1.0e-5");
152 B2DEBUG(20,
" --> PoissonSmearing: " << (
m_applyPoisson ?
"true" :
"false"));
153 B2DEBUG(20,
" --> +IntegrationWindow: " << (
m_applyWindow ?
"true" :
"false") <<
", size defined in xml");
155 B2DEBUG(20,
" --> ElectronGroupSize: " <<
m_elGroupSize <<
" e-");
156 B2DEBUG(20,
" --> ElectronStepTime: " <<
m_elStepTime <<
" ns");
158 B2DEBUG(20,
" --> ADU unit: " <<
m_eToADU <<
" e-/ADU");
185 m_gated = (*m_storePXDTiming).isGated();
190 B2DEBUG(20,
"Reading from PXDInjectionBGTiming: Event has TriggerGate " <<
m_triggerGate <<
" and is gated " <<
m_gated);
194 B2DEBUG(20,
"Gating start time " << gatingStartTime <<
"ns ");
197 if (gatingStartTime >= 0) {
203 if (lastGated >= firstGated) {
205 B2DEBUG(20,
"Gated channel interval from " << firstGated <<
" to " << lastGated);
209 B2DEBUG(20,
"Gated channel interval from " << firstGated <<
" to " <<
m_nGates);
210 B2DEBUG(20,
"Gated channel interval from " << 0 <<
" to " << lastGated);
225 for (Sensors::value_type& sensor :
m_sensors) {
226 sensor.second.clear();
237 RelationArray mcParticlesToSimHits(storeMCParticles, storeSimHits,
239 RelationArray mcParticlesToTrueHits(storeMCParticles, storeTrueHits);
246 RelationArray relDigitMCParticle(storeDigits, storeMCParticles,
248 relDigitMCParticle.
clear();
252 relDigitTrueHit.
clear();
254 unsigned int nSimHits = storeSimHits.
getEntries();
264 for (
unsigned int i = 0; i < nSimHits; ++i) {
281 "Could not find MCParticle which produced PXDSimhit " << i);
287 if (trueRel && trueRel->
weight > 0) {
298 "SensorInformation for Sensor " << sensorID <<
" not found, make sure that the geometry is set up correctly");
300 B2DEBUG(20,
"Sensor Parameters for Sensor " << sensorID <<
": " << endl
321 "Processing hit " << i <<
" in Sensor " << sensorID <<
", related to MCParticle " <<
m_currentParticle);
338 const ROOT::Math::XYZVector startPoint =
m_currentHit->getPosIn();
339 const ROOT::Math::XYZVector stopPoint =
m_currentHit->getPosOut();
348 int distanceToTriggerGate = 0;
358 "Checking if hit is in timeframe " << gateMinTime <<
" <= " << currentHitTime <<
359 " <= " << gateMaxTime);
371 "Checking if hit is in gated timeframe " << gatingStartTime <<
" <= " << currentHitTime <<
384 const ROOT::Math::XYZVector startPoint =
m_currentHit->getPosIn();
385 ROOT::Math::XYZVector stopPoint =
m_currentHit->getPosOut();
386 double dx = stopPoint.X() - startPoint.X();
387 double dy = stopPoint.Y() - startPoint.Y();
388 double dz = stopPoint.Z() - startPoint.Z();
389 double trackLength2 = dx * dx + dy * dy + dz * dz;
401 double lastFraction {0};
402 double lastElectrons {0};
404 for (
auto& segment : segments) {
407 const double f = (segment.first + lastFraction) * 0.5;
408 const double e = segment.second - lastElectrons;
410 std::tie(lastFraction, lastElectrons) = segment;
413 stopPoint.SetXYZ(startPoint.X() + f * dx, startPoint.Y() + f * dy, startPoint.Z() + f * dz);
419inline double pol3(
double x,
const double* c)
421 return c[0] + x * (c[1] + x * (c[2] + x * c[3]));
430 static const double cx[] = {0, 2.611995e-01, -1.948316e+01, 9.167064e+02};
431 static const double cu[] = {4.223817e-04, -1.041434e-03, 5.139596e-02, -1.180229e+00};
432 static const double cv[] = {4.085681e-04, -8.615459e-04, 5.706439e-02, -1.774319e+00};
433 static const double ct[] = {2.823743e-04, -1.138627e-06, 4.310888e-09, -1.559542e-11};
435 double dz = 0.5 * info.getThickness() - info.getGateDepth() - r.Z(), adz = fabs(dz);
439 double sigmaDrift_u = pol3(adz, cu);
440 double sigmaDrift_v = pol3(adz, cv);
441 double sigmaDiffus = pol3(info.getTemperature() - 300.0,
445 double groupCharge = electrons / nGroups;
447 double du = gRandom->Gaus(0.0, sigmaDrift_u), dv = gRandom->Gaus(0.0, sigmaDrift_v);
448 double uPos = r.X() + du + dx;
449 double vPos = r.Y() + dv;
452 int id = info.getTrappedID(uPos, vPos);
460 dv = gRandom->Gaus(0.0, sigmaDiffus), du = gRandom->Gaus(0.0, sigmaDiffus);
466 int iu = info.getUCellID(uPos, vPos, 1), iv = info.getVCellID(vPos, 1);
481 if (charge > (1000. *
Unit::e))
482 charge = gRandom->Gaus(charge,
sqrt(charge));
485 charge = gRandom->PoissonD(charge);
499 for (Sensors::value_type& sensor :
m_sensors) {
500 Sensor& s = sensor.second;
506 int nV = info.getVCells();
507 int nPixels = gRandom->Poisson(info.getNoiseFraction() * nU * nV);
513 for (
int i = 0; i < nPixels; ++i) {
514 Digit d(gRandom->Integer(nU), gRandom->Integer(nV));
523 [gate](
auto interval) {return gate >= interval.first && gate < interval.second;});
531 RelationArray relDigitMCParticle(storeDigits, storeMCParticles,
537 for (Sensors::value_type& sensor :
m_sensors) {
538 auto sensorID = sensor.first;
543 for (Sensor::value_type& digitAndValue : sensor.second) {
544 const Digit& d = digitAndValue.first;
557 charge = round(charge * gain /
m_eToADU + pedestal);
560 charge = std::min(charge, 255.0);
563 charge = charge - round(pedestal);
572 if (sensorID.getLayerNumber() == 1)
587 PXDDigit(sensorID, d.u(), d.v(), charge));
596 relDigitTrueHit.
add(digIndex, v.
truehits().begin(),
static const ParticleType photon
photon particle
static std::string relationName(const std::string &fromName, const std::string &toName, std::string const &namedRelation="")
Return storage name for a relation between two arrays of the given names.
static std::string arrayName(const TClass *t, const std::string &name)
Return the storage name for an object of the given TClass and name.
void setDescription(const std::string &description)
Sets the description of the module.
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
void setReturnValue(int value)
Sets the return value for this module as integer.
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Class representing the charge and particle contributions for one pixel.
const relations_map & truehits() const
Return the map containing all truehit contributions to the pixel charge.
const relations_map & particles() const
Return the map containing all particle contributions to the pixel charge.
double charge() const
Return the charge collected in the pixel.
Class to represent the coordinates of one pixel.
double m_pedestalRMS
RMS pedestal in ADU.
void processHit()
Process one PXDSimHit by dividing the step in smaller steps and drifting the charge.
double m_pxdIntegrationTime
Integration time for each gate of the PXD in ns.
double m_pedestalMean
Mean pedestal in ADU.
bool m_gatingWithoutReadout
Digits from gated rows not sent to DHH.
void initialize() override final
Initialize the module and check the parameters.
std::string m_relTrueHitSimHitName
Name of the relation between PXDTrueHits and PXDSimHits.
void saveDigits()
Save all digits to the datastore.
int m_currentParticle
Index of the particle which caused the current hit.
double m_ADCUnit
Slope of the linear ADC transfer curve in nA/ADU.
Sensors m_sensors
Structure containing all existing sensors.
double m_elStepTime
Timeframe for one random walk step.
std::string m_relDigitMCParticleName
Name of the relation between PXDDigits and MCParticles.
PXDDigitizerModule()
Constructor.
int m_triggerGate
PXD triggergate.
int m_elMaxSteps
Maximum number of random walks before abort.
bool m_applyWindow
Whether or not to apply a time window cut.
std::string m_storeTrueHitsName
Name of the collection for the PXDTrueHits.
double m_timePerGate
Time needed to sample and clear a readout gate.
double m_gatingWithoutReadoutTime
Time window during which digits from gated rows are not sent to DHH.
void driftCharge(const ROOT::Math::XYZVector &position, double electrons)
Drift the charge inside the silicon.
const PXDSimHit * m_currentHit
Pointer to the PXDSimhit currently digitized.
std::vector< float > m_gatingStartTimes
Vector of start times for gating.
std::string m_storeMCParticlesName
Name of the collection for the MCParticles.
double m_noiseFraction
Fraction of noisy pixels per sensor.
bool m_applyNoise
Whether or not to apply noise.
void addNoiseDigits()
Add pure noise digits to the Sensors.
void event() override final
Digitize one event.
int m_nGates
Number of readout gates (or total number of Switcher channels)
double m_chargeThreshold
Zero-suppression threshold in ADU.
bool m_applyPoisson
Whether or not to apply poission fluctuation of charge.
Sensor * m_currentSensor
Pointer to the sensor in which the current hit occurred.
double m_chargeThresholdElectrons
... and its equivalent in electrons
std::vector< std::pair< int, int > > m_gatedChannelIntervals
Vector of gated readout channels.
ROOT::Math::XYZVector m_currentBField
Current magnetic field.
double m_elNoise
Amount of noise to apply.
std::string m_storeDigitsName
Name of the collection for the PXDDigits.
const SensorInfo * m_currentSensorInfo
Pointer to the SensorInfo of the current sensor.
std::string m_relDigitTrueHitName
Name of the relation between PXDDigits and PXDTrueHits.
double m_gatingTime
Time window during which the PXD is not collecting charge.
void beginRun() override final
Initialize the list of existing PXD Sensors.
double m_eToADU
ENC equivalent of 1 ADU.
bool checkIfGated(int gate)
Check if gate was read while in gated mode.
double m_gq
g_q of a pixel in nA/electrons.
StoreObjPtr< PXDInjectionBGTiming > m_storePXDTiming
Input array for timings.
double m_hwdelay
Hardware delay between time of bunch crossing and switching on triggergate in ns.
int m_elGroupSize
Max number of electrons per random walk.
bool m_gated
Gated mode flag.
int m_currentTrueHit
Index of the TrueHit the current hit belongs to.
std::string m_relMCParticleSimHitName
Name of the relation between MCParticles and PXDSimHits.
double addNoise(double charge)
Calculate the noise contribution to one pixel with given charge.
double m_segmentLength
Max.
std::string m_storeSimHitsName
Name of the collection for the PXDSimhits.
float getGainCorrection(VxdID id, unsigned int uid, unsigned int vid) const
Get gain correction.
static PXDGainCalibrator & getInstance()
Main (and only) way to access the PXDGainCalibrator.
static PXDPixelMasker & getInstance()
Main (and only) way to access the PXDPixelMasker.
Specific implementation of SensorInfo for PXD Sensors which provides additional pixel specific inform...
Low-level class to create/modify relations between StoreArrays.
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
void clear() override
Clear all elements from the relation.
Provides access to fast ( O(log n) ) bi-directional lookups on a specified relation.
const Element * getFirstElementTo(const TO &to) const
Return a pointer to the first relation Element of the given object.
RelationIndexContainer< FROM, TO >::Element Element
Struct representing a single element in the index.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Accessor to arrays stored in the data store.
T * appendNew()
Construct a new T object at the end of the array.
int getEntries() const
Get the number of objects in the array.
void clear() override
Delete all entries in this array.
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
static const double mm
[millimeters]
static const double e
Standard of [electric charge].
static const double um
[micrometers]
static const double ns
Standard of [time].
static const double T
[tesla]
Class to facilitate easy access to sensor information of the VXD like coordinate transformations or p...
const std::set< Belle2::VxdID > getLayers(SensorInfoBase::SensorType sensortype=SensorInfoBase::VXD)
Return a set of all known Layers.
const SensorInfoBase & getSensorInfo(Belle2::VxdID id) const
Return a reference to the SensorInfo of a given SensorID.
const std::set< Belle2::VxdID > & getSensors(Belle2::VxdID ladder) const
Return a set of all sensor IDs belonging to a given ladder.
static GeoCache & getInstance()
Return a reference to the singleton instance.
const GeoTools * getGeoTools()
Return a raw pointer to a GeoTools object.
const std::set< Belle2::VxdID > & getLadders(Belle2::VxdID layer) const
Return a set of all ladder IDs belonging to a given layer.
int getUCells() const
Return number of pixel/strips in u direction.
Class to uniquely identify a any structure of the PXD and SVD.
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
double sqrt(double a)
sqrt for double
Namespace to encapsulate code needed for simulation and reconstrucion of the PXD.
std::map< Digit, DigitValue > Sensor
Map of all hits in one Sensor.
Abstract base class for different kinds of events.
RelationElement::index_type indexFrom
index of the element from which the relation points.
RelationElement::weight_type weight
weight of the relation.