Skip to content

Commit fc9e3de

Browse files
joachimckhjoachimckh
authored
[PWGJE] mc dev, fix bugs (AliceO2Group#10583)
1 parent 4054efb commit fc9e3de

File tree

1 file changed

+153
-69
lines changed

1 file changed

+153
-69
lines changed

‎PWGJE/Tasks/jetSpectraEseTask.cxx‎

Lines changed: 153 additions & 69 deletions
Original file line numberDiff line numberDiff line change
@@ -68,17 +68,23 @@ struct JetSpectraEseTask {
6868
Configurable<float> vertexZCut{"vertexZCut", 10.0, "vertex z cut"};
6969
Configurable<std::vector<float>> centRange{"centRange", {0, 90}, "centrality region of interest"};
7070
Configurable<bool> cfgSelCentrality{"cfgSelCentrality", true, "Flag for centrality selection"};
71-
Configurable<double> leadingTrackPtCut{"leadingTrackPtCut", 5.0, "leading jet pT cut"};
72-
Configurable<double> jetAreaFractionMin{"jetAreaFractionMin", 0.56, "used to make a cut on the jet areas"};
73-
Configurable<bool> fjetAreaCut{"fjetAreaCut", true, "Flag for jet area cut"};
71+
// Configurable<double> leadingTrackPtCut{"leadingTrackPtCut", 5.0, "leading jet pT cut"};
72+
Configurable<double> jetAreaFractionMin{"jetAreaFractionMin", -99, "used to make a cut on the jet areas"};
7473
Configurable<bool> cfgCentVariant{"cfgCentVariant", false, "Flag for centrality variant 1"};
7574
Configurable<bool> cfgisPbPb{"cfgisPbPb", false, "Flag for using MC centrality in PbPb"};
7675
Configurable<bool> cfgbkgSubMC{"cfgbkgSubMC", true, "Flag for MC background subtraction"};
76+
Configurable<bool> cfgUseMCEventWeights{"cfgUseMCEventWeights", false, "Flag for using MC event weights"};
77+
Configurable<float> leadingConstituentPtMin{"leadingConstituentPtMin", -99.0, "minimum pT selection on jet constituent"};
78+
Configurable<float> leadingConstituentPtMax{"leadingConstituentPtMax", 9999.0, "maximum pT selection on jet constituent"};
79+
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"};
80+
Configurable<float> pTHatMaxMCD{"pTHatMaxMCD", 999.0, "maximum fraction of hard scattering for jet acceptance in detector MC"};
81+
Configurable<float> pTHatMaxMCP{"pTHatMaxMCP", 999.0, "maximum fraction of hard scattering for jet acceptance in particle MC"};
82+
Configurable<float> pTHatExponent{"pTHatExponent", 6.0, "exponent of the event weight for the calculation of pTHat"};
7783

7884
Configurable<float> trackEtaMin{"trackEtaMin", -0.9, "minimum eta acceptance for tracks"};
7985
Configurable<float> trackEtaMax{"trackEtaMax", 0.9, "maximum eta acceptance for tracks"};
8086
Configurable<float> trackPtMin{"trackPtMin", 0.15, "minimum pT acceptance for tracks"};
81-
Configurable<float> trackPtMax{"trackPtMax", 100.0, "maximum pT acceptance for tracks"};
87+
Configurable<float> trackPtMax{"trackPtMax", 200.0, "maximum pT acceptance for tracks"};
8288

8389
Configurable<float> jetEtaMin{"jetEtaMin", -0.7, "minimum jet pseudorapidity"};
8490
Configurable<float> jetEtaMax{"jetEtaMax", 0.7, "maximum jet pseudorapidity"};
@@ -118,10 +124,11 @@ struct JetSpectraEseTask {
118124
std::vector<int> eventSelectionBits;
119125
int trackSelection{-1};
120126

127+
Filter trackCuts = (aod::jtrack::pt >= trackPtMin && aod::jtrack::pt < trackPtMax && aod::jtrack::eta > trackEtaMin && aod::jtrack::eta < trackEtaMax);
121128
Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f) && nabs(aod::jet::eta) < 0.9f - jetR;
122129
Filter colFilter = nabs(aod::jcollision::posZ) < vertexZCut;
123130
Filter mcCollisionFilter = nabs(aod::jmccollision::posZ) < vertexZCut;
124-
using ChargedMCDJets = soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>>;
131+
using ChargedMCDJets = soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets, aod::ChargedMCDetectorLevelJetEventWeights>>;
125132
Preslice<ChargedMCDJets> mcdjetsPerJCollision = o2::aod::jet::collisionId;
126133

127134
enum class DetID { FT0C,
@@ -247,7 +254,9 @@ struct JetSpectraEseTask {
247254
LOGF(info, "JetSpectraEseTask::init() - MC Charged Matched");
248255

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

252261
registry.add("mcm/hMCEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
253262
registry.add("mcm/hMCDMatchedEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
@@ -284,7 +293,7 @@ struct JetSpectraEseTask {
284293
}
285294

286295
void processESEDataCharged(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::Qvectors, aod::QPercentileFT0Cs>::iterator const& collision,
287-
soa::Filtered<aod::ChargedJets> const& jets,
296+
soa::Filtered<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents>> const& jets,
288297
aod::JetTracks const& tracks)
289298
{
290299
float counter{0.5f};
@@ -307,10 +316,6 @@ struct JetSpectraEseTask {
307316
if (qPerc[0] < 0)
308317
return;
309318
registry.fill(HIST("hEventCounter"), counter++);
310-
311-
if (!isAcceptedLeadTrack(tracks))
312-
return;
313-
314319
std::unique_ptr<TF1> rhoFit{nullptr};
315320
if (cfgrhoPhi) {
316321
rhoFit = fitRho<true>(collision, psi, tracks, jets);
@@ -322,10 +327,11 @@ struct JetSpectraEseTask {
322327
registry.fill(HIST("hRho"), centrality, collision.rho());
323328
registry.fill(HIST("hCentralityAnalyzed"), centrality);
324329
for (auto const& jet : jets) {
325-
if (fjetAreaCut && !isJetAreaAccepted(jet))
326-
continue;
327330
if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax))
328331
continue;
332+
if (!isAcceptedJet<aod::JetTracks>(jet)) {
333+
continue;
334+
}
329335
registry.fill(HIST("hJetPt"), jet.pt());
330336
registry.fill(HIST("hJetPt_bkgsub"), jet.pt() - collision.rho() * jet.area());
331337
registry.fill(HIST("hJetEta"), jet.eta());
@@ -405,7 +411,8 @@ struct JetSpectraEseTask {
405411

406412
void processMCParticleLevel(soa::Filtered<soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>>::iterator const& mcCollision,
407413
soa::SmallGroups<aod::JetCollisionsMCD> const& collisions,
408-
soa::Filtered<aod::ChargedMCParticleLevelJets> const& jets)
414+
soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetEventWeights, aod::ChargedMCParticleLevelJetConstituents>> const& jets,
415+
aod::JetParticles const&)
409416
{
410417
float counter{0.5f};
411418
registry.fill(HIST("mcp/hEventCounter"), counter++);
@@ -420,24 +427,31 @@ struct JetSpectraEseTask {
420427
registry.fill(HIST("mcp/hEventCounter"), counter++);
421428
auto centrality{-1};
422429
bool fOccupancy = true;
430+
bool eventSel = true;
423431
for (const auto& col : collisions) {
424432
if (cfgisPbPb)
425433
centrality = col.centrality();
426434
if (cfgEvSelOccupancy && !isOccupancyAccepted(col))
427435
fOccupancy = false;
436+
if (!jetderiveddatautilities::selectCollision(col, eventSelectionBits))
437+
eventSel = false;
428438
}
429439
if (cfgEvSelOccupancy && !fOccupancy)
430440
return;
441+
if (!(std::abs(mcCollision.posZ()) < vertexZCut)) {
442+
return;
443+
}
444+
if (!eventSel)
445+
return;
431446

432447
registry.fill(HIST("mcp/hEventCounter"), counter++);
433448

434449
registry.fill(HIST("mcp/hCentralitySel"), centrality);
435-
436-
jetLoop<JetType::MCP>(jets, centrality, mcCollision.rho());
450+
jetLoopMCP<aod::JetParticles>(jets, centrality, mcCollision.rho());
437451
}
438452
PROCESS_SWITCH(JetSpectraEseTask, processMCParticleLevel, "jets on particle level MC", false);
439453

440-
void processMCDetectorLevel(soa::Join<aod::JetCollisionsMCD, aod::BkgChargedRhos>::iterator const& collision,
454+
void processMCDetectorLevel(soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::BkgChargedRhos>>::iterator const& collision,
441455
ChargedMCDJets const& mcdjets,
442456
aod::JetTracks const&,
443457
aod::JetParticles const&)
@@ -451,11 +465,14 @@ struct JetSpectraEseTask {
451465
return;
452466
registry.fill(HIST("mcd/hEventCounter"), counter++);
453467

468+
if (!(std::abs(collision.posZ()) < vertexZCut)) {
469+
return;
470+
}
471+
454472
auto centrality = cfgisPbPb ? collision.centrality() : -1;
455473

456474
registry.fill(HIST("mcd/hCentralitySel"), centrality);
457-
458-
jetLoop<JetType::MCD>(mcdjets, centrality, collision.rho());
475+
jetLoopMCD<aod::JetTracks>(mcdjets, centrality, collision.rho());
459476
}
460477
PROCESS_SWITCH(JetSpectraEseTask, processMCDetectorLevel, "jets on detector level", false);
461478

@@ -472,6 +489,9 @@ struct JetSpectraEseTask {
472489
if (mcCol.size() < 1) {
473490
return;
474491
}
492+
if (collisions.size() != 1) {
493+
return;
494+
}
475495
registry.fill(HIST("mcm/hMCEventCounter"), counter++);
476496
if (!(std::abs(mcCol.posZ()) < vertexZCut)) {
477497
return;
@@ -500,44 +520,14 @@ struct JetSpectraEseTask {
500520
registry.fill(HIST("mcm/hMCDMatchedEventCounter"), secCount++);
501521

502522
registry.fill(HIST("mcm/hCentralityAnalyzed"), centrality);
503-
504-
jetLoop<JetType::MCM>(
505-
mcdjets.sliceBy(mcdjetsPerJCollision, collision.globalIndex()),
506-
centrality,
507-
collision.rho(),
508-
[](const auto& jet) { return jet.template matchedJetGeo_as<JetMCPTable>(); },
509-
mcCol.rho());
523+
matchedJetLoop<JetMCPTable, aod::JetTracks>(mcdjets.sliceBy(mcdjetsPerJCollision, collision.globalIndex()), centrality, collision.rho(), mcCol.rho());
510524

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

517-
template <typename T>
518-
bool isAcceptedLeadTrack(T const& tracks)
519-
{
520-
double leadingTrackPt = 0.0;
521-
for (const auto& track : tracks) {
522-
if (track.pt() > leadingTrackPtCut) {
523-
if (track.pt() > leadingTrackPt) {
524-
leadingTrackPt = track.pt();
525-
}
526-
}
527-
}
528-
if (leadingTrackPt == 0.0)
529-
return false;
530-
else
531-
return true;
532-
}
533-
template <typename Jet>
534-
bool isJetAreaAccepted(const Jet& jet)
535-
{
536-
if (jet.area() < jetAreaFractionMin * o2::constants::math::PI * (jet.r() / 100.0) * (jet.r() / 100.0))
537-
return false;
538-
else
539-
return true;
540-
}
541531
// template <bool FillAllPsi, bool FillHist, typename EPCol>
542532
template <EventPlaneFiller P, typename EPCol>
543533
EventPlane procEP(EPCol const& vec)
@@ -663,8 +653,8 @@ struct JetSpectraEseTask {
663653
return true;
664654
}
665655

666-
template <bool fillHist, typename C, typename T, typename J>
667-
std::unique_ptr<TF1> fitRho(const C& col, const EventPlane& ep, T const& tracks, J const& jets)
656+
template <bool fillHist, typename Col, typename TTracks, typename Jets>
657+
std::unique_ptr<TF1> fitRho(const Col& col, const EventPlane& ep, TTracks const& tracks, Jets const& jets)
668658
{
669659
float leadingJetPt = 0.0;
670660
float leadingJetEta = 0.0;
@@ -864,40 +854,134 @@ struct JetSpectraEseTask {
864854
return false;
865855
}
866856

867-
static constexpr std::string_view LevelJets[] = {"mcd/", "mcp/", "mcm/"};
868-
enum JetType { MCP = 0,
869-
MCD = 1,
870-
MCM = 2 };
871-
template <JetType jetLvl, typename Jets, typename matchJet = std::nullptr_t>
872-
void jetLoop(const Jets& jets, const float& centrality, const float& rho, matchJet matchjet = nullptr, const float& rho2 = 0)
857+
// static constexpr std::string_view LevelJets[] = {"mcd/", "mcp/"};
858+
// enum JetType { MCP = 0,
859+
// MCD = 1
860+
// };
861+
// template <JetType jetLvl, typename Jets>
862+
template <typename JTracks, typename Jets>
863+
void jetLoopMCD(const Jets& jets, const float& centrality, const float& rho)
864+
{
865+
float weight = 1.0;
866+
for (const auto& jet : jets) {
867+
if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
868+
continue;
869+
}
870+
if (!isAcceptedJet<JTracks>(jet)) {
871+
continue;
872+
}
873+
auto pt = jet.pt();
874+
if (cfgbkgSubMC) {
875+
pt = jet.pt() - (rho * jet.area());
876+
}
877+
if (cfgUseMCEventWeights) {
878+
weight = jet.eventWeight();
879+
}
880+
registry.fill(/*HIST(LevelJets[jetLvl]) +*/ HIST("mcd/hJetSparse"), centrality, pt, jet.eta(), jet.phi(), weight); /* detector level mcm*/
881+
}
882+
}
883+
884+
template <typename JTracks, typename Jets>
885+
void jetLoopMCP(const Jets& jets, const float& centrality, const float& rho)
873886
{
887+
bool mcLevelIsParticleLevel = true;
888+
float weight = 1.0;
889+
for (const auto& jet : jets) {
890+
if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
891+
continue;
892+
}
893+
if (!isAcceptedJet<JTracks>(jet, mcLevelIsParticleLevel)) {
894+
continue;
895+
}
896+
auto pt = jet.pt();
897+
if (cfgbkgSubMC) {
898+
pt = jet.pt() - (rho * jet.area());
899+
}
900+
if (cfgUseMCEventWeights) {
901+
weight = jet.eventWeight();
902+
}
903+
registry.fill(/*HIST(LevelJets[jetLvl]) +*/ HIST("mcp/hJetSparse"), centrality, pt, jet.eta(), jet.phi(), weight); /* detector level mcm*/
904+
}
905+
}
874906

907+
template <typename MCPTab, typename JTracks, typename Jets>
908+
void matchedJetLoop(const Jets& jets, const float& centrality, const float& rho, const float& rho2)
909+
{
910+
float weight = 1.0;
875911
for (const auto& jet : jets) {
912+
if (!jetfindingutilities::isInEtaAcceptance(jet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax)) {
913+
continue;
914+
}
915+
if (!isAcceptedJet<JTracks>(jet)) {
916+
continue;
917+
}
918+
float pTHat = 10. / (std::pow(weight, 1.0 / pTHatExponent));
919+
if (jet.pt() > pTHatMaxMCD * pTHat) {
920+
return;
921+
}
922+
876923
auto pt = jet.pt();
877924
if (cfgbkgSubMC) {
878925
pt = jet.pt() - (rho * jet.area());
879926
}
880-
// registry.fill(HIST(levelJets[jetLvl]) + HIST("hJetSparse"), centrality, jet.pt(), jet.eta(), jet.phi());
881-
registry.fill(HIST(LevelJets[jetLvl]) + HIST("hJetSparse"), centrality, pt, jet.eta(), jet.phi()); /* detector level mcm*/
927+
if (cfgUseMCEventWeights) {
928+
weight = jet.eventWeight();
929+
}
930+
registry.fill(HIST("mcm/hJetSparse"), centrality, pt, jet.eta(), jet.phi(), weight); /* detector level mcm*/
882931

883-
if constexpr (jetLvl == MCM && !std::is_same_v<matchJet, std::nullptr_t>) {
884-
for (const auto& matchedJet : matchjet(jet)) {
932+
if (jet.has_matchedJetGeo()) {
933+
registry.fill(HIST("mcm/hDetSparseMatch"), centrality, pt, jet.eta(), jet.phi(), weight);
934+
for (const auto& matchedJet : jet.template matchedJetGeo_as<MCPTab>()) {
935+
if (matchedJet.pt() > pTHatMaxMCD * pTHat)
936+
continue;
885937
auto matchedpt = matchedJet.pt();
886938
if (cfgbkgSubMC) {
887939
matchedpt = matchedJet.pt() - (rho2 * matchedJet.area());
888940
}
889-
registry.fill(HIST("mcm/hPartSparseMatch"), centrality, matchedpt, matchedJet.eta(), matchedJet.phi());
941+
registry.fill(HIST("mcm/hPartSparseMatch"), centrality, matchedpt, matchedJet.eta(), matchedJet.phi(), weight);
942+
registry.fill(HIST("mcm/hMatchedJetsPtDelta"), matchedJet.pt(), jet.pt() - matchedJet.pt(), weight);
943+
registry.fill(HIST("mcm/hMatchedJetsPhiDelta"), matchedJet.phi(), jet.phi() - matchedJet.phi(), weight);
944+
registry.fill(HIST("mcm/hMatchedJetsEtaDelta"), matchedJet.eta(), jet.eta() - matchedJet.eta(), weight);
945+
registry.fill(HIST("mcm/hGenMatchedJetsPtDeltadPt"), centrality, matchedpt, (pt - matchedpt) / matchedpt, weight);
946+
registry.fill(HIST("mcm/hRecoMatchedJetsPtDeltadPt"), centrality, pt, (pt - matchedpt) / pt, weight);
947+
registry.fill(HIST("mcm/hRespMcDMcPMatch"), centrality, pt, matchedpt, weight);
948+
}
949+
}
950+
}
951+
}
952+
953+
template <typename TTracks, typename TJets>
954+
bool isAcceptedJet(TJets const& jet, bool mcLevelIsParticleLevel = false)
955+
{
956+
if (jetAreaFractionMin > -98.0) {
957+
if (jet.area() < jetAreaFractionMin * o2::constants::math::PI * (jet.r() / 100.0) * (jet.r() / 100.0)) {
958+
return false;
959+
}
960+
}
961+
bool checkConstituentPt = true;
962+
bool checkConstituentMinPt = (leadingConstituentPtMin > -98.0);
963+
bool checkConstituentMaxPt = (leadingConstituentPtMax < 9998.0);
964+
if (!checkConstituentMinPt && !checkConstituentMaxPt) {
965+
checkConstituentPt = false;
966+
}
890967

891-
registry.fill(HIST("mcm/hMatchedJetsPtDelta"), matchedJet.pt(), jet.pt() - matchedJet.pt());
892-
registry.fill(HIST("mcm/hMatchedJetsPhiDelta"), matchedJet.phi(), jet.phi() - matchedJet.phi());
893-
registry.fill(HIST("mcm/hMatchedJetsEtaDelta"), matchedJet.eta(), jet.eta() - matchedJet.eta());
894-
registry.fill(HIST("mcm/hGenMatchedJetsPtDeltadPt"), centrality, matchedpt, (pt - matchedpt) / matchedpt);
895-
registry.fill(HIST("mcm/hRecoMatchedJetsPtDeltadPt"), centrality, pt, (pt - matchedpt) / pt);
968+
if (checkConstituentPt) {
969+
bool isMinLeadingConstituent = !checkConstituentMinPt;
970+
bool isMaxLeadingConstituent = true;
896971

897-
registry.fill(HIST("mcm/hRespMcDMcPMatch"), centrality, pt, matchedpt);
972+
for (const auto& constituent : jet.template tracks_as<TTracks>()) {
973+
double pt = constituent.pt();
974+
975+
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
976+
isMinLeadingConstituent = true;
977+
}
978+
if (checkConstituentMaxPt && pt > leadingConstituentPtMax) {
979+
isMaxLeadingConstituent = false;
898980
}
899981
}
982+
return isMinLeadingConstituent && isMaxLeadingConstituent;
900983
}
984+
return true;
901985
}
902986
};
903987
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<JetSpectraEseTask>(cfgc)}; }

0 commit comments

Comments
 (0)