Skip to content

Commit 10ea0e9

Browse files
sawenzelsawenzel
authored
[WIP] HMPID digitizer improvements (#1616)
This PR brings the following improvements to the HMPID digitization: - treatment of pad-crosstalk (necessitating some refactoring of Mathison functions in digit class) - treatment of labels - zero-suppression logic - correct calculation of total charge (ported from AliRoot, previously missing) - cleanup - using generic RootTreeWriter for digit and label output authors: @gvolpe79 and @sawenzel (coding sprint)
1 parent 212917a commit 10ea0e9

8 files changed

Lines changed: 220 additions & 113 deletions

File tree

‎Detectors/HMPID/base/include/HMPIDBase/Digit.h‎

Lines changed: 38 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -31,48 +31,60 @@ class Digit : public DigitBase
3131
Digit(float charge) : mQ(charge) {}
3232
float getCharge() const { return mQ; }
3333
int getPadID() const { return mPad; }
34+
// convenience conversion to x-y pad coordinates
35+
int getPx() const { return Param::A2X(mPad); }
36+
int getPy() const { return Param::A2Y(mPad); }
3437

35-
Digit(HitType const& hit)
38+
static void getPadAndTotalCharge(HitType const& hit, int& chamber, int& pc, int& px, int& py, float& totalcharge)
3639
{
37-
int pc, px, py;
3840
float localX;
3941
float localY;
40-
41-
// chamber number is in detID
42-
const auto chamber = hit.GetDetectorID();
43-
42+
chamber = hit.GetDetectorID();
4443
double tmp[3] = { hit.GetX(), hit.GetY(), hit.GetZ() };
4544
Param::Instance()->Mars2Lors(chamber, tmp, localX, localY);
46-
4745
Param::Lors2Pad(localX, localY, pc, px, py);
4846

49-
// TODO: check if this digit is valid
50-
// mark as invalid otherwise
51-
52-
// calculate pad id
53-
mPad = Param::Abs(chamber, pc, px, py);
47+
totalcharge = Digit::QdcTot(hit.GetEnergyLoss(), hit.GetTime(), pc, px, py, localX, localY);
48+
}
5449

55-
// calculate charge
56-
mQ = /*fQHit **/ InMathieson(localX, localY);
50+
static float getFractionalContributionForPad(HitType const& hit, int somepad)
51+
{
52+
float localX;
53+
float localY;
5754

58-
// what about time stamp??
55+
// chamber number is in detID
56+
const auto chamber = hit.GetDetectorID();
57+
double tmp[3] = { hit.GetX(), hit.GetY(), hit.GetZ() };
58+
// converting chamber id and hit coordiates to local coordinates
59+
Param::Instance()->Mars2Lors(chamber, tmp, localX, localY);
60+
// calculate charge fraction in given pad
61+
return Digit::InMathieson(localX, localY, somepad);
5962
}
6063

64+
Digit(int pad, float charge) : mPad(pad), mQ(charge) {}
65+
66+
// add charge to existing digit
67+
void addCharge(float q) { mQ += q; }
68+
6169
private:
6270
float mQ = 0.;
6371
int mPad = 0.; // -1 indicates invalid digit
6472

65-
float LorsX() const { return Param::LorsX(Param::A2P(mPad), Param::A2X(mPad)); } //center of the pad x, [cm]
66-
float LorsY() const { return Param::LorsY(Param::A2P(mPad), Param::A2Y(mPad)); } //center of the pad y, [cm]
73+
static float LorsX(int pad) { return Param::LorsX(Param::A2P(pad), Param::A2X(pad)); } //center of the pad x, [cm]
74+
static float LorsY(int pad) { return Param::LorsY(Param::A2P(pad), Param::A2Y(pad)); } //center of the pad y, [cm]
75+
76+
// determines the total charge created by a hit
77+
// might modify the localX, localY coordiates associated to the hit
78+
static float QdcTot(float e, float time, int pc, int px, int py, float& localX, float& localY);
6779

68-
float IntPartMathiX(float x) const
80+
static float IntPartMathiX(float x, int pad)
6981
{
7082
// Integration of Mathieson.
7183
// This is the answer to electrostatic problem of charge distrubution in MWPC described elsewhere. (NIM A370(1988)602-603)
7284
// Arguments: x,y- position of the center of Mathieson distribution
7385
// Returns: a charge fraction [0-1] imposed into the pad
74-
auto shift1 = -LorsX() + 0.5 * Param::SizePadX();
75-
auto shift2 = -LorsX() - 0.5 * Param::SizePadX();
86+
auto shift1 = -LorsX(pad) + 0.5 * Param::SizePadX();
87+
auto shift2 = -LorsX(pad) - 0.5 * Param::SizePadX();
7688

7789
auto ux1 = Param::SqrtK3x() * TMath::TanH(Param::K2x() * (x + shift1) / Param::PitchAnodeCathode());
7890
auto ux2 = Param::SqrtK3x() * TMath::TanH(Param::K2x() * (x + shift2) / o2::hmpid::Param::PitchAnodeCathode());
@@ -81,27 +93,27 @@ class Digit : public DigitBase
8193
}
8294
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8395

84-
Double_t IntPartMathiY(Double_t y) const
96+
static Double_t IntPartMathiY(Double_t y, int pad)
8597
{
8698
// Integration of Mathieson.
8799
// This is the answer to electrostatic problem of charge distrubution in MWPC described elsewhere. (NIM A370(1988)602-603)
88100
// Arguments: x,y- position of the center of Mathieson distribution
89101
// Returns: a charge fraction [0-1] imposed into the pad
90-
Double_t shift1 = -LorsY() + 0.5 * o2::hmpid::Param::SizePadY();
91-
Double_t shift2 = -LorsY() - 0.5 * o2::hmpid::Param::SizePadY();
102+
Double_t shift1 = -LorsY(pad) + 0.5 * o2::hmpid::Param::SizePadY();
103+
Double_t shift2 = -LorsY(pad) - 0.5 * o2::hmpid::Param::SizePadY();
92104

93105
Double_t uy1 = Param::SqrtK3y() * TMath::TanH(Param::K2y() * (y + shift1) / Param::PitchAnodeCathode());
94106
Double_t uy2 = Param::SqrtK3y() * TMath::TanH(Param::K2y() * (y + shift2) / Param::PitchAnodeCathode());
95107

96108
return Param::K4y() * (TMath::ATan(uy2) - TMath::ATan(uy1));
97109
}
98110

99-
float InMathieson(float localX, float localY) const
111+
static float InMathieson(float localX, float localY, int pad)
100112
{
101-
return 4. * IntPartMathiX(localX) * IntPartMathiY(localY);
113+
return 4. * Digit::IntPartMathiX(localX, pad) * Digit::IntPartMathiY(localY, pad);
102114
}
103115

104-
ClassDefNV(Digit, 1);
116+
ClassDefNV(Digit, 2);
105117
};
106118

