94 if (not
m_hits.isValid()) {
97 B2WARNING(
"StoreArray 'CDCDedxHits' does not exist, returning. Probably running on old cdst.");
99 B2WARNING(
"StoreArray 'CDCDedxHits' does not exist, returning. ...message will be suppressed now.");
117 double runGain = isData ?
m_DBRunGain->getRunGain() : 1.0;
126 B2ERROR(
"Scale factor from DB is zero! Will be set to one");
131 for (
const auto& track :
m_tracks) {
133 const auto* fitResult = track.getTrackFitResultWithClosestMass(
Const::pion);
135 B2WARNING(
"No related fit for track, skip it.");
140 const auto hits = track.getRelationsTo<
CDCDedxHit>();
141 if (hits.size() == 0)
continue;
144 const auto& trackMom = fitResult->getMomentum();
145 double theta = trackMom.Theta();
146 double cosTheta = std::cos(theta);
147 double sinTheta = std::sin(theta);
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;
155 const auto* mcParticle = isData ? nullptr : track.getRelated<
MCParticle>();
160 dedxTrack->m_track = track.getArrayIndex();
161 dedxTrack->m_charge = fitResult->getChargeSign();
162 dedxTrack->m_cosTheta = cosTheta;
163 dedxTrack->m_p = trackMom.R();
165 dedxTrack->m_injring =
m_TTDInfo->isHER();
166 dedxTrack->m_injtime =
m_TTDInfo->getTimeSinceLastInjectionInMicroSeconds();
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());
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;
188 std::map<int, DEDX> dedxWires;
189 for (
const auto& hit : hits) {
191 const auto& wireID = hit.getWireID();
192 int layer = wireID.getILayer();
193 int superlayer = wireID.getISuperLayer();
194 int currentLayer = (superlayer == 0) ? layer : (8 + (superlayer - 1) * 6 + layer);
196 lastLayer = currentLayer;
199 if (pCDC == 0) pCDC = hit.getPOCAMomentum().R();
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;
213 DedxPoint(bottomHalfWidth, -bottomHeight),
214 DedxPoint(-bottomHalfWidth, -bottomHeight));
217 double doca = hit.getSignedDOCAXY();
218 double entAng = hit.getEntranceAngle();
219 double celldx = cell.
dx(doca, entAng) / sinTheta;
220 if (not cell.
isValid())
continue;
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;
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);
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);
237 double onedcor = isData ?
m_DB1DCell->getMean(currentLayer, entAngRS) : 1.0;
238 double twodcor = isData ?
m_DB2DCell->getMean(currentLayer, normDocaRS, entAngRS) : 1.0;
241 double correction = runGain * cosCor * cosEdgeCor * timeGain * wiregain * twodcor * onedcor;
244 double adcCount = isData ?
m_DBNonlADC->getCorrectedADC(hit.getADCCount(), currentLayer) : hit.getADCCount();
245 double adcCalibrated = correction != 0 ? adcCount / scale / correction : 0;
248 if (correction != 0) dedxWires[iwire].add(hit, iwire, currentLayer, celldx, adcCalibrated);
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;
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());
271 std::map<int, DEDX> dedxLayers;
272 for (
const auto& dedxWire : dedxWires) {
273 const auto& dedx = dedxWire.second;
274 dedxLayers[dedx.cLayer].add(dedx);
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);
284 if (dedxTrack) dedxTrack->addDedx(dedx.nhits, dedx.cWire, dedx.cLayer, dedx.dx, dedxValues.back());
287 if (dedxValues.empty())
continue;
290 std::sort(dedxValues.begin(), dedxValues.end());
294 for (
auto x : dedxValues) mean += x;
295 mean /= dedxValues.size();
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;
303 for (
int i = lowEdgeTrunc; i < highEdgeTrunc; i++) {
304 double x = dedxValues[i];
306 sumOfSquares += x * x;
310 truncatedMean /= numValues;
312 truncatedMean = mean;
313 numValues = dedxValues.size();
315 double truncatedError = numValues > 1 ? std::sqrt(sumOfSquares / numValues - truncatedMean * truncatedMean) / (numValues - 1) : 0;
318 double correctedMean = isData ?
m_DBHadronCor->getCorrectedMean(truncatedMean, cosTheta) : truncatedMean;
322 double mass = mcParticle->getMass();
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);
329 if (dedxTrack) dedxTrack->m_simDedx = correctedMean;
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;
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];
354 track.addRelationTo(likelihoods);
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);