Skip to content

Commit 341158b

Browse files
nzardoshnzardosh
authored
restructuring PWGJE (#4512)
1 parent b33c799 commit 341158b

14 files changed

Lines changed: 921 additions & 140 deletions

File tree

‎Analysis/Core/CMakeLists.txt‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,3 +30,13 @@ o2_target_root_dictionary(AnalysisCore
3030
include/Analysis/TriggerAliases.h
3131
include/Analysis/MC.h
3232
LINKDEF src/AnalysisCoreLinkDef.h)
33+
34+
if(FastJet_FOUND)
35+
o2_add_library(AnalysisJets
36+
SOURCES src/JetFinder.cxx
37+
PUBLIC_LINK_LIBRARIES O2::AnalysisCore FastJet::FastJet FastJet::Contrib)
38+
39+
o2_target_root_dictionary(AnalysisJets
40+
HEADERS include/Analysis/JetFinder.h
41+
LINKDEF src/AnalysisJetsLinkDef.h)
42+
endif()

‎Analysis/Core/include/Analysis/JetFinder.h‎

Lines changed: 182 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,182 @@
1+
// Copyright CERN and copyright holders of ALICE O2. This software is
2+
// distributed under the terms of the GNU General Public License v3 (GPL
3+
// Version 3), copied verbatim in the file "COPYING".
4+
//
5+
// See http://alice-o2.web.cern.ch/license for full licensing information.
6+
//
7+
// In applying this license CERN does not waive the privileges and immunities
8+
// granted to it by virtue of its status as an Intergovernmental Organization
9+
// or submit itself to any jurisdiction.
10+
11+
// jet finder task
12+
//
13+
// Authors: Nima Zardoshti, Jochen Klein
14+
15+
#ifndef O2_ANALYSIS_JETFINDER_H
16+
#define O2_ANALYSIS_JETFINDER_H
17+
18+
#include <TDatabasePDG.h>
19+
#include <TPDGCode.h>
20+
#include <TMath.h>
21+
22+
#include "fastjet/PseudoJet.hh"
23+
#include "fastjet/ClusterSequenceArea.hh"
24+
#include "fastjet/AreaDefinition.hh"
25+
#include "fastjet/JetDefinition.hh"
26+
#include "fastjet/tools/JetMedianBackgroundEstimator.hh"
27+
#include "fastjet/tools/Subtractor.hh"
28+
#include "fastjet/contrib/ConstituentSubtractor.hh"
29+
30+
#include <vector>
31+
32+
class JetFinder
33+
{
34+
35+
public:
36+
enum class BkgSubMode { none,
37+
rhoAreaSub,
38+
constSub };
39+
BkgSubMode bkgSubMode;
40+
41+
void setBkgSubMode(BkgSubMode bSM) { bkgSubMode = bSM; }
42+
43+
/// Performs jet finding
44+
/// \note the input particle and jet lists are passed by reference
45+
/// \param inputParticles vector of input particles/tracks
46+
/// \param jets veector of jets to be filled
47+
/// \return ClusterSequenceArea object needed to access constituents
48+
// fastjet::ClusterSequenceArea findJets(std::vector<fastjet::PseudoJet> &inputParticles, std::vector<fastjet::PseudoJet> &jets);
49+
50+
static constexpr float mPion = 0.139; // TDatabasePDG::Instance()->GetParticle(211)->Mass(); //can be removed when pion mass becomes default for unidentified tracks
51+
52+
float phiMin;
53+
float phiMax;
54+
float etaMin;
55+
float etaMax;
56+
57+
float jetR;
58+
float jetPtMin;
59+
float jetPtMax;
60+
float jetPhiMin;
61+
float jetPhiMax;
62+
float jetEtaMin;
63+
float jetEtaMax;
64+
65+
float ghostEtaMin;
66+
float ghostEtaMax;
67+
float ghostArea;
68+
int ghostRepeatN;
69+
double ghostktMean;
70+
float gridScatter;
71+
float ktScatter;
72+
73+
float jetBkgR;
74+
float bkgPhiMin;
75+
float bkgPhiMax;
76+
float bkgEtaMin;
77+
float bkgEtaMax;
78+
79+
float constSubAlpha;
80+
float constSubRMax;
81+
82+
bool isReclustering;
83+
84+
fastjet::JetAlgorithm algorithm;
85+
fastjet::RecombinationScheme recombScheme;
86+
fastjet::Strategy strategy;
87+
fastjet::AreaType areaType;
88+
fastjet::GhostedAreaSpec ghostAreaSpec;
89+
fastjet::JetDefinition jetDef;
90+
fastjet::AreaDefinition areaDef;
91+
fastjet::Selector selJets;
92+
fastjet::Selector selGhosts;
93+
94+
fastjet::JetAlgorithm algorithmBkg;
95+
fastjet::RecombinationScheme recombSchemeBkg;
96+
fastjet::Strategy strategyBkg;
97+
fastjet::AreaType areaTypeBkg;
98+
fastjet::JetDefinition jetDefBkg;
99+
fastjet::AreaDefinition areaDefBkg;
100+
fastjet::Selector selRho;
101+
102+
/// Default constructor
103+
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),
104+
phiMax(phi_Max),
105+
etaMin(eta_Min),
106+
etaMax(eta_Max),
107+
jetR(0.4),
108+
jetPtMin(0.0),
109+
jetPtMax(1000.0),
110+
jetPhiMin(phi_Min),
111+
jetPhiMax(phi_Max),
112+
jetEtaMin(eta_Min),
113+
jetEtaMax(eta_Max),
114+
ghostEtaMin(eta_Min),
115+
ghostEtaMax(eta_Max),
116+
ghostArea(0.005),
117+
ghostRepeatN(1),
118+
ghostktMean(1e-100), //is float precise enough?
119+
gridScatter(1.0),
120+
ktScatter(0.1),
121+
jetBkgR(0.2),
122+
bkgPhiMin(phi_Min),
123+
bkgPhiMax(phi_Max),
124+
bkgEtaMin(eta_Min),
125+
bkgEtaMax(eta_Max),
126+
algorithm(fastjet::antikt_algorithm),
127+
recombScheme(fastjet::E_scheme),
128+
strategy(fastjet::Best),
129+
areaType(fastjet::active_area),
130+
algorithmBkg(fastjet::JetAlgorithm(fastjet::kt_algorithm)),
131+
recombSchemeBkg(fastjet::RecombinationScheme(fastjet::E_scheme)),
132+
strategyBkg(fastjet::Best),
133+
areaTypeBkg(fastjet::active_area),
134+
bkgSubMode(BkgSubMode::none),
135+
constSubAlpha(1.0),
136+
constSubRMax(0.6),
137+
isReclustering(false)
138+
139+
{
140+
141+
//default constructor
142+
}
143+
144+
/// Default destructor
145+
~JetFinder() = default;
146+
147+
/// Sets the jet finding parameters
148+
void setParams();
149+
150+
/// Sets the background subtraction estimater pointer
151+
void setBkgE();
152+
153+
/// Sets the background subtraction pointer
154+
void setSub();
155+
156+
/// Performs jet finding
157+
/// \note the input particle and jet lists are passed by reference
158+
/// \param inputParticles vector of input particles/tracks
159+
/// \param jets veector of jets to be filled
160+
/// \return ClusterSequenceArea object needed to access constituents
161+
fastjet::ClusterSequenceArea findJets(std::vector<fastjet::PseudoJet>& inputParticles, std::vector<fastjet::PseudoJet>& jets); //ideally find a way of passing the cluster sequence as a reeference
162+
163+
private:
164+
//void setParams();
165+
//void setBkgSub();
166+
std::unique_ptr<fastjet::BackgroundEstimatorBase> bkgE;
167+
std::unique_ptr<fastjet::Subtractor> sub;
168+
std::unique_ptr<fastjet::contrib::ConstituentSubtractor> constituentSub;
169+
170+
ClassDef(JetFinder, 1);
171+
};
172+
173+
//does this belong here?
174+
template <typename T>
175+
void fillConstituents(const T& constituent, std::vector<fastjet::PseudoJet>& constituents)
176+
{
177+
178+
auto energy = std::sqrt(constituent.p() * constituent.p() + JetFinder::mPion * JetFinder::mPion);
179+
constituents.emplace_back(constituent.px(), constituent.py(), constituent.pz(), energy);
180+
}
181+
182+
#endif

