diff --git a/PWGCF/TwoParticleCorrelations/Tasks/flowDecorrelation.cxx b/PWGCF/TwoParticleCorrelations/Tasks/flowDecorrelation.cxx index 89bb63cc0cd..692a6583f5c 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/flowDecorrelation.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/flowDecorrelation.cxx @@ -100,13 +100,24 @@ struct FlowDecorrelation { O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, bool, false, "Use local efficiency object") O2_DEFINE_CONFIGURABLE(cfgUseEventWeights, bool, false, "Use event weights for mixed event") O2_DEFINE_CONFIGURABLE(cfgDrawEtaPhiDis, bool, false, "draw eta-phi distribution for detectors in used") - O2_DEFINE_CONFIGURABLE(nClustersMftTrack, int, 5, "Minimum number of clusters for MFT track") - O2_DEFINE_CONFIGURABLE(cfgCutChi2Mft, float, -1.0f, "max chi2 of MFT track") - O2_DEFINE_CONFIGURABLE(cfgCutTrackTimeMft, float, -1.0f, "max deviation of MFT track wrt. bc in ns") - O2_DEFINE_CONFIGURABLE(cfgRejectFT0AInside, bool, false, "Rejection of inner ring channels of the FT0A detector") - O2_DEFINE_CONFIGURABLE(cfgRejectFT0AOutside, bool, false, "Rejection of outer ring channels of the FT0A detector") - O2_DEFINE_CONFIGURABLE(cfgRejectFT0CInside, bool, false, "Rejection of inner ring channels of the FT0C detector") - O2_DEFINE_CONFIGURABLE(cfgRejectFT0COutside, bool, false, "Rejection of outer ring channels of the FT0C detector") + struct : ConfigurableGroup { + O2_DEFINE_CONFIGURABLE(nClustersMftTrack, int, 5, "Minimum number of clusters for MFT track") + O2_DEFINE_CONFIGURABLE(cfgCutChi2Mft, float, -1.0f, "max chi2 of MFT track") + O2_DEFINE_CONFIGURABLE(cfgCutTrackTimeMft, float, -1.0f, "max deviation of MFT track wrt. bc in ns"); + } cfgMftCuts; + struct : ConfigurableGroup { + O2_DEFINE_CONFIGURABLE(cfgRejectFT0AInside, bool, false, "Rejection of inner ring channels of the FT0A detector") + O2_DEFINE_CONFIGURABLE(cfgRejectFT0AOutside, bool, false, "Rejection of outer ring channels of the FT0A detector") + O2_DEFINE_CONFIGURABLE(cfgRejectFT0CInside, bool, false, "Rejection of inner ring channels of the FT0C detector") + O2_DEFINE_CONFIGURABLE(cfgRejectFT0COutside, bool, false, "Rejection of outer ring channels of the FT0C detector"); + } cfgFt0RingRejections; + struct : ConfigurableGroup { + O2_DEFINE_CONFIGURABLE(cfgEtaTpcCut, float, 0.8f, "Eta cut of TPC MC particles") + O2_DEFINE_CONFIGURABLE(cfgMinEtaFt0cCut, float, -3.4f, "Min eta cut of FT0C MC particles") + O2_DEFINE_CONFIGURABLE(cfgMaxEtaFt0cCut, float, -2.0f, "Max eta cut of FT0C MC particles") + O2_DEFINE_CONFIGURABLE(cfgUseFt0cStructure, bool, true, "Use the true structure of FT0C in MC-true"); + O2_DEFINE_CONFIGURABLE(cfgUseCFStepAll, bool, true, "Use CFStepAll in addition to primry"); + } cfgMcTrue; struct : ConfigurableGroup { O2_DEFINE_CONFIGURABLE(cfgMultCentHighCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + 10.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); O2_DEFINE_CONFIGURABLE(cfgMultCentLowCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x - 3.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); @@ -147,8 +158,8 @@ struct FlowDecorrelation { ConfigurableAxis axisDeltaEtaTpcFt0c{"axisDeltaEtaTpcFt0c", {32, 1.2, 4.2}, "delta eta axis, 1.2~4.2 for TPC-FT0C"}; ConfigurableAxis axisDeltaEtaFt0aFt0c{"axisDeltaEtaFt0aFt0c", {32, 4.2, 8.2}, "delta eta axis, 4.2~8.2 for FT0A-FT0C"}; ConfigurableAxis axisDeltaEtaTpcMft{"axisDeltaEtaTpcMft", {32, 1.3, 4.8}, "delta eta axis, 1.3~4.8 for TPC-MFT"}; - ConfigurableAxis axisDeltaEtaTpcFv0{"axisDeltaEtaTpcFv0", {32, -1.2, -6.1}, "delta eta axis for TPC-FV0 histograms"}; - ConfigurableAxis axisEtaTrigger{"axisEtaTrigger", {VARIABLE_WIDTH, -3.3, -2.1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 3.5, 4.9}, "eta trigger axis for histograms"}; + ConfigurableAxis axisDeltaEtaTpcFv0{"axisDeltaEtaTpcFv0", {32, -6.1, -1.2}, "delta eta axis for TPC-FV0 histograms"}; + ConfigurableAxis axisEtaTrigger{"axisEtaTrigger", {VARIABLE_WIDTH, -3.4, -2.0, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 3.5, 4.9}, "eta trigger axis for histograms"}; ConfigurableAxis axisEtaAssoc{"axisEtaAssoc", {VARIABLE_WIDTH, -3.3, -2.1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6, 0.8, 3.5, 4.9}, "eta associated axis for histograms"}; ConfigurableAxis axisVtxMix{"axisVtxMix", {VARIABLE_WIDTH, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, "vertex axis for mixed event histograms"}; ConfigurableAxis axisMultMix{"axisMultMix", {VARIABLE_WIDTH, 0, 10, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260}, "multiplicity / centrality axis for mixed event histograms"}; @@ -171,6 +182,17 @@ struct FlowDecorrelation { using FilteredCollisions = soa::Filtered>; using FilteredTracks = soa::Filtered>; + Filter particleFilter = (nabs(aod::mcparticle::eta) < cfgMcTrue.cfgEtaTpcCut || ((aod::mcparticle::eta > cfgMcTrue.cfgMinEtaFt0cCut) && (aod::mcparticle::eta < cfgMcTrue.cfgMaxEtaFt0cCut))) && (aod::mcparticle::pt > cfgCutPtMin) && (aod::mcparticle::pt < cfgCutPtMax); + using FilteredMcParticles = soa::Filtered; + + // Filter for MCcollisions + Filter mccollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVtxZ; + using FilteredMcCollisions = soa::Filtered; + + using SmallGroupMcCollisions = soa::SmallGroups>; + + PresliceUnsorted collisionPerMCCollision = aod::mccollisionlabel::mcCollisionId; + // FT0 geometry o2::ft0::Geometry ft0Det; o2::fv0::Geometry* fv0Det{}; @@ -243,6 +265,7 @@ struct FlowDecorrelation { const AxisSpec axisPhi{72, 0.0, constants::math::TwoPI, "#varphi"}; const AxisSpec axisEta{40, -1., 1., "#eta"}; const AxisSpec axisEtaFull{90, -4., 5., "#eta"}; + const AxisSpec axisCentrality{20, 0., 100., "cent"}; ccdb->setURL("http://alice-ccdb.cern.ch"); ccdb->setCaching(true); @@ -268,7 +291,23 @@ struct FlowDecorrelation { registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(12, "MultCorrelation"); registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(13, "cfgEvSelV0AT0ACut"); } - + if (doprocessMcSameTpcFt0c) { + registry.add("MCTrue/MCeventcount", "MCeventcount", {HistType::kTH1F, {{5, 0, 5, "bin"}}}); // histogram to see how many events are in the same and mixed event + registry.get(HIST("MCTrue/MCeventcount"))->GetXaxis()->SetBinLabel(2, "same all"); + registry.get(HIST("MCTrue/MCeventcount"))->GetXaxis()->SetBinLabel(3, "same reco"); + registry.get(HIST("MCTrue/MCeventcount"))->GetXaxis()->SetBinLabel(4, "mixed all"); + registry.get(HIST("MCTrue/MCeventcount"))->GetXaxis()->SetBinLabel(5, "mixed reco"); + registry.add("MCTrue/MCCentrality", "cent", {HistType::kTH1D, {axisCentrality}}); + registry.add("MCTrue/MCNch", "N_{ch}", {HistType::kTH1D, {axisMultiplicity}}); + registry.add("MCTrue/MCzVtx", "MCzVtx", {HistType::kTH1D, {axisVertex}}); + registry.add("MCTrue/MCPhi", "MCPhi", {HistType::kTH1D, {axisPhi}}); + registry.add("MCTrue/MCEta", "MCEta", {HistType::kTH1D, {axisEtaFull}}); + registry.add("MCTrue/MCEtaTrueShape", "MCEta", {HistType::kTH1D, {axisEtaFull}}); + registry.add("MCTrue/MCpT", "MCpT", {HistType::kTH1D, {axisPtEfficiency}}); + registry.add("MCTrue/MCTrig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtEfficiency}}}); + registry.add("MCTrue/MCdeltaEta_deltaPhi_same", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFt0c}}); // check to see the delta eta and delta phi distribution + registry.add("MCTrue/MCdeltaEta_deltaPhi_mixed", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFt0c}}); + } if (cfgEvSelMultCorrelation) { cfgFuncParas.multT0CCutPars = cfgFuncParas.cfgMultT0CCutPars; cfgFuncParas.multPVT0CCutPars = cfgFuncParas.cfgMultPVT0CCutPars; @@ -421,6 +460,10 @@ struct FlowDecorrelation { same.setObject(new CorrelationContainer("sameEvent_TPC_FV0", "sameEvent_TPC_FV0", corrAxisTpcFv0, effAxis, userAxis)); mixed.setObject(new CorrelationContainer("mixedEvent_TPC_FV0", "mixedEvent_TPC_FV0", corrAxisTpcFv0, effAxis, userAxis)); } + if (doprocessMcSameTpcFt0c) { + same.setObject(new CorrelationContainer("sameEvent_TPC_FV0", "sameEvent_TPC_FV0", corrAxisTpcFt0c, effAxis, userAxis)); + mixed.setObject(new CorrelationContainer("mixedEvent_TPC_FV0", "mixedEvent_TPC_FV0", corrAxisTpcFt0c, effAxis, userAxis)); + } LOGF(info, "End of init"); } @@ -428,13 +471,13 @@ struct FlowDecorrelation { bool isAcceptedMftTrack(TTrackAssoc const& mftTrack) { // cut on the number of clusters of the reconstructed MFT track - if (mftTrack.nClusters() < nClustersMftTrack) + if (mftTrack.nClusters() < cfgMftCuts.nClustersMftTrack) return false; - if (cfgCutChi2Mft > 0. && mftTrack.chi2() > cfgCutChi2Mft) + if (cfgMftCuts.cfgCutChi2Mft > 0. && mftTrack.chi2() > cfgMftCuts.cfgCutChi2Mft) return false; - if (cfgCutTrackTimeMft > 0. && std::abs(mftTrack.trackTime()) > cfgCutTrackTimeMft) + if (cfgMftCuts.cfgCutTrackTimeMft > 0. && std::abs(mftTrack.trackTime()) > cfgMftCuts.cfgCutTrackTimeMft) return false; return true; @@ -729,11 +772,11 @@ struct FlowDecorrelation { float ampl = 0.; getChannel(ft0, iCh, chanelid, ampl, corType); if (corType == kFT0C) { - if ((cfgRejectFT0CInside && (chanelid >= kFT0CInnerRingMin && chanelid <= kFT0CInnerRingMax)) || (cfgRejectFT0COutside && (chanelid >= kFT0COuterRingMin && chanelid <= kFT0COuterRingMax))) { + if ((cfgFt0RingRejections.cfgRejectFT0CInside && (chanelid >= kFT0CInnerRingMin && chanelid <= kFT0CInnerRingMax)) || (cfgFt0RingRejections.cfgRejectFT0COutside && (chanelid >= kFT0COuterRingMin && chanelid <= kFT0COuterRingMax))) { continue; } } else if (corType == kFT0A) { - if ((cfgRejectFT0AInside && (chanelid >= kFT0AInnerRingMin && chanelid <= kFT0AInnerRingMax)) || (cfgRejectFT0AOutside && (chanelid >= kFT0AOuterRingMin && chanelid <= kFT0AOuterRingMax))) { + if ((cfgFt0RingRejections.cfgRejectFT0AInside && (chanelid >= kFT0AInnerRingMin && chanelid <= kFT0AInnerRingMax)) || (cfgFt0RingRejections.cfgRejectFT0AOutside && (chanelid >= kFT0AOuterRingMin && chanelid <= kFT0AOuterRingMax))) { continue; } } @@ -786,7 +829,7 @@ struct FlowDecorrelation { int channelIdA = 0; float amplA = 0.f; getChannel(ft0Trig, iChA, channelIdA, amplA, kFT0A); - if ((cfgRejectFT0AInside && (channelIdA >= kFT0AInnerRingMin && channelIdA <= kFT0AInnerRingMax)) || (cfgRejectFT0AOutside && (channelIdA >= kFT0AOuterRingMin && channelIdA <= kFT0AOuterRingMax))) { + if ((cfgFt0RingRejections.cfgRejectFT0AInside && (channelIdA >= kFT0AInnerRingMin && channelIdA <= kFT0AInnerRingMax)) || (cfgFt0RingRejections.cfgRejectFT0AOutside && (channelIdA >= kFT0AOuterRingMin && channelIdA <= kFT0AOuterRingMax))) { continue; } @@ -801,7 +844,7 @@ struct FlowDecorrelation { int channelIdC = 0; float amplC = 0.f; getChannel(ft0Assoc, iChC, channelIdC, amplC, kFT0C); - if ((cfgRejectFT0CInside && (channelIdC >= kFT0CInnerRingMin && channelIdC <= kFT0CInnerRingMax)) || (cfgRejectFT0COutside && (channelIdC >= kFT0COuterRingMin && channelIdC <= kFT0COuterRingMax))) { + if ((cfgFt0RingRejections.cfgRejectFT0CInside && (channelIdC >= kFT0CInnerRingMin && channelIdC <= kFT0CInnerRingMax)) || (cfgFt0RingRejections.cfgRejectFT0COutside && (channelIdC >= kFT0COuterRingMin && channelIdC <= kFT0COuterRingMax))) { continue; } @@ -983,6 +1026,24 @@ struct FlowDecorrelation { return 1; } + bool ft0cCollisionCourse(float eta, float phi, float vtxZ) + { + double theta = 2 * std::atan(std::exp(-eta)); + double vx = std::sin(theta) * std::cos(phi); // veloctiy component along x + double vy = std::sin(theta) * std::sin(phi); // veloctiy component along y + double vz = std::cos(theta); // veloctiy component along z + + double x = vx * ((-83.44 - vtxZ) / vz); // location of particle on x at FT0C distance from vertex + double y = vy * ((-83.44 - vtxZ) / vz); // location of particle on x at FT0C distance from vertex + + if (std::abs(x) < 24 && (std::abs(y) > 6.825 && std::abs(y) < 18.25)) + return true; + else if (std::abs(y) < 24 && (std::abs(x) > 6.675 && std::abs(x) < 18.175)) + return true; + else + return false; + } + void processSameTpcFt0a(FilteredCollisions::iterator const& collision, FilteredTracks const& tracks, aod::FT0s const&, aod::BCsWithTimestamps const&) { if (!collision.sel8()) @@ -1558,6 +1619,149 @@ struct FlowDecorrelation { } } PROCESS_SWITCH(FlowDecorrelation, processMixedTpcFv0, "Process mixed events for TPC-FV0 correlation", false); + + template + void fillMCCorrelations(TTracks tracks1, TTracksAssoc tracks2, float posZ, int system, float eventWeight) // function to fill the Output functions (sparse) and the delta eta and delta phi histograms + { + int fSampleIndex = gRandom->Uniform(0, cfgSampleSize); + + float triggerWeight = 1.0f; + float associatedWeight = 1.0f; + // loop over all FT0C tracks + for (auto const& track1 : tracks1) { + + if (!cfgMcTrue.cfgUseFt0cStructure) { + if (std::abs(track1.eta()) < 0.9) + continue; + } else if (!ft0cCollisionCourse(track1.eta(), track1.phi(), posZ)) { + continue; + } + + if (step >= CorrelationContainer::kCFStepTrackedOnlyPrim && !track1.isPhysicalPrimary()) + continue; + + if (step >= CorrelationContainer::kCFStepTrackedOnlyPrim && system == SameEvent) + registry.fill(HIST("MCTrue/MCEtaTrueShape"), track1.eta()); + + if (system == SameEvent && doprocessMcSameTpcFt0c) + registry.fill(HIST("MCTrue/MCTrig_hist"), fSampleIndex, posZ, track1.pt(), eventWeight * triggerWeight); + + // loop over all TPC tracks + for (auto const& track2 : tracks2) { + + if (std::abs(track2.eta()) > 0.9) + continue; + + if (step >= CorrelationContainer::kCFStepTrackedOnlyPrim && !track2.isPhysicalPrimary()) + continue; + + if (track1.globalIndex() == track2.globalIndex()) + continue; // For pt-differential correlations, skip if the trigger and associate are the same track + + float deltaPhi = RecoDecay::constrainAngle(track2.phi() - track1.phi(), -PIHalf); + float deltaEta = track2.eta() - track1.eta(); + + // fill the right sparse and histograms + if (system == SameEvent) { + same->getPairHist()->Fill(step, fSampleIndex, posZ, track2.eta(), track1.eta(), deltaPhi, deltaEta, eventWeight * triggerWeight * associatedWeight); + registry.fill(HIST("MCTrue/MCdeltaEta_deltaPhi_same"), deltaPhi, deltaEta, eventWeight * triggerWeight * associatedWeight); + } else if (system == MixedEvent) { + mixed->getPairHist()->Fill(step, fSampleIndex, posZ, track2.eta(), track1.eta(), deltaPhi, deltaEta, eventWeight * triggerWeight * associatedWeight); + registry.fill(HIST("MCTrue/MCdeltaEta_deltaPhi_mixed"), deltaPhi, deltaEta, eventWeight * triggerWeight * associatedWeight); + } + } + } + } + + void processMcSameTpcFt0c(FilteredMcCollisions::iterator const& mcCollision, FilteredMcParticles const& mcParticles, SmallGroupMcCollisions const& collisions) + { + float cent = -1; + if (!cfgCentTableUnavailable) { + for (const auto& collision : collisions) { + cent = getCentrality(collision); + } + } + + if (cfgSelCollByNch && (mcParticles.size() < cfgCutMultMin || mcParticles.size() >= cfgCutMultMax)) { + return; + } + if (!cfgSelCollByNch && !cfgCentTableUnavailable && (cent < cfgCutCentMin || cent >= cfgCutCentMax)) { + return; + } + + registry.fill(HIST("MCTrue/MCeventcount"), SameEvent); // because its same event i put it in the 1 bin + if (!cfgCentTableUnavailable) + registry.fill(HIST("MCTrue/MCCentrality"), cent); + registry.fill(HIST("MCTrue/MCNch"), mcParticles.size()); + registry.fill(HIST("MCTrue/MCzVtx"), mcCollision.posZ()); + for (const auto& mcParticle : mcParticles) { + if (mcParticle.isPhysicalPrimary()) { + registry.fill(HIST("MCTrue/MCPhi"), mcParticle.phi()); + registry.fill(HIST("MCTrue/MCEta"), mcParticle.eta()); + registry.fill(HIST("MCTrue/MCpT"), mcParticle.pt()); + } + } + if (cfgMcTrue.cfgUseCFStepAll) { + same->fillEvent(mcParticles.size(), CorrelationContainer::kCFStepAll); + fillMCCorrelations(mcParticles, mcParticles, mcCollision.posZ(), SameEvent, 1.0f); + } + + registry.fill(HIST("MCTrue/MCeventcount"), 2.5); + same->fillEvent(mcParticles.size(), CorrelationContainer::kCFStepTrackedOnlyPrim); + fillMCCorrelations(mcParticles, mcParticles, mcCollision.posZ(), SameEvent, 1.0f); + } + PROCESS_SWITCH(FlowDecorrelation, processMcSameTpcFt0c, "Process MC same event", false); + + void processMcMixed(FilteredMcCollisions const& mcCollisions, FilteredMcParticles const& mcParticles, SmallGroupMcCollisions const& collisions) + { + auto getTracksSize = [&mcParticles, this](FilteredMcCollisions::iterator const& mcCollision) { + auto associatedTracks = mcParticles.sliceByCached(o2::aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), this->cache); + auto mult = associatedTracks.size(); + return mult; + }; + + using MixedBinning = FlexibleBinningPolicy, o2::aod::mccollision::PosZ, decltype(getTracksSize)>; + + MixedBinning binningOnVtxAndMult{{getTracksSize}, {axisVtxMix, axisMultMix}, true}; + + auto tracksTuple = std::make_tuple(mcParticles, mcParticles); + Pair pairs{binningOnVtxAndMult, cfgMixEventNumMin, -1, mcCollisions, tracksTuple, &cache}; // -1 is the number of the bin to skip + for (auto it = pairs.begin(); it != pairs.end(); it++) { + auto& [collision1, tracks1, collision2, tracks2] = *it; + + if (cfgSelCollByNch && (tracks1.size() < cfgCutMultMin || tracks1.size() >= cfgCutMultMax)) + continue; + + if (cfgSelCollByNch && (tracks2.size() < cfgCutMultMin || tracks2.size() >= cfgCutMultMax)) + continue; + + auto groupedCollisions = collisions.sliceBy(collisionPerMCCollision, collision1.globalIndex()); + + float cent = -1; + if (!cfgCentTableUnavailable) { + for (const auto& collision : groupedCollisions) { + cent = getCentrality(collision); + } + } + + if (!cfgSelCollByNch && !cfgCentTableUnavailable && (cent < cfgCutCentMin || cent >= cfgCutCentMax)) + continue; + + registry.fill(HIST("MCTrue/MCeventcount"), MixedEvent); // fill the mixed event in the 3 bin + float eventWeight = 1.0f; + if (cfgUseEventWeights) { + eventWeight = 1.0f / it.currentWindowNeighbours(); + } + + if (cfgMcTrue.cfgUseCFStepAll) { + fillMCCorrelations(tracks1, tracks2, collision1.posZ(), MixedEvent, eventWeight); + } + + registry.fill(HIST("MCTrue/MCeventcount"), 4.5); + fillMCCorrelations(tracks1, tracks2, collision1.posZ(), MixedEvent, eventWeight); + } + } + PROCESS_SWITCH(FlowDecorrelation, processMcMixed, "Process MC mixed events", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)