[ALICE3] Add proto decayer

ALICE3/TableProducer/OTF/CMakeLists.txt

@@ -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.
 #
@@ -9,17 +9,17 @@
 # granted to it by virtue of its status as an Intergovernmental Organization
 # or submit itself to any jurisdiction.
 
 o2physics_add_dpl_workflow(onthefly-tracker
                     SOURCES onTheFlyTracker.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2::DCAFitter O2Physics::ALICE3Core O2Physics::FastTracker
                     COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(onthefly-tofpid
                     SOURCES onTheFlyTofPid.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core
                     COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(onthefly-richpid
                     SOURCES onTheFlyRichPid.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core
                     COMPONENT_NAME Analysis)
@@ -28,3 +28,9 @@
                     SOURCES onTheFlyTrackerPid.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2Physics::ALICE3Core
                     COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(on-the-fly-decayer
+                    SOURCES onTheFlyDecayer.cxx
+                    PUBLIC_LINK_LIBRARIES O2::Framework
+                    COMPONENT_NAME Analysis)
+

ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx

@@ -0,0 +1,131 @@
+// 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 onTheFlyDecayer.cxx
+///
+/// \brief Task to decay long lived particles and to propagate the information to other tasks
+///
+/// \author Nicolò Jacazio <nicolo.jacazio@cern.ch>, UniUPO
+///
+
+#include "Framework/O2DatabasePDGPlugin.h"
+
+using namespace o2;
+using namespace o2::framework;
+using std::array;
+
+struct OnTheFlyDecayer {
+  Service<o2::framework::O2DatabasePDG> pdgDB;
+
+  void init(o2::framework::InitContext&)
+  {
+  }
+
+  template <typename ParticleType>
+  bool decayParticle(const auto& particle)
+  {
+
+    bool canDecay = false;
+
+    switch (particle.pdgCode()) {
+      case 3312:
+        canDecay = true;
+    }
+    if (!canDecay) {
+      return false;
+    }
+    // Check that it does not have daughters
+    if (particle.hasDaughters()) {
+      LOG(fatal) << "Particle has daughters";
+    }
+
+    const auto& pdgInfo = pdgDB->GetParticle(particle.pdgCode());
+    if (!pdgInfo) {
+      LOG(fatal) << "PDG code " << particle.pdgCode() << " not found in the database";
+    }
+
+    const double u = rand.Uniform(0, 1);
+    double xi_mass = o2::constants::physics::MassXiMinus;
+    double la_mass = o2::constants::physics::MassLambda;
+    double pi_mass = o2::constants::physics::MassPionCharged;
+    double pr_mass = o2::constants::physics::MassProton;
+
+    double mass = 0.;
+    double tau = 0.;
+    // Compute channel
+    switch (particle.pdgCode()) {
+      case 3312:
+        mass = xi_mass;
+        tau = 4.91;
+        break;
+      case 3112:
+        mass = la_mass;
+        tau = 7.89;
+        break;
+    }
+
+    const double gamma = 1 / sqrt(1 + (particle.p() * particle.p()) / (mass * mass));
+    const double ctau = tau * gamma;
+    const double rxyz = (-ctau * log(1.0 - u));
+    // If the particle is charged, then propagate in the mag field
+    o2::math_utils::CircleXYf_t circle;
+    if (pdgInfo->Charge() != 0) {
+      float sna, csa;
+      track.getCircleParams(magneticField, circle, sna, csa);
+    } else { // Neutral particles
+    }
+    const double rxy = rxyz / sqrt(1. + track.getTgl() * track.getTgl());
+    const double theta = rxy / circle.rC;
+    const double newX = ((particle.vx() - circle.xC) * std::cos(theta) - (particle.vy() - circle.yC) * std::sin(theta)) + circle.xC;
+    const double newY = ((particle.vy() - circle.yC) * std::cos(theta) + (particle.vx() - circle.xC) * std::sin(theta)) + circle.yC;
+    const double newPx = particle.px() * std::cos(theta) - particle.py() * std::sin(theta);
+    const double newPy = particle.py() * std::cos(theta) + particle.px() * std::sin(theta);
+    xiDecayVertex.push_back(newX);
+    xiDecayVertex.push_back(newY);
+    xiDecayVertex.push_back(particle.vz() + rxyz * (particle.pz() / particle.p()));
+
+    std::vector<double> xiDaughters = {la_mass, pi_mass};
+    TLorentzVector xi(newPx, newPy, particle.pz(), particle.e());
+    TGenPhaseSpace xiDecay;
+    xiDecay.SetDecay(xi, 2, xiDaughters.data());
+    xiDecay.Generate();
+    decayDaughters.push_back(*xiDecay.GetDecay(1));
+    return true;
+
+    TLorentzVector la = *xiDecay.GetDecay(0);
+
+    double la_gamma = 1 / sqrt(1 + (la.P() * la.P()) / (la_mass * la_mass));
+    double la_ctau = 7.89 * la_gamma;
+    std::vector<double> laDaughters = {pi_mass, pr_mass};
+    double la_rxyz = (-la_ctau * log(1 - u));
+    laDecayVertex.push_back(xiDecayVertex[0] + la_rxyz * (xiDecay.GetDecay(0)->Px() / xiDecay.GetDecay(0)->P()));
+    laDecayVertex.push_back(xiDecayVertex[1] + la_rxyz * (xiDecay.GetDecay(0)->Py() / xiDecay.GetDecay(0)->P()));
+    laDecayVertex.push_back(xiDecayVertex[2] + la_rxyz * (xiDecay.GetDecay(0)->Pz() / xiDecay.GetDecay(0)->P()));
+
+    TGenPhaseSpace laDecay;
+    laDecay.SetDecay(la, 2, laDaughters.data());
+    laDecay.Generate();
+    decayDaughters.push_back(*laDecay.GetDecay(0));
+    decayDaughters.push_back(*laDecay.GetDecay(1));
+  }
+
+  void process(aod::McCollision const& mcCollision,
+               aod::McParticles const& mcParticles)
+  {
+    for (const auto& particle : mcParticles) {
+      decayParticle(particle);
+    }
+  };
+
+  WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+  {
+    return WorkflowSpec{adaptAnalysisTask<OnTheFlyDecayer>(cfgc)};
+  }

ALICE3/Tasks/alice3-dilepton.cxx

@@ -16,16 +16,17 @@
 
 #include "Math/Vector4D.h"
 // O2 includes
-#include "Framework/AnalysisTask.h"
-#include "Framework/ASoAHelpers.h"
-#include "Framework/runDataProcessing.h"
-#include "Framework/HistogramRegistry.h"
-#include "CommonConstants/PhysicsConstants.h"
-#include "Common/DataModel/TrackSelectionTables.h"
-#include "Framework/AnalysisDataModel.h"
-#include "ALICE3/DataModel/OTFTOF.h"
 #include "ALICE3/DataModel/OTFRICH.h"
+#include "ALICE3/DataModel/OTFTOF.h"
 #include "ALICE3/DataModel/tracksAlice3.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+
+#include "CommonConstants/PhysicsConstants.h"
+#include "Framework/ASoAHelpers.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/HistogramRegistry.h"
+#include "Framework/runDataProcessing.h"
 
 using namespace o2;
 using namespace o2::aod;
@@ -128,6 +129,15 @@ struct Alice3Dilepton {
     registry.addClone("Reconstructed/Pair/ULS/", "Reconstructed/Pair/LSpp/");
     registry.addClone("Reconstructed/Pair/ULS/", "Reconstructed/Pair/LSnn/");
 
+    registry.add("ReconstructedFiltered/Pair/ULS/Mass", "Pair Mass", kTH1F, {axisM});
+    registry.add("ReconstructedFiltered/Pair/ULS/Pt", "Pair Pt", kTH1F, {axisPt});
+    registry.add("ReconstructedFiltered/Pair/ULS/Eta", "Pair Eta", kTH1F, {axisEta});
+    registry.add("ReconstructedFiltered/Pair/ULS/Phi", "Pair Phi", kTH1F, {axisPhi});
+    registry.add("ReconstructedFiltered/Pair/ULS/Mass_Pt", "Pair Mass vs. Pt", kTH2F, {axisM, axisPt}, true);
+
+    registry.addClone("ReconstructedFiltered/Pair/ULS/", "ReconstructedFiltered/Pair/LSpp/");
+    registry.addClone("ReconstructedFiltered/Pair/ULS/", "ReconstructedFiltered/Pair/LSnn/");
+
     HistogramConfigSpec hs_rec{HistType::kTHnSparseF, {axisM, axisPt, axisDCAxy}, 3};
     registry.add("Reconstructed/Pair/ULS/hs_rec", "", hs_rec);
     registry.add("Reconstructed/Pair/LSpp/hs_rec", "", hs_rec);
@@ -287,6 +297,69 @@ struct Alice3Dilepton {
     return HFllType::kUndef;
   }
 
+  template <typename T1, typename T2>
+  ROOT::Math::PtEtaPhiMVector buildPairDCA(T1 const& t1, T2 const& t2, float& pair_dca_xy)
+  {
+
+    const float dcaXY_t1 = t1.dcaXY();
+    const float dcaXY_t2 = t2.dcaXY();
+    const float dcaXY_res_t1 = sqrt(t1.cYY());
+    const float dcaXY_res_t2 = sqrt(t2.cYY());
+
+    pair_dca_xy = sqrt((dcaXY_t2 * dcaXY_t2 / dcaXY_res_t2 + dcaXY_t1 * dcaXY_t1 / dcaXY_res_t1) / 2.);
+    ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
+    ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
+    ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
+    return v12;
+  }
+
+  template <PairType pairtype, typename TTracks, typename TMCTracks>
+  void FillPairRecWithPrefilter(TTracks const& tracks1, TTracks const& tracks2, TMCTracks const& mcParticles)
+  {
+    std::vector<uint64_t> prefilteredTracks;
+    if constexpr (pairtype == PairType::kULS) {
+      for (auto& [t1, t2] : combinations(soa::CombinationsFullIndexPolicy(tracks1, tracks2))) {
+        if (!t1.has_mcParticle() || !t2.has_mcParticle()) {
+          continue;
+        }
+
+        float pair_dca_xy = 999.f;
+        ROOT::Math::PtEtaPhiMVector v12 = buildPairDCA(t1, t2, pair_dca_xy);
+        // prefilter for low-mass pairs
+        if (v12.M() > 0.10) {
+          continue;
+        }
+        // prefilter small opening angle pairs
+        if (std::cos(t1.phi() - t2.phi()) < 0.99) {
+          continue;
+        }
+        prefilteredTracks.push_back(t1.globalIndex());
+        prefilteredTracks.push_back(t2.globalIndex());
+
+      } // end of unlike-sign pair loop
+
+      for (auto& [t1, t2] : combinations(soa::CombinationsFullIndexPolicy(tracks1, tracks2))) {
+        // Skipping tracks that are in the prefiltered list
+        if (std::find(prefilteredTracks.begin(), prefilteredTracks.end(), t1.globalIndex()) != prefilteredTracks.end()) {
+          continue;
+        }
+        if (std::find(prefilteredTracks.begin(), prefilteredTracks.end(), t2.globalIndex()) != prefilteredTracks.end()) {
+          continue;
+        }
+
+        float pair_dca_xy = 999.f;
+        ROOT::Math::PtEtaPhiMVector v12 = buildPairDCA(t1, t2, pair_dca_xy);
+
+        registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Mass"), v12.M());
+        registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Pt"), v12.Pt());
+        registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Eta"), v12.Eta());
+        registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Phi"), v12.Phi() < 0.f ? v12.Phi() + TMath::TwoPi() : v12.Phi());
+        registry.fill(HIST("ReconstructedFiltered/Pair/ULS/Mass_Pt"), v12.M(), v12.Pt());
+        registry.fill(HIST("ReconstructedFiltered/Pair/ULS/hs_rec"), v12.M(), v12.Pt(), pair_dca_xy);
+      }
+    }
+  }
+
   template <PairType pairtype, typename TTracks, typename TMCTracks>
   void FillPairRec(TTracks const& tracks1, TTracks const& tracks2, TMCTracks const& mcParticles)
   {
@@ -314,15 +387,8 @@ struct Alice3Dilepton {
         }
         // auto mother = mcparticles.iteratorAt(motherid);
 
-        // float dcaXY_t1 = t1.dcaXY();
-        // float dcaXY_t2 = t2.dcaXY();
-        // float dcaXY_res_t1 = sqrt(t1.cYY());
-        // float dcaXY_res_t2 = sqrt(t2.cYY());
-
-        float pair_dca_xy = sqrt((pow(t2.dcaXY() / sqrt(t2.cYY()), 2) + pow(t1.dcaXY() / sqrt(t1.cYY()), 2)) / 2.);
-        ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
-        ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
-        ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
+        float pair_dca_xy = 999.f;
+        ROOT::Math::PtEtaPhiMVector v12 = buildPairDCA(t1, t2, pair_dca_xy);
 
         registry.fill(HIST("Reconstructed/Pair/ULS/Mass"), v12.M());
         registry.fill(HIST("Reconstructed/Pair/ULS/Pt"), v12.Pt());
@@ -356,15 +422,8 @@ struct Alice3Dilepton {
         }
         // auto mother = mcparticles.iteratorAt(motherid);
 
-        // float dcaXY_t1 = t1.dcaXY();
-        // float dcaXY_t2 = t2.dcaXY();
-        // float dcaXY_res_t1 = sqrt(t1.cYY());
-        // float dcaXY_res_t2 = sqrt(t2.cYY());
-
-        float pair_dca_xy = sqrt((pow(t2.dcaXY() / sqrt(t2.cYY()), 2) + pow(t1.dcaXY() / sqrt(t1.cYY()), 2)) / 2.);
-        ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
-        ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), std::abs(pdg) == 11 ? o2::constants::physics::MassElectron : o2::constants::physics::MassMuon); // reconstructed pt/eta/phi
-        ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
+        float pair_dca_xy = 999.f;
+        ROOT::Math::PtEtaPhiMVector v12 = buildPairDCA(t1, t2, pair_dca_xy);
 
         if constexpr (pairtype == PairType::kLSpp) {
           registry.fill(HIST("Reconstructed/Pair/LSpp/Mass"), v12.M());
@@ -464,6 +523,7 @@ struct Alice3Dilepton {
       } // end of like-sign pair loop
     }
   }
+
   // Functions for pid
   template <typename TTrack>
   bool electronIDTOF(TTrack const& track)
@@ -545,11 +605,10 @@ struct Alice3Dilepton {
   Partition<MyFilteredTracksMC> posTracks = o2::aod::track::signed1Pt > 0.f;
   Partition<MyFilteredTracksMC> negTracks = o2::aod::track::signed1Pt < 0.f;
 
-  void processRec(
-    const o2::aod::Collisions& collisions,
-    MyFilteredTracksMC const& tracks,
-    const o2::aod::McCollisions&,
-    const aod::McParticles& mcParticles)
+  void processRec(const o2::aod::Collisions& collisions,
+                  MyFilteredTracksMC const& tracks,
+                  const o2::aod::McCollisions&,
+                  const aod::McParticles& mcParticles)
   {
     for (const auto& collision : collisions) {
       registry.fill(HIST("Reconstructed/Event/VtxX"), collision.posX());