add decay length related variables for MCtruth pairs

Author: ffionda

PWGDQ/Core/VarManager.h

@@ -1330,6 +1330,8 @@ class VarManager : public TObject
   static void FillTrackCollision(T const& track, C const& collision, float* values = nullptr);
   template <int candidateType, uint32_t fillMap, typename T1, typename T2, typename C>
   static void FillTrackCollisionMC(T1 const& track, T2 const& MotherTrack, C const& collision, float* values = nullptr);
+  template <int candidateType, typename T1>
+  static void FillTrackCollisionMC(T1 const& track, const std::array<double, 3>& collPos, float massHyp = -1., float* values = nullptr);
   template <uint32_t fillMap, typename T, typename C, typename M, typename P>
   static void FillTrackCollisionMatCorr(T const& track, C const& collision, M const& materialCorr, P const& propagator, float* values = nullptr);
   template <typename U, typename T>
@@ -3301,6 +3303,36 @@ void VarManager::FillTrackCollisionMC(T1 const& track, T2 const& MotherTrack, C
   values[kMCVertexingTauxyzProjected] = values[kMCVertexingLxyzProjected] * m / (MotherTrack.p());
 }
 
+template <int candidateType, typename T1>
+void VarManager::FillTrackCollisionMC(T1 const& track, const std::array<double, 3>& collPos, float massHyp, float* values){
+
+ if (!values) {
+    values = fgValues;
+  }
+
+  float m = o2::constants::physics::MassJPsi;
+  if(massHyp) m = massHyp;
+
+  // displaced vertex is computed with decay product (track) and vertex collPos
+  values[kMCVertexingLxy] = (collPos[0]-track.vx()) * (collPos[0]-track.vx()) +
+                            (collPos[1]-track.vy()) * (collPos[1]-track.vy());
+  values[kMCVertexingLz] = (collPos[2]-track.vz()) * (collPos[2]-track.vz());
+  values[kMCVertexingLxyz] = values[kMCVertexingLxy] + values[kMCVertexingLz];
+  values[kMCVertexingLxy] = std::sqrt(values[kMCVertexingLxy]);
+  values[kMCVertexingLz] = std::sqrt(values[kMCVertexingLz]);
+  values[kMCVertexingLxyz] = std::sqrt(values[kMCVertexingLxyz]);
+  values[kMCVertexingTauz] = (collPos[2]-track.vz()) * m / (TMath::Abs(track.pz()) * o2::constants::physics::LightSpeedCm2NS);
+  values[kMCVertexingTauxy] = values[kMCVertexingLxy] * m / (track.pt() * o2::constants::physics::LightSpeedCm2NS);
+
+  values[kMCVertexingLzProjected] = ((track.vz()-collPos[2]) * track.pz()) / TMath::Abs(track.pz());
+  values[kMCVertexingLxyProjected] = (((track.vx()-collPos[0]) * track.px()) + ((track.vy()-collPos[1]) * track.py())) / TMath::Abs(track.pt());
+  values[kMCVertexingLxyzProjected] = (((track.vx()-collPos[0]) * track.px()) + ((track.vy()-collPos[1]) * track.py()) + ((track.vz()-collPos[2]) * track.pz())) / track.p();
+  values[kMCVertexingTauxyProjected] = values[kMCVertexingLxyProjected] * m / (track.pt());
+  values[kMCVertexingTauzProjected] = values[kMCVertexingLzProjected] * m / TMath::Abs(track.pz());
+  values[kMCVertexingTauxyzProjected] = values[kMCVertexingLxyzProjected] * m / (track.p());
+}
+
+
 template <int pairType, typename T, typename T1>
 void VarManager::FillEnergyCorrelatorsMC(T const& track, T1 const& t1, float* values, float Translow, float Transhigh)
 {

PWGDQ/Tasks/dqEfficiency_withAssoc.cxx

@@ -40,6 +40,7 @@
 #include <Framework/Configurable.h>
 #include <Framework/InitContext.h>
 #include <Framework/runDataProcessing.h>
+#include <Framework/O2DatabasePDGPlugin.h>
 
 #include <THashList.h>
 #include <TList.h>
@@ -1397,6 +1398,9 @@ struct AnalysisSameEventPairing {
 
   Service<o2::ccdb::BasicCCDBManager> fCCDB;
 
+  // PDG database
+  Service<o2::framework::O2DatabasePDG> pdgDB;
+
   // Filter filterEventSelected = aod::dqanalysisflags::isEventSelected & uint32_t(1);
   Filter eventFilter = aod::dqanalysisflags::isEventSelected > static_cast<uint32_t>(0);
 
@@ -2275,6 +2279,16 @@ struct AnalysisSameEventPairing {
                 // cout << "      Signal matched!" << endl;
                 mcDecision |= (static_cast<uint32_t>(1) << isig);
                 VarManager::FillPairMC<TPairType>(t1_raw, t2_raw);
+		// check if t1_raw and t2_raw have same mother to compute decay length related variables
+		if(t1_raw.has_mothers() && t2_raw.has_mothers()){ 
+			auto motherMCParticle_t1 = t1_raw.template mothers_first_as<ReducedMCTracks>();
+                	auto motherMCParticle_t2 = t2_raw.template mothers_first_as<ReducedMCTracks>();
+                	if(motherMCParticle_t1 == motherMCParticle_t2){
+                	auto mcEvent = mcEvents.rawIteratorAt(motherMCParticle_t1.reducedMCeventId());
+                	std::array<double, 3> collVtxPos = {mcEvent.mcPosX(), mcEvent.mcPosY(), mcEvent.mcPosZ()};
+                	VarManager::FillTrackCollisionMC<TPairType>(motherMCParticle_t1,collVtxPos,pdgDB->Mass(motherMCParticle_t1.pdgCode()));
+               		}
+		}
                 // cout << "      Filled VarManager for the pair." << endl;
                 fHistMan->FillHistClass(Form("MCTruthGenPairSel_%s", sig->GetName()), VarManager::fgValues);
                 fHistMan->FillHistClass(Form("MCTruthGenPseudoPolPairSel_%s", sig->GetName()), VarManager::fgValues);