‎Analysis/Core/src/AnalysisCoreLinkDef.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,3 +26,5 @@
2626
#pragma link C++ class HistogramManager + ;
2727
#pragma link C++ class AnalysisCut + ;
2828
#pragma link C++ class AnalysisCompositeCut + ;
29+
30+
// #pragma link C++ class JetFinder+;

‎Analysis/Core/src/AnalysisJetsLinkDef.h‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
// Copyright CERN and copyright holders of ALICE O2. This software is
2+
// distributed under the terms of the GNU General Public License v3 (GPL
3+
// Version 3), copied verbatim in the file "COPYING".
4+
//
5+
// See http://alice-o2.web.cern.ch/license for full licensing information.
6+
//
7+
// In applying this license CERN does not waive the privileges and immunities
8+
// granted to it by virtue of its status as an Intergovernmental Organization
9+
// or submit itself to any jurisdiction.
10+
11+
#pragma link off all globals;
12+
#pragma link off all classes;
13+
#pragma link off all functions;
14+
15+
#pragma link C++ class JetFinder + ;

‎Analysis/Core/src/JetFinder.cxx‎

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
// Copyright CERN and copyright holders of ALICE O2. This software is
2+
// distributed under the terms of the GNU General Public License v3 (GPL
3+
// Version 3), copied verbatim in the file "COPYING".
4+
//
5+
// See http://alice-o2.web.cern.ch/license for full licensing information.
6+
//
7+
// In applying this license CERN does not waive the privileges and immunities
8+
// granted to it by virtue of its status as an Intergovernmental Organization
9+
// or submit itself to any jurisdiction.
10+
11+
// jet finder task
12+
//
13+
// Author: Jochen Klein, Nima Zardoshti
14+
#include "Analysis/JetFinder.h"
15+
#include "Framework/AnalysisTask.h"
16+
17+
/// Sets the jet finding parameters
18+
void JetFinder::setParams()
19+
{
20+
21+
if (!isReclustering) {
22+
jetEtaMin = etaMin + jetR; //in aliphysics this was (-etaMax + 0.95*jetR)
23+
jetEtaMax = etaMax - jetR;
24+
}
25+
26+
//selGhosts =fastjet::SelectorRapRange(ghostEtaMin,ghostEtaMax) && fastjet::SelectorPhiRange(phiMin,phiMax);
27+
//ghostAreaSpec=fastjet::GhostedAreaSpec(selGhosts,ghostRepeatN,ghostArea,gridScatter,ktScatter,ghostktMean);
28+
ghostAreaSpec = fastjet::GhostedAreaSpec(ghostEtaMax, ghostRepeatN, ghostArea, gridScatter, ktScatter, ghostktMean); //the first argument is rapidity not pseudorapidity, to be checked
29+
jetDef = fastjet::JetDefinition(algorithm, jetR, recombScheme, strategy);
30+
areaDef = fastjet::AreaDefinition(areaType, ghostAreaSpec);
31+
selJets = fastjet::SelectorPtRange(jetPtMin, jetPtMax) && fastjet::SelectorEtaRange(jetEtaMin, jetEtaMax) && fastjet::SelectorPhiRange(jetPhiMin, jetPhiMax);
32+
jetDefBkg = fastjet::JetDefinition(algorithmBkg, jetBkgR, recombSchemeBkg, strategyBkg);
33+
areaDefBkg = fastjet::AreaDefinition(areaTypeBkg, ghostAreaSpec);
34+
selRho = fastjet::SelectorRapRange(bkgEtaMin, bkgEtaMax) && fastjet::SelectorPhiRange(bkgPhiMin, bkgPhiMax); //&& !fastjet::SelectorNHardest(2) //here we have to put rap range, to be checked!
35+
}
36+
37+
/// Sets the background subtraction estimater pointer
38+
void JetFinder::setBkgE()
39+
{
40+
if (bkgSubMode == BkgSubMode::rhoAreaSub || bkgSubMode == BkgSubMode::constSub) {
41+
bkgE = decltype(bkgE)(new fastjet::JetMedianBackgroundEstimator(selRho, jetDefBkg, areaDefBkg));
42+
} else {
43+
if (bkgSubMode != BkgSubMode::none)
44+
LOGF(ERROR, "requested subtraction mode not implemented!");
45+
}
46+
}
47+
48+
/// Sets the background subtraction pointer
49+
void JetFinder::setSub()
50+
{
51+
//if rho < 1e-6 it is set to 1e-6 in AliPhysics
52+
if (bkgSubMode == BkgSubMode::rhoAreaSub) {
53+
sub = decltype(sub){new fastjet::Subtractor{bkgE.get()}};
54+
} else if (bkgSubMode == BkgSubMode::constSub) { //event or jetwise
55+
constituentSub = decltype(constituentSub){new fastjet::contrib::ConstituentSubtractor{bkgE.get()}};
56+
constituentSub->set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR);
57+
constituentSub->set_max_distance(constSubRMax);
58+
constituentSub->set_alpha(constSubAlpha);
59+
constituentSub->set_ghost_area(ghostArea);
60+
constituentSub->set_max_eta(bkgEtaMax);
61+
constituentSub->set_background_estimator(bkgE.get()); //what about rho_m
62+
} else {
63+
if (bkgSubMode != BkgSubMode::none)
64+
LOGF(ERROR, "requested subtraction mode not implemented!");
65+
}
66+
}
67+
68+
/// Performs jet finding
69+
/// \note the input particle and jet lists are passed by reference
70+
/// \param inputParticles vector of input particles/tracks
71+
/// \param jets veector of jets to be filled
72+
/// \return ClusterSequenceArea object needed to access constituents
73+
fastjet::ClusterSequenceArea JetFinder::findJets(std::vector<fastjet::PseudoJet>& inputParticles, std::vector<fastjet::PseudoJet>& jets) //ideally find a way of passing the cluster sequence as a reeference
74+
{
75+
setParams();
76+
setBkgE();
77+
jets.clear();
78+
79+
if (bkgE) {
80+
bkgE->set_particles(inputParticles);
81+
setSub();
82+
}
83+
if (constituentSub) {
84+
inputParticles = constituentSub->subtract_event(inputParticles);
85+
}
86+
fastjet::ClusterSequenceArea clusterSeq(inputParticles, jetDef, areaDef);
87+
jets = sub ? (*sub)(clusterSeq.inclusive_jets()) : clusterSeq.inclusive_jets();
88+
jets = selJets(jets);
89+
return clusterSeq;
90+
}

