diff --git a/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt b/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt index c8484a3dfa1..e95f85f0f87 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt +++ b/PWGCF/TwoParticleCorrelations/Tasks/CMakeLists.txt @@ -43,8 +43,8 @@ o2physics_add_dpl_workflow(dpt-dpt-per-run-extra-qc PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore O2Physics::AnalysisCCDB COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(corr-sparse - SOURCES corrSparse.cxx +o2physics_add_dpl_workflow(corr-reso + SOURCES corrReso.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore COMPONENT_NAME Analysis) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/corrFit.cxx b/PWGCF/TwoParticleCorrelations/Tasks/corrFit.cxx index b78d6f28463..3db192cc17c 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/corrFit.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/corrFit.cxx @@ -855,20 +855,20 @@ struct CorrFit { correctionsLoaded = true; } - bool getEfficiencyCorrection_Nch(float& weight_Nch, float pt) + bool getEfficiencyCorrectionNch(float& weightNch, float pt) { - float eff_Nch = 1.; + float effNch = 1.; if (mEfficiencyNch) { int ptBin = mEfficiencyNch->FindBin(pt); - eff_Nch = mEfficiencyNch->GetBinContent(ptBin); + effNch = mEfficiencyNch->GetBinContent(ptBin); } else { - eff_Nch = 1.0; + effNch = 1.0; } - if (eff_Nch == 0) + if (effNch == 0) return false; - weight_Nch = 1. / eff_Nch; + weightNch = 1. / effNch; return true; } @@ -910,7 +910,7 @@ struct CorrFit { void trackCounter(TTracks tracks, double& multiplicity) // function to count the number of tracks in the event and fill the histogram { double nTracksCorrected = 0; - float weight_Nch = 1.0f; + float weightNch = 1.0f; for (auto const& track : tracks) { if (cfgRefMultiplicity) { @@ -918,11 +918,11 @@ struct CorrFit { continue; } - if (!getEfficiencyCorrection_Nch(weight_Nch, track.pt())) { + if (!getEfficiencyCorrectionNch(weightNch, track.pt())) { continue; } - nTracksCorrected += weight_Nch; + nTracksCorrected += weightNch; } multiplicity = nTracksCorrected; } @@ -956,7 +956,7 @@ struct CorrFit { } template - void fillNsigmaBeforeCut(TTrack track1, Int_t pid) // function to fill the QA before Nsigma selection + void fillNsigmaBeforeCut(TTrack track1, int pid) // function to fill the QA before Nsigma selection { switch (pid) { case kPions: // For Pions @@ -1008,7 +1008,7 @@ struct CorrFit { } // end of fillNsigmaBeforeCut template - void fillNsigmaAfterCut(TTrack track1, Int_t pid) // function to fill the QA after Nsigma selection + void fillNsigmaAfterCut(TTrack track1, int pid) // function to fill the QA after Nsigma selection { switch (pid) { case kPions: // For Pions @@ -1221,7 +1221,7 @@ struct CorrFit { if (cfgPIDConfgs.cfgPIDParticle > kCharged) fillNsigmaAfterCut(track1, cfgPIDConfgs.cfgPIDParticle); - if (!getEfficiencyCorrection_Nch(triggerWeight, track1.pt())) + if (!getEfficiencyCorrectionNch(triggerWeight, track1.pt())) continue; if (system == SameEvent) { @@ -1233,7 +1233,7 @@ struct CorrFit { if (!trackSelected(track2)) continue; - if (!getEfficiencyCorrection_Nch(associateWeight, track2.pt())) + if (!getEfficiencyCorrectionNch(associateWeight, track2.pt())) continue; if (cfgRefpTt) { @@ -1406,8 +1406,8 @@ struct CorrFit { trackCounter(tracks1, multiplicity); } - if (cfgQaCheck) { - registry.fill(HIST("Nch_corrected"), multiplicity); + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; } const auto& ft0 = collision2.foundFT0(); @@ -1460,6 +1460,10 @@ struct CorrFit { registry.fill(HIST("Nch_corrected"), multiplicity); } + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; + } + fillCorrelationsTPCFT0(tracks, ft0, collision.posZ(), SameEvent, multiplicity, kFT0C, 1.0f); } PROCESS_SWITCH(CorrFit, processSameTpcFt0c, "Process same event for TPC-FT0C correlation", false); @@ -1514,8 +1518,8 @@ struct CorrFit { trackCounter(tracks, multiplicity); } - if (cfgQaCheck) { - registry.fill(HIST("Nch_corrected"), multiplicity); + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; } fillCorrelationsTPCFT0(tracks1, ft0, collision1.posZ(), MixedEvent, multiplicity, kFT0C, eventWeight); @@ -1570,6 +1574,10 @@ struct CorrFit { registry.fill(HIST("Nch_corrected"), multiplicity); } + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; + } + fillCorrelationsFT0AFT0C(ft0, ft0, collision.posZ(), SameEvent, multiplicity, eventWeight); } PROCESS_SWITCH(CorrFit, processSameFt0aFt0c, "Process same event for FT0A-FT0C correlation", true); @@ -1626,9 +1634,10 @@ struct CorrFit { trackCounter(tracks1, multiplicity); } - if (cfgQaCheck) { - registry.fill(HIST("Nch_corrected"), multiplicity); + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; } + registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin fillCorrelationsFT0AFT0C(ft0Col1, ft0Col2, collision1.posZ(), MixedEvent, multiplicity, eventWeight); @@ -1673,6 +1682,10 @@ struct CorrFit { registry.fill(HIST("Nch_corrected"), multiplicity); } + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; + } + fillCorrelations(tracks, tracks, collision.posZ(), SameEvent, multiplicity, getMagneticField(bc.timestamp())); } PROCESS_SWITCH(CorrFit, processSameTPC, "Process same event for TPC-TPC correlation", false); @@ -1723,6 +1736,10 @@ struct CorrFit { trackCounter(tracks1, multiplicity); } + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; + } + fillCorrelations(tracks1, tracks2, collision1.posZ(), MixedEvent, multiplicity, getMagneticField(collision1.bc_as().timestamp())); } } diff --git a/PWGCF/TwoParticleCorrelations/Tasks/corrReso.cxx b/PWGCF/TwoParticleCorrelations/Tasks/corrReso.cxx new file mode 100644 index 00000000000..7a133289785 --- /dev/null +++ b/PWGCF/TwoParticleCorrelations/Tasks/corrReso.cxx @@ -0,0 +1,1660 @@ +// 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 corrReso.cxx +/// \brief Ultra long range correlation using forward FIT detectors and TPC, with focus on resonances +/// \author Thor Jensen (thor.kjaersgaard.jensen@cern.ch) + +#include "PWGCF/Core/CorrelationContainer.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/CCDB/RCTSelectionFlags.h" +#include "Common/CCDB/TriggerAliases.h" +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseITS.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::aod::rctsel; +using namespace constants::math; + +#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; + +template +auto readMatrix(Array2D const& mat, P& array, int rowStart, int rowEnd, int colStart, int colEnd) +{ + for (auto i = rowStart; i < rowEnd; ++i) { + for (auto j = colStart; j < colEnd; ++j) { + array[i][j] = mat(i, j); + } + } + + return; +} + +static constexpr float LongArrayFloat[3][20] = {{1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}, {2.1, 2.2, 2.3, -2.1, -2.2, -2.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}, {3.1, 3.2, 3.3, -3.1, -3.2, -3.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}}; +static constexpr int LongArrayInt[3][20] = {{1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1}, {2, 2, 2, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1}, {3, 3, 3, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1}}; + +struct CorrReso { + o2::aod::ITSResponse itsResponse; + Service ccdb; + + O2_DEFINE_CONFIGURABLE(cfgUseAdditionalEventCut, bool, true, "Use additional event cut on mult correlations") + O2_DEFINE_CONFIGURABLE(cfgZVtxCut, float, 10.0f, "Accepted z-vertex range") + O2_DEFINE_CONFIGURABLE(cfgUseTransverseMomentum, bool, false, "Use transverse momentum for correlation container") + O2_DEFINE_CONFIGURABLE(cfgQaCheck, bool, true, "Enable QA histograms for event selection") + O2_DEFINE_CONFIGURABLE(cfgStrictTrackCounter, bool, false, "Strict track counter for multiplicity correlation cut, counts only tracks that pass all cuts and are used in the correlation") + O2_DEFINE_CONFIGURABLE(cfgRefpTt, bool, false, "Apply upper pT cut on reference tracks") + O2_DEFINE_CONFIGURABLE(cfgRefpTMax, float, 3.0f, "maximum pT for reference tracks if cfgRefpTt is true") + O2_DEFINE_CONFIGURABLE(cfgMinMultForCorrelations, int, 0, "minimum multiplicity for correlations") + O2_DEFINE_CONFIGURABLE(cfgMaxMultForCorrelations, int, 20, "maximum multiplicity for correlations") + O2_DEFINE_CONFIGURABLE(cfgRefMultiplicity, bool, false, "Use multiplicity of reference tracks for multiplicity correlation cut instead of Nch") + Configurable> cfgRunRemoveList{"cfgRunRemoveList", std::vector{-1}, "excluded run numbers"}; + + struct : ConfigurableGroup{ + O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.2f, "minimum accepted track pT") + O2_DEFINE_CONFIGURABLE(cfgPtCutMax, float, 10.0f, "maximum accepted track pT") + O2_DEFINE_CONFIGURABLE(cfgEtaCut, float, 0.8f, "Eta cut") + O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5f, "max chi2 per TPC clusters") + O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 50.0f, "minimum TPC clusters") + O2_DEFINE_CONFIGURABLE(cfgCutTPCCrossedRows, float, 70.0f, "minimum TPC crossed rows") + O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 5.0f, "minimum ITS clusters") + O2_DEFINE_CONFIGURABLE(cfgCutDCAz, float, 2.0f, "max DCA to vertex z")} cfgTrackCuts; + + struct : ConfigurableGroup{ + O2_DEFINE_CONFIGURABLE(cfgEvSelkNoSameBunchPileup, bool, false, "rejects collisions which are associated with the same found-by-T0 bunch crossing") + O2_DEFINE_CONFIGURABLE(cfgEvSelkNoITSROFrameBorder, bool, false, "reject events at ITS ROF border") + O2_DEFINE_CONFIGURABLE(cfgEvSelkNoTimeFrameBorder, bool, false, "reject events at TF border") + O2_DEFINE_CONFIGURABLE(cfgEvSelkIsGoodZvtxFT0vsPV, bool, false, "removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference, use this cut at low multiplicities with caution") + O2_DEFINE_CONFIGURABLE(cfgEvSelkNoCollInTimeRangeStandard, bool, false, "no collisions in specified time range") + O2_DEFINE_CONFIGURABLE(cfgEvSelkIsGoodITSLayer0123, bool, true, "cut time intervals with dead ITS layers 0,1,2,3") + O2_DEFINE_CONFIGURABLE(cfgEvSelkIsGoodITSLayersAll, bool, false, "cut time intervals with dead ITS staves") + O2_DEFINE_CONFIGURABLE(cfgEvSelkNoCollInRofStandard, bool, false, "no other collisions in this Readout Frame with per-collision multiplicity above threshold") + O2_DEFINE_CONFIGURABLE(cfgEvSelkNoHighMultCollInPrevRof, bool, false, "veto an event if FT0C amplitude in previous ITS ROF is above threshold") + O2_DEFINE_CONFIGURABLE(cfgEvSelMultCorrelation, bool, false, "Multiplicity correlation cut") + O2_DEFINE_CONFIGURABLE(cfgEvSelV0AT0ACut, bool, false, "V0A T0A 5 sigma cut") + O2_DEFINE_CONFIGURABLE(cfgEvSelOccupancy, bool, false, "Occupancy cut") + O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 2000, "High cut on TPC occupancy") + O2_DEFINE_CONFIGURABLE(cfgCutOccupancyLow, int, 0, "Low cut on TPC occupancy")} cfgEventSelection; + + O2_DEFINE_CONFIGURABLE(cfgMinMixEventNum, int, 5, "Minimum number of events to mix") + O2_DEFINE_CONFIGURABLE(cfgSampleSize, double, 10, "Sample size for mixed event") + O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object") + O2_DEFINE_CONFIGURABLE(cfgEfficiencyNch, std::string, "", "CCDB path to multiplicity dependent efficiency object") + O2_DEFINE_CONFIGURABLE(cfgCentralityWeight, std::string, "", "CCDB path to centrality weight object") + O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, bool, false, "Use local efficiency object") + O2_DEFINE_CONFIGURABLE(cfgLocalEfficiencyNch, bool, false, "Use local multiplicity dependent efficiency object"); + O2_DEFINE_CONFIGURABLE(cfgUseEventWeights, bool, false, "Use event weights for mixed event") + + 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"); + O2_DEFINE_CONFIGURABLE(cfgMultT0CCutEnabled, bool, false, "Enable Global multiplicity vs T0C centrality cut") + Configurable> cfgMultT0CCutPars{"cfgMultT0CCutPars", std::vector{143.04, -4.58368, 0.0766055, -0.000727796, 2.86153e-06, 23.3108, -0.36304, 0.00437706, -4.717e-05, 1.98332e-07}, "Global multiplicity vs T0C centrality cut parameter values"}; + O2_DEFINE_CONFIGURABLE(cfgMultPVT0CCutEnabled, bool, false, "Enable PV multiplicity vs T0C centrality cut") + Configurable> cfgMultPVT0CCutPars{"cfgMultPVT0CCutPars", std::vector{195.357, -6.15194, 0.101313, -0.000955828, 3.74793e-06, 30.0326, -0.43322, 0.00476265, -5.11206e-05, 2.13613e-07}, "PV multiplicity vs T0C centrality cut parameter values"}; + O2_DEFINE_CONFIGURABLE(cfgMultMultPVHighCutFunction, std::string, "[0]+[1]*x + 5.*([2]+[3]*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultMultPVLowCutFunction, std::string, "[0]+[1]*x - 5.*([2]+[3]*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultGlobalPVCutEnabled, bool, false, "Enable global multiplicity vs PV multiplicity cut") + Configurable> cfgMultGlobalPVCutPars{"cfgMultGlobalPVCutPars", std::vector{-0.140809, 0.734344, 2.77495, 0.0165935}, "PV multiplicity vs T0C centrality cut parameter values"}; + O2_DEFINE_CONFIGURABLE(cfgMultMultV0AHighCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + 4.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultMultV0ALowCutFunction, 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"); + O2_DEFINE_CONFIGURABLE(cfgMultMultV0ACutEnabled, bool, false, "Enable global multiplicity vs V0A multiplicity cut") + Configurable> cfgMultMultV0ACutPars{"cfgMultMultV0ACutPars", std::vector{534.893, 184.344, 0.423539, -0.00331436, 5.34622e-06, 871.239, 53.3735, -0.203528, 0.000122758, 5.41027e-07}, "Global multiplicity vs V0A multiplicity cut parameter values"}; + std::vector multT0CCutPars; + std::vector multPVT0CCutPars; + std::vector multGlobalPVCutPars; + std::vector multMultV0ACutPars; + TF1* fMultPVT0CCutLow = nullptr; + TF1* fMultPVT0CCutHigh = nullptr; + TF1* fMultT0CCutLow = nullptr; + TF1* fMultT0CCutHigh = nullptr; + TF1* fMultGlobalPVCutLow = nullptr; + TF1* fMultGlobalPVCutHigh = nullptr; + TF1* fMultMultV0ACutLow = nullptr; + TF1* fMultMultV0ACutHigh = nullptr; + TF1* fT0AV0AMean = nullptr; + TF1* fT0AV0ASigma = nullptr; + } cfgFuncParas; + + struct : ConfigurableGroup { + O2_DEFINE_CONFIGURABLE(cfgUseOnlyTPC, bool, true, "Use only TPC PID for daughter selection") + + Configurable> nSigmas{"nSigmas", {LongArrayFloat[0], 6, 3, {"UpCut_pi", "UpCut_ka", "UpCut_pr", "LowCut_pi", "LowCut_ka", "LowCut_pr"}, {"TPC", "TOF", "ITS"}}, "Labeled array for n-sigma values for TPC, TOF, ITS for pions, kaons, protons (positive and negative)"}; + Configurable> cfgResoCuts{"cfgResoCuts", {LongArrayFloat[0], 12, 3, {"cos_PAs", "massMin", "massMax", "PosTrackPt", "NegTrackPt", "DCAPosToPVMin", "DCANegToPVMin", "Lifetime", "RadiusMin", "RadiusMax", "Rapidity", "ArmPodMinVal"}, {"K0", "Lambda", "Phi"}}, "Labeled array (float) for various cuts on resonances"}; + Configurable> cfgResoSwitches{"cfgResoSwitches", {LongArrayInt[0], 6, 3, {"UseCosPA", "NMassBins", "DCABetDaug", "UseProperLifetime", "UseV0Radius", "UseArmPodCut"}, {"K0", "Lambda", "Phi"}}, "Labeled array (int) for various cuts on resonances"}; + + O2_DEFINE_CONFIGURABLE(cfgUseAntiLambda, bool, true, "Use AntiLambda candidates for analysis") + O2_DEFINE_CONFIGURABLE(cfgPIDUseRejection, bool, true, "True: use exclusion exclusion criteria for PID determination, false: don't use exclusion") + + O2_DEFINE_CONFIGURABLE(cfgTpcCut, float, 3.0f, "TPC N-sigma cut for pions, kaons, protons") + O2_DEFINE_CONFIGURABLE(cfgPIDParticle, int, 0, "4 = kshort, 5 = lambda, 6 = phi, 0 for no PID") + O2_DEFINE_CONFIGURABLE(cfgUseItsPID, bool, true, "Use ITS PID for particle identification") + O2_DEFINE_CONFIGURABLE(cfgTofPtCut, float, 0.4f, "Minimum pt to use TOF N-sigma") + } cfgPIDConfigs; + + Configurable cfgCutFV0{"cfgCutFV0", 50., "FV0A threshold"}; + Configurable cfgCutFT0A{"cfgCutFT0A", 150., "FT0A threshold"}; + Configurable cfgCutFT0C{"cfgCutFT0C", 50., "FT0C threshold"}; + Configurable cfgCutZDC{"cfgCutZDC", 10., "ZDC threshold"}; + + SliceCache cache; + + ConfigurableAxis axisVertex{"axisVertex", {10, -10, 10}, "vertex axis for histograms"}; + ConfigurableAxis axisMult{"axisMult", {10, 0, 100}, "multiplicity axis for histograms"}; + ConfigurableAxis axisEta{"axisEta", {40, -1., 1.}, "eta axis for histograms"}; + ConfigurableAxis axisPhi{"axisPhi", {72, 0.0, constants::math::TwoPI}, "phi axis for histograms"}; + ConfigurableAxis axisPtFiner{"axisPtFiner", {98, 0.2, 10.0}, "pt axis for histograms"}; + ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, PIHalf * 3}, "delta phi axis for histograms"}; + ConfigurableAxis axisDeltaEtaTpcFt0a{"axisDeltaEtaTpcFt0a", {32, -5.8, -2.6}, "delta eta axis, -5.8~-2.6 for TPC-FT0A,"}; + ConfigurableAxis axisDeltaEtaTpcFt0c{"axisDeltaEtaTpcFt0c", {32, 1.2, 4.2}, "delta eta axis, 1.2~4.2 for TPC-FT0C"}; + ConfigurableAxis axisDeltaEtaFt0aFt0c{"axisDeltaEtaFt0aFt0c", {32, -1.5, 3.0}, "delta eta axis"}; + ConfigurableAxis axisPtTrigger{"axisPtTrigger", {VARIABLE_WIDTH, 0.2, 0.5, 1, 1.5, 2, 3, 4, 6, 10}, "pt trigger 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"}; + ConfigurableAxis axisSample{"axisSample", {cfgSampleSize, 0, cfgSampleSize}, "sample axis for histograms"}; + ConfigurableAxis axisNch{"axisNch", {VARIABLE_WIDTH, 0, 10, 50, 70, 100}, "multiplicity axis for correlation container"}; + + ConfigurableAxis axisNsigmaTPC{"axisNsigmaTPC", {80, -5, 5}, "nsigmaTPC axis"}; + ConfigurableAxis axisNsigmaTOF{"axisNsigmaTOF", {80, -5, 5}, "nsigmaTOF axis"}; + + ConfigurableAxis axisVertexEfficiency{"axisVertexEfficiency", {10, -10, 10}, "vertex axis for efficiency histograms"}; + ConfigurableAxis axisEtaEfficiency{"axisEtaEfficiency", {20, -1.0, 1.0}, "eta axis for efficiency histograms"}; + ConfigurableAxis axisPtEfficiency{"axisPtEfficiency", {VARIABLE_WIDTH, 0.2, 0.5, 1, 1.5, 2, 3, 4, 6, 10}, "pt axis for efficiency histograms"}; + ConfigurableAxis axisAmplitudeFt0a{"axisAmplitudeFt0a", {5000, 0, 1000}, "FT0A amplitude"}; + ConfigurableAxis axisChannelFt0aAxis{"axisChannelFt0aAxis", {96, 0.0, 96.0}, "FT0A channel"}; + + Configurable cfgGainEqPath{"cfgGainEqPath", "Analysis/EventPlane/GainEq", "CCDB path for gain equalization constants"}; + Configurable cfgCorrLevel{"cfgCorrLevel", 0, "calibration step: 0 = no corr, 1 = gain corr"}; + ConfigurableAxis cfgaxisFITamp{"cfgaxisFITamp", {1000, 0, 5000}, ""}; + AxisSpec axisFit{cfgaxisFITamp, "fit amplitude"}; + AxisSpec axisChID = {220, 0, 220}; + + Filter collisionFilter = (nabs(aod::collision::posZ) < cfgZVtxCut); + Filter trackFilter = (nabs(aod::track::eta) < cfgTrackCuts.cfgEtaCut) && (cfgTrackCuts.cfgPtCutMin < aod::track::pt) && (cfgTrackCuts.cfgPtCutMax > aod::track::pt) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)); + + using FilteredCollisions = soa::Filtered>; + using FilteredTracks = soa::Filtered>; + using V0TrackCandidate = aod::V0Datas; + + // FT0 geometry + o2::ft0::Geometry ft0Det; + static constexpr uint64_t Ft0IndexA = 96; + std::vector* offsetFT0; + std::vector cstFT0RelGain{}; + + // Corrections + TH3D* mEfficiency = nullptr; + TH1D* mEfficiencyNch = nullptr; + TH1D* mCentralityWeight = nullptr; + bool correctionsLoaded = false; + + // Define the outputs + OutputObj sameTpcFt0a{"sameEvent_TPC_FT0A"}; + OutputObj mixedTpcFt0a{"mixedEvent_TPC_FT0A"}; + OutputObj sameTpcFt0c{"sameEvent_TPC_FT0C"}; + OutputObj mixedTpcFt0c{"mixedEvent_TPC_FT0C"}; + OutputObj sameFt0aFt0c{"sameEvent_FT0A_FT0C"}; + OutputObj mixedFt0aFt0c{"mixedEvent_FT0A_FT0C"}; + + HistogramRegistry registry{"registry"}; + + // define global variables + TRandom3* gRandom = new TRandom3(); + + enum EventType { + SameEvent = 1, + MixedEvent = 3 + }; + + enum FITIndex { + kFT0A = 0, + kFT0C = 1 + }; + + enum PIDIndex { + kCharged = 0, + kPions, + kKaons, + kProtons, + kK0, + kLambda, + kPhi + }; + + enum PiKpArrayIndex { + iPionUp = 0, + iKaonUp, + iProtonUp, + iPionLow, + iKaonLow, + iProtonLow + }; + enum ResoArrayIndex { + iK0 = 0, + iLambda = 1, + iPhi = 2, + NResoParticles = 3 + }; + enum ResoParticleCuts { + kCosPA = 0, + kMassMin, + kMassMax, + kPosTrackPt, + kNegTrackPt, + kDCAPosToPVMin, + kDCANegToPVMin, + kLifeTime, + kRadiusMin, + kRadiusMax, + kRapidity, + kArmPodMinVal, + kNParticleCuts + }; + enum ResoParticleSwitches { + kUseCosPA = 0, + kMassBins, + kDCABetDaug, + kUseProperLifetime, + kUseV0Radius, + kUseArmPodCut, + kNParticleSwitches + }; + enum DetectorType { + kTPC = 0, + kTOF, + kITS + }; + + std::array, 12> resoCutVals; + std::array, 7> resoSwitchVals; + std::array tofNsigmaCut; + std::array itsNsigmaCut; + std::array tpcNsigmaCut; + + RCTFlagsChecker rctChecker{"CBT"}; + + void init(InitContext&) + { + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + auto now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); + ccdb->setCreatedNotAfter(now); + + LOGF(info, "Starting init"); + + if (cfgPIDConfigs.cfgPIDParticle < kK0) { + LOGF(fatal, "The input number does not correspond to any particle"); + return; + } + + readMatrix(cfgPIDConfigs.cfgResoCuts->getData(), resoCutVals, kCosPA, kNParticleCuts, iK0, NResoParticles); + readMatrix(cfgPIDConfigs.cfgResoSwitches->getData(), resoSwitchVals, kUseCosPA, kNParticleSwitches, iK0, NResoParticles); + + tpcNsigmaCut[iPionUp] = cfgPIDConfigs.nSigmas->getData()[iPionUp][kTPC]; + tpcNsigmaCut[iKaonUp] = cfgPIDConfigs.nSigmas->getData()[iKaonUp][kTPC]; + tpcNsigmaCut[iProtonUp] = cfgPIDConfigs.nSigmas->getData()[iProtonUp][kTPC]; + tpcNsigmaCut[iPionLow] = cfgPIDConfigs.nSigmas->getData()[iPionLow][kTPC]; + tpcNsigmaCut[iKaonLow] = cfgPIDConfigs.nSigmas->getData()[iKaonLow][kTPC]; + tpcNsigmaCut[iProtonLow] = cfgPIDConfigs.nSigmas->getData()[iProtonLow][kTPC]; + + tofNsigmaCut[iPionUp] = cfgPIDConfigs.nSigmas->getData()[iPionUp][kTOF]; + tofNsigmaCut[iKaonUp] = cfgPIDConfigs.nSigmas->getData()[iKaonUp][kTOF]; + tofNsigmaCut[iProtonUp] = cfgPIDConfigs.nSigmas->getData()[iProtonUp][kTOF]; + tofNsigmaCut[iPionLow] = cfgPIDConfigs.nSigmas->getData()[iPionLow][kTOF]; + tofNsigmaCut[iKaonLow] = cfgPIDConfigs.nSigmas->getData()[iKaonLow][kTOF]; + tofNsigmaCut[iProtonLow] = cfgPIDConfigs.nSigmas->getData()[iProtonLow][kTOF]; + + itsNsigmaCut[iPionUp] = cfgPIDConfigs.nSigmas->getData()[iPionUp][kITS]; + itsNsigmaCut[iKaonUp] = cfgPIDConfigs.nSigmas->getData()[iKaonUp][kITS]; + itsNsigmaCut[iProtonUp] = cfgPIDConfigs.nSigmas->getData()[iProtonUp][kITS]; + itsNsigmaCut[iPionLow] = cfgPIDConfigs.nSigmas->getData()[iPionLow][kITS]; + itsNsigmaCut[iKaonLow] = cfgPIDConfigs.nSigmas->getData()[iKaonLow][kITS]; + itsNsigmaCut[iProtonLow] = cfgPIDConfigs.nSigmas->getData()[iProtonLow][kITS]; + + // Creating mass axis depending on particle - 4 = kshort, 5 = lambda, 6 = phi + AxisSpec axisInvMass = {10, 0, 1, "mass"}; + if (cfgPIDConfigs.cfgPIDParticle == kK0) + axisInvMass = {resoSwitchVals[kMassBins][iK0], resoCutVals[kMassMin][iK0], resoCutVals[kMassMax][iK0], "M_{#pi^{+}#pi^{-}} (GeV/c^{2})"}; + if (cfgPIDConfigs.cfgPIDParticle == kLambda) + axisInvMass = {resoSwitchVals[kMassBins][iLambda], resoCutVals[kMassMin][iLambda], resoCutVals[kMassMax][iLambda], "M_{p#pi^{-}} (GeV/c^{2})"}; + if (cfgPIDConfigs.cfgPIDParticle == kPhi) + axisInvMass = {resoSwitchVals[kMassBins][iPhi], resoCutVals[kMassMin][iPhi], resoCutVals[kMassMax][iPhi], "M_{K^{+}K^{-}} (GeV/c^{2})"}; + + if (doprocessSameTpcFt0a || doprocessSameTpcFt0c) { + if (cfgPIDConfigs.cfgPIDParticle == kK0) { // For K0 + registry.add("PiPlusTPC_K0", "", {HistType::kTH2D, {{axisPtFiner, axisNsigmaTPC}}}); + registry.add("PiMinusTPC_K0", "", {HistType::kTH2D, {{axisPtFiner, axisNsigmaTPC}}}); + registry.add("PiPlusTOF_K0", "", {HistType::kTH2D, {{axisPtFiner, axisNsigmaTOF}}}); + registry.add("PiMinusTOF_K0", "", {HistType::kTH2D, {{axisPtFiner, axisNsigmaTOF}}}); + registry.add("hK0Phi", "", {HistType::kTH1D, {axisPhi}}); + registry.add("hK0Eta", "", {HistType::kTH1D, {axisEta}}); + + registry.add("hK0Count", "Number of K0;; Count", {HistType::kTH1D, {{11, 0, 11}}}); + registry.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(1, "K0 candidates"); + registry.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(2, "Daughter pt"); + registry.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(3, "Mass cut"); + registry.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(4, "Rapidity cut"); + registry.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(5, "DCA to PV"); + registry.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(6, "DCA between daughters"); + registry.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(7, "V0radius"); + registry.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(8, "CosPA"); + registry.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(9, "Proper lifetime"); + registry.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(10, "ArmenterosPod"); + registry.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(11, "Daughter track selection"); + } + if (cfgPIDConfigs.cfgPIDParticle == kLambda) { // For Lambda + registry.add("PrPlusTPC_L", "", {HistType::kTH2D, {{axisPtFiner, axisNsigmaTPC}}}); + registry.add("PiMinusTPC_L", "", {HistType::kTH2D, {{axisPtFiner, axisNsigmaTPC}}}); + registry.add("PrPlusTOF_L", "", {HistType::kTH2D, {{axisPtFiner, axisNsigmaTOF}}}); + registry.add("PiMinusTOF_L", "", {HistType::kTH2D, {{axisPtFiner, axisNsigmaTOF}}}); + registry.add("hLambdaPhi", "", {HistType::kTH1D, {axisPhi}}); + registry.add("hLambdaEta", "", {HistType::kTH1D, {axisEta}}); + + registry.add("hLambdaCount", "Number of Lambda;; Count", {HistType::kTH1D, {{10, 0, 10}}}); + registry.get(HIST("hLambdaCount"))->GetXaxis()->SetBinLabel(1, "Lambda candidates"); + registry.get(HIST("hLambdaCount"))->GetXaxis()->SetBinLabel(2, "Daughter pt"); + registry.get(HIST("hLambdaCount"))->GetXaxis()->SetBinLabel(3, "Mass cut"); + registry.get(HIST("hLambdaCount"))->GetXaxis()->SetBinLabel(4, "Rapidity cut"); + registry.get(HIST("hLambdaCount"))->GetXaxis()->SetBinLabel(5, "DCA to PV"); + registry.get(HIST("hLambdaCount"))->GetXaxis()->SetBinLabel(6, "DCA between daughters"); + registry.get(HIST("hLambdaCount"))->GetXaxis()->SetBinLabel(7, "V0radius"); + registry.get(HIST("hLambdaCount"))->GetXaxis()->SetBinLabel(8, "CosPA"); + registry.get(HIST("hLambdaCount"))->GetXaxis()->SetBinLabel(9, "Proper lifetime"); + registry.get(HIST("hLambdaCount"))->GetXaxis()->SetBinLabel(10, "Daughter track selection"); + } + } + if (doprocessSameFt0aFt0c || doprocessSameTpcFt0a || doprocessSameTpcFt0c) { + registry.add("hEventCountRct", "Number of Event;; Count", {HistType::kTH1D, {{2, 0, 2}}}); + registry.get(HIST("hEventCountRct"))->GetXaxis()->SetBinLabel(1, "rct fail"); + registry.get(HIST("hEventCountRct"))->GetXaxis()->SetBinLabel(2, "rct pass"); + registry.add("hEventCountSpecific", "Number of Event;; Count", {HistType::kTH1D, {{13, 0, 13}}}); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(1, "after sel8"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(2, "kNoSameBunchPileup"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(3, "kNoITSROFrameBorder"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(4, "kNoTimeFrameBorder"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(5, "kIsGoodZvtxFT0vsPV"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(6, "kNoCollInTimeRangeStandard"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(7, "kIsGoodITSLayer0123"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(8, "kIsGoodITSLayersAll"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(9, "kNoCollInRofStandard"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(10, "kNoHighMultCollInPrevRof"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(11, "occupancy"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(12, "MultCorrelation"); + registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(13, "cfgEvSelV0AT0ACut"); + } + + if ((doprocessSameFt0aFt0c || doprocessSameTpcFt0a || doprocessSameTpcFt0c) && cfgQaCheck) { + registry.add("hPassedEventSelection", "Number of Event;; Count", {HistType::kTH1D, {{12, 0, 12}}}); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(1, "all tracks"); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(2, "after sel8"); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(3, "kNoSameBunchPileup"); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(4, "kNoTimeFrameBorder"); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(5, "kNoITSROFrameBorder"); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(6, "kIsGoodZvtxFT0vsPV"); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(7, "kNoCollInTimeRangeStandard"); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(8, "kIsGoodITSLayer0123"); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(9, "kIsGoodITSLayersAll"); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(10, "kNoCollInRofStandard"); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(11, "kNoHighMultCollInPrevRof"); + registry.get(HIST("hPassedEventSelection"))->GetXaxis()->SetBinLabel(12, "occupancy"); + } + + if (doprocessSameTpcFt0a || doprocessSameTpcFt0c || doprocessSameFt0aFt0c) { + registry.add("Phi", "Phi", {HistType::kTH1D, {axisPhi}}); + registry.add("Eta", "Eta", {HistType::kTH1D, {axisEta}}); + registry.add("EtaCorrected", "EtaCorrected", {HistType::kTH1D, {axisEta}}); + registry.add("pTFiner", "pTFiner", {HistType::kTH1D, {axisPtFiner}}); + registry.add("pTFinerCorrected", "pTFinerCorrected", {HistType::kTH1D, {axisPtFiner}}); + registry.add("Nch", "N_{ch}", {HistType::kTH1D, {axisMult}}); + registry.add("zVtx", "zVtx", {HistType::kTH1D, {axisVertex}}); + + registry.add("FT0Amp", "", {HistType::kTH2F, {axisChID, axisFit}}); + registry.add("FT0AmpCorrect", "", {HistType::kTH2F, {axisChID, axisFit}}); + } // end of single particle distribution histograms + + if (cfgQaCheck) { + registry.add("Nch_corrected", "N_{ch} corrected", {HistType::kTH1D, {axisMult}}); + } + + if (doprocessSameTpcFt0a) { + registry.add("deltaEta_deltaPhi_same_TPC_FT0A", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFt0a}}); // check to see the delta eta and delta phi distribution + registry.add("deltaEta_deltaPhi_mixed_TPC_FT0A", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFt0a}}); + registry.add("Assoc_amp_same_TPC_FT0A", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); + registry.add("Assoc_amp_mixed_TPC_FT0A", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); + registry.add("Trig_hist_TPC_FT0A", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger, axisInvMass}}}); + } + if (doprocessSameTpcFt0c) { + registry.add("deltaEta_deltaPhi_same_TPC_FT0C", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFt0c}}); // check to see the delta eta and delta phi distribution + registry.add("deltaEta_deltaPhi_mixed_TPC_FT0C", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFt0c}}); + registry.add("Assoc_amp_same_TPC_FT0C", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); + registry.add("Assoc_amp_mixed_TPC_FT0C", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); + registry.add("Trig_hist_TPC_FT0C", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger, axisInvMass}}}); + } + if (doprocessSameFt0aFt0c) { + registry.add("deltaEta_deltaPhi_same_FT0A_FT0C", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaFt0aFt0c}}); // check to see the delta eta and delta phi distribution + registry.add("deltaEta_deltaPhi_mixed_FT0A_FT0C", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaFt0aFt0c}}); + registry.add("Trig_hist_FT0A_FT0C", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); + } + + registry.add("eventcount", "bin", {HistType::kTH1F, {{4, 0, 4, "bin"}}}); // histogram to see how many events are in the same and mixed event + + LOGF(info, "Initializing correlation container"); + + // Initialize Nch-related histograms and containers + + std::vector corrAxisTpcFt0a = {{axisSample, "Sample"}, + {axisVertex, "z-vtx (cm)"}, + {axisPtTrigger, "p_{T} (GeV/c)"}, + {axisInvMass}, + {axisDeltaPhi, "#Delta#varphi (rad)"}, + {axisDeltaEtaTpcFt0a, "#Delta#eta"}}; + + std::vector corrAxisTpcFt0c = {{axisSample, "Sample"}, + {axisVertex, "z-vtx (cm)"}, + {axisPtTrigger, "p_{T} (GeV/c)"}, + {axisInvMass}, + {axisDeltaPhi, "#Delta#varphi (rad)"}, + {axisDeltaEtaTpcFt0c, "#Delta#eta"}}; + + std::vector corrAxisFt0aFt0c = {{axisSample, "Sample"}, + {axisVertex, "z-vtx (cm)"}, + {axisPtTrigger, "p_{T} (GeV/c)"}, + {axisNch, "N_{ch}"}, + {axisDeltaPhi, "#Delta#varphi (rad)"}, + {axisDeltaEtaFt0aFt0c, "#Delta#eta"}}; + + std::vector effAxis = { + {axisEtaEfficiency, "#eta"}, + {axisPtEfficiency, "p_{T} (GeV/c)"}, + {axisVertexEfficiency, "z-vtx (cm)"}, + }; + std::vector userAxis; + + if (doprocessSameTpcFt0a) { + sameTpcFt0a.setObject(new CorrelationContainer("sameEvent_TPC_FT0A", "sameEvent_TPC_FT0A", corrAxisTpcFt0a, effAxis, userAxis)); + mixedTpcFt0a.setObject(new CorrelationContainer("mixedEvent_TPC_FT0A", "mixedEvent_TPC_FT0A", corrAxisTpcFt0a, effAxis, userAxis)); + } + if (doprocessSameTpcFt0c) { + sameTpcFt0c.setObject(new CorrelationContainer("sameEvent_TPC_FT0C", "sameEvent_TPC_FT0C", corrAxisTpcFt0c, effAxis, userAxis)); + mixedTpcFt0c.setObject(new CorrelationContainer("mixedEvent_TPC_FT0C", "mixedEvent_TPC_FT0C", corrAxisTpcFt0c, effAxis, userAxis)); + } + if (doprocessSameFt0aFt0c) { + sameFt0aFt0c.setObject(new CorrelationContainer("sameEvent_FT0A_FT0C", "sameEvent_FT0A_FT0C", corrAxisFt0aFt0c, effAxis, userAxis)); + mixedFt0aFt0c.setObject(new CorrelationContainer("mixedEvent_FT0A_FT0C", "mixedEvent_FT0A_FT0C", corrAxisFt0aFt0c, effAxis, userAxis)); + } + + LOGF(info, "End of init"); + } + template + bool eventRct(TCollision const& collision, const bool fillCounter) + { + if (!rctChecker(collision)) { + if (fillCounter) + registry.fill(HIST("hEventCountRct"), 0.5); + + return 0; + } + if (fillCounter) + registry.fill(HIST("hEventCountRct"), 1.5); + + return 1; + } + + template + bool eventSelected(TCollision const& collision, const int multTrk, const bool fillCounter) + { + registry.fill(HIST("hEventCountSpecific"), 0.5); + if (cfgEventSelection.cfgEvSelkNoSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + // rejects collisions which are associated with the same "found-by-T0" bunch crossing + // https://indico.cern.ch/event/1396220/#1-event-selection-with-its-rof + return 0; + } + if (fillCounter && cfgEventSelection.cfgEvSelkNoSameBunchPileup) + registry.fill(HIST("hEventCountSpecific"), 1.5); + if (cfgEventSelection.cfgEvSelkNoITSROFrameBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { + return 0; + } + if (fillCounter && cfgEventSelection.cfgEvSelkNoITSROFrameBorder) + registry.fill(HIST("hEventCountSpecific"), 2.5); + if (cfgEventSelection.cfgEvSelkNoTimeFrameBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) { + return 0; + } + if (fillCounter && cfgEventSelection.cfgEvSelkNoTimeFrameBorder) + registry.fill(HIST("hEventCountSpecific"), 3.5); + if (cfgEventSelection.cfgEvSelkIsGoodZvtxFT0vsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + // removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference + // use this cut at low multiplicities with caution + return 0; + } + if (fillCounter && cfgEventSelection.cfgEvSelkIsGoodZvtxFT0vsPV) + registry.fill(HIST("hEventCountSpecific"), 4.5); + if (cfgEventSelection.cfgEvSelkNoCollInTimeRangeStandard && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { + // no collisions in specified time range + return 0; + } + if (fillCounter && cfgEventSelection.cfgEvSelkNoCollInTimeRangeStandard) + registry.fill(HIST("hEventCountSpecific"), 5.5); + + if (cfgEventSelection.cfgEvSelkIsGoodITSLayer0123 && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayer0123)) { + // from Jan 9 2025 AOT meeting + // cut time intervals with dead ITS staves + return 0; + } + if (fillCounter && cfgEventSelection.cfgEvSelkIsGoodITSLayer0123) + registry.fill(HIST("hEventCountSpecific"), 6.5); + + if (cfgEventSelection.cfgEvSelkIsGoodITSLayersAll && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { + // from Jan 9 2025 AOT meeting + // cut time intervals with dead ITS staves + return 0; + } + + if (fillCounter && cfgEventSelection.cfgEvSelkIsGoodITSLayersAll) + registry.fill(HIST("hEventCountSpecific"), 7.5); + + if (cfgEventSelection.cfgEvSelkNoCollInRofStandard && !collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard)) { + // no other collisions in this Readout Frame with per-collision multiplicity above threshold + return 0; + } + if (fillCounter && cfgEventSelection.cfgEvSelkNoCollInRofStandard) + registry.fill(HIST("hEventCountSpecific"), 8.5); + if (cfgEventSelection.cfgEvSelkNoHighMultCollInPrevRof && !collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof)) { + // veto an event if FT0C amplitude in previous ITS ROF is above threshold + return 0; + } + if (fillCounter && cfgEventSelection.cfgEvSelkNoHighMultCollInPrevRof) + registry.fill(HIST("hEventCountSpecific"), 9.5); + auto occupancy = collision.trackOccupancyInTimeRange(); + if (cfgEventSelection.cfgEvSelOccupancy && (occupancy < cfgEventSelection.cfgCutOccupancyLow || occupancy > cfgEventSelection.cfgCutOccupancyHigh)) + return 0; + if (fillCounter && cfgEventSelection.cfgEvSelOccupancy) + registry.fill(HIST("hEventCountSpecific"), 10.5); + + auto multNTracksPV = collision.multNTracksPV(); + + if (cfgFuncParas.cfgMultGlobalPVCutEnabled) { + if (multTrk < cfgFuncParas.fMultGlobalPVCutLow->Eval(multNTracksPV)) + return 0; + if (multTrk > cfgFuncParas.fMultGlobalPVCutHigh->Eval(multNTracksPV)) + return 0; + } + if (cfgFuncParas.cfgMultMultV0ACutEnabled) { + if (collision.multFV0A() < cfgFuncParas.fMultMultV0ACutLow->Eval(multTrk)) + return 0; + if (collision.multFV0A() > cfgFuncParas.fMultMultV0ACutHigh->Eval(multTrk)) + return 0; + } + + if (fillCounter && cfgEventSelection.cfgEvSelMultCorrelation) + registry.fill(HIST("hEventCountSpecific"), 11.5); + + // V0A T0A 5 sigma cut + float sigma = 5.0; + if (cfgEventSelection.cfgEvSelV0AT0ACut && (std::fabs(collision.multFV0A() - cfgFuncParas.fT0AV0AMean->Eval(collision.multFT0A())) > sigma * cfgFuncParas.fT0AV0ASigma->Eval(collision.multFT0A()))) + return 0; + if (fillCounter && cfgEventSelection.cfgEvSelV0AT0ACut) + registry.fill(HIST("hEventCountSpecific"), 12.5); + + return 1; + } + + template + void eventSelectedIndividually(TCollision const& collision) + { + + registry.fill(HIST("hPassedEventSelection"), 0.5); + + if (collision.sel8()) { + registry.fill(HIST("hPassedEventSelection"), 1.5); + } + + if (collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + registry.fill(HIST("hPassedEventSelection"), 2.5); + } + + if (collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) { + registry.fill(HIST("hPassedEventSelection"), 3.5); + } + + if (collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { + registry.fill(HIST("hPassedEventSelection"), 4.5); + } + + if (collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + registry.fill(HIST("hPassedEventSelection"), 5.5); + } + + if (collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { + registry.fill(HIST("hPassedEventSelection"), 6.5); + } + + if (collision.selection_bit(o2::aod::evsel::kIsGoodITSLayer0123)) { + registry.fill(HIST("hPassedEventSelection"), 7.5); + } + + if (collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { + registry.fill(HIST("hPassedEventSelection"), 8.5); + } + + if (collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard)) { + registry.fill(HIST("hPassedEventSelection"), 9.5); + } + + if (collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof)) { + registry.fill(HIST("hPassedEventSelection"), 10.5); + } + + auto occupancy = collision.trackOccupancyInTimeRange(); + if (cfgEventSelection.cfgEvSelOccupancy && (occupancy < cfgEventSelection.cfgCutOccupancyLow || occupancy > cfgEventSelection.cfgCutOccupancyHigh)) { + registry.fill(HIST("hPassedEventSelection"), 11.5); + } + } + + double getPhiFT0(uint64_t chno, int i) + { + // offsetFT0[0]: FT0A, offsetFT0[1]: FT0C + if (i > 1 || i < 0) { + LOGF(fatal, "kFIT Index %d out of range", i); + } + + ft0Det.calculateChannelCenter(); + auto chPos = ft0Det.getChannelCenter(chno); + return RecoDecay::phi(chPos.X() + (*offsetFT0)[i].getX(), chPos.Y() + (*offsetFT0)[i].getY()); + } + + double getEtaFT0(uint64_t chno, int i) + { + // offsetFT0[0]: FT0A, offsetFT0[1]: FT0C + if (i > 1 || i < 0) { + LOGF(fatal, "kFIT Index %d out of range", i); + } + ft0Det.calculateChannelCenter(); + auto chPos = ft0Det.getChannelCenter(chno); + auto x = chPos.X() + (*offsetFT0)[i].getX(); + auto y = chPos.Y() + (*offsetFT0)[i].getY(); + auto z = chPos.Z() + (*offsetFT0)[i].getZ(); + if (chno >= Ft0IndexA) { + z = -z; + } + auto r = std::sqrt(x * x + y * y); + auto theta = std::atan2(r, z); + return -std::log(std::tan(0.5 * theta)); + } + + void loadAlignParam(uint64_t timestamp) + { + offsetFT0 = ccdb->getForTimeStamp>("FT0/Calib/Align", timestamp); + if (offsetFT0 == nullptr) { + LOGF(fatal, "Could not load FT0/Calib/Align for timestamp %d", timestamp); + } + } + + template + bool trackSelected(TTrack const& track) + { + return ((track.tpcNClsFound() >= cfgTrackCuts.cfgCutTPCclu) && (track.tpcNClsCrossedRows() >= cfgTrackCuts.cfgCutTPCCrossedRows) && (track.itsNCls() >= cfgTrackCuts.cfgCutITSclu) && (track.tpcChi2NCl() < cfgTrackCuts.cfgCutChi2prTPCcls) && (track.dcaZ() < cfgTrackCuts.cfgCutDCAz)); + } + + void loadGain(aod::BCsWithTimestamps::iterator const& bc) + { + cstFT0RelGain.clear(); + cstFT0RelGain = {}; + std::string fullPath; + + auto timestamp = bc.timestamp(); + constexpr int ChannelsFT0 = 208; + if (cfgCorrLevel == 0) { + for (auto i{0u}; i < ChannelsFT0; i++) { + cstFT0RelGain.push_back(1.); + } + } else { + fullPath = cfgGainEqPath; + fullPath += "/FT0"; + const auto objft0Gain = ccdb->getForTimeStamp>(fullPath, timestamp); + if (!objft0Gain) { + for (auto i{0u}; i < ChannelsFT0; i++) { + cstFT0RelGain.push_back(1.); + } + } else { + cstFT0RelGain = *(objft0Gain); + } + } + } + + template + void getChannel(TFT0s const& ft0, std::size_t const& iCh, int& id, float& ampl, int fitType, int system) + { + if (fitType == kFT0C) { + id = ft0.channelC()[iCh]; + id = id + Ft0IndexA; + ampl = ft0.amplitudeC()[iCh]; + if (system == SameEvent) + registry.fill(HIST("FT0Amp"), id, ampl); + ampl = ampl / cstFT0RelGain[id]; + if (system == SameEvent) { + registry.fill(HIST("FT0AmpCorrect"), id, ampl); + } + } else if (fitType == kFT0A) { + id = ft0.channelA()[iCh]; + ampl = ft0.amplitudeA()[iCh]; + if (system == SameEvent) + registry.fill(HIST("FT0Amp"), id, ampl); + ampl = ampl / cstFT0RelGain[id]; + if (system == SameEvent) { + registry.fill(HIST("FT0AmpCorrect"), id, ampl); + } + } else { + LOGF(fatal, "Cor Index %d out of range", fitType); + } + } + + template + int getNsigmaPID(TTrack track) + { + // Computing Nsigma arrays for pion, kaon, and protons + std::array nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()}; + std::array nSigmaTOF = {track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()}; + std::array nSigmaITS = {itsResponse.nSigmaITS(track), itsResponse.nSigmaITS(track), itsResponse.nSigmaITS(track)}; + int pid = -1; // -1 = not identified, 1 = pion, 2 = kaon, 3 = proton + + std::array nSigmaToUse = cfgPIDConfigs.cfgUseItsPID ? nSigmaITS : nSigmaTPC; // Choose which nSigma to use: TPC or ITS + std::array detectorNsigmaCut = cfgPIDConfigs.cfgUseItsPID ? itsNsigmaCut : tpcNsigmaCut; // Choose which nSigma to use: TPC or ITS + + bool isPion, isKaon, isProton; + bool isDetectedPion = nSigmaToUse[iPionUp] < detectorNsigmaCut[iPionUp] && nSigmaToUse[iPionUp] > detectorNsigmaCut[iPionLow]; + bool isDetectedKaon = nSigmaToUse[iKaonUp] < detectorNsigmaCut[iKaonUp] && nSigmaToUse[iKaonUp] > detectorNsigmaCut[iKaonLow]; + bool isDetectedProton = nSigmaToUse[iProtonUp] < detectorNsigmaCut[iProtonUp] && nSigmaToUse[iProtonUp] > detectorNsigmaCut[iProtonLow]; + + bool isTofPion = nSigmaTOF[iPionUp] < tofNsigmaCut[iPionUp] && nSigmaTOF[iPionUp] > tofNsigmaCut[iPionLow]; + bool isTofKaon = nSigmaTOF[iKaonUp] < tofNsigmaCut[iKaonUp] && nSigmaTOF[iKaonUp] > tofNsigmaCut[iKaonLow]; + bool isTofProton = nSigmaTOF[iProtonUp] < tofNsigmaCut[iProtonUp] && nSigmaTOF[iProtonUp] > tofNsigmaCut[iProtonLow]; + + if (track.pt() > cfgPIDConfigs.cfgTofPtCut && !track.hasTOF()) { + return -1; + } else if (track.pt() > cfgPIDConfigs.cfgTofPtCut && track.hasTOF()) { + isPion = isTofPion && isDetectedPion; + isKaon = isTofKaon && isDetectedKaon; + isProton = isTofProton && isDetectedProton; + } else { + isPion = isDetectedPion; + isKaon = isDetectedKaon; + isProton = isDetectedProton; + } + + if (cfgPIDConfigs.cfgPIDUseRejection && ((isPion && isKaon) || (isPion && isProton) || (isKaon && isProton))) { + return -1; // more than one particle satisfy the criteria + } + + if (isPion) { + pid = kPions; + } else if (isKaon) { + pid = kKaons; + } else if (isProton) { + pid = kProtons; + } else { + return -1; // no particle satisfies the criteria + } + + return pid; // -1 = not identified, 1 = pion, 2 = kaon, 3 = proton + } + + template + bool selectionV0Daughter(TTrack const& track, int pid) + { + if (!(track.itsNCls() > cfgTrackCuts.cfgCutITSclu)) + return false; + if (!track.hasTPC()) + return false; + if (!(track.tpcNClsFound() >= cfgTrackCuts.cfgCutTPCclu)) + return false; + if (!(track.tpcNClsCrossedRows() >= cfgTrackCuts.cfgCutTPCCrossedRows)) + return false; + if (!(track.dcaZ() < cfgTrackCuts.cfgCutDCAz)) + return false; + + if (cfgPIDConfigs.cfgUseOnlyTPC) { + if (pid == kPions && std::abs(track.tpcNSigmaPi()) > cfgPIDConfigs.cfgTpcCut) + return false; + if (pid == kKaons && std::abs(track.tpcNSigmaKa()) > cfgPIDConfigs.cfgTpcCut) + return false; + if (pid == kProtons && std::abs(track.tpcNSigmaPr()) > cfgPIDConfigs.cfgTpcCut) + return false; + } else { + int partIndex = getNsigmaPID(track); + int pidIndex = partIndex; // 1 = pion, 2 = kaon, 3 = proton + if (pidIndex != pid) + return false; + } + + return true; + } + + void loadCorrection(uint64_t timestamp) + { + if (correctionsLoaded) { + return; + } + if (cfgEfficiency.value.empty() == false) { + if (cfgLocalEfficiency) { + TFile* fEfficiencyTrigger = TFile::Open(cfgEfficiency.value.c_str(), "READ"); + mEfficiency = reinterpret_cast(fEfficiencyTrigger->Get("ccdb_object")); + + } else { + mEfficiency = ccdb->getForTimeStamp(cfgEfficiency, timestamp); + } + if (mEfficiency == nullptr) { + LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgEfficiency.value.c_str()); + } + LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)mEfficiency); + } + if (cfgEfficiencyNch.value.empty() == false) { + if (cfgLocalEfficiencyNch) { + TFile* fEfficiencyTrigger = TFile::Open(cfgEfficiencyNch.value.c_str(), "READ"); + mEfficiencyNch = reinterpret_cast(fEfficiencyTrigger->Get("ccdb_object")); + + } else { + mEfficiencyNch = ccdb->getForTimeStamp(cfgEfficiencyNch, timestamp); + } + if (!mEfficiencyNch) { + LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgEfficiencyNch.value.c_str()); + } + LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiencyNch.value.c_str(), (void*)mEfficiencyNch); + } + if (cfgCentralityWeight.value.empty() == false) { + mCentralityWeight = ccdb->getForTimeStamp(cfgCentralityWeight, timestamp); + if (mCentralityWeight == nullptr) { + LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgCentralityWeight.value.c_str()); + } + LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgCentralityWeight.value.c_str(), (void*)mCentralityWeight); + } + correctionsLoaded = true; + } + + bool getEfficiencyCorrectionNch(float& weightNch, float pt) + { + float effNch = 1.; + if (mEfficiencyNch) { + + int ptBin = mEfficiencyNch->FindBin(pt); + effNch = mEfficiencyNch->GetBinContent(ptBin); + + } else { + effNch = 1.0; + } + if (effNch == 0) + return false; + weightNch = 1. / effNch; + return true; + } + + bool getEfficiencyCorrection(float& weight, float pt, float eta, float vertex) + { + float eff = 1.; + if (mEfficiency) { + + int etaBin = mEfficiency->GetXaxis()->FindBin(eta); // use the eta bin corresponding to eta=0 for the trigger particle efficiency + int ptBin = mEfficiency->GetYaxis()->FindBin(pt); + int vertexBin = mEfficiency->GetZaxis()->FindBin(vertex); // use the vertex bin corresponding to z=0 for the trigger particle efficiency + eff = mEfficiency->GetBinContent(etaBin, ptBin, vertexBin); + + } else { + eff = 1.0; + } + if (eff == 0) + return false; + weight = 1. / eff; + return true; + } + + template + void trackCounter(TTracks tracks, double& multiplicity) // function to count the number of tracks in the event and fill the histogram + { + double nTracksCorrected = 0; + float weightNch = 1.0f; + for (auto const& track : tracks) { + + if (cfgRefMultiplicity) { + if (track.pt() > cfgRefpTMax) + continue; + } + + if (!getEfficiencyCorrectionNch(weightNch, track.pt())) { + continue; + } + + nTracksCorrected += weightNch; + } + multiplicity = nTracksCorrected; + } + + template + void fillYield(TCollision collision, TTracks tracks) // function to fill the yield and etaphi histograms. + { + + float weff1 = 1.0; + float zvtx = collision.posZ(); + + for (auto const& track1 : tracks) { + + if (!trackSelected(track1)) { + continue; + } + + if (!getEfficiencyCorrection(weff1, track1.eta(), track1.pt(), zvtx)) { + continue; + } + + registry.fill(HIST("Phi"), RecoDecay::constrainAngle(track1.phi(), 0.0)); + registry.fill(HIST("Eta"), track1.eta()); + registry.fill(HIST("EtaCorrected"), track1.eta(), weff1); + registry.fill(HIST("pTFiner"), track1.pt()); + registry.fill(HIST("pTFinerCorrected"), track1.pt(), weff1); + } + } + + template + bool isSelectedK0(V0 const& candidate, float posZ, float posY, float posX) + { + double mk0 = candidate.mK0Short(); + + // separate the positive and negative V0 daughters + auto postrack = candidate.template posTrack_as(); + auto negtrack = candidate.template negTrack_as(); + + registry.fill(HIST("hK0Count"), 0.5); + if (postrack.pt() < resoCutVals[kPosTrackPt][iK0] || negtrack.pt() < resoCutVals[kNegTrackPt][iK0]) + return false; + registry.fill(HIST("hK0Count"), 1.5); + if (mk0 < resoCutVals[kMassMin][iK0] && mk0 > resoCutVals[kMassMax][iK0]) + return false; + registry.fill(HIST("hK0Count"), 2.5); + // Rapidity correction + if (candidate.yK0Short() > resoCutVals[kRapidity][iK0]) + return false; + registry.fill(HIST("hK0Count"), 3.5); + // DCA cuts for K0short + if (std::abs(candidate.dcapostopv()) < resoCutVals[kDCAPosToPVMin][iK0] || std::abs(candidate.dcanegtopv()) < resoCutVals[kDCANegToPVMin][iK0]) + return false; + registry.fill(HIST("hK0Count"), 4.5); + if (std::abs(candidate.dcaV0daughters()) > resoSwitchVals[kDCABetDaug][iK0]) + return false; + registry.fill(HIST("hK0Count"), 5.5); + // v0 radius cuts + if (resoSwitchVals[kUseV0Radius][iK0] && (candidate.v0radius() < resoCutVals[kRadiusMin][iK0] || candidate.v0radius() > resoCutVals[kRadiusMax][iK0])) + return false; + registry.fill(HIST("hK0Count"), 6.5); + // cosine pointing angle cuts + if (candidate.v0cosPA() < resoCutVals[kCosPA][iK0]) + return false; + registry.fill(HIST("hK0Count"), 7.5); + // Proper lifetime + float cTauK0 = candidate.distovertotmom(posX, posY, posZ) * massK0Short; + if (resoSwitchVals[kUseProperLifetime][iK0] && std::abs(cTauK0) > resoCutVals[kLifeTime][iK0]) + return false; + registry.fill(HIST("hK0Count"), 8.5); + // ArmenterosPodolanskiCut + if (resoSwitchVals[kUseArmPodCut][iK0] && (candidate.qtarm() / std::abs(candidate.alpha())) < resoCutVals[kArmPodMinVal][iK0]) + return false; + registry.fill(HIST("hK0Count"), 9.5); + // Selection on V0 daughters + if (!selectionV0Daughter(postrack, kPions) || !selectionV0Daughter(negtrack, kPions)) + return false; + registry.fill(HIST("hK0Count"), 10.5); + + registry.fill(HIST("hK0Phi"), candidate.phi()); + registry.fill(HIST("hK0Eta"), candidate.eta()); + registry.fill(HIST("PiPlusTPC_K0"), postrack.pt(), postrack.tpcNSigmaPi()); + registry.fill(HIST("PiPlusTOF_K0"), postrack.pt(), postrack.tofNSigmaPi()); + registry.fill(HIST("PiMinusTPC_K0"), negtrack.pt(), negtrack.tpcNSigmaPi()); + registry.fill(HIST("PiMinusTOF_K0"), negtrack.pt(), negtrack.tofNSigmaPi()); + + return true; + } + + template + bool isSelectedLambda(V0 const& candidate, float posZ, float posY, float posX) + { + bool isL = false; // Is lambda candidate + bool isAL = false; // Is anti-lambda candidate + + double mlambda = candidate.mLambda(); + double mantilambda = candidate.mAntiLambda(); + + // separate the positive and negative V0 daughters + auto postrack = candidate.template posTrack_as(); + auto negtrack = candidate.template negTrack_as(); + + registry.fill(HIST("hLambdaCount"), 0.5); + if (postrack.pt() < resoCutVals[kPosTrackPt][iLambda] || negtrack.pt() < resoCutVals[kNegTrackPt][iLambda]) + return false; + + registry.fill(HIST("hLambdaCount"), 1.5); + if (mlambda > resoCutVals[kMassMin][iLambda] && mlambda < resoCutVals[kMassMax][iLambda]) + isL = true; + if (mantilambda > resoCutVals[kMassMin][iLambda] && mantilambda < resoCutVals[kMassMax][iLambda]) + isAL = true; + + if (!isL && !isAL) { + return false; + } + registry.fill(HIST("hLambdaCount"), 2.5); + + // Rapidity correction + if (candidate.yLambda() > resoCutVals[kRapidity][iLambda]) + return false; + registry.fill(HIST("hLambdaCount"), 3.5); + // DCA cuts for lambda and antilambda + if (isL) { + if (std::abs(candidate.dcapostopv()) < resoCutVals[kDCAPosToPVMin][iLambda] || std::abs(candidate.dcanegtopv()) < resoCutVals[kDCANegToPVMin][iLambda]) + return false; + } + if (isAL) { + if (std::abs(candidate.dcapostopv()) < resoCutVals[kDCANegToPVMin][iLambda] || std::abs(candidate.dcanegtopv()) < resoCutVals[kDCAPosToPVMin][iLambda]) + return false; + } + registry.fill(HIST("hLambdaCount"), 4.5); + if (std::abs(candidate.dcaV0daughters()) > resoSwitchVals[kDCABetDaug][iLambda]) + return false; + registry.fill(HIST("hLambdaCount"), 5.5); + // v0 radius cuts + if (resoSwitchVals[kUseV0Radius][iLambda] && (candidate.v0radius() < resoCutVals[kRadiusMin][iLambda] || candidate.v0radius() > resoCutVals[kRadiusMax][iLambda])) + return false; + registry.fill(HIST("hLambdaCount"), 6.5); + // cosine pointing angle cuts + if (candidate.v0cosPA() < resoCutVals[kCosPA][iLambda]) + return false; + registry.fill(HIST("hLambdaCount"), 7.5); + // Proper lifetime + float cTauLambda = candidate.distovertotmom(posX, posY, posZ) * massLambda; + if (resoSwitchVals[kUseProperLifetime][iLambda] && cTauLambda > resoCutVals[kLifeTime][iLambda]) + return false; + registry.fill(HIST("hLambdaCount"), 8.5); + if (isL) { + if (!selectionV0Daughter(postrack, kProtons) || !selectionV0Daughter(negtrack, kPions)) + return false; + } + if (isAL) { + if (!selectionV0Daughter(postrack, kPions) || !selectionV0Daughter(negtrack, kProtons)) + return false; + } + registry.fill(HIST("hLambdaCount"), 9.5); + + if (!cfgPIDConfigs.cfgUseAntiLambda && isAL) { // Reject the track if it is antilambda + return false; + } + + registry.fill(HIST("hLambdaPhi"), candidate.phi()); + registry.fill(HIST("hLambdaEta"), candidate.eta()); + registry.fill(HIST("PrPlusTPC_L"), postrack.pt(), postrack.tpcNSigmaPr()); + registry.fill(HIST("PrPlusTOF_L"), postrack.pt(), postrack.tofNSigmaPr()); + registry.fill(HIST("PiMinusTPC_L"), negtrack.pt(), negtrack.tpcNSigmaPi()); + registry.fill(HIST("PiMinusTOF_L"), negtrack.pt(), negtrack.tofNSigmaPi()); + + return true; + } + + template + void fillCorrelationsTPCFT0(TV0Tracks tracks1, TFT0s const& ft0, float posZ, float posY, float posX, int system, int corType, 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; + // loop over all tracks + for (auto const& track1 : tracks1) { + double resoMass = -1; + + // 4 = kshort, 5 = lambda, 6 = phi + if (cfgPIDConfigs.cfgPIDParticle == kK0) { + if (!isSelectedK0(track1, posZ, posY, posX)) + continue; // Reject if called for K0 but V0 is not K0 + + resoMass = track1.mK0Short(); + } + + if (cfgPIDConfigs.cfgPIDParticle == kLambda) { + if (!isSelectedLambda(track1, posZ, posY, posX)) + continue; // Reject if called for Lambda but V0 is not lambda + + resoMass = track1.mLambda(); + } + + if (system == SameEvent) { + if (corType == kFT0C) { + registry.fill(HIST("Trig_hist_TPC_FT0C"), fSampleIndex, posZ, track1.pt(), resoMass, eventWeight * triggerWeight); + } else if (corType == kFT0A) { + registry.fill(HIST("Trig_hist_TPC_FT0A"), fSampleIndex, posZ, track1.pt(), resoMass, eventWeight * triggerWeight); + } + } + + std::size_t channelSize = 0; + if (corType == kFT0C) { + channelSize = ft0.channelC().size(); + } else if (corType == kFT0A) { + channelSize = ft0.channelA().size(); + } else { + LOGF(fatal, "Cor Index %d out of range", corType); + } + for (std::size_t iCh = 0; iCh < channelSize; iCh++) { + int chanelid = 0; + float ampl = 0.; + getChannel(ft0, iCh, chanelid, ampl, corType, system); + + auto phi = getPhiFT0(chanelid, corType); + auto eta = getEtaFT0(chanelid, corType); + + float deltaPhi = RecoDecay::constrainAngle(track1.phi() - phi, -PIHalf); + float deltaEta = track1.eta() - eta; + // fill the right sparse and histograms + if (system == SameEvent) { + if (corType == kFT0A) { + if (cfgQaCheck) { + registry.fill(HIST("Assoc_amp_same_TPC_FT0A"), chanelid, ampl); + registry.fill(HIST("deltaEta_deltaPhi_same_TPC_FT0A"), deltaPhi, deltaEta, ampl * eventWeight * triggerWeight); + } + sameTpcFt0a->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), resoMass, deltaPhi, deltaEta, ampl * eventWeight * triggerWeight); + } else if (corType == kFT0C) { + if (cfgQaCheck) { + registry.fill(HIST("Assoc_amp_same_TPC_FT0C"), chanelid, ampl); + registry.fill(HIST("deltaEta_deltaPhi_same_TPC_FT0C"), deltaPhi, deltaEta, ampl * eventWeight * triggerWeight); + } + sameTpcFt0c->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), resoMass, deltaPhi, deltaEta, ampl * eventWeight * triggerWeight); + } + } else if (system == MixedEvent) { + if (corType == kFT0A) { + if (cfgQaCheck) { + registry.fill(HIST("Assoc_amp_mixed_TPC_FT0A"), chanelid, ampl); + registry.fill(HIST("deltaEta_deltaPhi_mixed_TPC_FT0A"), deltaPhi, deltaEta, ampl * eventWeight * triggerWeight); + } + mixedTpcFt0a->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), resoMass, deltaPhi, deltaEta, ampl * eventWeight * triggerWeight); + } else if (corType == kFT0C) { + if (cfgQaCheck) { + registry.fill(HIST("Assoc_amp_mixed_TPC_FT0C"), chanelid, ampl); + registry.fill(HIST("deltaEta_deltaPhi_mixed_TPC_FT0C"), deltaPhi, deltaEta, ampl * eventWeight * triggerWeight); + } + mixedTpcFt0c->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), resoMass, deltaPhi, deltaEta, ampl * eventWeight * triggerWeight); + } + } + } + } + } + + template + void fillCorrelationsFT0AFT0C(TFT0s const& ft0Col1, TFT0s const& ft0Col2, float posZ, int system, int multiplicity, 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; + std::size_t channelASize = ft0Col1.channelA().size(); + std::size_t channelCSize = ft0Col2.channelC().size(); + // loop over all tracks + for (std::size_t iChA = 0; iChA < channelASize; iChA++) { + + int chanelAid = 0; + float amplA = 0.; + getChannel(ft0Col1, iChA, chanelAid, amplA, kFT0A, system); + auto phiA = getPhiFT0(chanelAid, kFT0A); + auto etaA = getEtaFT0(chanelAid, kFT0A); + + if (system == SameEvent) { + registry.fill(HIST("Trig_hist_FT0A_FT0C"), fSampleIndex, posZ, 0.5, eventWeight * amplA); + } + + for (std::size_t iChC = 0; iChC < channelCSize; iChC++) { + int chanelCid = 0; + float amplC = 0.; + getChannel(ft0Col2, iChC, chanelCid, amplC, kFT0C, system); + auto phiC = getPhiFT0(chanelCid, kFT0C); + auto etaC = getEtaFT0(chanelCid, kFT0C); + float deltaPhi = RecoDecay::constrainAngle(phiA - phiC, -PIHalf); + float deltaEta = etaA - etaC; + + // fill the right sparse and histograms + if (system == SameEvent) { + if (cfgQaCheck) { + registry.fill(HIST("deltaEta_deltaPhi_same_FT0A_FT0C"), deltaPhi, deltaEta, amplA * amplC * eventWeight * triggerWeight); + } + sameFt0aFt0c->getPairHist()->Fill(step, fSampleIndex, posZ, 0.5, multiplicity, deltaPhi, deltaEta, amplA * amplC * eventWeight * triggerWeight); + } else if (system == MixedEvent) { + if (cfgQaCheck) { + registry.fill(HIST("deltaEta_deltaPhi_mixed_FT0A_FT0C"), deltaPhi, deltaEta, amplA * amplC * eventWeight * triggerWeight); + } + mixedFt0aFt0c->getPairHist()->Fill(step, fSampleIndex, posZ, 0.5, multiplicity, deltaPhi, deltaEta, amplA * amplC * eventWeight * triggerWeight); + } + } + } + } + + bool isGoodRun(int runNumber) + { + for (const auto& ExcludedRun : cfgRunRemoveList.value) { + if (runNumber == ExcludedRun) { + return false; + } + } + + return true; + } + + double massKaPlus = o2::constants::physics::MassKPlus; // same as MassKMinus + double massLambda = o2::constants::physics::MassLambda; // same as MassLambda0 and MassLambda0Bar + double massK0Short = o2::constants::physics::MassK0Short; // same as o2::constants::physics::MassK0 and o2::constants::physics::MassK0Bar + + void processSameTpcFt0a(FilteredCollisions::iterator const& collision, FilteredTracks const& tracks, aod::FT0s const&, aod::BCsWithTimestamps const&, aod::V0Datas const& V0s) + { + if (cfgQaCheck) { + eventSelectedIndividually(collision); + } + + if (!collision.sel8()) + return; + + if (!eventRct(collision, true)) + return; + + auto bc = collision.bc_as(); + + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } + + if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), true)) + return; + + if (!collision.has_foundFT0()) + return; + + loadAlignParam(bc.timestamp()); + loadGain(bc); + loadCorrection(bc.timestamp()); + float eventWeight = 1.0f; + + registry.fill(HIST("zVtx"), collision.posZ()); + + registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin + + fillYield(collision, tracks); + + double multiplicity = tracks.size(); + + if (cfgQaCheck) + registry.fill(HIST("Nch"), multiplicity); + + if (cfgStrictTrackCounter) { + trackCounter(tracks, multiplicity); + } + + if (cfgQaCheck) { + registry.fill(HIST("Nch_corrected"), multiplicity); + } + + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; + } + + const auto& ft0 = collision.foundFT0(); + fillCorrelationsTPCFT0(V0s, ft0, collision.posZ(), collision.posY(), collision.posX(), SameEvent, kFT0A, eventWeight); + } + PROCESS_SWITCH(CorrReso, processSameTpcFt0a, "Process same event for TPC-FT0 correlation", false); + + void processMixedTpcFt0a(FilteredCollisions const& collisions, FilteredTracks const& tracks, aod::FT0s const&, aod::BCsWithTimestamps const&, aod::V0Datas const& V0s) + { + + auto getTracksSize = [&tracks, this](FilteredCollisions::iterator const& collision) { + auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); + auto mult = associatedTracks.size(); + return mult; + }; + + using MixedBinning = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getTracksSize)>; + + MixedBinning binningOnVtxAndMult{{getTracksSize}, {axisVtxMix, axisMultMix}, true}; + + auto tracksTuple = std::make_tuple(V0s, tracks); + Pair pairs{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip + for (auto it = pairs.begin(); it != pairs.end(); it++) { + auto& [collision1, v0s1, collision2, tracks2] = *it; + + if (!collision1.sel8() || !collision2.sel8()) + continue; + + if (!eventRct(collision1, false) || !eventRct(collision2, false)) + continue; + + auto tracks1 = tracks.sliceByCached(o2::aod::track::collisionId, collision1.globalIndex(), this->cache); + + if (cfgUseAdditionalEventCut && !eventSelected(collision1, tracks1.size(), false)) + continue; + if (cfgUseAdditionalEventCut && !eventSelected(collision2, tracks2.size(), false)) + continue; + + if (!(collision1.has_foundFT0() && collision2.has_foundFT0())) + continue; + + registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin + + auto bc = collision1.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } + + loadAlignParam(bc.timestamp()); + loadCorrection(bc.timestamp()); + float eventWeight = 1.0f; + + double multiplicity = tracks1.size(); + + if (cfgStrictTrackCounter) { + trackCounter(tracks1, multiplicity); + } + + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; + } + + const auto& ft0 = collision2.foundFT0(); + fillCorrelationsTPCFT0(v0s1, ft0, collision1.posZ(), collision1.posY(), collision1.posX(), MixedEvent, kFT0A, eventWeight); + } + } + PROCESS_SWITCH(CorrReso, processMixedTpcFt0a, "Process mixed events for TPC-FT0A correlation", false); + + void processSameTpcFt0c(FilteredCollisions::iterator const& collision, FilteredTracks const& tracks, aod::FT0s const&, aod::BCsWithTimestamps const&, aod::V0Datas const& V0s) + { + + if (cfgQaCheck) { + eventSelectedIndividually(collision); + } + + if (!collision.sel8()) + return; + + if (!eventRct(collision, true)) + return; + + auto bc = collision.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } + + if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), true)) + return; + if (!collision.has_foundFT0()) + return; + loadAlignParam(bc.timestamp()); + loadGain(bc); + loadCorrection(bc.timestamp()); + + registry.fill(HIST("zVtx"), collision.posZ()); + + registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin + + const auto& ft0 = collision.foundFT0(); + + double multiplicity = tracks.size(); + + if (cfgQaCheck) + registry.fill(HIST("Nch"), multiplicity); + + if (cfgStrictTrackCounter) { + trackCounter(tracks, multiplicity); + } + + if (cfgQaCheck) { + registry.fill(HIST("Nch_corrected"), multiplicity); + } + + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; + } + + fillCorrelationsTPCFT0(V0s, ft0, collision.posZ(), collision.posY(), collision.posX(), SameEvent, kFT0C, 1.0f); + } + PROCESS_SWITCH(CorrReso, processSameTpcFt0c, "Process same event for TPC-FT0C correlation", false); + + void processMixedTpcFt0c(FilteredCollisions const& collisions, FilteredTracks const& tracks, aod::FT0s const&, aod::BCsWithTimestamps const&, aod::V0Datas const& V0s) + { + + auto getTracksSize = [&tracks, this](FilteredCollisions::iterator const& collision) { + auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); + auto mult = associatedTracks.size(); + return mult; + }; + + using MixedBinning = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getTracksSize)>; + + MixedBinning binningOnVtxAndMult{{getTracksSize}, {axisVtxMix, axisMultMix}, true}; + + auto tracksTuple = std::make_tuple(V0s, tracks); + Pair pairs{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip + for (auto it = pairs.begin(); it != pairs.end(); it++) { + auto& [collision1, v0s1, collision2, tracks2] = *it; + if (!collision1.sel8() || !collision2.sel8()) + continue; + + if (!eventRct(collision1, false) || !eventRct(collision2, false)) + continue; + + auto tracks1 = tracks.sliceByCached(o2::aod::track::collisionId, collision1.globalIndex(), this->cache); + + if (cfgUseAdditionalEventCut && !eventSelected(collision1, tracks1.size(), false)) + continue; + + if (cfgUseAdditionalEventCut && !eventSelected(collision2, tracks2.size(), false)) + continue; + + if (!(collision1.has_foundFT0() && collision2.has_foundFT0())) + continue; + + registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin + auto bc = collision1.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } + loadAlignParam(bc.timestamp()); + loadCorrection(bc.timestamp()); + float eventWeight = 1.0f; + + const auto& ft0 = collision2.foundFT0(); + double multiplicity = tracks1.size(); + + if (cfgStrictTrackCounter) { + trackCounter(tracks, multiplicity); + } + + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; + } + + fillCorrelationsTPCFT0(v0s1, ft0, collision1.posZ(), collision1.posY(), collision1.posX(), MixedEvent, kFT0C, eventWeight); + } + } + PROCESS_SWITCH(CorrReso, processMixedTpcFt0c, "Process mixed events for TPC-FT0C correlation", false); + + void processSameFt0aFt0c(FilteredCollisions::iterator const& collision, FilteredTracks const& tracks, aod::FT0s const&, aod::BCsWithTimestamps const&) + { + + if (cfgQaCheck) { + eventSelectedIndividually(collision); + } + + if (!collision.sel8()) + return; + + if (!eventRct(collision, true)) + return; + + auto bc = collision.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } + + if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), true)) + return; + + if (!collision.has_foundFT0()) + return; + + loadAlignParam(bc.timestamp()); + loadGain(bc); + loadCorrection(bc.timestamp()); + float eventWeight = 1.0f; + + const auto& ft0 = collision.foundFT0(); + + fillYield(collision, tracks); + + registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin + + double multiplicity = tracks.size(); + + if (cfgQaCheck) + registry.fill(HIST("Nch"), multiplicity); + + if (cfgStrictTrackCounter) { + trackCounter(tracks, multiplicity); + } + + if (cfgQaCheck) { + registry.fill(HIST("Nch_corrected"), multiplicity); + } + + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; + } + + fillCorrelationsFT0AFT0C(ft0, ft0, collision.posZ(), SameEvent, multiplicity, eventWeight); + } + PROCESS_SWITCH(CorrReso, processSameFt0aFt0c, "Process same event for FT0A-FT0C correlation", true); + + void processMixedFt0aFt0c(FilteredCollisions const& collisions, FilteredTracks const& tracks, aod::FT0s const&, aod::BCsWithTimestamps const&) + { + + auto getTracksSize = [&tracks, this](FilteredCollisions::iterator const& collision) { + auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); + auto mult = associatedTracks.size(); + return mult; + }; + + using MixedBinning = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getTracksSize)>; + + MixedBinning binningOnVtxAndMult{{getTracksSize}, {axisVtxMix, axisMultMix}, true}; + + auto tracksTuple = std::make_tuple(tracks, tracks); + Pair pairs{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, 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; + + // should have the same event to TPC-FT0A/C correlations + if (!collision1.sel8() || !collision2.sel8()) + continue; + + if (!eventRct(collision1, false) || !eventRct(collision2, false)) + continue; + + if (cfgUseAdditionalEventCut && !eventSelected(collision1, tracks1.size(), false)) + continue; + if (cfgUseAdditionalEventCut && !eventSelected(collision2, tracks2.size(), false)) + continue; + + if (!(collision1.has_foundFT0() && collision2.has_foundFT0())) + continue; + + auto bc = collision1.bc_as(); + int currentRunNumber = bc.runNumber(); + if (!cfgRunRemoveList.value.empty()) { + if (!isGoodRun(currentRunNumber)) // Rejects runs if bad run number + return; + } + loadAlignParam(bc.timestamp()); + loadCorrection(bc.timestamp()); + float eventWeight = 1.0f; + + const auto& ft0Col1 = collision1.foundFT0(); + const auto& ft0Col2 = collision2.foundFT0(); + + double multiplicity = tracks1.size(); + + if (cfgStrictTrackCounter) { + trackCounter(tracks1, multiplicity); + } + + if (multiplicity > cfgMaxMultForCorrelations || multiplicity < cfgMinMultForCorrelations) { + return; + } + + registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin + + fillCorrelationsFT0AFT0C(ft0Col1, ft0Col2, collision1.posZ(), MixedEvent, multiplicity, eventWeight); + } + } + PROCESS_SWITCH(CorrReso, processMixedFt0aFt0c, "Process mixed events for FT0A-FT0C correlation", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + }; +} diff --git a/PWGCF/TwoParticleCorrelations/Tasks/corrSparse.cxx b/PWGCF/TwoParticleCorrelations/Tasks/corrSparse.cxx deleted file mode 100644 index 08a490c2e1d..00000000000 --- a/PWGCF/TwoParticleCorrelations/Tasks/corrSparse.cxx +++ /dev/null @@ -1,2125 +0,0 @@ -// 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 CorrSparse.cxx -/// \brief Provides a sparse with usefull two particle correlation info -/// \author Thor Jensen (thor.kjaersgaard.jensen@cern.ch) - -#include "PWGCF/Core/CorrelationContainer.h" -#include "PWGMM/Mult/DataModel/bestCollisionTable.h" - -#include "Common/CCDB/EventSelectionParams.h" -#include "Common/Core/RecoDecay.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/TrackSelectionTables.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace o2; -using namespace o2::framework; -using namespace o2::framework::expressions; -using namespace constants::math; - -namespace o2::aod -{ -namespace corrsparse -{ -DECLARE_SOA_COLUMN(Multiplicity, multiplicity, int); -} -DECLARE_SOA_TABLE(Multiplicity, "AOD", "MULTIPLICITY", - corrsparse::Multiplicity); - -} // namespace o2::aod - -// define the filtered collisions and tracks -#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; - -// static constexpr float LongArrayFloat[3][20] = {{1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}, {2.1, 2.2, 2.3, -2.1, -2.2, -2.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}, {3.1, 3.2, 3.3, -3.1, -3.2, -3.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2, 1.3, -1.1, -1.2, -1.3, 1.1, 1.2}}; -// static constexpr int LongArrayInt[3][20] = {{1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1}, {2, 2, 2, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1}, {3, 3, 3, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1}}; - -struct CorrSparse { - Service ccdb; - - O2_DEFINE_CONFIGURABLE(cfgUseAdditionalEventCut, bool, false, "Use additional event cut on mult correlations") - O2_DEFINE_CONFIGURABLE(cfgZVtxCut, float, 10.0f, "Accepted z-vertex range") - O2_DEFINE_CONFIGURABLE(cfgQaCheck, bool, false, "Fill QA histograms for multiplicity and zVtx for events used in the analysis") - - struct : ConfigurableGroup{ - O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.2f, "minimum accepted track pT") - O2_DEFINE_CONFIGURABLE(cfgPtCutMax, float, 10.0f, "maximum accepted track pT") - O2_DEFINE_CONFIGURABLE(cfgEtaCut, float, 0.8f, "Eta cut") - O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5f, "max chi2 per TPC clusters") - O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 50.0f, "minimum TPC clusters") - O2_DEFINE_CONFIGURABLE(cfgCutTPCCrossedRows, float, 70.0f, "minimum TPC crossed rows") - O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 5.0f, "minimum ITS clusters") - O2_DEFINE_CONFIGURABLE(cfgCutDCAz, float, 2.0f, "max DCA to vertex z") - - } cfgTrackCuts; - - struct : ConfigurableGroup{ - - O2_DEFINE_CONFIGURABLE(processFV0, bool, true, "Process FV0 correlations") - O2_DEFINE_CONFIGURABLE(processFT0A, bool, true, "Process FT0A correlations") - O2_DEFINE_CONFIGURABLE(processFT0C, bool, true, "Process FT0C correlations") - O2_DEFINE_CONFIGURABLE(processMFT, bool, true, "Process MFT correlations") - O2_DEFINE_CONFIGURABLE(withGain, bool, false, "Use gain for FT0A and FT0C") - - } cfgDetectorConfig; - - struct : ConfigurableGroup { - O2_DEFINE_CONFIGURABLE(cfgPtCutMinMFT, float, 0.5f, "minimum accepted MFT track pT") - O2_DEFINE_CONFIGURABLE(cfgPtCutMaxMFT, float, 10.0f, "maximum accepted MFT track pT") - O2_DEFINE_CONFIGURABLE(etaMftTrackMin, float, -3.6, "Minimum eta for MFT track") - O2_DEFINE_CONFIGURABLE(etaMftTrackMax, float, -2.5, "Maximum eta for MFT track") - O2_DEFINE_CONFIGURABLE(nClustersMftTrack, int, 5, "Minimum number of clusters for MFT track") - Configurable cutBestCollisionId{"cutBestCollisionId", 0, "cut on the best collision Id used in a filter"}; - Configurable etaMftTrackMaxFilter{"etaMftTrackMaxFilter", -2.0f, "Maximum value for the eta of MFT tracks when used in filter"}; - Configurable etaMftTrackMinFilter{"etaMftTrackMinFilter", -3.9f, "Minimum value for the eta of MFT tracks when used in filter"}; - Configurable mftMaxDCAxy{"mftMaxDCAxy", 2.0f, "Cut on dcaXY for MFT tracks"}; - Configurable mftMaxDCAz{"mftMaxDCAz", 2.0f, "Cut on dcaZ for MFT tracks"}; - } cfgMftConfig; - - struct : ConfigurableGroup{ - O2_DEFINE_CONFIGURABLE(cfgMinMult, int, 0, "Minimum multiplicity for collision") - O2_DEFINE_CONFIGURABLE(cfgMaxMult, int, 10, "Maximum multiplicity for collision") - O2_DEFINE_CONFIGURABLE(cfgEvSelkNoSameBunchPileup, bool, false, "rejects collisions which are associated with the same found-by-T0 bunch crossing") - O2_DEFINE_CONFIGURABLE(cfgEvSelkNoITSROFrameBorder, bool, false, "reject events at ITS ROF border") - O2_DEFINE_CONFIGURABLE(cfgEvSelkNoTimeFrameBorder, bool, false, "reject events at TF border") - O2_DEFINE_CONFIGURABLE(cfgEvSelkIsGoodZvtxFT0vsPV, bool, false, "removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference, use this cut at low multiplicities with caution") - O2_DEFINE_CONFIGURABLE(cfgEvSelkNoCollInTimeRangeStandard, bool, false, "no collisions in specified time range") - O2_DEFINE_CONFIGURABLE(cfgEvSelkIsGoodITSLayer0123, bool, true, "cut time intervals with dead ITS layers 0,1,2,3") - O2_DEFINE_CONFIGURABLE(cfgEvSelkIsGoodITSLayersAll, bool, true, "cut time intervals with dead ITS staves") - O2_DEFINE_CONFIGURABLE(cfgEvSelkNoCollInRofStandard, bool, false, "no other collisions in this Readout Frame with per-collision multiplicity above threshold") - O2_DEFINE_CONFIGURABLE(cfgEvSelkNoHighMultCollInPrevRof, bool, false, "veto an event if FT0C amplitude in previous ITS ROF is above threshold") - O2_DEFINE_CONFIGURABLE(cfgEvSelMultCorrelation, bool, true, "Multiplicity correlation cut") - O2_DEFINE_CONFIGURABLE(cfgEvSelV0AT0ACut, bool, true, "V0A T0A 5 sigma cut") - O2_DEFINE_CONFIGURABLE(cfgEvSelOccupancy, bool, true, "Occupancy cut") - O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 2000, "High cut on TPC occupancy") - O2_DEFINE_CONFIGURABLE(cfgCutOccupancyLow, int, 0, "Low cut on TPC occupancy") - - } cfgEventSelection; - - 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") - O2_DEFINE_CONFIGURABLE(cfgRemapFT0ADeadChannels, bool, false, "If true, remap FT0A channels 60-63 to amplitudes from 92-95 respectively") - O2_DEFINE_CONFIGURABLE(cfgRemapFT0CDeadChannels, bool, false, "If true, remap FT0C channels 177->145, 176->144, 178->146, 179->147, 139->115")} cfgFITConfig; - - O2_DEFINE_CONFIGURABLE(cfgMinMixEventNum, int, 5, "Minimum number of events to mix") - O2_DEFINE_CONFIGURABLE(cfgMergingCut, float, 0.02, "Merging cut on track merge") - O2_DEFINE_CONFIGURABLE(cfgApplyTwoTrackEfficiency, bool, true, "Apply two track efficiency for tpc tpc") - O2_DEFINE_CONFIGURABLE(cfgRadiusLow, float, 0.8, "Low radius for merging cut") - O2_DEFINE_CONFIGURABLE(cfgRadiusHigh, float, 2.5, "High radius for merging cut") - O2_DEFINE_CONFIGURABLE(cfgSampleSize, double, 10, "Sample size for mixed event") - O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object") - O2_DEFINE_CONFIGURABLE(cfgCentralityWeight, std::string, "", "CCDB path to centrality weight object") - O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, bool, false, "Use local efficiency object") - O2_DEFINE_CONFIGURABLE(cfgUseEventWeights, bool, false, "Use event weights for mixed event") - - 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"); - O2_DEFINE_CONFIGURABLE(cfgMultT0CCutEnabled, bool, false, "Enable Global multiplicity vs T0C centrality cut") - Configurable> cfgMultT0CCutPars{"cfgMultT0CCutPars", std::vector{143.04, -4.58368, 0.0766055, -0.000727796, 2.86153e-06, 23.3108, -0.36304, 0.00437706, -4.717e-05, 1.98332e-07}, "Global multiplicity vs T0C centrality cut parameter values"}; - O2_DEFINE_CONFIGURABLE(cfgMultPVT0CCutEnabled, bool, false, "Enable PV multiplicity vs T0C centrality cut") - Configurable> cfgMultPVT0CCutPars{"cfgMultPVT0CCutPars", std::vector{195.357, -6.15194, 0.101313, -0.000955828, 3.74793e-06, 30.0326, -0.43322, 0.00476265, -5.11206e-05, 2.13613e-07}, "PV multiplicity vs T0C centrality cut parameter values"}; - O2_DEFINE_CONFIGURABLE(cfgMultMultPVHighCutFunction, std::string, "[0]+[1]*x + 5.*([2]+[3]*x)", "Functional for multiplicity correlation cut"); - O2_DEFINE_CONFIGURABLE(cfgMultMultPVLowCutFunction, std::string, "[0]+[1]*x - 5.*([2]+[3]*x)", "Functional for multiplicity correlation cut"); - O2_DEFINE_CONFIGURABLE(cfgMultGlobalPVCutEnabled, bool, false, "Enable global multiplicity vs PV multiplicity cut") - Configurable> cfgMultGlobalPVCutPars{"cfgMultGlobalPVCutPars", std::vector{-0.140809, 0.734344, 2.77495, 0.0165935}, "PV multiplicity vs T0C centrality cut parameter values"}; - O2_DEFINE_CONFIGURABLE(cfgMultMultV0AHighCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + 4.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); - O2_DEFINE_CONFIGURABLE(cfgMultMultV0ALowCutFunction, 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"); - O2_DEFINE_CONFIGURABLE(cfgMultMultV0ACutEnabled, bool, false, "Enable global multiplicity vs V0A multiplicity cut") - Configurable> cfgMultMultV0ACutPars{"cfgMultMultV0ACutPars", std::vector{534.893, 184.344, 0.423539, -0.00331436, 5.34622e-06, 871.239, 53.3735, -0.203528, 0.000122758, 5.41027e-07}, "Global multiplicity vs V0A multiplicity cut parameter values"}; - std::vector multT0CCutPars; - std::vector multPVT0CCutPars; - std::vector multGlobalPVCutPars; - std::vector multMultV0ACutPars; - TF1* fMultPVT0CCutLow = nullptr; - TF1* fMultPVT0CCutHigh = nullptr; - TF1* fMultT0CCutLow = nullptr; - TF1* fMultT0CCutHigh = nullptr; - TF1* fMultGlobalPVCutLow = nullptr; - TF1* fMultGlobalPVCutHigh = nullptr; - TF1* fMultMultV0ACutLow = nullptr; - TF1* fMultMultV0ACutHigh = nullptr; - TF1* fT0AV0AMean = nullptr; - TF1* fT0AV0ASigma = nullptr; - } cfgFuncParas; - - Configurable cfgCutFV0{"cfgCutFV0", 50., "FV0A threshold"}; - Configurable cfgCutFT0A{"cfgCutFT0A", 150., "FT0A threshold"}; - Configurable cfgCutFT0C{"cfgCutFT0C", 50., "FT0C threshold"}; - Configurable cfgCutZDC{"cfgCutZDC", 10., "ZDC threshold"}; - - SliceCache cache; - SliceCache cacheNch; - - ConfigurableAxis axisVertex{"axisVertex", {10, -10, 10}, "vertex axis for histograms"}; - ConfigurableAxis axisEta{"axisEta", {40, -1., 1.}, "eta axis for histograms"}; - ConfigurableAxis axisPhi{"axisPhi", {72, 0.0, constants::math::TwoPI}, "phi axis for histograms"}; - ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt axis for histograms"}; - ConfigurableAxis axisAmbiguity{"axisAmbiguity", {100, 0, 100}, "MFT track ambiguity axis for histograms"}; - - ConfigurableAxis axisEtaMft{"axisEtaMFT", {40, -3.6, -2.4}, "eta axis for MFT tracks in histograms"}; - ConfigurableAxis axisEtaFt0a{"axisEtaFT0A", {40, 3.5, 4.9}, "eta axis for FT0A in histograms"}; - ConfigurableAxis axisEtaFt0c{"axisEtaFT0C", {40, -3.3, -2.1}, "eta axis for FT0C in histograms"}; - ConfigurableAxis axisEtaFv0{"axisEtaFV0", {40, 2.2, 5.1}, "eta axis for FV0 in histograms"}; - - ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, PIHalf * 3}, "delta phi axis for histograms"}; - ConfigurableAxis axisDeltaEta{"axisDeltaEta", {48, -1.6, 1.6}, "delta eta axis for histograms"}; - ConfigurableAxis axisDeltaEtaTpcMft{"axisDeltaEtaTPCMFT", {48, 1.6, 4.6}, "delta eta axis for TPC-MFT histograms"}; - ConfigurableAxis axisDeltaEtaTpcFv0{"axisDeltaEtaTPCFV0", {48, -1.7, -5.9}, "delta eta axis for TPC-FV0 histograms"}; - ConfigurableAxis axisDeltaEtaTpcFt0a{"axisDeltaEtaTPCFT0A", {48, -5.7, -2.7}, "delta eta axis for TPC-FT0A histograms"}; - ConfigurableAxis axisDeltaEtaTpcFt0c{"axisDeltaEtaTPCFT0C", {48, 1.3, 4.1}, "delta eta axis for TPC-FT0C histograms"}; - ConfigurableAxis axisDeltaEtaMftFt0c{"axisDeltaEtaMFTFT0C", {48, -2.0, 0.6}, "delta eta axis for MFT-FT0C histograms"}; - ConfigurableAxis axisDeltaEtaMftFt0a{"axisDeltaEtaMFTFT0A", {48, -8.5, -5.9}, "delta eta axis for MFT-FT0A histograms"}; - ConfigurableAxis axisDeltaEtaMftFv0{"axisDeltaEtaMFTFV0", {48, -8.6, -4.7}, "delta eta axis for MFT-FV0 histograms"}; - ConfigurableAxis axisDeltaEtaFt0s{"axisDeltaEtaFT0S", {48, -8.6, 5.0}, "delta eta axis for FT0S histograms"}; - - ConfigurableAxis axisPtTrigger{"axisPtTrigger", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt trigger axis for histograms"}; - ConfigurableAxis axisPtAssoc{"axisPtAssoc", {VARIABLE_WIDTH, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0, 10.0}, "pt associated axis for histograms"}; - ConfigurableAxis axisMultiplicity{"axisMultiplicity", {VARIABLE_WIDTH, 0, 5, 10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100}, "multiplicity / centrality axis for histograms"}; - ConfigurableAxis vtxMix{"vtxMix", {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 multMix{"multMix", {VARIABLE_WIDTH, 0, 5, 10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100}, "multiplicity / centrality axis for mixed event histograms"}; - ConfigurableAxis axisSample{"axisSample", {cfgSampleSize, 0, cfgSampleSize}, "sample axis for histograms"}; - - ConfigurableAxis axisVertexEfficiency{"axisVertexEfficiency", {10, -10, 10}, "vertex axis for efficiency histograms"}; - ConfigurableAxis axisEtaEfficiency{"axisEtaEfficiency", {20, -1.0, 1.0}, "eta axis for efficiency histograms"}; - ConfigurableAxis axisPtEfficiency{"axisPtEfficiency", {VARIABLE_WIDTH, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0}, "pt axis for efficiency histograms"}; - - ConfigurableAxis axisAmplitudeFt0a{"axisAmplitudeFt0a", {5000, 0, 1000}, "FT0A amplitude"}; - ConfigurableAxis axisChannelFt0aAxis{"axisChannelFt0aAxis", {96, 0.0, 96.0}, "FT0A channel"}; - - Configurable cfgGainEqPath{"cfgGainEqPath", "Analysis/EventPlane/GainEq", "CCDB path for gain equalization constants"}; - Configurable cfgCorrLevel{"cfgCorrLevel", 1, "calibration step: 0 = no corr, 1 = gain corr"}; - ConfigurableAxis cfgaxisFITamp{"cfgaxisFITamp", {1000, 0, 5000}, ""}; - AxisSpec axisFit{cfgaxisFITamp, "fit amplitude"}; - AxisSpec axisChID = {220, 0, 220}; - - // make the filters and cuts. - Filter collisionFilter = (nabs(aod::collision::posZ) < cfgZVtxCut); - Filter trackFilter = (nabs(aod::track::eta) < cfgTrackCuts.cfgEtaCut) && (cfgTrackCuts.cfgPtCutMin < aod::track::pt) && (cfgTrackCuts.cfgPtCutMax > aod::track::pt) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)) && (aod::track::tpcChi2NCl < cfgTrackCuts.cfgCutChi2prTPCcls) && (aod::track::dcaZ < cfgTrackCuts.cfgCutDCAz); - - Filter mftTrackEtaFilter = ((aod::fwdtrack::eta < cfgMftConfig.etaMftTrackMaxFilter) && (aod::fwdtrack::eta > cfgMftConfig.etaMftTrackMinFilter)); - - // Filters below will be used for uncertainties - Filter mftTrackCollisionIdFilter = (aod::fwdtrack::bestCollisionId >= 0); - Filter mftTrackDcaXYFilter = (nabs(aod::fwdtrack::bestDCAXY) < cfgMftConfig.mftMaxDCAxy); - // Filter mftTrackDcaZFilter = (nabs(aod::fwdtrack::bestDCAZ) < cfgMftConfig.mftMaxDCAz); - - // using AodCollisions = soa::Filtered>; // aod::CentFT0Cs - // using AodTracks = soa::Filtered>; - - using AodCollisions = soa::Filtered>; - using AodTracks = soa::Filtered>; - using FilteredMftTracks = soa::Filtered; - using Reassociated2DMftTracks = aod::BestCollisionsFwd; - - Preslice perColGlobal = aod::track::collisionId; - - // FT0 geometry - o2::ft0::Geometry ft0Det; - o2::fv0::Geometry* fv0Det{}; - static constexpr uint64_t Ft0IndexA = 96; - std::vector* offsetFT0; - std::vector* offsetFV0; - - std::vector cstFT0RelGain{}; - - // Corrections - TH3D* mEfficiency = nullptr; - TH1D* mCentralityWeight = nullptr; - bool correctionsLoaded = false; - - // Define the outputs - - OutputObj same{"sameEvent"}; - OutputObj mixed{"mixedEvent"}; - - HistogramRegistry registry{"registry"}; - - // Define Global variables - - enum EventCutTypes { - kFilteredEvents = 0, - kAfterSel8, - kUseNoTimeFrameBorder, - kUseNoITSROFrameBorder, - kUseNoSameBunchPileup, - kUseGoodZvtxFT0vsPV, - kUseNoCollInTimeRangeStandard, - kUseGoodITSLayersAll, - kUseGoodITSLayer0123, - kUseNoCollInRofStandard, - kUseNoHighMultCollInPrevRof, - kUseOccupancy, - kUseMultCorrCut, - kUseT0AV0ACut, - kUseVertexITSTPC, - kUseTVXinTRD, - kNEventCuts - }; - - enum EventType { - SameEvent = 1, - MixedEvent = 3 - }; - - enum FITIndex { - kFT0A = 0, - kFT0C, - kFV0 - }; - - enum DetectorChannels { - kFT0AInnerRingMin = 0, - kFT0AInnerRingMax = 31, - kFT0AOuterRingMin = 32, - kFT0AOuterRingMax = 95, - kFT0CInnerRingMin = 96, - kFT0CInnerRingMax = 143, - kFT0COuterRingMin = 144, - kFT0COuterRingMax = 207 - }; - - std::array, 16> eventCuts; - - void init(InitContext&) - { - - const AxisSpec axisPhi{72, 0.0, constants::math::TwoPI, "#varphi"}; - const AxisSpec axisEta{40, -1., 1., "#eta"}; - const AxisSpec axisEtaFull{90, -4., 5., "#eta"}; - - ccdb->setURL("http://alice-ccdb.cern.ch"); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - - auto now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); - ccdb->setCreatedNotAfter(now); - fv0Det = o2::fv0::Geometry::instance(o2::fv0::Geometry::eUninitialized); - - LOGF(info, "Starting init"); - - // Event Counter - if ((doprocessSameTpcFIT || doprocessSameTpcMft || doprocessSameTPC || doprocessSameMFTFIT || doprocessSameTpcMftReassociated2D || doprocessSameMftReassociated2DFIT || doprocessSameFT0s) && cfgUseAdditionalEventCut) { - registry.add("hEventCountSpecific", "Number of Event;; Count", {HistType::kTH1D, {{13, 0, 13}}}); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(1, "after sel8"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(2, "kNoSameBunchPileup"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(3, "kNoITSROFrameBorder"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(4, "kNoTimeFrameBorder"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(5, "kIsGoodZvtxFT0vsPV"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(6, "kNoCollInTimeRangeStandard"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(7, "kIsGoodITSLayer0123"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(8, "kIsGoodITSLayersAll"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(9, "kNoCollInRofStandard"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(10, "kNoHighMultCollInPrevRof"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(11, "occupancy"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(12, "MultCorrelation"); - registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(13, "cfgEvSelV0AT0ACut"); - } - - if (doprocessSameTpcMftReassociated2D || doprocessSameMftReassociated2DFIT) { - registry.add("hEventCountMftReassoc", "Number of Events;; Count", {HistType::kTH1D, {{4, 0, 4}}}); - registry.get(HIST("hEventCountMftReassoc"))->GetXaxis()->SetBinLabel(1, "all MFT tracks"); - registry.get(HIST("hEventCountMftReassoc"))->GetXaxis()->SetBinLabel(2, "MFT tracks after selection"); - registry.get(HIST("hEventCountMftReassoc"))->GetXaxis()->SetBinLabel(3, "ambiguous MFT tracks"); - registry.get(HIST("hEventCountMftReassoc"))->GetXaxis()->SetBinLabel(4, "non-ambiguous MFT tracks"); - - registry.add("ReassociatedMftTracks", "Reassociated MFT tracks", {HistType::kTH1D, {{2, 0, 2}}}); - registry.get(HIST("ReassociatedMftTracks"))->GetXaxis()->SetBinLabel(1, "Not Reassociated MFT tracks"); - registry.get(HIST("ReassociatedMftTracks"))->GetXaxis()->SetBinLabel(2, "Reassociated MFT tracks"); - } - - // Make histograms to check the distributions after cuts - if (doprocessSameTpcFIT || doprocessSameTpcMft || doprocessSameTPC || doprocessSameMFTFIT || doprocessSameTpcMftReassociated2D || doprocessSameMftReassociated2DFIT || doprocessSameFT0s) { - registry.add("Phi", "Phi", {HistType::kTH1D, {axisPhi}}); - if (doprocessSameMFTFIT) { - registry.add("Eta", "EtaMFT", {HistType::kTH1D, {axisEtaMft}}); - } - if (doprocessSameTpcFIT || doprocessSameTpcMft || doprocessSameTPC || doprocessSameTpcMftReassociated2D) { - registry.add("Eta", "Eta", {HistType::kTH1D, {axisEta}}); - } - registry.add("EtaCorrected", "EtaCorrected", {HistType::kTH1D, {axisEta}}); - registry.add("pT", "pT", {HistType::kTH1D, {axisPtTrigger}}); - registry.add("pTCorrected", "pTCorrected", {HistType::kTH1D, {axisPtTrigger}}); - registry.add("Nch", "N_{ch}", {HistType::kTH1D, {axisMultiplicity}}); - registry.add("Nch_used", "N_{ch}", {HistType::kTH1D, {axisMultiplicity}}); // histogram to see how many events are in the same and mixed event - registry.add("zVtx", "zVtx", {HistType::kTH1D, {axisVertex}}); - registry.add("zVtx_used", "zVtx_used", {HistType::kTH1D, {axisVertex}}); - - if (doprocessSameTpcFIT || doprocessSameMFTFIT || doprocessSameMftReassociated2DFIT || doprocessSameFT0s) { - registry.add("FT0Amp", "", {HistType::kTH2F, {axisChID, axisFit}}); - registry.add("FV0Amp", "", {HistType::kTH2F, {axisChID, axisFit}}); - registry.add("FT0AmpCorrect", "", {HistType::kTH2F, {axisChID, axisFit}}); - registry.add("EtaPhi", "", {HistType::kTH2F, {axisEtaFull, axisPhi}}); - } - } - - if (doprocessSameTpcFIT) { - - if (cfgDetectorConfig.processFT0A) { - registry.add("deltaEta_deltaPhi_same_TPC_FT0A", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFt0a}}); // check to see the delta eta and delta phi distribution - registry.add("deltaEta_deltaPhi_mixed_TPC_FT0A", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFt0a}}); - registry.add("Assoc_amp_same", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); - registry.add("Assoc_amp_mixed", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); - registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); - } - if (cfgDetectorConfig.processFT0C) { - registry.add("deltaEta_deltaPhi_same_TPC_FT0C", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFt0c}}); // check to see the delta eta and delta phi distribution - registry.add("deltaEta_deltaPhi_mixed_TPC_FT0C", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFt0c}}); - registry.add("Assoc_amp_same", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); - registry.add("Assoc_amp_mixed", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); - registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); - } - if (cfgDetectorConfig.processFV0) { - registry.add("deltaEta_deltaPhi_same_TPC_FV0", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFv0}}); // check to see the delta eta and delta phi distribution - registry.add("deltaEta_deltaPhi_same_TPC_FV0", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFv0}}); // check to see the delta eta and delta phi distribution - registry.add("deltaEta_deltaPhi_mixed_TPC_FV0", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcFv0}}); - registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); - } - } - - if (doprocessSameMFTFIT || doprocessSameMftReassociated2DFIT) { - - if (cfgDetectorConfig.processFT0A) { - registry.add("deltaEta_deltaPhi_same_MFT_FT0A", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaMftFt0a}}); // check to see the delta eta and delta phi distribution - registry.add("deltaEta_deltaPhi_mixed_MFT_FT0A", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaMftFt0a}}); - registry.add("Assoc_amp_same", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); - registry.add("Assoc_amp_mixed", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); - registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); - } - if (cfgDetectorConfig.processFT0C) { - registry.add("deltaEta_deltaPhi_same_MFT_FT0C", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaMftFt0c}}); // check to see the delta eta and delta phi distribution - registry.add("deltaEta_deltaPhi_mixed_MFT_FT0C", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaMftFt0c}}); - registry.add("Assoc_amp_same", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); - registry.add("Assoc_amp_mixed", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); - registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); - } - if (cfgDetectorConfig.processFV0) { - registry.add("deltaEta_deltaPhi_same_MFT_FV0", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaMftFv0}}); // check to see the delta eta and delta phi distribution - registry.add("deltaEta_deltaPhi_mixed_MFT_FV0", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaMftFv0}}); - registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); - registry.add("Assoc_amp_same", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); - registry.add("Assoc_amp_mixed", "", {HistType::kTH2D, {axisChannelFt0aAxis, axisAmplitudeFt0a}}); - } - } - - if (doprocessSameTpcMft || doprocessSameTpcMftReassociated2D) { - registry.add("deltaEta_deltaPhi_same_TPC_MFT", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcMft}}); // check to see the delta eta and delta phi distribution - registry.add("deltaEta_deltaPhi_mixed_TPC_MFT", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaTpcMft}}); - registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); - } - - if (doprocessSameTPC) { - registry.add("deltaEta_deltaPhi_same_TPC_TPC", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEta}}); // check to see the delta eta and delta phi distribution - registry.add("deltaEta_deltaPhi_mixed_TPC_TPC", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEta}}); - registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); - } - - if (doprocessSameFT0s) { - registry.add("deltaEta_deltaPhi_same_FT0A_FT0C", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaFt0s}}); // check to see the delta eta and delta phi distribution - registry.add("deltaEta_deltaPhi_mixed_FT0A_FT0C", "", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEtaFt0s}}); - registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); - } - - registry.add("eventcount", "bin", {HistType::kTH1F, {{4, 0, 4, "bin"}}}); // histogram to see how many events are in the same and mixed event - - LOGF(info, "Initializing correlation container"); - - std::vector corrAxisTpcFt0c = {{axisSample, "Sample"}, - {axisVertex, "z-vtx (cm)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisDeltaEtaTpcFt0c, "#Delta#eta"}}; - - std::vector effAxis = { - {axisEtaEfficiency, "#eta"}, - {axisPtEfficiency, "p_{T} (GeV/c)"}, - {axisVertexEfficiency, "z-vtx (cm)"}, - }; - - std::vector userAxis; - - // Correlation axis For TPC-FIT - - std::vector corrAxisTpcFt0a = {{axisSample, "Sample"}, - {axisVertex, "z-vtx (cm)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisDeltaEtaTpcFt0a, "#Delta#eta"}}; - std::vector corrAxisTpcFv0 = {{axisSample, "Sample"}, - {axisVertex, "z-vtx (cm)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisDeltaEtaTpcFv0, "#Delta#eta"}}; - - // correlation axis for TPC-TPC - std::vector corrAxisTpcTpc = {{axisSample, "Sample"}, - {axisVertex, "z-vtx (cm)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisDeltaEta, "#Delta#eta"}}; - - // Correlation axis For TPC-MFT - std::vector corrAxisTpcMft = {{axisSample, "Sample"}, - {axisVertex, "z-vtx (cm)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisDeltaEtaTpcMft, "#Delta#eta"}}; - - // Correlation axis For MFT-FIT - std::vector corrAxisMftFt0a = {{axisSample, "Sample"}, - {axisVertex, "z-vtx (cm)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisDeltaEtaMftFt0a, "#Delta#eta"}}; - std::vector corrAxisMftFt0c = {{axisSample, "Sample"}, - {axisVertex, "z-vtx (cm)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisDeltaEtaMftFt0c, "#Delta#eta"}}; - std::vector corrAxisMftFv0 = {{axisSample, "Sample"}, - {axisVertex, "z-vtx (cm)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisDeltaEtaMftFv0, "#Delta#eta"}}; - - std::vector corrAxisFT0s = {{axisSample, "Sample"}, - {axisVertex, "z-vtx (cm)"}, - {axisPtTrigger, "p_{T} (GeV/c)"}, - {axisPtAssoc, "p_{T} (GeV/c)"}, - {axisDeltaPhi, "#Delta#varphi (rad)"}, - {axisDeltaEtaFt0s, "#Delta#eta"}}; - - if (doprocessSameTpcFIT) { - if (cfgDetectorConfig.processFT0A) { - same.setObject(new CorrelationContainer("sameEvent_TPC_FT0A", "sameEvent_TPC_FT0A", corrAxisTpcFt0a, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_TPC_FT0A", "mixedEvent_TPC_FT0A", corrAxisTpcFt0a, effAxis, userAxis)); - } - if (cfgDetectorConfig.processFT0C) { - same.setObject(new CorrelationContainer("sameEvent_TPC_FT0C", "sameEvent_TPC_FT0C", corrAxisTpcFt0c, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_TPC_FT0C", "mixedEvent_TPC_FT0C", corrAxisTpcFt0c, effAxis, userAxis)); - } - if (cfgDetectorConfig.processFV0) { - 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 (doprocessSameMFTFIT) { - if (cfgDetectorConfig.processFT0A) { - same.setObject(new CorrelationContainer("sameEvent_MFT_FT0A", "sameEvent_MFT_FT0A", corrAxisMftFt0a, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_MFT_FT0A", "mixedEvent_MFT_FT0A", corrAxisMftFt0a, effAxis, userAxis)); - } - if (cfgDetectorConfig.processFT0C) { - same.setObject(new CorrelationContainer("sameEvent_MFT_FT0C", "sameEvent_MFT_FT0C", corrAxisMftFt0c, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_MFT_FT0C", "mixedEvent_MFT_FT0C", corrAxisMftFt0c, effAxis, userAxis)); - } - if (cfgDetectorConfig.processFV0) { - same.setObject(new CorrelationContainer("sameEvent_MFT_FV0", "sameEvent_MFT_FV0", corrAxisMftFv0, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_MFT_FV0", "mixedEvent_MFT_FV0", corrAxisMftFv0, effAxis, userAxis)); - } - } - - if (doprocessSameMftReassociated2DFIT) { - if (cfgDetectorConfig.processFT0A) { - same.setObject(new CorrelationContainer("sameEvent_MFT_Reassociated2D_FT0A", "sameEvent_MFT_Reassociated2D_FT0A", corrAxisMftFt0a, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_MFT_Reassociated2D_FT0A", "mixedEvent_MFT_Reassociated2D_FT0A", corrAxisMftFt0a, effAxis, userAxis)); - } - if (cfgDetectorConfig.processFT0C) { - same.setObject(new CorrelationContainer("sameEvent_MFT_Reassociated2D_FT0C", "sameEvent_MFT_Reassociated2D_FT0C", corrAxisMftFt0c, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_MFT_Reassociated2D_FT0C", "mixedEvent_MFT_Reassociated2D_FT0C", corrAxisMftFt0c, effAxis, userAxis)); - } - if (cfgDetectorConfig.processFV0) { - same.setObject(new CorrelationContainer("sameEvent_MFT_Reassociated2D_FV0", "sameEvent_MFT_Reassociated2D_FV0", corrAxisMftFv0, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_MFT_Reassociated2D_FV0", "mixedEvent_MFT_Reassociated2D_FV0", corrAxisMftFv0, effAxis, userAxis)); - } - } - - if (doprocessSameTPC) { - same.setObject(new CorrelationContainer("sameEvent_TPC_TPC", "sameEvent_TPC_TPC", corrAxisTpcTpc, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_TPC_TPC", "mixedEvent_TPC_TPC", corrAxisTpcTpc, effAxis, userAxis)); - } - - if (doprocessSameTpcMft) { - same.setObject(new CorrelationContainer("sameEvent_TPC_MFT", "sameEvent_TPC_MFT", corrAxisTpcMft, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_TPC_MFT", "mixedEvent_TPC_MFT", corrAxisTpcMft, effAxis, userAxis)); - } - - if (doprocessSameTpcMftReassociated2D) { - same.setObject(new CorrelationContainer("sameEvent_TPC_MFT_Reassociated2D", "sameEvent_TPC_MFT_Reassociated2D", corrAxisTpcMft, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_TPC_MFT_Reassociated2D", "mixedEvent_TPC_MFT_Reassociated2D", corrAxisTpcMft, effAxis, userAxis)); - } - - if (doprocessSameFT0s) { - same.setObject(new CorrelationContainer("sameEvent_FT0A_FT0C", "sameEvent_FT0A_FT0C", corrAxisTpcFt0c, effAxis, userAxis)); - mixed.setObject(new CorrelationContainer("mixedEvent_FT0A_FT0C", "mixedEvent_FT0A_FT0C", corrAxisTpcFt0c, effAxis, userAxis)); - } - LOGF(info, "End of init"); - } - - TRandom3* gRandom = new TRandom3(); - - template - bool eventSelected(TCollision collision, const int multTrk, const bool fillCounter) - { - registry.fill(HIST("hEventCountSpecific"), 0.5); - if (cfgEventSelection.cfgEvSelkNoSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { - // rejects collisions which are associated with the same "found-by-T0" bunch crossing - // https://indico.cern.ch/event/1396220/#1-event-selection-with-its-rof - return 0; - } - if (fillCounter && cfgEventSelection.cfgEvSelkNoSameBunchPileup) - registry.fill(HIST("hEventCountSpecific"), 1.5); - if (cfgEventSelection.cfgEvSelkNoITSROFrameBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { - return 0; - } - if (fillCounter && cfgEventSelection.cfgEvSelkNoITSROFrameBorder) - registry.fill(HIST("hEventCountSpecific"), 2.5); - if (cfgEventSelection.cfgEvSelkNoTimeFrameBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) { - return 0; - } - if (fillCounter && cfgEventSelection.cfgEvSelkNoTimeFrameBorder) - registry.fill(HIST("hEventCountSpecific"), 3.5); - if (cfgEventSelection.cfgEvSelkIsGoodZvtxFT0vsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { - // removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference - // use this cut at low multiplicities with caution - return 0; - } - if (fillCounter && cfgEventSelection.cfgEvSelkIsGoodZvtxFT0vsPV) - registry.fill(HIST("hEventCountSpecific"), 4.5); - if (cfgEventSelection.cfgEvSelkNoCollInTimeRangeStandard && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { - // no collisions in specified time range - return 0; - } - if (fillCounter && cfgEventSelection.cfgEvSelkNoCollInTimeRangeStandard) - registry.fill(HIST("hEventCountSpecific"), 5.5); - - if (cfgEventSelection.cfgEvSelkIsGoodITSLayer0123 && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayer0123)) { - // from Jan 9 2025 AOT meeting - // cut time intervals with dead ITS staves - return 0; - } - if (fillCounter && cfgEventSelection.cfgEvSelkIsGoodITSLayer0123) - registry.fill(HIST("hEventCountSpecific"), 6.5); - - if (cfgEventSelection.cfgEvSelkIsGoodITSLayersAll && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { - // from Jan 9 2025 AOT meeting - // cut time intervals with dead ITS staves - return 0; - } - - if (fillCounter && cfgEventSelection.cfgEvSelkIsGoodITSLayersAll) - registry.fill(HIST("hEventCountSpecific"), 7.5); - - if (cfgEventSelection.cfgEvSelkNoCollInRofStandard && !collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard)) { - // no other collisions in this Readout Frame with per-collision multiplicity above threshold - return 0; - } - if (fillCounter && cfgEventSelection.cfgEvSelkNoCollInRofStandard) - registry.fill(HIST("hEventCountSpecific"), 8.5); - if (cfgEventSelection.cfgEvSelkNoHighMultCollInPrevRof && !collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof)) { - // veto an event if FT0C amplitude in previous ITS ROF is above threshold - return 0; - } - if (fillCounter && cfgEventSelection.cfgEvSelkNoHighMultCollInPrevRof) - registry.fill(HIST("hEventCountSpecific"), 9.5); - auto occupancy = collision.trackOccupancyInTimeRange(); - if (cfgEventSelection.cfgEvSelOccupancy && (occupancy < cfgEventSelection.cfgCutOccupancyLow || occupancy > cfgEventSelection.cfgCutOccupancyHigh)) - return 0; - if (fillCounter && cfgEventSelection.cfgEvSelOccupancy) - registry.fill(HIST("hEventCountSpecific"), 10.5); - - auto multNTracksPV = collision.multNTracksPV(); - - if (cfgFuncParas.cfgMultGlobalPVCutEnabled) { - if (multTrk < cfgFuncParas.fMultGlobalPVCutLow->Eval(multNTracksPV)) - return 0; - if (multTrk > cfgFuncParas.fMultGlobalPVCutHigh->Eval(multNTracksPV)) - return 0; - } - if (cfgFuncParas.cfgMultMultV0ACutEnabled) { - if (collision.multFV0A() < cfgFuncParas.fMultMultV0ACutLow->Eval(multTrk)) - return 0; - if (collision.multFV0A() > cfgFuncParas.fMultMultV0ACutHigh->Eval(multTrk)) - return 0; - } - - if (fillCounter && cfgEventSelection.cfgEvSelMultCorrelation) - registry.fill(HIST("hEventCountSpecific"), 11.5); - - // V0A T0A 5 sigma cut - float sigma = 5.0; - if (cfgEventSelection.cfgEvSelV0AT0ACut && (std::fabs(collision.multFV0A() - cfgFuncParas.fT0AV0AMean->Eval(collision.multFT0A())) > sigma * cfgFuncParas.fT0AV0ASigma->Eval(collision.multFT0A()))) - return 0; - if (fillCounter && cfgEventSelection.cfgEvSelV0AT0ACut) - registry.fill(HIST("hEventCountSpecific"), 12.5); - - return 1; - } - - double getPhiFV0(uint64_t chno) - { - o2::fv0::Point3Dsimple chPos{}; - int const cellsInLeft[] = {0, 1, 2, 3, 8, 9, 10, 11, 16, 17, 18, 19, 24, 25, 26, 27, 32, 40, 33, 41, 34, 42, 35, 43}; - bool const isChnoInLeft = std::find(std::begin(cellsInLeft), std::end(cellsInLeft), chno) != std::end(cellsInLeft); - - if (isChnoInLeft) { - chPos = fv0Det->getReadoutCenter(chno); - return RecoDecay::phi(chPos.x + (*offsetFV0)[0].getX(), chPos.y + (*offsetFV0)[0].getY()); - } else { - chPos = fv0Det->getReadoutCenter(chno); - return RecoDecay::phi(chPos.x + (*offsetFV0)[1].getX(), chPos.y + (*offsetFV0)[1].getY()); - } - } - - double getPhiFT0(uint64_t chno, int i) - { - // offsetFT0[0]: FT0A, offsetFT0[1]: FT0C - if (i > 1 || i < 0) { - LOGF(fatal, "kFIT Index %d out of range", i); - } - - ft0Det.calculateChannelCenter(); - auto chPos = ft0Det.getChannelCenter(chno); - return RecoDecay::phi(chPos.X() + (*offsetFT0)[i].getX(), chPos.Y() + (*offsetFT0)[i].getY()); - } - - double getEtaFV0(uint64_t chno) - { - - int const cellsInLeft[] = {0, 1, 2, 3, 8, 9, 10, 11, 16, 17, 18, 19, 24, 25, 26, 27, 32, 40, 33, 41, 34, 42, 35, 43}; - bool const isChnoInLeft = std::find(std::begin(cellsInLeft), std::end(cellsInLeft), chno) != std::end(cellsInLeft); - - o2::fv0::Point3Dsimple chPos{}; - chPos = fv0Det->getReadoutCenter(chno); - - float offsetX, offsetY, offsetZ; - - if (isChnoInLeft) { - offsetX = (*offsetFV0)[0].getX(); - offsetY = (*offsetFV0)[0].getY(); - offsetZ = (*offsetFV0)[0].getZ(); - } else { - offsetX = (*offsetFV0)[1].getX(); - offsetY = (*offsetFV0)[1].getY(); - offsetZ = (*offsetFV0)[1].getZ(); - } - - auto x = chPos.x + offsetX; - auto y = chPos.y + offsetY; - auto z = chPos.z + offsetZ; - - auto r = std::sqrt(x * x + y * y); - auto theta = std::atan2(r, z); - - return -std::log(std::tan(0.5 * theta)); - } - - double getEtaFT0(uint64_t chno, int i) - { - // offsetFT0[0]: FT0A, offsetFT0[1]: FT0C - if (i > 1 || i < 0) { - LOGF(fatal, "kFIT Index %d out of range", i); - } - ft0Det.calculateChannelCenter(); - auto chPos = ft0Det.getChannelCenter(chno); - auto x = chPos.X() + (*offsetFT0)[i].getX(); - auto y = chPos.Y() + (*offsetFT0)[i].getY(); - auto z = chPos.Z() + (*offsetFT0)[i].getZ(); - if (chno >= Ft0IndexA) { - z = -z; - } - auto r = std::sqrt(x * x + y * y); - auto theta = std::atan2(r, z); - return -std::log(std::tan(0.5 * theta)); - } - - template - bool isAcceptedMftTrack(TTrack const& mftTrack) - { - // cut on the eta of MFT tracks - if (mftTrack.eta() < cfgMftConfig.etaMftTrackMin || mftTrack.eta() > cfgMftConfig.etaMftTrackMax) { - return false; - } - - // cut on the number of clusters of the reconstructed MFT track - if (mftTrack.nClusters() < cfgMftConfig.nClustersMftTrack) { - return false; - } - - if (mftTrack.pt() < cfgMftConfig.cfgPtCutMinMFT || mftTrack.pt() > cfgMftConfig.cfgPtCutMaxMFT) { - return false; - } - - return true; - } - - template - bool isAmbiguousMftTrack(TTrack const& mftTrack, bool fillHistogram) - { - if (mftTrack.ambDegree() > 1) { - if (fillHistogram) { - registry.fill(HIST("hEventCountMftReassoc"), 2.5); // fill histogram for events with at least one ambiguous track); - } - return false; - } - registry.fill(HIST("hEventCountMftReassoc"), 3.5); // fill histogram for events without ambiguous tracks - return true; - } - - void loadAlignParam(uint64_t timestamp) - { - offsetFT0 = ccdb->getForTimeStamp>("FT0/Calib/Align", timestamp); - offsetFV0 = ccdb->getForTimeStamp>("FV0/Calib/Align", timestamp); - if (offsetFT0 == nullptr) { - LOGF(fatal, "Could not load FT0/Calib/Align for timestamp %d", timestamp); - } - if (offsetFV0 == nullptr) { - LOGF(fatal, "Could not load FV0/Calib/Align for timestamp %d", timestamp); - } - } - - void loadGain(aod::BCsWithTimestamps::iterator const& bc) - { - cstFT0RelGain.clear(); - cstFT0RelGain = {}; - std::string fullPath; - - auto timestamp = bc.timestamp(); - constexpr int ChannelsFT0 = 208; - if (cfgCorrLevel == 0) { - for (auto i{0u}; i < ChannelsFT0; i++) { - cstFT0RelGain.push_back(1.); - } - } else { - fullPath = cfgGainEqPath; - fullPath += "/FT0"; - const auto objft0Gain = ccdb->getForTimeStamp>(fullPath, timestamp); - if (!objft0Gain) { - for (auto i{0u}; i < ChannelsFT0; i++) { - cstFT0RelGain.push_back(1.); - } - } else { - cstFT0RelGain = *(objft0Gain); - } - } - } - - template - void getChannelWithGain(TFT0s const& ft0, std::size_t const& iCh, int& id, float& ampl, int fitType) - { - if (fitType == kFT0C) { - id = ft0.channelC()[iCh]; - id = id + Ft0IndexA; - ampl = ft0.amplitudeC()[iCh]; - registry.fill(HIST("FT0Amp"), id, ampl); - ampl = ampl / cstFT0RelGain[id]; - registry.fill(HIST("FT0AmpCorrect"), id, ampl); - } else if (fitType == kFT0A) { - id = ft0.channelA()[iCh]; - ampl = ft0.amplitudeA()[iCh]; - registry.fill(HIST("FT0Amp"), id, ampl); - ampl = ampl / cstFT0RelGain[id]; - registry.fill(HIST("FT0AmpCorrect"), id, ampl); - } else { - LOGF(fatal, "Cor Index %d out of range", fitType); - } - } - - void loadCorrection(uint64_t timestamp) - { - if (correctionsLoaded) { - return; - } - if (cfgEfficiency.value.empty() == false) { - if (cfgLocalEfficiency > 0) { - TFile* fEfficiencyTrigger = TFile::Open(cfgEfficiency.value.c_str(), "READ"); - mEfficiency = reinterpret_cast(fEfficiencyTrigger->Get("ccdb_object")); - } else { - mEfficiency = ccdb->getForTimeStamp(cfgEfficiency, timestamp); - } - if (mEfficiency == nullptr) { - LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgEfficiency.value.c_str()); - } - LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)mEfficiency); - } - if (cfgCentralityWeight.value.empty() == false) { - mCentralityWeight = ccdb->getForTimeStamp(cfgCentralityWeight, timestamp); - if (mCentralityWeight == nullptr) { - LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgCentralityWeight.value.c_str()); - } - LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgCentralityWeight.value.c_str(), (void*)mCentralityWeight); - } - correctionsLoaded = true; - } - - bool getEfficiencyCorrection(float& weight_nue, float eta, float pt, float posZ) - { - float eff = 1.; - if (mEfficiency) { - int etaBin = mEfficiency->GetXaxis()->FindBin(eta); - int ptBin = mEfficiency->GetYaxis()->FindBin(pt); - int zBin = mEfficiency->GetZaxis()->FindBin(posZ); - eff = mEfficiency->GetBinContent(etaBin, ptBin, zBin); - } else { - eff = 1.0; - } - if (eff == 0) - return false; - weight_nue = 1. / eff; - return true; - } - - template - void getChannelFT0(TFT0s const& ft0, std::size_t const& iCh, int& id, float& ampl, int fitType) - { - if (fitType == kFT0C) { - id = ft0.channelC()[iCh]; - id = id + Ft0IndexA; - ampl = ft0.amplitudeC()[iCh]; - } else if (fitType == kFT0A) { - id = ft0.channelA()[iCh]; - ampl = ft0.amplitudeA()[iCh]; - } else { - LOGF(fatal, "Cor Index %d out of range", fitType); - } - registry.fill(HIST("FT0Amp"), id, ampl); - } - - template - void getChannelFV0(TFT0s const& fv0, std::size_t const& iCh, int& id, float& ampl) - { - id = fv0.channel()[iCh]; - ampl = fv0.amplitude()[iCh]; - registry.fill(HIST("FV0Amp"), id, ampl); - } - - template - bool trackSelected(TTrack track) - { - return ((track.tpcNClsFound() >= cfgTrackCuts.cfgCutTPCclu) && (track.tpcNClsCrossedRows() >= cfgTrackCuts.cfgCutTPCCrossedRows) && (track.itsNCls() >= cfgTrackCuts.cfgCutITSclu)); - } - - int getMagneticField(uint64_t timestamp) - { - // Get the magnetic field - static o2::parameters::GRPMagField* grpo = nullptr; - if (grpo == nullptr) { - grpo = ccdb->getForTimeStamp("/GLO/Config/GRPMagField", timestamp); - if (grpo == nullptr) { - LOGF(fatal, "GRP object not found for timestamp %llu", timestamp); - return 0; - } - LOGF(info, "Retrieved GRP for timestamp %llu with magnetic field of %d kG", timestamp, grpo->getNominalL3Field()); - } - return grpo->getNominalL3Field(); - } - - // fill multiple histograms - template - void fillYield(TCollision collision, TTracks tracks) // function to fill the yield and etaphi histograms. - { - registry.fill(HIST("zVtx"), collision.posZ()); - registry.fill(HIST("Nch"), tracks.size()); - - float weff1 = 1.0; - float zvtx = collision.posZ(); - - for (auto const& track1 : tracks) { - - if constexpr (std::is_same_v) { - if (!isAcceptedMftTrack(track1)) { - continue; - } - } else { - if (!trackSelected(track1)) { - continue; - } - if (!getEfficiencyCorrection(weff1, track1.eta(), track1.pt(), zvtx)) { - continue; - } - } - - registry.fill(HIST("Phi"), RecoDecay::constrainAngle(track1.phi(), 0.0)); - registry.fill(HIST("Eta"), track1.eta()); - registry.fill(HIST("EtaCorrected"), track1.eta(), weff1); - registry.fill(HIST("pT"), track1.pt()); - registry.fill(HIST("pTCorrected"), track1.pt(), weff1); - } - } - - template - float getDPhiStar(TTrack const& track1, TTrackAssoc const& track2, float radius, int magField) - { - float charge1 = track1.sign(); - float charge2 = track2.sign(); - - float phi1 = track1.phi(); - float phi2 = track2.phi(); - - float pt1 = track1.pt(); - float pt2 = track2.pt(); - - int fbSign = (magField > 0) ? 1 : -1; - - float dPhiStar = phi1 - phi2 - charge1 * fbSign * std::asin(0.075 * radius / pt1) + charge2 * fbSign * std::asin(0.075 * radius / pt2); - - if (dPhiStar > constants::math::PI) - dPhiStar = constants::math::TwoPI - dPhiStar; - if (dPhiStar < -constants::math::PI) - dPhiStar = -constants::math::TwoPI - dPhiStar; - - return dPhiStar; - } - - // Correlations for detectors and TPC - - ////////////////////////////// - //////////MFT/TPC-FIT//////// - //////////////////////////// - template - void fillCorrelationsFIT(TTracks tracks1, TTracksAssociated tracks2, FITs const&, float posZ, int system, int corType, float multiplicity) - { - - int fSampleIndex = gRandom->Uniform(0, cfgSampleSize); - float triggerWeight = 1.0f; - - // loop over all tracks - if (cfgQaCheck) { - - if (system == SameEvent) { - registry.fill(HIST("Nch_used"), multiplicity); - registry.fill(HIST("zVtx_used"), posZ); - } - } - - for (auto const& track1 : tracks1) { - - if constexpr (std::is_same_v) { - - if (!isAcceptedMftTrack(track1)) { - continue; - } - } else { - if (!trackSelected(track1)) { - continue; - } - - if (!getEfficiencyCorrection(triggerWeight, track1.eta(), track1.pt(), posZ)) - continue; - } - - if (system == SameEvent) { - registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, track1.pt(), triggerWeight); - } - - // if using FV0A for correlations / using FV0A as associated particles - if constexpr (std::is_same_v) { - - std::size_t channelSize = tracks2.channel().size(); - for (std::size_t iCh = 0; iCh < channelSize; iCh++) { - int channelID = 0; - float amplitude = 0.; - - getChannelFV0(tracks2, iCh, channelID, amplitude); - - auto phi = getPhiFV0(channelID); - auto eta = getEtaFV0(channelID); - - float deltaPhi = RecoDecay::constrainAngle(track1.phi() - phi, -PIHalf); - float deltaEta = track1.eta() - eta; - - if (system == SameEvent) { - if (cfgDetectorConfig.processMFT) { - registry.fill(HIST("deltaEta_deltaPhi_same_MFT_FV0"), deltaPhi, deltaEta, amplitude * triggerWeight); - } else { - registry.fill(HIST("deltaEta_deltaPhi_same_TPC_FV0"), deltaPhi, deltaEta, amplitude * triggerWeight); - } - - registry.fill(HIST("Assoc_amp_same"), channelID, amplitude); - same->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track1.pt(), deltaPhi, deltaEta, amplitude * triggerWeight); - - } else if (system == MixedEvent) { - if (cfgDetectorConfig.processMFT) { - registry.fill(HIST("deltaEta_deltaPhi_mixed_MFT_FV0"), deltaPhi, deltaEta, amplitude); - } else { - registry.fill(HIST("deltaEta_deltaPhi_mixed_TPC_FV0"), deltaPhi, deltaEta, amplitude); - } - registry.fill(HIST("Assoc_amp_mixed"), channelID, amplitude); - mixed->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track1.pt(), deltaPhi, deltaEta, amplitude); - } - } - } - - // if using FT0A and FT0C for correlations / using FT0A and FT0C as associated particles - if constexpr (std::is_same_v) { - - std::size_t channelSize = 0; - if (corType == kFT0C) { - channelSize = tracks2.channelC().size(); - } else if (corType == kFT0A) { - channelSize = tracks2.channelA().size(); - } else { - LOGF(fatal, "Cor Index %d out of range", corType); - } - - for (std::size_t iCh = 0; iCh < channelSize; iCh++) { - int channelID = 0; - float amplitude = 0.; - if (cfgDetectorConfig.withGain) { - getChannelWithGain(tracks2, iCh, channelID, amplitude, corType); - } else { - getChannelFT0(tracks2, iCh, channelID, amplitude, corType); - } - - // reject depending on FT0C/FT0A rings - if (corType == kFT0C) { - if ((cfgFITConfig.cfgRejectFT0CInside && (channelID >= kFT0CInnerRingMin && channelID <= kFT0CInnerRingMax)) || (cfgFITConfig.cfgRejectFT0COutside && (channelID >= kFT0COuterRingMin && channelID <= kFT0COuterRingMax))) - continue; - } - if (corType == kFT0A) { - if ((cfgFITConfig.cfgRejectFT0AInside && (channelID >= kFT0AInnerRingMin && channelID <= kFT0AInnerRingMax)) || (cfgFITConfig.cfgRejectFT0AOutside && (channelID >= kFT0AOuterRingMin && channelID <= kFT0AOuterRingMax))) - continue; - } - - auto phi = getPhiFT0(channelID, corType); - auto eta = getEtaFT0(channelID, corType); - - float deltaPhi = RecoDecay::constrainAngle(track1.phi() - phi, -PIHalf); - float deltaEta = track1.eta() - eta; - - if (system == SameEvent) { - if (corType == kFT0A) { - if (cfgDetectorConfig.processMFT) { - registry.fill(HIST("deltaEta_deltaPhi_same_MFT_FT0A"), deltaPhi, deltaEta, amplitude * triggerWeight); - } else { - registry.fill(HIST("deltaEta_deltaPhi_same_TPC_FT0A"), deltaPhi, deltaEta, amplitude * triggerWeight); - } - registry.fill(HIST("Assoc_amp_same"), channelID, amplitude); - same->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track1.pt(), deltaPhi, deltaEta, amplitude * triggerWeight); - } - if (corType == kFT0C) { - if (cfgDetectorConfig.processMFT) { - registry.fill(HIST("deltaEta_deltaPhi_same_MFT_FT0C"), deltaPhi, deltaEta, amplitude * triggerWeight); - } else { - registry.fill(HIST("deltaEta_deltaPhi_same_TPC_FT0C"), deltaPhi, deltaEta, amplitude * triggerWeight); - } - registry.fill(HIST("Assoc_amp_same"), channelID, amplitude); - same->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track1.pt(), deltaPhi, deltaEta, amplitude * triggerWeight); - } - } else if (system == MixedEvent) { - if (corType == kFT0A) { - if (cfgDetectorConfig.processMFT) { - registry.fill(HIST("deltaEta_deltaPhi_mixed_MFT_FT0A"), deltaPhi, deltaEta, amplitude); - } else { - registry.fill(HIST("deltaEta_deltaPhi_mixed_TPC_FT0A"), deltaPhi, deltaEta, amplitude); - } - registry.fill(HIST("Assoc_amp_mixed"), channelID, amplitude); - mixed->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track1.pt(), deltaPhi, deltaEta, amplitude); - } - if (corType == kFT0C) { - if (cfgDetectorConfig.processMFT) { - registry.fill(HIST("deltaEta_deltaPhi_mixed_MFT_FT0C"), deltaPhi, deltaEta, amplitude); - } else { - registry.fill(HIST("deltaEta_deltaPhi_mixed_TPC_FT0C"), deltaPhi, deltaEta, amplitude); - } - registry.fill(HIST("Assoc_amp_mixed"), channelID, amplitude); - mixed->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track1.pt(), deltaPhi, deltaEta, amplitude); - } - } - } - } - } - } - ////////////////////////// - //////////MFT-Reassociated///////// - ////////////////////////// - - template - void fillCorrelationsMftReassociatedFIT(TTracks tracks1, TTracksAssociated tracks2, FITs const&, float posZ, int system, int corType, float multiplicity, bool cutAmbiguousTracks) - { - int fSampleIndex = gRandom->Uniform(0, cfgSampleSize); - float triggerWeight = 1.0f; - - // loop over all tracks - if (cfgQaCheck) { - - if (system == SameEvent) { - registry.fill(HIST("Nch_used"), multiplicity); - registry.fill(HIST("zVtx_used"), posZ); - } - } - - for (auto const& track1 : tracks1) { - - auto reassociatedMftTrack = track1.template mfttrack_as(); - - if (!isAcceptedMftTrack(reassociatedMftTrack)) { - continue; - } - - if (isAmbiguousMftTrack(track1, false)) { - if (cutAmbiguousTracks) { - continue; - } - } - - if (system == SameEvent) { - - registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, reassociatedMftTrack.pt(), triggerWeight); - } - - if constexpr (std::is_same_v) { - - std::size_t channelSize = tracks2.channel().size(); - for (std::size_t iCh = 0; iCh < channelSize; iCh++) { - int channelID = 0; - float amplitude = 0.; - - getChannelFV0(tracks2, iCh, channelID, amplitude); - - auto phi = getPhiFV0(channelID); - auto eta = getEtaFV0(channelID); - - float deltaPhi = RecoDecay::constrainAngle(reassociatedMftTrack.phi() - phi, -PIHalf); - float deltaEta = reassociatedMftTrack.eta() - eta; - - if (system == SameEvent) { - registry.fill(HIST("Assoc_amp_same"), channelID, amplitude); - same->getPairHist()->Fill(step, fSampleIndex, posZ, reassociatedMftTrack.pt(), reassociatedMftTrack.pt(), deltaPhi, deltaEta, amplitude * triggerWeight); - registry.fill(HIST("deltaEta_deltaPhi_same_MFT_FV0"), deltaPhi, deltaEta, amplitude * triggerWeight); - } else if (system == MixedEvent) { - registry.fill(HIST("Assoc_amp_mixed"), channelID, amplitude); - mixed->getPairHist()->Fill(step, fSampleIndex, posZ, reassociatedMftTrack.pt(), reassociatedMftTrack.pt(), deltaPhi, deltaEta, amplitude); - registry.fill(HIST("deltaEta_deltaPhi_mixed_MFT_FV0"), deltaPhi, deltaEta, amplitude); - } - } - } - - // if using FT0A and FT0C for correlations / using FT0A and FT0C as associated particles - if constexpr (std::is_same_v) { - - std::size_t channelSize = 0; - if (corType == kFT0C) { - channelSize = tracks2.channelC().size(); - } else if (corType == kFT0A) { - channelSize = tracks2.channelA().size(); - } else { - LOGF(fatal, "Cor Index %d out of range", corType); - } - - for (std::size_t iCh = 0; iCh < channelSize; iCh++) { - int channelID = 0; - float amplitude = 0.; - getChannelWithGain(tracks2, iCh, channelID, amplitude, corType); - - // reject depending on FT0C/FT0A rings - if (corType == kFT0C) { - if ((cfgFITConfig.cfgRejectFT0CInside && (channelID >= kFT0CInnerRingMin && channelID <= kFT0CInnerRingMax)) || (cfgFITConfig.cfgRejectFT0COutside && (channelID >= kFT0COuterRingMin && channelID <= kFT0COuterRingMax))) - continue; - } - if (corType == kFT0A) { - if ((cfgFITConfig.cfgRejectFT0AInside && (channelID >= kFT0AInnerRingMin && channelID <= kFT0AInnerRingMax)) || (cfgFITConfig.cfgRejectFT0AOutside && (channelID >= kFT0AOuterRingMin && channelID <= kFT0AOuterRingMax))) - continue; - } - - auto phi = getPhiFT0(channelID, corType); - auto eta = getEtaFT0(channelID, corType); - - float deltaPhi = RecoDecay::constrainAngle(reassociatedMftTrack.phi() - phi, -PIHalf); - float deltaEta = reassociatedMftTrack.eta() - eta; - - if (system == SameEvent) { - if (corType == kFT0A) { - registry.fill(HIST("Assoc_amp_same"), channelID, amplitude); - same->getPairHist()->Fill(step, fSampleIndex, posZ, reassociatedMftTrack.pt(), reassociatedMftTrack.pt(), deltaPhi, deltaEta, amplitude * triggerWeight); - registry.fill(HIST("deltaEta_deltaPhi_same_MFT_FT0A"), deltaPhi, deltaEta, amplitude * triggerWeight); - } - if (corType == kFT0C) { - registry.fill(HIST("Assoc_amp_same"), channelID, amplitude); - same->getPairHist()->Fill(step, fSampleIndex, posZ, reassociatedMftTrack.pt(), reassociatedMftTrack.pt(), deltaPhi, deltaEta, amplitude * triggerWeight); - registry.fill(HIST("deltaEta_deltaPhi_same_MFT_FT0C"), deltaPhi, deltaEta, amplitude * triggerWeight); - } - } else if (system == MixedEvent) { - if (corType == kFT0A) { - registry.fill(HIST("Assoc_amp_mixed"), channelID, amplitude); - mixed->getPairHist()->Fill(step, fSampleIndex, posZ, reassociatedMftTrack.pt(), reassociatedMftTrack.pt(), deltaPhi, deltaEta, amplitude); - registry.fill(HIST("deltaEta_deltaPhi_mixed_MFT_FT0A"), deltaPhi, deltaEta, amplitude); - } - if (corType == kFT0C) { - registry.fill(HIST("Assoc_amp_mixed"), channelID, amplitude); - mixed->getPairHist()->Fill(step, fSampleIndex, posZ, reassociatedMftTrack.pt(), reassociatedMftTrack.pt(), deltaPhi, deltaEta, amplitude); - registry.fill(HIST("deltaEta_deltaPhi_mixed_MFT_FT0C"), deltaPhi, deltaEta, amplitude); - } - } - } - } - } - } - - template - void fillCorrelationsMftReassociatedTracks(TTracks tracks1, TTracksAssoc tracks2, float multiplicity, float posZ, int system, int magneticField, bool cutAmbiguousTracks) // 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; - - auto loopCounter = 0; - if (cfgQaCheck) { - - if (system == SameEvent) { - registry.fill(HIST("Nch_used"), multiplicity); - registry.fill(HIST("zVtx_used"), posZ); - } - } - // loop over all tracks - for (auto const& track1 : tracks1) { - - loopCounter++; - - if (!trackSelected(track1)) - continue; - - if (!getEfficiencyCorrection(triggerWeight, track1.eta(), track1.pt(), posZ)) - continue; - - if (system == SameEvent) { - registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, track1.pt(), triggerWeight); - } - - for (auto const& track2 : tracks2) { - - if (!cutAmbiguousTracks && system == SameEvent && (loopCounter == 1)) { - registry.fill(HIST("hEventCountMftReassoc"), 0.5); // fill histogram for events with at least one reassociated track); - } - - auto reassociatedMftTrack = track2.template mfttrack_as(); - - if (!isAcceptedMftTrack(reassociatedMftTrack)) { - continue; - } - - if (!cutAmbiguousTracks && system == SameEvent && (loopCounter == 1)) { - registry.fill(HIST("hEventCountMftReassoc"), 1.5); // fill histogram for events with at least one reassociated track after track selection); - } - - if (isAmbiguousMftTrack(track2, (!cutAmbiguousTracks && system == SameEvent && (loopCounter == 1)))) { - - if (SameEvent && (loopCounter == 1)) { - registry.fill(HIST("ReassociatedMftTracks"), 0.5); - } - if (cutAmbiguousTracks) { - continue; - } - } - - if (reassociatedMftTrack.collisionId() != track2.bestCollisionId()) { - if (SameEvent && (loopCounter == 1)) { - registry.fill(HIST("ReassociatedMftTracks"), 1.5); - } - } - - float deltaPhi = RecoDecay::constrainAngle(track1.phi() - reassociatedMftTrack.phi(), -PIHalf); - float deltaEta = track1.eta() - reassociatedMftTrack.eta(); - - if (cfgApplyTwoTrackEfficiency && std::abs(deltaEta) < cfgMergingCut) { - - double dPhiStarHigh = getDPhiStar(track1, reassociatedMftTrack, cfgRadiusHigh, magneticField); - double dPhiStarLow = getDPhiStar(track1, reassociatedMftTrack, cfgRadiusLow, magneticField); - - const double kLimit = 3.0 * cfgMergingCut; - - bool bIsBelow = false; - - if (std::abs(dPhiStarLow) < kLimit || std::abs(dPhiStarHigh) < kLimit || dPhiStarLow * dPhiStarHigh < 0) { - for (double rad(cfgRadiusLow); rad < cfgRadiusHigh; rad += 0.01) { - double dPhiStar = getDPhiStar(track1, reassociatedMftTrack, rad, magneticField); - if (std::abs(dPhiStar) < kLimit) { - bIsBelow = true; - break; - } - } - if (bIsBelow) - continue; - } - } - - // fill the right sparse and histograms - if (system == SameEvent) { - - same->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), reassociatedMftTrack.pt(), deltaPhi, deltaEta); - registry.fill(HIST("deltaEta_deltaPhi_same_TPC_MFT"), deltaPhi, deltaEta); - } else if (system == MixedEvent) { - - mixed->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), reassociatedMftTrack.pt(), deltaPhi, deltaEta); - registry.fill(HIST("deltaEta_deltaPhi_mixed_TPC_MFT"), deltaPhi, deltaEta); - } - } - } - } - - /////////////////////////////////////// - //////////TPC-TPC and TPC-MFT///////// - ///////////////////////////////////// - template - void fillCorrelations(TTracks tracks1, TTracksAssoc tracks2, float posZ, float multiplicity, int system, int magneticField) // 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; - - if (cfgQaCheck) { - if (system == SameEvent) { - registry.fill(HIST("Nch_used"), multiplicity); - registry.fill(HIST("zVtx_used"), posZ); - } - } - // loop over all tracks - for (auto const& track1 : tracks1) { - - if constexpr (std::is_same_v) { - if (!isAcceptedMftTrack(track1)) { - continue; - } - } else { - if (!trackSelected(track1)) - continue; - - if (!getEfficiencyCorrection(triggerWeight, track1.eta(), track1.pt(), posZ)) - continue; - } - - if (system == SameEvent) { - registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, track1.pt(), triggerWeight); - } - - for (auto const& track2 : tracks2) { - - if constexpr (std::is_same_v) { - if (!isAcceptedMftTrack(track2)) { - continue; - } - } else { - if (!trackSelected(track2)) - continue; - - if (track1.pt() <= track2.pt()) - continue; // skip if the trigger pt is less than the associate pt - } - - float deltaPhi = RecoDecay::constrainAngle(track1.phi() - track2.phi(), -PIHalf); - float deltaEta = track1.eta() - track2.eta(); - - if (cfgApplyTwoTrackEfficiency && std::abs(deltaEta) < cfgMergingCut) { - - double dPhiStarHigh = getDPhiStar(track1, track2, cfgRadiusHigh, magneticField); - double dPhiStarLow = getDPhiStar(track1, track2, cfgRadiusLow, magneticField); - - const double kLimit = 3.0 * cfgMergingCut; - - bool bIsBelow = false; - - if (std::abs(dPhiStarLow) < kLimit || std::abs(dPhiStarHigh) < kLimit || dPhiStarLow * dPhiStarHigh < 0) { - for (double rad(cfgRadiusLow); rad < cfgRadiusHigh; rad += 0.01) { - double dPhiStar = getDPhiStar(track1, track2, rad, magneticField); - if (std::abs(dPhiStar) < kLimit) { - bIsBelow = true; - break; - } - } - if (bIsBelow) - continue; - } - } - - // fill the right sparse and histograms - if (system == SameEvent) { - if (cfgDetectorConfig.processMFT) { - same->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta); - registry.fill(HIST("deltaEta_deltaPhi_same_TPC_MFT"), deltaPhi, deltaEta); - } else { - same->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta); - registry.fill(HIST("deltaEta_deltaPhi_same_TPC_TPC"), deltaPhi, deltaEta); - } - } else if (system == MixedEvent) { - if (cfgDetectorConfig.processMFT) { - mixed->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta); - registry.fill(HIST("deltaEta_deltaPhi_mixed_TPC_MFT"), deltaPhi, deltaEta); - } else { - mixed->getPairHist()->Fill(step, fSampleIndex, posZ, track1.pt(), track2.pt(), deltaPhi, deltaEta); - registry.fill(HIST("deltaEta_deltaPhi_mixed_TPC_TPC"), deltaPhi, deltaEta); - } - } - } - } - } - - template - void fillCorrelationsFT0s(TFT0s const& ft01, TFT0s const& ft02, float posZ, int multiplicity, 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); - - if (cfgQaCheck) { - if (system == SameEvent) { - registry.fill(HIST("zVtx_used"), posZ); - registry.fill(HIST("Nch"), multiplicity); - registry.fill(HIST("Nch_used_ft0c"), ft02.channelC().size()); - registry.fill(HIST("Nch_used_ft0a"), ft01.channelA().size()); - } - } - - float triggerWeight = 1.0f; - std::size_t channelASize = ft01.channelA().size(); - std::size_t channelCSize = ft02.channelC().size(); - // loop over all tracks - for (std::size_t iChA = 0; iChA < channelASize; iChA++) { - - int chanelAid = 0; - float amplA = 0.; - if (cfgDetectorConfig.withGain) { - getChannelWithGain(ft01, iChA, chanelAid, amplA, kFT0A); - } else { - getChannelFT0(ft01, iChA, chanelAid, amplA, kFT0A); - } - auto phiA = getPhiFT0(chanelAid, kFT0A); - auto etaA = getEtaFT0(chanelAid, kFT0A); - - if (system == SameEvent) { - registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, 0.5, eventWeight * amplA); - } - - for (std::size_t iChC = 0; iChC < channelCSize; iChC++) { - int chanelCid = 0; - float amplC = 0.; - if (cfgDetectorConfig.withGain) { - getChannelWithGain(ft02, iChC, chanelCid, amplC, kFT0C); - } else { - getChannelFT0(ft02, iChC, chanelCid, amplC, kFT0C); - } - auto phiC = getPhiFT0(chanelCid, kFT0C); - auto etaC = getEtaFT0(chanelCid, kFT0C); - float deltaPhi = RecoDecay::constrainAngle(phiA - phiC, -PIHalf); - float deltaEta = etaA - etaC; - - // fill the right sparse and histograms - if (system == SameEvent) { - registry.fill(HIST("deltaEta_deltaPhi_same_FT0A_FT0C"), deltaPhi, deltaEta, amplA * amplC * eventWeight * triggerWeight); - same->getPairHist()->Fill(step, fSampleIndex, posZ, 0.5, 0.5, deltaPhi, deltaEta, amplA * amplC * eventWeight * triggerWeight); - } else if (system == MixedEvent) { - registry.fill(HIST("deltaEta_deltaPhi_mixed_FT0A_FT0C"), deltaPhi, deltaEta, amplA * amplC * eventWeight * triggerWeight); - mixed->getPairHist()->Fill(step, fSampleIndex, posZ, 0.5, 0.5, deltaPhi, deltaEta, amplA * amplC * eventWeight * triggerWeight); - } - } - } - } - ////////////////////////////////////// - ////// Same event processing ///////// - ////////////////////////////////////// - - //////////////////////////////////// - /////////// Fwrd-Bwrd ///////////// - //////////////////////////////////// - - void processSameMFTFIT(AodCollisions::iterator const& collision, AodTracks const& tpctracks, aod::MFTTracks const& mfts, aod::FT0s const& ft0as, aod::FV0As const& fv0as, aod::BCsWithTimestamps const&) - { - - if (!collision.sel8()) - return; - - auto bc = collision.bc_as(); - - if (cfgUseAdditionalEventCut && !eventSelected(collision, tpctracks.size(), true)) - return; - - if (!collision.has_foundFT0()) - return; - loadAlignParam(bc.timestamp()); - if (cfgDetectorConfig.withGain) { - loadGain(bc); - } - - loadCorrection(bc.timestamp()); - - if ((tpctracks.size() < cfgEventSelection.cfgMinMult || tpctracks.size() >= cfgEventSelection.cfgMaxMult)) { - return; - } - if (mfts.size() == 0) { - return; - } - - registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin - fillYield(collision, mfts); - const auto& multiplicity = tpctracks.size(); - - if (cfgDetectorConfig.processFV0) { - if (collision.has_foundFV0()) { - same->fillEvent(mfts.size(), CorrelationContainer::kCFStepReconstructed); - const auto& fv0 = collision.foundFV0(); - fillCorrelationsFIT(mfts, fv0, fv0as, collision.posZ(), SameEvent, kFV0, multiplicity); - } - } - if (cfgDetectorConfig.processFT0C) { - if (collision.has_foundFT0()) { - same->fillEvent(mfts.size(), CorrelationContainer::kCFStepReconstructed); - const auto& ft0 = collision.foundFT0(); - fillCorrelationsFIT(mfts, ft0, ft0as, collision.posZ(), SameEvent, kFT0C, multiplicity); - } - } - if (cfgDetectorConfig.processFT0A) { - if (collision.has_foundFT0()) { - same->fillEvent(mfts.size(), CorrelationContainer::kCFStepReconstructed); - const auto& ft0 = collision.foundFT0(); - fillCorrelationsFIT(mfts, ft0, ft0as, collision.posZ(), SameEvent, kFT0A, multiplicity); - } - } - } - PROCESS_SWITCH(CorrSparse, processSameMFTFIT, "Process same event for MFT-FIT correlation", true); - - void processSameMftReassociated2DFIT(AodCollisions::iterator const& collision, AodTracks const& tpctracks, - soa::SmallGroups const& reassociatedMftTracks, - FilteredMftTracks const&, - aod::FT0s const& ft0as, aod::FV0As const& fv0as, aod::BCsWithTimestamps const&) - { - if (!collision.sel8()) - return; - - auto bc = collision.bc_as(); - - if (cfgUseAdditionalEventCut && !eventSelected(collision, tpctracks.size(), true)) - return; - - if (!collision.has_foundFT0()) - return; - loadAlignParam(bc.timestamp()); - loadGain(bc); - loadCorrection(bc.timestamp()); - - if ((tpctracks.size() < cfgEventSelection.cfgMinMult || tpctracks.size() >= cfgEventSelection.cfgMaxMult)) { - return; - } - if (reassociatedMftTracks.size() == 0) { - return; - } - - registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin - const auto& multiplicity = tpctracks.size(); - - if (cfgDetectorConfig.processFV0) { - if (collision.has_foundFV0()) { - same->fillEvent(reassociatedMftTracks.size(), CorrelationContainer::kCFStepReconstructed); - const auto& fv0 = collision.foundFV0(); - fillCorrelationsMftReassociatedFIT(reassociatedMftTracks, fv0, fv0as, collision.posZ(), SameEvent, kFV0, multiplicity, false); - } - } - if (cfgDetectorConfig.processFT0C) { - if (collision.has_foundFT0()) { - same->fillEvent(reassociatedMftTracks.size(), CorrelationContainer::kCFStepReconstructed); - const auto& ft0 = collision.foundFT0(); - fillCorrelationsMftReassociatedFIT(reassociatedMftTracks, ft0, ft0as, collision.posZ(), SameEvent, kFT0C, multiplicity, false); - } - } - if (cfgDetectorConfig.processFT0A) { - if (collision.has_foundFT0()) { - same->fillEvent(reassociatedMftTracks.size(), CorrelationContainer::kCFStepReconstructed); - const auto& ft0 = collision.foundFT0(); - fillCorrelationsMftReassociatedFIT(reassociatedMftTracks, ft0, ft0as, collision.posZ(), SameEvent, kFT0A, multiplicity, false); - } - } - } - PROCESS_SWITCH(CorrSparse, processSameMftReassociated2DFIT, "Process same event for MFT-FIT correlation with reassociated tracks", false); - - ///////////////////////// - ////////Mid-Mid////////// - //////////////////////// - - void processSameTPC(AodCollisions::iterator const& collision, AodTracks const& tracks, aod::BCsWithTimestamps const&) - { - - auto bc = collision.bc_as(); - - if (!collision.sel8()) - return; - - if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), true)) - return; - - loadCorrection(bc.timestamp()); - - fillYield(collision, tracks); - - if (tracks.size() < cfgEventSelection.cfgMinMult || tracks.size() >= cfgEventSelection.cfgMaxMult) { - return; - } - - registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin - - fillCorrelations(tracks, tracks, collision.posZ(), tracks.size(), SameEvent, getMagneticField(bc.timestamp())); - } - PROCESS_SWITCH(CorrSparse, processSameTPC, "Process same event for TPC-TPC correlation", false); - - ///////////////////// - ////////back-Mid-Fwrd////// - ///////////////////// - - void processSameTpcFIT(AodCollisions::iterator const& collision, AodTracks const& tracks, aod::FT0s const& ft0as, aod::FV0As const& fv0as, aod::BCsWithTimestamps const&) - { - if (!collision.sel8()) - return; - auto bc = collision.bc_as(); - if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), true)) - return; - if (!collision.has_foundFT0() && !collision.has_foundFV0()) - return; - - registry.fill(HIST("Nch"), tracks.size()); - registry.fill(HIST("zVtx"), collision.posZ()); - - if ((tracks.size() < cfgEventSelection.cfgMinMult || tracks.size() >= cfgEventSelection.cfgMaxMult)) { - return; - } - - loadAlignParam(bc.timestamp()); - if (cfgDetectorConfig.withGain) { - loadGain(bc); - } - - loadCorrection(bc.timestamp()); - - registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin - fillYield(collision, tracks); - - const auto& multiplicity = tracks.size(); - - if (cfgDetectorConfig.processFV0) { - same->fillEvent(tracks.size(), CorrelationContainer::kCFStepReconstructed); - const auto& fv0 = collision.foundFV0(); - fillCorrelationsFIT(tracks, fv0, fv0as, collision.posZ(), SameEvent, kFV0, multiplicity); - } - if (cfgDetectorConfig.processFT0C) { - same->fillEvent(tracks.size(), CorrelationContainer::kCFStepReconstructed); - const auto& ft0 = collision.foundFT0(); - fillCorrelationsFIT(tracks, ft0, ft0as, collision.posZ(), SameEvent, kFT0C, multiplicity); - } - if (cfgDetectorConfig.processFT0A) { - same->fillEvent(tracks.size(), CorrelationContainer::kCFStepReconstructed); - const auto& ft0 = collision.foundFT0(); - fillCorrelationsFIT(tracks, ft0, ft0as, collision.posZ(), SameEvent, kFT0A, multiplicity); - } - } - PROCESS_SWITCH(CorrSparse, processSameTpcFIT, "process for forward or backwards correlations with TPC", false); - - void processSameTpcMft(AodCollisions::iterator const& collision, AodTracks const& tracks, aod::MFTTracks const& mfts, aod::BCsWithTimestamps const&) - { - if (!collision.sel8()) - return; - - auto bc = collision.bc_as(); - - if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), true)) - return; - - loadCorrection(bc.timestamp()); - - registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin - - fillYield(collision, tracks); - - if (tracks.size() < cfgEventSelection.cfgMinMult || tracks.size() >= cfgEventSelection.cfgMaxMult) { - return; - } - - fillCorrelations(tracks, mfts, collision.posZ(), tracks.size(), SameEvent, getMagneticField(bc.timestamp())); - } - PROCESS_SWITCH(CorrSparse, processSameTpcMft, "Process same event for TPC-MFT correlation", false); - - void processSameTpcMftReassociated2D(AodCollisions::iterator const& collision, AodTracks const& tracks, - soa::SmallGroups const& reassociatedMftTracks, - FilteredMftTracks const&, - aod::BCsWithTimestamps const&) - { - if (!collision.sel8()) - return; - - auto bc = collision.bc_as(); - - if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), true)) - return; - - loadCorrection(bc.timestamp()); - - registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin - - fillYield(collision, tracks); - - if (tracks.size() < cfgEventSelection.cfgMinMult || tracks.size() >= cfgEventSelection.cfgMaxMult) { - return; - } - - fillCorrelationsMftReassociatedTracks(tracks, reassociatedMftTracks, tracks.size(), collision.posZ(), SameEvent, getMagneticField(bc.timestamp()), false); - } - PROCESS_SWITCH(CorrSparse, processSameTpcMftReassociated2D, "Process same event for TPC-MFT correlation with reassociated tracks", false); - - void processSameFT0s(AodCollisions::iterator const& collision, AodTracks const& tracks, aod::FT0s const&, aod::BCsWithTimestamps const&) - { - if (!collision.sel8()) - return; - - auto bc = collision.bc_as(); - - if (cfgUseAdditionalEventCut && !eventSelected(collision, 0, true)) - return; - - if (!collision.has_foundFT0()) - return; - - loadAlignParam(bc.timestamp()); - if (cfgDetectorConfig.withGain) { - loadGain(bc); - } - loadCorrection(bc.timestamp()); - - if ((tracks.size() < cfgEventSelection.cfgMinMult || tracks.size() >= cfgEventSelection.cfgMaxMult)) { - return; - } - - registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin - - const auto& ft0s = collision.foundFT0(); - auto multiplicity = tracks.size(); - - fillCorrelationsFT0s(ft0s, ft0s, collision.posZ(), multiplicity, SameEvent, 1.0f); - } - PROCESS_SWITCH(CorrSparse, processSameFT0s, "Process same event for FT0 correlation", false); - //////////////////////////////////// - ////// Mixed event processing ////// - //////////////////////////////////// - - ///////////////////////////////////////// - ////////////// Fwrd- Bwrd////////////// - //////////////////////////////////////// - - void processMixedMFTFIT(AodCollisions const& collisions, AodTracks const& tpctracks, aod::MFTTracks const& mfts, aod::FT0s const& ft0as, aod::FV0As const& fv0as, aod::BCsWithTimestamps const&) - { - - auto getTracksSize = [&tpctracks, this](AodCollisions::iterator const& collision) { - auto associatedTracks = tpctracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); - auto mult = associatedTracks.size(); - return mult; - }; - - using MixedBinning = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getTracksSize)>; - - MixedBinning binningOnVtxAndMult{{getTracksSize}, {vtxMix, multMix}, true}; - - auto tracksTuple = std::make_tuple(mfts, mfts); - Pair pair{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip - - for (auto it = pair.begin(); it != pair.end(); it++) { - auto& [collision1, tracks1, collision2, tracks2] = *it; - - if (!collision1.sel8() || !collision2.sel8()) - continue; - - if (collision1.globalIndex() == collision2.globalIndex()) { - continue; - } - - auto slicedtracks = tpctracks.sliceBy(perColGlobal, collision1.globalIndex()); - auto multiplicity = slicedtracks.size(); - - if ((multiplicity < cfgEventSelection.cfgMinMult || multiplicity >= cfgEventSelection.cfgMaxMult)) - continue; - - if (cfgUseAdditionalEventCut && !eventSelected(collision1, tracks1.size(), false)) - continue; - - if (cfgUseAdditionalEventCut && !eventSelected(collision2, tracks2.size(), false)) - continue; - - auto bc = collision1.bc_as(); - loadAlignParam(bc.timestamp()); - - loadCorrection(bc.timestamp()); - - if (cfgDetectorConfig.processFT0A) { - if (!collision1.has_foundFT0() && !collision2.has_foundFT0()) - continue; - - const auto& ft0 = collision2.foundFT0(); - fillCorrelationsFIT(tracks1, ft0, ft0as, collision1.posZ(), MixedEvent, kFT0A, multiplicity); - } - if (cfgDetectorConfig.processFT0C) { - if (!collision1.has_foundFT0() && !collision2.has_foundFT0()) - continue; - - const auto& ft0 = collision2.foundFT0(); - fillCorrelationsFIT(tracks1, ft0, ft0as, collision1.posZ(), MixedEvent, kFT0C, multiplicity); - } - if (cfgDetectorConfig.processFV0) { - if (collision1.has_foundFV0() && collision2.has_foundFV0()) { - const auto& fv0 = collision2.foundFV0(); - fillCorrelationsFIT(tracks1, fv0, fv0as, collision1.posZ(), MixedEvent, kFV0, multiplicity); - } - } - } - } - PROCESS_SWITCH(CorrSparse, processMixedMFTFIT, "Process mixed events for MFT-FIT correlation", true); - - ///////////////////////////// - //////////Mid-Mid/////////// - //////////////////////////// - - void processMixedTpcTpc(AodCollisions const& collisions, AodTracks const& tracks, aod::BCsWithTimestamps const&) - { - - auto getTracksSize = [&tracks, this](AodCollisions::iterator const& collision) { - auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); - auto mult = associatedTracks.size(); - return mult; - }; - - using MixedBinning = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getTracksSize)>; - - MixedBinning binningOnVtxAndMult{{getTracksSize}, {vtxMix, multMix}, true}; - - auto tracksTuple = std::make_tuple(tracks, tracks); - Pair pair{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip - for (auto const& [collision1, tracks1, collision2, tracks2] : pair) { - - if (!collision1.sel8() || !collision2.sel8()) - continue; - - if (collision1.globalIndex() == collision2.globalIndex()) { - continue; - } - - registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin - auto bc = collision1.bc_as(); - - loadCorrection(bc.timestamp()); - - if ((tracks1.size() < cfgEventSelection.cfgMinMult || tracks1.size() >= cfgEventSelection.cfgMaxMult)) - continue; - - if (cfgUseAdditionalEventCut && !eventSelected(collision1, tracks1.size(), false)) - continue; - if (cfgUseAdditionalEventCut && !eventSelected(collision2, tracks2.size(), false)) - continue; - - fillCorrelations(tracks1, tracks2, collision1.posZ(), tracks1.size(), MixedEvent, getMagneticField(bc.timestamp())); - } - } - PROCESS_SWITCH(CorrSparse, processMixedTpcTpc, "Process mixed events for TPC-TPC correlation", false); - - ////////////////////////////// - /////// back-Mid-Fwrd /////// - ///////////////////////////// - - void processMixedTpcFIT(AodCollisions const& collisions, AodTracks const& tracks, aod::FV0As const& fv0as, aod::FT0s const& ft0as, aod::BCsWithTimestamps const&) - { - auto getTracksSize = [&tracks, this](AodCollisions::iterator const& collision) { - auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); - auto mult = associatedTracks.size(); - return mult; - }; - - using MixedBinning = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getTracksSize)>; - - MixedBinning binningOnVtxAndMult{{getTracksSize}, {vtxMix, multMix}, true}; - - auto tracksTuple = std::make_tuple(tracks, tracks); - Pair pair{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip - - for (auto it = pair.begin(); it != pair.end(); it++) { - auto& [collision1, tracks1, collision2, tracks2] = *it; - - if (!collision1.sel8() || !collision2.sel8()) - continue; - - if (collision1.globalIndex() == collision2.globalIndex()) { - continue; - } - - if ((tracks1.size() < cfgEventSelection.cfgMinMult || tracks1.size() >= cfgEventSelection.cfgMaxMult)) - continue; - - if (cfgUseAdditionalEventCut && !eventSelected(collision1, tracks1.size(), false)) - continue; - if (cfgUseAdditionalEventCut && !eventSelected(collision2, tracks2.size(), false)) - continue; - - if (!(collision1.has_foundFT0() && collision2.has_foundFT0())) - continue; - - registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin - auto bc = collision1.bc_as(); - loadAlignParam(bc.timestamp()); - loadCorrection(bc.timestamp()); - - auto multiplicity = tracks1.size(); - - if (cfgDetectorConfig.processFT0A) { - const auto& ft0 = collision2.foundFT0(); - fillCorrelationsFIT(tracks1, ft0, ft0as, collision1.posZ(), MixedEvent, kFT0A, multiplicity); - } - if (cfgDetectorConfig.processFT0C) { - const auto& ft0 = collision2.foundFT0(); - fillCorrelationsFIT(tracks1, ft0, ft0as, collision1.posZ(), MixedEvent, kFT0C, multiplicity); - } - if (cfgDetectorConfig.processFV0) { - const auto& fv0 = collision2.foundFV0(); - fillCorrelationsFIT(tracks1, fv0, fv0as, collision1.posZ(), MixedEvent, kFV0, multiplicity); - } - } - } - PROCESS_SWITCH(CorrSparse, processMixedTpcFIT, "Process mixed events for TPC-FIT correlation", false); - - void processMixedTpcMFT(AodCollisions const& collisions, AodTracks const& tracks, aod::MFTTracks const& MFTtracks, aod::BCsWithTimestamps const&) - { - - auto getTracksSize = [&tracks, this](AodCollisions::iterator const& collision) { - auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); - auto mult = associatedTracks.size(); - return mult; - }; - - using MixedBinning = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getTracksSize)>; - - MixedBinning binningOnVtxAndMult{{getTracksSize}, {vtxMix, multMix}, true}; - - auto tracksTuple = std::make_tuple(tracks, MFTtracks); - Pair pair{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip - for (auto const& [collision1, tracks1, collision2, tracks2] : pair) { - - if (!collision1.sel8() || !collision2.sel8()) - continue; - - if (collision1.globalIndex() == collision2.globalIndex()) { - continue; - } - - registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin - auto bc = collision1.bc_as(); - - loadCorrection(bc.timestamp()); - - if ((tracks1.size() < cfgEventSelection.cfgMinMult || tracks1.size() >= cfgEventSelection.cfgMaxMult)) - continue; - - fillCorrelations(tracks1, tracks2, collision1.posZ(), tracks1.size(), MixedEvent, getMagneticField(bc.timestamp())); - } - } - PROCESS_SWITCH(CorrSparse, processMixedTpcMFT, "Process mixed events for TPC-MFT correlation", false); - - void processMixedFT0s(AodCollisions const& collisions, AodTracks const& tracks, aod::FT0s const&, aod::BCsWithTimestamps const&) - { - - auto getTracksSize = [&tracks, this](AodCollisions::iterator const& collision) { - auto associatedTracks = tracks.sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), this->cache); - auto mult = associatedTracks.size(); - return mult; - }; - - using MixedBinning = FlexibleBinningPolicy, aod::collision::PosZ, decltype(getTracksSize)>; - - MixedBinning binningOnVtxAndMult{{getTracksSize}, {vtxMix, multMix}, true}; - - auto tracksTuple = std::make_tuple(tracks, tracks); - Pair pair{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip - - for (auto it = pair.begin(); it != pair.end(); it++) { - auto& [collision1, tracks1, collision2, tracks2] = *it; - - if (!collision1.sel8() || !collision2.sel8()) - continue; - - if (collision1.globalIndex() == collision2.globalIndex()) { - continue; - } - - if ((tracks1.size() < cfgEventSelection.cfgMinMult || tracks1.size() >= cfgEventSelection.cfgMaxMult)) - continue; - - if (cfgUseAdditionalEventCut && !eventSelected(collision1, tracks1.size(), false)) - continue; - if (cfgUseAdditionalEventCut && !eventSelected(collision2, tracks2.size(), false)) - continue; - - if (!(collision1.has_foundFT0() && collision2.has_foundFT0())) - continue; - - registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin - auto bc = collision1.bc_as(); - loadAlignParam(bc.timestamp()); - loadCorrection(bc.timestamp()); - - auto multiplicity = tracks1.size(); - - const auto& ft0_1 = collision1.foundFT0(); - const auto& ft0_2 = collision2.foundFT0(); - fillCorrelationsFT0s(ft0_1, ft0_2, collision1.posZ(), multiplicity, MixedEvent, 1.0f); - } - } - PROCESS_SWITCH(CorrSparse, processMixedFT0s, "Process mixed events for FT0 correlation", false); -}; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc), - }; -}