diff --git a/PWGHF/HFC/Tasks/taskFlow.cxx b/PWGHF/HFC/Tasks/taskFlow.cxx index 50f04ea6d69..feff1396b0a 100644 --- a/PWGHF/HFC/Tasks/taskFlow.cxx +++ b/PWGHF/HFC/Tasks/taskFlow.cxx @@ -225,6 +225,7 @@ struct HfTaskFlow { Configurable nSamples{"nSamples", 10, "number of different samples for correlations"}; Configurable nameCorrelationContainer{"nameCorrelationContainer", "", "Add the possibility to the rename the correlation container as configurable"}; Configurable useEfficiencyCorrection{"useEfficiencyCorrection", false, "Choose to use or not use efficiency correction, if not used, weight is set to 1"}; + Configurable useOnlyPrimaryInMc{"useOnlyPrimaryInMc", false, "Choose to use or not use only primaries for McGen"}; } configTask; // configurables for collisions @@ -318,10 +319,10 @@ struct HfTaskFlow { Service pdg{}; Service ccdb{}; std::vector* offsetFT0{}; - std::vector* offsetFV0{}; + // std::vector* offsetFV0{}; o2::ccdb::CcdbApi ccdbApi; o2::ft0::Geometry ft0Det; - o2::fv0::Geometry* fv0Det{}; + // o2::fv0::Geometry* fv0Det{}; std::vector cstFT0RelGain{}; RCTFlagsChecker rctChecker; RCTFlagsChecker correlationAnalysisRctChecker{kFT0Bad, kITSBad, kTPCBadTracking, kTPCBadPID, kMFTBad, kITSLimAccMCRepr, kMFTLimAccMCRepr, kTPCLimAccMCRepr}; @@ -349,6 +350,8 @@ struct HfTaskFlow { // ========================= using SmallGroupMcCollisions = soa::SmallGroups>; + using FilteredMcCollisionsWMult = soa::Filtered>; + using FilteredMcParticles = soa::Filtered; // ========================= // Filters & partitions : DATA @@ -383,11 +386,19 @@ struct HfTaskFlow { Filter mftTrackEtaFilter = ((aod::fwdtrack::eta < configMft.etaMftTrackMaxFilter) && (aod::fwdtrack::eta > configMft.etaMftTrackMinFilter)); - // Filters below will be used for uncertainties Filter mftTrackCollisionIdFilter = (aod::fwdtrack::bestCollisionId >= 0); // Filter mftTrackDcaXYFilter = (nabs(aod::fwdtrack::bestDCAXY) < configMft.mftMaxDCAxy); // Filter mftTrackDcaZFilter = (nabs(aod::fwdtrack::bestDCAZ) < configMft.mftMaxDCAz); + // ========================= + // Filters & partitions : MC + // ========================= + + Filter mcParticleFilter = (((aod::mcparticle::eta > configTask.etaMcParticlesTriggerMin) && (aod::mcparticle::eta < configTask.etaMcParticlesTriggerMax)) || ((aod::mcparticle::eta > configTask.etaMcParticlesAssocMin) && (aod::mcparticle::eta < configTask.etaMcParticlesAssocMax)) || (nabs(aod::mcparticle::eta) < configCentral.etaCentralTrackMax)) && (aod::mcparticle::pt > configTask.ptMcParticlesTriggerMin) && (aod::mcparticle::pt < configTask.ptMcParticlesTriggerMax); + + // Filter for MCcollisions + Filter mcCollisionFilter = nabs(aod::mccollision::posZ) < configCollision.zVertexMax; + // ========================= // Preslice : DATA // ========================= @@ -396,7 +407,7 @@ struct HfTaskFlow { Preslice perColLcs = aod::track::collisionId; Preslice perColMftTracks = o2::aod::fwdtrack::collisionId; Preslice perColTracks = aod::track::collisionId; - Preslice perMcColMcParticles = aod::mcparticle::mcCollisionId; + // Preslice perMcColMcParticles = aod::mcparticle::mcCollisionId; PresliceUnsorted> perColReassociated2dTracks = o2::aod::fwdtrack::collisionId; PresliceUnsorted> perColReassociated3dTracks = o2::aod::fwdtrack::collisionId; @@ -524,12 +535,12 @@ struct HfTaskFlow { ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); LOGF(info, "Getting alignment offsets from the CCDB..."); offsetFT0 = ccdb->getForTimeStamp>("FT0/Calib/Align", configCcdb.noLaterThan.value); - offsetFV0 = ccdb->getForTimeStamp>("FV0/Calib/Align", configCcdb.noLaterThan.value); + // offsetFV0 = ccdb->getForTimeStamp>("FV0/Calib/Align", configCcdb.noLaterThan.value); LOGF(info, "Offset for FT0A: x = %.3f y = %.3f z = %.3f\n", (*offsetFT0)[0].getX(), (*offsetFT0)[0].getY(), (*offsetFT0)[0].getZ()); LOGF(info, "Offset for FT0C: x = %.3f y = %.3f z = %.3f\n", (*offsetFT0)[1].getX(), (*offsetFT0)[1].getY(), (*offsetFT0)[1].getZ()); - LOGF(info, "Offset for FV0-left: x = %.3f y = %.3f z = %.3f\n", (*offsetFV0)[0].getX(), (*offsetFV0)[0].getY(), (*offsetFV0)[0].getZ()); - LOGF(info, "Offset for FV0-right: x = %.3f y = %.3f z = %.3f\n", (*offsetFV0)[1].getX(), (*offsetFV0)[1].getY(), (*offsetFV0)[1].getZ()); - fv0Det = o2::fv0::Geometry::instance(o2::fv0::Geometry::eUninitialized); + // LOGF(info, "Offset for FV0-left: x = %.3f y = %.3f z = %.3f\n", (*offsetFV0)[0].getX(), (*offsetFV0)[0].getY(), (*offsetFV0)[0].getZ()); + // LOGF(info, "Offset for FV0-right: x = %.3f y = %.3f z = %.3f\n", (*offsetFV0)[1].getX(), (*offsetFV0)[1].getY(), (*offsetFV0)[1].getZ()); + // fv0Det = o2::fv0::Geometry::instance(o2::fv0::Geometry::eUninitialized); // ========================= // Event histograms @@ -650,6 +661,9 @@ struct HfTaskFlow { registry.add("Data/hEfficiencyTrigger", "", {HistType::kTH3D, {{configAxis.axisPtTrigger}, {configAxis.axisEtaTrigger}, {configAxis.axisVertex}}}); registry.add("Data/hEfficiencyAssociated", "", {HistType::kTH3D, {{configAxis.axisPtAssoc}, {configAxis.axisEtaAssociated}, {configAxis.axisVertex}}}); + registry.add("Data/hMultiplicity_uncorrected", "", {HistType::kTH1D, {configAxis.axisMultiplicity}}); + registry.add("Data/hMultiplicity_corrected", "", {HistType::kTH1D, {configAxis.axisMultiplicity}}); + if (!configTask.doEtaDependentFlow && !configTask.doVariationContainers) { registry.add("Trig_hist_TPC_MFT", "", {HistType::kTHnSparseF, {{configAxis.axisSamples, configAxis.axisVertex, configAxis.axisPtTrigger}}}); sameEvent.setObject(new CorrelationContainer("sameEvent", "sameEvent", corrAxis, effAxis, {})); @@ -852,7 +866,7 @@ struct HfTaskFlow { sameEvent.setObject(new CorrelationContainer("sameEvent", "sameEvent", corrAxis, effAxis, {})); mixedEvent.setObject(new CorrelationContainer("mixedEvent", "mixedEvent", corrAxis, effAxis, {})); } else if (configTask.doEtaDependentFlow) { - registry.add("Trig_hist_FT0A_FT0C ", "", {HistType::kTHnSparseF, {{configAxis.axisSamples, configAxis.axisVertex, configAxis.axisEtaTrigger}}}); + registry.add("Trig_hist_FT0A_FT0C", "", {HistType::kTHnSparseF, {{configAxis.axisSamples, configAxis.axisVertex, configAxis.axisEtaTrigger}}}); sameEvent.setObject(new CorrelationContainer("sameEvent", "sameEvent", corrAxisEta, effAxis, {})); mixedEvent.setObject(new CorrelationContainer("mixedEvent", "mixedEvent", corrAxisEta, effAxis, {})); } else { @@ -881,6 +895,8 @@ struct HfTaskFlow { addHistograms(); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::TpcFt0a)) { addHistograms(); + } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::MftFt0a)) { + addHistograms(); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::TpcFt0c)) { addHistograms(); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::Ft0aFt0c)) { @@ -1216,9 +1232,6 @@ struct HfTaskFlow { trackCounter += 1; } - // if (!getEfficiencyCorrection_Nch(weight_Nch, track.pt())) { - // continue; - // } auto efficiencyNch = 1.0f; if (configTask.useEfficiencyCorrection) { @@ -1241,24 +1254,6 @@ struct HfTaskFlow { return trackCounter; } - // bool getEfficiencyCorrectionNch(float& weightNch, float pt) - // { - // float efficiencyNch = 1.; - - // if (mEfficiencyNch) { - // int ptBin = mEfficiencyNch->FindBin(pt); - // efficiencyNch = mEfficiencyNch->GetBinContent(ptBin); - // } - - // if (efficiencyNch == 0 ) { - // return false; - // } - - // weightNch = 1. / efficiencyNch; - - // return true; - // } - // ========================= // Cuts with functions // ========================= @@ -2380,8 +2375,8 @@ struct HfTaskFlow { TTracksTrig const& tracks1, TTracksAssoc const& tracks2, float multiplicity, float posZ, bool sameEvent) { - auto triggerWeight = 1; - auto associatedWeight = 1; + auto triggerWeight = 1.0f; + auto associatedWeight = 1.0f; auto loopCounter = 0; // To avoid filling associated tracks QA many times, I fill it only for the first trigger track of the collision int sampleIndex = gRandom->Uniform(0, configTask.nSamples); bool fillQaPlots = false; @@ -2405,7 +2400,7 @@ struct HfTaskFlow { if (track1.pt() < configTask.ptMcParticlesTriggerMin || track1.pt() > configTask.ptMcParticlesTriggerMax) { continue; } - if (step >= CorrelationContainer::kCFStepTrackedOnlyPrim && !track1.isPhysicalPrimary()) { + if (configTask.useOnlyPrimaryInMc && !track1.isPhysicalPrimary()) { continue; } @@ -2431,6 +2426,8 @@ struct HfTaskFlow { fillTriggerQa(multiplicity, track1.eta(), track1.phi(), track1.pt()); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::TpcFt0a)) { fillTriggerQa(multiplicity, track1.eta(), track1.phi(), track1.pt()); + } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::MftFt0a)) { + fillTriggerQa(multiplicity, track1.eta(), track1.phi(), track1.pt()); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::TpcFt0c)) { fillTriggerQa(multiplicity, track1.eta(), track1.phi(), track1.pt()); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::Ft0aFt0c)) { @@ -2449,7 +2446,7 @@ struct HfTaskFlow { if (track2.pt() < configTask.ptMcParticlesAssocMin || track2.pt() > configTask.ptMcParticlesAssocMax) { continue; } - if (step >= CorrelationContainer::kCFStepTrackedOnlyPrim && !track2.isPhysicalPrimary()) { + if (configTask.useOnlyPrimaryInMc && !track2.isPhysicalPrimary()) { continue; } @@ -2471,19 +2468,21 @@ struct HfTaskFlow { if (sameEvent && fillQaPlots && (loopCounter == 1)) { registry.fill(HIST("MC/hEfficiencyAssociated"), track2.pt(), track2.eta(), posZ); if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::TpcTpc)) { - fillAssociatedQa(track2.eta(), track2.phi(), track2.pt()); + fillAssociatedQa(multiplicity, track2.eta(), track2.phi()); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::TpcMft)) { - fillAssociatedQa(track2.eta(), track2.phi(), track2.pt()); + fillAssociatedQa(multiplicity, track2.eta(), track2.phi()); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::TpcFv0a)) { - fillAssociatedQa(track2.eta(), track2.phi(), track2.pt()); + fillAssociatedQa(multiplicity, track2.eta(), track2.phi()); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::MftFv0a)) { - fillAssociatedQa(track2.eta(), track2.phi(), track2.pt()); + fillAssociatedQa(multiplicity, track2.eta(), track2.phi()); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::TpcFt0a)) { - fillAssociatedQa(track2.eta(), track2.phi(), track2.pt()); + fillAssociatedQa(multiplicity, track2.eta(), track2.phi()); + } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::MftFt0a)) { + fillAssociatedQa(multiplicity, track2.eta(), track2.phi()); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::TpcFt0c)) { - fillAssociatedQa(track2.eta(), track2.phi(), track2.pt()); + fillAssociatedQa(multiplicity, track2.eta(), track2.phi()); } else if (configTask.chooseCorrelationCase.value == static_cast(CorrelationCase::Ft0aFt0c)) { - fillAssociatedQa(track2.eta(), track2.phi(), track2.pt()); + fillAssociatedQa(multiplicity, track2.eta(), track2.phi()); } } // end of fill QA } // end of loop over track2 @@ -3020,6 +3019,12 @@ struct HfTaskFlow { multiplicity = getMultiplicityEstimator(collision, true); } + registry.fill(HIST("Data/hMultiplicity_uncorrected"), multiplicity); + if (configCollision.useMultiplicityFromTracksCorrected) { + multiplicity = getCorrectedMultiplicity(tracks); + registry.fill(HIST("Data/hMultiplicity_corrected"), multiplicity); + } + if (multiplicity < configCollision.minMultiplicity || multiplicity >= configCollision.maxMultiplicity) { return; } @@ -3809,14 +3814,10 @@ struct HfTaskFlow { // MONTE-CARLO // =================================================================================================================================================================================================================================================================== - void processSameMcGen(soa::Join::iterator const& mcCollision, - aod::McParticles const& mcParticles, + void processSameMcGen(FilteredMcCollisionsWMult::iterator const& mcCollision, + FilteredMcParticles const& mcParticles, SmallGroupMcCollisions const& collisions) { - if (!(std::abs(mcCollision.posZ()) <= configCollision.zVertexMax)) { - return; - } - auto multiplicity = 0; if (configCollision.useMultiplicityFromTracks) { for (const auto& track : mcParticles) { @@ -3832,15 +3833,12 @@ struct HfTaskFlow { return; } - sameEvent->fillEvent(mcParticles.size(), CorrelationContainer::kCFStepAll); - fillCorrelationsMonteCarlo(sameEvent, CorrelationContainer::CFStep::kCFStepAll, mcParticles, mcParticles, multiplicity, mcCollision.posZ(), true); - if (collisions.size() == 0) { return; } - sameEvent->fillEvent(mcParticles.size(), CorrelationContainer::kCFStepTrackedOnlyPrim); - fillCorrelationsMonteCarlo(sameEvent, CorrelationContainer::CFStep::kCFStepTrackedOnlyPrim, mcParticles, mcParticles, multiplicity, mcCollision.posZ(), true); + sameEvent->fillEvent(multiplicity, CorrelationContainer::kCFStepAll); + fillCorrelationsMonteCarlo(sameEvent, CorrelationContainer::CFStep::kCFStepAll, mcParticles, mcParticles, multiplicity, mcCollision.posZ(), true); } PROCESS_SWITCH(HfTaskFlow, processSameMcGen, "MC Gen : Process same-event correlations", false); @@ -4205,15 +4203,10 @@ struct HfTaskFlow { // MONTE-CARLO // =================================================================================================================================================================================================================================================================== - void processMixedMcGen(soa::Join const& mcCollisions, - aod::McParticles const& mcParticles, + void processMixedMcGen(FilteredMcCollisionsWMult const& mcCollisions, + FilteredMcParticles const& mcParticles, SmallGroupMcCollisions const& collisions) { - // auto getTracksSize = [&mcParticles, this](aod::McCollisions::iterator const& mcCollision) { - // auto associatedTracks = mcParticles.sliceByCached(o2::aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), this->cache); - // auto multiplicity = associatedTracks.size(); - // return multiplicity; - // }; auto getTracksSize = [&mcParticles, this](soa::Join::iterator const& mcCollision) { auto associatedTracks = mcParticles.sliceByCached(o2::aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), this->cache); @@ -4233,10 +4226,12 @@ struct HfTaskFlow { using MixedBinning = FlexibleBinningPolicy, aod::mccollision::PosZ, decltype(getTracksSize)>; MixedBinning const binningOnVtxAndMult{{getTracksSize}, {configAxis.binsMixingVertex, configAxis.binsMixingMultiplicity}, true}; - for (auto const& [collision1, collision2] : soa::selfCombinations(binningOnVtxAndMult, configTask.nMixedEvents, -1, mcCollisions, mcCollisions)) { + auto tracksTuple = std::make_tuple(mcParticles, mcParticles); + + Pair pairs{binningOnVtxAndMult, configTask.nMixedEvents, -1, mcCollisions, tracksTuple, &cache}; // -1 is the number of the bin to skip - auto tracks1 = mcParticles.sliceBy(perMcColMcParticles, collision1.globalIndex()); - auto tracks2 = mcParticles.sliceBy(perMcColMcParticles, collision2.globalIndex()); + for (auto it = pairs.begin(); it != pairs.end(); it++) { + auto& [collision1, tracks1, collision2, tracks2] = *it; auto multiplicityCollision1 = 0; auto multiplicityCollision2 = 0; @@ -4264,17 +4259,17 @@ struct HfTaskFlow { continue; } - auto groupedCollisions = collisions.sliceBy(collisionPerMcCollision, collision1.globalIndex()); - - mixedEvent->fillEvent(mcParticles.size(), CorrelationContainer::kCFStepAll); - fillCorrelationsMonteCarlo(mixedEvent, CorrelationContainer::CFStep::kCFStepAll, mcParticles, mcParticles, multiplicityCollision1, collision1.posZ(), false); - - if (groupedCollisions.size() == 0) { + auto groupedCollisions1 = collisions.sliceBy(collisionPerMcCollision, collision1.globalIndex()); + auto groupedCollisions2 = collisions.sliceBy(collisionPerMcCollision, collision2.globalIndex()); + if (groupedCollisions1.size() == 0) { + return; + } + if (groupedCollisions2.size() == 0) { return; } - mixedEvent->fillEvent(mcParticles.size(), CorrelationContainer::kCFStepTrackedOnlyPrim); - fillCorrelationsMonteCarlo(mixedEvent, CorrelationContainer::CFStep::kCFStepTrackedOnlyPrim, mcParticles, mcParticles, multiplicityCollision1, collision1.posZ(), false); + mixedEvent->fillEvent(multiplicityCollision1, CorrelationContainer::kCFStepAll); + fillCorrelationsMonteCarlo(mixedEvent, CorrelationContainer::CFStep::kCFStepAll, tracks1, tracks2, multiplicityCollision1, collision1.posZ(), false); } } PROCESS_SWITCH(HfTaskFlow, processMixedMcGen, "MC Gen : Process mixed-event correlations", false);