diff --git a/ALICE3/Core/FastTracker.cxx b/ALICE3/Core/FastTracker.cxx index a74858cfe02..be4ef7d2a12 100644 --- a/ALICE3/Core/FastTracker.cxx +++ b/ALICE3/Core/FastTracker.cxx @@ -288,7 +288,7 @@ float FastTracker::ProbGoodChiSqHit(float radius, float searchRadiusRPhi, float // function to provide a reconstructed track from a perfect input track // returns number of intercepts (generic for now) -int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack, const float nch) +int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack, const float nch, const float maxRadius) { dNdEtaCent = nch; // set the number of charged particles per unit rapidity hits.clear(); @@ -335,6 +335,14 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa continue; // this layer should not be attempted, but go ahead } + if (layers[il].getRadius() > maxRadius) { + if (lastLayerReached == -1) { + // This means that we didn't reach the first layer + return -9; + } + break; // could not reach + } + // check if layer is reached float targetX = 1e+3; inputTrack.getXatLabR(layers[il].getRadius(), targetX, magneticField); @@ -367,6 +375,7 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa break; } } + if (std::abs(inputTrack.getZ()) > layers[il].getZ() && mApplyZacceptance) { break; // out of acceptance bounds } @@ -405,8 +414,9 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa static constexpr float kLargeErr2Dir = 0.7 * 0.7; static constexpr float kLargeErr2PtI = 30.5 * 30.5; std::array largeCov = {0.}; - for (int ic = o2::track::kCovMatSize; ic--;) + for (int ic = o2::track::kCovMatSize; ic--;) { largeCov[ic] = 0.; + } largeCov[o2::track::CovLabels::kSigY2] = largeCov[o2::track::CovLabels::kSigZ2] = kLargeErr2Coord; largeCov[o2::track::CovLabels::kSigSnp2] = largeCov[o2::track::CovLabels::kSigTgl2] = kLargeErr2Dir; largeCov[o2::track::CovLabels::kSigQ2Pt2] = kLargeErr2PtI * trPars[o2::track::ParLabels::kQ2Pt] * trPars[o2::track::ParLabels::kQ2Pt]; @@ -442,8 +452,10 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa std::cos(alpha) * spacePoint[0] + std::sin(alpha) * spacePoint[1], -std::sin(alpha) * spacePoint[0] + std::cos(alpha) * spacePoint[1], spacePoint[2]}; - if (!inwardTrack.propagateTo(xyz1[0], magneticField)) + + if (!inwardTrack.propagateTo(xyz1[0], magneticField)) { continue; + } if (!layers[il].isInert()) { // only update covm for tracker hits const o2::track::TrackParametrization::dim2_t hitpoint = { @@ -474,13 +486,14 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa } } - if (layers[il].isSilicon()) + if (layers[il].isSilicon()) { nSiliconPoints++; // count silicon hits - if (layers[il].isGas()) + } + if (layers[il].isGas()) { nGasPoints++; // count TPC/gas hits + } hits.push_back(thisHit); - if (!layers[il].isInert()) { // good hit probability calculation float sigYCmb = o2::math_utils::sqrt(inwardTrack.getSigmaY2() + layers[il].getResolutionRPhi() * layers[il].getResolutionRPhi()); float sigZCmb = o2::math_utils::sqrt(inwardTrack.getSigmaZ2() + layers[il].getResolutionZ() * layers[il].getResolutionZ()); @@ -502,21 +515,24 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa } // only attempt to continue if intercepts are at least four - if (nIntercepts < 4) + if (nIntercepts < 4) { return nIntercepts; + } // generate efficiency float eff = 1.; for (size_t i = 0; i < layers.size(); i++) { float iGoodHit = goodHitProbability[i]; - if (iGoodHit <= 0) + if (iGoodHit <= 0) { continue; + } eff *= iGoodHit; } if (mApplyEffCorrection) { - if (gRandom->Uniform() > eff) + if (gRandom->Uniform() > eff) { return -8; + } } outputTrack.setCov(inwardTrack.getCov()); @@ -524,8 +540,9 @@ int FastTracker::FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackPa // Use covariance matrix based smearing std::array covMat = {0.}; - for (int ii = 0; ii < o2::track::kCovMatSize; ii++) + for (int ii = 0; ii < o2::track::kCovMatSize; ii++) { covMat[ii] = outputTrack.getCov()[ii]; + } TMatrixDSym m(5); double fcovm[5][5]; // double precision is needed for regularisation diff --git a/ALICE3/Core/FastTracker.h b/ALICE3/Core/FastTracker.h index 5b8fa150c4c..d361e261603 100644 --- a/ALICE3/Core/FastTracker.h +++ b/ALICE3/Core/FastTracker.h @@ -92,7 +92,7 @@ class FastTracker * @param nch Charged particle multiplicity (used for hit density calculations). * @return int i.e. number of intercepts (implementation-defined). */ - int FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack, const float nch); + int FastTrack(o2::track::TrackParCov inputTrack, o2::track::TrackParCov& outputTrack, const float nch, const float maxRadius = 100.f); // For efficiency calculation float Dist(float z, float radius); diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 814b55623fe..14ab69f0dab 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -181,6 +181,19 @@ struct OnTheFlyTracker { Configurable doV0QA{"doV0QA", false, "QA plots for when treating V0"}; } v0DecaySettings; + struct : ConfigurableGroup { + std::string prefix = "cfgFitter"; + Configurable propagateToPCA{"propagateToPCA", false, "create tracks version propagated to PCA"}; + Configurable maxR{"maxR", 200., "reject PCA's above this radius"}; + Configurable minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"}; + Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; + Configurable maxDZIni{"maxDZIni", 1e9, "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; + Configurable maxDXYIni{"maxDXYIni", 4, "reject (if>0) PCA candidate if tracks DXY exceeds threshold"}; + Configurable maxVtxChi2{"maxVtxChi2", 1e9, "reject (if>0) vtx. chi2 above this value"}; + Configurable useAbsDCA{"useAbsDCA", true, "Minimise abs. distance rather than chi2"}; + Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; + } cfgFitter; + using PVertex = o2::dataformats::PrimaryVertex; // for secondary vertex finding @@ -465,6 +478,7 @@ struct OnTheFlyTracker { getHist(TH1, histPath + "hFastTrackerQA")->GetXaxis()->SetBinLabel(6, "multiple scattering"); getHist(TH1, histPath + "hFastTrackerQA")->GetXaxis()->SetBinLabel(7, "energy loss"); getHist(TH1, histPath + "hFastTrackerQA")->GetXaxis()->SetBinLabel(8, "efficiency"); + getHist(TH1, histPath + "hFastTrackerQA")->GetXaxis()->SetBinLabel(9, "no layers hit"); } } @@ -491,6 +505,27 @@ struct OnTheFlyTracker { hCovMatOK->GetXaxis()->SetBinLabel(1, "Not OK"); hCovMatOK->GetXaxis()->SetBinLabel(2, "OK"); + auto hFitterStatusCode = histos.add("hFitterStatusCode", "hFitterStatusCode", kTH1D, {{15, -0.5, 14.5}}); + hFitterStatusCode->GetXaxis()->SetBinLabel(1, "None"); // no status set (should not be possible!) + + /* Good Conditions */ + hFitterStatusCode->GetXaxis()->SetBinLabel(2, "Converged"); // fit converged + hFitterStatusCode->GetXaxis()->SetBinLabel(3, "MaxIter"); // max iterations reached before fit convergence + + /* Error Conditions */ + hFitterStatusCode->GetXaxis()->SetBinLabel(4, "NoCrossing"); // no reasaonable crossing was found + hFitterStatusCode->GetXaxis()->SetBinLabel(5, "RejRadius"); // radius of crossing was not acceptable + hFitterStatusCode->GetXaxis()->SetBinLabel(6, "RejTrackX"); // one candidate track x was below the mimimum required radius + hFitterStatusCode->GetXaxis()->SetBinLabel(7, "RejTrackRoughZ"); // rejected by rough cut on tracks Z difference + hFitterStatusCode->GetXaxis()->SetBinLabel(8, "RejChi2Max"); // rejected by maximum chi2 cut + hFitterStatusCode->GetXaxis()->SetBinLabel(9, "FailProp"); // propagation of at least prong to PCA failed + hFitterStatusCode->GetXaxis()->SetBinLabel(10, "FailInvCov"); // inversion of cov.-matrix failed + hFitterStatusCode->GetXaxis()->SetBinLabel(11, "FailInvWeight"); // inversion of Ti weight matrix failed + hFitterStatusCode->GetXaxis()->SetBinLabel(12, "FailInv2ndDeriv"); // inversion of 2nd derivatives failed + hFitterStatusCode->GetXaxis()->SetBinLabel(13, "FailCorrTracks"); // correction of tracks to updated x failed + hFitterStatusCode->GetXaxis()->SetBinLabel(14, "FailCloserAlt"); // alternative PCA is closer + hFitterStatusCode->GetXaxis()->SetBinLabel(15, "NStatusesDefined"); + if (doExtraQA) { histos.add("h2dVerticesVsContributors", "h2dVerticesVsContributors;Multiplicity;N vertices", kTH2F, {axes.axisMultiplicity, axes.axisNVertices}); histos.add("h1dVerticesNotReco", "h1dVerticesNotReco;Multiplicity;Vertices Not Reco", kTH1F, {axes.axisMultiplicity}); @@ -515,6 +550,7 @@ struct OnTheFlyTracker { h->GetXaxis()->SetBinLabel(6, "multiple scattering"); h->GetXaxis()->SetBinLabel(7, "energy loss"); h->GetXaxis()->SetBinLabel(8, "efficiency"); + h->GetXaxis()->SetBinLabel(8, "no layers hit"); histPointers.insert({v0histPath + "hFastTrackerQA", h}); // K0s insertHist(v0histPath + "K0/hGen", "hGen", kTH2F, {axes.axisDecayRadius, axes.axisMomentum}); @@ -583,15 +619,16 @@ struct OnTheFlyTracker { o2::vertexing::PVertexerParams::Instance().printKeyValues(); // initialize O2 2-prong fitter - fitter.setPropagateToPCA(true); - fitter.setMaxR(200.); - fitter.setMinParamChange(1e-3); - fitter.setMinRelChi2Change(0.9); - fitter.setMaxDZIni(1e9); - fitter.setMaxDXYIni(4); - fitter.setMaxChi2(1e9); - fitter.setUseAbsDCA(true); - fitter.setWeightedFinalPCA(false); + fitter.setPropagateToPCA(cfgFitter.propagateToPCA); + fitter.setMaxR(cfgFitter.maxR); + fitter.setMinParamChange(cfgFitter.minParamChange); + fitter.setMinRelChi2Change(cfgFitter.minRelChi2Change); + fitter.setMaxDZIni(cfgFitter.maxDZIni); + fitter.setMaxDXYIni(cfgFitter.maxDXYIni); + fitter.setMaxChi2(cfgFitter.maxVtxChi2); + fitter.setUseAbsDCA(cfgFitter.useAbsDCA); + fitter.setWeightedFinalPCA(cfgFitter.useWeightedFinalPCA); + fitter.setBz(mMagneticField); fitter.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrNONE); // such a light detector here // Set seed for TGenPhaseSpace @@ -987,6 +1024,12 @@ struct OnTheFlyTracker { if (nCand == 0) { dcaFitterOK_V0 = false; } + + fitter.propagateTracksToVertex(); + if (!fitter.isPropagateTracksToVertexDone()) { + dcaFitterOK_V0 = false; + } + // V0 found successfully if (dcaFitterOK_V0) { if (cascadeDecaySettings.doXiQA) { @@ -1026,6 +1069,7 @@ struct OnTheFlyTracker { covV[MomInd[i]] = 1e-6; covV[i] = 1e-6; } + o2::track::TrackParCov v0Track = o2::track::TrackParCov( {pos[0], pos[1], pos[2]}, {posP[0] + negP[0], posP[1] + negP[1], posP[2] + negP[2]}, @@ -1046,6 +1090,13 @@ struct OnTheFlyTracker { dcaFitterOK_Cascade = false; } + fitter.propagateTracksToVertex(); + if (!fitter.isPropagateTracksToVertexDone()) { + dcaFitterOK_Cascade = false; + } + + const u_int8_t fitterStatusCode = fitter.getFitStatus(); + histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); // Cascade found successfully if (dcaFitterOK_Cascade) { if (cascadeDecaySettings.doXiQA) { @@ -1098,8 +1149,9 @@ struct OnTheFlyTracker { // find perfect intercept XYZ float targetX = 1e+3; trackParCov.getXatLabR(layer.getRadius(), targetX, mMagneticField); - if (targetX > 999) + if (targetX > 999) { continue; // failed to find intercept + } if (!trackParCov.propagateTo(targetX, mMagneticField)) { continue; // failed to propagate @@ -1118,8 +1170,9 @@ struct OnTheFlyTracker { posClusterCandidate[2] = gRandom->Gaus(posClusterCandidate[2], layer.getResolutionZ()); } - if (std::isnan(phi)) + if (std::isnan(phi)) { continue; // Catch when getXatLabR misses layer[i] + } // towards adding cluster: move to track alpha double alpha = cascadeTrack.getAlpha(); @@ -1128,8 +1181,10 @@ struct OnTheFlyTracker { -TMath::Sin(alpha) * posClusterCandidate[0] + TMath::Cos(alpha) * posClusterCandidate[1], posClusterCandidate[2]}; - if (!(cascadeTrack.propagateTo(xyz1[0], mMagneticField))) + if (!(cascadeTrack.propagateTo(xyz1[0], mMagneticField))) { continue; + } + const o2::track::TrackParametrization::dim2_t hitpoint = {static_cast(xyz1[1]), static_cast(xyz1[2])}; const o2::track::TrackParametrization::dim3_t hitpointcov = {layer.getResolutionRPhi() * layer.getResolutionRPhi(), 0.f, layer.getResolutionZ() * layer.getResolutionZ()}; if (layer.isInDeadPhiRegion(phi)) { @@ -1139,10 +1194,10 @@ struct OnTheFlyTracker { cascadeTrack.update(hitpoint, hitpointcov); thisCascade.foundClusters++; // add to findable } - } - if (thisCascade.foundClusters < cascadeDecaySettings.minStraTrackHits) { - continue; // We didn't find enough hits for strangeness tracking + if (thisCascade.foundClusters < cascadeDecaySettings.minStraTrackHits) { + continue; // We didn't find enough hits for strangeness tracking + } } // add cascade track @@ -1156,20 +1211,19 @@ struct OnTheFlyTracker { } // end cascade building if (isReco[0] && ((cascadeDecaySettings.doKinkReco == 1 && tryKinkReco) || cascadeDecaySettings.doKinkReco == 2)) { // mode 1 or 2 - o2::track::TrackParCov prefectCascadeTrack, trackedCasc; + o2::track::TrackParCov prefectCascadeTrack, trackedCascade; const o2::track::TrackParCov& trackedBach = xiDaughterTrackParCovsTracked[0]; o2::upgrade::convertMCParticleToO2Track(mcParticle, prefectCascadeTrack, pdgDB); // back track is already smeared prefectCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas - int nCascHits = fastTracker[icfg]->FastTrack(prefectCascadeTrack, trackedCasc, dNdEta); - reconstructedCascade = (fastTrackerSettings.minSiliconHitsForKinkReco < nCascHits) ? false : true; - + const int nCascHits = fastTracker[icfg]->FastTrack(prefectCascadeTrack, trackedCascade, dNdEta, xiDecayRadius2D); + reconstructedCascade = (fastTrackerSettings.minSiliconHitsForKinkReco < nCascHits) ? true : false; if (reconstructedCascade) { std::array pCasc; std::array pBach; std::array pV0; - trackedCasc.getPxPyPzGlo(pCasc); + trackedCascade.getPxPyPzGlo(pCasc); trackedBach.getPxPyPzGlo(pBach); for (size_t i = 0; i < pCasc.size(); ++i) { pV0[i] = pCasc[i] - pBach[i]; @@ -1192,7 +1246,7 @@ struct OnTheFlyTracker { int nCand = 0; bool kinkFitterOK = true; try { - nCand = fitter.process(trackedCasc, trackedBach); + nCand = fitter.process(trackedCascade, trackedBach); } catch (...) { kinkFitterOK = false; } @@ -1201,43 +1255,43 @@ struct OnTheFlyTracker { kinkFitterOK = false; } + fitter.propagateTracksToVertex(); + if (!fitter.isPropagateTracksToVertexDone()) { + kinkFitterOK = false; + } + + const u_int8_t fitterStatusCode = fitter.getFitStatus(); + histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); if (kinkFitterOK) { if (cascadeDecaySettings.doXiQA) { getHist(TH1, histPath + "hXiBuilding")->Fill(6.0f); } - } - fitter.propagateTracksToVertex(); // propagate e and K to D vertex - if (!fitter.isPropagateTracksToVertexDone()) { - kinkFitterOK = false; - } + o2::track::TrackParCov newCascadeTrack = fitter.getTrack(0); // (cascade) + std::array kinkVtx = {-999, -999, -999}; + kinkVtx = fitter.getPCACandidatePos(); + thisCascade.bachelorId = lastTrackIndex + tracksAlice3.size() - isReco.size(); + thisCascade.cascadeTrackId = lastTrackIndex + tracksAlice3.size(); // this should be ok + thisCascade.dcaV0dau = -1.f; // unknown + thisCascade.v0radius = -1.f; // unknown + thisCascade.dcacascdau = std::sqrt(fitter.getChi2AtPCACandidate()); + thisCascade.cascradius = std::hypot(kinkVtx[0], kinkVtx[1]); + thisCascade.cascradiusMC = xiDecayRadius2D; + thisCascade.mLambda = o2::constants::physics::MassLambda; + thisCascade.findableClusters = nCascHits; + thisCascade.foundClusters = nCascHits; + thisCascade.pt = newCascadeTrack.getPt(); + thisCascade.eta = newCascadeTrack.getEta(); + thisCascade.mXi = RecoDecay::m(std::array{std::array{pBach[0], pBach[1], pBach[2]}, std::array{pV0[0], pV0[1], pV0[2]}}, + std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda}); - o2::track::TrackParCov newCascadeTrack = fitter.getTrack(0); // (cascade) - std::array kinkVtx = {-999, -999, -999}; - kinkVtx = fitter.getPCACandidatePos(); - - thisCascade.bachelorId = lastTrackIndex + tracksAlice3.size() - isReco.size(); - thisCascade.cascadeTrackId = lastTrackIndex + tracksAlice3.size(); // this should be ok - thisCascade.dcaV0dau = -1.f; // unknown - thisCascade.v0radius = -1.f; // unknown - thisCascade.dcacascdau = std::sqrt(fitter.getChi2AtPCACandidate()); - thisCascade.cascradius = std::hypot(kinkVtx[0], kinkVtx[1]); - thisCascade.cascradiusMC = xiDecayRadius2D; - thisCascade.mLambda = o2::constants::physics::MassLambda; - thisCascade.findableClusters = nCascHits; - thisCascade.foundClusters = nCascHits; - thisCascade.pt = newCascadeTrack.getPt(); - thisCascade.eta = newCascadeTrack.getEta(); - thisCascade.mXi = RecoDecay::m(std::array{std::array{pBach[0], pBach[1], pBach[2]}, - std::array{pV0[0], pV0[1], pV0[2]}}, - std::array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda}); - - newCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas - tracksAlice3.push_back(TrackAlice3{newCascadeTrack, mcParticle.globalIndex(), time, timeResolutionUs, false, false, 1, thisCascade.foundClusters}); - - // add this cascade to vector (will fill cursor later with collision ID) - cascadesAlice3.push_back(thisCascade); - } + newCascadeTrack.setPID(pdgCodeToPID(PDG_t::kXiMinus)); // FIXME: not OK for omegas + tracksAlice3.push_back(TrackAlice3{newCascadeTrack, mcParticle.globalIndex(), time, timeResolutionUs, false, false, 1, thisCascade.foundClusters}); + + // add this cascade to vector (will fill cursor later with collision ID) + cascadesAlice3.push_back(thisCascade); + } // end fitter OK + } // end cascade found } // end cascade kink building // +-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+-~-+ @@ -1247,7 +1301,7 @@ struct OnTheFlyTracker { getHist(TH1, histPath + "hMassLambda")->Fill(thisCascade.mLambda); getHist(TH1, histPath + "hMassXi")->Fill(thisCascade.mXi); getHist(TH2, histPath + "h2dMassXi")->Fill(thisCascade.mXi, thisCascade.pt); - getHist(TH2, histPath + "h2dDeltaPtVsPt")->Fill(thisCascade.pt, mcParticle.pt() - thisCascade.pt); + getHist(TH2, histPath + "h2dDeltaPtVsPt")->Fill(thisCascade.pt, (mcParticle.pt() - thisCascade.pt) / thisCascade.pt); getHist(TH2, histPath + "h2dDeltaEtaVsPt")->Fill(thisCascade.pt, mcParticle.eta() - thisCascade.eta); getHist(TH2, histPath + "hFoundVsFindable")->Fill(thisCascade.findableClusters, thisCascade.foundClusters); } @@ -1382,6 +1436,14 @@ struct OnTheFlyTracker { if (nCand == 0) { dcaFitterOK_V0 = false; } + + fitter.propagateTracksToVertex(); + if (!fitter.isPropagateTracksToVertexDone()) { + dcaFitterOK_V0 = false; + } + + const u_int8_t fitterStatusCode = fitter.getFitStatus(); + histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); // V0 found successfully if (dcaFitterOK_V0) { if (v0DecaySettings.doV0QA) { diff --git a/ALICE3/TableProducer/alice3MulticharmFinder.cxx b/ALICE3/TableProducer/alice3MulticharmFinder.cxx index 79437c2a1a2..16e0e9571c7 100644 --- a/ALICE3/TableProducer/alice3MulticharmFinder.cxx +++ b/ALICE3/TableProducer/alice3MulticharmFinder.cxx @@ -87,7 +87,20 @@ struct Alice3MulticharmFinder { Configurable fillMCharmIdx{"fillMCharmIdx", true, "fill MCharmIdx[] tables (careful: memory)"}; Configurable fillMCharmCore{"fillMCharmCore", true, "fill MCharmCores[] tables (careful: memory)"}; Configurable fillMCharmExtra{"fillMCharmExtra", false, "fill MCharmExtra[] tables (careful: memory)"}; - } derivedTable; // allows for gap between peak and bg in case someone wants to + } derivedTable; + + struct : ConfigurableGroup { + std::string prefix = "cfgFitter"; + Configurable propagateToPCA{"propagateToPCA", false, "create tracks version propagated to PCA"}; + Configurable maxR{"maxR", 200., "reject PCA's above this radius"}; + Configurable minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"}; + Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; + Configurable maxDZIni{"maxDZIni", 1e9, "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; + Configurable maxDXYIni{"maxDXYIni", 4, "reject (if>0) PCA candidate if tracks DXY exceeds threshold"}; + Configurable maxVtxChi2{"maxVtxChi2", 1e9, "reject (if>0) vtx. chi2 above this value"}; + Configurable useAbsDCA{"useAbsDCA", true, "Minimise abs. distance rather than chi2"}; + Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; + } cfgFitter; Configurable cfgMagneticField{"cfgMagneticField", 20.0f, "Magnetic field (in kilogauss) if value not found from geo provider"}; Configurable doDCAplots{"doDCAplots", true, "do daughter prong DCA plots for D mesons"}; @@ -164,7 +177,6 @@ struct Alice3MulticharmFinder { // filter expressions for pions static constexpr uint32_t TrackSelectionPic = 1 << kInnerTOFPion | 1 << kOuterTOFPion | 1 << kRICHPion | 1 << kTruePiFromXiC; static constexpr uint32_t TrackSelectionPicc = 1 << kInnerTOFPion | 1 << kOuterTOFPion | 1 << kRICHPion | 1 << kTruePiFromXiCC; - float magneticField{}; // partitions Partition trueXi = aod::mcparticle::pdgCode == static_cast(PDG_t::kXiMinus); @@ -224,14 +236,19 @@ struct Alice3MulticharmFinder { } catch (...) { return false; } + + const u_int8_t fitterStatusCode = fitter.getFitStatus(); + histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); if (nCand == 0) { return false; } - const u_int8_t fitterStatusCode = fitter.getFitStatus(); - histos.fill(HIST("hFitterStatusCode"), fitterStatusCode); - //}-{}-{}-{}-{}-{}-{}-{}-{}-{} + fitter.propagateTracksToVertex(); + if (!fitter.isPropagateTracksToVertexDone()) { + return false; + } + //}-{}-{}-{}-{}-{}-{}-{}-{}-{} o2::track::TrackParCov t0new = fitter.getTrack(0); o2::track::TrackParCov t1new = fitter.getTrack(1); t0new.getPxPyPzGlo(thisXiccCandidate.prong0mom); @@ -298,12 +315,17 @@ struct Alice3MulticharmFinder { } catch (...) { return false; } + + const u_int8_t fitter3StatusCode = fitter3.getFitStatus(); + histos.fill(HIST("hFitter3StatusCode"), fitter3StatusCode); if (nCand == 0) { return false; } - const u_int8_t fitter3StatusCode = fitter3.getFitStatus(); - histos.fill(HIST("hFitter3StatusCode"), fitter3StatusCode); + fitter3.propagateTracksToVertex(); + if (!fitter3.isPropagateTracksToVertexDone()) { + return false; + } //}-{}-{}-{}-{}-{}-{}-{}-{}-{} t0 = fitter3.getTrack(0); @@ -405,38 +427,31 @@ struct Alice3MulticharmFinder { return returnValue; } - void init(o2::framework::InitContext& initContext) + void init(o2::framework::InitContext&) { - const bool foundMagneticField = common::core::getTaskOptionValue(initContext, "on-the-fly-detector-geometry-provider", "magneticField", magneticField, false); - if (!foundMagneticField) { - LOG(info) << "Could not retrieve magnetic field from geometry provider."; - LOG(info) << "Using value from configurable cfgMagneticField: " << cfgMagneticField; - magneticField = cfgMagneticField; - } else { - LOG(info) << "Using magnetic field form geometry provider with value: " << magneticField; - } - // initialize O2 2-prong fitter (only once) - fitter.setPropagateToPCA(true); - fitter.setMaxR(200.); - fitter.setMinParamChange(1e-3); - fitter.setMinRelChi2Change(0.9); - fitter.setMaxDZIni(1e9); - fitter.setMaxChi2(1e9); - fitter.setUseAbsDCA(true); - fitter.setWeightedFinalPCA(false); - fitter.setBz(magneticField); + fitter.setPropagateToPCA(cfgFitter.propagateToPCA); + fitter.setMaxR(cfgFitter.maxR); + fitter.setMinParamChange(cfgFitter.minParamChange); + fitter.setMinRelChi2Change(cfgFitter.minRelChi2Change); + fitter.setMaxDZIni(cfgFitter.maxDZIni); + fitter.setMaxDXYIni(cfgFitter.maxDXYIni); + fitter.setMaxChi2(cfgFitter.maxVtxChi2); + fitter.setUseAbsDCA(cfgFitter.useAbsDCA); + fitter.setWeightedFinalPCA(cfgFitter.useWeightedFinalPCA); + fitter.setBz(cfgMagneticField); fitter.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrNONE); - fitter3.setPropagateToPCA(true); - fitter3.setMaxR(200.); - fitter3.setMinParamChange(1e-3); - fitter3.setMinRelChi2Change(0.9); - fitter3.setMaxDZIni(1e9); - fitter3.setMaxChi2(1e9); - fitter3.setUseAbsDCA(true); - fitter3.setWeightedFinalPCA(false); - fitter3.setBz(magneticField); + fitter3.setPropagateToPCA(cfgFitter.propagateToPCA); + fitter3.setMaxR(cfgFitter.maxR); + fitter3.setMinParamChange(cfgFitter.minParamChange); + fitter3.setMinRelChi2Change(cfgFitter.minRelChi2Change); + fitter3.setMaxDZIni(cfgFitter.maxDZIni); + fitter3.setMaxDZIni(cfgFitter.maxDXYIni); + fitter3.setMaxChi2(cfgFitter.maxVtxChi2); + fitter3.setUseAbsDCA(cfgFitter.useAbsDCA); + fitter3.setWeightedFinalPCA(cfgFitter.useWeightedFinalPCA); + fitter3.setBz(cfgMagneticField); fitter3.setMatCorrType(o2::base::Propagator::MatCorrType::USEMatCorrNONE); auto hFitterStatusCode = histos.add("hFitterStatusCode", "hFitterStatusCode", kTH1D, {{15, -0.5, 14.5}}); @@ -713,7 +728,7 @@ struct Alice3MulticharmFinder { o2::vertexing::PVertex primaryVertex; primaryVertex.setXYZ(collision.posX(), collision.posY(), collision.posZ()); - if (xicTrackCopy.propagateToDCA(primaryVertex, magneticField, &dcaInfo)) { + if (xicTrackCopy.propagateToDCA(primaryVertex, cfgMagneticField, &dcaInfo)) { xicdcaXY = dcaInfo.getY(); xicdcaZ = dcaInfo.getZ(); } @@ -808,7 +823,7 @@ struct Alice3MulticharmFinder { GET_HIST(TH1, histPath + "hMultiCharmBuilding")->Fill(6.0f); GET_HIST(TH2, histPath + "hXicRadiusVsXiccRadius")->Fill(xicDecayRadius2D * ToMicrons, xiccDecayRadius2D * ToMicrons); float xiccdcaXY = 1e+10, xiccdcaZ = 1e+10; - if (xiccTrack.propagateToDCA(primaryVertex, magneticField, &dcaInfo)) { + if (xiccTrack.propagateToDCA(primaryVertex, cfgMagneticField, &dcaInfo)) { xiccdcaXY = dcaInfo.getY(); xiccdcaZ = dcaInfo.getZ(); }