diff --git a/PWGJE/Tasks/recoilJets.cxx b/PWGJE/Tasks/recoilJets.cxx index 17ac15cc6c2..f28605ef879 100644 --- a/PWGJE/Tasks/recoilJets.cxx +++ b/PWGJE/Tasks/recoilJets.cxx @@ -39,6 +39,7 @@ #include #include +#include #include #include #include @@ -84,6 +85,9 @@ struct RecoilJets { Configurable jetR{"jetR", 0.4, "Jet cone radius"}; Configurable maxJetConstituentPt{"maxJetConstituentPt", 100., "Remove jets with constituent above this pT cut"}; + Configurable randomConeR{"randomConeR", 0.4, "Size of random cone for estimating background fluctuations"}; + Configurable randomConeJetDeltaR{"randomConeJetDeltaR", 0.0, "Min distance between high pT jet axis and random cone axis; if default, min distance is set to R_Jet + R_RC "}; + Configurable triggerMasks{"triggerMasks", "", "Relevant trigger masks: fTrackLowPt,fTrackHighPt"}; Configurable skipMBGapEvents{"skipMBGapEvents", false, "flag to choose to reject min. bias gap events; jet-level rejection " @@ -122,19 +126,20 @@ struct RecoilJets { ConfigurableAxis multFT0MThreshPartLevel{"multFT0MThreshPartLevel", {VARIABLE_WIDTH, 0.0, 0.167, 0.267, 0.4, 0.567, 0.8, 1.067, 1.4, 1.833, 2.433, 3.667, 5.1, 6.433, 20.}, "Percentiles of scaled FT0M: 100%, 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 1%, 0.1%, 0.01%"}; // Auxiliary variables - TRandom3* rand = new TRandom3(0); + std::unique_ptr randGen = std::make_unique(0); // Declare filter on collision Z vertex - Filter jCollisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; - Filter jCollisionFilterMC = nabs(aod::jmccollision::posZ) < vertexZCut; - Filter collisionFilter = nabs(aod::collision::posZ) < vertexZCut; + Filter jCollisionFilter = nabs(aod::jcollision::posZ) < vertexZCut.node(); + Filter jCollisionFilterMC = nabs(aod::jmccollision::posZ) < vertexZCut.node(); + Filter collisionFilter = nabs(aod::collision::posZ) < vertexZCut.node(); // Declare filters on accepted tracks and MC particles (settings for jet reco are provided in the jet finder wagon) - Filter trackFilter = aod::jtrack::pt > trkPtMin&& aod::jtrack::pt < trkPtMax&& nabs(aod::jtrack::eta) < trkEtaCut; - Filter partFilter = nabs(aod::jmcparticle::eta) < trkEtaCut; + Filter trackFilter = aod::jtrack::pt > trkPtMin.node() && aod::jtrack::pt < trkPtMax.node() && nabs(aod::jtrack::eta) < trkEtaCut.node(); + Filter partFilter = nabs(aod::jmcparticle::eta) < trkEtaCut.node(); // Declare filter on jets Filter jetRadiusFilter = aod::jet::r == nround(jetR.node() * 100.); + Filter jetEtaFilter = nabs(aod::jet::eta) < trkEtaCut.node() - jetR.node(); // 0.5 in our analysis HistogramRegistry spectra; @@ -509,6 +514,23 @@ struct RecoilJets { spectra.add("hScaledFT0C_Correlation_LeadTrack_AssociatTracks", Form("Leading track #it{p}_{T} #in (%.2f, %.2f); Associated track #it{p}_{T} #in (%.2f, #it{p}_{T, lead. trk})", pTLeadTrack->at(0), pTLeadTrack->at(1), pTAssociatTrackMin.value), kTH2F, {{scaledFT0C}, {160, -1.28, 5.0, "#it{#varphi} (rad)"}}); spectra.add("hScaledFT0M_Correlation_LeadTrack_AssociatTracks", Form("Leading track #it{p}_{T} #in (%.2f, %.2f); Associated track #it{p}_{T} #in (%.2f, #it{p}_{T, lead. trk})", pTLeadTrack->at(0), pTLeadTrack->at(1), pTAssociatTrackMin.value), kTH2F, {{scaledFT0M}, {160, -1.28, 5.0, "#it{#varphi} (rad)"}}); } + + if (doprocessBkgdFluctuationsRC) { + spectra.add("hScaledFT0C_deltaPtRandomCone", Form("Bkgd fluctuations RC with #it{R} = %.1f vs. EA", randomConeR.value), kTH2F, {{scaledFT0C}, {400, -40., 60., "#delta#it{p}_{T} (GeV/#it{c})"}}); + spectra.add("hScaledFT0M_deltaPtRandomCone", Form("Bkgd fluctuations RC with #it{R} = %.1f vs. EA", randomConeR.value), kTH2F, {{scaledFT0M}, {400, -40., 60., "#delta#it{p}_{T} (GeV/#it{c})"}}); + + spectra.add("hScaledFT0C_deltaPtRandomConeAvoidLeadJet", Form("Bkgd fluctuations RC with #it{R} = %.1f avoid lead jet in vic. %.1f vs. EA", randomConeR.value, randomConeJetDeltaR.value), kTH2F, {{scaledFT0C}, {400, -40., 60., "#delta#it{p}_{T} (GeV/#it{c})"}}); + spectra.add("hScaledFT0M_deltaPtRandomConeAvoidLeadJet", Form("Bkgd fluctuations RC with #it{R} = %.1f avoid lead jet in vic. %.1f vs. EA", randomConeR.value, randomConeJetDeltaR.value), kTH2F, {{scaledFT0M}, {400, -40., 60., "#delta#it{p}_{T} (GeV/#it{c})"}}); + + spectra.add("hScaledFT0C_deltaPtPerpConeAvoidLeadJet", Form("Bkgd fluctuations PC with #it{R} = %.1f avoid lead jet in vic. %.1f vs. EA", randomConeR.value, randomConeJetDeltaR.value), kTH2F, {{scaledFT0C}, {400, -40., 60., "#delta#it{p}_{T} (GeV/#it{c})"}}); + spectra.add("hScaledFT0M_deltaPtPerpConeAvoidLeadJet", Form("Bkgd fluctuations PC with #it{R} = %.1f avoid lead jet in vic. %.1f vs. EA", randomConeR.value, randomConeJetDeltaR.value), kTH2F, {{scaledFT0M}, {400, -40., 60., "#delta#it{p}_{T} (GeV/#it{c})"}}); + + spectra.add("hScaledFT0C_deltaPtRandomConeAvoidLeadAndSubleadJet", Form("Bkgd fluctuations RC with #it{R} = %.1f avoid lead, sublead jet in vic. %.1f vs. EA", randomConeR.value, randomConeJetDeltaR.value), kTH2F, {{scaledFT0C}, {400, -40., 60., "#delta#it{p}_{T} (GeV/#it{c})"}}); + spectra.add("hScaledFT0M_deltaPtRandomConeAvoidLeadAndSubleadJet", Form("Bkgd fluctuations RC with #it{R} = %.1f avoid lead, sublead jet in vic. %.1f vs. EA", randomConeR.value, randomConeJetDeltaR.value), kTH2F, {{scaledFT0M}, {400, -40., 60., "#delta#it{p}_{T} (GeV/#it{c})"}}); + + spectra.add("hScaledFT0C_deltaPtRandomConeInEventTTSig", Form("Bkgd fluctuations RC with #it{R} = %.1f in events with TT{%.0f, %.0f} in vic. %.1f vs. EA", randomConeR.value, ptTTsig->at(0), ptTTsig->at(1), randomConeJetDeltaR.value), kTH2F, {{scaledFT0C}, {400, -40., 60., "#delta#it{p}_{T} (GeV/#it{c})"}}); + spectra.add("hScaledFT0M_deltaPtRandomConeInEventTTSig", Form("Bkgd fluctuations RC with #it{R} = %.1f in events with TT{%.0f, %.0f} in vic. %.1f vs. EA", randomConeR.value, ptTTsig->at(0), ptTTsig->at(1), randomConeJetDeltaR.value), kTH2F, {{scaledFT0M}, {400, -40., 60., "#delta#it{p}_{T} (GeV/#it{c})"}}); + } } // Fill histograms with raw or MC det. level data @@ -524,7 +546,7 @@ struct RecoilJets { float scaledFT0C = getScaledFT0(collision.multFT0C(), meanFT0C); float scaledFT0M = getScaledFT0M(getScaledFT0(collision.multFT0A(), meanFT0A), scaledFT0C); - auto dice = rand->Rndm(); + auto dice = randGen->Rndm(); if (dice < fracSig) bSigEv = true; @@ -660,7 +682,7 @@ struct RecoilJets { float scaledFT0C = getScaledFT0(collision.multFT0C(), meanFT0CPartLevel); float scaledFT0M = getScaledFT0M(getScaledFT0(collision.multFT0A(), meanFT0APartLevel), scaledFT0C); - auto dice = rand->Rndm(); + auto dice = randGen->Rndm(); if (dice < fracSig) bSigEv = true; @@ -675,14 +697,14 @@ struct RecoilJets { if (!pdgParticle) continue; - float particlePt = particle.pt(); - float particlePhi = particle.phi(); - // Need charge and physical primary particles bool bParticleNeutral = (static_cast(pdgParticle->Charge()) == 0); if (bParticleNeutral || !particle.isPhysicalPrimary()) continue; + float particlePt = particle.pt(); + float particlePhi = particle.phi(); + spectra.fill(HIST("hScaledFT0CPartPtEtaPhi"), scaledFT0C, particlePt, particle.eta(), particlePhi, weight); spectra.fill(HIST("hScaledFT0MPartPtEtaPhi"), scaledFT0M, particlePt, particle.eta(), particlePhi, weight); @@ -1027,7 +1049,7 @@ struct RecoilJets { int nLeadingTracks = vPhiOfLeadingTracks.size(); if (nLeadingTracks > 0) { - auto indexLeadTrack = rand->Integer(nLeadingTracks); + auto indexLeadTrack = randGen->Integer(nLeadingTracks); double phiLeadingTrack = vPhiOfLeadingTracks[indexLeadTrack]; double pTLeadingTrack = vPtOfLeadingTracks[indexLeadTrack]; @@ -1052,6 +1074,409 @@ struct RecoilJets { } } + // Background fluctuations in MC part. level + template + void fillBkgdFluctuationsRC(JCollision const& collision, + Jets const& jets, + JTracks const& tracks, + float weight = 1.) + { + //---------------------------------------------------------- + float rho = collision.rho(); + float scaledFT0C = getScaledFT0(collision.multFT0C(), meanFT0C); + float scaledFT0M = getScaledFT0M(getScaledFT0(collision.multFT0A(), meanFT0A), scaledFT0C); + + //---------------------------------------------------------- + // Study bkgd fluctuations in events with TTsig + std::vector vCandForTT; + + //---------------------------------------------------------- + // Place absolutely random cone (anywhere in an event) + float randomConeEta = randGen->Uniform(-trkEtaCut + randomConeR, trkEtaCut - randomConeR); + float randomConePhi = randGen->Uniform(0.0, constants::math::TwoPI); + float radiusRC2 = std::pow(randomConeR, 2); + float areaRC = constants::math::PI * radiusRC2; + float randomConePt = 0.0; + + uint64_t index = 0; + for (const auto& track : tracks) { + if (skipTrack(track)) + continue; + + float dEta = std::pow(randomConeEta - track.eta(), 2); + float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - track.phi(), -constants::math::PI), 2); + + if ((dEta + dPhi) < radiusRC2) // inside RC + { + randomConePt += track.pt(); + } + + // Search for TT candidate + if (track.pt() > ptTTsig->at(0) && track.pt() < ptTTsig->at(1)) { + vCandForTT.emplace_back(index); + } + ++index; + } + spectra.fill(HIST("hScaledFT0C_deltaPtRandomCone"), scaledFT0C, randomConePt - areaRC * rho, weight); + spectra.fill(HIST("hScaledFT0M_deltaPtRandomCone"), scaledFT0M, randomConePt - areaRC * rho, weight); + + //---------------------------------------------------------- + // Avoid leading jet (JE jet reconstruction sorts jets by pT) + + // square of distance to accept RC placement in events with leading jet + float dMinR2 = std::pow(jetR + randomConeR + randomConeJetDeltaR, 2); + + // max # of attempts to find a place for RC; to avoid possibility with infinite loop in While cycle + const int maxAttempts = 15000; + + if (jets.size() > 0) // at least 1 jet + { + float leadJetEta = jets.iteratorAt(0).eta(); + float leadJetPhi = jets.iteratorAt(0).phi(); + + float dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); + float dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); + + bool isTherePlaceForRC = false; + for (int attempt = 0; attempt < maxAttempts; ++attempt) { + if ((dPhiLeadJet + dEtaLeadJet) > dMinR2) { + isTherePlaceForRC = true; + break; + } + + randomConeEta = randGen->Uniform(-trkEtaCut + randomConeR, trkEtaCut - randomConeR); + randomConePhi = randGen->Uniform(0.0, constants::math::TwoPI); + + dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); + dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); + } + + if (isTherePlaceForRC) { + randomConePt = 0.0; + for (const auto& track : tracks) { + if (skipTrack(track)) + continue; + + float dEta = std::pow(randomConeEta - track.eta(), 2); + float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - track.phi(), -constants::math::PI), 2); + + if ((dEta + dPhi) < radiusRC2) // inside RC + { + randomConePt += track.pt(); + } + } + spectra.fill(HIST("hScaledFT0C_deltaPtRandomConeAvoidLeadJet"), scaledFT0C, randomConePt - areaRC * rho, weight); + spectra.fill(HIST("hScaledFT0M_deltaPtRandomConeAvoidLeadJet"), scaledFT0M, randomConePt - areaRC * rho, weight); + } + + //---------------------------------------------------------- + // Place cone perpendicular to the leading jet (perpendicular cone) + float perpConeEta = leadJetEta; + float perpConePhi = leadJetPhi + constants::math::PI / 2.; + + float perpConePt = 0.0; + for (const auto& track : tracks) { + if (skipTrack(track)) + continue; + + float dEta = std::pow(perpConeEta - track.eta(), 2); + float dPhi = std::pow(RecoDecay::constrainAngle(perpConePhi - track.phi(), -constants::math::PI), 2); + + if ((dEta + dPhi) < radiusRC2) // inside RC + { + perpConePt += track.pt(); + } + } + spectra.fill(HIST("hScaledFT0C_deltaPtPerpConeAvoidLeadJet"), scaledFT0C, perpConePt - areaRC * rho, weight); + spectra.fill(HIST("hScaledFT0M_deltaPtPerpConeAvoidLeadJet"), scaledFT0M, perpConePt - areaRC * rho, weight); + } + + //---------------------------------------------------------- + // Avoid leading and subleading jets + if (jets.size() > 1) // at least 2 jets in an event + { + + // Leading jet + float leadJetEta = jets.iteratorAt(0).eta(); + float leadJetPhi = jets.iteratorAt(0).phi(); + float dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); + float dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); + + // Subleading jet + float subleadJetEta = jets.iteratorAt(1).eta(); + float subleadJetPhi = jets.iteratorAt(1).phi(); + float dEtaSubleadJet = std::pow(subleadJetEta - randomConeEta, 2); + float dPhiSubleadJet = std::pow(RecoDecay::constrainAngle(subleadJetPhi - randomConePhi, -constants::math::PI), 2); + + // Try to add events with TTsig + bool keepEventWithTT = false; + if (vCandForTT.size() > 0) // at least 1 TT + { + auto randIndexTrack = randGen->Integer(vCandForTT.size()); + auto objTT = tracks.iteratorAt(vCandForTT[randIndexTrack]); + + // Skip events where TT is not a part of leading or subleading jets (mutlijet event, difficult to place RC and avoid hard jets) + if (IsTrackInJet(jets.iteratorAt(0), objTT) || IsTrackInJet(jets.iteratorAt(1), objTT)) { + keepEventWithTT = true; + } + } + + //---------------------------------------------------------- + bool isTherePlaceForRC = false; + for (int attempt = 0; attempt < maxAttempts; ++attempt) { + if ((dPhiLeadJet + dEtaLeadJet) > dMinR2 && (dEtaSubleadJet + dPhiSubleadJet) > dMinR2) { + isTherePlaceForRC = true; + break; + } + + randomConeEta = randGen->Uniform(-trkEtaCut + randomConeR, trkEtaCut - randomConeR); + randomConePhi = randGen->Uniform(0.0, constants::math::TwoPI); + + dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); + dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); + + dEtaSubleadJet = std::pow(subleadJetEta - randomConeEta, 2); + dPhiSubleadJet = std::pow(RecoDecay::constrainAngle(subleadJetPhi - randomConePhi, -constants::math::PI), 2); + } + + if (isTherePlaceForRC) { + randomConePt = 0.0; + for (const auto& track : tracks) { + if (skipTrack(track)) + continue; + + float dEta = std::pow(randomConeEta - track.eta(), 2); + float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - track.phi(), -constants::math::PI), 2); + + if ((dEta + dPhi) < radiusRC2) // inside RC + { + randomConePt += track.pt(); + } + } + spectra.fill(HIST("hScaledFT0C_deltaPtRandomConeAvoidLeadAndSubleadJet"), scaledFT0C, randomConePt - areaRC * rho, weight); + spectra.fill(HIST("hScaledFT0M_deltaPtRandomConeAvoidLeadAndSubleadJet"), scaledFT0M, randomConePt - areaRC * rho, weight); + + if (keepEventWithTT) { + spectra.fill(HIST("hScaledFT0C_deltaPtRandomConeInEventTTSig"), scaledFT0C, randomConePt - areaRC * rho, weight); + spectra.fill(HIST("hScaledFT0M_deltaPtRandomConeInEventTTSig"), scaledFT0M, randomConePt - areaRC * rho, weight); + } + } + } + } + + template + void fillBkgdFluctuationsRCMCP(JCollision const& collision, + Jets const& jets, + JParticles const& particles, + float weight = 1.) + { + //---------------------------------------------------------- + float rho = collision.rho(); + float scaledFT0C = getScaledFT0(collision.multFT0C(), meanFT0C); + float scaledFT0M = getScaledFT0M(getScaledFT0(collision.multFT0A(), meanFT0A), scaledFT0C); + + //---------------------------------------------------------- + // Study bkgd fluctuations in events with TTsig + std::vector vCandForTT; + + //---------------------------------------------------------- + // Place absolutely random cone (anywhere in an event) + float randomConeEta = randGen->Uniform(-trkEtaCut + randomConeR, trkEtaCut - randomConeR); + float randomConePhi = randGen->Uniform(0.0, constants::math::TwoPI); + float radiusRC2 = std::pow(randomConeR, 2); + float areaRC = constants::math::PI * radiusRC2; + float randomConePt = 0.0; + + uint64_t index = 0; + for (const auto& particle : particles) { + auto pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (!pdgParticle) + continue; + + // Need charge and physical primary particles + bool bParticleNeutral = (static_cast(pdgParticle->Charge()) == 0); + if (bParticleNeutral || !particle.isPhysicalPrimary()) + continue; + + float dEta = std::pow(randomConeEta - particle.eta(), 2); + float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - particle.phi(), -constants::math::PI), 2); + + if ((dEta + dPhi) < radiusRC2) // inside RC + { + randomConePt += particle.pt(); + } + + // Search for TT candidate + if (particle.pt() > ptTTsig->at(0) && particle.pt() < ptTTsig->at(1)) { + vCandForTT.emplace_back(index); + } + ++index; + } + spectra.fill(HIST("hScaledFT0C_deltaPtRandomCone"), scaledFT0C, randomConePt - areaRC * rho, weight); + spectra.fill(HIST("hScaledFT0M_deltaPtRandomCone"), scaledFT0M, randomConePt - areaRC * rho, weight); + + //---------------------------------------------------------- + // Avoid leading jet (JE jet reconstruction sorts jets by pT) + + // square of distance to accept RC placement in events with leading jet + float dMinR2 = std::pow(jetR + randomConeR + randomConeJetDeltaR, 2); + + // max # of attempts to find a place for RC; to avoid possibility with infinite loop in While cycle + const int maxAttempts = 15000; + + if (jets.size() > 0) // at least 1 jet + { + float leadJetEta = jets.iteratorAt(0).eta(); + float leadJetPhi = jets.iteratorAt(0).phi(); + + float dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); + float dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); + + bool isTherePlaceForRC = false; + for (int attempt = 0; attempt < maxAttempts; ++attempt) { + if ((dPhiLeadJet + dEtaLeadJet) > dMinR2) { + isTherePlaceForRC = true; + break; + } + + randomConeEta = randGen->Uniform(-trkEtaCut + randomConeR, trkEtaCut - randomConeR); + randomConePhi = randGen->Uniform(0.0, constants::math::TwoPI); + + dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); + dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); + } + + if (isTherePlaceForRC) { + randomConePt = 0.0; + for (const auto& particle : particles) { + auto pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (!pdgParticle) + continue; + + // Need charge and physical primary particles + bool bParticleNeutral = (static_cast(pdgParticle->Charge()) == 0); + if (bParticleNeutral || !particle.isPhysicalPrimary()) + continue; + + float dEta = std::pow(randomConeEta - particle.eta(), 2); + float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - particle.phi(), -constants::math::PI), 2); + + if ((dEta + dPhi) < radiusRC2) // inside RC + { + randomConePt += particle.pt(); + } + } + spectra.fill(HIST("hScaledFT0C_deltaPtRandomConeAvoidLeadJet"), scaledFT0C, randomConePt - areaRC * rho, weight); + spectra.fill(HIST("hScaledFT0M_deltaPtRandomConeAvoidLeadJet"), scaledFT0M, randomConePt - areaRC * rho, weight); + } + + //---------------------------------------------------------- + // Place cone perpendicular to the leading jet (perpendicular cone) + float perpConeEta = leadJetEta; + float perpConePhi = leadJetPhi + constants::math::PI / 2.; + + float perpConePt = 0.0; + for (const auto& particle : particles) { + auto pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (!pdgParticle) + continue; + + // Need charge and physical primary particles + bool bParticleNeutral = (static_cast(pdgParticle->Charge()) == 0); + if (bParticleNeutral || !particle.isPhysicalPrimary()) + continue; + + float dEta = std::pow(perpConeEta - particle.eta(), 2); + float dPhi = std::pow(RecoDecay::constrainAngle(perpConePhi - particle.phi(), -constants::math::PI), 2); + + if ((dEta + dPhi) < radiusRC2) // inside RC + { + perpConePt += particle.pt(); + } + } + spectra.fill(HIST("hScaledFT0C_deltaPtPerpConeAvoidLeadJet"), scaledFT0C, perpConePt - areaRC * rho, weight); + spectra.fill(HIST("hScaledFT0M_deltaPtPerpConeAvoidLeadJet"), scaledFT0M, perpConePt - areaRC * rho, weight); + } + + //---------------------------------------------------------- + // Avoid leading and subleading jets + if (jets.size() > 1) // at least 2 jets in an event + { + + // Leading jet + float leadJetEta = jets.iteratorAt(0).eta(); + float leadJetPhi = jets.iteratorAt(0).phi(); + float dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); + float dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); + + // Subleading jet + float subleadJetEta = jets.iteratorAt(1).eta(); + float subleadJetPhi = jets.iteratorAt(1).phi(); + float dEtaSubleadJet = std::pow(subleadJetEta - randomConeEta, 2); + float dPhiSubleadJet = std::pow(RecoDecay::constrainAngle(subleadJetPhi - randomConePhi, -constants::math::PI), 2); + + // Try to add events with TTsig + bool keepEventWithTT = false; + if (vCandForTT.size() > 0) // at least 1 TT + { + auto randIndexParticle = randGen->Integer(vCandForTT.size()); + auto objTT = particles.iteratorAt(vCandForTT[randIndexParticle]); + + // Skip events where TT is not a part of leading or subleading jets (mutlijet event, difficult to place RC and avoid hard jets) + if (IsTrackInJet(jets.iteratorAt(0), objTT) || IsTrackInJet(jets.iteratorAt(1), objTT)) { + keepEventWithTT = true; + } + } + + //---------------------------------------------------------- + bool isTherePlaceForRC = false; + for (int attempt = 0; attempt < maxAttempts; ++attempt) { + if ((dPhiLeadJet + dEtaLeadJet) > dMinR2 && (dEtaSubleadJet + dPhiSubleadJet) > dMinR2) { + isTherePlaceForRC = true; + break; + } + + randomConeEta = randGen->Uniform(-trkEtaCut + randomConeR, trkEtaCut - randomConeR); + randomConePhi = randGen->Uniform(0.0, constants::math::TwoPI); + + dEtaLeadJet = std::pow(leadJetEta - randomConeEta, 2); + dPhiLeadJet = std::pow(RecoDecay::constrainAngle(leadJetPhi - randomConePhi, -constants::math::PI), 2); + + dEtaSubleadJet = std::pow(subleadJetEta - randomConeEta, 2); + dPhiSubleadJet = std::pow(RecoDecay::constrainAngle(subleadJetPhi - randomConePhi, -constants::math::PI), 2); + } + + if (isTherePlaceForRC) { + randomConePt = 0.0; + for (const auto& particle : particles) { + auto pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (!pdgParticle) + continue; + + // Need charge and physical primary particles + bool bParticleNeutral = (static_cast(pdgParticle->Charge()) == 0); + if (bParticleNeutral || !particle.isPhysicalPrimary()) + continue; + + float dEta = std::pow(randomConeEta - particle.eta(), 2); + float dPhi = std::pow(RecoDecay::constrainAngle(randomConePhi - particle.phi(), -constants::math::PI), 2); + + if ((dEta + dPhi) < radiusRC2) // inside RC + { + randomConePt += particle.pt(); + } + } + spectra.fill(HIST("hScaledFT0C_deltaPtRandomConeAvoidLeadAndSubleadJet"), scaledFT0C, randomConePt - areaRC * rho, weight); + spectra.fill(HIST("hScaledFT0M_deltaPtRandomConeAvoidLeadAndSubleadJet"), scaledFT0M, randomConePt - areaRC * rho, weight); + + if (keepEventWithTT) { + spectra.fill(HIST("hScaledFT0C_deltaPtRandomConeInEventTTSig"), scaledFT0C, randomConePt - areaRC * rho, weight); + spectra.fill(HIST("hScaledFT0M_deltaPtRandomConeInEventTTSig"), scaledFT0M, randomConePt - areaRC * rho, weight); + } + } + } + } + //------------------------------------------------------------------------------ // Process functions void processData(FilteredColl const& collision, FilteredTracks const& tracks, @@ -1250,6 +1675,40 @@ struct RecoilJets { } PROCESS_SWITCH(RecoilJets, processLeadingAndAssociatedTracksTask, "process function for correlation between leading and associated tracks", false); + void processBkgdFluctuationsRC(FilteredColl const& collision, + FilteredTracks const& tracks, + FilteredJets const& jets) + { + if (skipEvent(collision)) + return; + + fillBkgdFluctuationsRC(collision, jets, tracks); + } + PROCESS_SWITCH(RecoilJets, processBkgdFluctuationsRC, "process data to estimate bkgd fluctuations using random cone (RC)", false); + + void processBkgdFluctuationsRCMCP(FilteredCollPartLevel const& collision, + FilteredParticles const& particles, + FilteredJetsPartLevel const& jets) + { + if (skipMBGapEvent(collision)) + return; + + fillBkgdFluctuationsRCMCP(collision, jets, particles); + } + PROCESS_SWITCH(RecoilJets, processBkgdFluctuationsRCMCP, "process data to estimate bkgd fluctuations using random cone (RC) in MC Particle Level", false); + + void processBkgdFluctuationsRCMCPWeighted(FilteredCollPartLevel const& collision, + FilteredParticles const& particles, + FilteredJetsPartLevel const& jets) + { + if (skipMBGapEvent(collision) || collision.isOutlier()) + return; + + auto weight = collision.weight(); + fillBkgdFluctuationsRCMCP(collision, jets, particles, weight); + } + PROCESS_SWITCH(RecoilJets, processBkgdFluctuationsRCMCPWeighted, "process data to estimate bkgd fluctuations using random cone (RC) in MC Particle Level weighted", false); + //------------------------------------------------------------------------------ // Auxiliary functions template @@ -1280,7 +1739,7 @@ struct RecoilJets { double getPhiTT(const std::vector& vPhiOfTT) { - auto iTrig = rand->Integer(vPhiOfTT.size()); + auto iTrig = randGen->Integer(vPhiOfTT.size()); return vPhiOfTT[iTrig]; } @@ -1307,6 +1766,17 @@ struct RecoilJets { return bIsJetWithHighPtConstituent; } + template + bool IsTrackInJet(Jet const& jet, Track const& track) + { + for (auto const& constituentId : jet.tracksIds()) { + if (constituentId == track.globalIndex()) { + return true; + } + } + return false; + } + template int getBinNumberOnYaxisForGivenRun(std::shared_ptr histogram, int runNumber) {