[PWGDQ] Add open charm MC analyses

PWGDQ/Core/CutsLibrary.cxx

@@ -1,4 +1,4 @@
 // 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.
 //
@@ -8,7 +8,7 @@
 // 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.
 //
 // Contact: iarsene@cern.ch, i.c.arsene@fys.uio.no
 //
 #include "PWGDQ/Core/CutsLibrary.h"
@@ -19,7 +19,7 @@
 #include "AnalysisCompositeCut.h"
 #include "VarManager.h"
 
 AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName)
 {
   //
   // define composie cuts, typically combinations of all the ingredients needed for a full cut
@@ -302,7 +302,7 @@
     cut->AddCut(GetAnalysisCut("electronStandardQualityForO2MCdebug"));
     cut->AddCut(GetAnalysisCut("pidCut_lowP_Corr"));
 
     AnalysisCompositeCut* pidCut_highP = new AnalysisCompositeCut("pidCut_highP", "pidCut_highP", kFALSE);
     pidCut_highP->AddCut(GetAnalysisCut("EleInclusion_highP_Corr"));
     pidCut_highP->AddCut(GetAnalysisCut("PionExclusion_highP_Corr"));
     cut->AddCut(pidCut_highP);
@@ -2569,10 +2569,10 @@
     cut->AddCut(GetAnalysisCut("trackQuality_compareDQEMframework"));
     cut->AddCut(GetAnalysisCut("trackDCA1cm"));
     AnalysisCompositeCut* cut_tpc_nSigma = new AnalysisCompositeCut("pid_TPCnSigma", "pid_TPCnSigma", kTRUE);
     cut_tpc_nSigma->AddCut(GetAnalysisCut("lmee_commonDQEM_PID_TPC"));
 
     AnalysisCompositeCut* cut_tof_nSigma = new AnalysisCompositeCut("pid_TOFnSigma", "pid_TOFnSigma", kTRUE);
     cut_tof_nSigma->AddCut(GetAnalysisCut("lmee_commonDQEM_PID_TOF"));
 
     AnalysisCompositeCut* cut_pid_OR = new AnalysisCompositeCut("e_NSigma", "e_NSigma", kFALSE);
     cut_pid_OR->AddCut(cut_tpc_nSigma);
@@ -3172,6 +3172,27 @@
     return cut;
   }
 
+  if (!nameStr.compare("pairCosPointingPos")) {
+    cut->AddCut(GetAnalysisCut("pairCosPointingPos"));
+    return cut;
+  }
+
+  if (!nameStr.compare("pairCosPointingNeg90")) {
+    cut->AddCut(GetAnalysisCut("pairCosPointingNeg90"));
+    return cut;
+  }
+
+  if (!nameStr.compare("pairCosPointingNeg85")) {
+    cut->AddCut(GetAnalysisCut("pairCosPointingNeg85"));
+    return cut;
+  }
+
+  if (!nameStr.compare("pairTauxyzProjectedCosPointing1")) {
+    cut->AddCut(GetAnalysisCut("pairCosPointingNeg"));
+    cut->AddCut(GetAnalysisCut("pairTauxyzProjected1"));
+    return cut;
+  }
+
   // -------------------------------------------------------------------------------------------------
   //
   // Below are a list of single electron single muon and in order or optimize the trigger
