jetSpectraEseTask.cxx

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.

/// \file jetSpectraEseTask.cxx
/// \brief jet spectra analysis framework with ESE (19/08/2024)
///
/// \author Joachim C. K. B. Hansen, Lund University

#include "PWGJE/Core/JetDerivedDataUtilities.h"
#include "PWGJE/Core/JetFindingUtilities.h"
#include "PWGJE/DataModel/Jet.h"
#include "PWGJE/DataModel/JetReducedData.h"
#include "PWGJE/DataModel/JetSubtraction.h"

#include "Common/Core/RecoDecay.h"
#include "Common/DataModel/EseTable.h"
#include "Common/DataModel/Qvectors.h"

#include <CCDB/BasicCCDBManager.h>
#include <CommonConstants/MathConstants.h>
#include <Framework/ASoA.h>
#include <Framework/AnalysisHelpers.h>
#include <Framework/AnalysisTask.h>
#include <Framework/BinningPolicy.h>
#include <Framework/Configurable.h>
#include <Framework/GroupedCombinations.h>
#include <Framework/HistogramRegistry.h>
#include <Framework/HistogramSpec.h>
#include <Framework/InitContext.h>
#include <Framework/Logger.h>
#include <Framework/O2DatabasePDGPlugin.h>
#include <Framework/OutputObjHeader.h>
#include <Framework/runDataProcessing.h>

#include <TF1.h>
#include <TH1.h>
#include <THnSparse.h>
#include <TMath.h>
#include <TRandom3.h>

#include <fmt/format.h>

#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <cstdlib>
#include <map>
#include <memory>
#include <string>
#include <string_view>
#include <vector>

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

struct JetSpectraEseTask {
  Configurable<std::string> cfgEfficiency{"cfgEfficiency", "", "CCDB path to efficiency"};
  Configurable<float> jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"};
  Configurable<float> jetR{"jetR", 0.2, "jet resolution parameter"};
  Configurable<float> randomConeR{"randomConeR", 0.4, "size of random Cone for estimating background fluctuations"};
  Configurable<float> randomConeLeadJetDeltaR{"randomConeLeadJetDeltaR", -99.0, "min distance between leading jet axis and random cone (RC) axis; if negative, min distance is set to automatic value of R_leadJet+R_RC "};
  Configurable<float> vertexZCut{"vertexZCut", 10.0, "vertex z cut"};
  // Configurable<double> leadingTrackPtCut{"leadingTrackPtCut", 5.0, "leading jet pT cut"};
  Configurable<double> jetAreaFractionMin{"jetAreaFractionMin", -99, "used to make a cut on the jet areas"};
  Configurable<bool> cfgCentVariant{"cfgCentVariant", false, "Flag for centrality variant 1"};
  Configurable<bool> cfgisPbPb{"cfgisPbPb", false, "Flag for using MC centrality in PbPb"};
  Configurable<bool> cfgbkgSubMC{"cfgbkgSubMC", true, "Flag for MC background subtraction"};
  Configurable<bool> cfgUseMCEventWeights{"cfgUseMCEventWeights", false, "Flag for using MC event weights"};
  Configurable<float> leadingConstituentPtMin{"leadingConstituentPtMin", -99.0, "minimum pT selection on jet constituent"};
  Configurable<float> leadingConstituentPtMax{"leadingConstituentPtMax", 9999.0, "maximum pT selection on jet constituent"};
  Configurable<bool> checkLeadConstituentMinPtForMcpJets{"checkLeadConstituentMinPtForMcpJets", false, "flag to choose whether particle level jets should have their lead track pt above leadingConstituentPtMin to be accepted; off by default, as leadingConstituentPtMin cut is only applied on MCD jets for the Pb-Pb analysis using pp MC anchored to Pb-Pb for the response matrix"};
  Configurable<float> pTHatMaxMCD{"pTHatMaxMCD", 999.0, "maximum fraction of hard scattering for jet acceptance in detector MC"};
  Configurable<float> pTHatMaxMCP{"pTHatMaxMCP", 999.0, "maximum fraction of hard scattering for jet acceptance in particle MC"};
  Configurable<float> pTHatExponent{"pTHatExponent", 6.0, "exponent of the event weight for the calculation of pTHat"};

  Configurable<bool> fLeadJetPtCut{"fLeadJetPtCut", false, "Flag for leading jet pT cut"};
  Configurable<float> leadJetPtMin{"leadJetPtMin", 20.0, "minimum leading jet pT cut"};

  Configurable<float> trackEtaMin{"trackEtaMin", -0.9, "minimum eta acceptance for tracks"};
  Configurable<float> trackEtaMax{"trackEtaMax", 0.9, "maximum eta acceptance for tracks"};
  Configurable<float> trackPtMin{"trackPtMin", 0.15, "minimum pT acceptance for tracks"};
  Configurable<float> trackPtMax{"trackPtMax", 200.0, "maximum pT acceptance for tracks"};
  Configurable<float> trackPtMinRhoPhi{"trackPtMinRhoPhi", 0.2, "minimum pT acceptance for tracks used in rho(phi) calculation"};
  Configurable<float> trackPtMaxRhoPhi{"trackPtMaxRhoPhi", 5.0, "maximum pT acceptance for tracks used in rho(phi) calculation"};

  Configurable<float> jetEtaMin{"jetEtaMin", -0.7, "minimum jet pseudorapidity"};
  Configurable<float> jetEtaMax{"jetEtaMax", 0.7, "maximum jet pseudorapidity"};

  Configurable<std::string> eventSelections{"eventSelections", "sel8FullPbPb", "choose event selection"};
  Configurable<std::string> trackSelections{"trackSelections", "globalTracks", "set track selections"};

  Configurable<bool> cfgEvSelOccupancy{"cfgEvSelOccupancy", true, "Flag for occupancy cut"};
  Configurable<std::vector<int>> cfgCutOccupancy{"cfgCutOccupancy", {0, 1000}, "Occupancy cut"};
  Configurable<std::vector<float>> cfgOccupancyPtCut{"cfgOccupancyPtCut", {0, 100}, "pT cut"};

  Configurable<bool> cfgrhoPhi{"cfgrhoPhi", true, "Flag for rho(phi)"};

  Configurable<int> cfgnTotalSystem{"cfgnTotalSystem", 7, "total qvector number // look in Qvector table for this number"};
  Configurable<int> cfgnCorrLevel{"cfgnCorrLevel", 3, "QVector step: 0 = no corr, 1 = rect, 2 = twist, 3 = full"};

  Configurable<std::string> cfgEPRefA{"cfgEPRefA", "FT0A", "EP reference A"};
  Configurable<std::string> cfgEPRefB{"cfgEPRefB", "TPCpos", "EP reference B"};
  Configurable<std::string> cfgEPRefC{"cfgEPRefC", "TPCneg", "EP reference C"};

  ConfigurableAxis assocTrackPt{"assocTrackPt", {VARIABLE_WIDTH, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0}, "pT of associated track"};
  ConfigurableAxis jetPtAxis{"jetPtAxis", {VARIABLE_WIDTH, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 70.0, 100.0, 150.0, 200.0}, "#it{p}_{T,jet} (GeV/#it{c})"};
  ConfigurableAxis dPhiAxis{"dPhiAxis", {60, -o2::constants::math::PI, o2::constants::math::PI}, "#Delta#phi"};
  ConfigurableAxis eseAxis{"eseAxis", {100, 0, 100}, "#it{q}_{2}"};
  Configurable<bool> cfgSelCentrality{"cfgSelCentrality", true, "Flag for centrality selection"};
  Configurable<std::vector<float>> centRange{"centRange", {0, 90}, "centrality region of interest"};
  ConfigurableAxis centAxis{"centAxis", {91, -1, 90}, "centrality"};
  ConfigurableAxis cosAxis{"cosAxis", {50, -1.0, 1.0}, "cos(2(#Psi_{2}^{A}-#Psi_{2}^{B}))"};
  ConfigurableAxis occAxis{"occAxis", {5000, 0, 25000}, "Occupancy"};
  ConfigurableAxis qvecAxis{"qvecAxis", {60, -3, 3}, "Q-vector"};
  ConfigurableAxis phiAxis{"phiAxis", {30, -1.0, 7.0}, "#phi"};
  ConfigurableAxis etaAxis{"etaAxis", {30, -1.0, 1.0}, "#eta"};
  ConfigurableAxis detaAxis{"detaAxis", {60, -1.6, 1.6}, "#Delta#eta"};
  ConfigurableAxis dphiAxis{"dphiAxis", {120, -o2::constants::math::PIHalf, 3 * o2::constants::math::PI / 2}, "#Delta#phi"};
  ConfigurableAxis trackPtAxis{"trackPtAxis", {100, 0, 100}, "#it{p}_{T} (GeV/#it{c})"};
  ConfigurableAxis fitAxis0{"fitAxis0", {100, 0, 50}, ""};
  ConfigurableAxis fitAxis13{"fitAxis13", {100, 0, 1.4}, ""};
  ConfigurableAxis fitAxis24{"fitAxis24", {100, 0, 10}, ""};
  ConfigurableAxis rhoAxis{"rhoAxis", {50, 0, 200}, "#rho"};
  ConfigurableAxis vertexZAxis{"vertexZAxis", {20, -10.0, 10.0}, "z vertex"};
  Configurable<int> numberEventsMixed{"numberEventsMixed", 5, "number of events mixed in ME process"};
  ConfigurableAxis binsCentrality{"binsCentrality", {VARIABLE_WIDTH, 0.0, 10., 30., 50, 70., 100.}, "Mixing bins - centrality"};
  ConfigurableAxis binsZVtx{"binsZVtx", {VARIABLE_WIDTH, -10.0f, -2.5f, 2.5f, 10.0f}, "Mixing bins - z-vertex"};

  HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false};

  std::vector<int> eventSelectionBits;
  int trackSelection{-1};

  static constexpr float Scaler = 100.0f;
  static constexpr float Acceptance = 0.9f;
  static constexpr float LowFT0Cut = 1e-8;

  Service<ccdb::BasicCCDBManager> ccdb;
  struct Efficiency {
    TH1D* hEff = nullptr;
    bool isLoaded = false;
  } cfg;

  Filter trackCuts = (aod::jtrack::pt >= trackPtMin && aod::jtrack::pt < trackPtMax && aod::jtrack::eta > trackEtaMin && aod::jtrack::eta < trackEtaMax);
  Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * Scaler) && nabs(aod::jet::eta) < Acceptance - jetR;
  Filter colFilter = nabs(aod::jcollision::posZ) < vertexZCut;
  Filter mcCollisionFilter = nabs(aod::jmccollision::posZ) < vertexZCut;
  using ChargedMCDJets = soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets, aod::ChargedMCDetectorLevelJetEventWeights>>;
  Preslice<ChargedMCDJets> mcdjetsPerJCollision = o2::aod::jet::collisionId;
  Preslice<aod::JetTracks> tracksPerJCollision = o2::aod::jtrack::collisionId;
  Preslice<aod::JetTracksMCD> mcdTracksPerJCollision = o2::aod::jtrack::collisionId;

  SliceCache cache;
  using BinningType = ColumnBinningPolicy<aod::jcollision::PosZ, aod::jcollision::CentFT0C>;
  BinningType corrBinning{{binsZVtx, binsCentrality}, true};
  Service<o2::framework::O2DatabasePDG> pdg;

  enum class DetID { FT0C,
                     FT0A,
                     FT0M,
                     FV0A,
                     TPCpos,
                     TPCneg,
                     TPCall };

  struct EventPlane {
    float psi2;
    float psi3;
  };

  struct EventPlaneFiller {
    bool psi;
    bool hist;
  };

  enum EventSelFlags {
    kFilteredInputEv = 1,
    kEventSel,
    kOccupancyCut,
    kCentCut,
    kEse,
    kRhoLocal,
    kLeadJetCut
  };

  static constexpr EventPlaneFiller PsiFillerEP = {true, true};
  static constexpr EventPlaneFiller PsiFillerEse = {true, false};
  static constexpr EventPlaneFiller PsiFillerFalse = {false, false};
  TRandom3* fRndm = new TRandom3(0);
  static constexpr int NumSubSmpl = 10;
  std::array<std::shared_ptr<THnSparse>, NumSubSmpl> hSameSub;

  void init(o2::framework::InitContext&)
  {
    eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>(eventSelections));
    trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast<std::string>(trackSelections));

    LOGF(info, "jetSpectraEse::init()");

    if (doprocessESEDataCharged) {
      LOGF(info, "JetSpectraEseTask::init() - ESE Data Process");
      registry.add("eventQA/hEventCounter", "event status;event status;entries", {HistType::kTH1F, {{20, 0.0, 20.0}}});
      registry.add("eventQA/hCentralityAnalyzed", ";Centrality;entries", {HistType::kTH1F, {{centAxis}}});
      registry.add("eventQA/hJetR", ";jet R;entries", {HistType::kTH1F, {{10, 0.0, 1.0}}});
      registry.add("trackQA/hRhoTrackCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
      registry.add("hJetPt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
      registry.add("hJetPt_bkgsub", "jet pT background sub;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
      registry.add("hJetEta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{etaAxis}}});
      registry.add("hJetPhi", "jet #phi;#phi_{jet};entries", {HistType::kTH1F, {{phiAxis}}});
      registry.add("eventQA/hRho", ";#rho;entries", {HistType::kTH2F, {{centAxis}, {rhoAxis}}});
      registry.add("eventQA/hRhoPhi", ";#rho;entries", {HistType::kTH2F, {{centAxis}, {rhoAxis}}});
      registry.add("hJetArea", ";area_{jet};entries", {HistType::kTH1F, {{50, 0, 5.}}});
      registry.add("hJetAreaRho", "", {HistType::kTH1F, {{50, 0, 10.}}});
      registry.add("hCentJetPtdPhiq2", "", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {dPhiAxis}, {eseAxis}}});
      registry.add("hCentJetPtdPhiq2RhoPhi", "", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {dPhiAxis}, {eseAxis}}});

      registry.add("eventQA/before/hVtxZ", ";z_{vtx} (cm);entries", {HistType::kTH1F, {{vertexZAxis}}});
      registry.add("eventQA/after/hVtxZ", ";z_{vtx} (cm);entries", {HistType::kTH1F, {{vertexZAxis}}});

      registry.get<TH1>(HIST("eventQA/hEventCounter"))->GetXaxis()->SetBinLabel(kFilteredInputEv, "Input filtered event");
      registry.get<TH1>(HIST("eventQA/hEventCounter"))->GetXaxis()->SetBinLabel(kEventSel, "Event selection");
      registry.get<TH1>(HIST("eventQA/hEventCounter"))->GetXaxis()->SetBinLabel(kOccupancyCut, "Occupancy cut");
      registry.get<TH1>(HIST("eventQA/hEventCounter"))->GetXaxis()->SetBinLabel(kCentCut, "Centrality cut");
      registry.get<TH1>(HIST("eventQA/hEventCounter"))->GetXaxis()->SetBinLabel(kEse, "ESE available");
      registry.get<TH1>(HIST("eventQA/hEventCounter"))->GetXaxis()->SetBinLabel(kRhoLocal, "rho(#phi) available");
      registry.get<TH1>(HIST("eventQA/hEventCounter"))->GetXaxis()->SetBinLabel(kLeadJetCut, "leading jet pT cut");
      registry.get<TH1>(HIST("trackQA/hRhoTrackCounter"))->GetXaxis()->SetBinLabel(1, "Input tracks");
      registry.get<TH1>(HIST("trackQA/hRhoTrackCounter"))->GetXaxis()->SetBinLabel(2, "Track selection");
      registry.get<TH1>(HIST("trackQA/hRhoTrackCounter"))->GetXaxis()->SetBinLabel(3, "Tracks for rho(phi)");

      registry.add("trackQA/hPhiPtsum", "jet sumpt;sum p_{T};entries", {HistType::kTH1F, {{40, 0., o2::constants::math::TwoPI}}});
      registry.add("eventQA/hfitPar0", "", {HistType::kTH2F, {{centAxis}, {fitAxis0}}});
      registry.add("eventQA/hfitPar1", "", {HistType::kTH2F, {{centAxis}, {fitAxis13}}});
      registry.add("eventQA/hfitPar2", "", {HistType::kTH2F, {{centAxis}, {fitAxis24}}});
      registry.add("eventQA/hfitPar3", "", {HistType::kTH2F, {{centAxis}, {fitAxis13}}});
      registry.add("eventQA/hfitPar4", "", {HistType::kTH2F, {{centAxis}, {fitAxis24}}});
      registry.add("eventQA/hPValueCentCDF", "p-value cDF vs centrality; centrality; p-value", {HistType::kTH2F, {{centAxis}, {40, 0, 1}}});
      registry.add("eventQA/hCentChi2Ndf", "Chi2 vs centrality; centrality; #tilde{#chi^{2}}", {HistType::kTH2F, {{centAxis}, {50, 0, 5}}});
      registry.add("eventQA/hCentPhi", "centrality vs #rho(#varphi); centrality;  #rho(#varphi) ", {HistType::kTH2F, {{centAxis}, {rhoAxis}}});
      registry.add("eventQA/hdPhiRhoPhi", "#varphi vs #rho(#varphi); #varphi - #Psi_{EP,2};  #rho(#varphi) ", {HistType::kTH2F, {{40, -o2::constants::math::PI, o2::constants::math::PI}, {rhoAxis}}});

      registry.add("jetQA/before/hLeadJetPt", ";Centrality;#it{p}_{T,lead jet} (GeV/#it{c});entries", {HistType::kTH2F, {{centAxis}, {jetPtAxis}}});
      registry.add("jetQA/before/hLeadJetPhi", ";Centrality;#phi_{lead jet};entries", {HistType::kTH2F, {{centAxis}, {phiAxis}}});
      registry.add("jetQA/before/hLeadJetEta", ";Centrality;#eta_{lead jet};entries", {HistType::kTH2F, {{centAxis}, {etaAxis}}});
      registry.add("jetQA/after/hLeadJetPt", ";Centrality;#it{p}_{T,lead jet} (GeV/#it{c});entries", {HistType::kTH2F, {{centAxis}, {jetPtAxis}}});
      registry.add("jetQA/after/hLeadJetPhi", ";Centrality;#phi_{lead jet};entries", {HistType::kTH2F, {{centAxis}, {phiAxis}}});
      registry.add("jetQA/after/hLeadJetEta", ";Centrality;#eta_{lead jet};entries", {HistType::kTH2F, {{centAxis}, {etaAxis}}});
      registry.add("h3CenttrPhiPsi2", ";Centrality;#phi_{track} - #Psi_{2};#it{q}_{2}", {HistType::kTH3F, {{centAxis}, {dPhiAxis}, {eseAxis}}});

      std::vector<o2::framework::AxisSpec> axes = {
        {centAxis},
        {jetPtAxis},
        {assocTrackPt},
        {detaAxis},
        {dphiAxis},
        {dPhiAxis},
        {eseAxis}};
      registry.add("thn_jethad_corr_same", "jet-had; centrality; #it{p}_{T,jet} - #rho_{local} * area_{jet} (GeV/#it{c}); #it{p}_{T,track} (GeV/#it{c}); #Delta#eta; #Delta#phi; #Delta#phi to EP; #it{q}_{2}", {HistType::kTHnSparseF, axes});
      for (int i = 0; i < NumSubSmpl; ++i) {
        std::string n = fmt::format("subsamples/thn_jethad_corr_same_subsample{}", i);
        hSameSub[i] = registry.add<THnSparse>(n, "same;...", o2::framework::HistType::kTHnSparseF, axes);
      }
      registry.add("hNtrig", "", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {dPhiAxis}, {eseAxis}}});

      registry.add("trackQA/before/hTrackPt", "", {HistType::kTH2F, {{centAxis}, {trackPtAxis}}});
      registry.add("trackQA/before/hTrackEta", "", {HistType::kTH2F, {{centAxis}, {etaAxis}}});
      registry.add("trackQA/before/hTrackPhi", "", {HistType::kTH2F, {{centAxis}, {phiAxis}}});
      registry.add("trackQA/after/hTrackPt", "", {HistType::kTH2F, {{centAxis}, {trackPtAxis}}});
      registry.add("trackQA/after/hTrackEta", "", {HistType::kTH2F, {{centAxis}, {etaAxis}}});
      registry.add("trackQA/after/hTrackPhi", "", {HistType::kTH2F, {{centAxis}, {phiAxis}}});
    }
    if (doprocessESEDataChargedMixed) {
      registry.add("eventQA/hEventCounterMixed", "event status;event status;entries", {HistType::kTH1F, {{20, 0.0, 20.0}}});
      registry.get<TH1>(HIST("eventQA/hEventCounterMixed"))->GetXaxis()->SetBinLabel(kFilteredInputEv, "Input filtered event");
      registry.get<TH1>(HIST("eventQA/hEventCounterMixed"))->GetXaxis()->SetBinLabel(kEventSel, "Event selection");
      registry.get<TH1>(HIST("eventQA/hEventCounterMixed"))->GetXaxis()->SetBinLabel(kOccupancyCut, "Occupancy cut");
      registry.get<TH1>(HIST("eventQA/hEventCounterMixed"))->GetXaxis()->SetBinLabel(kCentCut, "Centrality cut");
      registry.get<TH1>(HIST("eventQA/hEventCounterMixed"))->GetXaxis()->SetBinLabel(kEse, "ESE available");
      registry.get<TH1>(HIST("eventQA/hEventCounterMixed"))->GetXaxis()->SetBinLabel(kRhoLocal, "rho(#phi) available");
      registry.get<TH1>(HIST("eventQA/hEventCounterMixed"))->GetXaxis()->SetBinLabel(kLeadJetCut, "leading jet pT cut");
      registry.add("eventQA/before/hVtxZMixed", ";z_{vtx} (cm);entries", {HistType::kTH1F, {{vertexZAxis}}});
      registry.add("eventQA/after/hVtxZMixed", ";z_{vtx} (cm);entries", {HistType::kTH1F, {{vertexZAxis}}});
      registry.add("eventQA/before/hVtxZMixed2", ";z_{vtx} (cm);entries", {HistType::kTH1F, {{vertexZAxis}}});
      registry.add("eventQA/after/hVtxZMixed2", ";z_{vtx} (cm);entries", {HistType::kTH1F, {{vertexZAxis}}});
      registry.add("eventQA/hCentralityAnalyzedMixed", ";Centrality;entries", {HistType::kTH1F, {{centAxis}}});
      registry.add("eventQA/hCentralityAnalyzedMixed2", ";Centrality;entries", {HistType::kTH1F, {{centAxis}}});
      registry.add("eventQA/hDeltaCentMixed", ";Centrality difference;entries", {HistType::kTH1F, {{50, -100, 100}}});
      registry.add("eventQA/hDeltaVtxZMixed", ";z_{vtx} difference (cm);entries", {HistType::kTH1F, {{50, -20.0, 20.0}}});

      registry.add("trackQA/before/hTrackPtMixed", "", {HistType::kTH2F, {{centAxis}, {trackPtAxis}}});
      registry.add("trackQA/before/hTrackEtaMixed", "", {HistType::kTH2F, {{centAxis}, {etaAxis}}});
      registry.add("trackQA/before/hTrackPhiMixed", "", {HistType::kTH2F, {{centAxis}, {phiAxis}}});
      registry.add("trackQA/after/hTrackPtMixed", "", {HistType::kTH2F, {{centAxis}, {trackPtAxis}}});
      registry.add("trackQA/after/hTrackEtaMixed", "", {HistType::kTH2F, {{centAxis}, {etaAxis}}});
      registry.add("trackQA/after/hTrackPhiMixed", "", {HistType::kTH2F, {{centAxis}, {phiAxis}}});
      std::vector<o2::framework::AxisSpec> axes = {
        {centAxis},
        {jetPtAxis},
        {assocTrackPt},
        {detaAxis},
        {dphiAxis},
        {dPhiAxis},
        {eseAxis}};
      registry.add("thn_jethad_corr_mixed", "jet-had; centrality; #it{p}_{T,jet} - #rho_{local} * area_{jet} (GeV/#it{c}); #it{p}_{T,track} (GeV/#it{c}); #Delta#eta; #Delta#phi; #Delta#phi to EP; #it{q}_{2}", {HistType::kTHnSparseF, axes});
      registry.add("hNtrigMixed", "", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {dPhiAxis}, {eseAxis}}});
    }
    if (doprocessESEEPData) {
      LOGF(info, "JetSpectraEseTask::init() - Event Plane Process");
      registry.add("eventQA/hPsi2FT0C", ";Centrality; #Psi_{2}", {HistType::kTH2F, {{centAxis}, {150, -2.5, 2.5}}});
      registry.addClone("eventQA/hPsi2FT0C", "hPsi2FT0A");
      registry.addClone("eventQA/hPsi2FT0C", "hPsi2FV0A");
      registry.addClone("eventQA/hPsi2FT0C", "hPsi2TPCpos");
      registry.addClone("eventQA/hPsi2FT0C", "hPsi2TPCneg");
      registry.add("eventQA/hCosPsi2AmC", ";Centrality;cos(2(#Psi_{2}^{A}-#Psi_{2}^{B}));#it{q}_{2}", {HistType::kTH3F, {{centAxis}, {cosAxis}, {eseAxis}}});
      registry.addClone("eventQA/hCosPsi2AmC", "hCosPsi2AmB");
      registry.addClone("eventQA/hCosPsi2AmC", "hCosPsi2BmC");
      registry.add("eventQA/hQvecUncorV2", ";Centrality;Q_x;Q_y", {HistType::kTH3F, {{centAxis}, {qvecAxis}, {qvecAxis}}});
      registry.addClone("eventQA/hQvecUncorV2", "hQvecRectrV2");
      registry.addClone("eventQA/hQvecUncorV2", "hQvecTwistV2");
      registry.addClone("eventQA/hQvecUncorV2", "hQvecFinalV2");
      registry.addClone("eventQA/hPsi2FT0C", "hEPUncorV2");
      registry.addClone("eventQA/hPsi2FT0C", "hEPRectrV2");
      registry.addClone("eventQA/hPsi2FT0C", "hEPTwistV2");
    }
    if (doprocessESEBackground) {
      LOGF(info, "JetSpectraEseTask::init() - Background Process");
      registry.add("hCentRhoRandomCone", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, {800, -400.0, 400.0}}});
      registry.add("hCentRhoRandomConeRandomTrackDir", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTH2F, {{1100, 0., 110.}, {800, -400.0, 400.0}}});
      registry.add("hCentRhoRandomConewoLeadingJet", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTHnSparseF, {{centAxis}, {800, -400.0, 400.0}, {dPhiAxis}, {eseAxis}}});
      registry.add("hCentRhoRandomConeRndTrackDirwoOneLeadingJet", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTHnSparseF, {{centAxis}, {800, -400.0, 400.0}, {dPhiAxis}, {eseAxis}}});
      registry.add("hCentRhoRandomConeRndTrackDirwoTwoLeadingJet", "; centrality; #it{p}_{T,random cone} - #it{area, random cone} * #it{rho} (GeV/c);", {HistType::kTHnSparseF, {{centAxis}, {800, -400.0, 400.0}, {dPhiAxis}, {eseAxis}}});
    }
    if (doprocessMCParticleLevel) {
      LOGF(info, "JetSpectraEseTask::init() - MC Particle level");
      registry.add("mcp/hEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
      registry.add("mcp/hCentralitySel", ";centr;entries", {HistType::kTH1F, {{centAxis}}});
      // registry.add("mcp/hJetSparse", ";Centrality;#it{p}_{T,jet part}; #eta; #phi", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {etaAxis}, {phiAxis}}});
      registry.add("mcp/hJetSparse", ";Centrality;#it{p}_{T,jet part}; #eta; #phi", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {etaAxis}, {phiAxis}}});
      // registry.add("mcp/hJetPt", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
      // registry.add("mcp/hJetEta", "particle level jet #eta;#eta_{jet part};entries", {HistType::kTH1F, {{etaAxis}}});
      // registry.add("mcp/hJetPhi", "particle level jet #phi;#phi_{jet part};entries", {HistType::kTH1F, {{phiAxis}}});

      registry.get<TH1>(HIST("mcp/hEventCounter"))->GetXaxis()->SetBinLabel(1, "Input event");
      registry.get<TH1>(HIST("mcp/hEventCounter"))->GetXaxis()->SetBinLabel(2, "Collision size < 1");
      registry.get<TH1>(HIST("mcp/hEventCounter"))->GetXaxis()->SetBinLabel(3, "MCD size != 1");
      registry.get<TH1>(HIST("mcp/hEventCounter"))->GetXaxis()->SetBinLabel(4, "Occupancy cut");
      registry.get<TH1>(HIST("mcp/hEventCounter"))->GetXaxis()->SetBinLabel(5, "Vtx z cut");
      registry.get<TH1>(HIST("mcp/hEventCounter"))->GetXaxis()->SetBinLabel(6, "Event selection");
    }
    if (doprocessMCDetectorLevel) {
      LOGF(info, "JetSpectraEseTask::init() - MC Detector level");
      registry.add("mcd/hEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
      registry.add("mcd/hCentralitySel", ";centr;entries", {HistType::kTH1F, {{centAxis}}});
      // registry.add("mcd/hJetPt", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
      // registry.add("mcd/hJetSparse", ";Centrality;#it{p}_{T,jet det}; #eta; #phi", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {etaAxis}, {phiAxis}}});
      registry.add("mcd/hJetSparse", ";Centrality;#it{p}_{T,jet det}; #eta; #phi", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {etaAxis}, {phiAxis}}});
      // registry.add("mcd/hJetEta", "particle level jet #eta;#eta_{jet part};entries", {HistType::kTH1F, {{etaAxis}}});
      // registry.add("mcd/hJetPhi", "particle level jet #phi;#phi_{jet part};entries", {HistType::kTH1F, {{phiAxis}}});

      registry.get<TH1>(HIST("mcd/hEventCounter"))->GetXaxis()->SetBinLabel(1, "Input event");
      registry.get<TH1>(HIST("mcd/hEventCounter"))->GetXaxis()->SetBinLabel(1, "Event selection");
      registry.get<TH1>(HIST("mcd/hEventCounter"))->GetXaxis()->SetBinLabel(2, "Occupancy cut");
      registry.get<TH1>(HIST("mcd/hEventCounter"))->GetXaxis()->SetBinLabel(3, "Vtx z cut");
    }
    if (doprocessMCChargedMatched) {
      LOGF(info, "JetSpectraEseTask::init() - MC Charged Matched");

      registry.add("mcm/hJetSparse", ";Centrality;#it{p}_{T,jet det}; #eta; #phi", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {etaAxis}, {phiAxis}}}); /* detector level */
      // registry.add("mcm/hPartSparseMatch", ";Centrality;#it{p}_{T,jet part}; #eta; #phi", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {etaAxis}, {phiAxis}}});
      registry.addClone("mcm/hJetSparse", "mcm/hDetSparseMatch");
      registry.addClone("mcm/hJetSparse", "mcm/hPartSparseMatch");

      registry.add("mcm/hMCEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
      registry.add("mcm/hMCDMatchedEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
      registry.add("mcm/hCentralityAnalyzed", ";Centrality;entries", {HistType::kTH1F, {{centAxis}}});

      registry.add("mcm/hMatchedJetsPtDelta", "#it{p}_{T,jet part}; det - part", {HistType::kTH2F, {{jetPtAxis}, {100, -20., 20.0}}});
      registry.add("mcm/hGenMatchedJetsPtDeltadPt", "#it{p}_{T,jet part}; det - part / part", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {200, -20., 20.0}}});
      registry.add("mcm/hRecoMatchedJetsPtDeltadPt", "#it{p}_{T,jet det}; det - part / det", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {200, -20., 20.0}}});
      registry.add("mcm/hMatchedJetsEtaDelta", "#eta_{jet part}; det - part", {HistType::kTH2F, {{etaAxis}, {200, -0.8, 0.8}}});
      registry.add("mcm/hMatchedJetsPhiDelta", "#phi_{jet part}; det - part", {HistType::kTH2F, {{phiAxis}, {200, -10.0, 10.}}});

      registry.add("mcm/hRespMcDMcPMatch", ";Centrality,#it{p}_{T, jet det}; #it{p}_{T, jet part}", {HistType::kTHnSparseF, {{centAxis}, {jetPtAxis}, {jetPtAxis}}});

      registry.get<TH1>(HIST("mcm/hMCEventCounter"))->GetXaxis()->SetBinLabel(1, "Input event");
      registry.get<TH1>(HIST("mcm/hMCEventCounter"))->GetXaxis()->SetBinLabel(2, "Collision size < 1");
      registry.get<TH1>(HIST("mcm/hMCEventCounter"))->GetXaxis()->SetBinLabel(3, "Vertex cut");
      registry.get<TH1>(HIST("mcm/hMCEventCounter"))->GetXaxis()->SetBinLabel(4, "After analysis");
      registry.get<TH1>(HIST("mcm/hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(1, "Input event");
      registry.get<TH1>(HIST("mcm/hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(2, "Vertex cut");
      registry.get<TH1>(HIST("mcm/hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(3, "Event selection");
      registry.get<TH1>(HIST("mcm/hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(4, "Occupancy cut");
      registry.get<TH1>(HIST("mcm/hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(5, "Centrality cut1:cut2");
      registry.get<TH1>(HIST("mcm/hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(6, "After analysis");
    }
    if (doprocessESEOccupancy) {
      LOGF(info, "JetSpectraEseTask::init() - Occupancy QA");
      registry.add("hEventCounterOcc", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
      registry.add("hTrackPt", "track pT;#it{p}_{T,track} (GeV/#it{c});entries", {HistType::kTHnSparseF, {{centAxis}, {100, 0, 100}, {eseAxis}, {occAxis}}});
      registry.add("hTrackEta", "track #eta;#eta_{track};entries", {HistType::kTH3F, {{centAxis}, {etaAxis}, {occAxis}}});
      registry.add("hTrackPhi", "track #phi;#phi_{track};entries", {HistType::kTH3F, {{centAxis}, {phiAxis}, {occAxis}}});
      registry.add("hOccupancy", "Occupancy;Occupancy;entries", {HistType::kTH1F, {{occAxis}}});
      registry.add("hPsiOccupancy", "Occupancy;#Psi_{2};entries", {HistType::kTH3F, {{centAxis}, {150, -2.5, 2.5}, {occAxis}}});
    }
    if (doprocessMCGenTrack) {
      LOGF(info, "JetSpectraEseTask::init() - MCGen track");
      registry.add("mcgen/hTrackPtGen", "", {HistType::kTH1F, {{assocTrackPt}}});
      registry.add("mcgen/hTrackEtaGen", "", {HistType::kTH1F, {{etaAxis}}});
      registry.add("mcgen/hTrackPhiGen", "", {HistType::kTH1F, {{phiAxis}}});
    }
    if (doprocessMCRecoTrack) {
      LOGF(info, "JetSpectraEseTask::init() - MCRec track");
      registry.add("mcrec/hTrackPtReco", "", {HistType::kTH1F, {{assocTrackPt}}});
      registry.add("mcrec/hTrackEtaReco", "", {HistType::kTH1F, {{etaAxis}}});
      registry.add("mcrec/hTrackPhiReco", "", {HistType::kTH1F, {{phiAxis}}});
    }
  }

  template <typename TBC>
  void loadEfficiency(TBC const& bc)
  {
    uint64_t timestamp = bc.timestamp();
    if (cfg.isLoaded) {
      return;
    }
    if (!cfgEfficiency.value.empty()) {
      cfg.hEff = ccdb->getForTimeStamp<TH1D>(cfgEfficiency, timestamp);
      if (cfg.hEff == nullptr) {
        LOGF(fatal, "Could not load track efficiency from %s", cfgEfficiency.value.c_str());
      }
      LOGF(info, "Loaded tracking efficiency from %s (%p)", cfgEfficiency.value.c_str(), (void*)cfg.hEff);
    }
    cfg.isLoaded = true;
    return;
  }

  template <typename TTrack>
  double getEfficiency(TTrack track)
  {
    double eff{1.0};
    if (cfg.hEff)
      eff = cfg.hEff->GetBinContent(cfg.hEff->FindBin(track.pt()));
    if (eff == 0)
      return -1.;
    else
      return 1. / eff;
  }

  template <typename TCollision, typename TJets>
  float corr(const TCollision& collision, const TJets& jet)
  {
    return jet.pt() - collision.rho() * jet.area();
  }

  template <typename TCollision, typename TJets, typename TTracks>
  void jetSpectra(TCollision const& collision,
                  soa::Filtered<TJets> const& jets,
                  TTracks const& tracks)
  {

    auto centrality = cfgCentVariant ? collision.centFT0CVariant1() : collision.centFT0M();
    if (cfgSelCentrality && !isCentralitySelected(centrality))
      return;
    registry.fill(HIST("eventQA/hEventCounter"), kCentCut);

    const auto psi{procEP<PsiFillerEse>(collision)};
    const auto qPerc{collision.qPERCFT0C()};
    if (qPerc[0] < 0)
      return;
    registry.fill(HIST("eventQA/hEventCounter"), kEse);
    std::unique_ptr<TF1> rhoFit{nullptr};
    if (cfgrhoPhi) {
      rhoFit = fitRho<true>(collision, psi, tracks, jets);
      if (!rhoFit)
        return;
    }
    registry.fill(HIST("eventQA/hEventCounter"), kRhoLocal);
    if (fLeadJetPtCut) {
      if (!isAcceptedLeadingJet(collision, jets, centrality))
        return;
    }
    registry.fill(HIST("eventQA/hEventCounter"), kLeadJetCut);
    registry.fill(HIST("eventQA/after/hVtxZ"), collision.posZ());

    registry.fill(HIST("eventQA/hRho"), centrality, collision.rho());
    registry.fill(HIST("eventQA/hCentralityAnalyzed"), centrality);
    int lRndInd = fRndm->Integer(NumSubSmpl);

    auto corrL = [&](const auto& j) { return j.pt() - evalRho(rhoFit.get(), jetR, j.phi(), collision.rho()) * j.area(); };
    for (const auto& jet : jets) {
      if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax))
        continue;
      if (!isAcceptedJet<aod::JetTracks>(jet)) {
        continue;
      }
      auto vCorr = corr(collision, jet);
      auto vCorrL = cfgrhoPhi ? corrL(jet) : vCorr;
      registry.fill(HIST("hJetPt"), jet.pt());
      registry.fill(HIST("hJetPt_bkgsub"), vCorr);
      registry.fill(HIST("hJetEta"), jet.eta());
      registry.fill(HIST("hJetPhi"), jet.phi());
      registry.fill(HIST("hJetArea"), jet.area());
      registry.fill(HIST("hJetAreaRho"), jet.area() * collision.rho());
      registry.fill(HIST("eventQA/hJetR"), jet.r() / 100.0f);

      float dPhi{RecoDecay::constrainAngle(jet.phi() - psi.psi2, -o2::constants::math::PI)};
      registry.fill(HIST("hCentJetPtdPhiq2"), centrality, vCorr, dPhi, qPerc[0]);

      if (cfgrhoPhi) {
        auto rhoLocal = evalRho(rhoFit.get(), jetR, jet.phi(), collision.rho());
        registry.fill(HIST("eventQA/hRhoPhi"), centrality, rhoLocal);
        registry.fill(HIST("hCentJetPtdPhiq2RhoPhi"), centrality, vCorrL, dPhi, qPerc[0]);
        registry.fill(HIST("eventQA/hCentPhi"), centrality, rhoFit->Eval(jet.phi()));
        registry.fill(HIST("eventQA/hdPhiRhoPhi"), dPhi, rhoLocal);
      }

      registry.fill(HIST("hNtrig"), centrality, vCorrL, dPhi, qPerc[0]);
      for (const auto& track : tracks) {
        if (!jetderiveddatautilities::selectTrack(track, trackSelection))
          continue;
        double weff = getEfficiency(track);
        if (weff < 0)
          continue;
        auto deta = track.eta() - jet.eta();
        auto dphi = RecoDecay::constrainAngle(track.phi() - jet.phi(), -o2::constants::math::PIHalf);
        registry.fill(HIST("thn_jethad_corr_same"), centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0], weff);
        hSameSub[lRndInd]->Fill(centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0], weff);
      }
    }
    for (const auto& track : tracks) {
      double weff = getEfficiency(track);
      if (weff < 0)
        continue;
      registry.fill(HIST("trackQA/before/hTrackPt"), centrality, track.pt(), weff);
      registry.fill(HIST("trackQA/before/hTrackEta"), centrality, track.eta());
      registry.fill(HIST("trackQA/before/hTrackPhi"), centrality, track.phi());
      if (!jetderiveddatautilities::selectTrack(track, trackSelection))
        continue;
      registry.fill(HIST("trackQA/after/hTrackPt"), centrality, track.pt(), weff);
      registry.fill(HIST("trackQA/after/hTrackEta"), centrality, track.eta());
      registry.fill(HIST("trackQA/after/hTrackPhi"), centrality, track.phi());
      registry.fill(HIST("h3CenttrPhiPsi2"), centrality, RecoDecay::constrainAngle(track.phi() - psi.psi2, -o2::constants::math::PI), qPerc[0]);
    }
  }

  template <typename TCollisions, typename TJets, typename TTracks>
  void jetMixed(TCollisions const& collisions,
                TJets const& jets,
                TTracks const& tracks)
  {
    auto tracksTuple = std::make_tuple(jets, tracks);
    Pair<TCollisions, TJets, TTracks, BinningType> pairData{corrBinning, numberEventsMixed, -1, collisions, tracksTuple, &cache};

    for (const auto& [c1, jets1, c2, tracks2] : pairData) {
      auto bc = c2.template bc_as<aod::JBCs>();
      loadEfficiency(bc);
      auto c1Tracks = tracks.sliceBy(tracksPerJCollision, c1.globalIndex());
      registry.fill(HIST("eventQA/before/hVtxZMixed"), c1.posZ());
      registry.fill(HIST("eventQA/before/hVtxZMixed2"), c2.posZ());
      registry.fill(HIST("eventQA/hEventCounterMixed"), kFilteredInputEv);
      if (!jetderiveddatautilities::selectCollision(c1, eventSelectionBits))
        continue;
      if (!jetderiveddatautilities::selectCollision(c2, eventSelectionBits))
        continue;
      registry.fill(HIST("eventQA/hEventCounterMixed"), kEventSel);
      if (cfgEvSelOccupancy && !isOccupancyAccepted(c1))
        continue;
      if (cfgEvSelOccupancy && !isOccupancyAccepted(c2))
        continue;
      registry.fill(HIST("eventQA/hEventCounterMixed"), kOccupancyCut);

      auto centrality = cfgCentVariant ? c1.centFT0CVariant1() : c1.centFT0M();
      if (cfgSelCentrality && !isCentralitySelected(centrality))
        continue;
      auto centrality2 = cfgCentVariant ? c2.centFT0CVariant1() : c2.centFT0M();
      if (cfgSelCentrality && !isCentralitySelected(centrality2))
        continue;
      registry.fill(HIST("eventQA/hEventCounterMixed"), kCentCut);

      const auto psi{procEP<PsiFillerFalse>(c1)};
      const auto qPerc{c1.qPERCFT0C()};
      if (qPerc[0] < 0)
        continue;
      registry.fill(HIST("eventQA/hEventCounterMixed"), kEse);
      std::unique_ptr<TF1> rhoFit{nullptr};
      if (cfgrhoPhi) {
        rhoFit = fitRho<false>(c1, psi, c1Tracks, jets1);
        if (!rhoFit)
          continue;
      }
      registry.fill(HIST("eventQA/hEventCounterMixed"), kRhoLocal);
      if (fLeadJetPtCut) {
        if (!isAcceptedLeadingJet<false>(c1, jets, centrality))
          return;
      }
      registry.fill(HIST("eventQA/hEventCounterMixed"), kLeadJetCut);
      registry.fill(HIST("eventQA/after/hVtxZMixed"), c1.posZ());
      registry.fill(HIST("eventQA/after/hVtxZMixed2"), c2.posZ());
      registry.fill(HIST("eventQA/hCentralityAnalyzedMixed"), centrality);
      registry.fill(HIST("eventQA/hCentralityAnalyzedMixed2"), centrality2);
      registry.fill(HIST("eventQA/hDeltaCentMixed"), centrality - centrality2);
      registry.fill(HIST("eventQA/hDeltaVtxZMixed"), c1.posZ() - c2.posZ());

      auto corrL = [&](const auto& j) { return j.pt() - evalRho(rhoFit.get(), jetR, j.phi(), c1.rho()) * j.area(); };
      for (const auto& jet : jets1) {
        if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax))
          continue;
        if (!isAcceptedJet<aod::JetTracks>(jet)) {
          continue;
        }
        auto vCorr = corr(c1, jet);
        auto vCorrL = cfgrhoPhi ? corrL(jet) : vCorr;
        float dPhi{RecoDecay::constrainAngle(jet.phi() - psi.psi2, -o2::constants::math::PI)};

        registry.fill(HIST("hNtrigMixed"), centrality, vCorrL, dPhi, qPerc[0]);
        for (const auto& track : tracks2) {
          if (!jetderiveddatautilities::selectTrack(track, trackSelection))
            continue;
          double weff = getEfficiency(track);
          if (weff < 0)
            continue;
          auto deta = track.eta() - jet.eta();
          auto dphi = RecoDecay::constrainAngle(track.phi() - jet.phi(), -o2::constants::math::PIHalf);
          registry.fill(HIST("thn_jethad_corr_mixed"), centrality, vCorrL, track.pt(), deta, dphi, dPhi, qPerc[0], weff);
        }
      }

      for (const auto& track : tracks2) {
        double weff = getEfficiency(track);
        if (weff < 0)
          continue;
        registry.fill(HIST("trackQA/before/hTrackPtMixed"), centrality, track.pt(), weff);
        registry.fill(HIST("trackQA/before/hTrackEtaMixed"), centrality, track.eta());
        registry.fill(HIST("trackQA/before/hTrackPhiMixed"), centrality, track.phi());
        if (!jetderiveddatautilities::selectTrack(track, trackSelection))
          continue;
        registry.fill(HIST("trackQA/after/hTrackPtMixed"), centrality, track.pt(), weff);
        registry.fill(HIST("trackQA/after/hTrackEtaMixed"), centrality, track.eta());
        registry.fill(HIST("trackQA/after/hTrackPhiMixed"), centrality, track.phi());
      }
    }
  }

  void processESEDataCharged(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::Qvectors, aod::QPercentileFT0Cs>::iterator const& collision,
                             soa::Filtered<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents>> const& jets,
                             aod::JetTracks const& tracks, aod::JBCs const&)
  {
    registry.fill(HIST("eventQA/hEventCounter"), kFilteredInputEv);
    registry.fill(HIST("eventQA/before/hVtxZ"), collision.posZ());
    if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits))
      return;
    registry.fill(HIST("eventQA/hEventCounter"), kEventSel);

    if (cfgEvSelOccupancy && !isOccupancyAccepted(collision))
      return;
    registry.fill(HIST("eventQA/hEventCounter"), kOccupancyCut);

    auto bc = collision.bc_as<aod::JBCs>();
    loadEfficiency(bc);
    jetSpectra(collision, jets, tracks);
  }
  PROCESS_SWITCH(JetSpectraEseTask, processESEDataCharged, "process ese collisions", true);

  void processESEDataChargedMixed(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::Qvectors, aod::QPercentileFT0Cs> const& collisions,
                                  soa::Filtered<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents>> const& jets,
                                  aod::JetTracks const& tracks, aod::JBCs const&)
  {
    jetMixed(collisions, jets, tracks);
  }
  PROCESS_SWITCH(JetSpectraEseTask, processESEDataChargedMixed, "process ese mixed collisions", false);

  void processESEEPData(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::Qvectors, aod::QPercentileFT0Cs>::iterator const& collision,
                        soa::Filtered<aod::ChargedJets> const&,
                        aod::JetTracks const&)
  {

    if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits))
      return;

    if (cfgEvSelOccupancy && !isOccupancyAccepted(collision))
      return;

    [[maybe_unused]] const auto psi{procEP<PsiFillerEP>(collision)};
  }
  PROCESS_SWITCH(JetSpectraEseTask, processESEEPData, "process ese collisions for filling EP and EPR", false);

  void processESEBackground(soa::Filtered<soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::Qvectors, aod::QPercentileFT0Cs>>::iterator const& collision,
                            soa::Join<aod::ChargedJets, aod::ChargedJetConstituents> const& jets,
                            aod::JetTracks const& tracks)
  {
    bkgFluctuationsRandomCone(collision, jets, tracks);
  }
  PROCESS_SWITCH(JetSpectraEseTask, processESEBackground, "process random cones with event plane and ESE", false);

  void processESEOccupancy(soa::Join<aod::JetCollisions, aod::Qvectors, aod::QPercentileFT0Cs>::iterator const& collision,
                           soa::Join<aod::JetTracks, aod::JTrackPIs> const& tracks)
  {
    float count{0.5};
    registry.fill(HIST("hEventCounterOcc"), count++);
    const auto psi{procEP<PsiFillerFalse>(collision)};
    const auto qPerc{collision.qPERCFT0C()};

    auto occupancy{collision.trackOccupancyInTimeRange()};
    registry.fill(HIST("hPsiOccupancy"), collision.centFT0M(), psi.psi2, occupancy);
    registry.fill(HIST("hOccupancy"), occupancy);

    if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits))
      return;
    registry.fill(HIST("hEventCounterOcc"), count++);

    for (auto const& track : tracks) {
      if (!jetderiveddatautilities::selectTrack(track, trackSelection))
        continue;

      registry.fill(HIST("hTrackPt"), collision.centFT0M(), track.pt(), qPerc[0], occupancy);
      if (track.pt() < cfgOccupancyPtCut->at(0) || track.pt() > cfgOccupancyPtCut->at(1))
        continue;
      registry.fill(HIST("hTrackEta"), collision.centFT0M(), track.eta(), occupancy);
      registry.fill(HIST("hTrackPhi"), collision.centFT0M(), track.phi(), occupancy);
    }
  }
  PROCESS_SWITCH(JetSpectraEseTask, processESEOccupancy, "process occupancy", false);

  void processMCParticleLevel(soa::Filtered<soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>>::iterator const& mcCollision,
                              soa::SmallGroups<aod::JetCollisionsMCD> const& collisions,
                              soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetEventWeights, aod::ChargedMCParticleLevelJetConstituents>> const& jets,
                              aod::JetParticles const&)
  {
    float counter{0.5f};
    registry.fill(HIST("mcp/hEventCounter"), counter++);
    if (mcCollision.size() < 1) {
      return;
    }
    registry.fill(HIST("mcp/hEventCounter"), counter++);
    if (collisions.size() != 1) {
      return;
    }

    registry.fill(HIST("mcp/hEventCounter"), counter++);
    auto centrality{-1};
    bool fOccupancy = true;
    bool eventSel = true;
    for (const auto& col : collisions) {
      if (cfgisPbPb)
        centrality = col.centFT0M();
      if (cfgEvSelOccupancy && !isOccupancyAccepted(col))
        fOccupancy = false;
      if (!jetderiveddatautilities::selectCollision(col, eventSelectionBits))
        eventSel = false;
    }
    if (cfgEvSelOccupancy && !fOccupancy)
      return;
    registry.fill(HIST("mcp/hEventCounter"), counter++);
    if (!(std::abs(mcCollision.posZ()) < vertexZCut)) {
      return;
    }
    registry.fill(HIST("mcp/hEventCounter"), counter++);
    if (!eventSel)
      return;
    registry.fill(HIST("mcp/hEventCounter"), counter++);

    registry.fill(HIST("mcp/hCentralitySel"), centrality);
    jetLoopMCP<aod::JetParticles>(jets, centrality, mcCollision.rho());
  }
  PROCESS_SWITCH(JetSpectraEseTask, processMCParticleLevel, "jets on particle level MC", false);

  void processMCDetectorLevel(soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::BkgChargedRhos>>::iterator const& collision,
                              ChargedMCDJets const& mcdjets,
                              aod::JetTracks const&,
                              aod::JetParticles const&)
  {
    float counter{0.5f};
    registry.fill(HIST("mcd/hEventCounter"), counter++);
    if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits))
      return;
    registry.fill(HIST("mcd/hEventCounter"), counter++);
    if (cfgEvSelOccupancy && !isOccupancyAccepted(collision))
      return;
    registry.fill(HIST("mcd/hEventCounter"), counter++);

    if (!(std::abs(collision.posZ()) < vertexZCut)) {
      return;
    }
    registry.fill(HIST("mcd/hEventCounter"), counter++);

    auto centrality = cfgisPbPb ? collision.centFT0M() : -1;

    registry.fill(HIST("mcd/hCentralitySel"), centrality);
    jetLoopMCD<aod::JetTracks>(mcdjets, centrality, collision.rho());
  }
  PROCESS_SWITCH(JetSpectraEseTask, processMCDetectorLevel, "jets on detector level", false);

  using JetMCPTable = soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents, aod::ChargedMCParticleLevelJetsMatchedToChargedMCDetectorLevelJets>>;
  void processMCChargedMatched(soa::Filtered<soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>>::iterator const& mcCol,
                               soa::SmallGroups<soa::Join<aod::JetCollisionsMCD, aod::BkgChargedRhos>> const& collisions,
                               ChargedMCDJets const& mcdjets,
                               JetMCPTable const&,
                               aod::JetTracks const&,
                               aod::JetParticles const&)
  {
    float counter{0.5f};
    registry.fill(HIST("mcm/hMCEventCounter"), counter++);
    if (mcCol.size() < 1) {
      return;
    }
    if (collisions.size() != 1) {
      return;
    }
    registry.fill(HIST("mcm/hMCEventCounter"), counter++);
    if (!(std::abs(mcCol.posZ()) < vertexZCut)) {
      return;
    }
    registry.fill(HIST("mcm/hMCEventCounter"), counter++);

    for (const auto& collision : collisions) {
      float secCount{0.5f};
      registry.fill(HIST("mcm/hMCDMatchedEventCounter"), secCount++);
      if (!(std::abs(collision.posZ()) < vertexZCut)) {
        return;
      }
      registry.fill(HIST("mcm/hMCDMatchedEventCounter"), secCount++);

      if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits))
        return;
      registry.fill(HIST("mcm/hMCDMatchedEventCounter"), secCount++);

      if (cfgEvSelOccupancy && !isOccupancyAccepted(collision))
        return;
      registry.fill(HIST("mcm/hMCDMatchedEventCounter"), secCount++);

      auto centrality = cfgisPbPb ? collision.centFT0M() : -1;
      if (cfgisPbPb)
        if (centrality < centRange->at(0) || centrality > centRange->at(1))
          registry.fill(HIST("mcm/hMCDMatchedEventCounter"), secCount++);

      registry.fill(HIST("mcm/hCentralityAnalyzed"), centrality);
      matchedJetLoop<JetMCPTable, aod::JetTracks>(mcdjets.sliceBy(mcdjetsPerJCollision, collision.globalIndex()), centrality, collision.rho(), mcCol.rho());

      registry.fill(HIST("mcm/hMCDMatchedEventCounter"), secCount++);
    }
    registry.fill(HIST("mcm/hMCEventCounter"), counter++);
  }
  PROCESS_SWITCH(JetSpectraEseTask, processMCChargedMatched, "jet MC process: geometrically matched MCP and MCD for response matrix and efficiency", false);

  // process for gen/reco tracks for pT efficiency
  bool isChargedParticle(int pdgCode)
  {
    auto pdgParticle = pdg->GetParticle(pdgCode);
    return pdgParticle && pdgParticle->Charge() != 0.0;
  }

  void processMCGenTrack(soa::Filtered<soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>>::iterator const& mcCollision,
                         soa::SmallGroups<aod::JetCollisionsMCD> const& collisions,
                         aod::JetTracksMCD const&,
                         aod::JetParticles const& particles)
  {
    if (mcCollision.size() < 1) {
      return;
    }
    if (collisions.size() != 1) {
      return;
    }
    if (!(std::abs(mcCollision.posZ()) < vertexZCut)) {
      return;
    }

    auto collision = collisions.begin();
    if (!(std::abs(collision.posZ()) < vertexZCut)) {
      return;
    }
    if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
      return;
    }
    if (cfgEvSelOccupancy && !isOccupancyAccepted(collision)) {
      return;
    }

    auto centrality = cfgisPbPb ? collision.centFT0M() : -1.0f;
    if (cfgSelCentrality && cfgisPbPb && !isCentralitySelected(centrality)) {
      return;
    }

    for (const auto& particle : particles) {
      if (particle.mcCollisionId() != mcCollision.globalIndex()) {
        continue;
      }
      if (!isChargedParticle(particle.pdgCode())) {
        continue;
      }
      if (!particle.isPhysicalPrimary()) {
        continue;
      }
      if (particle.pt() < trackPtMin || particle.pt() >= trackPtMax) {
        continue;
      }
      if (particle.eta() <= trackEtaMin || particle.eta() >= trackEtaMax) {
        continue;
      }

      registry.fill(HIST("mcgen/hTrackPtGen"), particle.pt());
      registry.fill(HIST("mcgen/hTrackEtaGen"), particle.eta());
      registry.fill(HIST("mcgen/hTrackPhiGen"), particle.phi());
    }
  }
  PROCESS_SWITCH(JetSpectraEseTask, processMCGenTrack, "jet MC process: Generated track", false);
  void processMCRecoTrack(soa::Filtered<soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>>::iterator const& mcCollision,
                          soa::SmallGroups<aod::JetCollisionsMCD> const& collisions,
                          aod::JetTracksMCD const& tracks,
                          aod::JetParticles const&)
  {
    if (mcCollision.size() < 1) {
      return;
    }
    if (collisions.size() < 1) {
      return;
    }
    if (!(std::abs(mcCollision.posZ()) < vertexZCut)) {
      return;
    }

    std::vector<int> seenMcParticles;
    for (const auto& collision : collisions) {
      if (!(std::abs(collision.posZ()) < vertexZCut)) {
        continue;
      }
      if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
        continue;
      }
      if (cfgEvSelOccupancy && !isOccupancyAccepted(collision)) {
        continue;
      }

      auto centrality = cfgisPbPb ? collision.centFT0M() : -1.0f;
      if (cfgSelCentrality && cfgisPbPb && !isCentralitySelected(centrality)) {
        continue;
      }

      auto tracksInCollision = tracks.sliceBy(mcdTracksPerJCollision, collision.globalIndex());
      for (const auto& track : tracksInCollision) {
        if (!jetderiveddatautilities::selectTrack(track, trackSelection)) {
          continue;
        }
        if (!track.has_mcParticle()) {
          continue;
        }

        auto particle = track.mcParticle_as<aod::JetParticles>();
        if (!isChargedParticle(particle.pdgCode())) {
          continue;
        }
        if (!particle.isPhysicalPrimary()) {
          continue;
        }
        if (particle.pt() < trackPtMin || particle.pt() >= trackPtMax) {
          continue;
        }
        if (particle.eta() <= trackEtaMin || particle.eta() >= trackEtaMax) {
          continue;
        }
        if (std::find(seenMcParticles.begin(), seenMcParticles.end(), particle.globalIndex()) != seenMcParticles.end()) {
          continue;
        }
        seenMcParticles.push_back(particle.globalIndex());

        registry.fill(HIST("mcrec/hTrackPtReco"), track.pt());
        registry.fill(HIST("mcrec/hTrackEtaReco"), track.eta());
        registry.fill(HIST("mcrec/hTrackPhiReco"), track.phi());
      }
    }
  }
  PROCESS_SWITCH(JetSpectraEseTask, processMCRecoTrack, "jet MC process: Reconstructed track", false);

  static constexpr float InvalidValue = 999.;

  // template <bool FillAllPsi, bool FillHist, typename EPCol>
  template <EventPlaneFiller P, typename EPCol>
  EventPlane procEP(EPCol const& vec)
  {
    constexpr std::array<float, 2> AmpCut{LowFT0Cut, 0.0};
    auto computeEP = [&AmpCut](const std::vector<float>& vec, auto det, float n) { return vec[2] > AmpCut[det] ? (1.0 / n) * std::atan2(vec[1], vec[0]) : InvalidValue; };
    std::map<std::string, float> epMap;
    std::map<std::string, float> ep3Map;
    auto vec1{qVecNoESE<DetID::FT0A, P.hist>(vec)};
    epMap["FT0A"] = computeEP(vec1, 0, 2.0);
    ep3Map["FT0A"] = computeEP(vec1, 0, 3.0);
    if constexpr (P.psi) {
      auto vec2{qVecNoESE<DetID::FT0C, false>(vec)};
      epMap["FT0C"] = computeEP(vec2, 0, 2.0);
      ep3Map["FT0C"] = computeEP(vec2, 0, 3.0);
      epMap["FV0A"] = computeEP(qVecNoESE<DetID::FV0A, false>(vec), 0, 2.0);
      epMap["TPCpos"] = computeEP(qVecNoESE<DetID::TPCpos, false>(vec), 1, 2.0);
      epMap["TPCneg"] = computeEP(qVecNoESE<DetID::TPCneg, false>(vec), 1, 2.0);
      if constexpr (P.hist)
        fillEP(/*std::make_index_sequence<5>{},*/ vec, epMap);
      auto cosPsi = [](float psiX, float psiY) { return (static_cast<double>(psiX) == InvalidValue || static_cast<double>(psiY) == InvalidValue) ? InvalidValue : std::cos(2.0 * (psiX - psiY)); };
      std::array<float, 3> epCorrContainer{};
      epCorrContainer[0] = cosPsi(epMap.at(cfgEPRefA), epMap.at(cfgEPRefC));
      epCorrContainer[1] = cosPsi(epMap.at(cfgEPRefA), epMap.at(cfgEPRefB));
      epCorrContainer[2] = cosPsi(epMap.at(cfgEPRefB), epMap.at(cfgEPRefC));
      if constexpr (P.hist)
        fillEPCos(/*std::make_index_sequence<3>{},*/ vec, epCorrContainer);
    }
    EventPlane localPlane;
    localPlane.psi2 = epMap.at(cfgEPRefA);
    localPlane.psi3 = ep3Map.at(cfgEPRefA);
    return localPlane;
    // return epMap.at(cfgEPRefA);
  }
  template </*std::size_t... Idx,*/ typename collision>
  void fillEPCos(/*const std::index_sequence<Idx...>&,*/ const collision& col, const std::array<float, 3>& Corr)
  {
    // static constexpr std::string CosList[] = {"hCosPsi2AmC", "hCosPsi2AmB", "hCosPsi2BmC"};
    // (registry.fill(HIST(CosList[Idx]), col.centrality(), Corr[Idx], col.qPERCFT0C()[0]), ...);
    registry.fill(HIST("eventQA/hCosPsi2AmC"), col.centFT0M(), Corr[0], col.qPERCFT0C()[0]);
    registry.fill(HIST("eventQA/hCosPsi2AmB"), col.centFT0M(), Corr[1], col.qPERCFT0C()[0]);
    registry.fill(HIST("eventQA/hCosPsi2BmC"), col.centFT0M(), Corr[2], col.qPERCFT0C()[0]);
  }

  template </*std::size_t... Idx,*/ typename collision>
  void fillEP(/*const std::index_sequence<Idx...>&,*/ const collision& col, const std::map<std::string, float>& epMap)
  {
    // static constexpr std::string_view EpList[] = {"hPsi2FT0A", "hPsi2FV0A", "hPsi2FT0C", "hPsi2TPCpos", "hPsi2TPCneg"};
    // (registry.fill(HIST(EpList[Idx]), col.centrality(), epMap.at(std::string(RemovePrefix(EpList[Idx])))), ...);
    registry.fill(HIST("eventQA/hPsi2FT0A"), col.centFT0M(), epMap.at("FT0A"));
    registry.fill(HIST("eventQA/hPsi2FV0A"), col.centFT0M(), epMap.at("FV0A"));
    registry.fill(HIST("eventQA/hPsi2FT0C"), col.centFT0M(), epMap.at("FT0C"));
    registry.fill(HIST("eventQA/hPsi2TPCpos"), col.centFT0M(), epMap.at("TPCpos"));
    registry.fill(HIST("eventQA/hPsi2TPCneg"), col.centFT0M(), epMap.at("TPCneg"));
  }
  constexpr std::string_view RemovePrefix(std::string_view str)
  {
    constexpr std::string_view Prefix{"hPsi2"};
    return str.substr(Prefix.size());
  }

  constexpr int detIDN(const DetID id)
  {
    switch (id) {
      case DetID::FT0C:
        return 0;
      case DetID::FT0A:
        return 1;
      case DetID::FT0M:
        return 2;
      case DetID::FV0A:
        return 3;
      case DetID::TPCpos:
        return 4;
      case DetID::TPCneg:
        return 5;
      case DetID::TPCall:
        return 6;
    }
    return -1;
  }

  template <DetID id, bool fill, typename Col>
  std::vector<float> qVecNoESE(Col collision)
  {
    // const int nmode{2};
    int detId{detIDN(id)};
    int detInd{detId * 4 /*+ cfgnTotalSystem * 4 * (nmode - 2)*/};
    if constexpr (fill) {
      if (collision.qvecAmp()[detInd] > LowFT0Cut) {
        registry.fill(HIST("eventQA/hQvecUncorV2"), collision.centFT0M(), collision.qvecRe()[detInd], collision.qvecIm()[detInd]);
        registry.fill(HIST("eventQA/hQvecRectrV2"), collision.centFT0M(), collision.qvecRe()[detInd + 1], collision.qvecIm()[detInd + 1]);
        registry.fill(HIST("eventQA/hQvecTwistV2"), collision.centFT0M(), collision.qvecRe()[detInd + 2], collision.qvecIm()[detInd + 2]);
        registry.fill(HIST("eventQA/hQvecFinalV2"), collision.centFT0M(), collision.qvecRe()[detInd + 3], collision.qvecIm()[detInd + 3]);
        registry.fill(HIST("eventQA/hEPUncorV2"), collision.centFT0M(), 0.5 * std::atan2(collision.qvecIm()[detInd], collision.qvecRe()[detInd]));
        registry.fill(HIST("eventQA/hEPRectrV2"), collision.centFT0M(), 0.5 * std::atan2(collision.qvecIm()[detInd + 1], collision.qvecRe()[detInd + 1]));
        registry.fill(HIST("eventQA/hEPTwistV2"), collision.centFT0M(), 0.5 * std::atan2(collision.qvecIm()[detInd + 2], collision.qvecRe()[detInd + 2]));
      }
    }
    std::vector<float> qVec{};
    qVec.push_back(collision.qvecRe()[detInd + cfgnCorrLevel]);
    qVec.push_back(collision.qvecIm()[detInd + cfgnCorrLevel]);
    qVec.push_back(collision.qvecAmp()[detId]);
    return qVec;
  }

  template <typename col>
  bool isOccupancyAccepted(const col& collision)
  {
    auto occupancy{collision.trackOccupancyInTimeRange()};
    if (occupancy < cfgCutOccupancy->at(0) || occupancy > cfgCutOccupancy->at(1))
      return false;
    else
      return true;
  }

  template <typename Cent>
  bool isCentralitySelected(const Cent& centrality)
  {
    if (centrality < centRange->at(0) || centrality > centRange->at(1))
      return false;
    else
      return true;
  }

  template <bool fillHist, typename Col, typename TTracks, typename Jets>
  std::unique_ptr<TF1> fitRho(const Col& col, const EventPlane& ep, TTracks const& tracks, Jets const& jets)
  {
    float leadingJetPt = 0.0;
    float leadingJetEta = 0.0;
    // float leadingJetPhi = 0.0;
    for (const auto& jet : jets) {
      if (jet.pt() > leadingJetPt) {
        leadingJetPt = jet.pt();
        leadingJetEta = jet.eta();
        // leadingJetPhi = jet.phi();
      }
    }

    int nTrk{0};
    if (jets.size() > 0) {
      for (const auto& track : tracks) {
        if constexpr (fillHist)
          registry.fill(HIST("trackQA/hRhoTrackCounter"), 0.5);
        if (jetderiveddatautilities::selectTrack(track, trackSelection) && (std::fabs(track.eta() - leadingJetEta) > jetR) && track.pt() >= trackPtMinRhoPhi && track.pt() <= trackPtMaxRhoPhi) {
          nTrk++;
          if constexpr (fillHist)
            registry.fill(HIST("trackQA/hRhoTrackCounter"), 1.5);
        }
      }
    }
    if (nTrk < 1)
      return nullptr;

    auto hPhiPt = std::unique_ptr<TH1F>(new TH1F("h_ptsum_sumpt_fit", "h_ptsum_sumpt fit use", TMath::CeilNint(std::sqrt(nTrk)), 0., o2::constants::math::TwoPI));
    for (const auto& track : tracks) {
      if (jetderiveddatautilities::selectTrack(track, trackSelection) && (std::fabs(track.eta() - leadingJetEta) > jetR) && track.pt() >= trackPtMinRhoPhi && track.pt() <= trackPtMaxRhoPhi) {
        hPhiPt->Fill(track.phi(), track.pt());
        if constexpr (fillHist) {
          registry.fill(HIST("trackQA/hRhoTrackCounter"), 2.5);
          registry.fill(HIST("trackQA/hPhiPtsum"), track.phi(), track.pt());
        }
      }
    }
    auto modulationFit = std::unique_ptr<TF1>(new TF1("fit_rholoc", "[0] * (1. + 2. * ([1] * std::cos(2. * (x - [2])) + [3] * std::cos(3. * (x - [4]))))", 0, o2::constants::math::TwoPI));

    modulationFit->SetParameter(0, 1.0);
    modulationFit->SetParameter(1, 0.01);
    modulationFit->SetParameter(3, 0.01);

    modulationFit->FixParameter(2, (ep.psi2 < 0) ? RecoDecay::constrainAngle(ep.psi2) : ep.psi2);
    modulationFit->FixParameter(4, (ep.psi3 < 0) ? RecoDecay::constrainAngle(ep.psi3) : ep.psi3);

    hPhiPt->Fit(modulationFit.get(), "Q", "", 0, o2::constants::math::TwoPI);

    if constexpr (fillHist) {
      registry.fill(HIST("eventQA/hfitPar0"), col.centFT0M(), modulationFit->GetParameter(0));
      registry.fill(HIST("eventQA/hfitPar1"), col.centFT0M(), modulationFit->GetParameter(1));
      registry.fill(HIST("eventQA/hfitPar2"), col.centFT0M(), modulationFit->GetParameter(2));
      registry.fill(HIST("eventQA/hfitPar3"), col.centFT0M(), modulationFit->GetParameter(3));
      registry.fill(HIST("eventQA/hfitPar4"), col.centFT0M(), modulationFit->GetParameter(4));
    }

    if (modulationFit->GetParameter(0) <= 0)
      return nullptr;

    double chi2{0.};
    for (int i{0}; i < hPhiPt->GetXaxis()->GetNbins(); i++) {
      if (hPhiPt->GetBinContent(i + 1) <= 0.)
        continue;
      chi2 += std::pow((hPhiPt->GetBinContent(i + 1) - modulationFit->Eval(hPhiPt->GetXaxis()->GetBinCenter(1 + i))), 2) / hPhiPt->GetBinContent(i + 1);
    }

    int nDF{1};
    int numParams{2};
    nDF = static_cast<int>(modulationFit->GetXaxis()->GetNbins()) - numParams;
    if (nDF <= 0)
      return nullptr;

    auto cDF = 1. - TMath::Gamma(nDF, chi2);

    if constexpr (fillHist) {
      registry.fill(HIST("eventQA/hPValueCentCDF"), col.centFT0M(), cDF);
      registry.fill(HIST("eventQA/hCentChi2Ndf"), col.centFT0M(), chi2 / (static_cast<float>(nDF)));
    }

    return modulationFit;
  }

  template <typename phiT, typename rhoT>
  double evalRho(TF1* fit, const float& radius, phiT const& jetPhi, rhoT const& colRho)
  {
    double integralValue{fit->Integral(jetPhi - radius, jetPhi + radius)};
    double rhoLocal{colRho / (2 * radius * fit->GetParameter(0)) * integralValue};
    return rhoLocal;
  }

  template <typename TCollisions, typename TJets, typename TTracks>
  void bkgFluctuationsRandomCone(TCollisions const& collision, TJets const& jets, TTracks const& tracks)
  {
    if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
      return;
    }
    if (cfgEvSelOccupancy && !isOccupancyAccepted(collision))
      return;

    const auto psi{procEP<PsiFillerEse>(collision)};
    auto qPerc{collision.qPERCFT0C()};
    if (qPerc[0] < 0)
      return;

    TRandom3 randomNumber(0);
    float randomConeEta = randomNumber.Uniform(trackEtaMin + randomConeR, trackEtaMax - randomConeR);
    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, -o2::constants::math::PI);
        float dEta = track.eta() - randomConeEta;
        if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
          randomConePt += track.pt();
        }
      }
    }
    registry.fill(HIST("hCentRhoRandomCone"), collision.centFT0M(), randomConePt - o2::constants::math::PI * randomConeR * randomConeR * collision.rho());

    randomConePt = 0;
    for (auto const& track : tracks) {
      if (jetderiveddatautilities::selectTrack(track, trackSelection)) {
        float dPhi = RecoDecay::constrainAngle(randomNumber.Uniform(0.0, o2::constants::math::TwoPI) - randomConePhi, -o2::constants::math::PI);
        float dEta = randomNumber.Uniform(trackEtaMin, trackEtaMax) - randomConeEta;
        if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
          randomConePt += track.pt();
        }
      }
    }
    registry.fill(HIST("hCentRhoRandomConeRandomTrackDir"), collision.centFT0M(), randomConePt - o2::constants::math::PI * randomConeR * randomConeR * collision.rho());

    if (jets.size() > 0) {
      float dPhiLeadingJet = RecoDecay::constrainAngle(jets.iteratorAt(0).phi() - randomConePhi, -o2::constants::math::PI);
      float dEtaLeadingJet = jets.iteratorAt(0).eta() - randomConeEta;

      bool jetWasInCone = false;
      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, o2::constants::math::TwoPI);
        dPhiLeadingJet = RecoDecay::constrainAngle(jets.iteratorAt(0).phi() - randomConePhi, -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)) {
            float dPhi = RecoDecay::constrainAngle(track.phi() - randomConePhi, -o2::constants::math::PI);
            float dEta = track.eta() - randomConeEta;
            if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
              randomConePt += track.pt();
            }
          }
        }
      }
    }
    auto rho = collision.rho();
    std::unique_ptr<TF1> rhoFit{nullptr};
    if (cfgrhoPhi) {
      rhoFit = fitRho<false>(collision, psi, tracks, jets);
      if (!rhoFit)
        return;
      rho = evalRho(rhoFit.get(), randomConeR, randomConePhi, rho);
    }
    float dPhiRC{RecoDecay::constrainAngle(randomConePhi - psi.psi2, -o2::constants::math::PI)};
    registry.fill(HIST("hCentRhoRandomConewoLeadingJet"), collision.centFT0M(), randomConePt - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC, qPerc[0]);

    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, o2::constants::math::TwoPI) - randomConePhi, -o2::constants::math::PI);
          float dEta = randomNumber.Uniform(trackEtaMin, trackEtaMax) - randomConeEta;
          if (std::sqrt(dEta * dEta + dPhi * dPhi) < randomConeR) {
            if (!isTrackInJet(track, jets.iteratorAt(0))) {
              randomConePtWithoutOneLeadJet += track.pt();
              if (!isTrackInJet(track, jets.iteratorAt(1))) {
                randomConePtWithoutTwoLeadJet += track.pt();
              }
            }
          }
        }
      }
    }
    registry.fill(HIST("hCentRhoRandomConeRndTrackDirwoOneLeadingJet"), collision.centFT0M(), randomConePtWithoutOneLeadJet - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC, qPerc[0]);
    registry.fill(HIST("hCentRhoRandomConeRndTrackDirwoTwoLeadingJet"), collision.centFT0M(), randomConePtWithoutTwoLeadJet - o2::constants::math::PI * randomConeR * randomConeR * rho, dPhiRC, qPerc[0]);
  }
  template <typename TTracks, typename TJets>
  bool isTrackInJet(TTracks const& track, TJets const& jet)
  {
    for (auto const& constituentId : jet.tracksIds()) {
      if (constituentId == track.globalIndex()) {
        return true;
      }
    }
    return false;
  }

  // static constexpr std::string_view LevelJets[] = {"mcd/", "mcp/"};
  // enum JetType { MCP = 0,
  //                MCD = 1
  //               };
  // template <JetType jetLvl, typename Jets>
  template <typename JTracks, typename Jets>
  void jetLoopMCD(const Jets& jets, const float& centrality, const float& rho)
  {
    float weight = 1.0;
    for (const auto& jet : jets) {
      if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
        continue;
      }
      if (!isAcceptedJet<JTracks>(jet)) {
        continue;
      }
      auto pt = jet.pt();
      if (cfgbkgSubMC) {
        pt = jet.pt() - (rho * jet.area());
      }
      if (cfgUseMCEventWeights) {
        weight = jet.eventWeight();
      }
      registry.fill(/*HIST(LevelJets[jetLvl]) +*/ HIST("mcd/hJetSparse"), centrality, pt, jet.eta(), jet.phi(), weight); /* detector level mcm*/
    }
  }

  template <typename JTracks, typename Jets>
  void jetLoopMCP(const Jets& jets, const float& centrality, const float& rho)
  {
    bool mcLevelIsParticleLevel = true;
    float weight = 1.0;
    for (const auto& jet : jets) {
      if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
        continue;
      }
      if (!isAcceptedJet<JTracks>(jet, mcLevelIsParticleLevel)) {
        continue;
      }
      auto pt = jet.pt();
      if (cfgbkgSubMC) {
        pt = jet.pt() - (rho * jet.area());
      }
      if (cfgUseMCEventWeights) {
        weight = jet.eventWeight();
      }
      registry.fill(/*HIST(LevelJets[jetLvl]) +*/ HIST("mcp/hJetSparse"), centrality, pt, jet.eta(), jet.phi(), weight); /* detector level mcm*/
    }
  }

  template <typename MCPTab, typename JTracks, typename Jets>
  void matchedJetLoop(const Jets& jets, const float& centrality, const float& rho, const float& rho2)
  {
    float weight = 1.0;
    for (const auto& jet : jets) {
      if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
        continue;
      }
      if (!isAcceptedJet<JTracks>(jet)) {
        continue;
      }
      float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent));
      if (jet.pt() > pTHatMaxMCD * pTHat) {
        return;
      }

      auto pt = jet.pt();
      if (cfgbkgSubMC) {
        pt = jet.pt() - (rho * jet.area());
      }
      if (cfgUseMCEventWeights) {
        weight = jet.eventWeight();
      }
      registry.fill(HIST("mcm/hJetSparse"), centrality, pt, jet.eta(), jet.phi(), weight); /* detector level mcm*/

      if (jet.has_matchedJetGeo()) {
        registry.fill(HIST("mcm/hDetSparseMatch"), centrality, pt, jet.eta(), jet.phi(), weight);
        for (const auto& matchedJet : jet.template matchedJetGeo_as<MCPTab>()) {
          if (matchedJet.pt() > pTHatMaxMCD * pTHat)
            continue;
          auto matchedpt = matchedJet.pt();
          if (cfgbkgSubMC) {
            matchedpt = matchedJet.pt() - (rho2 * matchedJet.area());
          }
          registry.fill(HIST("mcm/hPartSparseMatch"), centrality, matchedpt, matchedJet.eta(), matchedJet.phi(), weight);
          registry.fill(HIST("mcm/hMatchedJetsPtDelta"), matchedJet.pt(), jet.pt() - matchedJet.pt(), weight);
          registry.fill(HIST("mcm/hMatchedJetsPhiDelta"), matchedJet.phi(), jet.phi() - matchedJet.phi(), weight);
          registry.fill(HIST("mcm/hMatchedJetsEtaDelta"), matchedJet.eta(), jet.eta() - matchedJet.eta(), weight);
          registry.fill(HIST("mcm/hGenMatchedJetsPtDeltadPt"), centrality, matchedpt, (pt - matchedpt) / matchedpt, weight);
          registry.fill(HIST("mcm/hRecoMatchedJetsPtDeltadPt"), centrality, pt, (pt - matchedpt) / pt, weight);
          registry.fill(HIST("mcm/hRespMcDMcPMatch"), centrality, pt, matchedpt, weight);
        }
      }
    }
  }
  // leadJet(collision, jets)
  template <bool fill = true, typename TCollisions, typename TJets>
  bool isAcceptedLeadingJet(TCollisions const& collision, TJets const& jets, const float& centrality)
  {
    float leadJetPt = 0.0;
    float leadJetPhi = -999;
    float leadJetEta = -999;
    bool hasLeadingJet = false;
    for (const auto& jet : jets) {
      if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax))
        continue;
      if (!isAcceptedJet<aod::JetTracks>(jet)) {
        continue;
      }
      auto jetPtBkg = jet.pt() - (collision.rho() * jet.area());

      if (!hasLeadingJet) {
        leadJetPt = jetPtBkg;
        leadJetPhi = jet.phi();
        leadJetEta = jet.eta();
        hasLeadingJet = true;
      }
      if (jetPtBkg > leadJetPt) {
        leadJetPt = jetPtBkg;
        leadJetPhi = jet.phi();
        leadJetEta = jet.eta();
      }
    }
    if (!hasLeadingJet) {
      return false;
    }
    if constexpr (fill) {
      registry.fill(HIST("jetQA/before/hLeadJetPt"), centrality, leadJetPt);
      registry.fill(HIST("jetQA/before/hLeadJetPhi"), centrality, leadJetPhi);
      registry.fill(HIST("jetQA/before/hLeadJetEta"), centrality, leadJetEta);
    }
    if (leadJetPt < leadJetPtMin) {
      return false;
    }
    if constexpr (fill) {
      registry.fill(HIST("jetQA/after/hLeadJetPt"), centrality, leadJetPt);
      registry.fill(HIST("jetQA/after/hLeadJetPhi"), centrality, leadJetPhi);
      registry.fill(HIST("jetQA/after/hLeadJetEta"), centrality, leadJetEta);
    }
    return true;
  }

  template <typename TTracks, typename TJets>
  bool isAcceptedJet(TJets const& jet, bool mcLevelIsParticleLevel = false)
  {
    float minCutoff = -98.0;
    float maxCutoff = 9998.0;
    if (jetAreaFractionMin > minCutoff) {
      if (jet.area() < jetAreaFractionMin * o2::constants::math::PI * (jet.r() / 100.0) * (jet.r() / 100.0)) {
        return false;
      }
    }
    bool checkConstituentPt = true;
    bool checkConstituentMinPt = (leadingConstituentPtMin > minCutoff);
    bool checkConstituentMaxPt = (leadingConstituentPtMax < maxCutoff);
    if (!checkConstituentMinPt && !checkConstituentMaxPt) {
      checkConstituentPt = false;
    }

    if (checkConstituentPt) {
      bool isMinLeadingConstituent = !checkConstituentMinPt;
      bool isMaxLeadingConstituent = true;

      for (const auto& constituent : jet.template tracks_as<TTracks>()) {
        double pt = constituent.pt();

        if ((!checkLeadConstituentMinPtForMcpJets && mcLevelIsParticleLevel) || (checkConstituentMinPt && pt >= leadingConstituentPtMin)) { // if the jet is mcp level and checkLeadConstituentMinPtForMcpJets is true, then the pt of the leading track of that jet does not need to be below the defined leadingConstituentPtMin cut
          isMinLeadingConstituent = true;
        }
        if (checkConstituentMaxPt && pt > leadingConstituentPtMax) {
          isMaxLeadingConstituent = false;
        }
      }
      return isMinLeadingConstituent && isMaxLeadingConstituent;
    }
    return true;
  }
};
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<JetSpectraEseTask>(cfgc)}; }