diff --git a/Tools/KFparticle/CMakeLists.txt b/Tools/KFparticle/CMakeLists.txt index eb0e6571b46..096d53fd09e 100644 --- a/Tools/KFparticle/CMakeLists.txt +++ b/Tools/KFparticle/CMakeLists.txt @@ -1,4 +1,3 @@ - # Copyright 2019-2020 CERN and copyright holders of ALICE O2. # See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. # All rights not expressly granted are reserved. @@ -11,11 +10,16 @@ # or submit itself to any jurisdiction. o2physics_add_dpl_workflow(qa-kfparticle - SOURCES qaKFParticle.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle - COMPONENT_NAME Analysis) + SOURCES qaKFParticle.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(qa-kfparticle-lc + SOURCES qaKFParticleLc.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle + COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(qa-kfeventtrack - SOURCES qaKFEventTrack.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle - COMPONENT_NAME Analysis) + SOURCES qaKFEventTrack.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle + COMPONENT_NAME Analysis) diff --git a/Tools/KFparticle/KFUtilities.h b/Tools/KFparticle/KFUtilities.h index 2b64c70bdb5..ce97274eab9 100644 --- a/Tools/KFparticle/KFUtilities.h +++ b/Tools/KFparticle/KFUtilities.h @@ -43,8 +43,8 @@ KFPVertex createKFPVertexFromCollision(const T& collision) kfpVertex.SetXYZ(collision.posX(), collision.posY(), collision.posZ()); kfpVertex.SetCovarianceMatrix(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()); kfpVertex.SetChi2(collision.chi2()); - kfpVertex.SetNDF(2 * collision.multNTracksPV() - 3); - kfpVertex.SetNContributors(collision.multNTracksPV()); + kfpVertex.SetNDF(2 * collision.numContrib() - 3); + kfpVertex.SetNContributors(collision.numContrib()); return kfpVertex; } diff --git a/Tools/KFparticle/qaKFEventTrack.cxx b/Tools/KFparticle/qaKFEventTrack.cxx index dba08e165b7..e6ba40adaeb 100644 --- a/Tools/KFparticle/qaKFEventTrack.cxx +++ b/Tools/KFparticle/qaKFEventTrack.cxx @@ -88,6 +88,7 @@ struct qaKFEventTrack { Configurable d_etaRange{"d_etaRange", 0.8, "eta Range for tracks"}; /// Option to write variables in a tree Configurable d_DwnSmplFact{"d_DwnSmplFact", 1., "Downsampling factor for tree"}; + Configurable writeHistograms{"writeHistograms", true, "write histograms"}; Configurable writeTree{"writeTree", false, "write daughter variables in a tree"}; // Define which track selection should be used: @@ -158,54 +159,56 @@ struct qaKFEventTrack { } runNumber = 0; - const AxisSpec axisVertexPosX{500, -0.05, 0.05, "X [cm]"}; - const AxisSpec axisVertexPosY{500, -0.05, 0.05, "Y [cm]"}; + const AxisSpec axisVertexPosX{100, -0.05, 0.05, "X [cm]"}; + const AxisSpec axisVertexPosY{100, -0.05, 0.05, "Y [cm]"}; const AxisSpec axisVertexPosZ{100, -20., 20., "Z [cm]"}; const AxisSpec axisVertexNumContrib{160, 0, 160, "Number Of contributors to the PV"}; const AxisSpec axisVertexCov{100, -0.00005, 0.00005}; - const AxisSpec axisParX{300, -0.1, 0.1, "#it{x} [cm]"}; - const AxisSpec axisParY{200, -0.1, 0.1, "#it{y} [cm]"}; - const AxisSpec axisParZ{200, -20., 20., "#it{z} [cm]"}; + const AxisSpec axisParX{100, -0.1, 0.1, "#it{x} [cm]"}; + const AxisSpec axisParY{100, -0.1, 0.1, "#it{y} [cm]"}; + const AxisSpec axisParZ{100, -20., 20., "#it{z} [cm]"}; const AxisSpec axisParPX{binsPt, "#it{p}_{x} [GeV/c]"}; const AxisSpec axisParPY{binsPt, "#it{p}_{y} [GeV/c]"}; const AxisSpec axisParPZ{binsPt, "#it{p}_{z} [GeV/c]"}; - /// collisions - histos.add("Events/covXX", ";Cov_{xx} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("Events/covXY", ";Cov_{xy} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("Events/covYY", ";Cov_{yy} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("Events/covXZ", ";Cov_{xz} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("Events/covYZ", ";Cov_{yz} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("Events/covZZ", ";Cov_{zz} [cm^{2}]", kTH1D, {axisVertexCov}); - - histos.add("EventsKF/posX", "", kTH1D, {axisVertexPosX}); - histos.add("EventsKF/posY", "", kTH1D, {axisVertexPosY}); - histos.add("EventsKF/posZ", "", kTH1D, {axisVertexPosZ}); - histos.add("EventsKF/posXY", "", kTH2D, {axisVertexPosX, axisVertexPosY}); - histos.add("EventsKF/nContrib", "", kTH1D, {axisVertexNumContrib}); - histos.add("EventsKF/vertexChi2", ";#chi^{2}", kTH1D, {{100, 0, 100}}); - histos.add("EventsKF/covXX", ";Cov_{xx} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("EventsKF/covXY", ";Cov_{xy} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("EventsKF/covYY", ";Cov_{yy} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("EventsKF/covXZ", ";Cov_{xz} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("EventsKF/covYZ", ";Cov_{yz} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("EventsKF/covZZ", ";Cov_{zz} [cm^{2}]", kTH1D, {axisVertexCov}); - - /// tracks - histos.add("Tracks/x", "track #it{x} position at dca in local coordinate system", kTH1D, {axisParX}); - histos.add("Tracks/y", "track #it{y} position at dca in local coordinate system", kTH1D, {axisParY}); - histos.add("Tracks/z", "track #it{z} position at dca in local coordinate system", kTH1D, {axisParZ}); - histos.add("Tracks/px", "track #it{p_{x}} momentum at dca in local coordinate system", kTH1D, {axisParPX}); - histos.add("Tracks/py", "track #it{p_{y}} momentum at dca in local coordinate system", kTH1D, {axisParPY}); - histos.add("Tracks/pz", "track #it{p_{z}} momentum at dca in local coordinate system", kTH1D, {axisParPZ}); - histos.add("Tracks/chi2perNDF", "Chi2/NDF of the track;#it{chi2/ndf};", kTH1D, {{200, 0., 25.}}); - histos.add("Tracks/dcaXYToPV", "distance of closest approach in #it{xy} plane;#it{dcaXY} [cm];", kTH1D, {{200, -0.15, 0.15}}); - histos.add("Tracks/dcaToPV", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, -0.15, 0.15}}); - histos.add("Tracks/dcaToPVLargeRange", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); - histos.add("Tracks/dcaToPVLargeRangeMCBeforeReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); - histos.add("Tracks/dcaToPVLargeRangeMCAfterReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); - histos.add("Tracks/deviationPiToPV", "deviation of Pi to PV", kTH1D, {{200, 0., 10.}}); + if (writeHistograms) { + /// collisions + histos.add("Events/covXX", ";Cov_{xx} [cm^{2}]", kTH1D, {axisVertexCov}); + histos.add("Events/covXY", ";Cov_{xy} [cm^{2}]", kTH1D, {axisVertexCov}); + histos.add("Events/covYY", ";Cov_{yy} [cm^{2}]", kTH1D, {axisVertexCov}); + histos.add("Events/covXZ", ";Cov_{xz} [cm^{2}]", kTH1D, {axisVertexCov}); + histos.add("Events/covYZ", ";Cov_{yz} [cm^{2}]", kTH1D, {axisVertexCov}); + histos.add("Events/covZZ", ";Cov_{zz} [cm^{2}]", kTH1D, {axisVertexCov}); + + histos.add("EventsKF/posX", "", kTH1D, {axisVertexPosX}); + histos.add("EventsKF/posY", "", kTH1D, {axisVertexPosY}); + histos.add("EventsKF/posZ", "", kTH1D, {axisVertexPosZ}); + histos.add("EventsKF/posXY", "", kTH2D, {axisVertexPosX, axisVertexPosY}); + histos.add("EventsKF/nContrib", "", kTH1D, {axisVertexNumContrib}); + histos.add("EventsKF/vertexChi2", ";#chi^{2}", kTH1D, {{100, 0, 100}}); + histos.add("EventsKF/covXX", ";Cov_{xx} [cm^{2}]", kTH1D, {axisVertexCov}); + histos.add("EventsKF/covXY", ";Cov_{xy} [cm^{2}]", kTH1D, {axisVertexCov}); + histos.add("EventsKF/covYY", ";Cov_{yy} [cm^{2}]", kTH1D, {axisVertexCov}); + histos.add("EventsKF/covXZ", ";Cov_{xz} [cm^{2}]", kTH1D, {axisVertexCov}); + histos.add("EventsKF/covYZ", ";Cov_{yz} [cm^{2}]", kTH1D, {axisVertexCov}); + histos.add("EventsKF/covZZ", ";Cov_{zz} [cm^{2}]", kTH1D, {axisVertexCov}); + + /// tracks + histos.add("Tracks/x", "track #it{x} position at dca in local coordinate system", kTH1D, {axisParX}); + histos.add("Tracks/y", "track #it{y} position at dca in local coordinate system", kTH1D, {axisParY}); + histos.add("Tracks/z", "track #it{z} position at dca in local coordinate system", kTH1D, {axisParZ}); + histos.add("Tracks/px", "track #it{p_{x}} momentum at dca in local coordinate system", kTH1D, {axisParPX}); + histos.add("Tracks/py", "track #it{p_{y}} momentum at dca in local coordinate system", kTH1D, {axisParPY}); + histos.add("Tracks/pz", "track #it{p_{z}} momentum at dca in local coordinate system", kTH1D, {axisParPZ}); + histos.add("Tracks/chi2perNDF", "Chi2/NDF of the track;#it{chi2/ndf};", kTH1D, {{100, 0., 25.}}); + histos.add("Tracks/dcaXYToPV", "distance of closest approach in #it{xy} plane;#it{dcaXY} [cm];", kTH1D, {{100, -0.15, 0.15}}); + histos.add("Tracks/dcaToPV", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{100, -0.15, 0.15}}); + histos.add("Tracks/dcaToPVLargeRange", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); + histos.add("Tracks/dcaToPVLargeRangeMCBeforeReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{100, 0., 20.}}); + histos.add("Tracks/dcaToPVLargeRangeMCAfterReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{100, 0., 20.}}); + histos.add("Tracks/deviationPiToPV", "deviation of Pi to PV", kTH1D, {{100, 0., 10.}}); + } auto hSelectionMC = histos.add("DZeroCandTopo/SelectionsMC", "Selections MC", kTH1D, {{5, 0.5, 5.5}}); hSelectionMC->GetXaxis()->SetBinLabel(hSelectionMC->FindBin(1), "All Tracks"); @@ -243,18 +246,20 @@ struct qaKFEventTrack { template void fillHistograms(const T1& kfpTrackPi, const T2& KFPion, const T2& KFPV) { - /// fill daughter track parameters - histos.fill(HIST("Tracks/x"), kfpTrackPi.GetX()); - histos.fill(HIST("Tracks/y"), kfpTrackPi.GetY()); - histos.fill(HIST("Tracks/z"), kfpTrackPi.GetZ()); - histos.fill(HIST("Tracks/px"), kfpTrackPi.GetPx()); - histos.fill(HIST("Tracks/py"), kfpTrackPi.GetPy()); - histos.fill(HIST("Tracks/pz"), kfpTrackPi.GetPz()); - histos.fill(HIST("Tracks/chi2perNDF"), kfpTrackPi.GetChi2perNDF()); - histos.fill(HIST("Tracks/dcaXYToPV"), KFPion.GetDistanceFromVertexXY(KFPV)); - histos.fill(HIST("Tracks/dcaToPV"), KFPion.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("Tracks/dcaToPVLargeRange"), KFPion.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("Tracks/deviationPiToPV"), KFPion.GetDeviationFromVertex(KFPV)); + if (writeHistograms) { + /// fill daughter track parameters + histos.fill(HIST("Tracks/x"), kfpTrackPi.GetX()); + histos.fill(HIST("Tracks/y"), kfpTrackPi.GetY()); + histos.fill(HIST("Tracks/z"), kfpTrackPi.GetZ()); + histos.fill(HIST("Tracks/px"), kfpTrackPi.GetPx()); + histos.fill(HIST("Tracks/py"), kfpTrackPi.GetPy()); + histos.fill(HIST("Tracks/pz"), kfpTrackPi.GetPz()); + histos.fill(HIST("Tracks/chi2perNDF"), kfpTrackPi.GetChi2perNDF()); + histos.fill(HIST("Tracks/dcaXYToPV"), KFPion.GetDistanceFromVertexXY(KFPV)); + histos.fill(HIST("Tracks/dcaToPV"), KFPion.GetDistanceFromVertex(KFPV)); + histos.fill(HIST("Tracks/dcaToPVLargeRange"), KFPion.GetDistanceFromVertex(KFPV)); + histos.fill(HIST("Tracks/deviationPiToPV"), KFPion.GetDeviationFromVertex(KFPV)); + } } template @@ -300,13 +305,14 @@ struct qaKFEventTrack { KFParticle::SetField(magneticField); #endif } - - histos.fill(HIST("Events/covXX"), collision.covXX()); - histos.fill(HIST("Events/covXY"), collision.covXY()); - histos.fill(HIST("Events/covXZ"), collision.covXZ()); - histos.fill(HIST("Events/covYY"), collision.covYY()); - histos.fill(HIST("Events/covYZ"), collision.covYZ()); - histos.fill(HIST("Events/covZZ"), collision.covZZ()); + if (writeHistograms) { + histos.fill(HIST("Events/covXX"), collision.covXX()); + histos.fill(HIST("Events/covXY"), collision.covXY()); + histos.fill(HIST("Events/covXZ"), collision.covXZ()); + histos.fill(HIST("Events/covYY"), collision.covYY()); + histos.fill(HIST("Events/covYZ"), collision.covYZ()); + histos.fill(HIST("Events/covZZ"), collision.covZZ()); + } /// Apply event selection if (!isSelectedCollision(collision)) { return; @@ -314,20 +320,21 @@ struct qaKFEventTrack { /// set KF primary vertex KFPVertex kfpVertex = createKFPVertexFromCollision(collision); KFParticle KFPV(kfpVertex); - - /// fill collision parameters - histos.fill(HIST("EventsKF/posX"), kfpVertex.GetX()); - histos.fill(HIST("EventsKF/posY"), kfpVertex.GetY()); - histos.fill(HIST("EventsKF/posZ"), kfpVertex.GetZ()); - histos.fill(HIST("EventsKF/posXY"), kfpVertex.GetX(), kfpVertex.GetY()); - histos.fill(HIST("EventsKF/nContrib"), kfpVertex.GetNContributors()); - histos.fill(HIST("EventsKF/vertexChi2"), kfpVertex.GetChi2()); - histos.fill(HIST("EventsKF/covXX"), kfpVertex.GetCovariance(0)); - histos.fill(HIST("EventsKF/covXY"), kfpVertex.GetCovariance(1)); - histos.fill(HIST("EventsKF/covYY"), kfpVertex.GetCovariance(2)); - histos.fill(HIST("EventsKF/covXZ"), kfpVertex.GetCovariance(3)); - histos.fill(HIST("EventsKF/covYZ"), kfpVertex.GetCovariance(4)); - histos.fill(HIST("EventsKF/covZZ"), kfpVertex.GetCovariance(5)); + if (writeHistograms) { + /// fill collision parameters + histos.fill(HIST("EventsKF/posX"), kfpVertex.GetX()); + histos.fill(HIST("EventsKF/posY"), kfpVertex.GetY()); + histos.fill(HIST("EventsKF/posZ"), kfpVertex.GetZ()); + histos.fill(HIST("EventsKF/posXY"), kfpVertex.GetX(), kfpVertex.GetY()); + histos.fill(HIST("EventsKF/nContrib"), kfpVertex.GetNContributors()); + histos.fill(HIST("EventsKF/vertexChi2"), kfpVertex.GetChi2()); + histos.fill(HIST("EventsKF/covXX"), kfpVertex.GetCovariance(0)); + histos.fill(HIST("EventsKF/covXY"), kfpVertex.GetCovariance(1)); + histos.fill(HIST("EventsKF/covYY"), kfpVertex.GetCovariance(2)); + histos.fill(HIST("EventsKF/covXZ"), kfpVertex.GetCovariance(3)); + histos.fill(HIST("EventsKF/covYZ"), kfpVertex.GetCovariance(4)); + histos.fill(HIST("EventsKF/covZZ"), kfpVertex.GetCovariance(5)); + } int ntracks = 0; @@ -374,7 +381,7 @@ struct qaKFEventTrack { using CollisionTableMC = soa::Join; using CollisionTableDataMult = soa::Join; using TrackTableMC = soa::Join; - Preslice perMcCollision = aod::mccollisionlabel::mcCollisionId; + // Preslice perMcCollision = aod::mccollisionlabel::mcCollisionId; void processMC(CollisionTableMC::iterator const& collision, CollisionTableMC const& collisions, soa::Filtered const& tracks, aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions, aod::BCsWithTimestamps const&) { auto bc = collision.bc_as(); @@ -391,12 +398,15 @@ struct qaKFEventTrack { if (!collision.has_mcCollision()) { return; } - histos.fill(HIST("Events/covXX"), collision.covXX()); - histos.fill(HIST("Events/covXY"), collision.covXY()); - histos.fill(HIST("Events/covXZ"), collision.covXZ()); - histos.fill(HIST("Events/covYY"), collision.covYY()); - histos.fill(HIST("Events/covYZ"), collision.covYZ()); - histos.fill(HIST("Events/covZZ"), collision.covZZ()); + if (writeHistograms) { + histos.fill(HIST("Events/covXX"), collision.covXX()); + histos.fill(HIST("Events/covXY"), collision.covXY()); + histos.fill(HIST("Events/covXZ"), collision.covXZ()); + histos.fill(HIST("Events/covYY"), collision.covYY()); + histos.fill(HIST("Events/covYZ"), collision.covYZ()); + histos.fill(HIST("Events/covZZ"), collision.covZZ()); + } + /// Apply event selection if (!isSelectedCollision(collision)) { return; @@ -407,20 +417,21 @@ struct qaKFEventTrack { KFPVertex kfpVertexDefault = createKFPVertexFromCollision(collision); KFParticle KFPVDefault(kfpVertexDefault); - - /// fill collision parameters - histos.fill(HIST("EventsKF/posX"), kfpVertex.GetX()); - histos.fill(HIST("EventsKF/posY"), kfpVertex.GetY()); - histos.fill(HIST("EventsKF/posZ"), kfpVertex.GetZ()); - histos.fill(HIST("EventsKF/posXY"), kfpVertex.GetX(), kfpVertex.GetY()); - histos.fill(HIST("EventsKF/nContrib"), kfpVertex.GetNContributors()); - histos.fill(HIST("EventsKF/vertexChi2"), kfpVertex.GetChi2()); - histos.fill(HIST("EventsKF/covXX"), kfpVertex.GetCovariance(0)); - histos.fill(HIST("EventsKF/covXY"), kfpVertex.GetCovariance(1)); - histos.fill(HIST("EventsKF/covYY"), kfpVertex.GetCovariance(2)); - histos.fill(HIST("EventsKF/covXZ"), kfpVertex.GetCovariance(3)); - histos.fill(HIST("EventsKF/covYZ"), kfpVertex.GetCovariance(4)); - histos.fill(HIST("EventsKF/covZZ"), kfpVertex.GetCovariance(5)); + if (writeHistograms) { + /// fill collision parameters + histos.fill(HIST("EventsKF/posX"), kfpVertex.GetX()); + histos.fill(HIST("EventsKF/posY"), kfpVertex.GetY()); + histos.fill(HIST("EventsKF/posZ"), kfpVertex.GetZ()); + histos.fill(HIST("EventsKF/posXY"), kfpVertex.GetX(), kfpVertex.GetY()); + histos.fill(HIST("EventsKF/nContrib"), kfpVertex.GetNContributors()); + histos.fill(HIST("EventsKF/vertexChi2"), kfpVertex.GetChi2()); + histos.fill(HIST("EventsKF/covXX"), kfpVertex.GetCovariance(0)); + histos.fill(HIST("EventsKF/covXY"), kfpVertex.GetCovariance(1)); + histos.fill(HIST("EventsKF/covYY"), kfpVertex.GetCovariance(2)); + histos.fill(HIST("EventsKF/covXZ"), kfpVertex.GetCovariance(3)); + histos.fill(HIST("EventsKF/covYZ"), kfpVertex.GetCovariance(4)); + histos.fill(HIST("EventsKF/covZZ"), kfpVertex.GetCovariance(5)); + } int ntracks = 0; @@ -440,51 +451,51 @@ struct qaKFEventTrack { continue; } - /// Check whether the track was assigned to the true MC PV - auto particle = track.mcParticle(); - auto collMC = particle.mcCollision(); - auto mcCollID_recoColl = track.collision_as().mcCollisionId(); - auto mcCollID_particle = particle.mcCollisionId(); - bool indexMatchOK = (mcCollID_recoColl == mcCollID_particle); - if (!indexMatchOK) { - histos.fill(HIST("DZeroCandTopo/SelectionsMC"), 4.f); - const auto matchedCollisions = collisions.sliceBy(perMcCollision, collMC.globalIndex()); - int i = 0; - std::array dcaZ{100, 100, 100, 100, 100}; - float min = 100; - for (auto matchedCollision : matchedCollisions) { - dcaZ[i] = abs(matchedCollision.posZ() - collMC.posZ()); - if (i == 0) { - min = dcaZ[i]; - } - if (i > 0) { - if (dcaZ[i] < dcaZ[i - 1]) { - min = dcaZ[i]; - } - } - - i = i + 1; - } - if (min > 10.) { - histos.fill(HIST("DZeroCandTopo/SelectionsMC"), 5.f); - } - int j = 0; - for (auto matchedCollision : matchedCollisions) { - if (i == 1) { - kfpVertex = createKFPVertexFromCollision(matchedCollision); - KFParticle KFPVNew(kfpVertex); - KFPV = KFPVNew; - } - if (i > 1) { - if (abs(matchedCollision.posZ() - collMC.posZ()) == min) { - kfpVertex = createKFPVertexFromCollision(matchedCollision); - KFParticle KFPVNew(kfpVertex); - KFPV = KFPVNew; - } - } - j = j + 1; - } - } + // /// Check whether the track was assigned to the true MC PV + // auto particle = track.mcParticle(); + // auto collMC = particle.mcCollision(); + // auto mcCollID_recoColl = track.collision_as().mcCollisionId(); + // auto mcCollID_particle = particle.mcCollisionId(); + // bool indexMatchOK = (mcCollID_recoColl == mcCollID_particle); + // if (!indexMatchOK) { + // histos.fill(HIST("DZeroCandTopo/SelectionsMC"), 4.f); + // const auto matchedCollisions = collisions.sliceBy(perMcCollision, collMC.globalIndex()); + // int i = 0; + // std::array dcaZ{100, 100, 100, 100, 100}; + // float min = 100; + // for (auto matchedCollision : matchedCollisions) { + // dcaZ[i] = abs(matchedCollision.posZ() - collMC.posZ()); + // if (i == 0) { + // min = dcaZ[i]; + // } + // if (i > 0) { + // if (dcaZ[i] < dcaZ[i - 1]) { + // min = dcaZ[i]; + // } + // } + + // i = i + 1; + // } + // if (min > 10.) { + // histos.fill(HIST("DZeroCandTopo/SelectionsMC"), 5.f); + // } + // int j = 0; + // for (auto matchedCollision : matchedCollisions) { + // if (i == 1) { + // kfpVertex = createKFPVertexFromCollision(matchedCollision); + // KFParticle KFPVNew(kfpVertex); + // KFPV = KFPVNew; + // } + // if (i > 1) { + // if (abs(matchedCollision.posZ() - collMC.posZ()) == min) { + // kfpVertex = createKFPVertexFromCollision(matchedCollision); + // KFParticle KFPVNew(kfpVertex); + // KFPV = KFPVNew; + // } + // } + // j = j + 1; + // } + // } /// Remove fake tracks if (!track.has_mcParticle()) { @@ -519,9 +530,10 @@ struct qaKFEventTrack { fillHistograms(kfpTrack, KFParticleTrack, KFPV); writeVarTree(kfpTrack, KFParticleTrack, KFPV, track, source, pVContrib, runNumber, collision); - - histos.fill(HIST("Tracks/dcaToPVLargeRangeMCBeforeReassignment"), KFParticleTrack.GetDistanceFromVertex(KFPVDefault)); - histos.fill(HIST("Tracks/dcaToPVLargeRangeMCAfterReassignment"), KFParticleTrack.GetDistanceFromVertex(KFPV)); + if (writeHistograms) { + histos.fill(HIST("Tracks/dcaToPVLargeRangeMCBeforeReassignment"), KFParticleTrack.GetDistanceFromVertex(KFPVDefault)); + histos.fill(HIST("Tracks/dcaToPVLargeRangeMCAfterReassignment"), KFParticleTrack.GetDistanceFromVertex(KFPV)); + } } } PROCESS_SWITCH(qaKFEventTrack, processMC, "process mc", false); @@ -541,6 +553,7 @@ struct qaKFEvent { o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; int runNumber; double magneticField = 0.; + int tfID = 0; /// option to select good events Configurable eventSelection{"eventSelection", true, "select good events"}; // currently only sel8 is defined for run3 Configurable writeTree{"writeTree", true, "write daughter variables in a tree"}; @@ -603,7 +616,7 @@ struct qaKFEvent { } template - void writeVarTreeColl(const T4& collision, double timeColl, double timestamp, double timeDiff) + void writeVarTreeColl(const T4& collision, uint64_t timeColl, double timestamp, double timeDiff, int tfID) { if (writeTree) { /// Filling the tree @@ -619,18 +632,25 @@ struct qaKFEvent { runNumber, timeColl, timestamp, - timeDiff); + timeDiff, + collision.bcId(), + tfID); } } /// Process function for data void processCollisions(CollisionTableData const& collisions, aod::BCsWithTimestamps const&) { - double timeColl = 0; - double timestamp = 0; - double timeDiff = 0; + uint64_t timeColl = 0; + tfID = tfID + 1; + LOGP(info, "processing TF {}", tfID); + int bc0 = 0; for (auto& collisionIndex : collisions) { auto bc = collisionIndex.bc_as(); + if (bc0 == 0) { + bc0 = bc.globalBC(); + } + // LOGP(info, "{} vs {}", collisionIndex.bcId(), bc.globalIndex()); if (runNumber != bc.runNumber()) { initMagneticFieldCCDB(bc, runNumber, ccdb, isRun3 ? ccdbPathGrpMag : ccdbPathGrp, lut, isRun3); magneticField = o2::base::Propagator::Instance()->getNominalBz(); @@ -639,10 +659,8 @@ struct qaKFEvent { if (!isSelectedCollision(collisionIndex)) { continue; } - timestamp = bc.timestamp() * 1.e6; - timeDiff = collisionIndex.collisionTime(); - timeColl = timestamp + timeDiff; - writeVarTreeColl(collisionIndex, timeColl, bc.timestamp(), timeDiff); + timeColl = uint64_t((bc.globalBC() - bc0) * o2::constants::lhc::LHCBunchSpacingNS) + collisionIndex.collisionTime(); + writeVarTreeColl(collisionIndex, timeColl, bc.timestamp(), collisionIndex.collisionTime(), tfID); } } PROCESS_SWITCH(qaKFEvent, processCollisions, "process collision", true); diff --git a/Tools/KFparticle/qaKFEventTrack.h b/Tools/KFparticle/qaKFEventTrack.h index 77d54b3c4f1..bec8a5e1dfb 100644 --- a/Tools/KFparticle/qaKFEventTrack.h +++ b/Tools/KFparticle/qaKFEventTrack.h @@ -44,9 +44,11 @@ DECLARE_SOA_COLUMN(ETA, Eta, float); DECLARE_SOA_COLUMN(PHI, Phi, float); DECLARE_SOA_COLUMN(TPCSIGNAL, Tpcsignal, float); DECLARE_SOA_COLUMN(RUNNUMBER, Runnumber, float); -DECLARE_SOA_COLUMN(TIMECOLL, TimeColl, double); +DECLARE_SOA_COLUMN(TIMECOLL, TimeColl, uint64_t); DECLARE_SOA_COLUMN(TIMESTAMP, TimeStamp, double); DECLARE_SOA_COLUMN(TIMEDIFF, TimeDiff, double); +DECLARE_SOA_COLUMN(BCID, BCid, int); +DECLARE_SOA_COLUMN(TFID, Tfid, int); DECLARE_SOA_COLUMN(XPV, Xpv, float); DECLARE_SOA_COLUMN(YPV, Ypv, float); DECLARE_SOA_COLUMN(ZPV, Zpv, float); @@ -94,7 +96,9 @@ DECLARE_SOA_TABLE(TreeCollisions, "AOD", "TREECOLLISIONS", kfeventtrack::RUNNUMBER, kfeventtrack::TIMECOLL, kfeventtrack::TIMESTAMP, - kfeventtrack::TIMEDIFF); + kfeventtrack::TIMEDIFF, + kfeventtrack::BCID, + kfeventtrack::TFID); } // namespace o2::aod #endif // TOOLS_KFPARTICLE_QAKFEVENTTRACK_H_ diff --git a/Tools/KFparticle/qaKFParticle.cxx b/Tools/KFparticle/qaKFParticle.cxx index 8cd437501b5..b52ae5f72c9 100644 --- a/Tools/KFparticle/qaKFParticle.cxx +++ b/Tools/KFparticle/qaKFParticle.cxx @@ -38,12 +38,10 @@ #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/TrackSelectionTables.h" #include "Common/DataModel/PIDResponse.h" -#include "Common/DataModel/Multiplicity.h" #include "Common/Core/trackUtilities.h" #include "Common/Core/TrackSelection.h" #include "Common/Core/TrackSelectionDefaults.h" #include "Common/Core/RecoDecay.h" -#include "Common/Core/TrackSelectorPID.h" #include "Tools/KFparticle/KFUtilities.h" /// includes KFParticle @@ -90,30 +88,23 @@ struct qaKFParticle { /// Particle Identification // TPC PID - Configurable ptPidTpcMinPi{"ptPidTpcMinPi", 0.15, "Lower bound of track pT for TPC PID"}; - Configurable ptPidTpcMaxPi{"ptPidTpcMaxPi", 5., "Upper bound of track pT for TPC PID"}; Configurable nSigmaTpcMaxPi{"nSigmaTpcMaxPi", 3., "Nsigma cut on TPC only"}; - Configurable nSigmaTpcCombinedMaxPi{"nSigmaTpcCombinedMaxPi", 5., "Nsigma cut on TPC combined with TOF"}; - Configurable ptPidTpcMinKa{"ptPidTpcMinKa", 0.15, "Lower bound of track pT for TPC PID"}; - Configurable ptPidTpcMaxKa{"ptPidTpcMaxKa", 5., "Upper bound of track pT for TPC PID"}; Configurable nSigmaTpcMaxKa{"nSigmaTpcMaxKa", 3., "Nsigma cut on TPC only"}; - Configurable nSigmaTpcCombinedMaxKa{"nSigmaTpcCombinedMaxKa", 5., "Nsigma cut on TPC combined with TOF"}; // TOF PID Configurable ptPidTofMinPi{"ptPidTofMinPi", 0.15, "Lower bound of track pT for TOF PID"}; - Configurable ptPidTofMaxPi{"ptPidTofMaxPi", 5., "Upper bound of track pT for TOF PID"}; Configurable nSigmaTofMaxPi{"nSigmaTofMaxPi", 3., "Nsigma cut on TOF only"}; - Configurable nSigmaTofCombinedMaxPi{"nSigmaTofCombinedMaxPi", 5., "Nsigma cut on TOF combined with TPC"}; Configurable ptPidTofMinKa{"ptPidTofMinKa", 0.15, "Lower bound of track pT for TOF PID"}; - Configurable ptPidTofMaxKa{"ptPidTofMaxKa", 5., "Upper bound of track pT for TOF PID"}; Configurable nSigmaTofMaxKa{"nSigmaTofMaxKa", 3., "Nsigma cut on TOF only"}; - Configurable nSigmaTofCombinedMaxKa{"nSigmaTofCombinedMaxKa", 5., "Nsigma cut on TOF combined with TPC"}; + // TPC & TOF Combined + Configurable nSigmaCombMaxPi{"nSigmaCombMaxPi", 3., "Nsigma cut on TPC & TOF"}; + Configurable nSigmaCombMaxKa{"nSigmaCombMaxKa", 3., "Nsigma cut on TPC & TOF"}; /// singe track selections Configurable d_pTMin{"d_pTMin", 0.3, "minimum momentum for tracks"}; Configurable d_etaRange{"d_etaRange", 0.8, "eta Range for tracks"}; Configurable d_dcaXYTrackPV{"d_dcaXYTrackPV", 2., "DCA XY of the daughter tracks to the PV"}; Configurable d_dcaZTrackPV{"d_dcaZTrackPV", 10., "DCA Z of the daughter tracks to the PV"}; - /// D0 selections - Configurable applySelectionDoWithTopoConst{"applySelectionDoWithTopoConst", true, "Apply selections on the D0 after constraining it to the PV"}; + // /// D0 selections + Configurable applySelectionDoWithTopoConst{"applySelectionDoWithTopoConst", false, "Apply selections on the D0 after constraining it to the PV"}; Configurable d_pTMinD0{"d_pTMinD0", 0., "minimum momentum for D0 candidates"}; Configurable d_pTMaxD0{"d_pTMaxD0", 36., "maximum momentum for D0 candidates"}; Configurable d_massMinD0{"d_massMinD0", 1.65, "minimum mass for D0"}; @@ -131,6 +122,7 @@ struct qaKFParticle { /// Option to write D0 variables in a tree Configurable d_DwnSmplFact{"d_DwnSmplFact", 1., "Downsampling factor for tree"}; Configurable writeTree{"writeTree", false, "write daughter variables in a tree"}; + Configurable writeHistograms{"writeHistograms", true, "write histograms"}; Configurable writeQAHistograms{"writeQAHistograms", false, "write all QA histograms"}; // Define which track selection should be used: @@ -151,7 +143,7 @@ struct qaKFParticle { ((trackSelection.node() == 4) && requireQualityTracksInFilter()) || ((trackSelection.node() == 5) && requireTrackCutInFilter(TrackSelectionFlags::kInAcceptanceTracks)); - using CollisionTableData = soa::Join; + using CollisionTableData = soa::Join; using BigTracks = soa::Join; using BigTracksExtended = soa::Join; using BigTracksPID = soa::Join; @@ -201,109 +193,80 @@ struct qaKFParticle { } runNumber = 0; - const AxisSpec axisVertexPosX{500, -0.05, 0.05, "X [cm]"}; - const AxisSpec axisVertexPosY{500, -0.05, 0.05, "Y [cm]"}; + const AxisSpec axisVertexPosX{100, -0.05, 0.05, "X [cm]"}; + const AxisSpec axisVertexPosY{100, -0.05, 0.05, "Y [cm]"}; const AxisSpec axisVertexPosZ{100, -20., 20., "Z [cm]"}; const AxisSpec axisVertexNumContrib{160, 0, 160, "Number Of contributors to the PV"}; const AxisSpec axisVertexCov{100, -0.00005, 0.00005}; - const AxisSpec axisParX{300, -0.1, 0.1, "#it{x} [cm]"}; - const AxisSpec axisParY{200, -0.1, 0.1, "#it{y} [cm]"}; - const AxisSpec axisParZ{200, -20., 20., "#it{z} [cm]"}; + const AxisSpec axisParX{100, -0.1, 0.1, "#it{x} [cm]"}; + const AxisSpec axisParY{100, -0.1, 0.1, "#it{y} [cm]"}; + const AxisSpec axisParZ{100, -20., 20., "#it{z} [cm]"}; const AxisSpec axisParPX{binsPt, "#it{p}_{x} [GeV/c]"}; const AxisSpec axisParPY{binsPt, "#it{p}_{y} [GeV/c]"}; const AxisSpec axisParPZ{binsPt, "#it{p}_{z} [GeV/c]"}; - - /// collisions - histos.add("Events/covXX", ";Cov_{xx} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("Events/covXY", ";Cov_{xy} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("Events/covYY", ";Cov_{yy} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("Events/covXZ", ";Cov_{xz} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("Events/covYZ", ";Cov_{yz} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("Events/covZZ", ";Cov_{zz} [cm^{2}]", kTH1D, {axisVertexCov}); - - histos.add("EventsKF/posX", "", kTH1D, {axisVertexPosX}); - histos.add("EventsKF/posY", "", kTH1D, {axisVertexPosY}); - histos.add("EventsKF/posZ", "", kTH1D, {axisVertexPosZ}); - histos.add("EventsKF/posXY", "", kTH2D, {axisVertexPosX, axisVertexPosY}); - histos.add("EventsKF/nContrib", "", kTH1D, {axisVertexNumContrib}); - histos.add("EventsKF/vertexChi2", ";#chi^{2}", kTH1D, {{100, 0, 100}}); - histos.add("EventsKF/covXX", ";Cov_{xx} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("EventsKF/covXY", ";Cov_{xy} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("EventsKF/covYY", ";Cov_{yy} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("EventsKF/covXZ", ";Cov_{xz} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("EventsKF/covYZ", ";Cov_{yz} [cm^{2}]", kTH1D, {axisVertexCov}); - histos.add("EventsKF/covZZ", ";Cov_{zz} [cm^{2}]", kTH1D, {axisVertexCov}); - - histos.add("Tracks/dcaXYToPV", "distance of closest approach in #it{xy} plane;#it{dcaXY} [cm];", kTH1D, {{200, -0.15, 0.15}}); - histos.add("Tracks/dcaZToPV", "distance of closest approach in #it{z} plane;#it{dcaXY} [cm];", kTH1D, {{200, -10., 10.}}); - - /// tracks - histos.add("TracksKFPi/x", "track #it{x} position at dca in local coordinate system", kTH1D, {axisParX}); - histos.add("TracksKFPi/y", "track #it{y} position at dca in local coordinate system", kTH1D, {axisParY}); - histos.add("TracksKFPi/z", "track #it{z} position at dca in local coordinate system", kTH1D, {axisParZ}); - histos.add("TracksKFPi/px", "track #it{p_{x}} momentum at dca in local coordinate system", kTH1D, {axisParPX}); - histos.add("TracksKFPi/py", "track #it{p_{y}} momentum at dca in local coordinate system", kTH1D, {axisParPY}); - histos.add("TracksKFPi/pz", "track #it{p_{z}} momentum at dca in local coordinate system", kTH1D, {axisParPZ}); - histos.add("TracksKFPi/chi2perNDF", "Chi2/NDF of the track;#it{chi2/ndf};", kTH1D, {{200, 0., 25.}}); - histos.add("TracksKFPi/dcaXYToPV", "distance of closest approach in #it{xy} plane;#it{dcaXY} [cm];", kTH1D, {{200, -0.15, 0.15}}); - histos.add("TracksKFPi/dcaToPV", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, -0.15, 0.15}}); - histos.add("TracksKFPi/dcaToPVLargeRange", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); - histos.add("TracksKFPi/dcaToPVLargeRangeMCBeforeReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); - histos.add("TracksKFPi/dcaToPVLargeRangeMCAfterReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); - histos.add("TracksKFPi/nSigmaTPC", "nSigmaTPC vs P", kTH2D, {{axisParPX}, {100, -6., 6.}}); // Add array in function - histos.add("TracksKFPi/deviationPiToSV", "deviation of Pi to SV", kTH1D, {{200, 0., 10.}}); - histos.add("TracksKFPi/deviationPiToPV", "deviation of Pi to PV", kTH1D, {{200, 0., 10.}}); - histos.add("TracksKFPi/dcaToSV", "distance of Pi to SV", kTH1D, {{100, -0.15, 0.15}}); - histos.add("TracksKFPi/dcaXYToSV", "distance xy of Pi to SV", kTH1D, {{100, -0.15, 0.15}}); - - histos.add("TracksKFKa/x", "track #it{x} position at dca in local coordinate system", kTH1D, {axisParX}); - histos.add("TracksKFKa/y", "track #it{y} position at dca in local coordinate system", kTH1D, {axisParY}); - histos.add("TracksKFKa/z", "track #it{z} position at dca in local coordinate system", kTH1D, {axisParZ}); - histos.add("TracksKFKa/px", "track #it{p_{x}} momentum at dca in local coordinate system", kTH1D, {axisParPX}); - histos.add("TracksKFKa/py", "track #it{p_{y}} momentum at dca in local coordinate system", kTH1D, {axisParPY}); - histos.add("TracksKFKa/pz", "track #it{p_{z}} momentum at dca in local coordinate system", kTH1D, {axisParPZ}); - histos.add("TracksKFKa/chi2perNDF", "Chi2/NDF of the track;#it{chi2/ndf};", kTH1D, {{200, 0., 25.}}); - histos.add("TracksKFKa/dcaXYToPV", "distance of closest approach in #it{xy} plane;#it{dcaXY} [cm];", kTH1D, {{200, -0.15, 0.15}}); - histos.add("TracksKFKa/dcaToPV", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, -0.15, 0.15}}); - histos.add("TracksKFKa/dcaToPVLargeRange", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); - histos.add("TracksKFKa/dcaToPVLargeRangeMCBeforeReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); - histos.add("TracksKFKa/dcaToPVLargeRangeMCAfterReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); - histos.add("TracksKFKa/nSigmaTPC", "nSigmaTPC vs P", kTH2D, {{axisParPX}, {100, -6., 6.}}); // Add Array in function - histos.add("TracksKFKa/deviationKaToSV", "deviation of Ka to SV", kTH1D, {{200, 0., 10.}}); - histos.add("TracksKFKa/deviationKaToPV", "deviation of Ka to PV", kTH1D, {{200, 0., 10.}}); - histos.add("TracksKFKa/dcaToSV", "distance of Ka to SV", kTH1D, {{100, -0.15, 0.15}}); - histos.add("TracksKFKa/dcaXYToSV", "distance xy of Ka to SV", kTH1D, {{100, -0.15, 0.15}}); - - histos.add("TracksDaughter/deviationDaugtherTracks", "chi2 in 3D of daughter tracks at the SV", kTH1D, {{200, 0., 0.2}}); - histos.add("TracksDaughter/dcaDaugtherTracks", "distance in 3D of daughter tracks at the SV", kTH1D, {{100, 0., 0.01}}); - histos.add("TracksDaughter/d0pid0ka", "product of impact parameters of daughters to the PV", kTH1D, {{100, -0.003, 0.003}}); - - /// D0 candidates - - histos.add("DZeroCandTopo/X", "X [cm]", kTH1D, {axisParX}); - histos.add("DZeroCandTopo/Y", "Y [cm]", kTH1D, {axisParY}); - histos.add("DZeroCandTopo/Z", "Z [cm]", kTH1D, {axisParZ}); - histos.add("DZeroCandTopo/E", "E", kTH1D, {{100, 0., 50.}}); - histos.add("DZeroCandTopo/Chi2", "Chi2", kTH1D, {{100, 0., 100.}}); - histos.add("DZeroCandTopo/NDF", "NDF", kTH1D, {{5, 0., 5.}}); - histos.add("DZeroCandTopo/Chi2OverNDF", "Chi2OverNDF", kTH1D, {{100, 0., 25.}}); - histos.add("DZeroCandTopo/p", "momentum", kTH1D, {axisParPX}); - histos.add("DZeroCandTopo/pt", "transverse momentum", kTH1D, {axisParPX}); - histos.add("DZeroCandTopo/eta", "eta", kTH1D, {{100, -2., 2.}}); - histos.add("DZeroCandTopo/phi", "phi", kTH1D, {{100, 0., 3.6}}); - histos.add("DZeroCandTopo/mass", "mass", kTH1D, {{430, 1.65, 2.08}}); - histos.add("DZeroCandTopo/massvspt", "mass vs pt", kTH2D, {{axisParPX}, {430, 1.65, 2.08}}); - histos.add("DZeroCandTopo/decayLength", "decay length [cm]", kTH1D, {{200, 0., 0.5}}); - histos.add("DZeroCandTopo/decayLengthXY", "decay length in xy plane [cm]", kTH1D, {{200, 0., 0.5}}); - histos.add("DZeroCandTopo/cosPA", "cosine of pointing angle", kTH1D, {{100, -1, 1.}}); - histos.add("DZeroCandTopo/lifetime", "life time", kTH1D, {{100, 0., 0.2}}); - histos.add("DZeroCandTopo/massErr", "error mass", kTH1D, {{100, 0., 0.1}}); - histos.add("DZeroCandTopo/decayLengthErr", "decay length error [cm]", kTH1D, {{200, 0., 0.05}}); - histos.add("DZeroCandTopo/dcaDToPV", "distance to PV", kTH1D, {{100, -0.01, 0.01}}); - histos.add("DZeroCandTopo/deviationDToPV", "deviation to PV", kTH1D, {{200, 0., 5.}}); - histos.add("DZeroCandTopo/dcaDToPVXY", "distance to PV in xy plane", kTH1D, {{100, -0.02, 0.02}}); - histos.add("DZeroCandTopo/deviationDToPVXY", "deviation to PV in xy plane", kTH1D, {{200, 0., 5.}}); + if (writeHistograms) { + /// tracks + histos.add("TracksKFPi/px", "track #it{p_{x}} momentum at dca in local coordinate system", kTH1D, {axisParPX}); + histos.add("TracksKFPi/py", "track #it{p_{y}} momentum at dca in local coordinate system", kTH1D, {axisParPY}); + histos.add("TracksKFPi/pz", "track #it{p_{z}} momentum at dca in local coordinate system", kTH1D, {axisParPZ}); + histos.add("TracksKFPi/chi2perNDF", "Chi2/NDF of the track;#it{chi2/ndf};", kTH1D, {{100, 0., 25.}}); + histos.add("TracksKFPi/dcaXYToPV", "distance of closest approach in #it{xy} plane;#it{dcaXY} [cm];", kTH1D, {{100, -0.15, 0.15}}); + histos.add("TracksKFPi/dcaToPV", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{100, -0.15, 0.15}}); + histos.add("TracksKFPi/dcaToPVLargeRange", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); + histos.add("TracksKFPi/dcaToPVLargeRangeMCBeforeReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{100, 0., 20.}}); + histos.add("TracksKFPi/dcaToPVLargeRangeMCAfterReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{100, 0., 20.}}); + histos.add("TracksKFPi/nSigmaTPC", "nSigmaTPC vs P", kTH2D, {{axisParPX}, {100, -6., 6.}}); // Add array in function + histos.add("TracksKFPi/deviationPiToSV", "deviation of Pi to SV", kTH1D, {{100, 0., 10.}}); + histos.add("TracksKFPi/deviationPiToPV", "deviation of Pi to PV", kTH1D, {{100, 0., 10.}}); + histos.add("TracksKFPi/dcaToSV", "distance of Pi to SV", kTH1D, {{100, -0.15, 0.15}}); + histos.add("TracksKFPi/dcaXYToSV", "distance xy of Pi to SV", kTH1D, {{100, -0.15, 0.15}}); + + histos.add("TracksKFKa/px", "track #it{p_{x}} momentum at dca in local coordinate system", kTH1D, {axisParPX}); + histos.add("TracksKFKa/py", "track #it{p_{y}} momentum at dca in local coordinate system", kTH1D, {axisParPY}); + histos.add("TracksKFKa/pz", "track #it{p_{z}} momentum at dca in local coordinate system", kTH1D, {axisParPZ}); + histos.add("TracksKFKa/chi2perNDF", "Chi2/NDF of the track;#it{chi2/ndf};", kTH1D, {{100, 0., 25.}}); + histos.add("TracksKFKa/dcaXYToPV", "distance of closest approach in #it{xy} plane;#it{dcaXY} [cm];", kTH1D, {{100, -0.15, 0.15}}); + histos.add("TracksKFKa/dcaToPV", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{100, -0.15, 0.15}}); + histos.add("TracksKFKa/dcaToPVLargeRange", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{200, 0., 20.}}); + histos.add("TracksKFKa/dcaToPVLargeRangeMCBeforeReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{100, 0., 20.}}); + histos.add("TracksKFKa/dcaToPVLargeRangeMCAfterReassignment", "distance of closest approach ;#it{dca} [cm];", kTH1D, {{100, 0., 20.}}); + histos.add("TracksKFKa/nSigmaTPC", "nSigmaTPC vs P", kTH2D, {{axisParPX}, {100, -6., 6.}}); // Add Array in function + histos.add("TracksKFKa/deviationKaToSV", "deviation of Ka to SV", kTH1D, {{100, 0., 10.}}); + histos.add("TracksKFKa/deviationKaToPV", "deviation of Ka to PV", kTH1D, {{100, 0., 10.}}); + histos.add("TracksKFKa/dcaToSV", "distance of Ka to SV", kTH1D, {{100, -0.15, 0.15}}); + histos.add("TracksKFKa/dcaXYToSV", "distance xy of Ka to SV", kTH1D, {{100, -0.15, 0.15}}); + + histos.add("TracksDaughter/deviationDaugtherTracks", "chi2 in 3D of daughter tracks at the SV", kTH1D, {{100, 0., 0.2}}); + histos.add("TracksDaughter/dcaDaugtherTracks", "distance in 3D of daughter tracks at the SV", kTH1D, {{100, 0., 0.01}}); + histos.add("TracksDaughter/d0pid0ka", "product of impact parameters of daughters to the PV", kTH1D, {{100, -0.003, 0.003}}); + + /// D0 candidates + + histos.add("DZeroCandTopo/X", "X [cm]", kTH1D, {axisParX}); + histos.add("DZeroCandTopo/Y", "Y [cm]", kTH1D, {axisParY}); + histos.add("DZeroCandTopo/Z", "Z [cm]", kTH1D, {axisParZ}); + histos.add("DZeroCandTopo/E", "E", kTH1D, {{100, 0., 50.}}); + histos.add("DZeroCandTopo/Chi2", "Chi2", kTH1D, {{100, 0., 100.}}); + histos.add("DZeroCandTopo/NDF", "NDF", kTH1D, {{5, 0., 5.}}); + histos.add("DZeroCandTopo/Chi2OverNDF", "Chi2OverNDF", kTH1D, {{100, 0., 25.}}); + histos.add("DZeroCandTopo/p", "momentum", kTH1D, {axisParPX}); + histos.add("DZeroCandTopo/pt", "transverse momentum", kTH1D, {axisParPX}); + histos.add("DZeroCandTopo/eta", "eta", kTH1D, {{100, -2., 2.}}); + histos.add("DZeroCandTopo/phi", "phi", kTH1D, {{100, 0., 3.6}}); + histos.add("DZeroCandTopo/mass", "mass", kTH1D, {{172, 1.65, 2.08}}); + histos.add("DZeroCandTopo/massvspt", "mass vs pt", kTH2D, {{axisParPX}, {172, 1.65, 2.08}}); + histos.add("DZeroCandTopo/decayLength", "decay length [cm]", kTH1D, {{100, 0., 0.5}}); + histos.add("DZeroCandTopo/decayLengthXY", "decay length in xy plane [cm]", kTH1D, {{100, 0., 0.5}}); + histos.add("DZeroCandTopo/cosPA", "cosine of pointing angle", kTH1D, {{100, -1, 1.}}); + histos.add("DZeroCandTopo/lifetime", "life time", kTH1D, {{100, 0., 0.2}}); + histos.add("DZeroCandTopo/massErr", "error mass", kTH1D, {{100, 0., 0.1}}); + histos.add("DZeroCandTopo/decayLengthErr", "decay length error [cm]", kTH1D, {{100, 0., 0.05}}); + histos.add("DZeroCandTopo/dcaDToPV", "distance to PV", kTH1D, {{100, -0.01, 0.01}}); + histos.add("DZeroCandTopo/deviationDToPV", "deviation to PV", kTH1D, {{100, 0., 5.}}); + histos.add("DZeroCandTopo/dcaDToPVXY", "distance to PV in xy plane", kTH1D, {{100, -0.02, 0.02}}); + histos.add("DZeroCandTopo/deviationDToPVXY", "deviation to PV in xy plane", kTH1D, {{100, 0., 5.}}); + } auto hSelection = histos.add("DZeroCandTopo/Selections", "Selections", kTH1D, {{25, 0.5, 25.5}}); hSelection->GetXaxis()->SetBinLabel(hSelection->FindBin(1), "All Collisions"); hSelection->GetXaxis()->SetBinLabel(hSelection->FindBin(2), "Cov < 0"); @@ -349,36 +312,18 @@ struct qaKFParticle { histos.add("DZeroCandGeo/pt", "transverse momentum", kTH1D, {axisParPX}); histos.add("DZeroCandGeo/eta", "eta", kTH1D, {{100, -2., 2.}}); histos.add("DZeroCandGeo/phi", "phi", kTH1D, {{100, 0., 3.6}}); - histos.add("DZeroCandGeo/mass", "mass", kTH1D, {{430, 1.65, 2.08}}); - histos.add("DZeroCandGeo/massvspt", "mass vs pt", kTH2D, {{axisParPX}, {430, 1.65, 2.08}}); + histos.add("DZeroCandGeo/mass", "mass", kTH1D, {{172, 1.65, 2.08}}); + histos.add("DZeroCandGeo/massvspt", "mass vs pt", kTH2D, {{axisParPX}, {172, 1.65, 2.08}}); histos.add("DZeroCandGeo/cosPA", "cosine of pointing angle", kTH1D, {{100, -1, 1.}}); histos.add("DZeroCandGeo/massErr", "error mass", kTH1D, {{100, 0., 0.1}}); histos.add("DZeroCandGeo/dcaDToPV", "distance to PV", kTH1D, {{100, -0.01, 0.01}}); - histos.add("DZeroCandGeo/deviationDToPV", "deviation to PV", kTH1D, {{200, 0., 5.}}); + histos.add("DZeroCandGeo/deviationDToPV", "deviation to PV", kTH1D, {{100, 0., 5.}}); histos.add("DZeroCandGeo/dcaDToPVXY", "distance to PV in xy plane", kTH1D, {{100, -0.02, 0.02}}); - histos.add("DZeroCandGeo/deviationDToPVXY", "deviation to PV in xy plane", kTH1D, {{200, 0., 5.}}); + histos.add("DZeroCandGeo/deviationDToPVXY", "deviation to PV in xy plane", kTH1D, {{100, 0., 5.}}); histos.add("DZeroCandGeo/cosThetaStar", "cosine theta star", kTH1D, {{100, -1., 1.}}); - - histos.add("DZeroCandTopoAtSV/X", "X [cm]", kTH1D, {axisParX}); - histos.add("DZeroCandTopoAtSV/Y", "Y [cm]", kTH1D, {axisParY}); - histos.add("DZeroCandTopoAtSV/Z", "Z [cm]", kTH1D, {axisParZ}); - histos.add("DZeroCandTopoAtSV/E", "E", kTH1D, {{100, 0., 50.}}); - histos.add("DZeroCandTopoAtSV/Chi2", "Chi2", kTH1D, {{100, 0., 100.}}); - histos.add("DZeroCandTopoAtSV/NDF", "NDF", kTH1D, {{5, 0., 5.}}); - histos.add("DZeroCandTopoAtSV/Chi2OverNDF", "Chi2OverNDF", kTH1D, {{100, 0., 25.}}); - histos.add("DZeroCandTopoAtSV/p", "momentum", kTH1D, {axisParPX}); - histos.add("DZeroCandTopoAtSV/pt", "transverse momentum", kTH1D, {axisParPX}); - histos.add("DZeroCandTopoAtSV/eta", "eta", kTH1D, {{100, -2., 2.}}); - histos.add("DZeroCandTopoAtSV/phi", "phi", kTH1D, {{100, 0., 3.6}}); - histos.add("DZeroCandTopoAtSV/mass", "mass", kTH1D, {{430, 1.65, 2.08}}); - histos.add("DZeroCandTopoAtSV/massvspt", "mass vs pt", kTH2D, {{axisParPX}, {430, 1.65, 2.08}}); - histos.add("DZeroCandTopoAtSV/cosPA", "cosine of pointing angle", kTH1D, {{100, -1, 1.}}); - histos.add("DZeroCandTopoAtSV/massErr", "error mass", kTH1D, {{100, 0., 0.1}}); - histos.add("DZeroCandTopoAtSV/dcaDToPV", "distance to PV", kTH1D, {{100, -0.01, 0.01}}); - histos.add("DZeroCandTopoAtSV/deviationDToPV", "deviation to PV", kTH1D, {{200, 0., 5.}}); - histos.add("DZeroCandTopoAtSV/dcaDToPVXY", "distance to PV in xy plane", kTH1D, {{100, -0.02, 0.02}}); - histos.add("DZeroCandTopoAtSV/deviationDToPVXY", "deviation to PV in xy plane", kTH1D, {{200, 0., 5.}}); } + LOGF(info, "End of init"); + histos.print(); } /// End init /// Function to select collisions @@ -541,69 +486,114 @@ struct qaKFParticle { return true; } + template + bool SelectPIDCombined(const T1& track, int particle) + { + switch (particle) { + case kPiPlus: { + if ((track.pt() <= ptPidTofMinPi) && track.hasTPC() && (abs(track.tpcNSigmaPi()) < nSigmaTpcMaxPi)) { + return true; + } else if ((track.pt() > ptPidTofMinPi) && track.hasTPC() && !track.hasTOF() && (abs(track.tpcNSigmaPi()) < nSigmaTpcMaxPi)) { + return true; + } else if ((track.pt() > ptPidTofMinPi) && !track.hasTPC() && track.hasTOF() && (abs(track.tofNSigmaPi()) < nSigmaTofMaxPi)) { + return true; + } else if ((track.pt() > ptPidTofMinPi) && track.hasTPC() && track.hasTOF()) { + float CombinednSigma = 1. / sqrt(2) * sqrt((track.tpcNSigmaPi() * track.tpcNSigmaPi()) + (track.tofNSigmaPi() * track.tofNSigmaPi())); + if (abs(CombinednSigma) < nSigmaCombMaxPi) { + return true; + } else { + return false; + } + } else { + return false; + } + break; + } + case kKPlus: { + if ((track.pt() <= ptPidTofMinKa) && track.hasTPC() && (abs(track.tpcNSigmaKa()) < nSigmaTpcMaxKa)) { + return true; + } else if ((track.pt() > ptPidTofMinKa) && track.hasTPC() && !track.hasTOF() && (abs(track.tpcNSigmaKa()) < nSigmaTpcMaxKa)) { + return true; + } else if ((track.pt() > ptPidTofMinKa) && !track.hasTPC() && track.hasTOF() && (abs(track.tofNSigmaKa()) < nSigmaTofMaxKa)) { + return true; + } else if ((track.pt() > ptPidTofMinKa) && track.hasTPC() && track.hasTOF()) { + float CombinednSigma = 1. / sqrt(2) * sqrt((track.tpcNSigmaKa() * track.tpcNSigmaKa()) + (track.tofNSigmaKa() * track.tofNSigmaKa())); + if (abs(CombinednSigma) < nSigmaCombMaxKa) { + return true; + } else { + return false; + } + } else { + return false; + } + break; + } + default: { + LOGF(error, "ERROR: Species is not implemented"); + return false; + } + } + } + template void fillHistograms(const T1& kfpTrackPi, const T1& kfpTrackKa, const T2& KFPion, const T2& KFKaon, const T2& KFDZero_PV, const T2& KFDZero, const T2& KFPV, const T2& KFDZero_DecayVtx, float cosThetaStar) { - /// fill daughter track parameters - histos.fill(HIST("TracksKFPi/x"), kfpTrackPi.GetX()); - histos.fill(HIST("TracksKFPi/y"), kfpTrackPi.GetY()); - histos.fill(HIST("TracksKFPi/z"), kfpTrackPi.GetZ()); - histos.fill(HIST("TracksKFPi/px"), kfpTrackPi.GetPx()); - histos.fill(HIST("TracksKFPi/py"), kfpTrackPi.GetPy()); - histos.fill(HIST("TracksKFPi/pz"), kfpTrackPi.GetPz()); - histos.fill(HIST("TracksKFPi/chi2perNDF"), kfpTrackPi.GetChi2perNDF()); - histos.fill(HIST("TracksKFPi/dcaXYToPV"), KFPion.GetDistanceFromVertexXY(KFPV)); - histos.fill(HIST("TracksKFPi/dcaToPV"), KFPion.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("TracksKFPi/dcaToPVLargeRange"), KFPion.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("TracksKFPi/deviationPiToSV"), KFPion.GetDeviationFromVertex(KFDZero)); - histos.fill(HIST("TracksKFPi/deviationPiToPV"), KFPion.GetDeviationFromVertex(KFPV)); - histos.fill(HIST("TracksKFPi/dcaToSV"), KFPion.GetDistanceFromVertex(KFDZero)); - histos.fill(HIST("TracksKFPi/dcaXYToSV"), KFPion.GetDistanceFromVertexXY(KFDZero)); - - histos.fill(HIST("TracksKFKa/x"), kfpTrackKa.GetX()); - histos.fill(HIST("TracksKFKa/y"), kfpTrackKa.GetY()); - histos.fill(HIST("TracksKFKa/z"), kfpTrackKa.GetZ()); - histos.fill(HIST("TracksKFKa/px"), kfpTrackKa.GetPx()); - histos.fill(HIST("TracksKFKa/py"), kfpTrackKa.GetPy()); - histos.fill(HIST("TracksKFKa/pz"), kfpTrackKa.GetPz()); - histos.fill(HIST("TracksKFKa/chi2perNDF"), kfpTrackKa.GetChi2perNDF()); - histos.fill(HIST("TracksKFKa/dcaXYToPV"), KFKaon.GetDistanceFromVertexXY(KFPV)); - histos.fill(HIST("TracksKFKa/dcaToPV"), KFKaon.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("TracksKFKa/dcaToPVLargeRange"), KFPion.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("TracksKFKa/deviationKaToSV"), KFKaon.GetDeviationFromVertex(KFDZero)); - histos.fill(HIST("TracksKFKa/deviationKaToPV"), KFKaon.GetDeviationFromVertex(KFPV)); - histos.fill(HIST("TracksKFKa/dcaToSV"), KFKaon.GetDistanceFromVertex(KFDZero)); - histos.fill(HIST("TracksKFKa/dcaXYToSV"), KFKaon.GetDistanceFromVertexXY(KFDZero)); - - histos.fill(HIST("TracksDaughter/deviationDaugtherTracks"), KFPion.GetDeviationFromParticle(KFKaon)); - histos.fill(HIST("TracksDaughter/dcaDaugtherTracks"), KFPion.GetDistanceFromParticle(KFKaon)); - float d0pid0ka = KFPion.GetDistanceFromVertexXY(KFPV) * KFKaon.GetDistanceFromVertexXY(KFPV); - histos.fill(HIST("TracksDaughter/d0pid0ka"), d0pid0ka); - - histos.fill(HIST("DZeroCandTopo/X"), KFDZero_PV.GetX()); - histos.fill(HIST("DZeroCandTopo/Y"), KFDZero_PV.GetY()); - histos.fill(HIST("DZeroCandTopo/Z"), KFDZero_PV.GetZ()); - histos.fill(HIST("DZeroCandTopo/E"), KFDZero_PV.GetE()); - histos.fill(HIST("DZeroCandTopo/Chi2"), KFDZero_PV.GetChi2()); - histos.fill(HIST("DZeroCandTopo/NDF"), KFDZero_PV.GetNDF()); - float chi2topo = KFDZero_PV.GetChi2() / KFDZero_PV.GetNDF(); - histos.fill(HIST("DZeroCandTopo/Chi2OverNDF"), chi2topo); - histos.fill(HIST("DZeroCandTopo/p"), KFDZero_PV.GetP()); - histos.fill(HIST("DZeroCandTopo/pt"), KFDZero_PV.GetPt()); - histos.fill(HIST("DZeroCandTopo/eta"), KFDZero_PV.GetEta()); - histos.fill(HIST("DZeroCandTopo/phi"), KFDZero_PV.GetPhi()); - histos.fill(HIST("DZeroCandTopo/mass"), KFDZero_PV.GetMass()); - histos.fill(HIST("DZeroCandTopo/massvspt"), KFDZero_PV.GetPt(), KFDZero_PV.GetMass()); - histos.fill(HIST("DZeroCandTopo/decayLength"), KFDZero_PV.GetDecayLength()); - histos.fill(HIST("DZeroCandTopo/decayLengthXY"), KFDZero_PV.GetDecayLengthXY()); - histos.fill(HIST("DZeroCandTopo/cosPA"), cpaFromKF(KFDZero_DecayVtx, KFPV)); - histos.fill(HIST("DZeroCandTopo/lifetime"), KFDZero_PV.GetLifeTime()); - histos.fill(HIST("DZeroCandTopo/massErr"), KFDZero_PV.GetErrMass()); - histos.fill(HIST("DZeroCandTopo/decayLengthErr"), KFDZero_PV.GetErrDecayLength()); - histos.fill(HIST("DZeroCandTopo/dcaDToPV"), KFDZero_PV.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("DZeroCandTopo/deviationDToPV"), KFDZero_PV.GetDeviationFromVertex(KFPV)); - histos.fill(HIST("DZeroCandTopo/dcaDToPVXY"), KFDZero_PV.GetDistanceFromVertexXY(KFPV)); - histos.fill(HIST("DZeroCandTopo/deviationDToPVXY"), KFDZero_PV.GetDeviationFromVertexXY(KFPV)); + if (writeHistograms) { + /// fill daughter track parameters + histos.fill(HIST("TracksKFPi/px"), kfpTrackPi.GetPx()); + histos.fill(HIST("TracksKFPi/py"), kfpTrackPi.GetPy()); + histos.fill(HIST("TracksKFPi/pz"), kfpTrackPi.GetPz()); + histos.fill(HIST("TracksKFPi/chi2perNDF"), kfpTrackPi.GetChi2perNDF()); + histos.fill(HIST("TracksKFPi/dcaXYToPV"), KFPion.GetDistanceFromVertexXY(KFPV)); + histos.fill(HIST("TracksKFPi/dcaToPV"), KFPion.GetDistanceFromVertex(KFPV)); + histos.fill(HIST("TracksKFPi/dcaToPVLargeRange"), KFPion.GetDistanceFromVertex(KFPV)); + histos.fill(HIST("TracksKFPi/deviationPiToSV"), KFPion.GetDeviationFromVertex(KFDZero)); + histos.fill(HIST("TracksKFPi/deviationPiToPV"), KFPion.GetDeviationFromVertex(KFPV)); + histos.fill(HIST("TracksKFPi/dcaToSV"), KFPion.GetDistanceFromVertex(KFDZero)); + histos.fill(HIST("TracksKFPi/dcaXYToSV"), KFPion.GetDistanceFromVertexXY(KFDZero)); + + histos.fill(HIST("TracksKFKa/px"), kfpTrackKa.GetPx()); + histos.fill(HIST("TracksKFKa/py"), kfpTrackKa.GetPy()); + histos.fill(HIST("TracksKFKa/pz"), kfpTrackKa.GetPz()); + histos.fill(HIST("TracksKFKa/chi2perNDF"), kfpTrackKa.GetChi2perNDF()); + histos.fill(HIST("TracksKFKa/dcaXYToPV"), KFKaon.GetDistanceFromVertexXY(KFPV)); + histos.fill(HIST("TracksKFKa/dcaToPV"), KFKaon.GetDistanceFromVertex(KFPV)); + histos.fill(HIST("TracksKFKa/dcaToPVLargeRange"), KFPion.GetDistanceFromVertex(KFPV)); + histos.fill(HIST("TracksKFKa/deviationKaToSV"), KFKaon.GetDeviationFromVertex(KFDZero)); + histos.fill(HIST("TracksKFKa/deviationKaToPV"), KFKaon.GetDeviationFromVertex(KFPV)); + histos.fill(HIST("TracksKFKa/dcaToSV"), KFKaon.GetDistanceFromVertex(KFDZero)); + histos.fill(HIST("TracksKFKa/dcaXYToSV"), KFKaon.GetDistanceFromVertexXY(KFDZero)); + + histos.fill(HIST("TracksDaughter/deviationDaugtherTracks"), KFPion.GetDeviationFromParticle(KFKaon)); + histos.fill(HIST("TracksDaughter/dcaDaugtherTracks"), KFPion.GetDistanceFromParticle(KFKaon)); + float d0pid0ka = KFPion.GetDistanceFromVertexXY(KFPV) * KFKaon.GetDistanceFromVertexXY(KFPV); + histos.fill(HIST("TracksDaughter/d0pid0ka"), d0pid0ka); + + histos.fill(HIST("DZeroCandTopo/X"), KFDZero_PV.GetX()); + histos.fill(HIST("DZeroCandTopo/Y"), KFDZero_PV.GetY()); + histos.fill(HIST("DZeroCandTopo/Z"), KFDZero_PV.GetZ()); + histos.fill(HIST("DZeroCandTopo/E"), KFDZero_PV.GetE()); + histos.fill(HIST("DZeroCandTopo/Chi2"), KFDZero_PV.GetChi2()); + histos.fill(HIST("DZeroCandTopo/NDF"), KFDZero_PV.GetNDF()); + float chi2topo = KFDZero_PV.GetChi2() / KFDZero_PV.GetNDF(); + histos.fill(HIST("DZeroCandTopo/Chi2OverNDF"), chi2topo); + histos.fill(HIST("DZeroCandTopo/p"), KFDZero_PV.GetP()); + histos.fill(HIST("DZeroCandTopo/pt"), KFDZero_PV.GetPt()); + histos.fill(HIST("DZeroCandTopo/eta"), KFDZero_PV.GetEta()); + histos.fill(HIST("DZeroCandTopo/phi"), KFDZero_PV.GetPhi()); + histos.fill(HIST("DZeroCandTopo/mass"), KFDZero_PV.GetMass()); + histos.fill(HIST("DZeroCandTopo/massvspt"), KFDZero_PV.GetPt(), KFDZero_PV.GetMass()); + histos.fill(HIST("DZeroCandTopo/decayLength"), KFDZero_PV.GetDecayLength()); + histos.fill(HIST("DZeroCandTopo/decayLengthXY"), KFDZero_PV.GetDecayLengthXY()); + histos.fill(HIST("DZeroCandTopo/cosPA"), cpaFromKF(KFDZero_DecayVtx, KFPV)); + histos.fill(HIST("DZeroCandTopo/lifetime"), KFDZero_PV.GetLifeTime()); + histos.fill(HIST("DZeroCandTopo/massErr"), KFDZero_PV.GetErrMass()); + histos.fill(HIST("DZeroCandTopo/decayLengthErr"), KFDZero_PV.GetErrDecayLength()); + histos.fill(HIST("DZeroCandTopo/dcaDToPV"), KFDZero_PV.GetDistanceFromVertex(KFPV)); + histos.fill(HIST("DZeroCandTopo/deviationDToPV"), KFDZero_PV.GetDeviationFromVertex(KFPV)); + histos.fill(HIST("DZeroCandTopo/dcaDToPVXY"), KFDZero_PV.GetDistanceFromVertexXY(KFPV)); + histos.fill(HIST("DZeroCandTopo/deviationDToPVXY"), KFDZero_PV.GetDeviationFromVertexXY(KFPV)); + } if (writeQAHistograms) { histos.fill(HIST("DZeroCandGeo/X"), KFDZero.GetX()); @@ -627,26 +617,6 @@ struct qaKFParticle { histos.fill(HIST("DZeroCandGeo/dcaDToPVXY"), KFDZero.GetDistanceFromVertexXY(KFPV)); histos.fill(HIST("DZeroCandGeo/deviationDToPVXY"), KFDZero.GetDeviationFromVertexXY(KFPV)); histos.fill(HIST("DZeroCandGeo/cosThetaStar"), cosThetaStar); - - histos.fill(HIST("DZeroCandTopoAtSV/X"), KFDZero_DecayVtx.GetX()); - histos.fill(HIST("DZeroCandTopoAtSV/Y"), KFDZero_DecayVtx.GetY()); - histos.fill(HIST("DZeroCandTopoAtSV/Z"), KFDZero_DecayVtx.GetZ()); - histos.fill(HIST("DZeroCandTopoAtSV/E"), KFDZero_DecayVtx.GetE()); - histos.fill(HIST("DZeroCandTopoAtSV/Chi2"), KFDZero_DecayVtx.GetChi2()); - histos.fill(HIST("DZeroCandTopoAtSV/NDF"), KFDZero_DecayVtx.GetNDF()); - float chi2Topo = KFDZero_DecayVtx.GetChi2() / KFDZero_DecayVtx.GetNDF(); - histos.fill(HIST("DZeroCandTopoAtSV/Chi2OverNDF"), chi2Topo); - histos.fill(HIST("DZeroCandTopoAtSV/p"), KFDZero_DecayVtx.GetP()); - histos.fill(HIST("DZeroCandTopoAtSV/pt"), KFDZero_DecayVtx.GetPt()); - histos.fill(HIST("DZeroCandTopoAtSV/eta"), KFDZero_DecayVtx.GetEta()); - histos.fill(HIST("DZeroCandTopoAtSV/phi"), KFDZero_DecayVtx.GetPhi()); - histos.fill(HIST("DZeroCandTopoAtSV/mass"), KFDZero_DecayVtx.GetMass()); - histos.fill(HIST("DZeroCandTopoAtSV/massvspt"), KFDZero_DecayVtx.GetPt(), KFDZero_DecayVtx.GetMass()); - histos.fill(HIST("DZeroCandTopoAtSV/massErr"), KFDZero_DecayVtx.GetErrMass()); - histos.fill(HIST("DZeroCandTopoAtSV/dcaDToPV"), KFDZero_DecayVtx.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("DZeroCandTopoAtSV/deviationDToPV"), KFDZero_DecayVtx.GetDeviationFromVertex(KFPV)); - histos.fill(HIST("DZeroCandTopoAtSV/dcaDToPVXY"), KFDZero_DecayVtx.GetDistanceFromVertexXY(KFPV)); - histos.fill(HIST("DZeroCandTopoAtSV/deviationDToPVXY"), KFDZero_DecayVtx.GetDeviationFromVertexXY(KFPV)); } } template @@ -703,7 +673,6 @@ struct qaKFParticle { /// Process function for data void processData(CollisionTableData::iterator const& collision, soa::Filtered const& tracks, aod::BCsWithTimestamps const&) { - auto bc = collision.bc_as(); if (runNumber != bc.runNumber()) { initMagneticFieldCCDB(bc, runNumber, ccdb, isRun3 ? ccdbPathGrpMag : ccdbPathGrp, lut, isRun3); @@ -713,13 +682,8 @@ struct qaKFParticle { KFParticle::SetField(magneticField); #endif } + histos.fill(HIST("DZeroCandTopo/Selections"), 1.f); - histos.fill(HIST("Events/covXX"), collision.covXX()); - histos.fill(HIST("Events/covXY"), collision.covXY()); - histos.fill(HIST("Events/covXZ"), collision.covXZ()); - histos.fill(HIST("Events/covYY"), collision.covYY()); - histos.fill(HIST("Events/covYZ"), collision.covYZ()); - histos.fill(HIST("Events/covZZ"), collision.covZZ()); /// Apply event selection if (!isSelectedCollision(collision)) { return; @@ -727,37 +691,6 @@ struct qaKFParticle { /// set KF primary vertex KFPVertex kfpVertex = createKFPVertexFromCollision(collision); KFParticle KFPV(kfpVertex); - - /// fill collision parameters - histos.fill(HIST("EventsKF/posX"), kfpVertex.GetX()); - histos.fill(HIST("EventsKF/posY"), kfpVertex.GetY()); - histos.fill(HIST("EventsKF/posZ"), kfpVertex.GetZ()); - histos.fill(HIST("EventsKF/posXY"), kfpVertex.GetX(), kfpVertex.GetY()); - histos.fill(HIST("EventsKF/nContrib"), kfpVertex.GetNContributors()); - histos.fill(HIST("EventsKF/vertexChi2"), kfpVertex.GetChi2()); - histos.fill(HIST("EventsKF/covXX"), kfpVertex.GetCovariance(0)); - histos.fill(HIST("EventsKF/covXY"), kfpVertex.GetCovariance(1)); - histos.fill(HIST("EventsKF/covYY"), kfpVertex.GetCovariance(2)); - histos.fill(HIST("EventsKF/covXZ"), kfpVertex.GetCovariance(3)); - histos.fill(HIST("EventsKF/covYZ"), kfpVertex.GetCovariance(4)); - histos.fill(HIST("EventsKF/covZZ"), kfpVertex.GetCovariance(5)); - - TrackSelectorPID selectorPion(kPiPlus); - selectorPion.setRangePtTPC(ptPidTpcMinPi, ptPidTpcMaxPi); - selectorPion.setRangeNSigmaTPC(-nSigmaTpcMaxPi, nSigmaTpcMaxPi); - selectorPion.setRangeNSigmaTPCCondTOF(-nSigmaTpcCombinedMaxPi, nSigmaTpcCombinedMaxPi); - selectorPion.setRangePtTOF(ptPidTofMinPi, ptPidTofMaxPi); - selectorPion.setRangeNSigmaTOF(-nSigmaTofMaxPi, nSigmaTofMaxPi); - selectorPion.setRangeNSigmaTOFCondTPC(-nSigmaTofCombinedMaxPi, nSigmaTofCombinedMaxPi); - - TrackSelectorPID selectorKaon(kKPlus); - selectorPion.setRangePtTPC(ptPidTpcMinKa, ptPidTpcMaxKa); - selectorPion.setRangeNSigmaTPC(-nSigmaTpcMaxKa, nSigmaTpcMaxKa); - selectorPion.setRangeNSigmaTPCCondTOF(-nSigmaTpcCombinedMaxKa, nSigmaTpcCombinedMaxKa); - selectorPion.setRangePtTOF(ptPidTofMinKa, ptPidTofMaxKa); - selectorPion.setRangeNSigmaTOF(-nSigmaTofMaxKa, nSigmaTofMaxKa); - selectorPion.setRangeNSigmaTOFCondTPC(-nSigmaTofCombinedMaxKa, nSigmaTofCombinedMaxKa); - for (auto& [track1, track2] : combinations(soa::CombinationsStrictlyUpperIndexPolicy(tracks, tracks))) { histos.fill(HIST("DZeroCandTopo/Selections"), 3.f); @@ -768,11 +701,6 @@ struct qaKFParticle { continue; } - histos.fill(HIST("Tracks/dcaXYToPV"), track1.dcaXY()); - histos.fill(HIST("Tracks/dcaXYToPV"), track2.dcaXY()); - histos.fill(HIST("Tracks/dcaZToPV"), track1.dcaZ()); - histos.fill(HIST("Tracks/dcaZToPV"), track2.dcaZ()); - KFPTrack kfpTrackPosPi; KFPTrack kfpTrackNegPi; KFPTrack kfpTrackPosKa; @@ -791,10 +719,10 @@ struct qaKFParticle { float TOFnSigmaNegKa = 0; int source = 0; - int pidKaonTr1 = selectorKaon.getStatusTrackPIDTpcOrTof(track1); - int pidKaonTr2 = selectorKaon.getStatusTrackPIDTpcOrTof(track2); - int pidPionTr1 = selectorPion.getStatusTrackPIDTpcOrTof(track1); - int pidPionTr2 = selectorPion.getStatusTrackPIDTpcOrTof(track2); + bool pidKaonTr1 = SelectPIDCombined(track1, kKPlus); + bool pidKaonTr2 = SelectPIDCombined(track2, kKPlus); + bool pidPionTr1 = SelectPIDCombined(track1, kPiPlus); + bool pidPionTr2 = SelectPIDCombined(track2, kPiPlus); if (track1.isPVContributor() || track2.isPVContributor()) { PVContributor = 1; @@ -805,7 +733,7 @@ struct qaKFParticle { } /// Select D0 and D0bar candidates - if (pidPionTr1 == TrackSelectorPID::Status::PIDAccepted && pidKaonTr2 == TrackSelectorPID::Status::PIDAccepted) { + if (pidPionTr1 && pidKaonTr2) { if (track1.sign() == 1 && track2.sign() == -1) { CandD0 = true; source = 1; @@ -828,10 +756,13 @@ struct qaKFParticle { continue; } } - if (pidKaonTr1 == TrackSelectorPID::Status::PIDAccepted && pidPionTr2 == TrackSelectorPID::Status::PIDAccepted) { + if (pidKaonTr1 && pidPionTr2) { if (track1.sign() == 1 && track2.sign() == -1) { CandD0bar = true; source = 2; + if (CandD0 == true) { + source = 3; + } kfpTrackNegPi = createKFPTrackFromTrack(track2); kfpTrackPosKa = createKFPTrackFromTrack(track1); TPCnSigmaNegPi = track2.tpcNSigmaPi(); @@ -841,6 +772,9 @@ struct qaKFParticle { } else if (track1.sign() == -1 && track2.sign() == 1) { CandD0 = true; source = 1; + if (CandD0bar == true) { + source = 3; + } kfpTrackPosPi = createKFPTrackFromTrack(track2); kfpTrackNegKa = createKFPTrackFromTrack(track1); TPCnSigmaPosPi = track2.tpcNSigmaPi(); @@ -855,9 +789,6 @@ struct qaKFParticle { histos.fill(HIST("DZeroCandTopo/Selections"), 9.f); continue; } - if (CandD0 && CandD0bar) { - source = 3; - } KFParticle KFPosPion(kfpTrackPosPi, 211); KFParticle KFNegPion(kfpTrackNegPi, 211); @@ -895,10 +826,11 @@ struct qaKFParticle { continue; } } - fillHistograms(kfpTrackPosPi, kfpTrackNegKa, KFPosPion, KFNegKaon, KFDZero_PV, KFDZero, KFPV, KFDZero_DecayVtx, cosThetaStar); - histos.fill(HIST("TracksKFPi/nSigmaTPC"), kfpTrackPosPi.GetPt(), TPCnSigmaPosPi); - histos.fill(HIST("TracksKFKa/nSigmaTPC"), kfpTrackNegKa.GetPt(), TPCnSigmaNegKa); + if (writeHistograms) { + histos.fill(HIST("TracksKFPi/nSigmaTPC"), kfpTrackPosPi.GetPt(), TPCnSigmaPosPi); + histos.fill(HIST("TracksKFKa/nSigmaTPC"), kfpTrackNegKa.GetPt(), TPCnSigmaNegKa); + } writeVarTree(kfpTrackPosPi, kfpTrackNegKa, KFPosPion, KFNegKaon, KFDZero_PV, KFDZero, KFPV, KFDZero_DecayVtx, TPCnSigmaPosPi, TOFnSigmaPosPi, TPCnSigmaNegKa, TOFnSigmaNegKa, cosThetaStar, track1, source); } if (CandD0bar) { @@ -930,8 +862,10 @@ struct qaKFParticle { } } fillHistograms(kfpTrackNegPi, kfpTrackPosKa, KFNegPion, KFPosKaon, KFDZeroBar_PV, KFDZeroBar, KFPV, KFDZeroBar_DecayVtx, cosThetaStar); - histos.fill(HIST("TracksKFPi/nSigmaTPC"), kfpTrackNegPi.GetPt(), TPCnSigmaNegPi); - histos.fill(HIST("TracksKFKa/nSigmaTPC"), kfpTrackPosKa.GetPt(), TPCnSigmaPosKa); + if (writeHistograms) { + histos.fill(HIST("TracksKFPi/nSigmaTPC"), kfpTrackNegPi.GetPt(), TPCnSigmaNegPi); + histos.fill(HIST("TracksKFKa/nSigmaTPC"), kfpTrackPosKa.GetPt(), TPCnSigmaPosKa); + } writeVarTree(kfpTrackNegPi, kfpTrackPosKa, KFNegPion, KFPosKaon, KFDZeroBar_PV, KFDZeroBar, KFPV, KFDZeroBar_DecayVtx, TPCnSigmaNegPi, TOFnSigmaNegPi, TPCnSigmaPosKa, TOFnSigmaPosKa, cosThetaStar, track1, source); } } @@ -940,9 +874,9 @@ struct qaKFParticle { /// Process function for MC using CollisionTableMC = soa::Join; - using CollisionTableDataMult = soa::Join; + using CollisionTableDataMult = soa::Join; using TrackTableMC = soa::Join; - Preslice perMcCollision = aod::mccollisionlabel::mcCollisionId; + // Preslice perMcCollision = o2::aod::mccollisionlabel::mcCollisionId; void processMC(CollisionTableMC::iterator const& collision, CollisionTableMC const& collisions, soa::Filtered const& tracks, aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions, aod::BCsWithTimestamps const&) { auto bc = collision.bc_as(); @@ -959,12 +893,6 @@ struct qaKFParticle { return; } histos.fill(HIST("DZeroCandTopo/Selections"), 1.f); - histos.fill(HIST("Events/covXX"), collision.covXX()); - histos.fill(HIST("Events/covXY"), collision.covXY()); - histos.fill(HIST("Events/covXZ"), collision.covXZ()); - histos.fill(HIST("Events/covYY"), collision.covYY()); - histos.fill(HIST("Events/covYZ"), collision.covYZ()); - histos.fill(HIST("Events/covZZ"), collision.covZZ()); /// Apply event selection if (!isSelectedCollision(collision)) { return; @@ -976,36 +904,6 @@ struct qaKFParticle { KFPVertex kfpVertexDefault = createKFPVertexFromCollision(collision); KFParticle KFPVDefault(kfpVertexDefault); - /// fill collision parameters - histos.fill(HIST("EventsKF/posX"), kfpVertex.GetX()); - histos.fill(HIST("EventsKF/posY"), kfpVertex.GetY()); - histos.fill(HIST("EventsKF/posZ"), kfpVertex.GetZ()); - histos.fill(HIST("EventsKF/posXY"), kfpVertex.GetX(), kfpVertex.GetY()); - histos.fill(HIST("EventsKF/nContrib"), kfpVertex.GetNContributors()); - histos.fill(HIST("EventsKF/vertexChi2"), kfpVertex.GetChi2()); - histos.fill(HIST("EventsKF/covXX"), kfpVertex.GetCovariance(0)); - histos.fill(HIST("EventsKF/covXY"), kfpVertex.GetCovariance(1)); - histos.fill(HIST("EventsKF/covYY"), kfpVertex.GetCovariance(2)); - histos.fill(HIST("EventsKF/covXZ"), kfpVertex.GetCovariance(3)); - histos.fill(HIST("EventsKF/covYZ"), kfpVertex.GetCovariance(4)); - histos.fill(HIST("EventsKF/covZZ"), kfpVertex.GetCovariance(5)); - - TrackSelectorPID selectorPion(kPiPlus); - selectorPion.setRangePtTPC(ptPidTpcMinPi, ptPidTpcMaxPi); - selectorPion.setRangeNSigmaTPC(-nSigmaTpcMaxPi, nSigmaTpcMaxPi); - selectorPion.setRangeNSigmaTPCCondTOF(-nSigmaTpcCombinedMaxPi, nSigmaTpcCombinedMaxPi); - selectorPion.setRangePtTOF(ptPidTofMinPi, ptPidTofMaxPi); - selectorPion.setRangeNSigmaTOF(-nSigmaTofMaxPi, nSigmaTofMaxPi); - selectorPion.setRangeNSigmaTOFCondTPC(-nSigmaTofCombinedMaxPi, nSigmaTofCombinedMaxPi); - - TrackSelectorPID selectorKaon(kKPlus); - selectorPion.setRangePtTPC(ptPidTpcMinKa, ptPidTpcMaxKa); - selectorPion.setRangeNSigmaTPC(-nSigmaTpcMaxKa, nSigmaTpcMaxKa); - selectorPion.setRangeNSigmaTPCCondTOF(-nSigmaTpcCombinedMaxKa, nSigmaTpcCombinedMaxKa); - selectorPion.setRangePtTOF(ptPidTofMinKa, ptPidTofMaxKa); - selectorPion.setRangeNSigmaTOF(-nSigmaTofMaxKa, nSigmaTofMaxKa); - selectorPion.setRangeNSigmaTOFCondTPC(-nSigmaTofCombinedMaxKa, nSigmaTofCombinedMaxKa); - for (auto& [track1, track2] : combinations(soa::CombinationsStrictlyUpperIndexPolicy(tracks, tracks))) { histos.fill(HIST("DZeroCandTopo/Selections"), 3.f); @@ -1018,52 +916,52 @@ struct qaKFParticle { continue; } - /// Check whether the track was assigned to the true MC PV + // /// Check whether the track was assigned to the true MC PV auto particle1 = track1.mcParticle(); auto particle2 = track2.mcParticle(); - auto collMC = particle1.mcCollision(); - auto mcCollID_recoColl = track1.collision_as().mcCollisionId(); - auto mcCollID_particle = particle1.mcCollisionId(); - bool indexMatchOK = (mcCollID_recoColl == mcCollID_particle); - if (!indexMatchOK) { - histos.fill(HIST("DZeroCandTopo/SelectionsMC"), 4.f); - const auto matchedCollisions = collisions.sliceBy(perMcCollision, collMC.globalIndex()); - int i = 0; - std::array dcaZ{100, 100, 100, 100, 100}; - float min = 100; - for (auto matchedCollision : matchedCollisions) { - dcaZ[i] = abs(matchedCollision.posZ() - collMC.posZ()); - if (i == 0) { - min = dcaZ[i]; - } - if (i > 0) { - if (dcaZ[i] < dcaZ[i - 1]) { - min = dcaZ[i]; - } - } - - i = i + 1; - } - if (min > 10.) { - histos.fill(HIST("DZeroCandTopo/SelectionsMC"), 5.f); - } - int j = 0; - for (auto matchedCollision : matchedCollisions) { - if (i == 1) { - kfpVertex = createKFPVertexFromCollision(matchedCollision); - KFParticle KFPVNew(kfpVertex); - KFPV = KFPVNew; - } - if (i > 1) { - if (abs(matchedCollision.posZ() - collMC.posZ()) == min) { - kfpVertex = createKFPVertexFromCollision(matchedCollision); - KFParticle KFPVNew(kfpVertex); - KFPV = KFPVNew; - } - } - j = j + 1; - } - } + // auto collMC = particle1.mcCollision(); + // auto mcCollID_recoColl = track1.collision_as().mcCollisionId(); + // auto mcCollID_particle = particle1.mcCollisionId(); + // bool indexMatchOK = (mcCollID_recoColl == mcCollID_particle); + // if (!indexMatchOK) { + // histos.fill(HIST("DZeroCandTopo/SelectionsMC"), 4.f); + // const auto matchedCollisions = collisions.sliceBy(perMcCollision, collMC.globalIndex()); + // int i = 0; + // std::array dcaZ{100, 100, 100, 100, 100}; + // float min = 100; + // for (auto matchedCollision : matchedCollisions) { + // dcaZ[i] = abs(matchedCollision.posZ() - collMC.posZ()); + // if (i == 0) { + // min = dcaZ[i]; + // } + // if (i > 0) { + // if (dcaZ[i] < dcaZ[i - 1]) { + // min = dcaZ[i]; + // } + // } + + // i = i + 1; + // } + // if (min > 10.) { + // histos.fill(HIST("DZeroCandTopo/SelectionsMC"), 5.f); + // } + // int j = 0; + // for (auto matchedCollision : matchedCollisions) { + // if (i == 1) { + // kfpVertex = createKFPVertexFromCollision(matchedCollision); + // KFParticle KFPVNew(kfpVertex); + // KFPV = KFPVNew; + // } + // if (i > 1) { + // if (abs(matchedCollision.posZ() - collMC.posZ()) == min) { + // kfpVertex = createKFPVertexFromCollision(matchedCollision); + // KFParticle KFPVNew(kfpVertex); + // KFPV = KFPVNew; + // } + // } + // j = j + 1; + // } + // } /// Remove fake tracks if (!track1.has_mcParticle() || !track2.has_mcParticle()) { @@ -1099,11 +997,6 @@ struct qaKFParticle { continue; } - histos.fill(HIST("Tracks/dcaXYToPV"), track1.dcaXY()); - histos.fill(HIST("Tracks/dcaXYToPV"), track2.dcaXY()); - histos.fill(HIST("Tracks/dcaZToPV"), track1.dcaZ()); - histos.fill(HIST("Tracks/dcaZToPV"), track2.dcaZ()); - KFPTrack kfpTrackPosPi; KFPTrack kfpTrackNegPi; KFPTrack kfpTrackPosKa; @@ -1120,10 +1013,10 @@ struct qaKFParticle { float TOFnSigmaPosKa = 0; float TOFnSigmaNegKa = 0; - int pidKaonTr1 = selectorKaon.getStatusTrackPIDTpcOrTof(track1); - int pidKaonTr2 = selectorKaon.getStatusTrackPIDTpcOrTof(track2); - int pidPionTr1 = selectorPion.getStatusTrackPIDTpcOrTof(track1); - int pidPionTr2 = selectorPion.getStatusTrackPIDTpcOrTof(track2); + bool pidKaonTr1 = SelectPIDCombined(track1, kKPlus); + bool pidKaonTr2 = SelectPIDCombined(track2, kKPlus); + bool pidPionTr1 = SelectPIDCombined(track1, kPiPlus); + bool pidPionTr2 = SelectPIDCombined(track2, kPiPlus); if (track1.isPVContributor() || track2.isPVContributor()) { PVContributor = 1; @@ -1135,7 +1028,7 @@ struct qaKFParticle { int pdgMother = mcParticles.rawIteratorAt(indexRec - mcParticles.offset()).pdgCode(); /// Select D0 and D0bar candidates - if (pidPionTr1 == TrackSelectorPID::Status::PIDAccepted && pidKaonTr2 == TrackSelectorPID::Status::PIDAccepted) { + if (pidPionTr1 && pidKaonTr2) { if (track1.sign() == 1 && track2.sign() == -1) { CandD0 = true; particle1 = track1.mcParticle(); @@ -1164,7 +1057,7 @@ struct qaKFParticle { continue; } } - if (pidKaonTr1 == TrackSelectorPID::Status::PIDAccepted && pidPionTr2 == TrackSelectorPID::Status::PIDAccepted) { + if (pidKaonTr1 && pidPionTr2) { if (track1.sign() == 1 && track2.sign() == -1) { CandD0bar = true; if (pdgMother == 421) { @@ -1234,12 +1127,14 @@ struct qaKFParticle { } fillHistograms(kfpTrackPosPi, kfpTrackNegKa, KFPosPion, KFNegKaon, KFDZero_PV, KFDZero, KFPV, KFDZero_DecayVtx, cosThetaStar); - histos.fill(HIST("TracksKFPi/dcaToPVLargeRangeMCBeforeReassignment"), KFPosPion.GetDistanceFromVertex(KFPVDefault)); - histos.fill(HIST("TracksKFPi/dcaToPVLargeRangeMCAfterReassignment"), KFPosPion.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("TracksKFKa/dcaToPVLargeRangeMCBeforeReassignment"), KFNegKaon.GetDistanceFromVertex(KFPVDefault)); - histos.fill(HIST("TracksKFKa/dcaToPVLargeRangeMCAfterReassignment"), KFNegKaon.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("TracksKFPi/nSigmaTPC"), kfpTrackPosPi.GetPt(), TPCnSigmaPosPi); - histos.fill(HIST("TracksKFKa/nSigmaTPC"), kfpTrackNegKa.GetPt(), TPCnSigmaNegKa); + if (writeHistograms) { + histos.fill(HIST("TracksKFPi/dcaToPVLargeRangeMCBeforeReassignment"), KFPosPion.GetDistanceFromVertex(KFPVDefault)); + histos.fill(HIST("TracksKFPi/dcaToPVLargeRangeMCAfterReassignment"), KFPosPion.GetDistanceFromVertex(KFPV)); + histos.fill(HIST("TracksKFKa/dcaToPVLargeRangeMCBeforeReassignment"), KFNegKaon.GetDistanceFromVertex(KFPVDefault)); + histos.fill(HIST("TracksKFKa/dcaToPVLargeRangeMCAfterReassignment"), KFNegKaon.GetDistanceFromVertex(KFPV)); + histos.fill(HIST("TracksKFPi/nSigmaTPC"), kfpTrackPosPi.GetPt(), TPCnSigmaPosPi); + histos.fill(HIST("TracksKFKa/nSigmaTPC"), kfpTrackNegKa.GetPt(), TPCnSigmaNegKa); + } writeVarTree(kfpTrackPosPi, kfpTrackNegKa, KFPosPion, KFNegKaon, KFDZero_PV, KFDZero, KFPV, KFDZero_DecayVtx, TPCnSigmaPosPi, TOFnSigmaPosPi, TPCnSigmaNegKa, TOFnSigmaNegKa, cosThetaStar, track1, sourceD0); } if (CandD0bar) { @@ -1271,12 +1166,14 @@ struct qaKFParticle { } } fillHistograms(kfpTrackNegPi, kfpTrackPosKa, KFNegPion, KFPosKaon, KFDZeroBar_PV, KFDZeroBar, KFPV, KFDZeroBar_DecayVtx, cosThetaStar); - histos.fill(HIST("TracksKFPi/dcaToPVLargeRangeMCBeforeReassignment"), KFNegPion.GetDistanceFromVertex(KFPVDefault)); - histos.fill(HIST("TracksKFPi/dcaToPVLargeRangeMCAfterReassignment"), KFNegPion.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("TracksKFKa/dcaToPVLargeRangeMCBeforeReassignment"), KFPosKaon.GetDistanceFromVertex(KFPVDefault)); - histos.fill(HIST("TracksKFKa/dcaToPVLargeRangeMCAfterReassignment"), KFPosKaon.GetDistanceFromVertex(KFPV)); - histos.fill(HIST("TracksKFPi/nSigmaTPC"), kfpTrackNegPi.GetPt(), TPCnSigmaNegPi); - histos.fill(HIST("TracksKFKa/nSigmaTPC"), kfpTrackPosKa.GetPt(), TPCnSigmaPosKa); + if (writeHistograms) { + histos.fill(HIST("TracksKFPi/dcaToPVLargeRangeMCBeforeReassignment"), KFNegPion.GetDistanceFromVertex(KFPVDefault)); + histos.fill(HIST("TracksKFPi/dcaToPVLargeRangeMCAfterReassignment"), KFNegPion.GetDistanceFromVertex(KFPV)); + histos.fill(HIST("TracksKFKa/dcaToPVLargeRangeMCBeforeReassignment"), KFPosKaon.GetDistanceFromVertex(KFPVDefault)); + histos.fill(HIST("TracksKFKa/dcaToPVLargeRangeMCAfterReassignment"), KFPosKaon.GetDistanceFromVertex(KFPV)); + histos.fill(HIST("TracksKFPi/nSigmaTPC"), kfpTrackNegPi.GetPt(), TPCnSigmaNegPi); + histos.fill(HIST("TracksKFKa/nSigmaTPC"), kfpTrackPosKa.GetPt(), TPCnSigmaPosKa); + } writeVarTree(kfpTrackPosPi, kfpTrackNegKa, KFPosPion, KFNegKaon, KFDZeroBar_PV, KFDZeroBar, KFPV, KFDZeroBar_DecayVtx, TPCnSigmaNegPi, TOFnSigmaNegPi, TPCnSigmaPosKa, TOFnSigmaPosKa, cosThetaStar, track1, sourceD0Bar); } } diff --git a/Tools/KFparticle/qaKFParticleLc.cxx b/Tools/KFparticle/qaKFParticleLc.cxx new file mode 100644 index 00000000000..334fa708590 --- /dev/null +++ b/Tools/KFparticle/qaKFParticleLc.cxx @@ -0,0 +1,627 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \file qaKFParticleLc.cxx +/// \author Annalena Kalteyer , GSI Darmstadt +/// \brief Task to test the performance of the KFParticle package on the Lc to pKpi decay +/// + +#include "Tools/KFparticle/qaKFParticleLc.h" +#include +#include +#include +#include +#include "TableHelper.h" + +/// includes O2 +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/DCA.h" +#include "ReconstructionDataFormats/Track.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" + +/// includes O2Physics +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/Core/trackUtilities.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/TrackSelectionDefaults.h" +#include "Common/Core/RecoDecay.h" +#include "Tools/KFparticle/KFUtilities.h" + +/// includes KFParticle +#include "KFParticle.h" +#include "KFPTrack.h" +#include "KFPVertex.h" +#include "KFParticleBase.h" +#include "KFVertex.h" + +#ifndef HomogeneousField + +#define HomogeneousField + +#endif + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::dataformats; + +struct qaKFParticleLc { + + /// general steering settings + Configurable isRun3{"isRun3", true, "Is Run3 dataset"}; + Configurable ccdbUrl{"ccdburl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable ccdbPathLut{"ccdbPathLut", "GLO/Param/MatLUT", "Path for LUT parametrization"}; + Configurable ccdbPathGeo{"ccdbPathGeo", "GLO/Config/GeometryAligned", "Path of the geometry file"}; + Configurable ccdbPathGrp{"ccdbPathGrp", "GLO/GRP/GRP", "Path of the grp file (Run 2)"}; + Configurable ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object (Run 3)"}; + Service ccdb; + o2::base::MatLayerCylSet* lut; + o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; + int runNumber; + double magneticField = 0.; + KFParticle KFPion, KFKaon, KFProton, KFLc, KFLc_PV; + + /// option to select good events + Configurable eventSelection{"eventSelection", true, "select good events"}; // currently only sel8 is defined for run3 + /// options to select only specific tracks + Configurable trackSelection{"trackSelection", 1, "Track selection: 0 -> No Cut, 1 -> kGlobalTrack, 2 -> kGlobalTrackWoPtEta, 3 -> kGlobalTrackWoDCA, 4 -> kQualityTracks, 5 -> kInAcceptanceTracks"}; + + /// Particle Identification + // TPC PID + Configurable nSigmaTpcMaxPi{"nSigmaTpcMaxPi", 3., "Nsigma cut on TPC only"}; + Configurable nSigmaTpcMaxKa{"nSigmaTpcMaxKa", 3., "Nsigma cut on TPC only"}; + Configurable nSigmaTpcMaxPr{"nSigmaTpcMaxPr", 3., "Nsigma cut on TPC only"}; + // TOF PID + Configurable ptPidTofMinPi{"ptPidTofMinPi", 0.15, "Lower bound of track pT for TOF PID"}; + Configurable nSigmaTofMaxPi{"nSigmaTofMaxPi", 3., "Nsigma cut on TOF only"}; + Configurable ptPidTofMinKa{"ptPidTofMinKa", 0.15, "Lower bound of track pT for TOF PID"}; + Configurable nSigmaTofMaxKa{"nSigmaTofMaxKa", 3., "Nsigma cut on TOF only"}; + Configurable ptPidTofMinPr{"ptPidTofMinPr", 0.15, "Lower bound of track pT for TOF PID"}; + Configurable nSigmaTofMaxPr{"nSigmaTofMaxPr", 3., "Nsigma cut on TOF only"}; + // TPC & TOF Combined + Configurable nSigmaCombMaxPi{"nSigmaCombMaxPi", 3., "Nsigma cut on TPC & TOF"}; + Configurable nSigmaCombMaxKa{"nSigmaCombMaxKa", 3., "Nsigma cut on TPC & TOF"}; + Configurable nSigmaCombMaxPr{"nSigmaCombMaxPr", 3., "Nsigma cut on TPC & TOF"}; + + /// singe track selections + Configurable d_pTMin{"d_pTMin", 0.3, "minimum momentum for tracks"}; + Configurable d_etaRange{"d_etaRange", 0.8, "eta Range for tracks"}; + Configurable d_dcaXYTrackPV{"d_dcaXYTrackPV", 2., "DCA XY of the daughter tracks to the PV"}; + Configurable d_dcaZTrackPV{"d_dcaZTrackPV", 10., "DCA Z of the daughter tracks to the PV"}; + /// Lc Daughter selections + Configurable d_PtMinPi{"d_PtMinPi", 0., "minimum momentum for Pi from Lc candidates"}; + Configurable d_PtMinKa{"d_PtMinKa", 0., "minimum momentum for Ka from Lc candidates"}; + Configurable d_PtMinPr{"d_PtMinPr", 0., "minimum momentum for Pr from Lc candidates"}; + Configurable d_dist3DSVDau{"d_dist3DSVDau", 1000., "maximum geometrical distance 3D daughter tracks at the SV"}; + /// Lc selection after geometrical fitting + Configurable d_pTMinLc{"d_pTMinLc", 0., "minimum momentum for Lc candidates"}; + Configurable d_pTMaxLc{"d_pTMaxLc", 36., "maximum momentum for Lc candidates"}; + Configurable d_massMin{"d_massMin", 1.65, "minimum mass"}; + Configurable d_massMax{"d_massMax", 2.08, "minimum mass"}; + Configurable d_cosPA{"d_cosPA", -1., "minimum cosine Pointing angle"}; + Configurable d_cosPAXY{"d_cosPAXY", -1., "minimum cosine Pointing angle"}; + Configurable d_distPVSV{"d_distPVSV", -1., "minimum distance between PV and SV"}; + Configurable d_chi2geo{"d_chi2geo", 1000., "maximum chi2 geometrical"}; + /// Lc selection after topological constrain to the PV + Configurable applySelectionWithTopoConst{"applySelectionWithTopoConst", true, "Apply selections constraining the mother to the PV"}; + Configurable d_decayLength{"d_decayLength", 0., "minimum decay length"}; + Configurable d_normdecayLength{"d_normdecayLength", 100., "minimum normalised decay length"}; + Configurable d_chi2topo{"d_chi2topo", 1000., "maximum chi2 topological"}; + /// Option to write D0 variables in a tree + Configurable d_DwnSmplFact{"d_DwnSmplFact", 1., "Downsampling factor for tree"}; + + // Define which track selection should be used: + // 0 -> No track selection is applied + // 1 kGlobalTrack = kQualityTracks | kPrimaryTracks | kInAcceptanceTracks + // kQualityTracks = kTrackType | kTPCNCls | kTPCCrossedRows | kTPCCrossedRowsOverNCls | kTPCChi2NDF | + // kTPCRefit | kITSNCls | kITSChi2NDF | kITSRefit | kITSHits + // kPrimaryTracks = kGoldenChi2 | kDCAxy | kDCAz + // kInAcceptanceTracks = kPtRange | kEtaRange + // 2 kGlobalTrackWoPtEta = kQualityTracks | kPrimaryTracks + // 3 kGlobalTrackWoDCA = kQualityTracks | kInAcceptanceTracks + // 4 kQualityTracks + // 5 kInAcceptanceTracks + Filter trackFilter = (trackSelection.node() == 0) || + ((trackSelection.node() == 1) && requireGlobalTrackInFilter()) || + ((trackSelection.node() == 2) && requireGlobalTrackWoPtEtaInFilter()) || + ((trackSelection.node() == 3) && requireGlobalTrackWoDCAInFilter()) || + ((trackSelection.node() == 4) && requireQualityTracksInFilter()) || + ((trackSelection.node() == 5) && requireTrackCutInFilter(TrackSelectionFlags::kInAcceptanceTracks)); + + using CollisionTableData = soa::Join; + using BigTracks = soa::Join; + using BigTracksExtended = soa::Join; + using BigTracksPID = soa::Join; + using TrackTableData = soa::Join; + + /// Table to be produced + Produces rowKFLc; + + void initMagneticFieldCCDB(o2::aod::BCsWithTimestamps::iterator const& bc, int& mRunNumber, + o2::framework::Service const& ccdb, std::string ccdbPathGrp, o2::base::MatLayerCylSet* lut, + bool isRun3) + { + + if (mRunNumber != bc.runNumber()) { + + LOGF(info, "====== initCCDB function called (isRun3==%d)", isRun3); + if (!isRun3) { // Run 2 GRP object + o2::parameters::GRPObject* grpo = ccdb->getForTimeStamp(ccdbPathGrp, bc.timestamp()); + if (grpo == nullptr) { + LOGF(fatal, "Run 2 GRP object (type o2::parameters::GRPObject) is not available in CCDB for run=%d at timestamp=%llu", bc.runNumber(), bc.timestamp()); + } + o2::base::Propagator::initFieldFromGRP(grpo); + o2::base::Propagator::Instance()->setMatLUT(lut); + LOGF(info, "Setting magnetic field to %d kG for run %d from its GRP CCDB object (type o2::parameters::GRPObject)", grpo->getNominalL3Field(), bc.runNumber()); + } else { // Run 3 GRP object + o2::parameters::GRPMagField* grpo = ccdb->getForTimeStamp(ccdbPathGrp, bc.timestamp()); + if (grpo == nullptr) { + LOGF(fatal, "Run 3 GRP object (type o2::parameters::GRPMagField) is not available in CCDB for run=%d at timestamp=%llu", bc.runNumber(), bc.timestamp()); + } + o2::base::Propagator::initFieldFromGRP(grpo); + o2::base::Propagator::Instance()->setMatLUT(lut); + LOGF(info, "Setting magnetic field to current %f A for run %d from its GRP CCDB object (type o2::parameters::GRPMagField)", grpo->getL3Current(), bc.runNumber()); + } + mRunNumber = bc.runNumber(); + } + } /// end initMagneticFieldCCDB + + void init(InitContext const&) + { + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(ccdbPathLut)); + if (!o2::base::GeometryManager::isGeometryLoaded()) { + ccdb->get(ccdbPathGeo); + } + runNumber = 0; + } /// End init + + /// Function to select collisions + template + bool isSelectedCollision(const T& collision) + { + /// Trigger selection + if (eventSelection && !(isRun3 ? collision.sel8() : collision.sel7())) { // currently only sel8 is defined for run3 + return false; + } + /// Reject collisions with negative covariance matrix elemts on the digonal + if (collision.covXX() < 0. || collision.covYY() < 0. || collision.covZZ() < 0.) { + return false; + } + return true; + } + + /// Function for single track selection + template + bool isSelectedTracks(const T& track1, const T& track2, const T& track3) + { + if ((track1.p() < d_pTMin) || (track2.p() < d_pTMin) || (track3.p() < d_pTMin)) { + return false; + } + /// Eta range + if ((abs(track1.eta()) > d_etaRange) || (abs(track2.eta()) > d_etaRange) || (abs(track3.eta()) > d_etaRange)) { + return false; + } + /// DCA XY of the daughter tracks to the primaty vertex + if ((track1.dcaXY() > d_dcaXYTrackPV) || (track2.dcaXY() > d_dcaXYTrackPV) || (track3.dcaXY() > d_dcaXYTrackPV)) { + return false; + } + /// DCA Z of the daughter tracks to the primaty vertex + if ((track1.dcaZ() > d_dcaZTrackPV) || (track2.dcaZ() > d_dcaZTrackPV) || (track3.dcaZ() > d_dcaZTrackPV)) { + return false; + } + /// reject if the tracks with sum of charge != 1 + if (abs(track1.sign() + track2.sign() + track3.sign()) != 1) { + return false; + } + return true; + } + + template + bool isSelectedLc(const T& trackKaon, const T& trackPion, const T& trackProton) + { + if (!(trackKaon.sign() == -1 && trackPion.sign() == 1 && trackProton.sign() == 1)) { + return false; + } + bool pidKaon = SelectPIDCombined(trackKaon, kKPlus); + bool pidPion = SelectPIDCombined(trackPion, kPiPlus); + bool pidProton = SelectPIDCombined(trackProton, kProton); + if (!(pidKaon && pidPion && pidProton)) { + return false; + } + return true; + } + template + bool isSelectedLcBar(const T& trackKaon, const T& trackPion, const T& trackProton) + { + if (!(trackKaon.sign() == 1 && trackPion.sign() == -1 && trackProton.sign() == -1)) { + return false; + } + bool pidKaon = SelectPIDCombined(trackKaon, kKPlus); + bool pidPion = SelectPIDCombined(trackPion, kPiPlus); + bool pidProton = SelectPIDCombined(trackProton, kProton); + if (!(pidKaon && pidPion && pidProton)) { + return false; + } + return true; + } + + template + bool ReconstructLc(const T& trackKaon, const T& trackPion, const T& trackProton, const T2& KFPV) + { + KFPTrack kfpTrackKa = createKFPTrackFromTrack(trackKaon); + KFPTrack kfpTrackPi = createKFPTrackFromTrack(trackPion); + KFPTrack kfpTrackPr = createKFPTrackFromTrack(trackProton); + KFParticle KFKa(kfpTrackKa, 321); + KFKaon = KFKa; + KFParticle KFPi(kfpTrackPi, 211); + KFPion = KFPi; + KFParticle KFPr(kfpTrackPr, 2212); + KFProton = KFPr; + const KFParticle* LcDaughters[3] = {&KFKaon, &KFPion, &KFProton}; + KFLc.SetConstructMethod(2); + KFLc.Construct(LcDaughters, 3); + if (!isSelectedDaughters(KFPion, KFKaon, KFProton)) { + return false; + } + if (!isSelectedGeo(KFLc, KFPV)) { + return false; + } + KFLc_PV = KFLc; + KFLc_PV.SetProductionVertex(KFPV); + if (applySelectionWithTopoConst) { + if (!isSelectedLcTopo(KFLc_PV, KFPV)) { + return false; + } + } + return true; + } + + template + bool isSelectedDaughters(const T& KFPion, const T& KFKaon, const T& KFProton) + { + /// Pt of daughters + if ((KFPion.GetPt() < d_PtMinPi) || (KFKaon.GetPt() < d_PtMinKa) || (KFProton.GetPt() < d_PtMinPr)) { + return false; + } + /// distance 3D daughter tracks at the secondary vertex + if ((KFPion.GetDistanceFromParticle(KFKaon) > d_dist3DSVDau) || (KFPion.GetDistanceFromParticle(KFProton) > d_dist3DSVDau) || (KFKaon.GetDistanceFromParticle(KFProton) > d_dist3DSVDau)) { + return false; + } + return true; + } + + template + bool isSelectedGeo(const T& KFLc, const T& KFPV) + { + /// Pt selection + if (KFLc.GetPt() < d_pTMinLc || KFLc.GetPt() > d_pTMaxLc) { + return false; + } + /// Mass window selection + if (KFLc.GetMass() < d_massMin || KFLc.GetMass() > d_massMax) { + return false; + } + /// cosine pointing angle selection + if (cpaFromKF(KFLc, KFPV) < d_cosPA) { + return false; + } + /// cosine pointing XY angle selection + if (cpaXYFromKF(KFLc, KFPV) < d_cosPAXY) { + return false; + } + /// Minimum distance between PV and SV + if (KFLc.GetDistanceFromVertex(KFPV) < d_distPVSV) { + return false; + } + /// chi2 geometrical + float chi2geo = KFLc.GetChi2() / KFLc.GetNDF(); + if (chi2geo > d_chi2geo) { + return false; + } + return true; + } + + template + bool isSelectedLcTopo(const T& KFLc_PV, const T& KFPV) + { + /// Pt selection + if (KFLc_PV.GetPt() < d_pTMinLc || KFLc_PV.GetPt() > d_pTMaxLc) { + return false; + } + /// Mass window selection + if (KFLc_PV.GetMass() < d_massMin || KFLc_PV.GetMass() > d_massMax) { + return false; + } + /// decay length selection + if (KFLc_PV.GetDecayLength() < d_decayLength) { + return false; + } + /// decay length error selection + float normdecayLength = KFLc_PV.GetDecayLength() / KFLc_PV.GetErrDecayLength(); + if (normdecayLength < d_normdecayLength) { + return false; + } + /// chi2 topological + float chi2topo = KFLc_PV.GetChi2() / KFLc_PV.GetNDF(); + if (chi2topo > d_chi2topo) { + return false; + } + return true; + } + + template + bool SelectPIDCombined(const T1& track, int particle) + { + switch (particle) { + case kPiPlus: { + if ((track.pt() <= ptPidTofMinPi) && track.hasTPC() && (abs(track.tpcNSigmaPi()) < nSigmaTpcMaxPi)) { + return true; + } else if ((track.pt() > ptPidTofMinPi) && track.hasTPC() && !track.hasTOF() && (abs(track.tpcNSigmaPi()) < nSigmaTpcMaxPi)) { + return true; + } else if ((track.pt() > ptPidTofMinPi) && !track.hasTPC() && track.hasTOF() && (abs(track.tofNSigmaPi()) < nSigmaTofMaxPi)) { + return true; + } else if ((track.pt() > ptPidTofMinPi) && track.hasTPC() && track.hasTOF()) { + float CombinednSigma = 1. / sqrt(2) * sqrt((track.tpcNSigmaPi() * track.tpcNSigmaPi()) + (track.tofNSigmaPi() * track.tofNSigmaPi())); + if (abs(CombinednSigma) < nSigmaCombMaxPi) { + return true; + } else { + return false; + } + } else { + return false; + } + break; + } + case kKPlus: { + if ((track.pt() <= ptPidTofMinKa) && track.hasTPC() && (abs(track.tpcNSigmaKa()) < nSigmaTpcMaxKa)) { + return true; + } else if ((track.pt() > ptPidTofMinKa) && track.hasTPC() && !track.hasTOF() && (abs(track.tpcNSigmaKa()) < nSigmaTpcMaxKa)) { + return true; + } else if ((track.pt() > ptPidTofMinKa) && !track.hasTPC() && track.hasTOF() && (abs(track.tofNSigmaKa()) < nSigmaTofMaxKa)) { + return true; + } else if ((track.pt() > ptPidTofMinKa) && track.hasTPC() && track.hasTOF()) { + float CombinednSigma = 1. / sqrt(2) * sqrt((track.tpcNSigmaKa() * track.tpcNSigmaKa()) + (track.tofNSigmaKa() * track.tofNSigmaKa())); + if (abs(CombinednSigma) < nSigmaCombMaxKa) { + return true; + } else { + return false; + } + } else { + return false; + } + break; + } + case kProton: { + if ((track.pt() <= ptPidTofMinPr) && track.hasTPC() && (abs(track.tpcNSigmaPr()) < nSigmaTpcMaxPr)) { + return true; + } else if ((track.pt() > ptPidTofMinPr) && track.hasTPC() && !track.hasTOF() && (abs(track.tpcNSigmaPr()) < nSigmaTpcMaxPr)) { + return true; + } else if ((track.pt() > ptPidTofMinPr) && !track.hasTPC() && track.hasTOF() && (abs(track.tofNSigmaPr()) < nSigmaTofMaxPr)) { + return true; + } else if ((track.pt() > ptPidTofMinPr) && track.hasTPC() && track.hasTOF()) { + float CombinednSigma = 1. / sqrt(2) * sqrt((track.tpcNSigmaPr() * track.tpcNSigmaPr()) + (track.tofNSigmaPr() * track.tofNSigmaPr())); + if (abs(CombinednSigma) < nSigmaCombMaxPr) { + return true; + } else { + return false; + } + } else { + return false; + } + break; + } + default: { + LOGF(error, "ERROR: Species is not implemented"); + return false; + } + } + } + + template + void writeVarTree(const T2& KFPion, const T2& KFKaon, const T2& KFProton, const T2& KFLc, const T2& KFLc_PV, const T2& KFPV, const T1& trackKa, const T1& trackPi, const T1& trackPr, const int source) + { + + float chi2geo = KFLc.GetChi2() / KFLc.GetNDF(); + float normdecayLength = KFLc_PV.GetDecayLength() / KFLc_PV.GetErrDecayLength(); + float chi2topo = KFLc_PV.GetChi2() / KFLc_PV.GetNDF(); + const double pseudoRndm = trackKa.pt() * 1000. - (int64_t)(trackKa.pt() * 1000); + if (pseudoRndm < d_DwnSmplFact) { + /// Filling the D0 tree + rowKFLc(runNumber, + KFPion.GetPt(), + KFKaon.GetPt(), + KFProton.GetPt(), + KFPion.GetDistanceFromVertexXY(KFPV), + KFKaon.GetDistanceFromVertexXY(KFPV), + KFProton.GetDistanceFromVertexXY(KFPV), + KFPion.GetDistanceFromVertex(KFPV), + KFKaon.GetDistanceFromVertex(KFPV), + KFProton.GetDistanceFromVertex(KFPV), + trackPi.tpcNSigmaPi(), + trackKa.tpcNSigmaKa(), + trackPr.tpcNSigmaPr(), + trackPi.tofNSigmaPi(), + trackKa.tofNSigmaKa(), + trackPr.tofNSigmaPr(), + KFPion.GetDistanceFromVertexXY(KFLc), + KFKaon.GetDistanceFromVertexXY(KFLc), + KFProton.GetDistanceFromVertexXY(KFLc), + KFLc.GetPt(), + KFLc.GetMass(), + cpaFromKF(KFLc, KFPV), + cpaXYFromKF(KFLc, KFPV), + KFLc.GetDistanceFromVertex(KFPV), + KFLc.GetDistanceFromVertexXY(KFPV), + chi2geo, + KFLc_PV.GetPt(), + KFLc_PV.GetMass(), + KFLc_PV.GetDecayLength(), + KFLc_PV.GetDecayLengthXY(), + cpaFromKF(KFLc_PV, KFPV), + KFLc_PV.GetLifeTime(), + normdecayLength, + chi2topo, + source); + } + } + + /// Process function for data + void processData(CollisionTableData::iterator const& collision, soa::Filtered const& tracks, aod::BCsWithTimestamps const&) + { + + auto bc = collision.bc_as(); + if (runNumber != bc.runNumber()) { + initMagneticFieldCCDB(bc, runNumber, ccdb, isRun3 ? ccdbPathGrpMag : ccdbPathGrp, lut, isRun3); + magneticField = o2::base::Propagator::Instance()->getNominalBz(); +/// Set magnetic field for KF vertexing +#ifdef HomogeneousField + KFParticle::SetField(magneticField); +#endif + } + /// Apply event selection + if (!isSelectedCollision(collision)) { + return; + } + /// set KF primary vertex + KFPVertex kfpVertex = createKFPVertexFromCollision(collision); + KFParticle KFPV(kfpVertex); + + for (auto& [track1, track2, track3] : combinations(soa::CombinationsStrictlyUpperIndexPolicy(tracks, tracks, tracks))) { + + /// Apply single track selection + if (!isSelectedTracks(track1, track2, track3)) { + continue; + } + + bool CandLc = false; + bool CandLcbar = false; + int source = 0; + + if (isSelectedLc(track1, track2, track3)) { + CandLc = true; + source = 1; + bool LcReconstructed = ReconstructLc(track1, track2, track3, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track1, track2, track3, source); + } + } + if (isSelectedLc(track1, track3, track2)) { + CandLc = true; + source = 1; + bool LcReconstructed = ReconstructLc(track1, track3, track2, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track1, track3, track2, source); + } + } + if (isSelectedLc(track2, track1, track3)) { + CandLc = true; + source = 1; + bool LcReconstructed = ReconstructLc(track2, track1, track3, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track2, track1, track3, source); + } + } + if (isSelectedLc(track2, track3, track1)) { + CandLc = true; + source = 1; + bool LcReconstructed = ReconstructLc(track2, track3, track1, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track2, track3, track1, source); + } + } + if (isSelectedLc(track3, track1, track2)) { + CandLc = true; + source = 1; + bool LcReconstructed = ReconstructLc(track3, track1, track2, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track3, track1, track2, source); + } + } + if (isSelectedLc(track3, track2, track1)) { + CandLc = true; + source = 1; + bool LcReconstructed = ReconstructLc(track3, track2, track1, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track3, track2, track1, source); + } + } + if (isSelectedLcBar(track1, track2, track3)) { + CandLcbar = true; + source = 2; + bool LcReconstructed = ReconstructLc(track1, track2, track3, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track1, track2, track3, source); + } + } + if (isSelectedLcBar(track1, track3, track2)) { + CandLcbar = true; + source = 2; + bool LcReconstructed = ReconstructLc(track1, track3, track2, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track1, track3, track2, source); + } + } + if (isSelectedLcBar(track2, track1, track3)) { + CandLcbar = true; + source = 2; + bool LcReconstructed = ReconstructLc(track2, track1, track3, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track2, track1, track3, source); + } + } + if (isSelectedLcBar(track2, track3, track1)) { + CandLcbar = true; + source = 2; + bool LcReconstructed = ReconstructLc(track2, track3, track1, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track2, track3, track1, source); + } + } + if (isSelectedLcBar(track3, track1, track2)) { + CandLcbar = true; + source = 2; + bool LcReconstructed = ReconstructLc(track3, track1, track2, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track3, track1, track2, source); + } + } + if (isSelectedLcBar(track3, track2, track1)) { + CandLcbar = true; + source = 2; + bool LcReconstructed = ReconstructLc(track3, track2, track1, KFPV); + if (LcReconstructed) { + writeVarTree(KFPion, KFKaon, KFProton, KFLc, KFLc_PV, KFPV, track3, track2, track1, source); + } + } + if (!CandLc && !CandLcbar) { + continue; + } + } + } + PROCESS_SWITCH(qaKFParticleLc, processData, "process data", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/Tools/KFparticle/qaKFParticleLc.h b/Tools/KFparticle/qaKFParticleLc.h new file mode 100644 index 00000000000..8df0e08aba3 --- /dev/null +++ b/Tools/KFparticle/qaKFParticleLc.h @@ -0,0 +1,111 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file qaKFParticleLc.h +/// \author Annalena Kalteyer + +#ifndef TOOLS_KFPARTICLE_QAKFPARTICLELC_H_ +#define TOOLS_KFPARTICLE_QAKFPARTICLELC_H_ +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Common/Core/trackUtilities.h" +using namespace o2; +using namespace o2::framework; +using namespace o2::track; + +enum Source { + kPrompt = BIT(0), + kNonPrompt = BIT(1), + kReflection = BIT(2) +}; + +namespace o2::aod +{ +namespace kfparticleLc +{ +DECLARE_SOA_COLUMN(RUNNUMBER, Runnumber, int); +DECLARE_SOA_COLUMN(PtPi, ptpi, float); +DECLARE_SOA_COLUMN(PtKa, ptka, float); +DECLARE_SOA_COLUMN(PtPr, ptpr, float); +DECLARE_SOA_COLUMN(DCAXYPiToPV, dcaxypvpi, float); +DECLARE_SOA_COLUMN(DCAXYKaToPV, dcaxypvka, float); +DECLARE_SOA_COLUMN(DCAXYPrToPV, dcaxypvpr, float); +DECLARE_SOA_COLUMN(DCAPiToPV, dcapvpi, float); +DECLARE_SOA_COLUMN(DCAKaToPV, dcapvka, float); +DECLARE_SOA_COLUMN(DCAPrToPV, dcapvpr, float); +DECLARE_SOA_COLUMN(TPCNSigmaPi, tpcnsigmapi, float); +DECLARE_SOA_COLUMN(TPCNSigmaKA, tpcnsigmaka, float); +DECLARE_SOA_COLUMN(TPCNSigmaPr, tpcnsigmapr, float); +DECLARE_SOA_COLUMN(TOFNSigmaPi, tofnsigmapi, float); +DECLARE_SOA_COLUMN(TOFNSigmaKA, tofnsigmaka, float); +DECLARE_SOA_COLUMN(TOFNSigmaPr, tofnsigmapr, float); +DECLARE_SOA_COLUMN(DCAXYPiToSV, distpisv, float); +DECLARE_SOA_COLUMN(DCAXYKaToSV, distkasv, float); +DECLARE_SOA_COLUMN(DCAXYPrToSV, distprsv, float); +DECLARE_SOA_COLUMN(PtGeo, ptgeo, float); +DECLARE_SOA_COLUMN(MassGeo, massgeo, float); +DECLARE_SOA_COLUMN(CosPaGeo, cospageo, float); +DECLARE_SOA_COLUMN(CosPaXYGeo, cospaxygeo, float); +DECLARE_SOA_COLUMN(DCAPVGeo, distpvgeo, float); +DECLARE_SOA_COLUMN(DCAPVXYGeo, distpvxygeo, float); +DECLARE_SOA_COLUMN(Chi2Geo, chi2geo, float); +DECLARE_SOA_COLUMN(PtTopo, pttopo, float); +DECLARE_SOA_COLUMN(MassTopo, masstopo, float); +DECLARE_SOA_COLUMN(DecayLength, decaylength, float); +DECLARE_SOA_COLUMN(DecayLengthXY, decaylengthxy, float); +DECLARE_SOA_COLUMN(CosPaTopo, cospatopo, float); +DECLARE_SOA_COLUMN(Lifetime, lifetime, float); +DECLARE_SOA_COLUMN(NormDecayLength, normdecaylength, float); +DECLARE_SOA_COLUMN(Chi2Topo, chi2topo, float); +DECLARE_SOA_COLUMN(PartSource, source, float); + +} // namespace kfparticleLc + +DECLARE_SOA_TABLE(TreeKFLc, "AOD", "TREEKFLc", + kfparticleLc::RUNNUMBER, + kfparticleLc::PtPi, + kfparticleLc::PtKa, + kfparticleLc::PtPr, + kfparticleLc::DCAXYPiToPV, + kfparticleLc::DCAXYKaToPV, + kfparticleLc::DCAXYPrToPV, + kfparticleLc::DCAPiToPV, + kfparticleLc::DCAKaToPV, + kfparticleLc::DCAPrToPV, + kfparticleLc::TPCNSigmaPi, + kfparticleLc::TPCNSigmaKA, + kfparticleLc::TPCNSigmaPr, + kfparticleLc::TOFNSigmaPi, + kfparticleLc::TOFNSigmaKA, + kfparticleLc::TOFNSigmaPr, + kfparticleLc::DCAXYPiToSV, + kfparticleLc::DCAXYKaToSV, + kfparticleLc::DCAXYPrToSV, + kfparticleLc::PtGeo, + kfparticleLc::MassGeo, + kfparticleLc::CosPaGeo, + kfparticleLc::CosPaXYGeo, + kfparticleLc::DCAPVGeo, + kfparticleLc::DCAPVXYGeo, + kfparticleLc::Chi2Geo, + kfparticleLc::PtTopo, + kfparticleLc::MassTopo, + kfparticleLc::DecayLength, + kfparticleLc::DecayLengthXY, + kfparticleLc::CosPaTopo, + kfparticleLc::Lifetime, + kfparticleLc::NormDecayLength, + kfparticleLc::Chi2Topo, + kfparticleLc::PartSource); + +} // namespace o2::aod + +#endif // TOOLS_KFPARTICLE_QAKFPARTICLELC_H_