Belle II Software development
CDCDedxCorrectionModule.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/CDCDedxCorrection/CDCDedxCorrectionModule.h>
10
11#include <cdc/geometry/CDCGeometryPar.h>
12
13#include <TMath.h>
14
15using namespace Belle2;
16using namespace CDC;
17
18REG_MODULE(CDCDedxCorrection);
19
21{
22
23 setDescription("Apply hit level corrections to the dE/dx measurements.");
24 addParam("relativeCorrections", m_relative, "If true, apply corrections relative", true);
25 addParam("momentumCor", m_momCor, "Boolean to apply momentum correction", false);
26 addParam("momentumCorFromDB", m_useDBMomCor, "Boolean to apply momentum correction from DB", false);
27 addParam("scaleCor", m_scaleCor, "Boolean to apply scale correction", false);
28 addParam("cosineCor", m_cosineCor, "Boolean to apply cosine correction", false);
29 addParam("wireGain", m_wireGain, "Boolean to apply wire gains", false);
30 addParam("runGain", m_runGain, "Boolean to apply run gain", false);
31 addParam("timeGain", m_timeGain, "Boolean to apply time gain", false);
32 addParam("twoDCell", m_twoDCell, "Boolean to apply 2D correction", false);
33 addParam("oneDCell", m_oneDCell, "Boolean to apply 1D correction", false);
34 addParam("cosineEdge", m_cosineEdge, "Boolean to apply cosine edge correction", true);
35 addParam("nonlADC", m_nonlADC, "Boolean to apply non-linear adc correction", true);
36 addParam("removeLowest", m_removeLowest, "portion of events with low dE/dx", double(0.05));
37 addParam("removeHighest", m_removeHighest, "portion of events with high dE/dx", double(0.25));
38}
39
41
43{
44 // register in datastore
45 m_cdcDedxTracks.isRequired();
46
47 // make sure the calibration constants are reasonable
48 // run gains
49 if (m_DBRunGain->getRunGain() == 0) {
50 B2WARNING("Run gain is zero for this run");
51 }
52
53 // wire gains
54 for (unsigned int i = 0; i < c_nSenseWires; ++i) {
55 if (m_DBWireGains->getWireGain(i) == 0)
56 B2WARNING("Wire gain is zero for this wire: " << i);
57 }
58
59 // cosine correction (store the bin edges for extrapolation)
60 unsigned int ncosbins = m_DBCosineCor->getSize();
61 for (unsigned int i = 0; i < ncosbins; ++i) {
62 double gain = m_DBCosineCor->getMean(i);
63 if (gain == 0)
64 B2ERROR("Cosine gain is zero for this bin " << i);
65 }
66
67 // get the hadron correction parameters
68 if (!m_DBHadronCor || m_DBHadronCor->getSize() == 0) {
69 B2WARNING("No hadron correction parameters!");
70 for (int i = 0; i < 4; ++i)
71 m_hadronpars.push_back(0.0);
72 m_hadronpars.push_back(1.0);
73 } else m_hadronpars = m_DBHadronCor->getHadronPars();
74
75 int jwire = -1;
76 B2INFO("Creating CDCGeometryPar object");
78
79 for (unsigned int il = 0; il < c_maxNSenseLayers; il++) {
80 int activewires = 0;
81 m_lgainavg[il] = 0.0;
82
83 for (unsigned int iw = 0; iw < cdcgeo.nWiresInLayer(il); ++iw) {
84 jwire++;
85 if (m_DBWireGains->getWireGain(jwire) > 0.) {
86 //active wire only
87 m_lgainavg[il] += m_DBWireGains->getWireGain(jwire);
88 activewires++;
89 }
90 }
91 if (activewires > 0) m_lgainavg[il] /= activewires;
92 else m_lgainavg[il] = 1.0;
93 }
94}
95
96
98{
99
100 // **************************************************
101 // LOOP OVER EACH DEDX MEASUREMENT (TRACK LEVEL)
102 // **************************************************
103 for (auto& dedxTrack : m_cdcDedxTracks) {
104 if (dedxTrack.size() == 0) {
105 B2WARNING("No good hits on this track...");
106 continue;
107 }
108
109 // **************************************************
110 // LOOP OVER EACH DEDX MEASUREMENT (HIT LEVEL)
111 // **************************************************
112 // hit level
113 int nhits = dedxTrack.size();
114 double costh = dedxTrack.getCosTheta();
115 std::vector<double> newLayerHits;
116 double newLayerDe = 0, newLayerDx = 0;
117
118 if (costh < TMath::Cos(150.0 * TMath::DegToRad()))continue; //-0.866
119 if (costh > TMath::Cos(17.0 * TMath::DegToRad())) continue; //0.95
120
121
122 double injtime = dedxTrack.getInjectionTime();
123 double injring = dedxTrack.getInjectionRing();
124
125 //--------REMOVEBAL--------
126 //when CDST are reproduced with injection time
127 if (injtime == -1 || injring == -1) {
128 if (m_TTDInfo.isValid() && m_TTDInfo->hasInjection()) {
129 injring = m_TTDInfo->isHER();
130 injtime = m_TTDInfo->getTimeSinceLastInjectionInMicroSeconds();
131 }
132 }
133 dedxTrack.m_injring = injring;
134 dedxTrack.m_injtime = injtime;
135 //--------REMOVEBAL--------
136
137 for (int i = 0; i < nhits; ++i) {
138
139 //pay attention to deadwire or gain uses
140 //getADCount is already corrected w/ non linear ADC payload
141 //getADCbasecount is now incorrect ADC
142 int jadcbase = dedxTrack.getADCBaseCount(i);
143 double jLayer = dedxTrack.getHitLayer(i);
144 double jWire = dedxTrack.getWire(i);
145 double jNDocaRS = dedxTrack.getDocaRS(i) / dedxTrack.getCellHalfWidth(i);
146 double jEntaRS = dedxTrack.getEntaRS(i);
147 double jPath = dedxTrack.getPath(i);
148
149 double correction = dedxTrack.m_scale * dedxTrack.m_runGain * dedxTrack.m_timeGain * dedxTrack.m_cosCor * dedxTrack.m_cosEdgeCor *
150 dedxTrack.getTwoDCorrection(i) * dedxTrack.getOneDCorrection(i) * dedxTrack.getNonLADCCorrection(i);
151 if (dedxTrack.getWireGain(i) > 0)correction *= dedxTrack.getWireGain(i); //also keep dead wire
152
153 //Modify hit level dedx
154 double newhitdedx = (m_relative) ? 1 / correction : 1.0;
155 StandardCorrection(jadcbase, jLayer, jWire, jNDocaRS, jEntaRS, jPath, costh, injring, injtime, newhitdedx);
156 dedxTrack.setDedx(i, newhitdedx);
157
158 // do track level dedx and modify after loop over hits
159 // rel const -> upto 6 from calibrated GT and 2 are direct from dedx track (no rel cal for them now)
160 // abs const -> upto 6 from calibrated GT and 2 are direct from default GT
161 if (m_relative) {
162 //prewire gains need old tracks (old bad wire) and post need track new wg (new dead wire)
163 //get same base adc + rel correction factor
164 correction *= GetCorrection(jadcbase, jLayer, jWire, jNDocaRS, jEntaRS, costh, injring, injtime);
165 if (!m_DBWireGains && dedxTrack.getWireGain(i) == 0)correction = 0;
166 } else {
167 //get modified adc + abs correction factor
168 correction = GetCorrection(jadcbase, jLayer, jWire, jNDocaRS, jEntaRS, costh, injring, injtime);
169 }
170
171 // combine hits across layers
172 if (correction != 0) {
173 newLayerDe += jadcbase / correction;
174 newLayerDx += jPath;
175 }
176
177 if (i + 1 < nhits && dedxTrack.getHitLayer(i + 1) == jLayer) {
178 continue;
179 } else {
180 if (newLayerDx != 0)newLayerHits.push_back(newLayerDe / newLayerDx * std::sqrt(1 - costh * costh));
181 newLayerDe = 0;
182 newLayerDx = 0;
183 }
184 }
185
186 // recalculate the truncated means
187 calculateMeans(&(dedxTrack.m_dedxAvg),
188 &(dedxTrack.m_dedxAvgTruncatedNoSat),
189 &(dedxTrack.m_dedxAvgTruncatedErr),
190 newLayerHits);
191
192 dedxTrack.m_dedxAvgTruncated = dedxTrack.m_dedxAvgTruncatedNoSat;
193 HadronCorrection(costh, dedxTrack.m_dedxAvgTruncated);
194 } // end loop over tracks
195}
196
198{
199
200 B2INFO("CDCDedxCorrectionModule exiting...");
201}
202
204{
205 double gain = m_DBRunGain->getRunGain();
206 if (gain != 0) {
207 dedx = dedx / gain;
208 } else dedx = 0;
209}
210
211void CDCDedxCorrectionModule::TimeGainCorrection(double& dedx, double ring, double time) const
212{
213 double gain = m_DBInjectTime->getCorrection("mean", ring, time);
214 if (gain != 0) {
215 dedx = dedx / gain;
216 } else dedx = 0;
217}
218
219void CDCDedxCorrectionModule::WireGainCorrection(int wireID, double& dedx, int layer) const
220{
221 double gain = m_DBWireGains->getWireGain(wireID);
222 if (gain != 0) dedx = dedx / gain;
223 else {
224 //rel-abs method needs all wire for cal but w/ this method post calis (e.g.final RG)
225 //will also see all hitdedx but that is not an issue for track level calibration
226 if (m_relative)dedx = 0;
227 else {
228 if (m_lgainavg.at(layer) > 0)dedx = dedx / m_lgainavg.at(layer);
229 }
230 }
231}
232
233void CDCDedxCorrectionModule::TwoDCorrection(int layer, double doca, double enta, double& dedx) const
234{
235
236 double gain = (m_DB2DCell) ? m_DB2DCell->getMean(layer, doca, enta) : 1.0;
237 if (gain != 0) dedx = dedx / gain;
238 else dedx = 0;
239}
240
241
242void CDCDedxCorrectionModule::OneDCorrection(int layer, double enta, double& dedx) const
243{
244 double gain = (m_DB1DCell) ? m_DB1DCell->getMean(layer, enta) : 1.0;
245 if (gain != 0) dedx = dedx / gain;
246 else dedx = 0;
247}
248
249void CDCDedxCorrectionModule::CosineCorrection(double costh, double& dedx) const
250{
251 double coscor = m_DBCosineCor->getMean(costh);
252 if (coscor != 0) dedx = dedx / coscor;
253 else dedx = 0;
254}
255
256void CDCDedxCorrectionModule::CosineEdgeCorrection(double costh, double& dedx) const
257{
258
259 double cosedgecor = m_DBCosEdgeCor->getMean(costh);
260 if (cosedgecor != 0) dedx = dedx / cosedgecor;
261 else dedx = 0;
262}
263
264
265void CDCDedxCorrectionModule::HadronCorrection(double costheta, double& dedx) const
266{
267 dedx = D2I(costheta, I2D(costheta, 1.00) / 1.00 * dedx);
268}
269
270void CDCDedxCorrectionModule::StandardCorrection(int adc, int layer, int wireID, double doca, double enta, double length,
271 double costheta, double ring, double time, double& dedx) const
272{
273
274 if (!m_relative && m_nonlADC)
275 adc = m_DBNonlADC->getCorrectedADC(adc, layer);
276
277 dedx *= adc * std::sqrt(1 - costheta * costheta) / length;
278
279 if (m_scaleCor) {
280 double scale = m_DBScaleFactor->getScaleFactor();
281 if (scale != 0) dedx = dedx / scale;
282 else dedx = 0;
283 }
284
285 if (m_runGain)
286 RunGainCorrection(dedx);
287
288 if (m_timeGain)
289 TimeGainCorrection(dedx, ring, time);
290
291 if (m_wireGain)
292 WireGainCorrection(wireID, dedx, layer);
293
294 if (m_cosineCor)
295 CosineCorrection(costheta, dedx);
296
297 //only if constants are abs are for specific costh
298 if (!m_relative && m_cosineEdge && (costheta <= -0.850 || costheta >= 0.950))
299 CosineEdgeCorrection(costheta, dedx);
300
301 if (m_twoDCell)
302 TwoDCorrection(layer, doca, enta, dedx);
303
304 if (m_oneDCell)
305 OneDCorrection(layer, enta, dedx);
306}
307
308
309double CDCDedxCorrectionModule::GetCorrection(int& adc, int layer, int wireID, double doca, double enta, double costheta,
310 double ring, double time) const
311{
312 double correction = 1.0;
313 if (m_scaleCor) correction *= m_DBScaleFactor->getScaleFactor();
314 if (m_runGain) correction *= m_DBRunGain->getRunGain();
315 if (m_timeGain) correction *= m_DBInjectTime->getCorrection("mean", ring, time);
316 if (m_wireGain) correction *= m_DBWireGains->getWireGain(wireID);
317 if (m_twoDCell) correction *= m_DB2DCell->getMean(layer, doca, enta);
318 if (m_oneDCell) correction *= m_DB1DCell->getMean(layer, enta);
319 if (m_cosineCor) correction *= m_DBCosineCor->getMean(costheta);
320
321 //these last two are only for abs constant
322 if (!m_relative) {
323 if (m_nonlADC)adc = m_DBNonlADC->getCorrectedADC(adc, layer);
324 if (m_cosineEdge && (costheta <= -0.850 || costheta >= 0.950)) {
325 correction *= m_DBCosEdgeCor->getMean(costheta);
326 }
327 }
328 return correction;
329}
330
331
332double CDCDedxCorrectionModule::D2I(const double cosTheta, const double D) const
333{
334 double absCosTheta = fabs(cosTheta);
335 double projection = pow(absCosTheta, m_hadronpars[3]) + m_hadronpars[2];
336 if (projection == 0) {
337 B2WARNING("Something wrong with dE/dx hadron constants!");
338 return D;
339 }
340
341 double chargeDensity = D / projection;
342 double numerator = 1 + m_hadronpars[0] * chargeDensity;
343 double denominator = 1 + m_hadronpars[1] * chargeDensity;
344
345 if (denominator == 0) {
346 B2WARNING("Something wrong with dE/dx hadron constants!");
347 return D;
348 }
349
350 double I = D * m_hadronpars[4] * numerator / denominator;
351 return I;
352}
353
354double CDCDedxCorrectionModule::I2D(const double cosTheta, const double I) const
355{
356 double absCosTheta = fabs(cosTheta);
357 double projection = pow(absCosTheta, m_hadronpars[3]) + m_hadronpars[2];
358
359 if (projection == 0 || m_hadronpars[4] == 0) {
360 B2WARNING("Something wrong with dE/dx hadron constants!");
361 return I;
362 }
363
364 double a = m_hadronpars[0] / projection;
365 double b = 1 - m_hadronpars[1] / projection * (I / m_hadronpars[4]);
366 double c = -1.0 * I / m_hadronpars[4];
367
368 if (b == 0 && a == 0) {
369 B2WARNING("both a and b coefficiants for hadron correction are 0");
370 return I;
371 }
372
373 double discr = b * b - 4.0 * a * c;
374 if (discr < 0) {
375 B2WARNING("negative discriminant; return uncorrectecd value");
376 return I;
377 }
378
379 double D = (a != 0) ? (-b + sqrt(discr)) / (2.0 * a) : -c / b;
380 if (D < 0) {
381 B2WARNING("D is less 0! will try another solution");
382 D = (a != 0) ? (-b - sqrt(discr)) / (2.0 * a) : -c / b;
383 if (D < 0) {
384 B2WARNING("D is still less 0! just return uncorrectecd value");
385 return I;
386 }
387 }
388
389 return D;
390}
391
392
393void CDCDedxCorrectionModule::calculateMeans(double* mean, double* truncatedMean, double* truncatedMeanErr,
394 const std::vector<double>& dedx) const
395{
396 // Calculate the truncated average by skipping the lowest & highest
397 // events in the array of dE/dx values
398 std::vector<double> sortedDedx = dedx;
399 std::sort(sortedDedx.begin(), sortedDedx.end());
400 sortedDedx.erase(std::remove(sortedDedx.begin(), sortedDedx.end(), 0), sortedDedx.end());
401 sortedDedx.shrink_to_fit();
402
403 double truncatedMeanTmp = 0.0;
404 double meanTmp = 0.0;
405 double sumOfSquares = 0.0;
406 int numValuesTrunc = 0;
407 const int numDedx = sortedDedx.size();
408
409 // add a factor of 0.51 here to make sure we are rounding appropriately...
410 const int lowEdgeTrunc = int(numDedx * m_removeLowest + 0.51);
411 const int highEdgeTrunc = int(numDedx * (1 - m_removeHighest) + 0.51);
412 for (int i = 0; i < numDedx; i++) {
413 meanTmp += sortedDedx[i];
414 if (i >= lowEdgeTrunc and i < highEdgeTrunc) {
415 truncatedMeanTmp += sortedDedx[i];
416 sumOfSquares += sortedDedx[i] * sortedDedx[i];
417 numValuesTrunc++;
418 }
419 }
420
421 if (numDedx != 0) {
422 meanTmp /= numDedx;
423 }
424 if (numValuesTrunc != 0) {
425 truncatedMeanTmp /= numValuesTrunc;
426 } else {
427 truncatedMeanTmp = meanTmp;
428 }
429
430 *mean = meanTmp;
431 *truncatedMean = truncatedMeanTmp;
432
433 if (numValuesTrunc > 1) {
434 *truncatedMeanErr = sqrt(sumOfSquares / double(numValuesTrunc) - truncatedMeanTmp * truncatedMeanTmp) / double(
435 numValuesTrunc - 1);
436 } else {
437 *truncatedMeanErr = 0;
438 }
439}
bool m_twoDCell
boolean to apply 2D correction
double D2I(const double cosTheta, const double D) const
Saturation correction: convert the measured ionization (D) to actual ionization (I)
bool m_runGain
boolean to apply run gains
StoreObjPtr< EventLevelTriggerTimeInfo > m_TTDInfo
Store Object Ptr: EventLevelTriggerTimeInfo.
std::vector< double > m_hadronpars
hadron saturation parameters
bool m_cosineEdge
boolean to apply cosine edge
DBObjPtr< CDCDedxRunGain > m_DBRunGain
Run gain DB object.
bool m_useDBMomCor
boolean to apply momentum correction from DB
void TwoDCorrection(int layer, double doca, double enta, double &dedx) const
Perform a 2D correction.
bool m_momCor
boolean to apply momentum correction
virtual void initialize() override
Initialize the Module.
bool m_scaleCor
boolean to apply scale factor
DBObjPtr< CDCDedxHadronCor > m_DBHadronCor
hadron saturation parameters
virtual void event() override
This method is called for each event.
double m_removeHighest
upper bound for truncated mean
DBObjPtr< CDCDedxCosineEdge > m_DBCosEdgeCor
cosine edge calibration
double I2D(const double cosTheta, const double I) const
Saturation correction: convert the actual ionization (I) to measured ionization (D)
std::array< double, 56 > m_lgainavg
average calibration factor for the layer
void OneDCorrection(int layer, double enta, double &dedx) const
Perform a wire gain correction.
DBObjPtr< CDCDedxADCNonLinearity > m_DBNonlADC
hadron saturation non linearity
DBObjPtr< CDCDedx1DCell > m_DB1DCell
1D correction DB object
virtual void terminate() override
End of the event processing.
bool m_oneDCell
boolean to apply 1D correction
void HadronCorrection(double costheta, double &dedx) const
Perform a hadron saturation correction.
void calculateMeans(double *mean, double *truncatedMean, double *truncatedMeanErr, const std::vector< double > &dedx) const
Recalculate the dE/dx mean values after corrections.
void StandardCorrection(int adc, int layer, int wireID, double doca, double enta, double length, double costheta, double ring, double time, double &dedx) const
Perform a standard set of corrections.
DBObjPtr< CDCDedx2DCell > m_DB2DCell
2D correction DB object
bool m_nonlADC
boolean to apply non linear ADC
bool m_relative
boolean to apply relative or absolute correction
DBObjPtr< CDCDedxCosineCor > m_DBCosineCor
Electron saturation correction DB object.
StoreArray< CDCDedxTrack > m_cdcDedxTracks
Store array: CDCDedxTrack.
bool m_timeGain
boolean to apply injection time gains
void CosineEdgeCorrection(double costh, double &dedx) const
Perform the cosine edge correction.
DBObjPtr< CDCDedxWireGain > m_DBWireGains
Wire gain DB object.
void CosineCorrection(double costheta, double &dedx) const
Perform the cosine correction.
void WireGainCorrection(int wireID, double &dedx, int layer) const
Perform a wire gain correction.
bool m_cosineCor
boolean to apply cosine correction
void RunGainCorrection(double &dedx) const
Perform a run gain correction.
double GetCorrection(int &adc, int layer, int wireID, double doca, double enta, double costheta, double ring, double time) const
Get the standard set of corrections.
double m_removeLowest
lower bound for truncated mean
DBObjPtr< CDCDedxInjectionTime > m_DBInjectTime
time gain/reso DB object
DBObjPtr< CDCDedxScaleFactor > m_DBScaleFactor
Scale factor to make electrons ~1.
void TimeGainCorrection(double &dedx, double ring, double time) const
Perform a injection time gain correction.
bool m_wireGain
boolean to apply wire gains
The Class for CDC Geometry Parameters.
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
Module()
Constructor.
Definition Module.cc:30
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
Abstract base class for different kinds of events.