From 79c5685e599eef31164a7c0030db5e84ae5bff31 Mon Sep 17 00:00:00 2001 From: nzardosh Date: Wed, 16 Sep 2020 18:21:56 +0200 Subject: [PATCH 1/2] restructuring PWGJE restructuring jets --- Analysis/DataModel/include/Analysis/Jet.h | 56 ++++- .../DataModel/include/Analysis/JetFinder.h | 234 ++++++++++++++++++ Analysis/Tasks/CMakeLists.txt | 9 +- Analysis/Tasks/PWGJE/CMakeLists.txt | 43 ++++ Analysis/Tasks/PWGJE/jetfinder.cxx | 102 ++++++++ .../Tasks/PWGJE/jetfinderhadronrecoil.cxx | 126 ++++++++++ Analysis/Tasks/PWGJE/jetfinderhf.cxx | 112 +++++++++ Analysis/Tasks/PWGJE/jetskimming.cxx | 69 ++++++ Analysis/Tasks/PWGJE/jetsubstructure.cxx | 124 ++++++++++ Analysis/Tasks/jetfinder.cxx | 125 ---------- 10 files changed, 860 insertions(+), 140 deletions(-) create mode 100644 Analysis/DataModel/include/Analysis/JetFinder.h create mode 100644 Analysis/Tasks/PWGJE/CMakeLists.txt create mode 100644 Analysis/Tasks/PWGJE/jetfinder.cxx create mode 100644 Analysis/Tasks/PWGJE/jetfinderhadronrecoil.cxx create mode 100644 Analysis/Tasks/PWGJE/jetfinderhf.cxx create mode 100644 Analysis/Tasks/PWGJE/jetskimming.cxx create mode 100644 Analysis/Tasks/PWGJE/jetsubstructure.cxx delete mode 100644 Analysis/Tasks/jetfinder.cxx diff --git a/Analysis/DataModel/include/Analysis/Jet.h b/Analysis/DataModel/include/Analysis/Jet.h index 0d08276852a7d..6f36198fb58af 100644 --- a/Analysis/DataModel/include/Analysis/Jet.h +++ b/Analysis/DataModel/include/Analysis/Jet.h @@ -10,7 +10,7 @@ // table definitions for jets // -// Author: Jochen Klein +// Author: Jochen Klein, Nima Zardoshti #pragma once @@ -21,17 +21,31 @@ namespace o2::aod namespace jet { DECLARE_SOA_INDEX_COLUMN(Collision, collision); +DECLARE_SOA_COLUMN(Pt, pt, float); DECLARE_SOA_COLUMN(Eta, eta, float); DECLARE_SOA_COLUMN(Phi, phi, float); -DECLARE_SOA_COLUMN(Pt, pt, float); -DECLARE_SOA_COLUMN(Area, area, float); DECLARE_SOA_COLUMN(Energy, energy, float); DECLARE_SOA_COLUMN(Mass, mass, float); +DECLARE_SOA_COLUMN(Area, area, float); +DECLARE_SOA_DYNAMIC_COLUMN(Px, px, [](float pt, float phi) { return pt * TMath::Cos(phi); }); +DECLARE_SOA_DYNAMIC_COLUMN(Py, py, [](float pt, float phi) { return pt * TMath::Sin(phi); }); +DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, [](float pt, float eta) { return pt * TMath::SinH(eta); }); +DECLARE_SOA_DYNAMIC_COLUMN(P, p, [](float pt, float eta) { return pt * TMath::CosH(eta); }); //absolute p } // namespace jet DECLARE_SOA_TABLE(Jets, "AOD", "JET", - o2::soa::Index<>, jet::CollisionId, - jet::Eta, jet::Phi, jet::Pt, jet::Area, jet::Energy, jet::Mass); + o2::soa::Index<>, + jet::CollisionId, + jet::Pt, + jet::Eta, + jet::Phi, + jet::Energy, + jet::Mass, + jet::Area, + jet::Px, + jet::Py, + jet::Pz, + jet::P); using Jet = Jets::iterator; @@ -43,8 +57,36 @@ DECLARE_SOA_INDEX_COLUMN(Jet, jet); DECLARE_SOA_INDEX_COLUMN(Track, track); } // namespace constituents -DECLARE_SOA_TABLE(JetConstituents, "AOD", "JETCONSTITUENTS", - constituents::JetId, constituents::TrackId); +DECLARE_SOA_TABLE(JetConstituents, "AOD", "CONSTITUENTS", + constituents::JetId, + constituents::TrackId); using JetConstituent = JetConstituents::iterator; + +namespace constituentssub +{ +DECLARE_SOA_INDEX_COLUMN(Jet, jet); +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(Eta, eta, float); +DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(Energy, energy, float); +DECLARE_SOA_COLUMN(Mass, mass, float); +DECLARE_SOA_DYNAMIC_COLUMN(Px, px, [](float pt, float phi) { return pt * TMath::Cos(phi); }); +DECLARE_SOA_DYNAMIC_COLUMN(Py, py, [](float pt, float phi) { return pt * TMath::Sin(phi); }); +DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, [](float pt, float eta) { return pt * TMath::SinH(eta); }); +DECLARE_SOA_DYNAMIC_COLUMN(P, p, [](float pt, float eta) { return pt * TMath::CosH(eta); }); //absolute p +} //namespace constituentssub +DECLARE_SOA_TABLE(JetConstituentsSub, "AOD", "CONSTITUENTSSUB", + constituentssub::JetId, + constituentssub::Pt, + constituentssub::Eta, + constituentssub::Phi, + constituentssub::Energy, + constituentssub::Mass, + constituentssub::Px, + constituentssub::Py, + constituentssub::Pz, + constituentssub::P); +using JetConstituentSub = JetConstituentsSub::iterator; + } // namespace o2::aod diff --git a/Analysis/DataModel/include/Analysis/JetFinder.h b/Analysis/DataModel/include/Analysis/JetFinder.h new file mode 100644 index 0000000000000..7ff2aa5cc233b --- /dev/null +++ b/Analysis/DataModel/include/Analysis/JetFinder.h @@ -0,0 +1,234 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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. + +// jet finder task +// +// Authors: Nima Zardoshti, Jochen Klein + +#include +#include +#include + +#include "fastjet/PseudoJet.hh" +#include "fastjet/ClusterSequenceArea.hh" +#include "fastjet/AreaDefinition.hh" +#include "fastjet/JetDefinition.hh" +#include "fastjet/tools/JetMedianBackgroundEstimator.hh" +#include "fastjet/tools/Subtractor.hh" +#include "fastjet/contrib/ConstituentSubtractor.hh" + +#include + +float mPion = TDatabasePDG::Instance()->GetParticle(211)->Mass(); //can be removed when pion mass becomes default for unidentified tracks + +class JetFinder +{ + + public: + enum class BkgSubMode { none, + rhoAreaSub, + constSub }; + BkgSubMode bkgSubMode; + + void setBkgSubMode(BkgSubMode bSM) { bkgSubMode = bSM; } + + /// Performs jet finding + /// \note the input particle and jet lists are passed by reference + /// \param inputParticles vector of input particles/tracks + /// \param jets veector of jets to be filled + /// \return ClusterSequenceArea object needed to access constituents + // fastjet::ClusterSequenceArea findJets(std::vector &inputParticles, std::vector &jets); + + float phiMin; + float phiMax; + float etaMin; + float etaMax; + + float jetR; + float jetPtMin; + float jetPtMax; + float jetPhiMin; + float jetPhiMax; + float jetEtaMin; + float jetEtaMax; + + float ghostEtaMin; + float ghostEtaMax; + float ghostArea; + int ghostRepeatN; + double ghostktMean; + float gridScatter; + float ktScatter; + + float jetBkgR; + float bkgPhiMin; + float bkgPhiMax; + float bkgEtaMin; + float bkgEtaMax; + + float constSubAlpha; + float constSubRMax; + + bool isReclustering; + + fastjet::JetAlgorithm algorithm; + fastjet::RecombinationScheme recombScheme; + fastjet::Strategy strategy; + fastjet::AreaType areaType; + fastjet::GhostedAreaSpec ghostAreaSpec; + fastjet::JetDefinition jetDef; + fastjet::AreaDefinition areaDef; + fastjet::Selector selJets; + fastjet::Selector selGhosts; + + fastjet::JetAlgorithm algorithmBkg; + fastjet::RecombinationScheme recombSchemeBkg; + fastjet::Strategy strategyBkg; + fastjet::AreaType areaTypeBkg; + fastjet::JetDefinition jetDefBkg; + fastjet::AreaDefinition areaDefBkg; + fastjet::Selector selRho; + + /// Default constructor + JetFinder(float eta_Min = -0.9, float eta_Max = 0.9, float phi_Min = 0.0, float phi_Max = 2 * TMath::Pi()) : phiMin(phi_Min), + phiMax(phi_Max), + etaMin(eta_Min), + etaMax(eta_Max), + jetR(0.4), + jetPtMin(0.0), + jetPtMax(1000.0), + jetPhiMin(phi_Min), + jetPhiMax(phi_Max), + jetEtaMin(eta_Min), + jetEtaMax(eta_Max), + ghostEtaMin(eta_Min), + ghostEtaMax(eta_Max), + ghostArea(0.005), + ghostRepeatN(1), + ghostktMean(1e-100), //is float precise enough? + gridScatter(1.0), + ktScatter(0.1), + jetBkgR(0.2), + bkgPhiMin(phi_Min), + bkgPhiMax(phi_Max), + bkgEtaMin(eta_Min), + bkgEtaMax(eta_Max), + algorithm(fastjet::antikt_algorithm), + recombScheme(fastjet::E_scheme), + strategy(fastjet::Best), + areaType(fastjet::active_area), + algorithmBkg(fastjet::JetAlgorithm(fastjet::kt_algorithm)), + recombSchemeBkg(fastjet::RecombinationScheme(fastjet::E_scheme)), + strategyBkg(fastjet::Best), + areaTypeBkg(fastjet::active_area), + bkgSubMode(BkgSubMode::none), + constSubAlpha(1.0), + constSubRMax(0.6), + isReclustering(false) + + { + + //default constructor + } + + /// Default destructor + ~JetFinder() = default; + + /// Sets the jet finding parameters + void setParams() + { + + if (!isReclustering) { + jetEtaMin = etaMin + jetR; //in aliphysics this was (-etaMax + 0.95*jetR) + jetEtaMax = etaMax - jetR; + } + + //selGhosts =fastjet::SelectorRapRange(ghostEtaMin,ghostEtaMax) && fastjet::SelectorPhiRange(phiMin,phiMax); + //ghostAreaSpec=fastjet::GhostedAreaSpec(selGhosts,ghostRepeatN,ghostArea,gridScatter,ktScatter,ghostktMean); + ghostAreaSpec = fastjet::GhostedAreaSpec(ghostEtaMax, ghostRepeatN, ghostArea, gridScatter, ktScatter, ghostktMean); //the first argument is rapidity not pseudorapidity, to be checked + jetDef = fastjet::JetDefinition(algorithm, jetR, recombScheme, strategy); + areaDef = fastjet::AreaDefinition(areaType, ghostAreaSpec); + selJets = fastjet::SelectorPtRange(jetPtMin, jetPtMax) && fastjet::SelectorEtaRange(jetEtaMin, jetEtaMax) && fastjet::SelectorPhiRange(jetPhiMin, jetPhiMax); + jetDefBkg = fastjet::JetDefinition(algorithmBkg, jetBkgR, recombSchemeBkg, strategyBkg); + areaDefBkg = fastjet::AreaDefinition(areaTypeBkg, ghostAreaSpec); + selRho = fastjet::SelectorRapRange(bkgEtaMin, bkgEtaMax) && fastjet::SelectorPhiRange(bkgPhiMin, bkgPhiMax); //&& !fastjet::SelectorNHardest(2) //here we have to put rap range, to be checked! + } + + /// Sets the background subtraction estimater pointer + void setBkgE() + { + if (bkgSubMode == BkgSubMode::rhoAreaSub || bkgSubMode == BkgSubMode::constSub) { + bkgE = decltype(bkgE)(new fastjet::JetMedianBackgroundEstimator(selRho, jetDefBkg, areaDefBkg)); + } else { + if (bkgSubMode != BkgSubMode::none) + LOGF(ERROR, "requested subtraction mode not implemented!"); + } + } + + /// Sets the background subtraction pointer + void setSub() + { + //if rho < 1e-6 it is set to 1e-6 in AliPhysics + if (bkgSubMode == BkgSubMode::rhoAreaSub) { + sub = decltype(sub){new fastjet::Subtractor{bkgE.get()}}; + } else if (bkgSubMode == BkgSubMode::constSub) { //event or jetwise + constituentSub = decltype(constituentSub){new fastjet::contrib::ConstituentSubtractor{bkgE.get()}}; + constituentSub->set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR); + constituentSub->set_max_distance(constSubRMax); + constituentSub->set_alpha(constSubAlpha); + constituentSub->set_ghost_area(ghostArea); + constituentSub->set_max_eta(bkgEtaMax); + constituentSub->set_background_estimator(bkgE.get()); //what about rho_m + } else { + if (bkgSubMode != BkgSubMode::none) + LOGF(ERROR, "requested subtraction mode not implemented!"); + } + } + + /// Performs jet finding + /// \note the input particle and jet lists are passed by reference + /// \param inputParticles vector of input particles/tracks + /// \param jets veector of jets to be filled + /// \return ClusterSequenceArea object needed to access constituents + fastjet::ClusterSequenceArea findJets(std::vector& inputParticles, std::vector& jets) //ideally find a way of passing the cluster sequence as a reeference + { + setParams(); + setBkgE(); + jets.clear(); + + if (bkgE) { + bkgE->set_particles(inputParticles); + setSub(); + } + if (constituentSub) { + inputParticles = constituentSub->subtract_event(inputParticles); + } + fastjet::ClusterSequenceArea clusterSeq(inputParticles, jetDef, areaDef); + jets = sub ? (*sub)(clusterSeq.inclusive_jets()) : clusterSeq.inclusive_jets(); + jets = selJets(jets); + return clusterSeq; + } + + private: + //void setParams(); + //void setBkgSub(); + std::unique_ptr bkgE; + std::unique_ptr sub; + std::unique_ptr constituentSub; +}; + +//does this belong here? +template +void fillConstituents(const T& constituent, std::vector& constituents) +{ + + auto energy = std::sqrt(constituent.p() * constituent.p() + mPion * mPion); + constituents.emplace_back(constituent.px(), constituent.py(), constituent.pz(), energy); +} diff --git a/Analysis/Tasks/CMakeLists.txt b/Analysis/Tasks/CMakeLists.txt index db35ed2080e91..9be71512c5e7e 100644 --- a/Analysis/Tasks/CMakeLists.txt +++ b/Analysis/Tasks/CMakeLists.txt @@ -11,6 +11,7 @@ add_subdirectory(PWGCF) add_subdirectory(PWGDQ) add_subdirectory(PWGHF) +add_subdirectory(PWGJE) add_subdirectory(PWGLF) add_subdirectory(PWGUD) @@ -44,14 +45,6 @@ o2_add_dpl_workflow(validation PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase COMPONENT_NAME Analysis) -if(FastJet_FOUND) -o2_add_dpl_workflow(jetfinder - SOURCES jetfinder.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisCore O2::AnalysisDataModel - FastJet::FastJet - COMPONENT_NAME Analysis) -endif() - o2_add_dpl_workflow(event-selection SOURCES eventSelection.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsBase O2::CCDB diff --git a/Analysis/Tasks/PWGJE/CMakeLists.txt b/Analysis/Tasks/PWGJE/CMakeLists.txt new file mode 100644 index 0000000000000..889fef51ba988 --- /dev/null +++ b/Analysis/Tasks/PWGJE/CMakeLists.txt @@ -0,0 +1,43 @@ +# Copyright CERN and copyright holders of ALICE O2. This software is distributed +# under the terms of the GNU General Public License v3 (GPL Version 3), copied +# verbatim in the file "COPYING". +# +# See http://alice-o2.web.cern.ch/license for full licensing information. +# +# 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. + + +if(FastJet_FOUND) +o2_add_dpl_workflow(jet-finder + SOURCES jetfinder.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisCore O2::AnalysisDataModel + FastJet::FastJet FastJet::ConstituentSubtractor + COMPONENT_NAME Analysis) + +o2_add_dpl_workflow(jet-finder-hf + SOURCES jetfinderhf.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisCore O2::AnalysisDataModel + FastJet::FastJet FastJet::ConstituentSubtractor + COMPONENT_NAME Analysis) + +o2_add_dpl_workflow(jet-finder-hadron-recoil + SOURCES jetfinderhadronrecoil.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisCore O2::AnalysisDataModel + FastJet::FastJet FastJet::ConstituentSubtractor + COMPONENT_NAME Analysis) + +o2_add_dpl_workflow(jet-substructure + SOURCES jetsubstructure.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisCore O2::AnalysisDataModel + FastJet::FastJet FastJet::ConstituentSubtractor + COMPONENT_NAME Analysis) + +o2_add_dpl_workflow(jet-skimmer + SOURCES jetskimming.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisCore O2::AnalysisDataModel + FastJet::FastJet FastJet::ConstituentSubtractor + COMPONENT_NAME Analysis) +endif() + diff --git a/Analysis/Tasks/PWGJE/jetfinder.cxx b/Analysis/Tasks/PWGJE/jetfinder.cxx new file mode 100644 index 0000000000000..187b130a961e6 --- /dev/null +++ b/Analysis/Tasks/PWGJE/jetfinder.cxx @@ -0,0 +1,102 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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. + +// jet finder task +// +// Author: Jochen Klein, Nima Zardoshti + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoA.h" + +#include "fastjet/PseudoJet.hh" +#include "fastjet/ClusterSequenceArea.hh" + +#include "Analysis/Jet.h" +#include "Analysis/JetFinder.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct JetFinderTask { + Produces jetsTable; + Produces constituentsTable; + Produces constituentsSubTable; + OutputObj hJetPt{"h_jet_pt"}; + OutputObj hJetPhi{"h_jet_phi"}; + OutputObj hJetEta{"h_jet_eta"}; + OutputObj hJetN{"h_jet_n"}; + + Configurable b_DoRhoAreaSub{"b_DoRhoAreaSub", false, "do rho area subtraction"}; + Configurable b_DoConstSub{"b_DoConstSub", false, "do constituent subtraction"}; + + Filter trackCuts = aod::track::pt >= 0.15f && aod::track::eta >= -0.9f && aod::track::eta <= 0.9f; + + std::vector jets; + std::vector inputParticles; + JetFinder jetFinder; + + void init(InitContext const&) + { + hJetPt.setObject(new TH1F("h_jet_pt", "jet p_{T};p_{T} (GeV/#it{c})", + 100, 0., 100.)); + hJetPhi.setObject(new TH1F("h_jet_phi", "jet #phi;#phi", + 80, -1., 7.)); + hJetEta.setObject(new TH1F("h_jet_eta", "jet #eta;#eta", + 70, -0.7, 0.7)); + hJetN.setObject(new TH1F("h_jet_n", "jet n;n constituents", + 30, 0., 30.)); + if (b_DoRhoAreaSub) + jetFinder.setBkgSubMode(JetFinder::BkgSubMode::rhoAreaSub); + if (b_DoConstSub) + jetFinder.setBkgSubMode(JetFinder::BkgSubMode::constSub); + } + + void process(aod::Collision const& collision, + soa::Filtered const& tracks) + { + + jets.clear(); + inputParticles.clear(); + + for (auto& track : tracks) { + /* auto energy = std::sqrt(track.p() * track.p() + mPion*mPion); + inputParticles.emplace_back(track.px(), track.py(), track.pz(), energy); + inputParticles.back().set_user_index(track.globalIndex());*/ + fillConstituents(track, inputParticles); + inputParticles.back().set_user_index(track.globalIndex()); + } + + fastjet::ClusterSequenceArea clusterSeq(jetFinder.findJets(inputParticles, jets)); + + for (const auto& jet : jets) { + jetsTable(collision, jet.pt(), jet.eta(), jet.phi(), + jet.E(), jet.m(), jet.area()); + hJetPt->Fill(jet.pt()); + hJetPhi->Fill(jet.phi()); + hJetEta->Fill(jet.eta()); + hJetN->Fill(jet.constituents().size()); + for (const auto& constituent : jet.constituents()) { //event or jetwise + if (b_DoConstSub) + constituentsSubTable(jetsTable.lastIndex(), constituent.pt(), constituent.eta(), constituent.phi(), + constituent.E(), constituent.m()); + constituentsTable(jetsTable.lastIndex(), constituent.user_index()); + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const&) +{ + return WorkflowSpec{ + adaptAnalysisTask("jet-finder")}; +} diff --git a/Analysis/Tasks/PWGJE/jetfinderhadronrecoil.cxx b/Analysis/Tasks/PWGJE/jetfinderhadronrecoil.cxx new file mode 100644 index 0000000000000..d85576b3be959 --- /dev/null +++ b/Analysis/Tasks/PWGJE/jetfinderhadronrecoil.cxx @@ -0,0 +1,126 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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. + +// jet finder task +// +// Author: Jochen Klein + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoA.h" + +#include "fastjet/PseudoJet.hh" +#include "fastjet/ClusterSequenceArea.hh" + +#include "Analysis/Jet.h" +#include "Analysis/JetFinder.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct JetFinderHadronRecoilTask { + OutputObj hJetPt{"h_jet_pt"}; + OutputObj hHadronPt{"h_hadron_pt"}; + OutputObj hJetHadronDeltaPhi{"h_jet_hadron_deltaphi"}; + + Configurable f_trackTTMin{"f_trackTTMin", 8.0, "TT hadron min pT"}; + Configurable f_trackTTMax{"f_trackTTMax", 9.0, "TT hadron max pT"}; + Configurable f_recoilWindow{"f_recoilWindow", 0.6, "jet finding reecoil window phi distance from TT hadrom"}; + + Filter trackCuts = aod::track::pt >= 0.15f && aod::track::eta > -0.9f && aod::track::eta < 0.9f; + int collisionSplit = 0; //can we partition the collisions? + //can we also directly filter the collision based on the max track::pt? + + std::vector trackTTPt; + std::vector jets; + std::vector inputParticles; + JetFinder jetFinder; + + template + T relativePhi(T phi1, T phi2) + { + if (phi1 < -TMath::Pi()) { + phi1 += (2 * TMath::Pi()); + } else if (phi1 > TMath::Pi()) { + phi1 -= (2 * TMath::Pi()); + } + if (phi2 < -TMath::Pi()) { + phi2 += (2 * TMath::Pi()); + } else if (phi2 > TMath::Pi()) { + phi2 -= (2 * TMath::Pi()); + } + T deltaPhi = phi2 - phi1; + if (deltaPhi < -TMath::Pi()) { + deltaPhi += (2 * TMath::Pi()); + } else if (deltaPhi > TMath::Pi()) { + deltaPhi -= (2 * TMath::Pi()); + } + return deltaPhi; + } + + void init(InitContext const&) + { + hJetPt.setObject(new TH1F("h_jet_pt", "jet p_{T};jet p_{T} (GeV/#it{c})", + 100, 0., 100.)); + hHadronPt.setObject(new TH1F("h_hadron_pt", "hadron p_{T};hadron p_{T} (GeV/#it{c})", + 120, 0., 60.)); + hJetHadronDeltaPhi.setObject(new TH1F("h_jet_hadron_deltaphi", "jet #eta;#eta", + 40, 0.0, 4.)); + } + + void process(aod::Collision const& collision, + soa::Filtered const& tracks) + { + + jets.clear(); + inputParticles.clear(); + auto trackTTPhi = 0.0; + auto trackTTPt = 0.0; + auto comparisonPt = 0.0; + bool isTT = false; + for (auto& track : tracks) { + if (track.pt() >= f_trackTTMin && track.pt() < f_trackTTMax) { //can this also go into a partition? + isTT = true; + if (track.pt() >= comparisonPt) { //currently take highest pT but later to randomise + comparisonPt = track.pt(); + trackTTPt = track.pt(); + trackTTPhi = track.phi(); + } + } + auto energy = std::sqrt(track.p() * track.p() + mPion * mPion); + inputParticles.emplace_back(track.px(), track.py(), track.pz(), energy); + inputParticles.back().set_user_index(track.globalIndex()); + } + if (!isTT) + return; + hHadronPt->Fill(trackTTPt); + + // you can set phi selector here for jets + fastjet::ClusterSequenceArea clusterSeq(jetFinder.findJets(inputParticles, jets)); + + for (const auto& jet : jets) { + auto deltaPhi = TMath::Abs(relativePhi(jet.phi(), trackTTPhi)); + if (deltaPhi >= (TMath::Pi() - f_recoilWindow)) { + hJetPt->Fill(jet.pt()); + } + if (deltaPhi >= TMath::Pi() / 2.0 && deltaPhi <= TMath::Pi()) { + hJetHadronDeltaPhi->Fill(deltaPhi); + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const&) +{ + return WorkflowSpec{ + adaptAnalysisTask("jet-finder-hadron-recoil")}; +} diff --git a/Analysis/Tasks/PWGJE/jetfinderhf.cxx b/Analysis/Tasks/PWGJE/jetfinderhf.cxx new file mode 100644 index 0000000000000..7c25428caa3ec --- /dev/null +++ b/Analysis/Tasks/PWGJE/jetfinderhf.cxx @@ -0,0 +1,112 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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. + +// jet finder task +// +// Author: Nima Zardoshti + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoA.h" + +#include "Analysis/SecondaryVertexHF.h" +#include "Analysis/CandidateSelectionTables.h" + +#include "fastjet/PseudoJet.hh" +#include "fastjet/ClusterSequenceArea.hh" + +#include "Analysis/Jet.h" +#include "Analysis/JetFinder.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct JetFinderHFTask { + Produces jetsTable; + Produces constituents; + OutputObj hJetPt{"h_jet_pt"}; + OutputObj hD0Pt{"h_D0_pt"}; + + std::vector jets; + std::vector inputParticles; + JetFinder jetFinder; + + void init(InitContext const&) + { + hJetPt.setObject(new TH1F("h_jet_pt", "jet p_{T};p_{T} (GeV/#it{c})", + 100, 0., 100.)); + hD0Pt.setObject(new TH1F("h_D0_pt", "jet p_{T,D};p_{T,D} (GeV/#it{c})", + 60, 0., 60.)); + } + + Configurable d_selectionFlagD0{"d_selectionFlagD0", 1, "Selection Flag for D0"}; + Configurable d_selectionFlagD0bar{"d_selectionFlagD0bar", 1, "Selection Flag for D0bar"}; + + //need enum as configurable + enum pdgCode { pdgD0 = 421 }; + + Filter trackCuts = (aod::track::pt > 0.15f && aod::track::eta > -0.9f && aod::track::eta < 0.9f); + Filter seltrack = (aod::hfselcandidate::isSelD0 >= d_selectionFlagD0 || aod::hfselcandidate::isSelD0bar >= d_selectionFlagD0bar); + + void process(aod::Collision const& collision, + soa::Filtered const& tracks, soa::Filtered> const& candidates) + { + std::cout << "Per Event" << std::endl; + // TODO: retrieve pion mass from somewhere + bool isHFJet; + + //this loop should be made more efficient + for (auto& candidate : candidates) { + jets.clear(); + inputParticles.clear(); + for (auto& track : tracks) { + auto energy = std::sqrt(track.p() * track.p() + mPion * mPion); + if (candidate.index0().globalIndex() == track.globalIndex() || candidate.index1().globalIndex() == track.globalIndex()) { //is it global index? + continue; + } + inputParticles.emplace_back(track.px(), track.py(), track.pz(), energy); + inputParticles.back().set_user_index(track.globalIndex()); + } + inputParticles.emplace_back(candidate.px(), candidate.py(), candidate.pz(), candidate.e(RecoDecay::getMassPDG(pdgD0))); + inputParticles.back().set_user_index(1); + + fastjet::ClusterSequenceArea clusterSeq(jetFinder.findJets(inputParticles, jets)); + + for (const auto& jet : jets) { + isHFJet = false; + for (const auto& constituent : jet.constituents()) { + if (constituent.user_index() == 1 && (candidate.isSelD0() == 1 || candidate.isSelD0bar() == 1)) { + isHFJet = true; + break; + } + } + if (isHFJet) { + jetsTable(collision, jet.eta(), jet.phi(), jet.pt(), + jet.area(), jet.E(), jet.m()); + for (const auto& constituent : jet.constituents()) { + constituents(jetsTable.lastIndex(), constituent.user_index()); + } + hJetPt->Fill(jet.pt()); + std::cout << "Filling" << std::endl; + hD0Pt->Fill(candidate.pt()); + break; + } + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const&) +{ + return WorkflowSpec{ + adaptAnalysisTask("jet-finder-hf")}; +} diff --git a/Analysis/Tasks/PWGJE/jetskimming.cxx b/Analysis/Tasks/PWGJE/jetskimming.cxx new file mode 100644 index 0000000000000..490ce1e284e3b --- /dev/null +++ b/Analysis/Tasks/PWGJE/jetskimming.cxx @@ -0,0 +1,69 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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. + +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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. + +// jet skimming task +// +// Author: Nima Zardoshti + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoA.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +namespace o2::aod +{ +namespace jetskim +{ +DECLARE_SOA_INDEX_COLUMN(Collision, collision); +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(Eta, eta, float); +DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(Energy, energy, float); +} // namespace jetskim +DECLARE_SOA_TABLE(JetSkim, "AOD", "JETSKIM1", + jetskim::CollisionId, + jetskim::Pt, jetskim::Eta, jetskim::Phi, jetskim::Energy); +} // namespace o2::aod + +struct JetSkimmingTask1 { + Produces skim; + + Filter trackCuts = aod::track::pt > 0.15f; + float mPionSquared = 0.139 * 0.139; + + void process(aod::Collision const& collision, + soa::Filtered const& tracks) + { + for (auto& track : tracks) { + float energy = std::sqrt(track.p() * track.p() + mPionSquared); + skim(collision, track.pt(), track.eta(), track.phi(), energy); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const&) +{ + return WorkflowSpec{ + adaptAnalysisTask("jet-skimmer")}; +} diff --git a/Analysis/Tasks/PWGJE/jetsubstructure.cxx b/Analysis/Tasks/PWGJE/jetsubstructure.cxx new file mode 100644 index 0000000000000..5a51b6e739e03 --- /dev/null +++ b/Analysis/Tasks/PWGJE/jetsubstructure.cxx @@ -0,0 +1,124 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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. + +// jet analysis tasks (subscribing to jet finder task) +// +// Author: Nima Zardoshti +// + +#include "TH1F.h" +#include "TTree.h" + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoA.h" + +#include "Analysis/Jet.h" +#include "Analysis/JetFinder.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +namespace o2::aod +{ +namespace jetsubstructure +{ +DECLARE_SOA_COLUMN(Zg, zg, float); +DECLARE_SOA_COLUMN(Rg, rg, float); +DECLARE_SOA_COLUMN(Nsd, nsd, float); +} // namespace jetsubstructure +DECLARE_SOA_TABLE(JetSubtructure, "AOD", "JETSUBSTRUCTURE", jetsubstructure::Zg, jetsubstructure::Rg, jetsubstructure::Nsd); +} // namespace o2::aod + +struct JetSubstructure { + Produces jetSubstructure; + OutputObj hZg{"h_jet_zg"}; + OutputObj hRg{"h_jet_rg"}; + OutputObj hNsd{"h_jet_nsd"}; + + Configurable f_jetPtMin{"f_jetPtMin", 0.0, "minimum jet pT cut"}; + Configurable f_zCut{"f_zCut", 0.1, "soft drop z cut"}; + Configurable f_beta{"f_beta", 0.0, "soft drop beta"}; + Configurable f_jetR{"f_jetR", 0.4, "jer resolution parameter"}; //possible to get configurable from another task? jetR + Configurable b_DoConstSub{"b_DoConstSub", false, "do constituent subtraction"}; + + std::vector jetConstituents; + std::vector jetReclustered; + JetFinder jetReclusterer; + + void init(InitContext const&) + { + hZg.setObject(new TH1F("h_jet_zg", "zg ;zg", + 10, 0.0, 0.5)); + hRg.setObject(new TH1F("h_jet_rg", "rg ;rg", + 10, 0.0, 0.5)); + hNsd.setObject(new TH1F("h_jet_nsd", "nsd ;nsd", + 7, -0.5, 6.5)); + jetReclusterer.isReclustering = true; + jetReclusterer.algorithm = fastjet::JetAlgorithm::cambridge_algorithm; + jetReclusterer.jetR = f_jetR * 2.5; + } + + //Filter jetCuts = aod::jet::pt > f_jetPtMin; //how does this work? + + void process(aod::Jet const& jet, + aod::Tracks const& tracks, + aod::JetConstituents const& constituents, + aod::JetConstituentsSub const& constituentsSub) + { + jetConstituents.clear(); + jetReclustered.clear(); + if (b_DoConstSub) { + for (const auto constituent : constituentsSub) { + fillConstituents(constituent, jetConstituents); + } + } else { + for (const auto constituentIndex : constituents) { + auto constituent = constituentIndex.track(); + fillConstituents(constituent, jetConstituents); + } + } + fastjet::ClusterSequenceArea clusterSeq(jetReclusterer.findJets(jetConstituents, jetReclustered)); + jetReclustered = sorted_by_pt(jetReclustered); + fastjet::PseudoJet daughterSubJet = jetReclustered[0]; + fastjet::PseudoJet parentSubJet1; + fastjet::PseudoJet parentSubJet2; + bool softDropped = false; + int nsd = 0.0; + auto zg = -1.0; + auto rg = -1.0; + while (daughterSubJet.has_parents(parentSubJet1, parentSubJet2)) { + if (parentSubJet1.perp() < parentSubJet2.perp()) + std::swap(parentSubJet1, parentSubJet2); + auto z = parentSubJet2.perp() / (parentSubJet1.perp() + parentSubJet2.perp()); + auto r = parentSubJet1.delta_R(parentSubJet2); + if (z >= f_zCut * TMath::Power(r / f_jetR, f_beta)) { + if (!softDropped) { + zg = z; + rg = r; + hZg->Fill(zg); + hRg->Fill(rg); + softDropped = true; + } + nsd++; + } + daughterSubJet = parentSubJet1; + } + hNsd->Fill(nsd); + jetSubstructure(zg, rg, nsd); + } +}; +WorkflowSpec defineDataProcessing(ConfigContext const&) +{ + return WorkflowSpec{ + adaptAnalysisTask("jet-substructure")}; +} diff --git a/Analysis/Tasks/jetfinder.cxx b/Analysis/Tasks/jetfinder.cxx deleted file mode 100644 index 7145cd58541c0..0000000000000 --- a/Analysis/Tasks/jetfinder.cxx +++ /dev/null @@ -1,125 +0,0 @@ -// Copyright CERN and copyright holders of ALICE O2. This software is -// distributed under the terms of the GNU General Public License v3 (GPL -// Version 3), copied verbatim in the file "COPYING". -// -// See http://alice-o2.web.cern.ch/license for full licensing information. -// -// 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. - -// jet finder task -// -// Author: Jochen Klein - -#include "Framework/runDataProcessing.h" -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/ASoA.h" - -#include "fastjet/PseudoJet.hh" -#include "fastjet/ClusterSequenceArea.hh" -#include "fastjet/AreaDefinition.hh" -#include "fastjet/JetDefinition.hh" -#include "fastjet/tools/JetMedianBackgroundEstimator.hh" -#include "fastjet/tools/Subtractor.hh" - -#include "Analysis/Jet.h" - -using namespace o2; -using namespace o2::framework; -using namespace o2::framework::expressions; - -struct JetFinderTask { - Produces jets; - Produces constituents; - - // options for jet finding - Configurable rParam{"rParam", 0.4, "jet radius"}; - Configurable ghostEtamax{"ghostEtamax", 1.0, "eta max for ghosts"}; - Configurable minJetPt{"minJetPt", 0., "minimum jet pt"}; - // TODO: initialize with 0.9 - rParam (requires lazy evaluation) - Configurable maxJetEta{"maxJetEta", 0.5, "eta range for jets"}; - // TODO: use configurables also for enums - fastjet::JetAlgorithm algorithm{fastjet::antikt_algorithm}; - fastjet::RecombinationScheme recombScheme{fastjet::E_scheme}; - fastjet::Strategy strategy{fastjet::Best}; - fastjet::AreaType areaType{fastjet::passive_area}; - fastjet::GhostedAreaSpec ghostSpec{ghostEtamax}; - fastjet::AreaDefinition areaDef{areaType, ghostSpec}; - fastjet::JetDefinition jetDef{algorithm, rParam, recombScheme, strategy}; - fastjet::Selector selJet = fastjet::SelectorPtMin(minJetPt) && - fastjet::SelectorAbsRapMax(maxJetEta); - - // options for background subtraction - enum class BkgMode { none, - rhoArea }; - BkgMode bkgMode = BkgMode::none; - Configurable rParamBkg{"rParamBkg", 0.2, "jet radius for background"}; - Configurable rapBkg{"rapBkg", .9, "rapidity range for background"}; - // TODO: use configurables also for enums - fastjet::JetAlgorithm algorithmBkg{fastjet::kt_algorithm}; - fastjet::RecombinationScheme recombSchemeBkg{fastjet::E_scheme}; - fastjet::JetDefinition jetDefBkg{algorithmBkg, rParamBkg, recombSchemeBkg, strategy}; - fastjet::AreaDefinition areaDefBkg{areaType, ghostSpec}; - fastjet::Selector selBkg = fastjet::SelectorAbsRapMax(rapBkg); - - // TODO: use values from configurables - // TODO: add eta cuts - Filter trackCuts = aod::track::pt > 0.1f; - - std::unique_ptr bge; - std::unique_ptr sub; - - std::vector pJets; - - void init(InitContext const&) - { - if (bkgMode == BkgMode::none) { - } else if (bkgMode == BkgMode::rhoArea) { - bge = decltype(bge)(new fastjet::JetMedianBackgroundEstimator(selBkg, jetDefBkg, areaDefBkg)); - sub = decltype(sub){new fastjet::Subtractor{bge.get()}}; - } else { - LOGF(ERROR, "requested subtraction mode not implemented!"); - } - } - - void process(aod::Collision const& collision, - soa::Filtered const& fullTracks) - { - // TODO: retrieve pion mass from somewhere - const float mPionSquared = 0.139 * 0.139; - - pJets.clear(); - - std::vector inputParticles; - for (auto& track : fullTracks) { - auto energy = std::sqrt(track.p() * track.p() + mPionSquared); - inputParticles.emplace_back(track.px(), track.py(), track.pz(), energy); - inputParticles.back().set_user_index(track.globalIndex()); - } - - fastjet::ClusterSequenceArea clust_seq(inputParticles, jetDef, areaDef); - if (bge) - bge->set_particles(inputParticles); - pJets = sub ? (*sub)(clust_seq.inclusive_jets()) : clust_seq.inclusive_jets(); - - pJets = selJet(pJets); - - for (const auto& pjet : pJets) { - jets(collision, pjet.eta(), pjet.phi(), pjet.pt(), - pjet.area(), pjet.Et(), pjet.m()); - for (const auto& track : pjet.constituents()) { - LOGF(DEBUG, "jet %d constituent %d: %f %f %f", jets.lastIndex(), - track.user_index(), track.eta(), track.phi(), track.pt()); - constituents(jets.lastIndex(), track.user_index()); - } - } - } -}; - -WorkflowSpec defineDataProcessing(ConfigContext const&) -{ - return WorkflowSpec{ - adaptAnalysisTask("jet-finder")}; -} From 614857fbafda5bd37501ef32b40cd7f28f53c85a Mon Sep 17 00:00:00 2001 From: Jochen Klein Date: Tue, 6 Oct 2020 13:07:10 +0200 Subject: [PATCH 2/2] Refactor PWGJE code --- Analysis/Core/CMakeLists.txt | 10 +++ .../include/Analysis/JetFinder.h | 80 +++-------------- Analysis/Core/src/AnalysisCoreLinkDef.h | 2 + Analysis/Core/src/AnalysisJetsLinkDef.h | 15 ++++ Analysis/Core/src/JetFinder.cxx | 90 +++++++++++++++++++ Analysis/Tasks/PWGJE/CMakeLists.txt | 15 ++-- .../Tasks/PWGJE/jetfinderhadronrecoil.cxx | 6 +- Analysis/Tasks/PWGJE/jetfinderhf.cxx | 11 +-- 8 files changed, 145 insertions(+), 84 deletions(-) rename Analysis/{DataModel => Core}/include/Analysis/JetFinder.h (72%) create mode 100644 Analysis/Core/src/AnalysisJetsLinkDef.h create mode 100644 Analysis/Core/src/JetFinder.cxx diff --git a/Analysis/Core/CMakeLists.txt b/Analysis/Core/CMakeLists.txt index 03b0393c6029e..145b41d3ed9ce 100644 --- a/Analysis/Core/CMakeLists.txt +++ b/Analysis/Core/CMakeLists.txt @@ -30,3 +30,13 @@ o2_target_root_dictionary(AnalysisCore include/Analysis/TriggerAliases.h include/Analysis/MC.h LINKDEF src/AnalysisCoreLinkDef.h) + +if(FastJet_FOUND) +o2_add_library(AnalysisJets + SOURCES src/JetFinder.cxx + PUBLIC_LINK_LIBRARIES O2::AnalysisCore FastJet::FastJet FastJet::Contrib) + +o2_target_root_dictionary(AnalysisJets + HEADERS include/Analysis/JetFinder.h + LINKDEF src/AnalysisJetsLinkDef.h) +endif() diff --git a/Analysis/DataModel/include/Analysis/JetFinder.h b/Analysis/Core/include/Analysis/JetFinder.h similarity index 72% rename from Analysis/DataModel/include/Analysis/JetFinder.h rename to Analysis/Core/include/Analysis/JetFinder.h index 7ff2aa5cc233b..486d11daedf5e 100644 --- a/Analysis/DataModel/include/Analysis/JetFinder.h +++ b/Analysis/Core/include/Analysis/JetFinder.h @@ -12,6 +12,9 @@ // // Authors: Nima Zardoshti, Jochen Klein +#ifndef O2_ANALYSIS_JETFINDER_H +#define O2_ANALYSIS_JETFINDER_H + #include #include #include @@ -26,8 +29,6 @@ #include -float mPion = TDatabasePDG::Instance()->GetParticle(211)->Mass(); //can be removed when pion mass becomes default for unidentified tracks - class JetFinder { @@ -46,6 +47,8 @@ class JetFinder /// \return ClusterSequenceArea object needed to access constituents // fastjet::ClusterSequenceArea findJets(std::vector &inputParticles, std::vector &jets); + static constexpr float mPion = 0.139; // TDatabasePDG::Instance()->GetParticle(211)->Mass(); //can be removed when pion mass becomes default for unidentified tracks + float phiMin; float phiMax; float etaMin; @@ -142,79 +145,20 @@ class JetFinder ~JetFinder() = default; /// Sets the jet finding parameters - void setParams() - { - - if (!isReclustering) { - jetEtaMin = etaMin + jetR; //in aliphysics this was (-etaMax + 0.95*jetR) - jetEtaMax = etaMax - jetR; - } - - //selGhosts =fastjet::SelectorRapRange(ghostEtaMin,ghostEtaMax) && fastjet::SelectorPhiRange(phiMin,phiMax); - //ghostAreaSpec=fastjet::GhostedAreaSpec(selGhosts,ghostRepeatN,ghostArea,gridScatter,ktScatter,ghostktMean); - ghostAreaSpec = fastjet::GhostedAreaSpec(ghostEtaMax, ghostRepeatN, ghostArea, gridScatter, ktScatter, ghostktMean); //the first argument is rapidity not pseudorapidity, to be checked - jetDef = fastjet::JetDefinition(algorithm, jetR, recombScheme, strategy); - areaDef = fastjet::AreaDefinition(areaType, ghostAreaSpec); - selJets = fastjet::SelectorPtRange(jetPtMin, jetPtMax) && fastjet::SelectorEtaRange(jetEtaMin, jetEtaMax) && fastjet::SelectorPhiRange(jetPhiMin, jetPhiMax); - jetDefBkg = fastjet::JetDefinition(algorithmBkg, jetBkgR, recombSchemeBkg, strategyBkg); - areaDefBkg = fastjet::AreaDefinition(areaTypeBkg, ghostAreaSpec); - selRho = fastjet::SelectorRapRange(bkgEtaMin, bkgEtaMax) && fastjet::SelectorPhiRange(bkgPhiMin, bkgPhiMax); //&& !fastjet::SelectorNHardest(2) //here we have to put rap range, to be checked! - } + void setParams(); /// Sets the background subtraction estimater pointer - void setBkgE() - { - if (bkgSubMode == BkgSubMode::rhoAreaSub || bkgSubMode == BkgSubMode::constSub) { - bkgE = decltype(bkgE)(new fastjet::JetMedianBackgroundEstimator(selRho, jetDefBkg, areaDefBkg)); - } else { - if (bkgSubMode != BkgSubMode::none) - LOGF(ERROR, "requested subtraction mode not implemented!"); - } - } + void setBkgE(); /// Sets the background subtraction pointer - void setSub() - { - //if rho < 1e-6 it is set to 1e-6 in AliPhysics - if (bkgSubMode == BkgSubMode::rhoAreaSub) { - sub = decltype(sub){new fastjet::Subtractor{bkgE.get()}}; - } else if (bkgSubMode == BkgSubMode::constSub) { //event or jetwise - constituentSub = decltype(constituentSub){new fastjet::contrib::ConstituentSubtractor{bkgE.get()}}; - constituentSub->set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR); - constituentSub->set_max_distance(constSubRMax); - constituentSub->set_alpha(constSubAlpha); - constituentSub->set_ghost_area(ghostArea); - constituentSub->set_max_eta(bkgEtaMax); - constituentSub->set_background_estimator(bkgE.get()); //what about rho_m - } else { - if (bkgSubMode != BkgSubMode::none) - LOGF(ERROR, "requested subtraction mode not implemented!"); - } - } + void setSub(); /// Performs jet finding /// \note the input particle and jet lists are passed by reference /// \param inputParticles vector of input particles/tracks /// \param jets veector of jets to be filled /// \return ClusterSequenceArea object needed to access constituents - fastjet::ClusterSequenceArea findJets(std::vector& inputParticles, std::vector& jets) //ideally find a way of passing the cluster sequence as a reeference - { - setParams(); - setBkgE(); - jets.clear(); - - if (bkgE) { - bkgE->set_particles(inputParticles); - setSub(); - } - if (constituentSub) { - inputParticles = constituentSub->subtract_event(inputParticles); - } - fastjet::ClusterSequenceArea clusterSeq(inputParticles, jetDef, areaDef); - jets = sub ? (*sub)(clusterSeq.inclusive_jets()) : clusterSeq.inclusive_jets(); - jets = selJets(jets); - return clusterSeq; - } + fastjet::ClusterSequenceArea findJets(std::vector& inputParticles, std::vector& jets); //ideally find a way of passing the cluster sequence as a reeference private: //void setParams(); @@ -222,6 +166,8 @@ class JetFinder std::unique_ptr bkgE; std::unique_ptr sub; std::unique_ptr constituentSub; + + ClassDef(JetFinder, 1); }; //does this belong here? @@ -229,6 +175,8 @@ template void fillConstituents(const T& constituent, std::vector& constituents) { - auto energy = std::sqrt(constituent.p() * constituent.p() + mPion * mPion); + auto energy = std::sqrt(constituent.p() * constituent.p() + JetFinder::mPion * JetFinder::mPion); constituents.emplace_back(constituent.px(), constituent.py(), constituent.pz(), energy); } + +#endif diff --git a/Analysis/Core/src/AnalysisCoreLinkDef.h b/Analysis/Core/src/AnalysisCoreLinkDef.h index 5eac2f6879fa5..bc79d27b607be 100644 --- a/Analysis/Core/src/AnalysisCoreLinkDef.h +++ b/Analysis/Core/src/AnalysisCoreLinkDef.h @@ -26,3 +26,5 @@ #pragma link C++ class HistogramManager + ; #pragma link C++ class AnalysisCut + ; #pragma link C++ class AnalysisCompositeCut + ; + +// #pragma link C++ class JetFinder+; diff --git a/Analysis/Core/src/AnalysisJetsLinkDef.h b/Analysis/Core/src/AnalysisJetsLinkDef.h new file mode 100644 index 0000000000000..8de8bce16da64 --- /dev/null +++ b/Analysis/Core/src/AnalysisJetsLinkDef.h @@ -0,0 +1,15 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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. + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class JetFinder + ; diff --git a/Analysis/Core/src/JetFinder.cxx b/Analysis/Core/src/JetFinder.cxx new file mode 100644 index 0000000000000..fb1498d911e80 --- /dev/null +++ b/Analysis/Core/src/JetFinder.cxx @@ -0,0 +1,90 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// 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. + +// jet finder task +// +// Author: Jochen Klein, Nima Zardoshti +#include "Analysis/JetFinder.h" +#include "Framework/AnalysisTask.h" + +/// Sets the jet finding parameters +void JetFinder::setParams() +{ + + if (!isReclustering) { + jetEtaMin = etaMin + jetR; //in aliphysics this was (-etaMax + 0.95*jetR) + jetEtaMax = etaMax - jetR; + } + + //selGhosts =fastjet::SelectorRapRange(ghostEtaMin,ghostEtaMax) && fastjet::SelectorPhiRange(phiMin,phiMax); + //ghostAreaSpec=fastjet::GhostedAreaSpec(selGhosts,ghostRepeatN,ghostArea,gridScatter,ktScatter,ghostktMean); + ghostAreaSpec = fastjet::GhostedAreaSpec(ghostEtaMax, ghostRepeatN, ghostArea, gridScatter, ktScatter, ghostktMean); //the first argument is rapidity not pseudorapidity, to be checked + jetDef = fastjet::JetDefinition(algorithm, jetR, recombScheme, strategy); + areaDef = fastjet::AreaDefinition(areaType, ghostAreaSpec); + selJets = fastjet::SelectorPtRange(jetPtMin, jetPtMax) && fastjet::SelectorEtaRange(jetEtaMin, jetEtaMax) && fastjet::SelectorPhiRange(jetPhiMin, jetPhiMax); + jetDefBkg = fastjet::JetDefinition(algorithmBkg, jetBkgR, recombSchemeBkg, strategyBkg); + areaDefBkg = fastjet::AreaDefinition(areaTypeBkg, ghostAreaSpec); + selRho = fastjet::SelectorRapRange(bkgEtaMin, bkgEtaMax) && fastjet::SelectorPhiRange(bkgPhiMin, bkgPhiMax); //&& !fastjet::SelectorNHardest(2) //here we have to put rap range, to be checked! +} + +/// Sets the background subtraction estimater pointer +void JetFinder::setBkgE() +{ + if (bkgSubMode == BkgSubMode::rhoAreaSub || bkgSubMode == BkgSubMode::constSub) { + bkgE = decltype(bkgE)(new fastjet::JetMedianBackgroundEstimator(selRho, jetDefBkg, areaDefBkg)); + } else { + if (bkgSubMode != BkgSubMode::none) + LOGF(ERROR, "requested subtraction mode not implemented!"); + } +} + +/// Sets the background subtraction pointer +void JetFinder::setSub() +{ + //if rho < 1e-6 it is set to 1e-6 in AliPhysics + if (bkgSubMode == BkgSubMode::rhoAreaSub) { + sub = decltype(sub){new fastjet::Subtractor{bkgE.get()}}; + } else if (bkgSubMode == BkgSubMode::constSub) { //event or jetwise + constituentSub = decltype(constituentSub){new fastjet::contrib::ConstituentSubtractor{bkgE.get()}}; + constituentSub->set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR); + constituentSub->set_max_distance(constSubRMax); + constituentSub->set_alpha(constSubAlpha); + constituentSub->set_ghost_area(ghostArea); + constituentSub->set_max_eta(bkgEtaMax); + constituentSub->set_background_estimator(bkgE.get()); //what about rho_m + } else { + if (bkgSubMode != BkgSubMode::none) + LOGF(ERROR, "requested subtraction mode not implemented!"); + } +} + +/// Performs jet finding +/// \note the input particle and jet lists are passed by reference +/// \param inputParticles vector of input particles/tracks +/// \param jets veector of jets to be filled +/// \return ClusterSequenceArea object needed to access constituents +fastjet::ClusterSequenceArea JetFinder::findJets(std::vector& inputParticles, std::vector& jets) //ideally find a way of passing the cluster sequence as a reeference +{ + setParams(); + setBkgE(); + jets.clear(); + + if (bkgE) { + bkgE->set_particles(inputParticles); + setSub(); + } + if (constituentSub) { + inputParticles = constituentSub->subtract_event(inputParticles); + } + fastjet::ClusterSequenceArea clusterSeq(inputParticles, jetDef, areaDef); + jets = sub ? (*sub)(clusterSeq.inclusive_jets()) : clusterSeq.inclusive_jets(); + jets = selJets(jets); + return clusterSeq; +} diff --git a/Analysis/Tasks/PWGJE/CMakeLists.txt b/Analysis/Tasks/PWGJE/CMakeLists.txt index 889fef51ba988..d0afc844d3375 100644 --- a/Analysis/Tasks/PWGJE/CMakeLists.txt +++ b/Analysis/Tasks/PWGJE/CMakeLists.txt @@ -12,32 +12,27 @@ if(FastJet_FOUND) o2_add_dpl_workflow(jet-finder SOURCES jetfinder.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisCore O2::AnalysisDataModel - FastJet::FastJet FastJet::ConstituentSubtractor + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisJets O2::AnalysisDataModel COMPONENT_NAME Analysis) o2_add_dpl_workflow(jet-finder-hf SOURCES jetfinderhf.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisCore O2::AnalysisDataModel - FastJet::FastJet FastJet::ConstituentSubtractor + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisJets O2::AnalysisDataModel COMPONENT_NAME Analysis) o2_add_dpl_workflow(jet-finder-hadron-recoil SOURCES jetfinderhadronrecoil.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisCore O2::AnalysisDataModel - FastJet::FastJet FastJet::ConstituentSubtractor + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisJets O2::AnalysisDataModel COMPONENT_NAME Analysis) o2_add_dpl_workflow(jet-substructure SOURCES jetsubstructure.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisCore O2::AnalysisDataModel - FastJet::FastJet FastJet::ConstituentSubtractor + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisJets O2::AnalysisDataModel COMPONENT_NAME Analysis) o2_add_dpl_workflow(jet-skimmer SOURCES jetskimming.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisCore O2::AnalysisDataModel - FastJet::FastJet FastJet::ConstituentSubtractor + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisJets O2::AnalysisDataModel COMPONENT_NAME Analysis) endif() diff --git a/Analysis/Tasks/PWGJE/jetfinderhadronrecoil.cxx b/Analysis/Tasks/PWGJE/jetfinderhadronrecoil.cxx index d85576b3be959..d5eda89cbca71 100644 --- a/Analysis/Tasks/PWGJE/jetfinderhadronrecoil.cxx +++ b/Analysis/Tasks/PWGJE/jetfinderhadronrecoil.cxx @@ -10,7 +10,7 @@ // jet finder task // -// Author: Jochen Klein +// Author: Nima Zardoshti #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" @@ -34,7 +34,7 @@ struct JetFinderHadronRecoilTask { Configurable f_trackTTMin{"f_trackTTMin", 8.0, "TT hadron min pT"}; Configurable f_trackTTMax{"f_trackTTMax", 9.0, "TT hadron max pT"}; - Configurable f_recoilWindow{"f_recoilWindow", 0.6, "jet finding reecoil window phi distance from TT hadrom"}; + Configurable f_recoilWindow{"f_recoilWindow", 0.6, "jet finding phi window reecoilling from hadron"}; Filter trackCuts = aod::track::pt >= 0.15f && aod::track::eta > -0.9f && aod::track::eta < 0.9f; int collisionSplit = 0; //can we partition the collisions? @@ -96,7 +96,7 @@ struct JetFinderHadronRecoilTask { trackTTPhi = track.phi(); } } - auto energy = std::sqrt(track.p() * track.p() + mPion * mPion); + auto energy = std::sqrt(track.p() * track.p() + JetFinder::mPion * JetFinder::mPion); inputParticles.emplace_back(track.px(), track.py(), track.pz(), energy); inputParticles.back().set_user_index(track.globalIndex()); } diff --git a/Analysis/Tasks/PWGJE/jetfinderhf.cxx b/Analysis/Tasks/PWGJE/jetfinderhf.cxx index 7c25428caa3ec..d32a7018e94b4 100644 --- a/Analysis/Tasks/PWGJE/jetfinderhf.cxx +++ b/Analysis/Tasks/PWGJE/jetfinderhf.cxx @@ -17,8 +17,8 @@ #include "Framework/AnalysisDataModel.h" #include "Framework/ASoA.h" -#include "Analysis/SecondaryVertexHF.h" -#include "Analysis/CandidateSelectionTables.h" +#include "Analysis/HFSecondaryVertex.h" +#include "Analysis/HFCandidateSelectionTables.h" #include "fastjet/PseudoJet.hh" #include "fastjet/ClusterSequenceArea.hh" @@ -55,10 +55,11 @@ struct JetFinderHFTask { enum pdgCode { pdgD0 = 421 }; Filter trackCuts = (aod::track::pt > 0.15f && aod::track::eta > -0.9f && aod::track::eta < 0.9f); - Filter seltrack = (aod::hfselcandidate::isSelD0 >= d_selectionFlagD0 || aod::hfselcandidate::isSelD0bar >= d_selectionFlagD0bar); + Filter seltrack = (aod::hf_selcandidate::isSelD0 >= d_selectionFlagD0 || aod::hf_selcandidate::isSelD0bar >= d_selectionFlagD0bar); void process(aod::Collision const& collision, - soa::Filtered const& tracks, soa::Filtered> const& candidates) + soa::Filtered const& tracks, + soa::Filtered> const& candidates) { std::cout << "Per Event" << std::endl; // TODO: retrieve pion mass from somewhere @@ -69,7 +70,7 @@ struct JetFinderHFTask { jets.clear(); inputParticles.clear(); for (auto& track : tracks) { - auto energy = std::sqrt(track.p() * track.p() + mPion * mPion); + auto energy = std::sqrt(track.p() * track.p() + JetFinder::mPion * JetFinder::mPion); if (candidate.index0().globalIndex() == track.globalIndex() || candidate.index1().globalIndex() == track.globalIndex()) { //is it global index? continue; }