‎Analysis/DataModel/include/Analysis/Jet.h‎

Lines changed: 49 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@
1010

1111
// table definitions for jets
1212
//
13-
// Author: Jochen Klein
13+
// Author: Jochen Klein, Nima Zardoshti
1414

1515
#pragma once
1616

@@ -21,17 +21,31 @@ namespace o2::aod
2121
namespace jet
2222
{
2323
DECLARE_SOA_INDEX_COLUMN(Collision, collision);
24+
DECLARE_SOA_COLUMN(Pt, pt, float);
2425
DECLARE_SOA_COLUMN(Eta, eta, float);
2526
DECLARE_SOA_COLUMN(Phi, phi, float);
26-
DECLARE_SOA_COLUMN(Pt, pt, float);
27-
DECLARE_SOA_COLUMN(Area, area, float);
2827
DECLARE_SOA_COLUMN(Energy, energy, float);
2928
DECLARE_SOA_COLUMN(Mass, mass, float);
29+
DECLARE_SOA_COLUMN(Area, area, float);
30+
DECLARE_SOA_DYNAMIC_COLUMN(Px, px, [](float pt, float phi) { return pt * TMath::Cos(phi); });
31+
DECLARE_SOA_DYNAMIC_COLUMN(Py, py, [](float pt, float phi) { return pt * TMath::Sin(phi); });
32+
DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, [](float pt, float eta) { return pt * TMath::SinH(eta); });
33+
DECLARE_SOA_DYNAMIC_COLUMN(P, p, [](float pt, float eta) { return pt * TMath::CosH(eta); }); //absolute p
3034
} // namespace jet
3135