@@ -4553,7 +4574,7 @@
     return cut;
   }
 
   if (!nameStr.compare("lmee_commonDQEM_PID_TPC")) { // cut setting to check least common factor between reduced data sets of PWGEM and PWGDQ
     cut->AddCut(VarManager::kTPCnSigmaEl, -2.5, 3., false, VarManager::kPin, 0.0, 1e+10, false);
     cut->AddCut(VarManager::kTPCnSigmaPi, -1e12, 3.5, true, VarManager::kPin, 0.0, 1e+10, false);
     cut->AddCut(VarManager::kTPCnSigmaKa, -3., 3., true, VarManager::kPin, 0.0, 1e+10, false);
@@ -4561,7 +4582,7 @@
     return cut;
   }
 
   if (!nameStr.compare("lmee_commonDQEM_PID_TOF")) { // cut setting to check least common factor between reduced data sets of PWGEM and PWGDQ
     cut->AddCut(VarManager::kTPCnSigmaEl, -2.5, 3., false, VarManager::kPin, 0.0, 1e+10, false);
     cut->AddCut(VarManager::kTPCnSigmaPi, -3., 3.5, true, VarManager::kPin, 0.0, 1e+10, false);
     cut->AddCut(VarManager::kTOFnSigmaEl, -3., 3., false, VarManager::kPin, 0.3, 1e+10, false);
@@ -5962,6 +5983,21 @@
     return cut;
   }
 
+  if (!nameStr.compare("pairCosPointingPos")) {
+    cut->AddCut(VarManager::kCosPointingAngle, 0.9, 1000.);
+    return cut;
+  }
+
+  if (!nameStr.compare("pairCosPointingNeg90")) {
+    cut->AddCut(VarManager::kCosPointingAngle, -1000., -0.9);
+    return cut;
+  }
+
+  if (!nameStr.compare("pairCosPointingNeg85")) {
+    cut->AddCut(VarManager::kCosPointingAngle, -1000., -0.85);
+    return cut;
+  }
+
   // -------------------------------------------------------------------------------------------------
   //
   // Below are a list of single electron single muon and pair selection in order or optimize the trigger

PWGDQ/Core/HistogramsLibrary.cxx

@@ -941,6 +941,7 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
         hm->AddHistogram(histClass, "LxyProj_Pt", "", false, 10, 0.0, 20.0, VarManager::kPt, 1000, -1.0, 1.0, VarManager::kVertexingLxyProjected);
         hm->AddHistogram(histClass, "LxyProj_Mass_Pt", "", false, 50, 2.0, 4.0, VarManager::kMass, 10, 0.0, 20.0, VarManager::kPt, 1000, -1.0, 1.0, VarManager::kVertexingLxyProjected);
         hm->AddHistogram(histClass, "LzProj_Mass_Pt", "", false, 50, 2.0, 4.0, VarManager::kMass, 10, 0.0, 20.0, VarManager::kPt, 1000, -1.0, 1.0, VarManager::kVertexingLzProjected);
+        hm->AddHistogram(histClass, "CosPointingAngle", "", false, 200, -1.0, 1.0, VarManager::kCosPointingAngle);
       }
 
       if (subGroupStr.Contains("kalman-filter")) {
@@ -1378,6 +1379,7 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
         hm->AddHistogram(histClass, "MassD0region_Rapidity", "", false, 140, 1.5, 2.2, VarManager::kMass, 10, -0.8, 0.8, VarManager::kRap);
         hm->AddHistogram(histClass, "MassD0region_eta", "", false, 140, 1.5, 2.2, VarManager::kMass, 40, -2., 2., VarManager::kEta);
         hm->AddHistogram(histClass, "MassD0region_TauxyzProj", "", false, 140, 1.5, 2.2, VarManager::kMass, 200, -0.03, 0.03, VarManager::kVertexingTauxyzProjected);
+        hm->AddHistogram(histClass, "MassD0region_CosPointing", "", false, 140, 1.5, 2.2, VarManager::kMass, 200, -1.0, 1.0, VarManager::kCosPointingAngle);
         hm->AddHistogram(histClass, "MassD0region_VtxNContribReal", "", false, 140, 1.5, 2.2, VarManager::kMass, 50, 0, 50, VarManager::kVtxNcontribReal);
         hm->AddHistogram(histClass, "MassD0region_Rapidity_AveragePt", "", true, 140, 1.5, 2.2, VarManager::kMass, 10, -0.8, 0.8, VarManager::kRap, 150, 0.0, 30.0, VarManager::kPt);
       }

