Skip to content
Closed
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
85 changes: 52 additions & 33 deletions PWGJE/Tasks/jetBackgroundAnalysis.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -8,11 +8,12 @@
// 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 QA task
//
/// \author Aimeric Landou <aimeric.landou@cern.ch>
/// \author Nima Zardoshti <nima.zardoshti@cern.ch>
/// \author Yubiao Wang <yubiao.wang@cern.ch>
/// \file jetBackgroundAnalysis.cxx
/// \brief This file contains the implementation for the Charged Jet v2 analysis in the ALICE experiment

#include "RecoDecay.h"

Expand All @@ -35,20 +36,20 @@
#include <cmath>
#include <string>
#include <vector>

#include <math.h>

using namespace o2;
using namespace o2::framework;
using namespace o2::framework::expressions;

struct JetBackgroundAnalysisTask {
struct JetBackgroundAnalysis {
HistogramRegistry registry;

Configurable<std::string> eventSelections{"eventSelections", "sel8", "choose event selection"};
Configurable<float> vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range for collisions"};
Configurable<float> centralityMin{"centralityMin", -999.0, "minimum centrality for collisions"};
Configurable<float> centralityMax{"centralityMax", 999.0, "maximum centrality for collisions"};
Configurable<bool> checkCentFT0M{"checkCentFT0M", false, "0: centFT0C as default, 1: use centFT0M estimator"};
Configurable<int> trackOccupancyInTimeRangeMax{"trackOccupancyInTimeRangeMax", 999999, "maximum track occupancy of collisions in neighbouring collisions in a given time range; only applied to reconstructed collisions (data and mcd jets), not mc collisions (mcp jets)"};
Configurable<int> trackOccupancyInTimeRangeMin{"trackOccupancyInTimeRangeMin", -999999, "minimum track occupancy of collisions in neighbouring collisions in a given time range; only applied to reconstructed collisions (data and mcd jets), not mc collisions (mcp jets)"};
Configurable<bool> skipMBGapEvents{"skipMBGapEvents", false, "flag to choose to reject min. bias gap events"};
Expand Down Expand Up @@ -99,7 +100,7 @@ struct JetBackgroundAnalysisTask {
}

Filter trackCuts = (aod::jtrack::pt >= trackPtMin && aod::jtrack::pt < trackPtMax && aod::jtrack::eta > trackEtaMin && aod::jtrack::eta < trackEtaMax);
Filter eventCuts = (nabs(aod::jcollision::posZ) < vertexZCut && aod::jcollision::centFT0M >= centralityMin && aod::jcollision::centFT0M < centralityMax);
Filter eventCuts = (nabs(aod::jcollision::posZ) < vertexZCut);

template <typename TTracks, typename TJets>
bool trackIsInJet(TTracks const& track, TJets const& jet)
Expand All @@ -113,73 +114,73 @@ struct JetBackgroundAnalysisTask {
}

template <typename TCollisions, typename TJets, typename TTracks>
void bkgFluctuationsRandomCone(TCollisions const& collision, TJets const& jets, TTracks const& tracks)
void bkgFluctuationsRandomCone(TCollisions const& collision, TJets const& jets, TTracks const& tracks, float centrality)
{
TRandom3 randomNumber(0);
float randomConeEta = randomNumber.Uniform(trackEtaMin + randomConeR, trackEtaMax - randomConeR);
float randomConePhi = randomNumber.Uniform(0.0, 2 * M_PI);
float randomConePhi = randomNumber.Uniform(0.0, o2::constants::math::TwoPI);
float randomConePt = 0;
for (auto const& track : tracks) {
if (jetderiveddatautilities::selectTrack(track, trackSelection)) {
float dPhi = RecoDecay::constrainAngle(track.phi() - randomConePhi, static_cast<float>(-M_PI));
float dPhi = RecoDecay::constrainAngle(track.phi() - randomConePhi, static_cast<float>(-o2::constants::math::PI));
float dEta = track.eta() - randomConeEta;
if (TMath::Sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
randomConePt += track.pt();
}
}
}
registry.fill(HIST("h2_centrality_rhorandomcone"), collision.centFT0M(), randomConePt - M_PI * randomConeR * randomConeR * collision.rho());
registry.fill(HIST("h2_centrality_rhorandomcone"), centrality, randomConePt - o2::constants::math::PI * randomConeR * randomConeR * collision.rho());

// randomised eta,phi for tracks, to assess part of fluctuations coming from statistically independently emitted particles
randomConePt = 0;
for (auto const& track : tracks) {
if (jetderiveddatautilities::selectTrack(track, trackSelection)) {
float dPhi = RecoDecay::constrainAngle(randomNumber.Uniform(0.0, 2 * M_PI) - randomConePhi, static_cast<float>(-M_PI)); // ignores actual phi of track
float dPhi = RecoDecay::constrainAngle(randomNumber.Uniform(0.0, o2::constants::math::TwoPI) - randomConePhi, static_cast<float>(-o2::constants::math::PI)); // ignores actual phi of track
float dEta = randomNumber.Uniform(trackEtaMin, trackEtaMax) - randomConeEta; // ignores actual eta of track
if (TMath::Sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
randomConePt += track.pt();
}
}
}
registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirection"), collision.centFT0M(), randomConePt - M_PI * randomConeR * randomConeR * collision.rho());
registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirection"), centrality, randomConePt - o2::constants::math::PI * randomConeR * randomConeR * collision.rho());

// removing the leading jet from the random cone
if (jets.size() > 0) { // if there are no jets in the acceptance (from the jetfinder cuts) then there can be no leading jet
float dPhiLeadingJet = RecoDecay::constrainAngle(jets.iteratorAt(0).phi() - randomConePhi, static_cast<float>(-M_PI));
float dPhiLeadingJet = RecoDecay::constrainAngle(jets.iteratorAt(0).phi() - randomConePhi, static_cast<float>(-o2::constants::math::PI));
float dEtaLeadingJet = jets.iteratorAt(0).eta() - randomConeEta;

bool jetWasInCone = false;
while ((randomConeLeadJetDeltaR <= 0 && (TMath::Sqrt(dEtaLeadingJet * dEtaLeadingJet + dPhiLeadingJet * dPhiLeadingJet) < jets.iteratorAt(0).r() / 100.0 + randomConeR)) || (randomConeLeadJetDeltaR > 0 && (TMath::Sqrt(dEtaLeadingJet * dEtaLeadingJet + dPhiLeadingJet * dPhiLeadingJet) < randomConeLeadJetDeltaR))) {
while ((randomConeLeadJetDeltaR <= 0 && (std::sqrt(dEtaLeadingJet * dEtaLeadingJet + dPhiLeadingJet * dPhiLeadingJet) < jets.iteratorAt(0).r() / 100.0 + randomConeR)) || (randomConeLeadJetDeltaR > 0 && (std::sqrt(dEtaLeadingJet * dEtaLeadingJet + dPhiLeadingJet * dPhiLeadingJet) < randomConeLeadJetDeltaR))) {
jetWasInCone = true;
randomConeEta = randomNumber.Uniform(trackEtaMin + randomConeR, trackEtaMax - randomConeR);
randomConePhi = randomNumber.Uniform(0.0, 2 * M_PI);
dPhiLeadingJet = RecoDecay::constrainAngle(jets.iteratorAt(0).phi() - randomConePhi, static_cast<float>(-M_PI));
randomConePhi = randomNumber.Uniform(0.0, o2::constants::math::TwoPI);
dPhiLeadingJet = RecoDecay::constrainAngle(jets.iteratorAt(0).phi() - randomConePhi, static_cast<float>(-o2::constants::math::PI));
dEtaLeadingJet = jets.iteratorAt(0).eta() - randomConeEta;
}
if (jetWasInCone) {
randomConePt = 0.0;
for (auto const& track : tracks) {
if (jetderiveddatautilities::selectTrack(track, trackSelection)) { // if track selection is uniformTrack, dcaXY and dcaZ cuts need to be added as they aren't in the selection so that they can be studied here
float dPhi = RecoDecay::constrainAngle(track.phi() - randomConePhi, static_cast<float>(-M_PI));
float dPhi = RecoDecay::constrainAngle(track.phi() - randomConePhi, static_cast<float>(-o2::constants::math::PI));
float dEta = track.eta() - randomConeEta;
if (TMath::Sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
randomConePt += track.pt();
}
}
}
}
}
registry.fill(HIST("h2_centrality_rhorandomconewithoutleadingjet"), collision.centFT0M(), randomConePt - M_PI * randomConeR * randomConeR * collision.rho());
registry.fill(HIST("h2_centrality_rhorandomconewithoutleadingjet"), centrality, randomConePt - o2::constants::math::PI * randomConeR * randomConeR * collision.rho());

// randomised eta,phi for tracks, to assess part of fluctuations coming from statistically independently emitted particles, removing tracks from 2 leading jets
double randomConePtWithoutOneLeadJet = 0;
double randomConePtWithoutTwoLeadJet = 0;
if (jets.size() > 1) { // if there are no jets, or just one, in the acceptance (from the jetfinder cuts) then one cannot find 2 leading jets
for (auto const& track : tracks) {
if (jetderiveddatautilities::selectTrack(track, trackSelection)) {
float dPhi = RecoDecay::constrainAngle(randomNumber.Uniform(0.0, 2 * M_PI) - randomConePhi, static_cast<float>(-M_PI)); // ignores actual phi of track
float dPhi = RecoDecay::constrainAngle(randomNumber.Uniform(0.0, o2::constants::math::TwoPI) - randomConePhi, static_cast<float>(-o2::constants::math::PI)); // ignores actual phi of track
float dEta = randomNumber.Uniform(trackEtaMin, trackEtaMax) - randomConeEta; // ignores actual eta of track
if (TMath::Sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
if (!trackIsInJet(track, jets.iteratorAt(0))) {
randomConePtWithoutOneLeadJet += track.pt();
if (!trackIsInJet(track, jets.iteratorAt(1))) {
Expand All @@ -190,8 +191,8 @@ struct JetBackgroundAnalysisTask {
}
}
}
registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirectionwithoutoneleadingjets"), collision.centFT0M(), randomConePtWithoutOneLeadJet - M_PI * randomConeR * randomConeR * collision.rho());
registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirectionwithouttwoleadingjets"), collision.centFT0M(), randomConePtWithoutTwoLeadJet - M_PI * randomConeR * randomConeR * collision.rho());
registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirectionwithoutoneleadingjets"), centrality, randomConePtWithoutOneLeadJet - o2::constants::math::PI * randomConeR * randomConeR * collision.rho());
registry.fill(HIST("h2_centrality_rhorandomconerandomtrackdirectionwithouttwoleadingjets"), centrality, randomConePtWithoutTwoLeadJet - o2::constants::math::PI * randomConeR * randomConeR * collision.rho());
}

void processRho(soa::Filtered<soa::Join<aod::JetCollisions, aod::BkgChargedRhos>>::iterator const& collision, soa::Filtered<aod::JetTracks> const& tracks)
Expand All @@ -202,19 +203,24 @@ struct JetBackgroundAnalysisTask {
if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) {
return;
}
float centrality = checkCentFT0M ? collision.centFT0M() : collision.centFT0C();
if (centrality < centralityMin || centralityMax < centrality) {
return;
}

int nTracks = 0;
for (auto const& track : tracks) {
if (jetderiveddatautilities::selectTrack(track, trackSelection)) {
nTracks++;
}
}
registry.fill(HIST("h2_centrality_ntracks"), collision.centFT0M(), nTracks);
registry.fill(HIST("h2_centrality_ntracks"), centrality, nTracks);
registry.fill(HIST("h2_ntracks_rho"), nTracks, collision.rho());
registry.fill(HIST("h2_ntracks_rhom"), nTracks, collision.rhoM());
registry.fill(HIST("h2_centrality_rho"), collision.centFT0M(), collision.rho());
registry.fill(HIST("h2_centrality_rhom"), collision.centFT0M(), collision.rhoM());
registry.fill(HIST("h2_centrality_rho"), centrality, collision.rho());
registry.fill(HIST("h2_centrality_rhom"), centrality, collision.rhoM());
}
PROCESS_SWITCH(JetBackgroundAnalysisTask, processRho, "QA for rho-area subtracted jets", false);
PROCESS_SWITCH(JetBackgroundAnalysis, processRho, "QA for rho-area subtracted jets", false);

void processBkgFluctuationsData(soa::Filtered<soa::Join<aod::JetCollisions, aod::BkgChargedRhos>>::iterator const& collision, soa::Join<aod::ChargedJets, aod::ChargedJetConstituents> const& jets, soa::Filtered<aod::JetTracks> const& tracks)
{
Expand All @@ -224,9 +230,14 @@ struct JetBackgroundAnalysisTask {
if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) {
return;
}
bkgFluctuationsRandomCone(collision, jets, tracks);
float centrality = checkCentFT0M ? collision.centFT0M() : collision.centFT0C();
if (centrality < centralityMin || centralityMax < centrality) {
return;
}

bkgFluctuationsRandomCone(collision, jets, tracks, centrality);
}
PROCESS_SWITCH(JetBackgroundAnalysisTask, processBkgFluctuationsData, "QA for random cone estimation of background fluctuations in data", false);
PROCESS_SWITCH(JetBackgroundAnalysis, processBkgFluctuationsData, "QA for random cone estimation of background fluctuations in data", false);

void processBkgFluctuationsMCD(soa::Filtered<soa::Join<aod::JetCollisions, aod::BkgChargedRhos>>::iterator const& collision, soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents> const& jets, soa::Filtered<aod::JetTracks> const& tracks)
{
Expand All @@ -236,9 +247,17 @@ struct JetBackgroundAnalysisTask {
if (collision.trackOccupancyInTimeRange() < trackOccupancyInTimeRangeMin || trackOccupancyInTimeRangeMax < collision.trackOccupancyInTimeRange()) {
return;
}
bkgFluctuationsRandomCone(collision, jets, tracks);
float centrality = checkCentFT0M ? collision.centFT0M() : collision.centFT0C();
if (centrality < centralityMin || centralityMax < centrality) {
return;
}

bkgFluctuationsRandomCone(collision, jets, tracks, centrality);
}
PROCESS_SWITCH(JetBackgroundAnalysisTask, processBkgFluctuationsMCD, "QA for random cone estimation of background fluctuations in mcd", false);
PROCESS_SWITCH(JetBackgroundAnalysis, processBkgFluctuationsMCD, "QA for random cone estimation of background fluctuations in mcd", false);
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<JetBackgroundAnalysisTask>(cfgc, TaskName{"jet-background-analysis"})}; }
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
{
return WorkflowSpec{adaptAnalysisTask<JetBackgroundAnalysis>(cfgc)};
}
Loading