@@ -349,11 +349,13 @@ class VarManager : public TObject
349349 static void FillEvent (T const & event, float * values = nullptr );
350350 template <uint32_t fillMap, typename T>
351351 static void FillTrack (T const & track, float * values = nullptr );
352- template <int pairType, typename T1 , typename T2 >
352+ template <int pairType, uint32_t fillMap, typename T1 , typename T2 >
353353 static void FillPair (T1 const & t1, T2 const & t2, float * values = nullptr );
354+ template <int pairType, typename T1 , typename T2 >
355+ static void FillPairME (T1 const & t1, T2 const & t2, float * values = nullptr );
354356 template <typename T1 , typename T2 >
355357 static void FillPairMC (T1 const & t1, T2 const & t2, float * values = nullptr , PairCandidateType pairType = kJpsiToEE );
356- template <int pairType, typename C, typename T>
358+ template <int pairType, uint32_t collFillMap, uint32_t fillMap, typename C, typename T>
357359 static void FillPairVertexing (C const & collision, T const & t1, T const & t2, float * values = nullptr );
358360 template <typename T1 , typename T2 >
359361 static void FillDileptonHadron (T1 const & dilepton, T2 const & hadron, float * values = nullptr , float hadronMass = 0 .0f );
@@ -419,7 +421,7 @@ void VarManager::FillEvent(T const& event, float* values)
419421 }
420422
421423 if constexpr ((fillMap & Collision) > 0 ) {
422- // TODO: trigger info from the evenet selection requires a separate flag
424+ // TODO: trigger info from the event selection requires a separate flag
423425 // so that it can be switched off independently of the rest of Collision variables (e.g. if event selection is not available)
424426 if (fgUsedVars[kIsINT7 ]) {
425427 values[kIsINT7 ] = (event.alias ()[kINT7 ] > 0 );
@@ -759,7 +761,7 @@ void VarManager::FillTrack(T const& track, float* values)
759761 FillTrackDerived (values);
760762}
761763
762- template <int pairType, typename T1 , typename T2 >
764+ template <int pairType, uint32_t fillMap, typename T1 , typename T2 >
763765void VarManager::FillPair (T1 const & t1, T2 const & t2, float * values)
764766{
765767 if (!values) {
@@ -791,15 +793,11 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values)
791793 ROOT ::Math::XYZVectorF v2_CM{(boostv12 (v2).Vect ()).Unit ()};
792794 ROOT ::Math::XYZVectorF zaxis{(v12.Vect ()).Unit ()};
793795
794- double cosTheta = 0 ;
795- if (t1.sign () > 0 ) {
796- cosTheta = zaxis.Dot (v1_CM);
797- } else {
798- cosTheta = zaxis.Dot (v2_CM);
796+ if (fgUsedVars[kCosThetaHE ]) {
797+ values[kCosThetaHE ] = (t1.sign () > 0 ? zaxis.Dot (v1_CM) : zaxis.Dot (v2_CM));
799798 }
800- values[kCosThetaHE ] = cosTheta;
801799
802- if constexpr (pairType == kJpsiToEE ) {
800+ if constexpr (( pairType == kJpsiToEE ) && ((fillMap & TrackCov) > 0 || (fillMap & ReducedTrackBarrelCov) > 0 ) ) {
803801
804802 if (fgUsedVars[kQuadDCAabsXY ] || fgUsedVars[kQuadDCAsigXY ]) {
805803 // Quantities based on the barrel tables
@@ -814,6 +812,37 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values)
814812 }
815813}
816814
815+ template <int pairType, typename T1 , typename T2 >
816+ void VarManager::FillPairME (T1 const & t1, T2 const & t2, float * values)
817+ {
818+ //
819+ // Lightweight fill function called from the innermost event mixing loop
820+ //
821+ if (!values) {
822+ values = fgValues;
823+ }
824+
825+ float m1 = fgkElectronMass;
826+ float m2 = fgkElectronMass;
827+ if constexpr (pairType == kJpsiToMuMu ) {
828+ m1 = fgkMuonMass;
829+ m2 = fgkMuonMass;
830+ }
831+
832+ if constexpr (pairType == kElectronMuon ) {
833+ m2 = fgkMuonMass;
834+ }
835+
836+ ROOT ::Math::PtEtaPhiMVector v1 (t1.pt (), t1.eta (), t1.phi (), m1);
837+ ROOT ::Math::PtEtaPhiMVector v2 (t2.pt (), t2.eta (), t2.phi (), m2);
838+ ROOT ::Math::PtEtaPhiMVector v12 = v1 + v2;
839+ values[kMass ] = v12.M ();
840+ values[kPt ] = v12.Pt ();
841+ values[kEta ] = v12.Eta ();
842+ values[kPhi ] = v12.Phi ();
843+ values[kRap ] = -v12.Rapidity ();
844+ }
845+
817846template <typename T1 , typename T2 >
818847void VarManager::FillPairMC (T1 const & t1, T2 const & t2, float * values, PairCandidateType pairType)
819848{
@@ -843,9 +872,14 @@ void VarManager::FillPairMC(T1 const& t1, T2 const& t2, float* values, PairCandi
843872 values[kRap ] = -v12.Rapidity ();
844873}
845874
846- template <int pairType, typename C, typename T>
875+ template <int pairType, uint32_t collFillMap, uint32_t fillMap, typename C, typename T>
847876void VarManager::FillPairVertexing (C const & collision, T const & t1, T const & t2, float * values)
848877{
878+ // check at compile time that the event and cov matrix have the cov matrix
879+ constexpr bool eventHasVtxCov = ((collFillMap & Collision) > 0 || (collFillMap & ReducedEventVtxCov) > 0 );
880+ constexpr bool trackHasCov = ((fillMap & TrackCov) > 0 || (fillMap & ReducedTrackBarrelCov) > 0 );
881+ constexpr bool muonHasCov = ((fillMap & MuonCov) > 0 || (fillMap & ReducedMuonCov) > 0 );
882+
849883 if (!values) {
850884 values = fgValues;
851885 }
@@ -855,7 +889,8 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2,
855889 // TODO: use trackUtilities functions to initialize the various matrices to avoid code duplication
856890 // auto pars1 = getTrackParCov(t1);
857891 // auto pars2 = getTrackParCov(t2);
858- if constexpr (pairType == kJpsiToEE ) {
892+ // We need to hide the cov data members from the cases when no cov table is provided
893+ if constexpr ((pairType == kJpsiToEE ) && trackHasCov) {
859894 std::array<float , 5 > t1pars = {t1.y (), t1.z (), t1.snp (), t1.tgl (), t1.signed1Pt ()};
860895 std::array<float , 15 > t1covs = {t1.cYY (), t1.cZY (), t1.cZZ (), t1.cSnpY (), t1.cSnpZ (),
861896 t1.cSnpSnp (), t1.cTglY (), t1.cTglZ (), t1.cTglSnp (), t1.cTglTgl (),
@@ -867,7 +902,7 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2,
867902 t2.c1PtY (), t2.c1PtZ (), t2.c1PtSnp (), t2.c1PtTgl (), t2.c1Pt21Pt2 ()};
868903 o2::track::TrackParCov pars2{t2.x (), t2.alpha (), t2pars, t2covs};
869904 procCode = fgFitterTwoProng.process (pars1, pars2);
870- } else if constexpr (pairType == kJpsiToMuMu ) {
905+ } else if constexpr (( pairType == kJpsiToMuMu ) && muonHasCov ) {
871906 // Initialize track parameters for forward
872907 double chi21 = t1.chi2 ();
873908 double chi22 = t2.chi2 ();
@@ -884,6 +919,8 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2,
884919 SMatrix55 t2covs (v2.begin (), v2.end ());
885920 o2::track::TrackParCovFwd pars2{t2.z (), t2pars, t2covs, chi22};
886921 procCode = FwdfgFitterTwoProng.process (pars1, pars2);
922+ } else {
923+ return ;
887924 }
888925
889926 values[kVertexingProcCode ] = procCode;
@@ -901,92 +938,91 @@ void VarManager::FillPairVertexing(C const& collision, T const& t1, T const& t2,
901938 values[kVertexingTauz ] = -999 .;
902939 values[kVertexingTauxyErr ] = -999 .;
903940 values[kVertexingTauzErr ] = -999 .;
904-
905941 return ;
906942 }
907943
908944 float m1 = fgkElectronMass;
909945 float m2 = fgkElectronMass;
910-
911946 Vec3D secondaryVertex;
912947 float bz = 0 ;
913948 std::array<float , 3 > pvec0;
914949 std::array<float , 3 > pvec1;
915950
916- double * covMatrixPCA;
917-
918- o2::dataformats::DCA impactParameter0;
919- o2::dataformats::DCA impactParameter1;
920- // get track impact parameters
921- // This modifies track momenta!
922- o2::math_utils::Point3D<float > vtxXYZ (collision.posX (), collision.posY (), collision.posZ ());
923- std::array<float , 6 > vtxCov{collision.covXX (), collision.covXY (), collision.covYY (), collision.covXZ (), collision.covYZ (), collision.covZZ ()};
924- o2::dataformats::VertexBase primaryVertex = {std::move (vtxXYZ), std::move (vtxCov)};
925- // auto primaryVertex = getPrimaryVertex(collision);
926- auto covMatrixPV = primaryVertex.getCov ();
927-
928- if constexpr (pairType == kJpsiToEE ) {
929- secondaryVertex = fgFitterTwoProng.getPCACandidate ();
930- bz = fgFitterTwoProng.getBz ();
931- covMatrixPCA = fgFitterTwoProng.calcPCACovMatrix ().Array ();
932- auto chi2PCA = fgFitterTwoProng.getChi2AtPCACandidate ();
933- auto trackParVar0 = fgFitterTwoProng.getTrack (0 );
934- auto trackParVar1 = fgFitterTwoProng.getTrack (1 );
935- values[kVertexingChi2PCA ] = chi2PCA;
936- trackParVar0.getPxPyPzGlo (pvec0);
937- trackParVar1.getPxPyPzGlo (pvec1);
938- trackParVar0.propagateToDCA (primaryVertex, bz, &impactParameter0);
939- trackParVar1.propagateToDCA (primaryVertex, bz, &impactParameter1);
940- } else if constexpr (pairType == kJpsiToMuMu ) {
941- // Get pca candidate from forward DCA fitter
942- m1 = fgkMuonMass;
943- m2 = fgkMuonMass;
944-
945- secondaryVertex = FwdfgFitterTwoProng.getPCACandidate ();
946- bz = FwdfgFitterTwoProng.getBz ();
947- covMatrixPCA = FwdfgFitterTwoProng.calcPCACovMatrix ().Array ();
948- auto chi2PCA = FwdfgFitterTwoProng.getChi2AtPCACandidate ();
949- auto trackParVar0 = FwdfgFitterTwoProng.getTrack (0 );
950- auto trackParVar1 = FwdfgFitterTwoProng.getTrack (1 );
951- values[kVertexingChi2PCA ] = chi2PCA;
952- pvec0[0 ] = trackParVar0.getPx ();
953- pvec0[1 ] = trackParVar0.getPy ();
954- pvec0[2 ] = trackParVar0.getPz ();
955- pvec1[0 ] = trackParVar1.getPx ();
956- pvec1[1 ] = trackParVar1.getPy ();
957- pvec1[2 ] = trackParVar1.getPz ();
958- trackParVar0.propagateToZlinear (primaryVertex.getZ ());
959- trackParVar1.propagateToZlinear (primaryVertex.getZ ());
960- }
961- // get uncertainty of the decay length
962- // double phi, theta;
963- // getPointDirection(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertex, phi, theta);
964- double phi = std::atan2 (secondaryVertex[1 ] - collision.posY (), secondaryVertex[0 ] - collision.posX ());
965- double theta = std::atan2 (secondaryVertex[2 ] - collision.posZ (),
966- std::sqrt ((secondaryVertex[0 ] - collision.posX ()) * (secondaryVertex[0 ] - collision.posX ()) +
967- (secondaryVertex[1 ] - collision.posY ()) * (secondaryVertex[1 ] - collision.posY ())));
968-
969- values[kVertexingLxyzErr ] = std::sqrt (getRotatedCovMatrixXX (covMatrixPV, phi, theta) + getRotatedCovMatrixXX (covMatrixPCA, phi, theta));
970- values[kVertexingLxyErr ] = std::sqrt (getRotatedCovMatrixXX (covMatrixPV, phi, 0 .) + getRotatedCovMatrixXX (covMatrixPCA, phi, 0 .));
971- values[kVertexingLzErr ] = std::sqrt (getRotatedCovMatrixXX (covMatrixPV, 0 , theta) + getRotatedCovMatrixXX (covMatrixPCA, 0 , theta));
972-
973- values[kVertexingLxy ] = (collision.posX () - secondaryVertex[0 ]) * (collision.posX () - secondaryVertex[0 ]) +
974- (collision.posY () - secondaryVertex[1 ]) * (collision.posY () - secondaryVertex[1 ]);
975- values[kVertexingLz ] = (collision.posZ () - secondaryVertex[2 ]) * (collision.posZ () - secondaryVertex[2 ]);
976- values[kVertexingLxyz ] = values[kVertexingLxy ] + values[kVertexingLz ];
977- values[kVertexingLxy ] = std::sqrt (values[kVertexingLxy ]);
978- values[kVertexingLz ] = std::sqrt (values[kVertexingLz ]);
979- values[kVertexingLxyz ] = std::sqrt (values[kVertexingLxyz ]);
980-
981- ROOT ::Math::PtEtaPhiMVector v1 (t1.pt (), t1.eta (), t1.phi (), m1);
982- ROOT ::Math::PtEtaPhiMVector v2 (t2.pt (), t2.eta (), t2.phi (), m2);
983- ROOT ::Math::PtEtaPhiMVector v12 = v1 + v2;
951+ if constexpr (eventHasVtxCov) {
952+ double * covMatrixPCA;
953+ o2::dataformats::DCA impactParameter0;
954+ o2::dataformats::DCA impactParameter1;
955+ // get track impact parameters
956+ // This modifies track momenta!
957+ o2::math_utils::Point3D<float > vtxXYZ (collision.posX (), collision.posY (), collision.posZ ());
958+ std::array<float , 6 > vtxCov{collision.covXX (), collision.covXY (), collision.covYY (), collision.covXZ (), collision.covYZ (), collision.covZZ ()};
959+ o2::dataformats::VertexBase primaryVertex = {std::move (vtxXYZ), std::move (vtxCov)};
960+ // auto primaryVertex = getPrimaryVertex(collision);
961+ auto covMatrixPV = primaryVertex.getCov ();
962+
963+ if constexpr (pairType == kJpsiToEE && trackHasCov ) {
964+ secondaryVertex = fgFitterTwoProng.getPCACandidate ();
965+ bz = fgFitterTwoProng.getBz ();
966+ covMatrixPCA = fgFitterTwoProng.calcPCACovMatrix ().Array ();
967+ auto chi2PCA = fgFitterTwoProng.getChi2AtPCACandidate ();
968+ auto trackParVar0 = fgFitterTwoProng.getTrack (0 );
969+ auto trackParVar1 = fgFitterTwoProng.getTrack (1 );
970+ values[kVertexingChi2PCA ] = chi2PCA;
971+ trackParVar0.getPxPyPzGlo (pvec0);
972+ trackParVar1.getPxPyPzGlo (pvec1);
973+ trackParVar0.propagateToDCA (primaryVertex, bz, &impactParameter0);
974+ trackParVar1.propagateToDCA (primaryVertex, bz, &impactParameter1);
975+ } else if constexpr (pairType == kJpsiToMuMu && muonHasCov ) {
976+ // Get pca candidate from forward DCA fitter
977+ m1 = fgkMuonMass;
978+ m2 = fgkMuonMass;
979+
980+ secondaryVertex = FwdfgFitterTwoProng.getPCACandidate ();
981+ bz = FwdfgFitterTwoProng.getBz ();
982+ covMatrixPCA = FwdfgFitterTwoProng.calcPCACovMatrix ().Array ();
983+ auto chi2PCA = FwdfgFitterTwoProng.getChi2AtPCACandidate ();
984+ auto trackParVar0 = FwdfgFitterTwoProng.getTrack (0 );
985+ auto trackParVar1 = FwdfgFitterTwoProng.getTrack (1 );
986+ values[kVertexingChi2PCA ] = chi2PCA;
987+ pvec0[0 ] = trackParVar0.getPx ();
988+ pvec0[1 ] = trackParVar0.getPy ();
989+ pvec0[2 ] = trackParVar0.getPz ();
990+ pvec1[0 ] = trackParVar1.getPx ();
991+ pvec1[1 ] = trackParVar1.getPy ();
992+ pvec1[2 ] = trackParVar1.getPz ();
993+ trackParVar0.propagateToZlinear (primaryVertex.getZ ());
994+ trackParVar1.propagateToZlinear (primaryVertex.getZ ());
995+ }
996+ // get uncertainty of the decay length
997+ // double phi, theta;
998+ // getPointDirection(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertex, phi, theta);
999+ double phi = std::atan2 (secondaryVertex[1 ] - collision.posY (), secondaryVertex[0 ] - collision.posX ());
1000+ double theta = std::atan2 (secondaryVertex[2 ] - collision.posZ (),
1001+ std::sqrt ((secondaryVertex[0 ] - collision.posX ()) * (secondaryVertex[0 ] - collision.posX ()) +
1002+ (secondaryVertex[1 ] - collision.posY ()) * (secondaryVertex[1 ] - collision.posY ())));
1003+
1004+ values[kVertexingLxyzErr ] = std::sqrt (getRotatedCovMatrixXX (covMatrixPV, phi, theta) + getRotatedCovMatrixXX (covMatrixPCA, phi, theta));
1005+ values[kVertexingLxyErr ] = std::sqrt (getRotatedCovMatrixXX (covMatrixPV, phi, 0 .) + getRotatedCovMatrixXX (covMatrixPCA, phi, 0 .));
1006+ values[kVertexingLzErr ] = std::sqrt (getRotatedCovMatrixXX (covMatrixPV, 0 , theta) + getRotatedCovMatrixXX (covMatrixPCA, 0 , theta));
1007+
1008+ values[kVertexingLxy ] = (collision.posX () - secondaryVertex[0 ]) * (collision.posX () - secondaryVertex[0 ]) +
1009+ (collision.posY () - secondaryVertex[1 ]) * (collision.posY () - secondaryVertex[1 ]);
1010+ values[kVertexingLz ] = (collision.posZ () - secondaryVertex[2 ]) * (collision.posZ () - secondaryVertex[2 ]);
1011+ values[kVertexingLxyz ] = values[kVertexingLxy ] + values[kVertexingLz ];
1012+ values[kVertexingLxy ] = std::sqrt (values[kVertexingLxy ]);
1013+ values[kVertexingLz ] = std::sqrt (values[kVertexingLz ]);
1014+ values[kVertexingLxyz ] = std::sqrt (values[kVertexingLxyz ]);
1015+
1016+ ROOT ::Math::PtEtaPhiMVector v1 (t1.pt (), t1.eta (), t1.phi (), m1);
1017+ ROOT ::Math::PtEtaPhiMVector v2 (t2.pt (), t2.eta (), t2.phi (), m2);
1018+ ROOT ::Math::PtEtaPhiMVector v12 = v1 + v2;
9841019
985- values[kVertexingTauz ] = (collision.posZ () - secondaryVertex[2 ]) * v12.M () / (TMath::Abs (v12.Pz ()) * o2::constants::physics::LightSpeedCm2NS);
986- values[kVertexingTauxy ] = values[kVertexingLxy ] * v12.M () / (v12.P () * o2::constants::physics::LightSpeedCm2NS);
1020+ values[kVertexingTauz ] = (collision.posZ () - secondaryVertex[2 ]) * v12.M () / (TMath::Abs (v12.Pz ()) * o2::constants::physics::LightSpeedCm2NS);
1021+ values[kVertexingTauxy ] = values[kVertexingLxy ] * v12.M () / (v12.P () * o2::constants::physics::LightSpeedCm2NS);
9871022
988- values[kVertexingTauzErr ] = values[kVertexingLzErr ] * v12.M () / (TMath::Abs (v12.Pz ()) * o2::constants::physics::LightSpeedCm2NS);
989- values[kVertexingTauxyErr ] = values[kVertexingLxyErr ] * v12.M () / (v12.P () * o2::constants::physics::LightSpeedCm2NS);
1023+ values[kVertexingTauzErr ] = values[kVertexingLzErr ] * v12.M () / (TMath::Abs (v12.Pz ()) * o2::constants::physics::LightSpeedCm2NS);
1024+ values[kVertexingTauxyErr ] = values[kVertexingLxyErr ] * v12.M () / (v12.P () * o2::constants::physics::LightSpeedCm2NS);
1025+ }
9901026}
9911027
9921028template <typename T1 , typename T2 >
0 commit comments