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/PXDSimHit.h>
25#include <pxd/dataobjects/PXDTrueHit.h>
26#include <pxd/dataobjects/PXDDigit.h>
27#include <pxd/dataobjects/PXDInjectionBGTiming.h>
57 "Noise added by the electronics, set in ENC", 150.0);
60 "MCParticle collection name",
string(
""));
68 "Apply Poisson smearing of electrons collected on pixels?",
true);
75 "Split Signalpoints in smaller groups of N electrons (in e)", 100);
77 "Time step for tracking electron groups in readout plane (in ns)",
80 "Maximum number of steps when propagating electrons", 200);
81 addParam(
"Gq",
m_gq,
"Gq of a pixel in nA/electron", 0.6);
82 addParam(
"ADCUnit",
m_ADCUnit,
"Slope of the linear ADC transfer curve in nA/ADU", 130.0);
86 "Time window during which the PXD is not collecting charge in nano seconds", 400.0);
88 "Digits from gated rows not sent to DHH ",
true);
90 "Time window during which digits from gated rows are not sent to DHH in nano seconds", 1400.0);
92 "Constant time delay between bunch crossing and switching on triggergate in nano seconds", 0.0);
130 m_nGates = gTools->getNumberOfPXDReadoutGates();
137 B2DEBUG(20,
"PXDDigitizer Parameters (in system units, *=calculated +=set in xml):");
138 B2DEBUG(20,
" --> ElectronicEffects: " << (
m_applyNoise ?
"true" :
"false"));
139 B2DEBUG(20,
" --> ElectronicNoise: " <<
m_elNoise <<
" e-");
140 B2DEBUG(20,
" --> +ChargeThreshold: " <<
"set in xml by sensor, nominal 4 ADU");
141 B2DEBUG(20,
" --> *NoiseFraction: " <<
"set in xml by sensor, nominal 1.0e-5");
154 B2DEBUG(20,
" --> PoissonSmearing: " << (
m_applyPoisson ?
"true" :
"false"));
155 B2DEBUG(20,
" --> +IntegrationWindow: " << (
m_applyWindow ?
"true" :
"false") <<
", size defined in xml");
157 B2DEBUG(20,
" --> ElectronGroupSize: " <<
m_elGroupSize <<
" e-");
158 B2DEBUG(20,
" --> ElectronStepTime: " <<
m_elStepTime <<
" ns");
160 B2DEBUG(20,
" --> ADU unit: " <<
m_eToADU <<
" e-/ADU");
187 m_gated = (*m_storePXDTiming).isGated();
192 B2DEBUG(20,
"Reading from PXDInjectionBGTiming: Event has TriggerGate " <<
m_triggerGate <<
" and is gated " <<
m_gated);
196 B2DEBUG(20,
"Gating start time " << gatingStartTime <<
"ns ");
199 if (gatingStartTime >= 0) {
205 if (lastGated >= firstGated) {
207 B2DEBUG(20,
"Gated channel interval from " << firstGated <<
" to " << lastGated);
211 B2DEBUG(20,
"Gated channel interval from " << firstGated <<
" to " <<
m_nGates);
212 B2DEBUG(20,
"Gated channel interval from " << 0 <<
" to " << lastGated);
227 for (Sensors::value_type& sensor :
m_sensors) {
228 sensor.second.clear();
239 RelationArray mcParticlesToSimHits(storeMCParticles, storeSimHits,
241 RelationArray mcParticlesToTrueHits(storeMCParticles, storeTrueHits);
248 RelationArray relDigitMCParticle(storeDigits, storeMCParticles,
250 relDigitMCParticle.
clear();
254 relDigitTrueHit.
clear();
256 unsigned int nSimHits = storeSimHits.
getEntries();
266 for (
unsigned int i = 0; i < nSimHits; ++i) {
283 "Could not find MCParticle which produced PXDSimhit " << i);
289 if (trueRel && trueRel->
weight > 0) {
300 "SensorInformation for Sensor " << sensorID <<
" not found, make sure that the geometry is set up correctly");
302 B2DEBUG(20,
"Sensor Parameters for Sensor " << sensorID <<
": " << endl
323 "Processing hit " << i <<
" in Sensor " << sensorID <<
", related to MCParticle " <<
m_currentParticle);
340 const ROOT::Math::XYZVector startPoint =
m_currentHit->getPosIn();
341 const ROOT::Math::XYZVector stopPoint =
m_currentHit->getPosOut();
350 int distanceToTriggerGate = 0;
360 "Checking if hit is in timeframe " << gateMinTime <<
" <= " << currentHitTime <<
361 " <= " << gateMaxTime);
373 "Checking if hit is in gated timeframe " << gatingStartTime <<
" <= " << currentHitTime <<
386 const ROOT::Math::XYZVector startPoint =
m_currentHit->getPosIn();
387 ROOT::Math::XYZVector stopPoint =
m_currentHit->getPosOut();
388 double dx = stopPoint.X() - startPoint.X();
389 double dy = stopPoint.Y() - startPoint.Y();
390 double dz = stopPoint.Z() - startPoint.Z();
391 double trackLength2 = dx * dx + dy * dy + dz * dz;
403 double lastFraction {0};
404 double lastElectrons {0};
406 for (
auto& segment : segments) {
409 const double f = (segment.first + lastFraction) * 0.5;
410 const double e = segment.second - lastElectrons;
412 std::tie(lastFraction, lastElectrons) = segment;
415 stopPoint.SetXYZ(startPoint.X() + f * dx, startPoint.Y() + f * dy, startPoint.Z() + f * dz);
421inline double pol3(
double x,
const double* c)
423 return c[0] + x * (c[1] + x * (c[2] + x * c[3]));
432 static const double cx[] = {0, 2.611995e-01, -1.948316e+01, 9.167064e+02};
433 static const double cu[] = {4.223817e-04, -1.041434e-03, 5.139596e-02, -1.180229e+00};
434 static const double cv[] = {4.085681e-04, -8.615459e-04, 5.706439e-02, -1.774319e+00};
435 static const double ct[] = {2.823743e-04, -1.138627e-06, 4.310888e-09, -1.559542e-11};
437 double dz = 0.5 * info.getThickness() - info.getGateDepth() - r.Z(), adz = fabs(dz);
441 double sigmaDrift_u = pol3(adz, cu);
442 double sigmaDrift_v = pol3(adz, cv);
443 double sigmaDiffus = pol3(info.getTemperature() - 300.0,
447 double groupCharge = electrons / nGroups;
449 double du = gRandom->Gaus(0.0, sigmaDrift_u), dv = gRandom->Gaus(0.0, sigmaDrift_v);
450 double uPos = r.X() + du + dx;
451 double vPos = r.Y() + dv;
454 int id = info.getTrappedID(uPos, vPos);
462 dv = gRandom->Gaus(0.0, sigmaDiffus), du = gRandom->Gaus(0.0, sigmaDiffus);
468 int iu = info.getUCellID(uPos, vPos, 1), iv = info.getVCellID(vPos, 1);
483 if (charge > (1000. *
Unit::e))
484 charge = gRandom->Gaus(charge,
sqrt(charge));
487 charge = gRandom->PoissonD(charge);
501 for (Sensors::value_type& sensor :
m_sensors) {
502 Sensor& s = sensor.second;
508 int nV = info.getVCells();
509 int nPixels = gRandom->Poisson(info.getNoiseFraction() * nU * nV);
515 for (
int i = 0; i < nPixels; ++i) {
516 Digit d(gRandom->Integer(nU), gRandom->Integer(nV));
525 [gate](
auto interval) {return gate >= interval.first && gate < interval.second;});
533 RelationArray relDigitMCParticle(storeDigits, storeMCParticles,
539 for (Sensors::value_type& sensor :
m_sensors) {
540 auto sensorID = sensor.first;
545 for (Sensor::value_type& digitAndValue : sensor.second) {
546 const Digit& d = digitAndValue.first;
559 charge = round(charge * gain /
m_eToADU + pedestal);
562 charge = std::min(charge, 255.0);
565 charge = charge - round(pedestal);
574 if (sensorID.getLayerNumber() == 1)
589 PXDDigit(sensorID, d.u(), d.v(), charge));
598 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.