Skip to content

Commit 9758ce6

Browse files
hahassan7hahassan7
authored
PWGJE: Implementing background subtraction in the JE framework (AliceO2Group#3587)
* Implementing background subtraction in the JE framework * RhoAreaSub applied on single jet instead of vector * Fix eventConstSub
1 parent 18ddf3b commit 9758ce6

11 files changed

+483
-151
lines changed

‎PWGJE/Core/CMakeLists.txt

+4-1
Original file line numberDiff line numberDiff line change
@@ -11,13 +11,16 @@
1111

1212
if(FastJet_FOUND)
1313
o2physics_add_library(PWGJECore
14-
SOURCES JetFinder.cxx
14+
SOURCES FastJetUtilities.cxx
15+
JetFinder.cxx
16+
JetBkgSubUtils.cxx
1517
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore FastJet::FastJet FastJet::Contrib)
1618

1719
o2physics_target_root_dictionary(PWGJECore
1820
HEADERS JetFinder.h
1921
JetUtilities.h
2022
FastJetUtilities.h
2123
JetTaggingUtilities.h
24+
JetBkgSubUtils.h
2225
LINKDEF PWGJECoreLinkDef.h)
2326
endif()

‎PWGJE/Core/FastJetUtilities.cxx

+38
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
#include "FastJetUtilities.h"
13+
14+
void FastJetUtilities::setFastJetUserInfo(std::vector<fastjet::PseudoJet>& constituents, int index, int status)
15+
{
16+
fastjet_user_info* user_info = new fastjet_user_info(status, index); // FIXME: can setting this as a pointer be avoided?
17+
constituents.back().set_user_info(user_info);
18+
if (index != -99999999) { // FIXME: needed for constituent subtraction as user_info is not propagated, but need to be quite careful to make sure indices dont overlap between tracks, clusters and HF candidates. Current solution might not be optimal
19+
int i = index;
20+
if (status == static_cast<int>(JetConstituentStatus::track)) {
21+
i = i + 1;
22+
}
23+
if (status == static_cast<int>(JetConstituentStatus::cluster)) {
24+
i = -1 * (i + 1);
25+
}
26+
if (status == static_cast<int>(JetConstituentStatus::candidateHF)) {
27+
i = 0;
28+
}
29+
constituents.back().set_user_index(i); // FIXME: needed for constituent subtraction, but need to be quite careful to make sure indices dont overlap between tracks, clusters and HF candidates. Current solution might not be optimal
30+
}
31+
}
32+
33+
// Selector of HF candidates
34+
fastjet::Selector FastJetUtilities::SelectorIsHFCand()
35+
{
36+
// This method is applied on particles or jet constituents only !!!
37+
return fastjet::Selector(new SW_IsHFCand());
38+
}

‎PWGJE/Core/FastJetUtilities.h

+29-19
Original file line numberDiff line numberDiff line change
@@ -22,18 +22,22 @@
2222
#include <numeric>
2323
#include <tuple>
2424
#include <vector>
25+
#include <string>
2526

26-
#include "PWGJE/Core/JetFinder.h"
27+
#include "fastjet/PseudoJet.hh"
28+
#include "fastjet/Selector.hh"
2729

2830
enum class JetConstituentStatus {
2931
track = 0,
3032
cluster = 1,
31-
candidateHF = 2,
33+
candidateHF = 2
3234
};
3335

3436
namespace FastJetUtilities
3537
{
3638

39+
static constexpr float mPion = 0.139; // TDatabasePDG::Instance()->GetParticle(211)->Mass(); //can be removed when pion mass becomes default for unidentified tracks
40+
3741
// Class defined to store additional info which is passed to the FastJet object
3842
class fastjet_user_info : public fastjet::PseudoJet::UserInfoBase
3943
{
@@ -66,24 +70,30 @@ class fastjet_user_info : public fastjet::PseudoJet::UserInfoBase
6670
* @param status status of constituent type
6771
*/
6872

69-
void setFastJetUserInfo(std::vector<fastjet::PseudoJet>& constituents, int index = -99999999, int status = static_cast<int>(JetConstituentStatus::track))
73+
void setFastJetUserInfo(std::vector<fastjet::PseudoJet>& constituents, int index = -99999999, int status = static_cast<int>(JetConstituentStatus::track));
74+
75+
// Class defined to select the HF candidate particle
76+
// This method is applied on particles or jet constituents only !!!
77+
class SW_IsHFCand : public fastjet::SelectorWorker
7078
{
71-
fastjet_user_info* user_info = new fastjet_user_info(status, index); // FIXME: can setting this as a pointer be avoided?
72-
constituents.back().set_user_info(user_info);
73-
if (index != -99999999) { // FIXME: needed for constituent subtraction as user_info is not propagated, but need to be quite careful to make sure indices dont overlap between tracks, clusters and HF candidates. Current solution might not be optimal
74-
int i = index;
75-
if (status == static_cast<int>(JetConstituentStatus::track)) {
76-
i = i + 1;
77-
}
78-
if (status == static_cast<int>(JetConstituentStatus::cluster)) {
79-
i = -1 * (i + 1);
80-
}
81-
if (status == static_cast<int>(JetConstituentStatus::candidateHF)) {
82-
i = 0;
83-
}
84-
constituents.back().set_user_index(i); // FIXME: needed for constituent subtraction, but need to be quite careful to make sure indices dont overlap between tracks, clusters and HF candidates. Current solution might not be optimal
79+
public:
80+
// default ctor
81+
SW_IsHFCand() {}
82+
83+
// the selector's description
84+
std::string description() const
85+
{
86+
return "HF candidate selector";
8587
}
86-
}
88+
89+
bool pass(const fastjet::PseudoJet& p) const
90+
{
91+
return (p.user_info<fastjet_user_info>().getStatus() == static_cast<int>(JetConstituentStatus::candidateHF));
92+
}
93+
};
94+
95+
// Selector of HF candidates
96+
fastjet::Selector SelectorIsHFCand();
8797

8898
/**
8999
* Add track as a pseudojet object to the fastjet vector
@@ -96,7 +106,7 @@ void setFastJetUserInfo(std::vector<fastjet::PseudoJet>& constituents, int index
96106
*/
97107

98108
template <typename T>
99-
void fillTracks(const T& constituent, std::vector<fastjet::PseudoJet>& constituents, int index = -99999999, int status = static_cast<int>(JetConstituentStatus::track), double mass = JetFinder::mPion)
109+
void fillTracks(const T& constituent, std::vector<fastjet::PseudoJet>& constituents, int index = -99999999, int status = static_cast<int>(JetConstituentStatus::track), double mass = mPion)
100110
{
101111
if (status == static_cast<int>(JetConstituentStatus::track) || status == static_cast<int>(JetConstituentStatus::candidateHF)) {
102112
auto p = std::sqrt((constituent.px() * constituent.px()) + (constituent.py() * constituent.py()) + (constituent.pz() * constituent.pz()));

‎PWGJE/Core/JetBkgSubUtils.cxx

+211
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,211 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
// jet finder task
13+
//
14+
// Author: Hadi Hassan, Universiy of Jväskylä, hadi.hassan@cern.ch
15+
#include <memory>
16+
#include <tuple>
17+
#include "Framework/Logger.h"
18+
#include "Common/Core/RecoDecay.h"
19+
#include "PWGJE/Core/JetUtilities.h"
20+
#include "PWGJE/Core/JetBkgSubUtils.h"
21+
#include "PWGJE/Core/FastJetUtilities.h"
22+
23+
JetBkgSubUtils::JetBkgSubUtils(float jetBkgR_out, float constSubAlpha_out, float constSubRMax_out, float bkgEtaMin_out, float bkgEtaMax_out, float bkgPhiMin_out, float bkgPhiMax_out, fastjet::GhostedAreaSpec ghostAreaSpec_out, int nHardReject) : jetBkgR(jetBkgR_out),
24+
constSubAlpha(constSubAlpha_out),
25+
constSubRMax(constSubRMax_out),
26+
bkgEtaMin(bkgEtaMin_out),
27+
bkgEtaMax(bkgEtaMax_out),
28+
bkgPhiMin(bkgPhiMin_out),
29+
bkgPhiMax(bkgPhiMax_out),
30+
ghostAreaSpec(ghostAreaSpec_out)
31+
32+
{
33+
// Note: if you are using the PerpCone method you should jetBkgR to be the same as the anit_kt jets R, otherwise use R=0.2
34+
jetDefBkg = fastjet::JetDefinition(algorithmBkg, jetBkgR, recombSchemeBkg, fastjet::Best);
35+
areaDefBkg = fastjet::AreaDefinition(fastjet::active_area_explicit_ghosts, ghostAreaSpec);
36+
selRho = fastjet::SelectorRapRange(bkgEtaMin, bkgEtaMax) && fastjet::SelectorPhiRange(bkgPhiMin, bkgPhiMax) && !fastjet::SelectorNHardest(nHardReject); // here we have to put rap range, to be checked!
37+
}
38+
39+
std::tuple<double, double> JetBkgSubUtils::estimateRhoAreaMedian(const std::vector<fastjet::PseudoJet>& inputParticles, bool doSparseSub)
40+
{
41+
42+
if (inputParticles.size() == 0) {
43+
return std::make_tuple(0.0, 0.0);
44+
}
45+
46+
// cluster the kT jets
47+
fastjet::ClusterSequenceArea clusterSeq(removeHFCand ? selRemoveHFCand(inputParticles) : inputParticles, jetDefBkg, areaDefBkg);
48+
49+
// select jets in detector acceptance
50+
std::vector<fastjet::PseudoJet> alljets = selRho(clusterSeq.inclusive_jets());
51+
52+
double totaljetAreaPhys(0), totalAreaCovered(0);
53+
std::vector<double> rhovector;
54+
std::vector<double> rhoMdvector;
55+
56+
// Fill a vector for pT/area to be used for the median
57+
for (auto& ijet : alljets) {
58+
59+
// Physical area/ Physical jets (no ghost)
60+
if (!clusterSeq.is_pure_ghost(ijet)) {
61+
rhovector.push_back(ijet.perp() / ijet.area());
62+
rhoMdvector.push_back(getMd(ijet) / ijet.area());
63+
64+
totaljetAreaPhys += ijet.area();
65+
}
66+
// Full area
67+
totalAreaCovered += ijet.area();
68+
}
69+
// calculate Rho as the median of the jet pT / jet area
70+
double rho = TMath::Median<double>(rhovector.size(), rhovector.data());
71+
double rhoM = TMath::Median<double>(rhoMdvector.size(), rhoMdvector.data());
72+
73+
if (doSparseSub) {
74+
// calculate The ocupancy factor, which the ratio of covered area / total area
75+
double occupancyFactor = totalAreaCovered > 0 ? totaljetAreaPhys / totalAreaCovered : 1.;
76+
rho *= occupancyFactor;
77+
rhoM *= occupancyFactor;
78+
}
79+
80+
return std::make_tuple(rho, rhoM);
81+
}
82+
83+
std::tuple<double, double> JetBkgSubUtils::estimateRhoPerpCone(const std::vector<fastjet::PseudoJet>& inputParticles, const std::vector<fastjet::PseudoJet>& jets)
84+
{
85+
86+
if (inputParticles.size() == 0 || jets.size() == 0) {
87+
return std::make_tuple(0.0, 0.0);
88+
}
89+
90+
// Select a list of particles without the HF candidate
91+
std::vector<fastjet::PseudoJet> inputPartnoHF = removeHFCand ? selRemoveHFCand(inputParticles) : inputParticles;
92+
93+
double perpPtDensity1 = 0;
94+
double perpPtDensity2 = 0;
95+
double perpMdDensity1 = 0;
96+
double perpMdDensity2 = 0;
97+
98+
fastjet::Selector selectJet = fastjet::SelectorEtaRange(bkgEtaMin, bkgEtaMax) && fastjet::SelectorPhiRange(bkgPhiMin, bkgPhiMax);
99+
100+
std::vector<fastjet::PseudoJet> selectedJets = fastjet::sorted_by_pt(selectJet(jets));
101+
102+
if (selectedJets.size() == 0) {
103+
return std::make_tuple(0.0, 0.0);
104+
}
105+
106+
fastjet::PseudoJet leadingJet = selectedJets[0];
107+
108+
double dPhi1 = 999.;
109+
double dPhi2 = 999.;
110+
double dEta = 999.;
111+
double PerpendicularConeAxisPhi1 = 999., PerpendicularConeAxisPhi2 = 999.;
112+
// build 2 perp cones in phi around the leading jet (right and left of the jet)
113+
PerpendicularConeAxisPhi1 = RecoDecay::constrainAngle<double, double>(leadingJet.phi() + (M_PI / 2.)); // This will contrain the angel between 0-2Pi
114+
PerpendicularConeAxisPhi2 = RecoDecay::constrainAngle<double, double>(leadingJet.phi() - (M_PI / 2.)); // This will contrain the angel between 0-2Pi
115+
116+
for (auto& particle : inputPartnoHF) {
117+
// sum the momentum of all paricles that fill the two cones
118+
dPhi1 = particle.phi() - PerpendicularConeAxisPhi1;
119+
dPhi1 = RecoDecay::constrainAngle<double, double>(dPhi1, -M_PI); // This will contrain the angel between -pi & Pi
120+
dPhi2 = particle.phi() - PerpendicularConeAxisPhi2;
121+
dPhi2 = RecoDecay::constrainAngle<double, double>(dPhi2, -M_PI); // This will contrain the angel between -pi & Pi
122+
dEta = leadingJet.eta() - particle.eta(); // The perp cone eta is the same as the leading jet since the cones are perpendicular only in phi
123+
if (TMath::Sqrt(dPhi1 * dPhi1 + dEta * dEta) <= jetBkgR) {
124+
perpPtDensity1 += particle.perp();
125+
perpMdDensity1 += TMath::Sqrt(particle.m() * particle.m() + particle.pt() * particle.pt()) - particle.pt();
126+
}
127+
128+
if (TMath::Sqrt(dPhi2 * dPhi2 + dEta * dEta) <= jetBkgR) {
129+
perpPtDensity2 += particle.perp();
130+
perpMdDensity2 += TMath::Sqrt(particle.m() * particle.m() + particle.pt() * particle.pt()) - particle.pt();
131+
}
132+
}
133+
134+
// Caculate rho as the ratio of average pT of the two cones / the cone area
135+
double perpPtDensity = (perpPtDensity1 + perpPtDensity2) / (2 * M_PI * jetBkgR * jetBkgR);
136+
double perpMdDensity = (perpMdDensity1 + perpMdDensity2) / (2 * M_PI * jetBkgR * jetBkgR);
137+
138+
return std::make_tuple(perpPtDensity, perpMdDensity);
139+
}
140+
141+
fastjet::PseudoJet JetBkgSubUtils::doRhoAreaSub(fastjet::PseudoJet& jet, double& rhoParam, double& rhoMParam)
142+
{
143+
144+
fastjet::Subtractor sub = fastjet::Subtractor(rhoParam, rhoMParam);
145+
if (doRhoMassSub) {
146+
sub.set_safe_mass();
147+
}
148+
return sub(jet);
149+
}
150+
151+
std::vector<fastjet::PseudoJet> JetBkgSubUtils::doEventConstSub(std::vector<fastjet::PseudoJet>& inputParticles, double& rhoParam, double& rhoMParam)
152+
{
153+
154+
fastjet::contrib::ConstituentSubtractor constituentSub(rhoParam, rhoMParam);
155+
constituentSub.set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR); /// deltaR=sqrt((y_i-y_j)^2+(phi_i-phi_j)^2)), longitudinal Lorentz invariant
156+
constituentSub.set_max_distance(constSubRMax);
157+
constituentSub.set_alpha(constSubAlpha);
158+
constituentSub.set_ghost_area(ghostAreaSpec.ghost_area());
159+
constituentSub.set_max_eta(maxEtaEvent);
160+
if (removeHFCand) {
161+
constituentSub.set_particle_selector(&selRemoveHFCand);
162+
}
163+
164+
// by default, the masses of all particles are set to zero. With this flag the jet mass will also be subtracted
165+
if (doRhoMassSub) {
166+
constituentSub.set_do_mass_subtraction();
167+
}
168+
169+
return constituentSub.subtract_event(inputParticles, maxEtaEvent);
170+
}
171+
172+
std::vector<fastjet::PseudoJet> JetBkgSubUtils::doJetConstSub(std::vector<fastjet::PseudoJet>& jets, double& rhoParam, double& rhoMParam)
173+
{
174+
175+
if (jets.size() == 0) {
176+
return std::vector<fastjet::PseudoJet>();
177+
}
178+
179+
// FIXME, this method works only if the input jets "jets" are reconstructed with area def "active_area_explicit_ghosts"
180+
// because it needs the ghosts to estimate the backgeound
181+
fastjet::contrib::ConstituentSubtractor constituentSub(rhoParam, rhoMParam);
182+
constituentSub.set_distance_type(fastjet::contrib::ConstituentSubtractor::deltaR); /// deltaR=sqrt((y_i-y_j)^2+(phi_i-phi_j)^2)), longitudinal Lorentz invariant
183+
constituentSub.set_max_distance(constSubRMax);
184+
constituentSub.set_alpha(constSubAlpha);
185+
constituentSub.set_ghost_area(ghostAreaSpec.ghost_area());
186+
constituentSub.set_max_eta(bkgEtaMax);
187+
if (removeHFCand) {
188+
constituentSub.set_particle_selector(&selRemoveHFCand);
189+
}
190+
191+
// by default, the masses of all particles are set to zero. With this flag the jet mass will also be subtracted
192+
if (doRhoMassSub) {
193+
constituentSub.set_do_mass_subtraction();
194+
}
195+
196+
// FIXME, This method doesn't propagate the area information, since after constituent subtraction
197+
// the jet structure will change, so it no longer has the same area. fastjet developers said calculatig the area
198+
// information will difficult
199+
return constituentSub(jets);
200+
}
201+
202+
double JetBkgSubUtils::getMd(fastjet::PseudoJet jet) const
203+
{
204+
// Refere to https://arxiv.org/abs/1211.2811 for the rhoM caclulation
205+
double sum(0);
206+
for (auto constituent : jet.constituents()) {
207+
sum += TMath::Sqrt(constituent.m() * constituent.m() + constituent.pt() * constituent.pt()) - constituent.pt();
208+
}
209+
210+
return sum;
211+
}

0 commit comments

Comments
 (0)