PWGDQ/Core/MCSignalLibrary.cxx

@@ -1185,6 +1185,126 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name)
     return signal;
   }
 
+  //------------------------------------------------------------------------------------
+
+  if (!nameStr.compare("D0")) {
+    MCProng prong(1, {Pdg::kD0}, {true}, {false}, {0}, {0}, {false});
+    signal = new MCSignal(name, "D0", {prong}, {-1});
+    return signal;
+  }
+  if (!nameStr.compare("KPiFromD0")) {
+    MCProng prongKaon(2, {321, Pdg::kD0}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    MCProng prongPion(2, {211, Pdg::kD0}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    signal = new MCSignal(name, "Kaon and pion pair from D0", {prongKaon, prongPion}, {1, 1});
+    return signal;
+  }
+  if (!nameStr.compare("Dcharged")) {
+    MCProng prong(1, {Pdg::kDPlus}, {true}, {false}, {0}, {0}, {false});
+    signal = new MCSignal(name, "D+/-", {prong}, {-1});
+    return signal;
+  }
+  if (!nameStr.compare("Dplus")) {
+    MCProng prong(1, {Pdg::kDPlus}, {false}, {false}, {0}, {0}, {false});
+    signal = new MCSignal(name, "D+", {prong}, {-1});
+    return signal;
+  }
+  if (!nameStr.compare("Dminus")) {
+    MCProng prong(1, {-Pdg::kDPlus}, {false}, {false}, {0}, {0}, {false});
+    signal = new MCSignal(name, "D+", {prong}, {-1});
+    return signal;
+  }
+  if (!nameStr.compare("KPiPiFromDcharged")) {
+    MCProng prongKaon(2, {321, Pdg::kDPlus}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    MCProng prongPion(2, {211, Pdg::kDPlus}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    signal = new MCSignal(name, "Kaon pion pion triplet from D+/-", {prongKaon, prongPion, prongPion}, {1, 1, 1});
+    return signal;
+  }
+  if (!nameStr.compare("KPiPiFromDplus")) {
+    MCProng prongKaon(2, {-321, Pdg::kDPlus}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    MCProng prongPion(2, {211, Pdg::kDPlus}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    signal = new MCSignal(name, "Kaon pion pion triplet from D+", {prongKaon, prongPion, prongPion}, {1, 1, 1});
+    return signal;
+  }
+  if (!nameStr.compare("KPiPiFromDminus")) {
+    MCProng prongKaon(2, {321, -Pdg::kDPlus}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    MCProng prongPion(2, {-211, -Pdg::kDPlus}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    signal = new MCSignal(name, "Kaon pion pion triplet from D-", {prongKaon, prongPion, prongPion}, {1, 1, 1});
+    return signal;
+  }
+  if (!nameStr.compare("Dstar")) {
+    MCProng prong(1, {Pdg::kDStar}, {true}, {false}, {0}, {0}, {false});
+    signal = new MCSignal(name, "D*", {prong}, {-1});
+    return signal;
+  }
+  if (!nameStr.compare("DstarPlus")) {
+    MCProng prong(1, {Pdg::kDStar}, {false}, {false}, {0}, {0}, {false});
+    signal = new MCSignal(name, "D*+", {prong}, {-1});
+    return signal;
+  }
+  if (!nameStr.compare("DstarMinus")) {
+    MCProng prong(1, {-Pdg::kDStar}, {false}, {false}, {0}, {0}, {false});
+    signal = new MCSignal(name, "D*-", {prong}, {-1});
+    return signal;
+  }
+  if (!nameStr.compare("pionFromDstar")) {
+    MCProng prong(2, {211, Pdg::kDStar}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    signal = new MCSignal(name, "Pions from D* decays", {prong}, {1});
+    return signal;
+  }
+  if (!nameStr.compare("D0FromDstar")) {
+    MCProng prong(2, {Pdg::kD0, Pdg::kDStar}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    signal = new MCSignal(name, "D0 from D* decays", {prong}, {1});
+    return signal;
+  }
+  if (!nameStr.compare("KFromD0FromDstar")) {
+    MCProng prong(3, {321, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
+    signal = new MCSignal(name, "Kaons from D0 from D* decays", {prong}, {1});
+    return signal;
+  }
+  if (!nameStr.compare("PiFromD0FromDstar")) {
+    MCProng prong(3, {211, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
+    signal = new MCSignal(name, "Pions from D0 from D* decays", {prong}, {1});
+    return signal;
+  }
+  if (!nameStr.compare("KPiFromD0FromDstar")) {
+    MCProng prongKaon(3, {321, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
+    MCProng prongPion(3, {321, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
+    signal = new MCSignal(name, "Kaon and pion pair from D0 from D* decay", {prongKaon, prongPion}, {1, 1});
+    return signal;
+  }
+  if (!nameStr.compare("KPiPiFromD0FromDstar")) {
+    MCProng prongKaon(3, {321, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
+    MCProng prongPionSecondary(3, {211, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
+    MCProng prongPion(2, {211, Pdg::kDStar}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    signal = new MCSignal(name, "Kaon pion pion triplet from D*", {prongKaon, prongPionSecondary, prongPion}, {2, 2, 1});
+    return signal;
+  }
+  if (!nameStr.compare("KPiPiFromD0FromDstarPlus")) {
+    MCProng prongKaon(3, {-321, Pdg::kD0, Pdg::kDStar}, {false, false, false}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
+    MCProng prongPionSecondary(3, {211, Pdg::kD0, Pdg::kDStar}, {false, false, false}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
+    MCProng prongPion(2, {211, Pdg::kDStar}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    signal = new MCSignal(name, "Kaon pion pion triplet from D*+", {prongKaon, prongPionSecondary, prongPion}, {2, 2, 1});
+    return signal;
+  }
+  if (!nameStr.compare("KPiPiFromD0FromDstarMinus")) {
+    MCProng prongKaon(3, {321, Pdg::kD0, Pdg::kDStar}, {false, false, false}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
+    MCProng prongPionSecondary(3, {-211, Pdg::kD0, Pdg::kDStar}, {false, false, false}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false});
+    MCProng prongPion(2, {-211, Pdg::kDStar}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false});
+    signal = new MCSignal(name, "Kaon pion pion triplet from D*-", {prongKaon, prongPionSecondary, prongPion}, {2, 2, 1});
+    return signal;
+  }
+  if (!nameStr.compare("KFromDplus")) {
+    MCProng prong(2, {321, Pdg::kDPlus}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {502}, {true});
+    prong.SetSourceBit(0, MCProng::kPhysicalPrimary);
+    signal = new MCSignal(name, "Kaons from D+/- decays", {prong}, {-1});
+    return signal;
+  }
+  if (!nameStr.compare("LambdaC")) {
+    MCProng prong(1, {Pdg::kLambdaCPlus}, {true}, {false}, {0}, {0}, {false});
+    signal = new MCSignal(name, "Lambda_c", {prong}, {-1});
+    return signal;
+  }
+
   //--------------------------------------------------------------------------------
 
   if (!nameStr.compare("JpsiFromChic0")) {

PWGDQ/Core/VarManager.cxx

@@ -642,6 +642,8 @@ void VarManager::SetDefaultVarNames()
   fgVariableUnits[kVertexingLzOverErr] = "";
   fgVariableNames[kVertexingLxyzOverErr] = "Pair Lxyz/DLxyz";
   fgVariableUnits[kVertexingLxyzOverErr] = "";
+  fgVariableNames[kCosPointingAngle] = "Cos #theta_{pointing}";
+  fgVariableUnits[kCosPointingAngle] = "";
   fgVariableNames[kKFTrack0DCAxyz] = "Daughter0 DCAxyz";
   fgVariableUnits[kKFTrack0DCAxyz] = "cm";
   fgVariableNames[kKFTrack1DCAxyz] = "Daughter1 DCAxyz";

PWGDQ/Core/VarManager.h

@@ -146,7 +146,6 @@ class VarManager : public TObject
     kDecayToKPi,                // e.g. D0           -> K+ pi- or cc.
     kTripleCandidateToKPiPi,    // e.g. D+ -> K- pi+ pi+
     kTripleCandidateToPKPi,     // e.g. Lambda_c -> p K- pi+
-    kTripleCandidateToKKPi,     // e.g. D_s -> K+ K- pi+
     kNMaxCandidateTypes
   };
 
@@ -2885,21 +2884,6 @@ void VarManager::FillTriple(T1 const& t1, T2 const& t2, T3 const& t3, float* val
     values[kPhi] = v123.Phi();
     values[kRap] = -v123.Rapidity();
   }
-
-  if (pairType == kTripleCandidateToKKPi) {
-    float m1 = o2::constants::physics::MassKaonCharged;
-    float m2 = o2::constants::physics::MassPionCharged;
-
-    ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), m1);
-    ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), m1);
-    ROOT::Math::PtEtaPhiMVector v3(t3.pt(), t3.eta(), t3.phi(), m2);
-    ROOT::Math::PtEtaPhiMVector v123 = v1 + v2 + v3;
-    values[kMass] = v123.M();
-    values[kPt] = v123.Pt();
-    values[kEta] = v123.Eta();
-    values[kPhi] = v123.Phi();
-    values[kRap] = -v123.Rapidity();
-  }
 }
 
 template <uint32_t fillMap, int pairType, typename T1, typename T2>
@@ -3063,6 +3047,21 @@ void VarManager::FillTripleMC(T1 const& t1, T2 const& t2, T3 const& t3, float* v
     values[kPt1] = t1.pt();
     values[kPt2] = t2.pt();
   }
+
+  if (pairType == kTripleCandidateToKPiPi) {
+    float m1 = o2::constants::physics::MassKaonCharged;
+    float m2 = o2::constants::physics::MassPionCharged;
+
+    ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), m1);
+    ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), m2);
+    ROOT::Math::PtEtaPhiMVector v3(t3.pt(), t3.eta(), t3.phi(), m2);
+    ROOT::Math::PtEtaPhiMVector v123 = v1 + v2 + v3;
+    values[kMass] = v123.M();
+    values[kPt] = v123.Pt();
+    values[kEta] = v123.Eta();
+    values[kPhi] = v123.Phi();
+    values[kRap] = -v123.Rapidity();
+  }
 }
 
 template <int pairType, uint32_t collFillMap, uint32_t fillMap, typename C, typename T>
@@ -3574,6 +3573,10 @@ void VarManager::FillTripletVertexing(C const& collision, T const& t1, T const&
       values[kVertexingTauzErr] = values[kVertexingLzErr] * v123.M() / (TMath::Abs(v123.Pz()) * o2::constants::physics::LightSpeedCm2NS);
       values[kVertexingTauxyErr] = values[kVertexingLxyErr] * v123.M() / (v123.Pt() * o2::constants::physics::LightSpeedCm2NS);
 
+      values[kCosPointingAngle] = ((collision.posX() - secondaryVertex[0]) * v123.Px() +
+                                   (collision.posY() - secondaryVertex[1]) * v123.Py() +
+                                   (collision.posZ() - secondaryVertex[2]) * v123.Pz()) /
+                                  (v123.P() * values[VarManager::kVertexingLxyz]);
       // run 2 definitions: Decay length projected onto the momentum vector of the candidate
       values[kVertexingLzProjected] = (secondaryVertex[2] - collision.posZ()) * v123.Pz();
       values[kVertexingLzProjected] = values[kVertexingLzProjected] / TMath::Sqrt(v123.Pz() * v123.Pz());