From 10be228ccd09017a8f80a661e50ad9765ce723b3 Mon Sep 17 00:00:00 2001 From: sarahherrmann <83958698+sarahherrmann@users.noreply.github.com> Date: Thu, 28 Jul 2022 16:25:06 +0200 Subject: [PATCH 1/4] Commented added template and cut of DCAxy --- PWGMM/Mult/Tasks/vertexing-fwd.cxx | 229 +++++++---------------------- 1 file changed, 53 insertions(+), 176 deletions(-) diff --git a/PWGMM/Mult/Tasks/vertexing-fwd.cxx b/PWGMM/Mult/Tasks/vertexing-fwd.cxx index 94305acc19d..a21905f847c 100644 --- a/PWGMM/Mult/Tasks/vertexing-fwd.cxx +++ b/PWGMM/Mult/Tasks/vertexing-fwd.cxx @@ -8,6 +8,10 @@ // 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. + +//Code written by Robin Caron and Sarah Herrmann +//This code loops over every ambiguous MFT tracks and associates them to a collision that has the smallest DCAxy + #include #include "ReconstructionDataFormats/TrackFwd.h" #include "Math/MatrixFunctions.h" @@ -36,7 +40,7 @@ namespace o2::aod { DECLARE_SOA_TABLE(AmbiguousTracksMFT, "AOD", "AMBIGUOUSTRMFT", //! Table for MFT tracks which are not uniquely associated with a collision o2::soa::Index<>, o2::aod::ambiguous::MFTTrackId, o2::aod::ambiguous::BCIdSlice, o2::soa::Marker<2>); -} +} // namespace o2::aod struct vertexingfwd { @@ -44,10 +48,10 @@ struct vertexingfwd { /// into different std::vector to easily handle them later outside the loops std::vector vecCollForAmb; // vector for collisions associated to an ambiguous track std::vector vecDCACollForAmb; // vector for dca collision associated to an ambiguous track - std::vector vecAmbTrack; // vector for ambiguous track quantities + std::vector vecAmbTrack; // vector for ambiguous track quantities like chi2 and z std::vector vecZposCollForAmb; // vector for z vertex of collisions associated to an ambiguous track - Configurable maxDCAXY{"maxDCAXY", 3.0, "max allowed transverse DCA"}; // To be used when associating ambitrack to collision using best DCA + Configurable maxDCAXY{"maxDCAXY", 6.0, "max allowed transverse DCA"}; // To be used when associating ambitrack to collision using best DCA HistogramRegistry registry{ "registry", @@ -55,17 +59,15 @@ struct vertexingfwd { {"TracksDCAX", "; DCA_{x} (cm); counts", {HistType::kTH1F, {{100, -10, 10}}}}, {"TracksDCAY", "; DCA_{y} (cm); counts", {HistType::kTH1F, {{100, -10, 10}}}}, {"AmbiguousTracksStatus", "; Status; counts", {HistType::kTH1F, {{6, -0.5, 5.5}}}}, - {"NbCollComp", "; NbCollComp", {HistType::kTH1F, {{10, -0.5, 9.5}}}}, + {"NbCollComp", "; NbCollComp", {HistType::kTH1F, {{20, -0.5, 19.5}}}}, {"NumberOfContributors", "; N_{tr} for vertexing; counts", {HistType::kTH1F, {{100, 0, 100}}}}, {"CollisionsMatchIndicesMC", "; Rec. minDCA ambitrack coll.ID; Gen. ambitrack coll.ID", {HistType::kTH2F, {{401, -0.5, 1000.5}, {401, -0.5, 1000.5}}}}, {"TracksDCAXYBest", "; DCA_{xy}^{best} (cm); counts", {HistType::kTH1F, {{100, -1, 10}}}}, - //{"TracksDCAXYOther", "; DCA_{xy}^{best, false} (cm); counts", {HistType::kTH1F, {{100, -1, 10}}}}, {"TracksDCAXYBestFalse", "; DCA_{xy}^{best, false} (cm); counts", {HistType::kTH1F, {{100, -1, 10}}}}, {"EfficiencyZvtx", "; z vertex; Efficiency", {HistType::kTProfile, {{100, -30, 30}}}}, {"DeltaZvtx", "; #delta z (cm); counts", {HistType::kTH1F, {{400, -20, 20}}}}, {"DeltaZvtxBest", "; #delta z = z_{best} - z_{true} (cm); counts", {HistType::kTH1F, {{400, -20, 20}}}}, - //{"DeltaDCAminNcoll", "; N_{coll}; #delta DCA = DCA_{xy}^{min} - DCA_{xy}^{other} (cm); counts", {HistType::kTH2F, {{10, -0.5, 9.5}, {100, -0.5, 9.5}}}}, - {"CorrectMatch", "; Matching value; counts", {HistType::kTH1F, {{4, -0.5, 3.5}}}}}}; + {"CorrectMatch", "; Matching value; counts", {HistType::kTH1F, {{5, -0.5, 4.5}}}}}}; void init(InitContext&) { @@ -73,17 +75,18 @@ struct vertexingfwd { auto* x1 = hstat->GetXaxis(); x1->SetBinLabel(1, "Incorrect match"); x1->SetBinLabel(2, "Correct match"); - x1->SetBinLabel(3, "Total w N_{coll} > 0"); - x1->SetBinLabel(4, "Total w N_{coll} #geq 0"); + x1->SetBinLabel(3, "N_{ambitrack} associable"); + x1->SetBinLabel(4, "N_{ambitrack} w N_{coll} > 0"); + x1->SetBinLabel(5, "N_{ambitrack} total"); auto hstatus = registry.get(HIST("AmbiguousTracksStatus")); auto* x2 = hstatus->GetXaxis(); x2->SetBinLabel(1, "MFT tracks "); x2->SetBinLabel(2, "MFT ambiguous tracks "); - x2->SetBinLabel(3, "is from primary (all) "); - x2->SetBinLabel(4, "is from primary "); - x2->SetBinLabel(5, "not from primary (all) "); - x2->SetBinLabel(6, "not from primary "); + x2->SetBinLabel(3, "All ambiguous primary"); + x2->SetBinLabel(4, "Orphan + primary"); + x2->SetBinLabel(5, "All ambiguous secondary"); + x2->SetBinLabel(6, "Orphan + secondary"); } int getIndexBestCollision(std::vector vecOfDCA, int method = 0) @@ -98,9 +101,9 @@ struct vertexingfwd { return indice; } - void processNew(aod::AmbiguousMFTTracks const& ambitracks, aod::BCs const& bcs, soa::Join const& tracks, FullCollision const& collisions, aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions) + template + void doProcess(T const& ambitracks, aod::BCs const& bcs, MFTTracksLabeled const& tracks, FullCollision const& collisions, aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions) { - int ntracks = tracks.size(); int nambitracks = ambitracks.size(); @@ -117,8 +120,8 @@ struct vertexingfwd { double zVtxMCAmbi = 0; // z vertex associated to the mc collision int mcCollAmbiID = -1; // mc value for the collision containing the ambiguous track - auto track = ambitrack.mfttrack_as(); // Obtain the MFT ambiguous track with the MC labels - // auto extAmbiTrackid = ambitrack.mfttrackId(); // Global index of the MFT ambiguous track + //auto track = ambitrack.mfttrack_as(); // Obtain the MFT ambiguous track with the MC labels + auto track = ambitrack.template mfttrack_as(); //This is the synthax when calling a templated funcction on a template if (!track.has_mcParticle()) { LOGF(warning, "No MC particle for ambiguous track, skip..."); @@ -156,8 +159,10 @@ struct vertexingfwd { if (bcambi.globalBC() != mostProbableBC) { continue; } + //here the bc of the ambitrack is the bc of the collision we are looking at collCounter++; + // We compute the DCAxy of this track wrt the primary vertex of the current collision SMatrix5 tpars(vecAmbTrack[0], vecAmbTrack[1], vecAmbTrack[2], vecAmbTrack[3], vecAmbTrack[4]); // std::vector v1{extAmbiTrack.cXX(), extAmbiTrack.cXY(), extAmbiTrack.cYY(), extAmbiTrack.cPhiX(), extAmbiTrack.cPhiY(), // extAmbiTrack.cPhiPhi(), extAmbiTrack.cTglX(), extAmbiTrack.cTglY(), extAmbiTrack.cTglPhi(), extAmbiTrack.cTglTgl(), @@ -179,6 +184,10 @@ struct vertexingfwd { registry.fill(HIST("TracksDCAY"), dcaY); registry.fill(HIST("NumberOfContributors"), collision.numContrib()); + if (dcaXY > maxDCAXY) { + continue; + } + vecDCACollForAmb.push_back(dcaXY); if (!collision.has_mcCollision()) { @@ -191,188 +200,45 @@ struct vertexingfwd { vecZposCollForAmb.push_back(collision.mcCollision().posZ()); registry.fill(HIST("DeltaZvtx"), collision.mcCollision().posZ() - zVtxMCAmbi); - break; + break; //once we found a bc that corresponds to the bc of the collision we are working on, go out of the bcambi loop + //We then look at the next collision } } registry.fill(HIST("NbCollComp"), collCounter); - registry.fill(HIST("CorrectMatch"), 3.0); // counting for ambiguous track with N collisions >=0 + registry.fill(HIST("CorrectMatch"), 4.0); // counting for ambiguous track with N collisions >=0 - int indexMinDCA = getIndexBestCollision(vecDCACollForAmb, 0); // obtain min value in the stored vector of DCAs - int indexMCcoll = -1; - if (indexMinDCA == -1) { - continue; - } - indexMCcoll = vecCollForAmb[indexMinDCA]; - - registry.fill(HIST("CollisionsMatchIndicesMC"), mcCollAmbiID, indexMCcoll); - - if (collCounter == 1) { - printf("strange ambiguous track of mfttrackId %d\n", ambitrack.mfttrackId()); - if (mcCollAmbiID == indexMCcoll) { - printf("and this is a correct match for the ambiguous track of mfttrackid %d\n", ambitrack.mfttrackId()); - } - } - - if (vecDCACollForAmb.size() == 0) { // do not use the vector with no collisions - if (!particle.isPhysicalPrimary()) { + if (collCounter == 0) { //these are orphan tracks + if (!particle.isPhysicalPrimary()) { //orphan and secondary registry.fill(HIST("AmbiguousTracksStatus"), 5.0); } - if (particle.isPhysicalPrimary()) { + if (particle.isPhysicalPrimary()) { //orphan and primary registry.fill(HIST("AmbiguousTracksStatus"), 3.0); } continue; } - registry.fill(HIST("DeltaZvtxBest"), vecZposCollForAmb[indexMinDCA] - zVtxMCAmbi); - // for (auto& dca : vecDCACollForAmb) { - // if (dca != vecDCACollForAmb[indexMinDCA]) { - // registry.fill(HIST("DeltaDCAminNcoll"), vecCollForAmb.size(), std::abs(vecDCACollForAmb[indexMinDCA] - dca)); - // registry.fill(HIST("TracksDCAXYOther"), dca); - // } - // } - if (mcCollAmbiID == indexMCcoll) { - value = 1.0; - // LOGF(info, " --> Ambitrack correctly associated to collision, dca= %f", vecDCACollForAmb[indexMinDCA]); - } - registry.fill(HIST("TracksDCAXYBest"), vecDCACollForAmb[indexMinDCA]); - registry.fill(HIST("CorrectMatch"), value); - registry.fill(HIST("EfficiencyZvtx"), zVtxMCAmbi, value); - registry.fill(HIST("CorrectMatch"), 2.0); // Counting for amibuous track with N collisions > 0 - - if (value == 0.0) { - registry.fill(HIST("TracksDCAXYBestFalse"), vecDCACollForAmb[indexMinDCA]); // Incorrect association with min DCA - } - - } // ambitracks loop - } - - PROCESS_SWITCH(vertexingfwd, processNew, "Process ambiguous track DCA", true); - - void processOld(aod::AmbiguousTracksMFT const& ambitracks, aod::BCs const& bcs, soa::Join const& tracks, FullCollision const& collisions, aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions) - { - - int ntracks = tracks.size(); - int nambitracks = ambitracks.size(); - - registry.fill(HIST("AmbiguousTracksStatus"), 0.0, ntracks); - registry.fill(HIST("AmbiguousTracksStatus"), 1.0, nambitracks); - - for (auto& ambitrack : ambitracks) { - vecCollForAmb.clear(); - vecDCACollForAmb.clear(); - vecAmbTrack.clear(); - vecZposCollForAmb.clear(); - - double value = 0.0; // matching value for collision association to an ambiguous track - double zVtxMCAmbi = 0; // z vertex associated to the mc collision - int mcCollAmbiID = -1; // mc value for the collision containing the ambiguous track - - auto track = ambitrack.mfttrack_as(); // Obtain the MFT ambiguous track with the MC labels - // auto extAmbiTrackid = ambitrack.mfttrackId(); // Global index of the MFT ambiguous track - - if (!track.has_mcParticle()) { - LOGF(warning, "No MC particle for ambiguous track, skip..."); - continue; - } - auto particle = track.mcParticle(); - mcCollAmbiID = particle.mcCollisionId(); - zVtxMCAmbi = particle.mcCollision().posZ(); - if (particle.isPhysicalPrimary()) { - registry.fill(HIST("AmbiguousTracksStatus"), 2.0); - } else { - registry.fill(HIST("AmbiguousTracksStatus"), 4.0); - } - - // Fill the std::vector for ambiguous tracks with the quantities needed - vecAmbTrack.push_back(track.x()); - vecAmbTrack.push_back(track.y()); - vecAmbTrack.push_back(track.phi()); - vecAmbTrack.push_back(track.tgl()); - vecAmbTrack.push_back(track.signed1Pt()); - vecAmbTrack.push_back(track.z()); - vecAmbTrack.push_back(track.chi2()); - - auto bcambis = ambitrack.bc(); - - int collCounter = 0; - for (auto& collision : collisions) { - uint64_t mostProbableBC = collision.bc().globalBC(); - //uint64_t meanBC = mostProbableBC - std::lround(collision.collisionTime() / (o2::constants::lhc::LHCBunchSpacingNS / 1000)); - //int deltaBC = std::ceil(collision.collisionTimeRes() / (o2::constants::lhc::LHCBunchSpacingNS / 1000) * 4); - - for (auto& bcambi : bcambis) { - if (bcambi.globalBC() != mostProbableBC) { - continue; - } - collCounter++; - - SMatrix5 tpars(vecAmbTrack[0], vecAmbTrack[1], vecAmbTrack[2], vecAmbTrack[3], vecAmbTrack[4]); - // std::vector v1{extAmbiTrack.cXX(), extAmbiTrack.cXY(), extAmbiTrack.cYY(), extAmbiTrack.cPhiX(), extAmbiTrack.cPhiY(), - // extAmbiTrack.cPhiPhi(), extAmbiTrack.cTglX(), extAmbiTrack.cTglY(), extAmbiTrack.cTglPhi(), extAmbiTrack.cTglTgl(), - // extAmbiTrack.c1PtX(), extAmbiTrack.c1PtY(), extAmbiTrack.c1PtPhi(), extAmbiTrack.c1PtTgl(), extAmbiTrack.c1Pt21Pt2()}; - - std::vector v1; // Temporary null vector for the computation of the covariance matrix - SMatrix55 tcovs(v1.begin(), v1.end()); - o2::track::TrackParCovFwd pars1{vecAmbTrack[5], tpars, tcovs, vecAmbTrack[6]}; - - // o2::track::TrackParCovFwd pars1{extAmbiTrack.z(), tpars, tcovs, chi2}; - pars1.propagateToZlinear(collision.posZ()); // track parameters propagation to the position of the z vertex - - const auto dcaX(pars1.getX() - collision.posX()); - const auto dcaY(pars1.getY() - collision.posY()); - auto dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - - registry.fill(HIST("TracksDCAXY"), dcaXY); - registry.fill(HIST("TracksDCAX"), dcaX); - registry.fill(HIST("TracksDCAY"), dcaY); - registry.fill(HIST("NumberOfContributors"), collision.numContrib()); - - vecDCACollForAmb.push_back(dcaXY); - - if (!collision.has_mcCollision()) { - continue; - } - - int mcCollindex = collision.mcCollision().globalIndex(); - vecCollForAmb.push_back(mcCollindex); - - vecZposCollForAmb.push_back(collision.mcCollision().posZ()); - - registry.fill(HIST("DeltaZvtx"), collision.mcCollision().posZ() - zVtxMCAmbi); - break; - } - } - - registry.fill(HIST("NbCollComp"), collCounter); + registry.fill(HIST("CorrectMatch"), 3.0); // counting for ambiguous track with N collisions >0 int indexMinDCA = getIndexBestCollision(vecDCACollForAmb, 0); // obtain min value in the stored vector of DCAs int indexMCcoll = -1; if (indexMinDCA == -1) { - continue; + continue; //if no DCAxy was < maxDCAXY } indexMCcoll = vecCollForAmb[indexMinDCA]; registry.fill(HIST("CollisionsMatchIndicesMC"), mcCollAmbiID, indexMCcoll); - registry.fill(HIST("CorrectMatch"), 3.0); // counting for amibuous track with N collisions >=0 - if (vecDCACollForAmb.size() == 0) { // do not use the vector with no collisions - if (!particle.isPhysicalPrimary()) { - registry.fill(HIST("AmbiguousTracksStatus"), 5.0); - } - if (particle.isPhysicalPrimary()) { - registry.fill(HIST("AmbiguousTracksStatus"), 3.0); + if (collCounter == 1) { //This shouldn't happen after Ruben's fix + printf("strange ambiguous track of mfttrackId %d\n", ambitrack.mfttrackId()); + if (mcCollAmbiID == indexMCcoll) { + printf("and this is a correct match for the ambiguous track of mfttrackid %d\n", ambitrack.mfttrackId()); } - continue; } + registry.fill(HIST("DeltaZvtxBest"), vecZposCollForAmb[indexMinDCA] - zVtxMCAmbi); - // for (auto& dca : vecDCACollForAmb) { - // if (dca != vecDCACollForAmb[indexMinDCA]) { - // registry.fill(HIST("DeltaDCAminNcoll"), vecCollForAmb.size(), std::abs(vecDCACollForAmb[indexMinDCA] - dca)); - // registry.fill(HIST("TracksDCAXYOther"), dca); - // } - // } - if (mcCollAmbiID == indexMCcoll) { + + if (mcCollAmbiID == indexMCcoll) { //correct match value = 1.0; // LOGF(info, " --> Ambitrack correctly associated to collision, dca= %f", vecDCACollForAmb[indexMinDCA]); } @@ -380,13 +246,24 @@ struct vertexingfwd { registry.fill(HIST("CorrectMatch"), value); registry.fill(HIST("EfficiencyZvtx"), zVtxMCAmbi, value); registry.fill(HIST("CorrectMatch"), 2.0); // Counting for amibuous track with N collisions > 0 - if (value == 0.0) { + + if (value == 0.0) { //incorrect match registry.fill(HIST("TracksDCAXYBestFalse"), vecDCACollForAmb[indexMinDCA]); // Incorrect association with min DCA } } // ambitracks loop + } //end of doProcess + + void processNew(aod::AmbiguousMFTTracks const& ambitracks, aod::BCs const& bcs, MFTTracksLabeled const& tracks, FullCollision const& collisions, aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions) + { + doProcess(ambitracks, bcs, tracks, collisions, mcParticles, mcCollisions); } + PROCESS_SWITCH(vertexingfwd, processNew, "Process ambiguous track DCA", true); + void processOld(aod::AmbiguousTracksMFT const& ambitracks, aod::BCs const& bcs, MFTTracksLabeled const& tracks, FullCollision const& collisions, aod::McParticles const& mcParticles, aod::McCollisions const& mcCollisions) + { + doProcess(ambitracks, bcs, tracks, collisions, mcParticles, mcCollisions); + } PROCESS_SWITCH(vertexingfwd, processOld, "Process ambiguous track DCA", false); }; From 6c4ebf2d960fd4adf88d843eaae41a7d4c0fd848 Mon Sep 17 00:00:00 2001 From: sarahherrmann <83958698+sarahherrmann@users.noreply.github.com> Date: Thu, 28 Jul 2022 17:45:28 +0200 Subject: [PATCH 2/4] Update vertexing-fwd.cxx --- PWGMM/Mult/Tasks/vertexing-fwd.cxx | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/PWGMM/Mult/Tasks/vertexing-fwd.cxx b/PWGMM/Mult/Tasks/vertexing-fwd.cxx index a21905f847c..a0c9889fa8b 100644 --- a/PWGMM/Mult/Tasks/vertexing-fwd.cxx +++ b/PWGMM/Mult/Tasks/vertexing-fwd.cxx @@ -9,8 +9,13 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -//Code written by Robin Caron and Sarah Herrmann -//This code loops over every ambiguous MFT tracks and associates them to a collision that has the smallest DCAxy +/// +/// \file vertexing-fwd.cxx +/// \author Robin Caron +/// \author Sarah Herrmann +/// +/// \brief This code loops over every ambiguous MFT tracks and associates +/// them to a collision that has the smallest DCAxy #include #include "ReconstructionDataFormats/TrackFwd.h" From 51058adb71562a9ef60f750a67560b0e2c018535 Mon Sep 17 00:00:00 2001 From: sarahherrmann <83958698+sarahherrmann@users.noreply.github.com> Date: Thu, 28 Jul 2022 17:47:08 +0200 Subject: [PATCH 3/4] Update vertexing-fwd.cxx --- PWGMM/Mult/Tasks/vertexing-fwd.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/PWGMM/Mult/Tasks/vertexing-fwd.cxx b/PWGMM/Mult/Tasks/vertexing-fwd.cxx index a0c9889fa8b..13ac5f69de6 100644 --- a/PWGMM/Mult/Tasks/vertexing-fwd.cxx +++ b/PWGMM/Mult/Tasks/vertexing-fwd.cxx @@ -16,7 +16,6 @@ /// /// \brief This code loops over every ambiguous MFT tracks and associates /// them to a collision that has the smallest DCAxy - #include #include "ReconstructionDataFormats/TrackFwd.h" #include "Math/MatrixFunctions.h" From a1c7b3fca703bcdb9084e99cac7ef850c84f7d45 Mon Sep 17 00:00:00 2001 From: sarahherrmann <83958698+sarahherrmann@users.noreply.github.com> Date: Thu, 28 Jul 2022 17:49:00 +0200 Subject: [PATCH 4/4] clang fix --- PWGMM/Mult/Tasks/vertexing-fwd.cxx | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/PWGMM/Mult/Tasks/vertexing-fwd.cxx b/PWGMM/Mult/Tasks/vertexing-fwd.cxx index 13ac5f69de6..722dae90505 100644 --- a/PWGMM/Mult/Tasks/vertexing-fwd.cxx +++ b/PWGMM/Mult/Tasks/vertexing-fwd.cxx @@ -9,13 +9,13 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// -/// \file vertexing-fwd.cxx -/// \author Robin Caron -/// \author Sarah Herrmann -/// -/// \brief This code loops over every ambiguous MFT tracks and associates -/// them to a collision that has the smallest DCAxy +// \file vertexing-fwd.cxx +// \author Robin Caron +// \author Sarah Herrmann +// +// \brief This code loops over every ambiguous MFT tracks and associates +// them to a collision that has the smallest DCAxy + #include #include "ReconstructionDataFormats/TrackFwd.h" #include "Math/MatrixFunctions.h"