3236
DECLARE_SOA_TABLE(Jets, "AOD", "JET",
33-
o2::soa::Index<>, jet::CollisionId,
34-
jet::Eta, jet::Phi, jet::Pt, jet::Area, jet::Energy, jet::Mass);
37+
o2::soa::Index<>,
38+
jet::CollisionId,
39+
jet::Pt,
40+
jet::Eta,
41+
jet::Phi,
42+
jet::Energy,
43+
jet::Mass,
44+
jet::Area,
45+
jet::Px<jet::Pt, jet::Phi>,
46+
jet::Py<jet::Pt, jet::Phi>,
47+
jet::Pz<jet::Pt, jet::Eta>,
48+
jet::P<jet::Pt, jet::Eta>);
3549

3650
using Jet = Jets::iterator;
3751

@@ -43,8 +57,36 @@ DECLARE_SOA_INDEX_COLUMN(Jet, jet);
4357
DECLARE_SOA_INDEX_COLUMN(Track, track);
4458
} // namespace constituents
4559

46-
DECLARE_SOA_TABLE(JetConstituents, "AOD", "JETCONSTITUENTS",
47-
constituents::JetId, constituents::TrackId);
60+
DECLARE_SOA_TABLE(JetConstituents, "AOD", "CONSTITUENTS",
61+
constituents::JetId,
62+
constituents::TrackId);
4863

