195 for (
const auto& track :
m_tracks) {
196 std::shared_ptr<CDCDedxTrack> dedxTrack = std::make_shared<CDCDedxTrack>();
197 dedxTrack->m_track = mtrack++;
203 B2WARNING(
"No related fit for track ...");
217 dedxTrack->m_pdg = mcpart->
getPDG();
218 dedxTrack->m_mcmass = mcpart->
getMass();
220 dedxTrack->m_motherPDG = mother ? mother->
getPDG() : 0;
222 const ROOT::Math::XYZVector trueMomentum = mcpart->
getMomentum();
223 dedxTrack->m_pTrue = trueMomentum.R();
224 dedxTrack->m_cosThetaTrue = cos(trueMomentum.Theta());
227 dedxTrack->m_pdg = -999;
231 const ROOT::Math::XYZVector& trackMom = fitResult->
getMomentum();
232 dedxTrack->m_p = trackMom.R();
233 bool nomom = (dedxTrack->m_p != dedxTrack->m_p);
234 double costh = std::cos(std::atan(1 / fitResult->
getCotTheta()));
237 costh = cos(trackMom.Theta());
240 dedxTrack->m_cosTheta = costh;
241 dedxTrack->m_charge = charge;
243 double injring = -1.0, injtime = -1.0;
245 if (m_TTDInfo.
isValid() && m_TTDInfo->hasInjection()) {
246 injring = m_TTDInfo->isHER();
247 injtime = m_TTDInfo->getTimeSinceLastInjectionInMicroSeconds();
249 dedxTrack->m_injring = injring;
250 dedxTrack->m_injtime = injtime;
255 B2WARNING(
"No related track for this fit...");
261 B2ERROR(
"GFTrack is pruned, please run CDCDedxPID only on unpruned tracks! Skipping this track.");
278 if ((abs(costh + 0.860) < 0.010) || (abs(costh - 0.955) <= 0.005))isEdge =
true;
281 bool isvalidTime =
true;
282 if (injtime < 0 || injring < 0)isvalidTime =
false;
287 double layerdE = 0.0;
288 double layerdx = 0.0;
290 int nhitscombined = 0;
291 int wirelongesthit = 0;
292 double longesthit = 0;
297 const std::vector< genfit::AbsTrackRep* >& gftrackRepresentations = recoTrack->
getRepresentations();
299 for (std::vector< genfit::TrackPoint* >::const_iterator tp = gftrackPoints.begin();
300 tp != gftrackPoints.end(); ++tp) {
303 genfit::AbsMeasurement* aAbsMeasurementPtr = (*tp)->getRawMeasurement(0);
305 if (!cdcRecoHit)
continue;
307 if (!cdcHit)
continue;
311 const genfit::AbsFitterInfo* fi = (*tp)->getFitterInfo();
313 B2DEBUG(29,
"No fitter info, skipping...");
324 int foundByTrackFinder = 0;
329 double weightPionHypo = 0;
330 double weightProtHypo = 0;
331 double weightKaonHypo = 0;
333 for (std::vector<genfit::AbsTrackRep* >::const_iterator trep = gftrackRepresentations.begin();
334 trep != gftrackRepresentations.end(); ++trep) {
335 const int pdgCode = TMath::Abs((*trep)->getPDG());
340 const genfit::KalmanFitterInfo* kalmanFitterInfo = (*tp)->getKalmanFitterInfo(*trep);
341 if (kalmanFitterInfo == NULL) {
347 std::vector<double> weights = kalmanFitterInfo->getWeights();
348 const double maxWeight = weights[0] > weights[1] ? weights[0] : weights[1];
349 if (pdgCode ==
Const::pion.getPDGCode()) weightPionHypo = maxWeight;
350 else if (pdgCode ==
Const::kaon.getPDGCode()) weightKaonHypo = maxWeight;
351 else if (pdgCode ==
Const::proton.getPDGCode()) weightProtHypo = maxWeight;
356 int currentLayer = (superlayer == 0) ? layer : (8 + (superlayer - 1) * 6 + layer);
357 int nextLayer = currentLayer;
360 const int iwire = (superlayer == 0) ? 160 * layer + wire :
m_nLayerWires[superlayer - 1] + (160 + 32 *
361 (superlayer - 1)) * layer + wire;
364 const bool lastHit = (tp + 1 == gftrackPoints.end());
365 bool lastHitInCurrentLayer = lastHit;
368 genfit::AbsMeasurement* aAbsMeasurementPtrNext = (*(tp + 1))->getRawMeasurement(0);
371 if (!nextcdcRecoHit || !(cdcRecoHit->
getCDCHit()) || !((*(tp + 1))->getFitterInfo())) {
375 const int nextILayer = nextcdcHit->
getILayer();
377 nextLayer = (nextSuperlayer == 0) ? nextILayer : (8 + (nextSuperlayer - 1) * 6 + nextILayer);
378 lastHitInCurrentLayer = (nextLayer != currentLayer);
382 const ROOT::Math::XYZVector& wirePosF = cdcgeo.
wireForwardPosition(wireID, CDCGeometryPar::c_Aligned);
390 double topHeight = outer - wirePosF.Rho();
391 double bottomHeight = wirePosF.Rho() - inner;
392 double cellHeight = topHeight + bottomHeight;
393 double topHalfWidth = M_PI * outer / nWires;
394 double bottomHalfWidth = M_PI * inner / nWires;
395 double cellHalfWidth = M_PI * wirePosF.Rho() / nWires;
406 const genfit::MeasuredStateOnPlane& mop = fi->getFittedState();
409 ROOT::Math::XYZVector fittedPoca = ROOT::Math::XYZVector(mop.getPos());
410 const ROOT::Math::XYZVector& pocaMom = ROOT::Math::XYZVector(mop.getMom());
411 if (tp == gftrackPoints.begin() || cdcMom == 0) {
412 cdcMom = pocaMom.R();
413 dedxTrack->m_pCDC = cdcMom;
416 dedxTrack->m_p = cdcMom;
425 ROOT::Math::XYZVector pocaOnWire = ROOT::Math::XYZVector(mop.getPlane()->getO());
428 ROOT::Math::XYZVector B2WireDoca = fittedPoca - pocaOnWire;
431 double doca = B2WireDoca.Rho();
432 double phidiff = fittedPoca.Phi() - pocaOnWire.Phi();
435 if (phidiff > -3.1416 && (phidiff < 0 || phidiff > 3.1416)) doca *= -1;
438 const double px = pocaMom.X();
439 const double py = pocaMom.Y();
440 const double wx = pocaOnWire.X();
441 const double wy = pocaOnWire.Y();
442 const double cross = wx * py - wy * px;
443 const double dot = wx * px + wy * py;
444 double entAng = atan2(cross,
dot);
447 double cellR = 2 * cellHalfWidth / cellHeight;
449 if (std::abs(2 *
atan(1) - std::abs(entAng)) < 0.01)tana = 100 * (entAng / std::abs(entAng));
450 else tana = std::tan(entAng);
451 double docaRS = doca * std::sqrt((1 + cellR * cellR * tana * tana) / (1 + tana * tana));
452 double entAngRS = std::atan(tana / cellR);
457 && numMCParticles == 0) ?
m_DBNonlADC->getCorrectedADC(adcbaseCount, currentLayer) : adcbaseCount;
459 double hitCharge = translator.
getCharge(adcCount, wireID,
false, pocaOnWire.Z(), pocaMom.Phi());
463 double dadcCount = 1.0 * adcCount;
464 if (dedxTrack->m_scale != 0) dadcCount /= dedxTrack->m_scale;
467 double driftDRealistic = realistictdc.
getDriftLength(driftT, wireID, 0,
true, pocaOnWire.Z(), pocaMom.Phi(), pocaMom.Theta());
468 double driftDRealisticRes = realistictdc.
getDriftLengthResolution(driftDRealistic, wireID,
true, pocaOnWire.Z(), pocaMom.Phi(),
472 double celldx = c.dx(doca, entAng);
478 double normDocaRS = docaRS / cellHalfWidth;
480 && numMCParticles == 0) ?
m_DB2DCell->getMean(currentLayer, normDocaRS, entAngRS) : 1.0;
488 double correction = dedxTrack->m_runGain * dedxTrack->m_cosCor * dedxTrack->m_cosEdgeCor * dedxTrack->m_timeGain * wiregain *
494 if (correction != 0) dadcCount = dadcCount / correction;
497 double cellDedx = (dadcCount / celldx);
500 if (nomom) cellDedx *= std::sin(std::atan(1 / fitResult->
getCotTheta()));
501 else cellDedx *= std::sin(trackMom.Theta());
504 dedxTrack->addHit(wire, iwire, currentLayer, doca, docaRS, entAng, entAngRS,
505 adcCount, adcbaseCount, hitCharge, celldx, cellDedx, cellHeight, cellHalfWidth, driftT,
506 driftDRealistic, driftDRealisticRes, wiregain, twodcor, onedcor,
507 foundByTrackFinder, weightPionHypo, weightKaonHypo, weightProtHypo);
512 if (correction != 0) {
514 layerdE += dadcCount;
517 if (celldx > longesthit) {
519 wirelongesthit = iwire;
524 }
catch (genfit::Exception&) {
525 B2WARNING(
"Track: " << mtrack <<
": genfit::MeasuredStateOnPlane exception...");
530 if (lastHitInCurrentLayer) {
532 double totalDistance;
533 if (nomom) totalDistance = layerdx / std::sin(std::atan(1 / fitResult->
getCotTheta()));
534 else totalDistance = layerdx / std::sin(trackMom.Theta());
535 double layerDedx = layerdE / totalDistance;
539 dedxTrack->addDedx(nhitscombined, wirelongesthit, currentLayer, totalDistance, layerDedx);
543 saveLookupLogl(dedxTrack->m_cdcLogl, dedxTrack->m_pCDC, layerDedx);
560 if (dedxTrack->m_lDedx.empty()) {
561 B2DEBUG(29,
"Found track with no hits, ignoring.");
565 const int numDedx = dedxTrack->m_lDedx.size();
569 dedxTrack->m_lNHitsUsed = highEdgeTrunc - lowEdgeTrunc;
579 &(dedxTrack->m_dedxAvgTruncatedNoSat),
580 &(dedxTrack->m_dedxAvgTruncatedErr),
584 if (numMCParticles == 0)
585 dedxTrack->m_dedxAvgTruncated =
D2I(costh,
I2D(costh, 1.00) / 1.00 * dedxTrack->m_dedxAvgTruncatedNoSat);
586 else dedxTrack->m_dedxAvgTruncated = dedxTrack->m_dedxAvgTruncatedNoSat;
590 if (numMCParticles != 0 && dedxTrack->m_mcmass > 0 && dedxTrack->m_pTrue != 0) {
592 double mean =
getMean(dedxTrack->m_pCDC / dedxTrack->m_mcmass);
593 double sigma =
getSigma(mean, dedxTrack->m_lNHitsUsed,
594 dedxTrack->m_cosTheta, dedxTrack->m_timeReso);
595 dedxTrack->m_simDedx = gRandom->Gaus(mean, sigma);
596 while (dedxTrack->m_simDedx < 0)
597 dedxTrack->m_simDedx = gRandom->Gaus(mean, sigma);
599 dedxTrack->m_dedxAvgTruncated = dedxTrack->m_simDedx;
603 saveChiValue(dedxTrack->m_cdcChi, dedxTrack->m_predmean, dedxTrack->m_predres, dedxTrack->m_pCDC, dedxTrack->m_dedxAvgTruncated,
604 dedxTrack->m_cosTheta, dedxTrack->m_lNHitsUsed, dedxTrack->m_timeReso);
606 saveLookupLogl(dedxTrack->m_cdcLogl, dedxTrack->m_pCDC, dedxTrack->m_dedxAvgTruncated);
614 pidvalues[pdgIter.getIndex()] = -0.5 * dedxTrack->m_cdcChi[pdgIter.getIndex()] * dedxTrack->m_cdcChi[pdgIter.getIndex()];
618 pidvalues[pdgIter.getIndex()] = dedxTrack->m_cdcLogl[pdgIter.getIndex()];
623 track.addRelationTo(likelihoodObj);
628 track.addRelationTo(newCDCDedxTrack);