11 #include <reconstruction/modules/CDCDedxCorrection/CDCDedxCorrectionModule.h>
12 #include <reconstruction/dataobjects/CDCDedxTrack.h>
13 #include <reconstruction/dataobjects/DedxConstants.h>
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));
48 m_cdcDedxTracks.isRequired();
52 if (m_DBRunGain->getRunGain() == 0) {
53 B2WARNING(
"Run gain is zero for this run");
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);
63 unsigned int ncosbins = m_DBCosineCor->getSize();
64 for (
unsigned int i = 0; i < ncosbins; ++i) {
65 double gain = m_DBCosineCor->getMean(i);
67 B2ERROR(
"Cosine gain is zero for this bin " << i);
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();
88 for (
auto& dedxTrack : m_cdcDedxTracks) {
89 if (dedxTrack.size() == 0) {
90 B2WARNING(
"No good hits on this track...");
98 int nhits = dedxTrack.size();
99 double costh = dedxTrack.getCosTheta();
100 std::vector<double> newLayerHits;
101 double newLayerDe = 0, newLayerDx = 0;
103 if (costh < TMath::Cos(150.0 * TMath::DegToRad()))
continue;
104 if (costh > TMath::Cos(17.0 * TMath::DegToRad()))
continue;
106 for (
int i = 0; i < nhits; ++i) {
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);
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);
123 double newhitdedx = (m_relative) ? 1 / correction : 1.0;
124 StandardCorrection(jadcbase, jLayer, jWire, jNDocaRS, jEntaRS, jPath, costh, newhitdedx);
125 dedxTrack.setDedx(i, newhitdedx);
133 correction *= GetCorrection(jadcbase, jLayer, jWire, jNDocaRS, jEntaRS, costh);
134 if (!m_DBWireGains && dedxTrack.getWireGain(i) == 0)correction = 0;
137 correction = GetCorrection(jadcbase, jLayer, jWire, jNDocaRS, jEntaRS, costh);
141 if (correction != 0) {
142 newLayerDe += jadcbase / correction;
146 if (i + 1 < nhits && dedxTrack.getHitLayer(i + 1) == jLayer) {
149 if (newLayerDx != 0)newLayerHits.push_back(newLayerDe / newLayerDx * std::sqrt(1 - costh * costh));
156 calculateMeans(&(dedxTrack.m_dedxAvg),
157 &(dedxTrack.m_dedxAvgTruncatedNoSat),
158 &(dedxTrack.m_dedxAvgTruncatedErr),
161 dedxTrack.m_dedxAvgTruncated = dedxTrack.m_dedxAvgTruncatedNoSat;
162 HadronCorrection(costh, dedxTrack.m_dedxAvgTruncated);
169 B2INFO(
"CDCDedxCorrectionModule exiting...");
175 double gain = m_DBRunGain->getRunGain();
183 double gain = m_DBWireGains->getWireGain(wireID);
184 if (gain != 0) dedx = dedx / gain;
188 if (m_relative)dedx = 0;
195 double gain = (m_DB2DCell) ? m_DB2DCell->getMean(layer, doca, enta) : 1.0;
196 if (gain != 0) dedx = dedx / gain;
203 double gain = (m_DB1DCell) ? m_DB1DCell->getMean(layer, enta) : 1.0;
204 if (gain != 0) dedx = dedx / gain;
210 double coscor = m_DBCosineCor->getMean(costh);
211 if (coscor != 0) dedx = dedx / coscor;
218 double cosedgecor = m_DBCosEdgeCor->getMean(costh);
219 if (cosedgecor != 0) dedx = dedx / cosedgecor;
226 dedx = D2I(costheta, I2D(costheta, 1.00) / 1.00 * dedx);
230 double costheta,
double& dedx)
const
233 if (!m_relative && m_nonlADC)
234 adc = m_DBNonlADC->getCorrectedADC(adc, layer);
236 dedx *= adc * std::sqrt(1 - costheta * costheta) / length;
239 double scale = m_DBScaleFactor->getScaleFactor();
240 if (scale != 0) dedx = dedx / scale;
245 RunGainCorrection(dedx);
248 WireGainCorrection(wireID, dedx);
251 CosineCorrection(costheta, dedx);
254 if (!m_relative && m_cosineEdge && (costheta <= -0.850 || costheta >= 0.950))
255 CosineEdgeCorrection(costheta, dedx);
258 TwoDCorrection(layer, doca, enta, dedx);
261 OneDCorrection(layer, enta, dedx);
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);
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);
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!");
296 double chargeDensity = D / projection;
297 double numerator = 1 + m_hadronpars[0] * chargeDensity;
298 double denominator = 1 + m_hadronpars[1] * chargeDensity;
300 if (denominator == 0) {
301 B2WARNING(
"Something wrong with dE/dx hadron constants!");
305 double I = D * m_hadronpars[4] * numerator / denominator;
311 double absCosTheta = fabs(cosTheta);
312 double projection = pow(absCosTheta, m_hadronpars[3]) + m_hadronpars[2];
314 if (projection == 0 || m_hadronpars[4] == 0) {
315 B2WARNING(
"Something wrong with dE/dx hadron constants!");
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];
323 if (b == 0 && a == 0) {
324 B2WARNING(
"both a and b coefficiants for hadron correction are 0");
328 double D = (a != 0) ? (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a) : -c / b;
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;
333 B2WARNING(
"D is still less 0! just return uncorrectecd value");
343 const std::vector<double>& dedx)
const
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();
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();
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];
373 if (numValuesTrunc != 0) {
374 truncatedMeanTmp /= numValuesTrunc;
376 truncatedMeanTmp = meanTmp;
380 *truncatedMean = truncatedMeanTmp;
382 if (numValuesTrunc > 1) {
383 *truncatedMeanErr = sqrt(sumOfSquares /
double(numValuesTrunc) - truncatedMeanTmp * truncatedMeanTmp) / double(
386 *truncatedMeanErr = 0;