Belle II Software  release-05-02-19
CDCDedxCorrectionModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2012 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Jake Bennett
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <reconstruction/modules/CDCDedxCorrection/CDCDedxCorrectionModule.h>
12 #include <reconstruction/dataobjects/CDCDedxTrack.h>
13 #include <reconstruction/dataobjects/DedxConstants.h>
14 
15 #include <cmath>
16 #include <algorithm>
17 #include "TMath.h"
18 
19 using namespace Belle2;
20 using namespace Dedx;
21 
22 REG_MODULE(CDCDedxCorrection)
23 
25 {
26 
27  setDescription("Apply hit level corrections to the dE/dx measurements.");
28  addParam("relativeCorrections", m_relative, "If true, apply corrections relative to those used in production", true);
29  addParam("momentumCor", m_momCor, "Boolean to apply momentum correction", false);
30  addParam("momentumCorFromDB", m_useDBMomCor, "Boolean to apply momentum correction from DB", false);
31  addParam("scaleCor", m_scaleCor, "Boolean to apply scale correction", false);
32  addParam("cosineCor", m_cosineCor, "Boolean to apply cosine correction", false);
33  addParam("wireGain", m_wireGain, "Boolean to apply wire gains", false);
34  addParam("runGain", m_runGain, "Boolean to apply run gain", false);
35  addParam("twoDCell", m_twoDCell, "Boolean to apply 2D correction", false);
36  addParam("oneDCell", m_oneDCell, "Boolean to apply 1D correction", false);
37  addParam("cosineEdge", m_cosineEdge, "Boolean to apply cosine edge correction", true);
38  addParam("nonlADC", m_nonlADC, "Boolean to apply non-linear adc correction", true);
39  addParam("removeLowest", m_removeLowest, "portion of events with low dE/dx that should be discarded", double(0.05));
40  addParam("removeHighest", m_removeHighest, "portion of events with high dE/dx that should be discarded", double(0.25));
41 }
42 
44 
46 {
47  // register in datastore
48  m_cdcDedxTracks.isRequired();
49 
50  // make sure the calibration constants are reasonable
51  // run gains
52  if (m_DBRunGain->getRunGain() == 0) {
53  B2WARNING("Run gain is zero for this run");
54  }
55 
56  // wire gains
57  for (unsigned int i = 0; i < 14336; ++i) {
58  if (m_DBWireGains->getWireGain(i) == 0)
59  B2WARNING("Wire gain is zero for this wire: " << i);
60  }
61 
62  // cosine correction (store the bin edges for extrapolation)
63  unsigned int ncosbins = m_DBCosineCor->getSize();
64  for (unsigned int i = 0; i < ncosbins; ++i) {
65  double gain = m_DBCosineCor->getMean(i);
66  if (gain == 0)
67  B2ERROR("Cosine gain is zero for this bin " << i);
68  }
69 
70  // get the hadron correction parameters
71  if (!m_DBHadronCor || m_DBHadronCor->getSize() == 0) {
72  B2WARNING("No hadron correction parameters!");
73  for (int i = 0; i < 4; ++i)
74  m_hadronpars.push_back(0.0);
75  m_hadronpars.push_back(1.0);
76  } else m_hadronpars = m_DBHadronCor->getHadronPars();
77 }
78 
80 {
81 
82  // outputs
83  StoreArray<CDCDedxTrack> dedxArray;
84 
85  // **************************************************
86  // LOOP OVER EACH DEDX MEASUREMENT (TRACK LEVEL)
87  // **************************************************
88  for (auto& dedxTrack : m_cdcDedxTracks) {
89  if (dedxTrack.size() == 0) {
90  B2WARNING("No good hits on this track...");
91  continue;
92  }
93 
94  // **************************************************
95  // LOOP OVER EACH DEDX MEASUREMENT (HIT LEVEL)
96  // **************************************************
97  // hit level
98  int nhits = dedxTrack.size();
99  double costh = dedxTrack.getCosTheta();
100  std::vector<double> newLayerHits;
101  double newLayerDe = 0, newLayerDx = 0;
102 
103  if (costh < TMath::Cos(150.0 * TMath::DegToRad()))continue; //-0.866
104  if (costh > TMath::Cos(17.0 * TMath::DegToRad())) continue; //0.95
105 
106  for (int i = 0; i < nhits; ++i) {
107 
108  //pay attention to deadwire or gain uses
109  //getADCount is already corrected w/ non linear ADC payload
110  //getADCbasecount is now uncorrect ADC
111  int jadcbase = dedxTrack.getADCBaseCount(i);
112  double jLayer = dedxTrack.getHitLayer(i);
113  double jWire = dedxTrack.getWire(i);
114  double jNDocaRS = dedxTrack.getDocaRS(i) / dedxTrack.getCellHalfWidth(i);
115  double jEntaRS = dedxTrack.getEntaRS(i);
116  double jPath = dedxTrack.getPath(i);
117 
118  double correction = dedxTrack.m_scale * dedxTrack.m_runGain * dedxTrack.m_cosCor * dedxTrack.m_cosEdgeCor *
119  dedxTrack.getTwoDCorrection(i) * dedxTrack.getOneDCorrection(i) * dedxTrack.getNonLADCCorrection(i);
120  if (dedxTrack.getWireGain(i) > 0)correction *= dedxTrack.getWireGain(i); //also keep dead wire
121 
122  //Modify hit level dedx
123  double newhitdedx = (m_relative) ? 1 / correction : 1.0;
124  StandardCorrection(jadcbase, jLayer, jWire, jNDocaRS, jEntaRS, jPath, costh, newhitdedx);
125  dedxTrack.setDedx(i, newhitdedx);
126 
127  // do track level dedx and modifiy after loop over hits
128  // rel const -> upto 6 from calibrated GT and 2 are direct from dedx track (no rel cal for them now)
129  // abs const -> upto 6 from calibrated GT and 2 are direct from default GT
130  if (m_relative) {
131  //prewire gains need old tracks (old bad wire) and post need track new wg (new dead wire)
132  //get same base adc + rel correction factor
133  correction *= GetCorrection(jadcbase, jLayer, jWire, jNDocaRS, jEntaRS, costh);
134  if (!m_DBWireGains && dedxTrack.getWireGain(i) == 0)correction = 0;
135  } else {
136  //get modifed adc + abs correction factor
137  correction = GetCorrection(jadcbase, jLayer, jWire, jNDocaRS, jEntaRS, costh);
138  }
139 
140  // combine hits accross layers
141  if (correction != 0) {
142  newLayerDe += jadcbase / correction;
143  newLayerDx += jPath;
144  }
145 
146  if (i + 1 < nhits && dedxTrack.getHitLayer(i + 1) == jLayer) {
147  continue;
148  } else {
149  if (newLayerDx != 0)newLayerHits.push_back(newLayerDe / newLayerDx * std::sqrt(1 - costh * costh));
150  newLayerDe = 0;
151  newLayerDx = 0;
152  }
153  }
154 
155  // recalculate the truncated means
156  calculateMeans(&(dedxTrack.m_dedxAvg),
157  &(dedxTrack.m_dedxAvgTruncatedNoSat),
158  &(dedxTrack.m_dedxAvgTruncatedErr),
159  newLayerHits);
160 
161  dedxTrack.m_dedxAvgTruncated = dedxTrack.m_dedxAvgTruncatedNoSat;
162  HadronCorrection(costh, dedxTrack.m_dedxAvgTruncated);
163  } // end loop over tracks
164 }
165 
167 {
168 
169  B2INFO("CDCDedxCorrectionModule exiting...");
170 }
171 
173 {
174 
175  double gain = m_DBRunGain->getRunGain();
176  if (gain != 0) {
177  dedx = dedx / gain;
178  } else dedx = 0;
179 }
180 
181 void CDCDedxCorrectionModule::WireGainCorrection(int wireID, double& dedx) const
182 {
183  double gain = m_DBWireGains->getWireGain(wireID);
184  if (gain != 0) dedx = dedx / gain;
185  else {
186  //rel-abs method needs all wire for cal but w/ this method post calis (e.g.final RG)
187  //will also see all hitdedx but that is not an issue for track level calibration
188  if (m_relative)dedx = 0;
189  }
190 }
191 
192 void CDCDedxCorrectionModule::TwoDCorrection(int layer, double doca, double enta, double& dedx) const
193 {
194 
195  double gain = (m_DB2DCell) ? m_DB2DCell->getMean(layer, doca, enta) : 1.0;
196  if (gain != 0) dedx = dedx / gain;
197  else dedx = 0;
198 }
199 
200 
201 void CDCDedxCorrectionModule::OneDCorrection(int layer, double enta, double& dedx) const
202 {
203  double gain = (m_DB1DCell) ? m_DB1DCell->getMean(layer, enta) : 1.0;
204  if (gain != 0) dedx = dedx / gain;
205  else dedx = 0;
206 }
207 
208 void CDCDedxCorrectionModule::CosineCorrection(double costh, double& dedx) const
209 {
210  double coscor = m_DBCosineCor->getMean(costh);
211  if (coscor != 0) dedx = dedx / coscor;
212  else dedx = 0;
213 }
214 
215 void CDCDedxCorrectionModule::CosineEdgeCorrection(double costh, double& dedx) const
216 {
217 
218  double cosedgecor = m_DBCosEdgeCor->getMean(costh);
219  if (cosedgecor != 0) dedx = dedx / cosedgecor;
220  else dedx = 0;
221 }
222 
223 
224 void CDCDedxCorrectionModule::HadronCorrection(double costheta, double& dedx) const
225 {
226  dedx = D2I(costheta, I2D(costheta, 1.00) / 1.00 * dedx);
227 }
228 
229 void CDCDedxCorrectionModule::StandardCorrection(int adc, int layer, int wireID, double doca, double enta, double length,
230  double costheta, double& dedx) const
231 {
232 
233  if (!m_relative && m_nonlADC)
234  adc = m_DBNonlADC->getCorrectedADC(adc, layer);
235 
236  dedx *= adc * std::sqrt(1 - costheta * costheta) / length;
237 
238  if (m_scaleCor) {
239  double scale = m_DBScaleFactor->getScaleFactor();
240  if (scale != 0) dedx = dedx / scale;
241  else dedx = 0;
242  }
243 
244  if (m_runGain)
245  RunGainCorrection(dedx);
246 
247  if (m_wireGain)
248  WireGainCorrection(wireID, dedx);
249 
250  if (m_cosineCor)
251  CosineCorrection(costheta, dedx);
252 
253  //only if constants are abs are for specific costh
254  if (!m_relative && m_cosineEdge && (costheta <= -0.850 || costheta >= 0.950))
255  CosineEdgeCorrection(costheta, dedx);
256 
257  if (m_twoDCell)
258  TwoDCorrection(layer, doca, enta, dedx);
259 
260  if (m_oneDCell)
261  OneDCorrection(layer, enta, dedx);
262 
263 }
264 
265 
266 double CDCDedxCorrectionModule::GetCorrection(int& adc, int layer, int wireID, double doca, double enta, double costheta) const
267 {
268  double correction = 1.0;
269  if (m_scaleCor) correction *= m_DBScaleFactor->getScaleFactor();
270  if (m_runGain) correction *= m_DBRunGain->getRunGain();
271  if (m_wireGain) correction *= m_DBWireGains->getWireGain(wireID);
272  if (m_twoDCell) correction *= m_DB2DCell->getMean(layer, doca, enta);
273  if (m_oneDCell) correction *= m_DB1DCell->getMean(layer, enta);
274  if (m_cosineCor) correction *= m_DBCosineCor->getMean(costheta);
275 
276  //these last two are only for abs constant
277  if (!m_relative) {
278  if (m_nonlADC)adc = m_DBNonlADC->getCorrectedADC(adc, layer);
279  if (m_cosineEdge && (costheta <= -0.850 || costheta >= 0.950)) {
280  correction *= m_DBCosEdgeCor->getMean(costheta);
281  }
282  }
283  return correction;
284 }
285 
286 
287 double CDCDedxCorrectionModule::D2I(const double cosTheta, const double D) const
288 {
289  double absCosTheta = fabs(cosTheta);
290  double projection = pow(absCosTheta, m_hadronpars[3]) + m_hadronpars[2];
291  if (projection == 0) {
292  B2WARNING("Something wrong with dE/dx hadron constants!");
293  return D;
294  }
295 
296  double chargeDensity = D / projection;
297  double numerator = 1 + m_hadronpars[0] * chargeDensity;
298  double denominator = 1 + m_hadronpars[1] * chargeDensity;
299 
300  if (denominator == 0) {
301  B2WARNING("Something wrong with dE/dx hadron constants!");
302  return D;
303  }
304 
305  double I = D * m_hadronpars[4] * numerator / denominator;
306  return I;
307 }
308 
309 double CDCDedxCorrectionModule::I2D(const double cosTheta, const double I) const
310 {
311  double absCosTheta = fabs(cosTheta);
312  double projection = pow(absCosTheta, m_hadronpars[3]) + m_hadronpars[2];
313 
314  if (projection == 0 || m_hadronpars[4] == 0) {
315  B2WARNING("Something wrong with dE/dx hadron constants!");
316  return I;
317  }
318 
319  double a = m_hadronpars[0] / projection;
320  double b = 1 - m_hadronpars[1] / projection * (I / m_hadronpars[4]);
321  double c = -1.0 * I / m_hadronpars[4];
322 
323  if (b == 0 && a == 0) {
324  B2WARNING("both a and b coefficiants for hadron correction are 0");
325  return I;
326  }
327 
328  double D = (a != 0) ? (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a) : -c / b;
329  if (D < 0) {
330  B2WARNING("D is less 0! will try another solution");
331  D = (a != 0) ? (-b - sqrt(b * b + 4.0 * a * c)) / (2.0 * a) : -c / b;
332  if (D < 0) {
333  B2WARNING("D is still less 0! just return uncorrectecd value");
334  return I;
335  }
336  }
337 
338  return D;
339 }
340 
341 
342 void CDCDedxCorrectionModule::calculateMeans(double* mean, double* truncatedMean, double* truncatedMeanErr,
343  const std::vector<double>& dedx) const
344 {
345  // Calculate the truncated average by skipping the lowest & highest
346  // events in the array of dE/dx values
347  std::vector<double> sortedDedx = dedx;
348  std::sort(sortedDedx.begin(), sortedDedx.end());
349  sortedDedx.erase(std::remove(sortedDedx.begin(), sortedDedx.end(), 0), sortedDedx.end());
350  sortedDedx.shrink_to_fit();
351 
352  double truncatedMeanTmp = 0.0;
353  double meanTmp = 0.0;
354  double sumOfSquares = 0.0;
355  int numValuesTrunc = 0;
356  const int numDedx = sortedDedx.size();
357 
358  // add a factor of 0.51 here to make sure we are rounding appropriately...
359  const int lowEdgeTrunc = int(numDedx * m_removeLowest + 0.51);
360  const int highEdgeTrunc = int(numDedx * (1 - m_removeHighest) + 0.51);
361  for (int i = 0; i < numDedx; i++) {
362  meanTmp += sortedDedx[i];
363  if (i >= lowEdgeTrunc and i < highEdgeTrunc) {
364  truncatedMeanTmp += sortedDedx[i];
365  sumOfSquares += sortedDedx[i] * sortedDedx[i];
366  numValuesTrunc++;
367  }
368  }
369 
370  if (numDedx != 0) {
371  meanTmp /= numDedx;
372  }
373  if (numValuesTrunc != 0) {
374  truncatedMeanTmp /= numValuesTrunc;
375  } else {
376  truncatedMeanTmp = meanTmp;
377  }
378 
379  *mean = meanTmp;
380  *truncatedMean = truncatedMeanTmp;
381 
382  if (numValuesTrunc > 1) {
383  *truncatedMeanErr = sqrt(sumOfSquares / double(numValuesTrunc) - truncatedMeanTmp * truncatedMeanTmp) / double(
384  numValuesTrunc - 1);
385  } else {
386  *truncatedMeanErr = 0;
387  }
388 }
Belle2::CDCDedxCorrectionModule::calculateMeans
void calculateMeans(double *mean, double *truncatedMean, double *truncatedMeanErr, const std::vector< double > &dedx) const
Recalculate the dE/dx mean values after corrections.
Definition: CDCDedxCorrectionModule.cc:342
Belle2::CDCDedxCorrectionModule::CosineEdgeCorrection
void CosineEdgeCorrection(double costh, double &dedx) const
Perform the cosine edge correction.
Definition: CDCDedxCorrectionModule.cc:215
Belle2::CDCDedxCorrectionModule::HadronCorrection
void HadronCorrection(double costheta, double &dedx) const
Perform a hadron saturation correction.
Definition: CDCDedxCorrectionModule.cc:224
Belle2::CDCDedxCorrectionModule::TwoDCorrection
void TwoDCorrection(int layer, double doca, double enta, double &dedx) const
Perform a 2D correction.
Definition: CDCDedxCorrectionModule.cc:192
Belle2::CDCDedxCorrectionModule::StandardCorrection
void StandardCorrection(int adc, int layer, int wireID, double doca, double enta, double length, double costheta, double &dedx) const
Perform a standard set of corrections.
Definition: CDCDedxCorrectionModule.cc:229
Belle2::CDCDedxCorrectionModule::~CDCDedxCorrectionModule
virtual ~CDCDedxCorrectionModule()
Destructor.
Definition: CDCDedxCorrectionModule.cc:43
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::CDCDedxCorrectionModule::WireGainCorrection
void WireGainCorrection(int wireID, double &dedx) const
Perform a wire gain correction.
Definition: CDCDedxCorrectionModule.cc:181
Belle2::CDCDedxCorrectionModule::terminate
virtual void terminate() override
End of the event processing.
Definition: CDCDedxCorrectionModule.cc:166
Belle2::CDCDedxCorrectionModule::OneDCorrection
void OneDCorrection(int layer, double enta, double &dedx) const
Perform a wire gain correction.
Definition: CDCDedxCorrectionModule.cc:201
Belle2::CDCDedxCorrectionModule::event
virtual void event() override
This method is called for each event.
Definition: CDCDedxCorrectionModule.cc:79
Belle2::CDCDedxCorrectionModule::initialize
virtual void initialize() override
Initialize the Module.
Definition: CDCDedxCorrectionModule.cc:45
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::CDCDedxCorrectionModule::CosineCorrection
void CosineCorrection(double costheta, double &dedx) const
Perform the cosine correction.
Definition: CDCDedxCorrectionModule.cc:208
Belle2::CDCDedxCorrectionModule::GetCorrection
double GetCorrection(int &adc, int layer, int wireID, double doca, double enta, double costheta) const
Get the standard set of corrections.
Definition: CDCDedxCorrectionModule.cc:266
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDCDedxCorrectionModule::RunGainCorrection
void RunGainCorrection(double &dedx) const
Perform a run gain correction.
Definition: CDCDedxCorrectionModule.cc:172
Belle2::CDCDedxCorrectionModule
This module may be used to apply the corrections to dE/dx per the calibration constants.
Definition: CDCDedxCorrectionModule.h:53
Belle2::CDCDedxCorrectionModule::I2D
double I2D(const double cosTheta, const double I) const
Saturation correction: convert the actural ionization (I) to measured ionization (D)
Definition: CDCDedxCorrectionModule.cc:309
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::CDCDedxCorrectionModule::D2I
double D2I(const double cosTheta, const double D) const
Saturation correction: convert the measured ionization (D) to actual ionization (I)
Definition: CDCDedxCorrectionModule.cc:287