107119
} // namespace hmpid

‎Detectors/HMPID/base/src/Digit.cxx‎

Lines changed: 36 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,8 +9,43 @@
99
// or submit itself to any jurisdiction.
1010

1111
#include "HMPIDBase/Digit.h"
12-
#include <iostream>
12+
#include "HMPIDBase/Param.h"
13+
#include "TRandom.h"
14+
#include "TMath.h"
1315

1416
using namespace o2::hmpid;
1517

1618
ClassImp(o2::hmpid::Digit);
19+
20+
float Digit::QdcTot(float e, float time, int pc, int px, int py, float& localX, float& localY)
21+
{
22+
// Samples total charge associated to a hit
23+
// Arguments: e- hit energy [GeV] for mip Eloss for photon Etot
24+
// Returns: total QDC
25+
float Q = 0;
26+
if (time > 1.2e-6) {
27+
Q = 0;
28+
}
29+
if (py < 0) {
30+
return 0;
31+
} else {
32+
float y = Param::LorsY(pc, py);
33+
localY = ((y - localY) > 0) ? y - 0.2 : y + 0.2; //shift to the nearest anod wire
34+
35+
float x = (localX > 66.6) ? localX - 66.6 : localX; //sagita is for PC (0-64) and not for chamber
36+
float qdcEle = 34.06311 + 0.2337070 * x + 5.807476e-3 * x * x - 2.956471e-04 * x * x * x + 2.310001e-06 * x * x * x * x; //reparametrised from DiMauro
37+
38+
int iNele = int((e / 26e-9) * 0.8);
39+
if (iNele < 1) {
40+
iNele = 1; //number of electrons created by hit, if photon e=0 implies iNele=1
41+
}
42+
for (Int_t i = 1; i <= iNele; i++) {
43+
double rnd = gRandom->Rndm();
44+
if (rnd == 0) {
45+
rnd = 1e-12; //1e-12 is a protection against 0 from rndm
46+
}
47+
Q -= qdcEle * TMath::Log(rnd);
48+
}
49+
}
50+
return Q;
51+
}

