Belle II Software development
CDCDedxPIDCreatorModule.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
9#include <cdc/modules/CDCDedxPID/CDCDedxPIDCreatorModule.h>
10#include <cdc/geometry/CDCGeometryPar.h>
11#include <cdc/translators/LinearGlobalADCCountTranslator.h>
12#include <cdc/translators/RealisticTDCCountTranslator.h>
13#include <cdc/modules/CDCDedxPID/LineHelper.h>
14#include <framework/core/Environment.h>
15#include <TRandom.h>
16#include <cmath>
17#include <algorithm>
18#include <map>
19#include <vector>
20
21using namespace std;
22
23namespace Belle2 {
28
29 using namespace CDC;
30
31 //-----------------------------------------------------------------
33 //-----------------------------------------------------------------
34
35 REG_MODULE(CDCDedxPIDCreator);
36
37 //-----------------------------------------------------------------
38 // Implementation
39 //-----------------------------------------------------------------
40
42
43 {
44 // set module description
45 setDescription("Module that creates PID likelihoods from CDC hit information stored in CDCDedxHits "
46 "using parameterized means and resolutions.");
48
49 addParam("removeLowest", m_removeLowest,
50 "Portion of events with low dE/dx that should be discarded", double(0.05));
51 addParam("removeHighest", m_removeHighest,
52 "Portion of events with high dE/dx that should be discarded", double(0.25));
53 addParam("useBackHalfCurlers", m_useBackHalfCurlers,
54 "Whether to use the back half of curlers", false);
55 addParam("trackLevel", m_trackLevel,
56 "ONLY USEFUL FOR MC: Use track-level MC (generate truncated mean from predicted mean and sigma using MC truth). "
57 "If false, use hit-level MC (use truncated mean determined from hits)", true);
58 addParam("enableDebugOutput", m_enableDebugOutput,
59 "Option to write out debugging information to CDCDedxTracks", true);
60 addParam("likelihoodsName", m_likelihoodsName,
61 "name of CDCDedxLikelihood collection", string(""));
62 addParam("dedxTracksName", m_dedxTracksName,
63 "name of CDCDedxTrack collection", string(""));
64 }
65
69
71 {
72 m_tracks.isRequired();
73 m_hits.isOptional(); // in order to run also with old cdst's where this collection doesn't exist
74 m_mcParticles.isOptional();
75 m_TTDInfo.isOptional();
76 m_likelihoods.registerInDataStore(m_likelihoodsName);
77 m_tracks.registerRelationTo(m_likelihoods);
78 m_dedxTracks.registerInDataStore(m_dedxTracksName);
79 m_tracks.registerRelationTo(m_dedxTracks);
80
81 m_nLayerWires[0] = 1280;
82 for (int i = 1; i < 9; ++i) {
83 m_nLayerWires[i] = m_nLayerWires[i - 1] + 6 * (160 + (i - 1) * 32);
84 }
85
86 if (not m_trackLevel)
87 B2WARNING("Hit-level MC still needs a precise calibration to perform well! Until then please use track-level MC.");
88
89 }
90
92 {
93 // check if CDCDedxHits are present; return if not.
94 if (not m_hits.isValid()) {
96 if (m_warnCount < 10) {
97 B2WARNING("StoreArray 'CDCDedxHits' does not exist, returning. Probably running on old cdst.");
98 } else if (m_warnCount == 10) {
99 B2WARNING("StoreArray 'CDCDedxHits' does not exist, returning. ...message will be suppressed now.");
100 }
101 return;
102 }
103
104 // clear output collections
105 m_likelihoods.clear();
106 m_dedxTracks.clear();
107
108 // CDC geometry parameters and translators
109 const auto& cdcgeo = CDCGeometryPar::Instance();
110 LinearGlobalADCCountTranslator adcTranslator;
111 RealisticTDCCountTranslator tdcTranslator;
112
113 // is data or MC ?
114 bool isData = not Environment::Instance().isMC();
115
116 // track independent calibration constants
117 double runGain = isData ? m_DBRunGain->getRunGain() : 1.0;
118 double timeGain = 1;
119 double timeReso = 1; // this is multiplicative constant
120 if (isData and m_TTDInfo.isValid() and m_TTDInfo->hasInjection()) {
121 timeGain = m_DBInjectTime->getCorrection("mean", m_TTDInfo->isHER(), m_TTDInfo->getTimeSinceLastInjectionInMicroSeconds());
122 timeReso = m_DBInjectTime->getCorrection("reso", m_TTDInfo->isHER(), m_TTDInfo->getTimeSinceLastInjectionInMicroSeconds());
123 }
124 double scale = m_DBScaleFactor->getScaleFactor(); // scale factor to make electron dE/dx ~ 1
125 if (scale == 0) {
126 B2ERROR("Scale factor from DB is zero! Will be set to one");
127 scale = 1;
128 }
129
130 // loop over tracks
131 for (const auto& track : m_tracks) {
132 // track fit result
133 const auto* fitResult = track.getTrackFitResultWithClosestMass(Const::pion);
134 if (not fitResult) {
135 B2WARNING("No related fit for track, skip it.");
136 continue;
137 }
138
139 // hits of this track
140 const auto hits = track.getRelationsTo<CDCDedxHit>();
141 if (hits.size() == 0) continue;
142
143 // track momentum
144 const auto& trackMom = fitResult->getMomentum();
145 double theta = trackMom.Theta();
146 double cosTheta = std::cos(theta);
147 double sinTheta = std::sin(theta);
148
149 // track dependent calibration constants
150 double cosCor = isData ? m_DBCosineCor->getMean(cosTheta) : 1.0;
151 bool isEdge = std::abs(cosTheta + 0.860) < 0.010 or std::abs(cosTheta - 0.955) <= 0.005;
152 double cosEdgeCor = (isData and isEdge) ? m_DBCosEdgeCor->getMean(cosTheta) : 1.0;
153
154 // MC particle
155 const auto* mcParticle = isData ? nullptr : track.getRelated<MCParticle>();
156
157 // debug output
158 CDCDedxTrack* dedxTrack = m_enableDebugOutput ? m_dedxTracks.appendNew() : nullptr;
159 if (dedxTrack) {
160 dedxTrack->m_track = track.getArrayIndex();
161 dedxTrack->m_charge = fitResult->getChargeSign();
162 dedxTrack->m_cosTheta = cosTheta;
163 dedxTrack->m_p = trackMom.R();
164 if (isData and m_TTDInfo.isValid() and m_TTDInfo->hasInjection()) {
165 dedxTrack->m_injring = m_TTDInfo->isHER();
166 dedxTrack->m_injtime = m_TTDInfo->getTimeSinceLastInjectionInMicroSeconds();
167 }
168 if (mcParticle) {
169 dedxTrack->m_pdg = mcParticle->getPDG();
170 dedxTrack->m_mcmass = mcParticle->getMass();
171 const auto* mother = mcParticle->getMother();
172 dedxTrack->m_motherPDG = mother ? mother->getPDG() : 0;
173 const auto& trueMom = mcParticle->getMomentum();
174 dedxTrack->m_pTrue = trueMom.R();
175 dedxTrack->m_cosThetaTrue = std::cos(trueMom.Theta());
176 }
177 dedxTrack->m_scale = scale;
178 dedxTrack->m_cosCor = cosCor;
179 dedxTrack->m_cosEdgeCor = cosEdgeCor;
180 dedxTrack->m_runGain = runGain;
181 dedxTrack->m_timeGain = timeGain;
182 dedxTrack->m_timeReso = timeReso;
183 }
184
185 // loop over hits
186 int lastLayer = -1;
187 double pCDC = 0;
188 std::map<int, DEDX> dedxWires;
189 for (const auto& hit : hits) {
190 // wire numbering: layer and superlayer
191 const auto& wireID = hit.getWireID();
192 int layer = wireID.getILayer(); // layer within superlayer
193 int superlayer = wireID.getISuperLayer();
194 int currentLayer = (superlayer == 0) ? layer : (8 + (superlayer - 1) * 6 + layer); // continuous layer number
195 if (not m_useBackHalfCurlers and currentLayer < lastLayer) break;
196 lastLayer = currentLayer;
197
198 // track momentum at the first hit
199 if (pCDC == 0) pCDC = hit.getPOCAMomentum().R();
200
201 // drift cell
202 double innerRadius = cdcgeo.innerRadiusWireLayer()[currentLayer];
203 double outerRadius = cdcgeo.outerRadiusWireLayer()[currentLayer];
204 const ROOT::Math::XYZVector& wirePosF = cdcgeo.wireForwardPosition(wireID, CDCGeometryPar::c_Aligned);
205 double wireRadius = wirePosF.Rho();
206 int nWires = cdcgeo.nWiresInLayer(currentLayer);
207 double topHeight = outerRadius - wireRadius;
208 double bottomHeight = wireRadius - innerRadius;
209 double topHalfWidth = M_PI * outerRadius / nWires;
210 double bottomHalfWidth = M_PI * innerRadius / nWires;
211 DedxDriftCell cell(DedxPoint(-topHalfWidth, topHeight),
212 DedxPoint(topHalfWidth, topHeight),
213 DedxPoint(bottomHalfWidth, -bottomHeight),
214 DedxPoint(-bottomHalfWidth, -bottomHeight));
215
216 // length of a track within the drift cell
217 double doca = hit.getSignedDOCAXY();
218 double entAng = hit.getEntranceAngle();
219 double celldx = cell.dx(doca, entAng) / sinTheta; // length of a track in the cell
220 if (not cell.isValid()) continue;
221
222 // wire gain calibration (iwire is a continuous wire number)
223 int wire = wireID.getIWire();
224 int iwire = (superlayer == 0) ? 160 * layer + wire : m_nLayerWires[superlayer - 1] + (160 + 32 * (superlayer - 1)) * layer + wire;
225 double wiregain = isData ? m_DBWireGains->getWireGain(iwire) : 1.0;
226
227 // re-scaled (RS) doca and entAng variable: map to square cell
228 double cellHalfWidth = M_PI * wireRadius / nWires;
229 double cellHeight = topHeight + bottomHeight;
230 double cellR = 2 * cellHalfWidth / cellHeight;
231 double tana = std::max(std::min(std::tan(entAng), 1e10), -1e10); // this fixes bug in CDCDedxPIDModule near +-pi/2
232 double docaRS = doca * std::sqrt((1 + cellR * cellR * tana * tana) / (1 + tana * tana));
233 double normDocaRS = docaRS / cellHalfWidth;
234 double entAngRS = std::atan(tana / cellR);
235
236 // one and two dimensional corrections
237 double onedcor = isData ? m_DB1DCell->getMean(currentLayer, entAngRS) : 1.0;
238 double twodcor = isData ? m_DB2DCell->getMean(currentLayer, normDocaRS, entAngRS) : 1.0;
239
240 // total correction
241 double correction = runGain * cosCor * cosEdgeCor * timeGain * wiregain * twodcor * onedcor;
242
243 // calibrated ADC count
244 double adcCount = isData ? m_DBNonlADC->getCorrectedADC(hit.getADCCount(), currentLayer) : hit.getADCCount();
245 double adcCalibrated = correction != 0 ? adcCount / scale / correction : 0;
246
247 // merge dEdx measurements on single wires; take active wires only
248 if (correction != 0) dedxWires[iwire].add(hit, iwire, currentLayer, celldx, adcCalibrated);
249
250 // debug output
251 if (dedxTrack) {
252 dedxTrack->m_pCDC = pCDC;
253 const auto& pocaMom = hit.getPOCAMomentum();
254 double pocaPhi = pocaMom.Phi();
255 double pocaTheta = pocaMom.Theta();
256 double pocaZ = hit.getPOCAOnWire().Z();
257 double hitCharge = adcTranslator.getCharge(adcCount, wireID, false, pocaZ, pocaPhi);
258 double driftDRealistic = tdcTranslator.getDriftLength(hit.getTDCCount(), wireID, 0, true, pocaZ, pocaPhi, pocaTheta);
259 double driftDRealisticRes = tdcTranslator.getDriftLengthResolution(driftDRealistic, wireID, true, pocaZ, pocaPhi, pocaTheta);
260 double cellDedx = adcCalibrated / celldx;
261
262 dedxTrack->addHit(wire, iwire, currentLayer, doca, docaRS, entAng, entAngRS,
263 adcCount, hit.getADCCount(), hitCharge, celldx * sinTheta, cellDedx, cellHeight, cellHalfWidth,
264 hit.getTDCCount(), driftDRealistic, driftDRealisticRes, wiregain, twodcor, onedcor,
265 hit.getFoundByTrackFinder(), hit.getWeightPionHypo(), hit.getWeightKaonHypo(), hit.getWeightProtonHypo());
266 }
267
268 } // end of loop over hits
269
270 // merge dEdx measurements in layers
271 std::map<int, DEDX> dedxLayers;
272 for (const auto& dedxWire : dedxWires) {
273 const auto& dedx = dedxWire.second;
274 dedxLayers[dedx.cLayer].add(dedx);
275 }
276
277 // push dEdx values to a vector
278 std::vector<double> dedxValues;
279 for (const auto& dedxLayer : dedxLayers) {
280 const auto& dedx = dedxLayer.second;
281 if (dedx.dx > 0 and dedx.dE > 0) {
282 dedxValues.push_back(dedx.dE / dedx.dx);
283 // debug output
284 if (dedxTrack) dedxTrack->addDedx(dedx.nhits, dedx.cWire, dedx.cLayer, dedx.dx, dedxValues.back());
285 }
286 }
287 if (dedxValues.empty()) continue;
288
289 // sort dEdx values
290 std::sort(dedxValues.begin(), dedxValues.end());
291
292 // calculate mean
293 double mean = 0;
294 for (auto x : dedxValues) mean += x;
295 mean /= dedxValues.size();
296
297 // calculate truncated mean and error
298 int lowEdgeTrunc = int(dedxValues.size() * m_removeLowest + 0.51);
299 int highEdgeTrunc = int(dedxValues.size() * (1 - m_removeHighest) + 0.51);
300 double truncatedMean = 0;
301 double sumOfSquares = 0;
302 int numValues = 0;
303 for (int i = lowEdgeTrunc; i < highEdgeTrunc; i++) {
304 double x = dedxValues[i];
305 truncatedMean += x;
306 sumOfSquares += x * x;
307 numValues++;
308 }
309 if (numValues > 0) {
310 truncatedMean /= numValues;
311 } else {
312 truncatedMean = mean;
313 numValues = dedxValues.size();
314 }
315 double truncatedError = numValues > 1 ? std::sqrt(sumOfSquares / numValues - truncatedMean * truncatedMean) / (numValues - 1) : 0;
316
317 // apply the saturation correction only to data (the so called "hadron correction")
318 double correctedMean = isData ? m_DBHadronCor->getCorrectedMean(truncatedMean, cosTheta) : truncatedMean;
319
320 // track level MC (e.g. replacing truncated mean with a generated one)
321 if (m_trackLevel and mcParticle) {
322 double mass = mcParticle->getMass();
323 if (mass > 0) {
324 double mcMean = m_DBMeanPars->getMean(pCDC / mass);
325 double mcSigma = m_DBSigmaPars->getSigma(mcMean, numValues, cosTheta, timeReso);
326 correctedMean = gRandom->Gaus(mcMean, mcSigma);
327 while (correctedMean < 0) correctedMean = gRandom->Gaus(mcMean, mcSigma);
328 // debug output
329 if (dedxTrack) dedxTrack->m_simDedx = correctedMean;
330 }
331 }
332
333 // calculate log likelihoods
334 double cdcLogL[Const::ChargedStable::c_SetSize] = {0};
335 for (const auto& chargedStable : Const::chargedStableSet) {
336 double betagamma = pCDC / chargedStable.getMass();
337 double predictedMean = m_DBMeanPars->getMean(betagamma);
338 double predictedSigma = m_DBSigmaPars->getSigma(predictedMean, numValues, cosTheta, timeReso);
339 if (predictedSigma <= 0) B2ERROR("Predicted sigma is not positive for PDG = " << chargedStable.getPDGCode());
340 double chi = (correctedMean - predictedMean) / predictedSigma;
341 int index = chargedStable.getIndex();
342 cdcLogL[index] = -0.5 * chi * chi;
343 // debug output
344 if (dedxTrack) {
345 dedxTrack->m_predmean[index] = predictedMean;
346 dedxTrack->m_predres[index] = predictedSigma;
347 dedxTrack->m_cdcChi[index] = chi;
348 dedxTrack->m_cdcLogl[index] = cdcLogL[index];
349 }
350 }
351
352 // save log likelihoods
353 auto* likelihoods = m_likelihoods.appendNew(cdcLogL);
354 track.addRelationTo(likelihoods);
355
356 // debug output
357 if (dedxTrack) {
358 double fullLength = 0;
359 for (const auto& dedxLayer : dedxLayers) fullLength += dedxLayer.second.dx;
360 dedxTrack->m_length = fullLength;
361 dedxTrack->m_dedxAvg = mean;
362 dedxTrack->m_dedxAvgTruncatedNoSat = truncatedMean;
363 dedxTrack->m_dedxAvgTruncatedErr = truncatedError;
364 dedxTrack->m_dedxAvgTruncated = correctedMean;
365 dedxTrack->m_lNHitsUsed = numValues;
366 track.addRelationTo(dedxTrack);
367 }
368
369 } // end of loop over tracks
370
371 }
372
374} // end Belle2 namespace
375
Class to store CDC hit information needed for dedx.
Definition CDCDedxHit.h:26
bool m_useBackHalfCurlers
whether to use the back half of curlers
StoreObjPtr< EventLevelTriggerTimeInfo > m_TTDInfo
injection time info
DBObjPtr< CDCDedxRunGain > m_DBRunGain
Run gain DB object.
DBObjPtr< CDCDedxMeanPars > m_DBMeanPars
dE/dx mean parameters
DBObjPtr< CDCDedxHadronCor > m_DBHadronCor
hadron saturation parameters
std::string m_dedxTracksName
name of collection of debug output
double m_removeHighest
portion of events with high dE/dx to discard
DBObjPtr< CDCDedxCosineEdge > m_DBCosEdgeCor
non-linearly ACD correction DB object
DBObjPtr< CDCDedxADCNonLinearity > m_DBNonlADC
non-linearly ACD correction DB object
DBObjPtr< CDCDedx1DCell > m_DB1DCell
1D correction DB object
StoreArray< CDCDedxLikelihood > m_likelihoods
collection of PID likelihoods
DBObjPtr< CDCDedx2DCell > m_DB2DCell
2D correction DB object
DBObjPtr< CDCDedxCosineCor > m_DBCosineCor
Electron saturation correction DB object.
StoreArray< Track > m_tracks
collection of tracks
int m_nLayerWires[9]
lookup table for number of wires per superlayer (indexed by superlayer)
DBObjPtr< CDCDedxSigmaPars > m_DBSigmaPars
dE/dx resolution parameters
DBObjPtr< CDCDedxWireGain > m_DBWireGains
Wire gain DB object.
bool m_trackLevel
whether to use track-level or hit-level MC
StoreArray< CDCDedxHit > m_hits
collection of hits
StoreArray< MCParticle > m_mcParticles
collection of MC particles
double m_removeLowest
portion of events with low dE/dx to discard
StoreArray< CDCDedxTrack > m_dedxTracks
collection of debug output
bool m_enableDebugOutput
option to write out debugging information to CDCDedxTracks
std::string m_likelihoodsName
name of collection of PID likelihoods
DBObjPtr< CDCDedxInjectionTime > m_DBInjectTime
time gain/reso DB object
DBObjPtr< CDCDedxScaleFactor > m_DBScaleFactor
Scale factor to make electrons ~1.
Debug output for CDCDedxPID module.
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
This class simply assumes a linear translation through (0,0)
float getCharge(unsigned short adcCount, const WireID &, bool, float, float)
just multiply with the conversion factor and return.
Translator mirroring the realistic Digitization.
double getDriftLength(unsigned short tdcCount, const WireID &wireID=WireID(), double timeOfFlightEstimator=0, bool leftRight=false, double z=0, double alpha=0, double theta=static_cast< double >(TMath::Pi()/2.), unsigned short adcCount=0) override
Get Drift length.
double getDriftLengthResolution(double driftLength, const WireID &wireID=WireID(), bool leftRight=false, double z=0, double alpha=0, double=static_cast< double >(TMath::Pi()/2.)) override
Get position resolution^2 corresponding to the drift length from getDriftLength of this class.
static const unsigned int c_SetSize
Number of elements (for use in array bounds etc.)
Definition Const.h:615
static const ParticleSet chargedStableSet
set of charged stable particles
Definition Const.h:618
static const ChargedStable pion
charged pion particle
Definition Const.h:661
A class to hold the geometry of a cell.
Definition LineHelper.h:186
double dx(const DedxPoint &poca, double entAng)
Calculate the path length through this cell for a track with a given DedxPoint Of Closest Approach (p...
Definition LineHelper.h:203
bool isValid()
Check if this is a valid calculation (number of intersections = 2)
Definition LineHelper.h:199
A collection of classes that are useful for making a simple path length correction to the dE/dx measu...
Definition LineHelper.h:29
bool isMC() const
Do we have generated, not real data?
static Environment & Instance()
Static method to get a reference to the Environment instance.
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
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
virtual void initialize() override
Initialize the module.
virtual void event() override
This method is called for each event.
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
Abstract base class for different kinds of events.
STL namespace.