diff --git a/ALICE3/TableProducer/OTF/CMakeLists.txt b/ALICE3/TableProducer/OTF/CMakeLists.txt index 69a4d2ec722..260dc179a40 100644 --- a/ALICE3/TableProducer/OTF/CMakeLists.txt +++ b/ALICE3/TableProducer/OTF/CMakeLists.txt @@ -28,3 +28,9 @@ o2physics_add_dpl_workflow(on-the-fly-tracker-pid 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) + diff --git a/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx b/ALICE3/TableProducer/OTF/onTheFlyDecayer.cxx new file mode 100644 index 00000000000..41a7881e368 --- /dev/null +++ b/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 , UniUPO +/// + +#include "Framework/O2DatabasePDGPlugin.h" + +using namespace o2; +using namespace o2::framework; +using std::array; + +struct OnTheFlyDecayer { + Service pdgDB; + + void init(o2::framework::InitContext&) + { + } + + template + 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 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 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(cfgc)}; + } diff --git a/ALICE3/Tasks/alice3-dilepton.cxx b/ALICE3/Tasks/alice3-dilepton.cxx index e3632fc5a02..455de21fc34 100644 --- a/ALICE3/Tasks/alice3-dilepton.cxx +++ b/ALICE3/Tasks/alice3-dilepton.cxx @@ -14,18 +14,21 @@ /// \author s.scheid@cern.ch, daiki.sekihata@cern.ch /// -#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 +#include +#include +#include +#include +#include + +#include + +#include using namespace o2; using namespace o2::aod; @@ -128,6 +131,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 +299,69 @@ struct Alice3Dilepton { return HFllType::kUndef; } + template + 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 + void FillPairRecWithPrefilter(TTracks const& tracks1, TTracks const& tracks2, TMCTracks const& /*mcParticles*/) + { + std::vector 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 void FillPairRec(TTracks const& tracks1, TTracks const& tracks2, TMCTracks const& mcParticles) { @@ -314,15 +389,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 +424,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 +525,7 @@ struct Alice3Dilepton { } // end of like-sign pair loop } } + // Functions for pid template bool electronIDTOF(TTrack const& track) @@ -545,11 +607,10 @@ struct Alice3Dilepton { Partition posTracks = o2::aod::track::signed1Pt > 0.f; Partition 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());