‎Detectors/HMPID/simulation/include/HMPIDSimulation/HMPIDDigitizer.h‎

Lines changed: 28 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,8 @@
1313

1414
#include "HMPIDBase/Digit.h"
1515
#include "HMPIDSimulation/Detector.h" // for the hit
16+
#include "SimulationDataFormat/MCTruthContainer.h"
17+
#include "SimulationDataFormat/MCCompLabel.h"
1618
#include <vector>
1719

1820
namespace o2
@@ -27,10 +29,24 @@ class HMPIDDigitizer
2729
void setEventID(int eventID) { mEventID = eventID; }
2830
void setSrcID(int sID) { mSrcID = sID; }
2931

32+
// user can pass a label container to be filled ... this activates the label mechanism
33+
// the passed label container can be readout after call to process
34+
void setLabelContainer(o2::dataformats::MCTruthContainer<o2::MCCompLabel>* labels)
35+
{
36+
mRegisteredLabelContainer = labels;
37+
}
38+
3039
// this will process hits and fill the digit vector with digits which are finalized
3140
void process(std::vector<o2::hmpid::HitType> const&, std::vector<o2::hmpid::Digit>& digit);
3241

3342
private:
43+
void zeroSuppress(std::vector<o2::hmpid::Digit> const& digits, std::vector<o2::hmpid::Digit>& newdigits,
44+
o2::dataformats::MCTruthContainer<o2::MCCompLabel> const& labels,
45+
o2::dataformats::MCTruthContainer<o2::MCCompLabel>* newlabels);
46+
47+
float getThreshold(o2::hmpid::Digit const&) const; // gives back threshold to apply for a certain digit
48+
// (using noise and other tables for pad)
49+
3450
double mTime = 0.;
3551
int mEventID = 0;
3652
int mSrcID = 0;
@@ -39,7 +55,18 @@ class HMPIDDigitizer
3955
std::vector<o2::hmpid::Digit> mSummable;
4056
std::vector<o2::hmpid::Digit> mFinal;
4157

42-
// other stuff needed for digitizaton
58+
std::vector<o2::hmpid::Digit> mDigits; // internal store for digits
59+
60+
//static constexpr int HMPID_NUMBEROFPADS = 161280;
61+
//std::array<short, HMPID_NUMBEROFPADS> mIndexForPad = { -1 }; //! mapping of pad to digit index
62+
63+
std::map<int, short> mIndexForPad; //! logarithmic mapping of pad to digit index
64+
65+
std::vector<int> mInvolvedPads; //! list of pads where digits created
66+
67+
// other stuff needed for digitization
68+
o2::dataformats::MCTruthContainer<o2::MCCompLabel> mTmpLabelContainer; // temp label container as workspace
69+
o2::dataformats::MCTruthContainer<o2::MCCompLabel>* mRegisteredLabelContainer = nullptr; // label container to be filled
4370

4471
ClassDefNV(HMPIDDigitizer, 1);
4572
};

‎Detectors/HMPID/simulation/src/HMPIDDigitizer.cxx‎

Lines changed: 90 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -15,13 +15,99 @@ using namespace o2::hmpid;
1515

1616
ClassImp(HMPIDDigitizer);
1717