4964
using JetConstituent = JetConstituents::iterator;
65+
66+
namespace constituentssub
67+
{
68+
DECLARE_SOA_INDEX_COLUMN(Jet, jet);
69+
DECLARE_SOA_COLUMN(Pt, pt, float);
70+
DECLARE_SOA_COLUMN(Eta, eta, float);
71+
DECLARE_SOA_COLUMN(Phi, phi, float);
72+
DECLARE_SOA_COLUMN(Energy, energy, float);
73+
DECLARE_SOA_COLUMN(Mass, mass, float);
74+
DECLARE_SOA_DYNAMIC_COLUMN(Px, px, [](float pt, float phi) { return pt * TMath::Cos(phi); });
75+
DECLARE_SOA_DYNAMIC_COLUMN(Py, py, [](float pt, float phi) { return pt * TMath::Sin(phi); });
76+
DECLARE_SOA_DYNAMIC_COLUMN(Pz, pz, [](float pt, float eta) { return pt * TMath::SinH(eta); });
77+
DECLARE_SOA_DYNAMIC_COLUMN(P, p, [](float pt, float eta) { return pt * TMath::CosH(eta); }); //absolute p
78+
} //namespace constituentssub
79+
DECLARE_SOA_TABLE(JetConstituentsSub, "AOD", "CONSTITUENTSSUB",
80+
constituentssub::JetId,
81+
constituentssub::Pt,
82+
constituentssub::Eta,
83+
constituentssub::Phi,
84+
constituentssub::Energy,
85+
constituentssub::Mass,
86+
constituentssub::Px<constituentssub::Pt, constituentssub::Phi>,
87+
constituentssub::Py<constituentssub::Pt, constituentssub::Phi>,
88+
constituentssub::Pz<constituentssub::Pt, constituentssub::Eta>,
89+
constituentssub::P<constituentssub::Pt, constituentssub::Eta>);
90+
using JetConstituentSub = JetConstituentsSub::iterator;
91+
5092
} // namespace o2::aod

0 commit comments

Comments
 (0)