18+
float HMPIDDigitizer::getThreshold(o2::hmpid::Digit const& digiti) const
19+
{
20+
// TODO: implement like in AliRoot some thresholding depending on conditions ...
21+
return 4.;
22+
}
23+
24+
// applies threshold to digits; removes the ones below a certain charge threshold
25+
void HMPIDDigitizer::zeroSuppress(std::vector<o2::hmpid::Digit> const& digits, std::vector<o2::hmpid::Digit>& newdigits,
26+
o2::dataformats::MCTruthContainer<o2::MCCompLabel> const& labels,
27+
o2::dataformats::MCTruthContainer<o2::MCCompLabel>* newlabels)
28+
{
29+
int index = 0;
30+
for (auto& digit : digits) {
31+
if (digit.getCharge() >= getThreshold(digit)) {
32+
newdigits.push_back(digit);
33+
34+
if (newlabels) {
35+
// copy the labels to the new place with the right new index
36+
newlabels->addElements(newdigits.size() - 1, labels.getLabels(index));
37+
}
38+
}
39+
index++;
40+
}
41+
}
42+
1843
// this will process hits and fill the digit vector with digits which are finalized
1944
void HMPIDDigitizer::process(std::vector<o2::hmpid::HitType> const& hits, std::vector<o2::hmpid::Digit>& digits)
2045
{
21-
// this is a very simple variant that creates one digit from one hit
22-
// conversion is done in the digit constructor
23-
// TODO: introduce cross-talk, pile-up, etc.
46+
mIndexForPad.clear();
47+
mInvolvedPads.clear();
48+
2449
for (auto& hit : hits) {
25-
digits.emplace_back(hit);
50+
int chamber, pc, px, py;
51+
float totalQ;
52+
// retrieves center pad and the total charge
53+
Digit::getPadAndTotalCharge(hit, chamber, pc, px, py, totalQ);
54+
55+
if (px < 0 || py < 0) {
56+
continue;
57+
}
58+
59+
// determine which pads to loop over
60+
std::array<int, 9> allpads;
61+
int counter = 0;
62+
for (int nx = -1; nx <= 1; ++nx) {
63+
for (int ny = -1; ny <= 1; ++ny) {
64+
allpads[counter] = Param::Abs(chamber, pc, px + nx, py + ny);
65+
counter++;
66+
}
67+
}
68+
69+
for (auto& pad : allpads) {
70+
auto iter = mIndexForPad.find(pad);
71+
int index = -1;
72+
if (iter != mIndexForPad.end()) {
73+
index = iter->second;
74+
}
75+
// auto index = mIndexForPad[pad];
76+
float fraction = Digit::getFractionalContributionForPad(hit, pad);
77+
// LOG(INFO) << "FRACTION ON PAD " << pad << " IS " << fraction;
78+
if (index != -1) {
79+
// digit exists ... reuse
80+
auto& digit = mDigits[index];
81+
digit.addCharge(totalQ * fraction);
82+
83+
if (mRegisteredLabelContainer) {
84+
auto labels = mTmpLabelContainer.getLabels(index);
85+
o2::MCCompLabel newlabel(hit.GetTrackID(), mEventID, mSrcID);
86+
bool newlabelneeded = true;
87+
for (auto& l : labels) {
88+
if (l == newlabel) {
89+
newlabelneeded = false;
90+
break;
91+
}
92+
}
93+
if (newlabelneeded) {
94+
mTmpLabelContainer.addElementRandomAccess(index, newlabel);
95+
}
96+
}
97+
} else {
98+
// create digit ... and register
99+
mDigits.emplace_back(pad, totalQ * fraction);
100+
mIndexForPad[pad] = mDigits.size() - 1;
101+
mInvolvedPads.emplace_back(pad);
102+
103+
if (mRegisteredLabelContainer) {
104+
// add label for this digit
105+
mTmpLabelContainer.addElement(mDigits.size() - 1, o2::MCCompLabel(hit.GetTrackID(), mEventID, mSrcID));
106+
}
107+
}
108+
}
26109
}
110+
111+
// apply zero suppression (removes digits in the noise)
112+
zeroSuppress(mDigits, digits, mTmpLabelContainer, mRegisteredLabelContainer);
27113
}

‎Steer/DigitizerWorkflow/CMakeLists.txt‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,6 @@ O2_GENERATE_EXECUTABLE(
3737
src/EMCALDigitizerSpec.cxx
3838
src/EMCALDigitWriterSpec.cxx
3939
src/HMPIDDigitizerSpec.cxx
40-
src/HMPIDDigitWriterSpec.cxx
4140
src/MCHDigitizerSpec.cxx
4241
src/TRDDigitizerSpec.cxx
4342
src/GRPUpdaterSpec.cxx

‎Steer/DigitizerWorkflow/src/HMPIDDigitWriterSpec.cxx‎

Lines changed: 0 additions & 73 deletions
This file was deleted.

0 commit comments

Comments
 (0)