diff --git a/Common/DataModel/Qvectors.h b/Common/DataModel/Qvectors.h index d723d659bc5..f62e8851a2b 100644 --- a/Common/DataModel/Qvectors.h +++ b/Common/DataModel/Qvectors.h @@ -155,6 +155,95 @@ using QvectorBNegVec = QvectorBNegVecs::iterator; using QvectorBTotVec = QvectorBTotVecs::iterator; ///////////////////////////////////////////////////////////////// +namespace ese_qvec +{ +DECLARE_SOA_COLUMN(IsCalibrated, isCalibrated, bool); + +DECLARE_SOA_COLUMN(EseQvecRe, eseQvecRe, std::vector); +DECLARE_SOA_COLUMN(EseQvecIm, eseQvecIm, std::vector); +DECLARE_SOA_COLUMN(EseQvecAmp, eseQvecAmp, std::vector); + +DECLARE_SOA_COLUMN(EseQvecFT0CReVec, eseQvecFT0CReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFT0CImVec, eseQvecFT0CImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFT0AReVec, eseQvecFT0AReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFT0AImVec, eseQvecFT0AImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFT0MReVec, eseQvecFT0MReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFT0MImVec, eseQvecFT0MImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFV0AReVec, eseQvecFV0AReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecFV0AImVec, eseQvecFV0AImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCposReVec, eseQvecTPCposReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCposImVec, eseQvecTPCposImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCnegReVec, eseQvecTPCnegReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCnegImVec, eseQvecTPCnegImVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCallReVec, eseQvecTPCallReVec, std::vector); +DECLARE_SOA_COLUMN(EseQvecTPCallImVec, eseQvecTPCallImVec, std::vector); + +DECLARE_SOA_COLUMN(EseQvecFT0CRe, eseQvecFT0CRe, float); +DECLARE_SOA_COLUMN(EseQvecFT0CIm, eseQvecFT0CIm, float); +DECLARE_SOA_COLUMN(EseQvecFT0ARe, eseQvecFT0ARe, float); +DECLARE_SOA_COLUMN(EseQvecFT0AIm, eseQvecFT0AIm, float); +DECLARE_SOA_COLUMN(EseQvecFT0MRe, eseQvecFT0MRe, float); +DECLARE_SOA_COLUMN(EseQvecFT0MIm, eseQvecFT0MIm, float); +DECLARE_SOA_COLUMN(EseQvecFV0ARe, eseQvecFV0ARe, float); +DECLARE_SOA_COLUMN(EseQvecFV0AIm, eseQvecFV0AIm, float); +DECLARE_SOA_COLUMN(EseQvecTPCposRe, eseQvecTPCposRe, float); +DECLARE_SOA_COLUMN(EseQvecTPCposIm, eseQvecTPCposIm, float); +DECLARE_SOA_COLUMN(EseQvecTPCnegRe, eseQvecTPCnegRe, float); +DECLARE_SOA_COLUMN(EseQvecTPCnegIm, eseQvecTPCnegIm, float); +DECLARE_SOA_COLUMN(EseQvecTPCallRe, eseQvecTPCallRe, float); +DECLARE_SOA_COLUMN(EseQvecTPCallIm, eseQvecTPCallIm, float); + +DECLARE_SOA_COLUMN(EseRedQFT0C, eseRedQFT0C, float); +DECLARE_SOA_COLUMN(EseRedQFT0A, eseRedQFT0A, float); +DECLARE_SOA_COLUMN(EseRedQFT0M, eseRedQFT0M, float); +DECLARE_SOA_COLUMN(EseRedQFV0A, eseRedQFV0A, float); +///////////////////////////////////////////////////////////////// +} // namespace ese_qvec + +DECLARE_SOA_TABLE(EseQvectors, "AOD", "ESEQVECTORDEVS", //! Table with all Qvectors. + qvec::Cent, ese_qvec::IsCalibrated, ese_qvec::EseQvecRe, ese_qvec::EseQvecIm, ese_qvec::EseQvecAmp); +using Qvector = Qvectors::iterator; + +DECLARE_SOA_TABLE(EseQvecFT0Cs, "AOD", "ESEQVECFT0C", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0CRe, ese_qvec::EseQvecFT0CIm, qvec::SumAmplFT0C); +DECLARE_SOA_TABLE(EseQvecFT0As, "AOD", "ESEQVECFT0A", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0ARe, ese_qvec::EseQvecFT0AIm, qvec::SumAmplFT0A); +DECLARE_SOA_TABLE(EseQvecFT0Ms, "AOD", "ESEQVECFT0M", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0MRe, ese_qvec::EseQvecFT0MIm, qvec::SumAmplFT0M); +DECLARE_SOA_TABLE(EseQvecFV0As, "AOD", "ESEQVECFV0A", ese_qvec::IsCalibrated, ese_qvec::EseQvecFV0ARe, ese_qvec::EseQvecFV0AIm, qvec::SumAmplFV0A); +DECLARE_SOA_TABLE(EseQvecTPCposs, "AOD", "ESEQVECTPCPOS", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCposRe, ese_qvec::EseQvecTPCposIm, qvec::NTrkTPCpos, qvec::LabelsTPCpos); +DECLARE_SOA_TABLE(EseQvecTPCnegs, "AOD", "ESEQVECTPCNEG", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCnegRe, ese_qvec::EseQvecTPCnegIm, qvec::NTrkTPCneg, qvec::LabelsTPCneg); +DECLARE_SOA_TABLE(EseQvecTPCalls, "AOD", "ESEQVECTPCALL", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCallRe, ese_qvec::EseQvecTPCallIm, qvec::NTrkTPCall, qvec::LabelsTPCall); + +DECLARE_SOA_TABLE(EseQvecFT0CVecs, "AOD", "ESEQVECFT0CVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0CReVec, ese_qvec::EseQvecFT0CImVec, qvec::SumAmplFT0C); +DECLARE_SOA_TABLE(EseQvecFT0AVecs, "AOD", "ESEQVECFT0AVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0AReVec, ese_qvec::EseQvecFT0AImVec, qvec::SumAmplFT0A); +DECLARE_SOA_TABLE(EseQvecFT0MVecs, "AOD", "ESEQVECFT0MVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecFT0MReVec, ese_qvec::EseQvecFT0MImVec, qvec::SumAmplFT0M); +DECLARE_SOA_TABLE(EseQvecFV0AVecs, "AOD", "ESEQVECFV0AVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecFV0AReVec, ese_qvec::EseQvecFV0AImVec, qvec::SumAmplFV0A); +DECLARE_SOA_TABLE(EseQvecTPCposVecs, "AOD", "ESEQVECTPCPVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCposReVec, ese_qvec::EseQvecTPCposImVec, qvec::NTrkTPCpos, qvec::LabelsTPCpos); +DECLARE_SOA_TABLE(EseQvecTPCnegVecs, "AOD", "ESEQVECTPCNVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCnegReVec, ese_qvec::EseQvecTPCnegImVec, qvec::NTrkTPCneg, qvec::LabelsTPCneg); +DECLARE_SOA_TABLE(EseQvecTPCallVecs, "AOD", "ESEQVECTPCAVEC", ese_qvec::IsCalibrated, ese_qvec::EseQvecTPCallReVec, ese_qvec::EseQvecTPCallImVec, qvec::NTrkTPCall, qvec::LabelsTPCall); + +DECLARE_SOA_TABLE(EseQvecPercs, "AOD", "ESEQVECPERC", ese_qvec::EseQvecFT0CRe, ese_qvec::EseQvecFT0CIm, qvec::SumAmplFT0C, + ese_qvec::EseQvecFT0ARe, ese_qvec::EseQvecFT0AIm, qvec::SumAmplFT0A, + ese_qvec::EseQvecFT0MRe, ese_qvec::EseQvecFT0MIm, qvec::SumAmplFT0M, + ese_qvec::EseQvecFV0ARe, ese_qvec::EseQvecFV0AIm, qvec::SumAmplFV0A, + ese_qvec::EseQvecTPCposRe, ese_qvec::EseQvecTPCposIm, qvec::NTrkTPCpos, + ese_qvec::EseQvecTPCnegRe, ese_qvec::EseQvecTPCnegIm, qvec::NTrkTPCneg, + ese_qvec::EseQvecTPCallRe, ese_qvec::EseQvecTPCallIm, qvec::NTrkTPCall); + +using EseQvectorFT0C = EseQvecFT0Cs::iterator; +using EseQvectorFT0A = EseQvecFT0As::iterator; +using EseQvectorFT0M = EseQvecFT0Ms::iterator; +using EseQvectorFV0A = EseQvecFV0As::iterator; +using EseQvectorTPCpos = EseQvecTPCposs::iterator; +using EseQvectorTPCneg = EseQvecTPCnegs::iterator; +using EseQvectorTPCall = EseQvecTPCalls::iterator; + +using EseQvectorFT0CVec = EseQvecFT0CVecs::iterator; +using EseQvectorFT0AVec = EseQvecFT0AVecs::iterator; +using EseQvectorFT0MVec = EseQvecFT0MVecs::iterator; +using EseQvectorFV0AVec = EseQvecFV0AVecs::iterator; +using EseQvectorTPCposVec = EseQvecTPCposVecs::iterator; +using EseQvectorTPCnegVec = EseQvecTPCnegVecs::iterator; +using EseQvectorTPCallVec = EseQvecTPCallVecs::iterator; +///////////////////////////////////////////////////////////////// } // namespace o2::aod #endif // COMMON_DATAMODEL_QVECTORS_H_ diff --git a/Common/TableProducer/qVectorsTable.cxx b/Common/TableProducer/qVectorsTable.cxx index 273394a8e67..574c370a16a 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -50,6 +50,7 @@ #include #include +#include #include #include #include @@ -65,14 +66,28 @@ using MyCollisions = soa::Join; struct qVectorsTable { - enum { + enum Detectors { kFT0C = 0, kFT0A = 1, kFT0M, kFV0A, - kTPCpos, - kTPCneg, - kTPCall + kTPCPos, + kTPCNeg, + kTPCAll, + kNDetectors + }; + enum Corrections { + kNoCorr = 0, + kRecenter, + kTwist, + kRescale, + kNCorrections + }; + enum MultNorms { + kNoNorm = 0, + kScalarProd, + kEsE, + kMultNormTypes }; // Configurables. @@ -98,7 +113,7 @@ struct qVectorsTable { Configurable useCorrectionForRun{"useCorrectionForRun", true, "Get Qvector corrections based on run number instead of timestamp"}; Configurable cfgGainEqPath{"cfgGainEqPath", "Users/j/junlee/Qvector/GainEq", "CCDB path for gain equalization constants"}; - Configurable cfgQvecCalibPath{"cfgQvecCalibPath", "Analysis/EventPlane/QVecCorrections", "CCDB pasth for Q-vecteor calibration constants"}; + Configurable cfgQvecCalibPath{"cfgQvecCalibPath", "Analysis/EventPlane/QVecCorrections", "CCDB path for Q-vector calibration constants"}; Configurable cfgShiftCorr{"cfgShiftCorr", false, "configurable flag for shift correction"}; Configurable cfgShiftPath{"cfgShiftPath", "", "CCDB path for shift correction"}; @@ -109,9 +124,10 @@ struct qVectorsTable { Configurable cfgUseFT0A{"cfgUseFT0A", false, "Initial value for using FT0A. By default obtained from DataModel."}; Configurable cfgUseFT0M{"cfgUseFT0M", false, "Initial value for using FT0M. By default obtained from DataModel."}; Configurable cfgUseFV0A{"cfgUseFV0A", false, "Initial value for using FV0A. By default obtained from DataModel."}; - Configurable cfgUseTPCpos{"cfgUseTPCpos", false, "Initial value for using TPCpos. By default obtained from DataModel."}; - Configurable cfgUseTPCneg{"cfgUseTPCneg", false, "Initial value for using TPCneg. By default obtained from DataModel."}; - Configurable cfgUseTPCall{"cfgUseTPCall", false, "Initial value for using TPCall. By default obtained from DataModel."}; + Configurable cfgUseTPCpos{"cfgUseTPCpos", false, "Initial value for using TPCPos. By default obtained from DataModel."}; + Configurable cfgUseTPCneg{"cfgUseTPCneg", false, "Initial value for using TPCNeg. By default obtained from DataModel."}; + Configurable cfgUseTPCall{"cfgUseTPCall", false, "Initial value for using TPCAll. By default obtained from DataModel."}; + Configurable cfgProduceRedQVecs{"cfgProduceRedQVecs", false, "Produce reduced Q-vectors for Event-Shape Engineering"}; // Table. Produces qVector; @@ -119,17 +135,35 @@ struct qVectorsTable { Produces qVectorFT0A; Produces qVectorFT0M; Produces qVectorFV0A; - Produces qVectorTPCpos; - Produces qVectorTPCneg; - Produces qVectorTPCall; + Produces qVectorTPCPos; + Produces qVectorTPCNeg; + Produces qVectorTPCAll; Produces qVectorFT0CVec; Produces qVectorFT0AVec; Produces qVectorFT0MVec; Produces qVectorFV0AVec; - Produces qVectorTPCposVec; - Produces qVectorTPCnegVec; - Produces qVectorTPCallVec; + Produces qVectorTPCPosVec; + Produces qVectorTPCNegVec; + Produces qVectorTPCAllVec; + + Produces eseQVector; + Produces eseQVectorFT0C; + Produces eseQVectorFT0A; + Produces eseQVectorFT0M; + Produces eseQVectorFV0A; + Produces eseQVectorTPCPos; + Produces eseQVectorTPCNeg; + Produces eseQVectorTPCAll; + + Produces eseQVectorFT0CVec; + Produces eseQVectorFT0AVec; + Produces eseQVectorFT0MVec; + Produces eseQVectorFV0AVec; + Produces eseQVectorTPCPosVec; + Produces eseQVectorTPCNegVec; + Produces eseQVectorTPCAllVec; + Produces eseQVectorPerc; std::vector FT0RelGainConst{}; std::vector FV0RelGainConst{}; @@ -147,11 +181,14 @@ struct qVectorsTable { HistogramRegistry histosQA{"histosQA", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + const float minAmplitude = 1.0e-8f; int runNumber{-1}; float cent; - std::vector objQvec{}; - std::vector shiftprofile{}; + std::vector corrsQvecSp{}; + std::vector corrsQvecEse{}; + std::vector shiftProfileSp{}; + std::vector shiftProfileEse{}; // Deprecated, will be removed in future after transition time // Configurable cfgUseBPos{"cfgUseBPos", false, "Initial value for using BPos. By default obtained from DataModel."}; @@ -193,8 +230,8 @@ struct qVectorsTable { goto allDetectorsInUse; // Added to break from nested loop if all detectors are in use. } for (auto const& det : useDetector) { - std::string table_name_with_vector = det.first; // for replacing s with Vecs at the end. - if (input.matcher.binding == det.first || input.matcher.binding == table_name_with_vector.replace(table_name_with_vector.size() - 1, 1, "Vecs")) { + std::string tableNameWithVector = det.first; // for replacing s with Vecs at the end. + if (input.matcher.binding == det.first || input.matcher.binding == tableNameWithVector.replace(tableNameWithVector.size() - 1, 1, "Vecs")) { useDetector[det.first.data()] = true; LOGF(info, Form("Using detector: %s.", det.first.data())); } @@ -257,38 +294,66 @@ struct qVectorsTable { LOGF(fatal, "Could not get the alignment parameters for FV0."); } - objQvec.clear(); + corrsQvecSp.clear(); for (std::size_t i = 0; i < cfgnMods->size(); i++) { int ind = cfgnMods->at(i); fullPath = cfgQvecCalibPath; fullPath += "/v"; fullPath += std::to_string(ind); - auto objqvec = getForTsOrRun(fullPath, timestamp, runnumber); - if (!objqvec) { + auto modeCorrQvecSp = getForTsOrRun(fullPath, timestamp, runnumber); + if (!modeCorrQvecSp) { fullPath = cfgQvecCalibPath; fullPath += "/v2"; - objqvec = getForTsOrRun(fullPath, timestamp, runnumber); + modeCorrQvecSp = getForTsOrRun(fullPath, timestamp, runnumber); } - objQvec.push_back(objqvec); + corrsQvecSp.push_back(modeCorrQvecSp); + } + + corrsQvecEse.clear(); + for (std::size_t i = 0; i < cfgnMods->size(); i++) { + int ind = cfgnMods->at(i); + fullPath = cfgQvecCalibPath; + fullPath += "/eseq"; + fullPath += std::to_string(ind); + auto modeCorrQvecEse = getForTsOrRun(fullPath, timestamp, runnumber); + if (!modeCorrQvecEse) { + fullPath = cfgQvecCalibPath; + fullPath += "/eseq2"; + modeCorrQvecEse = getForTsOrRun(fullPath, timestamp, runnumber); + } + corrsQvecEse.push_back(modeCorrQvecEse); } if (cfgShiftCorr) { - shiftprofile.clear(); + shiftProfileSp.clear(); for (std::size_t i = 0; i < cfgnMods->size(); i++) { int ind = cfgnMods->at(i); fullPath = cfgShiftPath; fullPath += "/v"; fullPath += std::to_string(ind); auto objshift = getForTsOrRun(fullPath, timestamp, runnumber); - shiftprofile.push_back(objshift); + shiftProfileSp.push_back(objshift); + } + + if (cfgProduceRedQVecs) { + shiftProfileEse.clear(); + for (std::size_t i = 0; i < cfgnMods->size(); i++) { + int ind = cfgnMods->at(i); + fullPath = cfgShiftPath; + fullPath += "/eseq"; + fullPath += std::to_string(ind); + auto objshift = getForTsOrRun(fullPath, timestamp, runnumber); + shiftProfileEse.push_back(objshift); + } } } fullPath = cfgGainEqPath; fullPath += "/FT0"; + const int nPixelsFT0 = 208; const auto objft0Gain = getForTsOrRun>(fullPath, timestamp, runnumber); if (!objft0Gain || cfgCorrLevel == 0) { - for (auto i{0u}; i < 208; i++) { + for (auto i{0u}; i < nPixelsFT0; i++) { FT0RelGainConst.push_back(1.); } } else { @@ -297,9 +362,10 @@ struct qVectorsTable { fullPath = cfgGainEqPath; fullPath += "/FV0"; + const int nChannelsFV0 = 48; const auto objfv0Gain = getForTsOrRun>(fullPath, timestamp, runnumber); if (!objfv0Gain || cfgCorrLevel == 0) { - for (auto i{0u}; i < 48; i++) { + for (auto i{0u}; i < nChannelsFV0; i++) { FV0RelGainConst.push_back(1.); } } else { @@ -308,7 +374,7 @@ struct qVectorsTable { } template - bool SelTrack(const TrackType track) + bool selTrack(const TrackType track) { if (track.pt() < cfgMinPtOnTPC) return false; @@ -347,19 +413,177 @@ struct qVectorsTable { } } + /// Function to normalize the q-vectors + /// \param qVecReNorm is the vector with the normalized real part of the q-vector for each detector and correction step + /// \param qVecImNorm is the vector with the normalized imaginary part of the q-vector for each detector and correction step + /// \param qVecReRaw is the vector with the raw real part of the q-vector for each detector and correction step + /// \param qVecImRaw is the vector with the raw imaginary part of the q-vector for each detector and correction step + /// \param qVecAmp is the vector with the amplitude of the q-vector for each detector and correction step + /// \param normType is the type of normalization to apply to the q-vectors + void normalizeQVec(std::vector& qVecReNorm, + std::vector& qVecImNorm, + std::vector const& qVecReRaw, + std::vector const& qVecImRaw, + std::vector const& qVecAmp, + MultNorms normType) + { + for (std::size_t i = 0; i < kNDetectors; i++) { + float qVecDetReNorm{-999.}, qVecDetImNorm{-999.}; + if (qVecAmp[i] > minAmplitude) { + switch (normType) { + case MultNorms::kScalarProd: + qVecDetReNorm = qVecReRaw[i] / qVecAmp[i]; + qVecDetImNorm = qVecImRaw[i] / qVecAmp[i]; + break; + case MultNorms::kEsE: + qVecDetReNorm = qVecReRaw[i] / std::sqrt(qVecAmp[i]); + qVecDetImNorm = qVecImRaw[i] / std::sqrt(qVecAmp[i]); + break; + default: + LOGP(fatal, "Undefined normalization type for Q-vector amplitude. Check the configuration."); + break; + } + } + for (int iCorr = 0; iCorr < kNCorrections; iCorr++) { + qVecReNorm.push_back(qVecDetReNorm); + qVecImNorm.push_back(qVecDetImNorm); + } + } + } + + /// Function to calculate the un-normalized q-vectors + /// \param cent is the collision centrality + /// \param qVecRe is the vector with the real part of the q-vector for each detector and correction step + /// \param qVecIm is the vector with the imaginary part of the q-vector for each detector and correction step + /// \param histsCorrs is the vector with the histograms with the correction constants for each detector and correction step + /// \param nMode is the modulation of interest + void correctQVec(float cent, std::vector& qVecRe, std::vector& qVecIm, TH3F* histsCorrs, std::vector& shiftProfile, int nMode) + { + int nCorrections = static_cast(kNCorrections); + if (cent < cfgMaxCentrality) { + for (auto i{0u}; i < kTPCAll + 1; i++) { + int idxDet = i * kNCorrections; + helperEP.DoRecenter(qVecRe[idxDet + kRecenter], qVecIm[idxDet + kRecenter], + histsCorrs->GetBinContent(static_cast(cent) + 1, 1, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 2, i + 1)); + + helperEP.DoRecenter(qVecRe[idxDet + kTwist], qVecIm[idxDet + kTwist], + histsCorrs->GetBinContent(static_cast(cent) + 1, 1, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 2, i + 1)); + helperEP.DoTwist(qVecRe[idxDet + kTwist], qVecIm[idxDet + kTwist], + histsCorrs->GetBinContent(static_cast(cent) + 1, 3, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 4, i + 1)); + + helperEP.DoRecenter(qVecRe[idxDet + kRescale], qVecIm[idxDet + kRescale], + histsCorrs->GetBinContent(static_cast(cent) + 1, 1, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 2, i + 1)); + helperEP.DoTwist(qVecRe[idxDet + kRescale], qVecIm[idxDet + kRescale], + histsCorrs->GetBinContent(static_cast(cent) + 1, 3, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 4, i + 1)); + helperEP.DoRescale(qVecRe[idxDet + kRescale], qVecIm[idxDet + kRescale], + histsCorrs->GetBinContent(static_cast(cent) + 1, 5, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 6, i + 1)); + } + if (cfgShiftCorr) { + auto deltaPsiFT0C = 0.0; + auto deltaPsiFT0A = 0.0; + auto deltaPsiFT0M = 0.0; + auto deltaPsiFV0A = 0.0; + auto deltaPsiTPCPos = 0.0; + auto deltaPsiTPCNeg = 0.0; + auto deltaPsiTPCAll = 0.0; + + auto psiDefFT0C = std::atan2(qVecIm[kFT0C * nCorrections + kRescale], qVecRe[kFT0C * nCorrections + kRescale]) / static_cast(nMode); + auto psiDefFT0A = std::atan2(qVecIm[kFT0A * nCorrections + kRescale], qVecRe[kFT0A * nCorrections + kRescale]) / static_cast(nMode); + auto psiDefFT0M = std::atan2(qVecIm[kFT0M * nCorrections + kRescale], qVecRe[kFT0M * nCorrections + kRescale]) / static_cast(nMode); + auto psiDefFV0A = std::atan2(qVecIm[kFV0A * nCorrections + kRescale], qVecRe[kFV0A * nCorrections + kRescale]) / static_cast(nMode); + auto psiDefTPCPos = std::atan2(qVecIm[kTPCPos * nCorrections + kRescale], qVecRe[kTPCPos * nCorrections + kRescale]) / static_cast(nMode); + auto psiDefTPCNeg = std::atan2(qVecIm[kTPCNeg * nCorrections + kRescale], qVecRe[kTPCNeg * nCorrections + kRescale]) / static_cast(nMode); + auto psiDefTPCAll = std::atan2(qVecIm[kTPCAll * nCorrections + kRescale], qVecRe[kTPCAll * nCorrections + kRescale]) / static_cast(nMode); + + for (int iShift = 1; iShift <= 10; iShift++) { + auto coeffShiftXFT0C = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kFT0C, iShift - 0.5)); + auto coeffShiftYFT0C = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kFT0C + 1, iShift - 0.5)); + auto coeffShiftXFT0A = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kFT0A, iShift - 0.5)); + auto coeffShiftYFT0A = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kFT0A + 1, iShift - 0.5)); + auto coeffShiftXFT0M = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kFT0M, iShift - 0.5)); + auto coeffShiftYFT0M = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kFT0M + 1, iShift - 0.5)); + auto coeffShiftXFV0A = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kFV0A, iShift - 0.5)); + auto coeffShiftYFV0A = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kFV0A + 1, iShift - 0.5)); + auto coeffShiftXTPCPos = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kTPCPos, iShift - 0.5)); + auto coeffShiftYTPCPos = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kTPCPos + 1, iShift - 0.5)); + auto coeffShiftXTPCNeg = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kTPCNeg, iShift - 0.5)); + auto coeffShiftYTPCNeg = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kTPCNeg + 1, iShift - 0.5)); + auto coeffShiftXTPCAll = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kTPCAll, iShift - 0.5)); + auto coeffShiftYTPCAll = shiftProfile.at(nMode - 2)->GetBinContent(shiftProfile.at(nMode - 2)->FindBin(cent, 2 * kTPCAll + 1, iShift - 0.5)); + + deltaPsiFT0C += ((2. / (1.0 * iShift)) * (-coeffShiftXFT0C * std::cos(iShift * static_cast(nMode) * psiDefFT0C) + coeffShiftYFT0C * std::sin(iShift * static_cast(nMode) * psiDefFT0C))) / static_cast(nMode); + deltaPsiFT0A += ((2. / (1.0 * iShift)) * (-coeffShiftXFT0A * std::cos(iShift * static_cast(nMode) * psiDefFT0A) + coeffShiftYFT0A * std::sin(iShift * static_cast(nMode) * psiDefFT0A))) / static_cast(nMode); + deltaPsiFT0M += ((2. / (1.0 * iShift)) * (-coeffShiftXFT0M * std::cos(iShift * static_cast(nMode) * psiDefFT0M) + coeffShiftYFT0M * std::sin(iShift * static_cast(nMode) * psiDefFT0M))) / static_cast(nMode); + deltaPsiFV0A += ((2. / (1.0 * iShift)) * (-coeffShiftXFV0A * std::cos(iShift * static_cast(nMode) * psiDefFV0A) + coeffShiftYFV0A * std::sin(iShift * static_cast(nMode) * psiDefFV0A))) / static_cast(nMode); + deltaPsiTPCPos += ((2. / (1.0 * iShift)) * (-coeffShiftXTPCPos * std::cos(iShift * static_cast(nMode) * psiDefTPCPos) + coeffShiftYTPCPos * std::sin(iShift * static_cast(nMode) * psiDefTPCPos))) / static_cast(nMode); + deltaPsiTPCNeg += ((2. / (1.0 * iShift)) * (-coeffShiftXTPCNeg * std::cos(iShift * static_cast(nMode) * psiDefTPCNeg) + coeffShiftYTPCNeg * std::sin(iShift * static_cast(nMode) * psiDefTPCNeg))) / static_cast(nMode); + deltaPsiTPCAll += ((2. / (1.0 * iShift)) * (-coeffShiftXTPCAll * std::cos(iShift * static_cast(nMode) * psiDefTPCAll) + coeffShiftYTPCAll * std::sin(iShift * static_cast(nMode) * psiDefTPCAll))) / static_cast(nMode); + } + + deltaPsiFT0C *= static_cast(nMode); + deltaPsiFT0A *= static_cast(nMode); + deltaPsiFT0M *= static_cast(nMode); + deltaPsiFV0A *= static_cast(nMode); + deltaPsiTPCPos *= static_cast(nMode); + deltaPsiTPCNeg *= static_cast(nMode); + deltaPsiTPCAll *= static_cast(nMode); + + float qVecReShiftedFT0C = qVecRe[kFT0C * nCorrections + kRescale] * std::cos(deltaPsiFT0C) - qVecIm[kFT0C * nCorrections + kRescale] * std::sin(deltaPsiFT0C); + float qVecImShiftedFT0C = qVecRe[kFT0C * nCorrections + kRescale] * std::sin(deltaPsiFT0C) + qVecIm[kFT0C * nCorrections + kRescale] * std::cos(deltaPsiFT0C); + float qVecReShiftedFT0A = qVecRe[kFT0A * nCorrections + kRescale] * std::cos(deltaPsiFT0A) - qVecIm[kFT0A * nCorrections + kRescale] * std::sin(deltaPsiFT0A); + float qVecImShiftedFT0A = qVecRe[kFT0A * nCorrections + kRescale] * std::sin(deltaPsiFT0A) + qVecIm[kFT0A * nCorrections + kRescale] * std::cos(deltaPsiFT0A); + float qVecReShiftedFT0M = qVecRe[kFT0M * nCorrections + kRescale] * std::cos(deltaPsiFT0M) - qVecIm[kFT0M * nCorrections + kRescale] * std::sin(deltaPsiFT0M); + float qVecImShiftedFT0M = qVecRe[kFT0M * nCorrections + kRescale] * std::sin(deltaPsiFT0M) + qVecIm[kFT0M * nCorrections + kRescale] * std::cos(deltaPsiFT0M); + float qVecReShiftedFV0A = qVecRe[kFV0A * nCorrections + kRescale] * std::cos(deltaPsiFV0A) - qVecIm[kFV0A * nCorrections + kRescale] * std::sin(deltaPsiFV0A); + float qVecImShiftedFV0A = qVecRe[kFV0A * nCorrections + kRescale] * std::sin(deltaPsiFV0A) + qVecIm[kFV0A * nCorrections + kRescale] * std::cos(deltaPsiFV0A); + float qVecReShiftedTPCPos = qVecRe[kTPCPos * nCorrections + kRescale] * std::cos(deltaPsiTPCPos) - qVecIm[kTPCPos * nCorrections + kRescale] * std::sin(deltaPsiTPCPos); + float qVecImShiftedTPCPos = qVecRe[kTPCPos * nCorrections + kRescale] * std::sin(deltaPsiTPCPos) + qVecIm[kTPCPos * nCorrections + kRescale] * std::cos(deltaPsiTPCPos); + float qVecReShiftedTPCNeg = qVecRe[kTPCNeg * nCorrections + kRescale] * std::cos(deltaPsiTPCNeg) - qVecIm[kTPCNeg * nCorrections + kRescale] * std::sin(deltaPsiTPCNeg); + float qVecImShiftedTPCNeg = qVecRe[kTPCNeg * nCorrections + kRescale] * std::sin(deltaPsiTPCNeg) + qVecIm[kTPCNeg * nCorrections + kRescale] * std::cos(deltaPsiTPCNeg); + float qVecReShiftedTPCAll = qVecRe[kTPCAll * nCorrections + kRescale] * std::cos(deltaPsiTPCAll) - qVecIm[kTPCAll * nCorrections + kRescale] * std::sin(deltaPsiTPCAll); + float qVecImShiftedTPCAll = qVecRe[kTPCAll * nCorrections + kRescale] * std::sin(deltaPsiTPCAll) + qVecIm[kTPCAll * nCorrections + kRescale] * std::cos(deltaPsiTPCAll); + + qVecRe[kFT0C * nCorrections + kRescale] = qVecReShiftedFT0C; + qVecIm[kFT0C * nCorrections + kRescale] = qVecImShiftedFT0C; + qVecRe[kFT0A * nCorrections + kRescale] = qVecReShiftedFT0A; + qVecIm[kFT0A * nCorrections + kRescale] = qVecImShiftedFT0A; + qVecRe[kFT0M * nCorrections + kRescale] = qVecReShiftedFT0M; + qVecIm[kFT0M * nCorrections + kRescale] = qVecImShiftedFT0M; + qVecRe[kFV0A * nCorrections + kRescale] = qVecReShiftedFV0A; + qVecIm[kFV0A * nCorrections + kRescale] = qVecImShiftedFV0A; + qVecRe[kTPCPos * nCorrections + kRescale] = qVecReShiftedTPCPos; + qVecIm[kTPCPos * nCorrections + kRescale] = qVecImShiftedTPCPos; + qVecRe[kTPCNeg * nCorrections + kRescale] = qVecReShiftedTPCNeg; + qVecIm[kTPCNeg * nCorrections + kRescale] = qVecImShiftedTPCNeg; + qVecRe[kTPCAll * nCorrections + kRescale] = qVecReShiftedTPCAll; + qVecIm[kTPCAll * nCorrections + kRescale] = qVecImShiftedTPCAll; + } + } + } + + /// Function to calculate the un-normalized q-vectors + /// \param nMode is the harmonic number of the q-vector + /// \param coll is the collision object + /// \param tracks are the tracks associated to the collision + /// \param qVecRe is the vector with the real part of the q-vector for each detector + /// \param qVecIm is the vector with the imaginary part of the q-vector for each detector + /// \param qVecAmp is the vector with the amplitude of the signal in each detector + /// \param trkTPCPosLabel is the vector with the number of TPC tracks with positive eta + /// \param trkTPCNegLabel is the vector with the number of TPC tracks with negative eta + /// \param trkTPCAllLabel is the vector with the number of TPC tracks with any eta template - void CalQvec(const Nmode nmode, const CollType& coll, const TrackType& track, std::vector& QvecRe, std::vector& QvecIm, std::vector& QvecAmp, std::vector& TrkTPCposLabel, std::vector& TrkTPCnegLabel, std::vector& TrkTPCallLabel) + void calcQVec(const Nmode nMode, const CollType& coll, const TrackType& tracks, std::vector& qVecRe, std::vector& qVecIm, std::vector& qVecAmp, std::vector& trkTPCPosLabel, std::vector& trkTPCNegLabel, std::vector& trkTPCAllLabel) { - float qVectFT0A[2] = {0.}; - float qVectFT0C[2] = {0.}; - float qVectFT0M[2] = {0.}; - float qVectFV0A[2] = {0.}; - float qVectTPCpos[2] = {0.}; - float qVectTPCneg[2] = {0.}; - float qVectTPCall[2] = {0.}; - - TComplex QvecDet(0); - TComplex QvecFT0M(0); + float qVectFT0A[2] = {-999., -999.}; + float qVectFT0C[2] = {-999., -999.}; + float qVectFT0M[2] = {-999., -999.}; + float qVectFV0A[2] = {-999., -999.}; + float qVectTPCPos[2] = {0., 0.}; // Always computed + float qVectTPCNeg[2] = {0., 0.}; // Always computed + float qVectTPCAll[2] = {0., 0.}; // Always computed + + TComplex qVecDet(0); + TComplex qVecFT0M(0); float sumAmplFT0A = 0.; float sumAmplFT0C = 0.; float sumAmplFT0M = 0.; @@ -376,21 +600,17 @@ struct qVectorsTable { histosQA.fill(HIST("FT0Amp"), ampl, FT0AchId); histosQA.fill(HIST("FT0AmpCor"), ampl / FT0RelGainConst[FT0AchId], FT0AchId); - helperEP.SumQvectors(0, FT0AchId, ampl / FT0RelGainConst[FT0AchId], nmode, QvecDet, sumAmplFT0A, ft0geom, fv0geom); - helperEP.SumQvectors(0, FT0AchId, ampl / FT0RelGainConst[FT0AchId], nmode, QvecFT0M, sumAmplFT0M, ft0geom, fv0geom); + helperEP.SumQvectors(0, FT0AchId, ampl / FT0RelGainConst[FT0AchId], nMode, qVecDet, sumAmplFT0A, ft0geom, fv0geom); + helperEP.SumQvectors(0, FT0AchId, ampl / FT0RelGainConst[FT0AchId], nMode, qVecFT0M, sumAmplFT0M, ft0geom, fv0geom); } - if (sumAmplFT0A > 1e-8) { - QvecDet /= sumAmplFT0A; - qVectFT0A[0] = QvecDet.Re(); - qVectFT0A[1] = QvecDet.Im(); + if (sumAmplFT0A > minAmplitude) { + qVectFT0A[0] = qVecDet.Re(); + qVectFT0A[1] = qVecDet.Im(); } - } else { - qVectFT0A[0] = 999.; - qVectFT0A[1] = 999.; } if (useDetector["QvectorFT0Cs"]) { - QvecDet = TComplex(0., 0.); + qVecDet = TComplex(0., 0.); for (std::size_t iChC = 0; iChC < ft0.channelC().size(); iChC++) { float ampl = ft0.amplitudeC()[iChC]; int FT0CchId = ft0.channelC()[iChC] + 96; @@ -398,73 +618,47 @@ struct qVectorsTable { histosQA.fill(HIST("FT0Amp"), ampl, FT0CchId); histosQA.fill(HIST("FT0AmpCor"), ampl / FT0RelGainConst[FT0CchId], FT0CchId); - helperEP.SumQvectors(0, FT0CchId, ampl / FT0RelGainConst[FT0CchId], nmode, QvecDet, sumAmplFT0C, ft0geom, fv0geom); - helperEP.SumQvectors(0, FT0CchId, ampl / FT0RelGainConst[FT0CchId], nmode, QvecFT0M, sumAmplFT0M, ft0geom, fv0geom); + helperEP.SumQvectors(0, FT0CchId, ampl / FT0RelGainConst[FT0CchId], nMode, qVecDet, sumAmplFT0C, ft0geom, fv0geom); + helperEP.SumQvectors(0, FT0CchId, ampl / FT0RelGainConst[FT0CchId], nMode, qVecFT0M, sumAmplFT0M, ft0geom, fv0geom); } - if (sumAmplFT0C > 1e-8) { - QvecDet /= sumAmplFT0C; - qVectFT0C[0] = QvecDet.Re(); - qVectFT0C[1] = QvecDet.Im(); - } else { - qVectFT0C[0] = 999.; - qVectFT0C[1] = 999.; + if (sumAmplFT0C > minAmplitude) { + qVectFT0C[0] = qVecDet.Re(); + qVectFT0C[1] = qVecDet.Im(); + } + if (sumAmplFT0M > minAmplitude && useDetector["QvectorFT0Ms"]) { + qVectFT0M[0] = qVecFT0M.Re(); + qVectFT0M[1] = qVecFT0M.Im(); } - } else { - qVectFT0C[0] = -999.; - qVectFT0C[1] = -999.; - } - - if (sumAmplFT0M > 1e-8 && useDetector["QvectorFT0Ms"]) { - QvecFT0M /= sumAmplFT0M; - qVectFT0M[0] = QvecFT0M.Re(); - qVectFT0M[1] = QvecFT0M.Im(); - } else { - qVectFT0M[0] = 999.; - qVectFT0M[1] = 999.; } - } else { - qVectFT0A[0] = -999.; - qVectFT0A[1] = -999.; - qVectFT0C[0] = -999.; - qVectFT0C[1] = -999.; - qVectFT0M[0] = -999.; - qVectFT0M[1] = -999.; - } - QvecDet = TComplex(0., 0.); - sumAmplFV0A = 0; - if (coll.has_foundFV0() && useDetector["QvectorFV0As"]) { - auto fv0 = coll.foundFV0(); + qVecDet = TComplex(0., 0.); + sumAmplFV0A = 0; + if (coll.has_foundFV0() && useDetector["QvectorFV0As"]) { + auto fv0 = coll.foundFV0(); - for (std::size_t iCh = 0; iCh < fv0.channel().size(); iCh++) { - float ampl = fv0.amplitude()[iCh]; - int FV0AchId = fv0.channel()[iCh]; - histosQA.fill(HIST("FV0Amp"), ampl, FV0AchId); - histosQA.fill(HIST("FV0AmpCor"), ampl / FV0RelGainConst[FV0AchId], FV0AchId); + for (std::size_t iCh = 0; iCh < fv0.channel().size(); iCh++) { + float ampl = fv0.amplitude()[iCh]; + int FV0AchId = fv0.channel()[iCh]; + histosQA.fill(HIST("FV0Amp"), ampl, FV0AchId); + histosQA.fill(HIST("FV0AmpCor"), ampl / FV0RelGainConst[FV0AchId], FV0AchId); - helperEP.SumQvectors(1, FV0AchId, ampl / FV0RelGainConst[FV0AchId], nmode, QvecDet, sumAmplFV0A, ft0geom, fv0geom); - } + helperEP.SumQvectors(1, FV0AchId, ampl / FV0RelGainConst[FV0AchId], nMode, qVecDet, sumAmplFV0A, ft0geom, fv0geom); + } - if (sumAmplFV0A > 1e-8) { - QvecDet /= sumAmplFV0A; - qVectFV0A[0] = QvecDet.Re(); - qVectFV0A[1] = QvecDet.Im(); - } else { - qVectFV0A[0] = 999.; - qVectFV0A[1] = 999.; + if (sumAmplFV0A > minAmplitude) { + qVectFV0A[0] = qVecDet.Re(); + qVectFV0A[1] = qVecDet.Im(); + } } - } else { - qVectFV0A[0] = -999.; - qVectFV0A[1] = -999.; } - int nTrkTPCpos = 0; - int nTrkTPCneg = 0; - int nTrkTPCall = 0; + int nTrkTPCPos = 0; + int nTrkTPCNeg = 0; + int nTrkTPCAll = 0; - for (auto const& trk : track) { - if (!SelTrack(trk)) { + for (auto const& trk : tracks) { + if (!selTrack(trk)) { continue; } histosQA.fill(HIST("ChTracks"), trk.pt(), trk.eta(), trk.phi(), cent); @@ -474,110 +668,90 @@ struct qVectorsTable { if (trk.eta() < cfgEtaMin) { continue; } - qVectTPCall[0] += trk.pt() * std::cos(trk.phi() * nmode); - qVectTPCall[1] += trk.pt() * std::sin(trk.phi() * nmode); - TrkTPCallLabel.push_back(trk.globalIndex()); - nTrkTPCall++; + qVectTPCAll[0] += trk.pt() * std::cos(trk.phi() * nMode); + qVectTPCAll[1] += trk.pt() * std::sin(trk.phi() * nMode); + trkTPCAllLabel.push_back(trk.globalIndex()); + nTrkTPCAll++; if (std::abs(trk.eta()) < 0.1) { continue; } if (trk.eta() > 0 && (useDetector["QvectorTPCposs"] || useDetector["QvectorBPoss"])) { - qVectTPCpos[0] += trk.pt() * std::cos(trk.phi() * nmode); - qVectTPCpos[1] += trk.pt() * std::sin(trk.phi() * nmode); - TrkTPCposLabel.push_back(trk.globalIndex()); - nTrkTPCpos++; + qVectTPCPos[0] += trk.pt() * std::cos(trk.phi() * nMode); + qVectTPCPos[1] += trk.pt() * std::sin(trk.phi() * nMode); + trkTPCPosLabel.push_back(trk.globalIndex()); + nTrkTPCPos++; } else if (trk.eta() < 0 && (useDetector["QvectorTPCnegs"] || useDetector["QvectorBNegs"])) { - qVectTPCneg[0] += trk.pt() * std::cos(trk.phi() * nmode); - qVectTPCneg[1] += trk.pt() * std::sin(trk.phi() * nmode); - TrkTPCnegLabel.push_back(trk.globalIndex()); - nTrkTPCneg++; + qVectTPCNeg[0] += trk.pt() * std::cos(trk.phi() * nMode); + qVectTPCNeg[1] += trk.pt() * std::sin(trk.phi() * nMode); + trkTPCNegLabel.push_back(trk.globalIndex()); + nTrkTPCNeg++; } } - if (nTrkTPCpos > 0) { - qVectTPCpos[0] /= nTrkTPCpos; - qVectTPCpos[1] /= nTrkTPCpos; - } else { - qVectTPCpos[0] = 999.; - qVectTPCpos[1] = 999.; - } - - if (nTrkTPCneg > 0) { - qVectTPCneg[0] /= nTrkTPCneg; - qVectTPCneg[1] /= nTrkTPCneg; - } else { - qVectTPCneg[0] = 999.; - qVectTPCneg[1] = 999.; - } - - if (nTrkTPCall > 0) { - qVectTPCall[0] /= nTrkTPCall; - qVectTPCall[1] /= nTrkTPCall; - } else { - qVectTPCall[0] = 999.; - qVectTPCall[1] = 999.; - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectFT0C[0]); - QvecIm.push_back(qVectFT0C[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectFT0A[0]); - QvecIm.push_back(qVectFT0A[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectFT0M[0]); - QvecIm.push_back(qVectFT0M[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectFV0A[0]); - QvecIm.push_back(qVectFV0A[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectTPCpos[0]); - QvecIm.push_back(qVectTPCpos[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectTPCneg[0]); - QvecIm.push_back(qVectTPCneg[1]); - } - for (auto i{0u}; i < 4; i++) { - QvecRe.push_back(qVectTPCall[0]); - QvecIm.push_back(qVectTPCall[1]); - } - - QvecAmp.push_back(sumAmplFT0C); - QvecAmp.push_back(sumAmplFT0A); - QvecAmp.push_back(sumAmplFT0M); - QvecAmp.push_back(sumAmplFV0A); - QvecAmp.push_back(static_cast(nTrkTPCpos)); - QvecAmp.push_back(static_cast(nTrkTPCneg)); - QvecAmp.push_back(static_cast(nTrkTPCall)); + qVecRe.push_back(qVectFT0C[0]); + qVecIm.push_back(qVectFT0C[1]); + qVecRe.push_back(qVectFT0A[0]); + qVecIm.push_back(qVectFT0A[1]); + qVecRe.push_back(qVectFT0M[0]); + qVecIm.push_back(qVectFT0M[1]); + qVecRe.push_back(qVectFV0A[0]); + qVecIm.push_back(qVectFV0A[1]); + qVecRe.push_back(qVectTPCPos[0]); + qVecIm.push_back(qVectTPCPos[1]); + qVecRe.push_back(qVectTPCNeg[0]); + qVecIm.push_back(qVectTPCNeg[1]); + qVecRe.push_back(qVectTPCAll[0]); + qVecIm.push_back(qVectTPCAll[1]); + + qVecAmp.push_back(sumAmplFT0C); + qVecAmp.push_back(sumAmplFT0A); + qVecAmp.push_back(sumAmplFT0M); + qVecAmp.push_back(sumAmplFV0A); + qVecAmp.push_back(static_cast(nTrkTPCPos)); + qVecAmp.push_back(static_cast(nTrkTPCNeg)); + qVecAmp.push_back(static_cast(nTrkTPCAll)); } void process(MyCollisions::iterator const& coll, aod::BCsWithTimestamps const&, aod::FT0s const&, aod::FV0As const&, MyTracks const& tracks) { - std::vector TrkTPCposLabel{}; - std::vector TrkTPCnegLabel{}; - std::vector TrkTPCallLabel{}; - std::vector qvecRe{}; - std::vector qvecIm{}; - std::vector qvecAmp{}; - - std::vector qvecReFT0C{}; - std::vector qvecImFT0C{}; - std::vector qvecReFT0A{}; - std::vector qvecImFT0A{}; - std::vector qvecReFT0M{}; - std::vector qvecImFT0M{}; - std::vector qvecReFV0A{}; - std::vector qvecImFV0A{}; - std::vector qvecReTPCpos{}; - std::vector qvecImTPCpos{}; - std::vector qvecReTPCneg{}; - std::vector qvecImTPCneg{}; - std::vector qvecReTPCall{}; - std::vector qvecImTPCall{}; + std::vector trkTPCPosLabel{}; + std::vector trkTPCNegLabel{}; + std::vector trkTPCAllLabel{}; + std::vector qVecAmp{}; + + std::vector qVecReSp{}; + std::vector qVecImSp{}; + std::vector qVecReFT0CSp{}; + std::vector qVecImFT0CSp{}; + std::vector qVecReFT0ASp{}; + std::vector qVecImFT0ASp{}; + std::vector qVecReFT0MSp{}; + std::vector qVecImFT0MSp{}; + std::vector qVecReFV0ASp{}; + std::vector qVecImFV0ASp{}; + std::vector qVecReTPCPosSp{}; + std::vector qVecImTPCPosSp{}; + std::vector qVecReTPCNegSp{}; + std::vector qVecImTPCNegSp{}; + std::vector qVecReTPCAllSp{}; + std::vector qVecImTPCAllSp{}; + + std::vector qVecReEse{}; + std::vector qVecImEse{}; + std::vector qVecReFT0CEse{}; + std::vector qVecImFT0CEse{}; + std::vector qVecReFT0AEse{}; + std::vector qVecImFT0AEse{}; + std::vector qVecReFT0MEse{}; + std::vector qVecImFT0MEse{}; + std::vector qVecReFV0AEse{}; + std::vector qVecImFV0AEse{}; + std::vector qVecReTPCPosEse{}; + std::vector qVecImTPCPosEse{}; + std::vector qVecReTPCNegEse{}; + std::vector qVecImTPCNegEse{}; + std::vector qVecReTPCAllEse{}; + std::vector qVecImTPCAllEse{}; auto bc = coll.bc_as(); int currentRun = bc.runNumber(); @@ -590,167 +764,138 @@ struct qVectorsTable { coll.centFT0M(), coll.centFT0A(), coll.centFT0C(), coll.centFV0A()}; cent = centAllEstim[cfgCentEsti]; - bool IsCalibrated = true; + bool isCalibrated = true; if (cent < 0. || cent > cfgMaxCentrality) { cent = 110.; - IsCalibrated = false; + isCalibrated = false; } - for (std::size_t id = 0; id < cfgnMods->size(); id++) { - int nmode = cfgnMods->at(id); - CalQvec(nmode, coll, tracks, qvecRe, qvecIm, qvecAmp, TrkTPCposLabel, TrkTPCnegLabel, TrkTPCallLabel); - if (cent < cfgMaxCentrality) { - for (auto i{0u}; i < kTPCall + 1; i++) { - helperEP.DoRecenter(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 1], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 1], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 1, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 2, i + 1)); - - helperEP.DoRecenter(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 2], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 2], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 1, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 2, i + 1)); - helperEP.DoTwist(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 2], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 2], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 3, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 4, i + 1)); - - helperEP.DoRecenter(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 3], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 3], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 1, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 2, i + 1)); - helperEP.DoTwist(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 3], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 3], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 3, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 4, i + 1)); - helperEP.DoRescale(qvecRe[(kTPCall + 1) * 4 * id + i * 4 + 3], qvecIm[(kTPCall + 1) * 4 * id + i * 4 + 3], - objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 5, i + 1), objQvec.at(id)->GetBinContent(static_cast(cent) + 1, 6, i + 1)); - } - if (cfgShiftCorr) { - auto deltapsiFT0C = 0.0; - auto deltapsiFT0A = 0.0; - auto deltapsiFT0M = 0.0; - auto deltapsiFV0A = 0.0; - auto deltapsiTPCpos = 0.0; - auto deltapsiTPCneg = 0.0; - auto deltapsiTPCall = 0.0; - - auto psidefFT0C = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3]) / static_cast(nmode); - auto psidefFT0A = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3]) / static_cast(nmode); - auto psidefFT0M = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3]) / static_cast(nmode); - auto psidefFV0A = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3]) / static_cast(nmode); - auto psidefTPCpos = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3]) / static_cast(nmode); - auto psidefTPCneg = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3]) / static_cast(nmode); - auto psidefTPCall = TMath::ATan2(qvecIm[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3], qvecRe[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3]) / static_cast(nmode); - - for (int ishift = 1; ishift <= 10; ishift++) { - auto coeffshiftxFT0C = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0C, ishift - 0.5)); - auto coeffshiftyFT0C = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0C + 1, ishift - 0.5)); - auto coeffshiftxFT0A = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0A, ishift - 0.5)); - auto coeffshiftyFT0A = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0A + 1, ishift - 0.5)); - auto coeffshiftxFT0M = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0M, ishift - 0.5)); - auto coeffshiftyFT0M = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFT0M + 1, ishift - 0.5)); - auto coeffshiftxFV0A = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFV0A, ishift - 0.5)); - auto coeffshiftyFV0A = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kFV0A + 1, ishift - 0.5)); - auto coeffshiftxTPCpos = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCpos, ishift - 0.5)); - auto coeffshiftyTPCpos = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCpos + 1, ishift - 0.5)); - auto coeffshiftxTPCneg = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCneg, ishift - 0.5)); - auto coeffshiftyTPCneg = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCneg + 1, ishift - 0.5)); - auto coeffshiftxTPCall = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCall, ishift - 0.5)); - auto coeffshiftyTPCall = shiftprofile.at(nmode - 2)->GetBinContent(shiftprofile.at(nmode - 2)->FindBin(cent, 2 * kTPCall + 1, ishift - 0.5)); - - deltapsiFT0C += ((2. / (1.0 * ishift)) * (-coeffshiftxFT0C * TMath::Cos(ishift * static_cast(nmode) * psidefFT0C) + coeffshiftyFT0C * TMath::Sin(ishift * static_cast(nmode) * psidefFT0C))) / static_cast(nmode); - deltapsiFT0A += ((2. / (1.0 * ishift)) * (-coeffshiftxFT0A * TMath::Cos(ishift * static_cast(nmode) * psidefFT0A) + coeffshiftyFT0A * TMath::Sin(ishift * static_cast(nmode) * psidefFT0A))) / static_cast(nmode); - deltapsiFT0M += ((2. / (1.0 * ishift)) * (-coeffshiftxFT0M * TMath::Cos(ishift * static_cast(nmode) * psidefFT0M) + coeffshiftyFT0M * TMath::Sin(ishift * static_cast(nmode) * psidefFT0M))) / static_cast(nmode); - deltapsiFV0A += ((2. / (1.0 * ishift)) * (-coeffshiftxFV0A * TMath::Cos(ishift * static_cast(nmode) * psidefFV0A) + coeffshiftyFV0A * TMath::Sin(ishift * static_cast(nmode) * psidefFV0A))) / static_cast(nmode); - deltapsiTPCpos += ((2. / (1.0 * ishift)) * (-coeffshiftxTPCpos * TMath::Cos(ishift * static_cast(nmode) * psidefTPCpos) + coeffshiftyTPCpos * TMath::Sin(ishift * static_cast(nmode) * psidefTPCpos))) / static_cast(nmode); - deltapsiTPCneg += ((2. / (1.0 * ishift)) * (-coeffshiftxTPCneg * TMath::Cos(ishift * static_cast(nmode) * psidefTPCneg) + coeffshiftyTPCneg * TMath::Sin(ishift * static_cast(nmode) * psidefTPCneg))) / static_cast(nmode); - deltapsiTPCall += ((2. / (1.0 * ishift)) * (-coeffshiftxTPCall * TMath::Cos(ishift * static_cast(nmode) * psidefTPCall) + coeffshiftyTPCall * TMath::Sin(ishift * static_cast(nmode) * psidefTPCall))) / static_cast(nmode); - } - deltapsiFT0C *= static_cast(nmode); - deltapsiFT0A *= static_cast(nmode); - deltapsiFT0M *= static_cast(nmode); - deltapsiFV0A *= static_cast(nmode); - deltapsiTPCpos *= static_cast(nmode); - deltapsiTPCneg *= static_cast(nmode); - deltapsiTPCall *= static_cast(nmode); - - float qvecReShiftedFT0C = qvecRe[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] * TMath::Cos(deltapsiFT0C) - qvecIm[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] * TMath::Sin(deltapsiFT0C); - float qvecImShiftedFT0C = qvecRe[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] * TMath::Sin(deltapsiFT0C) + qvecIm[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] * TMath::Cos(deltapsiFT0C); - float qvecReShiftedFT0A = qvecRe[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] * TMath::Cos(deltapsiFT0A) - qvecIm[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] * TMath::Sin(deltapsiFT0A); - float qvecImShiftedFT0A = qvecRe[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] * TMath::Sin(deltapsiFT0A) + qvecIm[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] * TMath::Cos(deltapsiFT0A); - float qvecReShiftedFT0M = qvecRe[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] * TMath::Cos(deltapsiFT0M) - qvecIm[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] * TMath::Sin(deltapsiFT0M); - float qvecImShiftedFT0M = qvecRe[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] * TMath::Sin(deltapsiFT0M) + qvecIm[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] * TMath::Cos(deltapsiFT0M); - float qvecReShiftedFV0A = qvecRe[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] * TMath::Cos(deltapsiFV0A) - qvecIm[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] * TMath::Sin(deltapsiFV0A); - float qvecImShiftedFV0A = qvecRe[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] * TMath::Sin(deltapsiFV0A) + qvecIm[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] * TMath::Cos(deltapsiFV0A); - float qvecReShiftedTPCpos = qvecRe[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] * TMath::Cos(deltapsiTPCpos) - qvecIm[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] * TMath::Sin(deltapsiTPCpos); - float qvecImShiftedTPCpos = qvecRe[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] * TMath::Sin(deltapsiTPCpos) + qvecIm[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] * TMath::Cos(deltapsiTPCpos); - float qvecReShiftedTPCneg = qvecRe[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] * TMath::Cos(deltapsiTPCneg) - qvecIm[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] * TMath::Sin(deltapsiTPCneg); - float qvecImShiftedTPCneg = qvecRe[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] * TMath::Sin(deltapsiTPCneg) + qvecIm[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] * TMath::Cos(deltapsiTPCneg); - float qvecReShiftedTPCall = qvecRe[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] * TMath::Cos(deltapsiTPCall) - qvecIm[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] * TMath::Sin(deltapsiTPCall); - float qvecImShiftedTPCall = qvecRe[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] * TMath::Sin(deltapsiTPCall) + qvecIm[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] * TMath::Cos(deltapsiTPCall); - - qvecRe[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] = qvecReShiftedFT0C; - qvecIm[(kTPCall + 1) * 4 * id + kFT0C * 4 + 3] = qvecImShiftedFT0C; - qvecRe[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] = qvecReShiftedFT0A; - qvecIm[(kTPCall + 1) * 4 * id + kFT0A * 4 + 3] = qvecImShiftedFT0A; - qvecRe[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] = qvecReShiftedFT0M; - qvecIm[(kTPCall + 1) * 4 * id + kFT0M * 4 + 3] = qvecImShiftedFT0M; - qvecRe[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] = qvecReShiftedFV0A; - qvecIm[(kTPCall + 1) * 4 * id + kFV0A * 4 + 3] = qvecImShiftedFV0A; - qvecRe[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] = qvecReShiftedTPCpos; - qvecIm[(kTPCall + 1) * 4 * id + kTPCpos * 4 + 3] = qvecImShiftedTPCpos; - qvecRe[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] = qvecReShiftedTPCneg; - qvecIm[(kTPCall + 1) * 4 * id + kTPCneg * 4 + 3] = qvecImShiftedTPCneg; - qvecRe[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] = qvecReShiftedTPCall; - qvecIm[(kTPCall + 1) * 4 * id + kTPCall * 4 + 3] = qvecImShiftedTPCall; - } + for (std::size_t id = 0; id < cfgnMods->size(); id++) { + int nMode = cfgnMods->at(id); + + // Raw Q-vectors, no multiplicity normalization and no corrections + std::vector qVecReRaw{}; + std::vector qVecImRaw{}; + calcQVec(nMode, coll, tracks, qVecReRaw, qVecImRaw, qVecAmp, trkTPCPosLabel, trkTPCNegLabel, trkTPCAllLabel); + + // Scalar Product Q-vectors, normalization by multiplicity/amplitude + std::vector nModeQVecReSp{}; + std::vector nModeQVecImSp{}; + normalizeQVec(nModeQVecReSp, nModeQVecImSp, qVecReRaw, qVecImRaw, qVecAmp, MultNorms::kScalarProd); + correctQVec(cent, nModeQVecReSp, nModeQVecImSp, corrsQvecSp[id], shiftProfileSp, nMode); + // Add to summary vector + qVecReSp.insert(qVecReSp.end(), nModeQVecReSp.begin(), nModeQVecReSp.end()); + qVecImSp.insert(qVecImSp.end(), nModeQVecImSp.begin(), nModeQVecImSp.end()); + + // Pick the desired correction level for the Q-vectors to be stored in the analysis table + int corrLevel = cfgCorrLevel == 0 ? 0 : cfgCorrLevel - 1; + int nCorrections = static_cast(kNCorrections); + + qVecReFT0CSp.push_back(nModeQVecReSp[kFT0C * nCorrections + corrLevel]); + qVecImFT0CSp.push_back(nModeQVecImSp[kFT0C * nCorrections + corrLevel]); + qVecReFT0ASp.push_back(nModeQVecReSp[kFT0A * nCorrections + corrLevel]); + qVecImFT0ASp.push_back(nModeQVecImSp[kFT0A * nCorrections + corrLevel]); + qVecReFT0MSp.push_back(nModeQVecReSp[kFT0M * nCorrections + corrLevel]); + qVecImFT0MSp.push_back(nModeQVecImSp[kFT0M * nCorrections + corrLevel]); + qVecReFV0ASp.push_back(nModeQVecReSp[kFV0A * nCorrections + corrLevel]); + qVecImFV0ASp.push_back(nModeQVecImSp[kFV0A * nCorrections + corrLevel]); + qVecReTPCPosSp.push_back(nModeQVecReSp[kTPCPos * nCorrections + corrLevel]); + qVecImTPCPosSp.push_back(nModeQVecImSp[kTPCPos * nCorrections + corrLevel]); + qVecReTPCNegSp.push_back(nModeQVecReSp[kTPCNeg * nCorrections + corrLevel]); + qVecImTPCNegSp.push_back(nModeQVecImSp[kTPCNeg * nCorrections + corrLevel]); + qVecReTPCAllSp.push_back(nModeQVecReSp[kTPCAll * nCorrections + corrLevel]); + qVecImTPCAllSp.push_back(nModeQVecImSp[kTPCAll * nCorrections + corrLevel]); + + if (cfgProduceRedQVecs) { + + // Ese Q-vectors, normalization by sqrt(multiplicity/amplitude) + std::vector nModeQVecReEse{}; + std::vector nModeQVecImEse{}; + normalizeQVec(nModeQVecReEse, nModeQVecImEse, qVecReRaw, qVecImRaw, qVecAmp, MultNorms::kEsE); + correctQVec(cent, nModeQVecReEse, nModeQVecImEse, corrsQvecEse[id], shiftProfileEse, nMode); + // Add to summary vector + qVecReEse.insert(qVecReEse.end(), nModeQVecReEse.begin(), nModeQVecReEse.end()); + qVecImEse.insert(qVecImEse.end(), nModeQVecImEse.begin(), nModeQVecImEse.end()); + + qVecReFT0CEse.push_back(nModeQVecReEse[kFT0C * nCorrections + corrLevel]); + qVecImFT0CEse.push_back(nModeQVecImEse[kFT0C * nCorrections + corrLevel]); + qVecReFT0AEse.push_back(nModeQVecReEse[kFT0A * nCorrections + corrLevel]); + qVecImFT0AEse.push_back(nModeQVecImEse[kFT0A * nCorrections + corrLevel]); + qVecReFT0MEse.push_back(nModeQVecReEse[kFT0M * nCorrections + corrLevel]); + qVecImFT0MEse.push_back(nModeQVecImEse[kFT0M * nCorrections + corrLevel]); + qVecReFV0AEse.push_back(nModeQVecReEse[kFV0A * nCorrections + corrLevel]); + qVecImFV0AEse.push_back(nModeQVecImEse[kFV0A * nCorrections + corrLevel]); + qVecReTPCPosEse.push_back(nModeQVecReEse[kTPCPos * nCorrections + corrLevel]); + qVecImTPCPosEse.push_back(nModeQVecImEse[kTPCPos * nCorrections + corrLevel]); + qVecReTPCNegEse.push_back(nModeQVecReEse[kTPCNeg * nCorrections + corrLevel]); + qVecImTPCNegEse.push_back(nModeQVecImEse[kTPCNeg * nCorrections + corrLevel]); + qVecReTPCAllEse.push_back(nModeQVecReEse[kTPCAll * nCorrections + corrLevel]); + qVecImTPCAllEse.push_back(nModeQVecImEse[kTPCAll * nCorrections + corrLevel]); } - int CorrLevel = cfgCorrLevel == 0 ? 0 : cfgCorrLevel - 1; - qvecReFT0C.push_back(qvecRe[(kTPCall + 1) * 4 * id + kFT0C * 4 + CorrLevel]); - qvecImFT0C.push_back(qvecIm[(kTPCall + 1) * 4 * id + kFT0C * 4 + CorrLevel]); - qvecReFT0A.push_back(qvecRe[(kTPCall + 1) * 4 * id + kFT0A * 4 + CorrLevel]); - qvecImFT0A.push_back(qvecIm[(kTPCall + 1) * 4 * id + kFT0A * 4 + CorrLevel]); - qvecReFT0M.push_back(qvecRe[(kTPCall + 1) * 4 * id + kFT0M * 4 + CorrLevel]); - qvecImFT0M.push_back(qvecIm[(kTPCall + 1) * 4 * id + kFT0M * 4 + CorrLevel]); - qvecReFV0A.push_back(qvecRe[(kTPCall + 1) * 4 * id + kFV0A * 4 + CorrLevel]); - qvecImFV0A.push_back(qvecIm[(kTPCall + 1) * 4 * id + kFV0A * 4 + CorrLevel]); - qvecReTPCpos.push_back(qvecRe[(kTPCall + 1) * 4 * id + kTPCpos * 4 + CorrLevel]); - qvecImTPCpos.push_back(qvecIm[(kTPCall + 1) * 4 * id + kTPCpos * 4 + CorrLevel]); - qvecReTPCneg.push_back(qvecRe[(kTPCall + 1) * 4 * id + kTPCneg * 4 + CorrLevel]); - qvecImTPCneg.push_back(qvecIm[(kTPCall + 1) * 4 * id + kTPCneg * 4 + CorrLevel]); - qvecReTPCall.push_back(qvecRe[(kTPCall + 1) * 4 * id + kTPCall * 4 + CorrLevel]); - qvecImTPCall.push_back(qvecIm[(kTPCall + 1) * 4 * id + kTPCall * 4 + CorrLevel]); } // Fill the columns of the Qvectors table. - qVector(cent, IsCalibrated, qvecRe, qvecIm, qvecAmp); + qVector(cent, isCalibrated, qVecReSp, qVecImSp, qVecAmp); if (useDetector["QvectorFT0Cs"]) - qVectorFT0C(IsCalibrated, qvecReFT0C.at(0), qvecImFT0C.at(0), qvecAmp[kFT0C]); + qVectorFT0C(isCalibrated, qVecReFT0CSp.at(0), qVecImFT0CSp.at(0), qVecAmp[kFT0C]); if (useDetector["QvectorFT0As"]) - qVectorFT0A(IsCalibrated, qvecReFT0A.at(0), qvecImFT0A.at(0), qvecAmp[kFT0A]); + qVectorFT0A(isCalibrated, qVecReFT0ASp.at(0), qVecImFT0ASp.at(0), qVecAmp[kFT0A]); if (useDetector["QvectorFT0Ms"]) - qVectorFT0M(IsCalibrated, qvecReFT0M.at(0), qvecImFT0M.at(0), qvecAmp[kFT0M]); + qVectorFT0M(isCalibrated, qVecReFT0MSp.at(0), qVecImFT0MSp.at(0), qVecAmp[kFT0M]); if (useDetector["QvectorFV0As"]) - qVectorFV0A(IsCalibrated, qvecReFV0A.at(0), qvecImFV0A.at(0), qvecAmp[kFV0A]); + qVectorFV0A(isCalibrated, qVecReFV0ASp.at(0), qVecImFV0ASp.at(0), qVecAmp[kFV0A]); if (useDetector["QvectorTPCposs"]) - qVectorTPCpos(IsCalibrated, qvecReTPCpos.at(0), qvecImTPCpos.at(0), qvecAmp[kTPCpos], TrkTPCposLabel); + qVectorTPCPos(isCalibrated, qVecReTPCPosSp.at(0), qVecImTPCPosSp.at(0), qVecAmp[kTPCPos], trkTPCPosLabel); if (useDetector["QvectorTPCnegs"]) - qVectorTPCneg(IsCalibrated, qvecReTPCneg.at(0), qvecImTPCneg.at(0), qvecAmp[kTPCneg], TrkTPCnegLabel); + qVectorTPCNeg(isCalibrated, qVecReTPCNegSp.at(0), qVecImTPCNegSp.at(0), qVecAmp[kTPCNeg], trkTPCNegLabel); if (useDetector["QvectorTPCalls"]) - qVectorTPCall(IsCalibrated, qvecReTPCall.at(0), qvecImTPCall.at(0), qvecAmp[kTPCall], TrkTPCallLabel); + qVectorTPCAll(isCalibrated, qVecReTPCAllSp.at(0), qVecImTPCAllSp.at(0), qVecAmp[kTPCAll], trkTPCAllLabel); - qVectorFT0CVec(IsCalibrated, qvecReFT0C, qvecImFT0C, qvecAmp[kFT0C]); - qVectorFT0AVec(IsCalibrated, qvecReFT0A, qvecImFT0A, qvecAmp[kFT0A]); - qVectorFT0MVec(IsCalibrated, qvecReFT0M, qvecImFT0M, qvecAmp[kFT0M]); - qVectorFV0AVec(IsCalibrated, qvecReFV0A, qvecImFV0A, qvecAmp[kFV0A]); - qVectorTPCposVec(IsCalibrated, qvecReTPCpos, qvecImTPCpos, qvecAmp[kTPCpos], TrkTPCposLabel); - qVectorTPCnegVec(IsCalibrated, qvecReTPCneg, qvecImTPCneg, qvecAmp[kTPCneg], TrkTPCnegLabel); - qVectorTPCallVec(IsCalibrated, qvecReTPCall, qvecImTPCall, qvecAmp[kTPCall], TrkTPCallLabel); + qVectorFT0CVec(isCalibrated, qVecReFT0CSp, qVecImFT0CSp, qVecAmp[kFT0C]); + qVectorFT0AVec(isCalibrated, qVecReFT0ASp, qVecImFT0ASp, qVecAmp[kFT0A]); + qVectorFT0MVec(isCalibrated, qVecReFT0MSp, qVecImFT0MSp, qVecAmp[kFT0M]); + qVectorFV0AVec(isCalibrated, qVecReFV0ASp, qVecImFV0ASp, qVecAmp[kFV0A]); + qVectorTPCPosVec(isCalibrated, qVecReTPCPosSp, qVecImTPCPosSp, qVecAmp[kTPCPos], trkTPCPosLabel); + qVectorTPCNegVec(isCalibrated, qVecReTPCNegSp, qVecImTPCNegSp, qVecAmp[kTPCNeg], trkTPCNegLabel); + qVectorTPCAllVec(isCalibrated, qVecReTPCAllSp, qVecImTPCAllSp, qVecAmp[kTPCAll], trkTPCAllLabel); // Deprecated, will be removed in future after transition time // if (useDetector["QvectorBPoss"]) - qVectorBPos(IsCalibrated, qvecReTPCpos.at(0), qvecImTPCpos.at(0), qvecAmp[kTPCpos], TrkTPCposLabel); + qVectorBPos(isCalibrated, qVecReTPCPosSp.at(0), qVecImTPCPosSp.at(0), qVecAmp[kTPCPos], trkTPCPosLabel); if (useDetector["QvectorBNegs"]) - qVectorBNeg(IsCalibrated, qvecReTPCneg.at(0), qvecImTPCneg.at(0), qvecAmp[kTPCneg], TrkTPCnegLabel); + qVectorBNeg(isCalibrated, qVecReTPCNegSp.at(0), qVecImTPCNegSp.at(0), qVecAmp[kTPCNeg], trkTPCNegLabel); if (useDetector["QvectorBTots"]) - qVectorBTot(IsCalibrated, qvecReTPCall.at(0), qvecImTPCall.at(0), qvecAmp[kTPCall], TrkTPCallLabel); + qVectorBTot(isCalibrated, qVecReTPCAllSp.at(0), qVecImTPCAllSp.at(0), qVecAmp[kTPCAll], trkTPCAllLabel); - qVectorBPosVec(IsCalibrated, qvecReTPCpos, qvecImTPCpos, qvecAmp[kTPCpos], TrkTPCposLabel); - qVectorBNegVec(IsCalibrated, qvecReTPCneg, qvecImTPCneg, qvecAmp[kTPCneg], TrkTPCnegLabel); - qVectorBTotVec(IsCalibrated, qvecReTPCall, qvecImTPCall, qvecAmp[kTPCall], TrkTPCallLabel); + qVectorBPosVec(isCalibrated, qVecReTPCPosSp, qVecImTPCPosSp, qVecAmp[kTPCPos], trkTPCPosLabel); + qVectorBNegVec(isCalibrated, qVecReTPCNegSp, qVecImTPCNegSp, qVecAmp[kTPCNeg], trkTPCNegLabel); + qVectorBTotVec(isCalibrated, qVecReTPCAllSp, qVecImTPCAllSp, qVecAmp[kTPCAll], trkTPCAllLabel); ///////////////////////////////////////////////////////////////// + if (cfgProduceRedQVecs) { + eseQVector(cent, isCalibrated, qVecReEse, qVecImEse, qVecAmp); + eseQVectorFT0C(isCalibrated, qVecReFT0CEse.at(0), qVecImFT0CEse.at(0), qVecAmp[kFT0C]); + eseQVectorFT0A(isCalibrated, qVecReFT0AEse.at(0), qVecImFT0AEse.at(0), qVecAmp[kFT0A]); + eseQVectorFT0M(isCalibrated, qVecReFT0MEse.at(0), qVecImFT0MEse.at(0), qVecAmp[kFT0M]); + eseQVectorFV0A(isCalibrated, qVecReFV0AEse.at(0), qVecImFV0AEse.at(0), qVecAmp[kFV0A]); + eseQVectorTPCPos(isCalibrated, qVecReTPCPosEse.at(0), qVecImTPCPosEse.at(0), qVecAmp[kTPCPos], trkTPCPosLabel); + eseQVectorTPCNeg(isCalibrated, qVecReTPCNegEse.at(0), qVecImTPCNegEse.at(0), qVecAmp[kTPCNeg], trkTPCNegLabel); + eseQVectorTPCAll(isCalibrated, qVecReTPCAllEse.at(0), qVecImTPCAllEse.at(0), qVecAmp[kTPCAll], trkTPCAllLabel); + eseQVectorFT0CVec(isCalibrated, qVecReFT0CEse, qVecImFT0CEse, qVecAmp[kFT0C]); + eseQVectorFT0AVec(isCalibrated, qVecReFT0AEse, qVecImFT0AEse, qVecAmp[kFT0A]); + eseQVectorFT0MVec(isCalibrated, qVecReFT0MEse, qVecImFT0MEse, qVecAmp[kFT0M]); + eseQVectorFV0AVec(isCalibrated, qVecReFV0AEse, qVecImFV0AEse, qVecAmp[kFV0A]); + eseQVectorTPCPosVec(isCalibrated, qVecReTPCPosEse, qVecImTPCPosEse, qVecAmp[kTPCPos], trkTPCPosLabel); + eseQVectorTPCNegVec(isCalibrated, qVecReTPCNegEse, qVecImTPCNegEse, qVecAmp[kTPCNeg], trkTPCNegLabel); + eseQVectorTPCAllVec(isCalibrated, qVecReTPCAllEse, qVecImTPCAllEse, qVecAmp[kTPCAll], trkTPCAllLabel); + eseQVectorPerc(qVecReFT0CEse.at(0), qVecImFT0CEse.at(0), qVecAmp[kFT0C], + qVecReFT0AEse.at(0), qVecImFT0AEse.at(0), qVecAmp[kFT0A], + qVecReFT0MEse.at(0), qVecImFT0MEse.at(0), qVecAmp[kFT0M], + qVecReFV0AEse.at(0), qVecImFV0AEse.at(0), qVecAmp[kFV0A], + qVecReTPCPosEse.at(0), qVecImTPCPosEse.at(0), qVecAmp[kTPCPos], + qVecReTPCNegEse.at(0), qVecImTPCNegEse.at(0), qVecAmp[kTPCNeg], + qVecReTPCAllEse.at(0), qVecImTPCAllEse.at(0), qVecAmp[kTPCAll]); + } } // End process. }; diff --git a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx index 78462e10a73..14d9c2c190f 100644 --- a/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx +++ b/PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx @@ -75,6 +75,7 @@ DECLARE_SOA_COLUMN(MlScore0, mlScore0, float); //! ML score of the first con DECLARE_SOA_COLUMN(MlScore1, mlScore1, float); //! ML score of the second configured index DECLARE_SOA_COLUMN(ScalarProd, scalarProd, float); //! Scalar product DECLARE_SOA_COLUMN(Cent, cent, float); //! Centrality +DECLARE_SOA_COLUMN(RedQVec, redQVec, float); //! Reduced Q-vector } // namespace full DECLARE_SOA_TABLE(HfCandMPtInfos, "AOD", "HFCANDMPTINFO", full::M, @@ -89,6 +90,7 @@ DECLARE_SOA_TABLE(HfCandFlowInfos, "AOD", "HFCANDFLOWINFO", full::MlScore1, full::ScalarProd, full::Cent); +DECLARE_SOA_TABLE(HfRedQVecEsEs, "AOD", "HFREDQVECESE", full::RedQVec); } // namespace o2::aod enum DecayChannel { DplusToPiKPi = 0, @@ -106,15 +108,18 @@ enum DecayChannel { DplusToPiKPi = 0, struct HfTaskFlowCharmHadrons { Produces rowCandMassPtMl; Produces rowCandMassPtMlSpCent; + Produces rowRedQVecEsE; Configurable harmonic{"harmonic", 2, "harmonic number"}; Configurable qVecDetector{"qVecDetector", 3, "Detector for Q vector estimation (FV0A: 0, FT0M: 1, FT0A: 2, FT0C: 3, TPC Pos: 4, TPC Neg: 5, TPC Tot: 6)"}; + Configurable qVecRedDetector{"qVecRedDetector", 6, "Detector for Q vector estimation (FT0C: 3, TPC Pos: 4, TPC Neg: 5, TPC Tot: 6)"}; Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3, FV0A: 4, NTracksPV: 5, FT0CVariant2: 6)"}; Configurable selectionFlag{"selectionFlag", 1, "Selection Flag for hadron (e.g. 1 for skimming, 3 for topo. and kine., 7 for PID)"}; Configurable centralityMin{"centralityMin", 0., "Minimum centrality accepted in SP/EP computation (not applied in resolution process)"}; Configurable centralityMax{"centralityMax", 100., "Maximum centrality accepted in SP/EP computation (not applied in resolution process)"}; Configurable storeEP{"storeEP", false, "Flag to store EP-related axis"}; Configurable storeMl{"storeMl", false, "Flag to store ML scores"}; + Configurable storeRedQVec{"storeRedQVec", false, "Flag to store reduced Q-vectors for ESE"}; Configurable fillMassPtMlTree{"fillMassPtMlTree", false, "Flag to fill mass, pt and ML scores tree"}; Configurable fillMassPtMlSpCentTree{"fillMassPtMlSpCentTree", false, "Flag to fill mass, pt, ML scores, SP and centrality tree"}; Configurable fillSparse{"fillSparse", true, "Flag to fill sparse"}; @@ -148,7 +153,7 @@ struct HfTaskFlowCharmHadrons { using CandXic0DataWMl = soa::Filtered>; using CandD0DataWMl = soa::Filtered>; using CandD0Data = soa::Filtered>; - using CollsWithQvecs = soa::Join; + using CollsWithQvecs = soa::Join; using TracksWithExtra = soa::Join; Filter filterSelectDsCandidates = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag || aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag; @@ -198,6 +203,7 @@ struct HfTaskFlowCharmHadrons { ConfigurableAxis thnConfigAxisResoFV0aTPCtot{"thnConfigAxisResoFV0aTPCtot", {160, -8, 8}, ""}; ConfigurableAxis thnConfigAxisCandidateEta{"thnConfigAxisCandidateEta", {100, -5, 5}, ""}; ConfigurableAxis thnConfigAxisSign{"thnConfigAxisSign", {6, -3.0, 3.0}, ""}; + ConfigurableAxis thnConfigAxisRedQVec{"thnConfigAxisRedQVec", {1000, 0, 100}, ""}; HistogramRegistry registry{"registry", {}}; @@ -225,6 +231,7 @@ struct HfTaskFlowCharmHadrons { const AxisSpec thnAxisNoCollInTimeRangeNarrow{thnConfigAxisNoCollInTimeRangeNarrow, "NoCollInTimeRangeNarrow"}; const AxisSpec thnAxisNoCollInTimeRangeStandard{thnConfigAxisNoCollInTimeRangeStandard, "NoCollInTimeRangeStandard"}; const AxisSpec thnAxisNoCollInRofStandard{thnConfigAxisNoCollInRofStandard, "NoCollInRofStandard"}; + const AxisSpec thnAxisRedQVec{thnConfigAxisRedQVec, "Reduced Q-vector"}; // TODO: currently only the Q vector of FT0c FV0a and TPCtot are considered const AxisSpec thnAxisResoFT0cFV0a{thnConfigAxisResoFT0cFV0a, "Q_{FT0c} #bullet Q_{FV0a}"}; const AxisSpec thnAxisResoFT0cTPCtot{thnConfigAxisResoFT0cTPCtot, "Q_{FT0c} #bullet Q_{TPCtot}"}; @@ -252,6 +259,9 @@ struct HfTaskFlowCharmHadrons { thnAxisNoCollInTimeRangeNarrow, thnAxisNoCollInTimeRangeStandard, thnAxisNoCollInRofStandard}); } } + if (storeRedQVec) { + axes.insert(axes.end(), {thnAxisRedQVec}); + } registry.add("hSparseFlowCharm", "THn for SP", HistType::kTHnSparseF, axes); registry.add("hCentEventWithCand", "Centrality distributions with charm candidates;Cent;entries", HistType::kTH1F, {{100, 0.f, 100.f}}); registry.add("hCentEventWithCandInSigRegion", "Centrality distributions with charm candidates in signal range;Cent;entries", HistType::kTH1F, {{100, 0.f, 100.f}}); @@ -288,6 +298,14 @@ struct HfTaskFlowCharmHadrons { registry.add("spReso/hSpResoFV0aTPCtot", "hSpResoFV0aTPCtot; centrality; Q_{FV0a} #bullet Q_{TPCtot}", {HistType::kTH2F, {thnAxisCent, thnAxisScalarProd}}); registry.add("spReso/hSpResoTPCposTPCneg", "hSpResoTPCposTPCneg; centrality; Q_{TPCpos} #bullet Q_{TPCneg}", {HistType::kTH2F, {thnAxisCent, thnAxisScalarProd}}); + if (storeRedQVec) { + registry.add("redQVecs/hSparseRedQVecs", "hSpResoRedQVec; centrality; Q_{red} #bullet Q_{TPCtot}", {HistType::kTHnSparseF, {thnAxisCent, thnAxisRedQVec, thnAxisRedQVec, thnAxisRedQVec, thnAxisRedQVec}}); + registry.add("redQVecs/hRedQVecFT0C", "hRedQVecFT0C; centrality; q_{FT0c}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}}); + registry.add("redQVecs/hRedQVecTpcPos", "hRedQVecTpcPos; centrality; q_{TPCpos}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}}); + registry.add("redQVecs/hRedQVecTpcNeg", "hRedQVecTpcNeg; centrality; q_{TPCneg}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}}); + registry.add("redQVecs/hRedQVecTpcAll", "hRedQVecTpcAll; centrality; q_{TPCall}", {HistType::kTH2F, {thnAxisCent, thnAxisRedQVec}}); + } + if (saveEpResoHisto) { registry.add("epReso/hEpResoFT0cFT0a", "hEpResoFT0cFT0a; centrality; #Delta#Psi_{sub}", {HistType::kTH2F, {thnAxisCent, thnAxisCosNPhi}}); registry.add("epReso/hEpResoFT0cFV0a", "hEpResoFT0cFV0a; centrality; #Delta#Psi_{sub}", {HistType::kTH2F, {thnAxisCent, thnAxisCosNPhi}}); @@ -430,6 +448,7 @@ struct HfTaskFlowCharmHadrons { /// \param outputMl are the ML scores /// \param occupancy is the occupancy of the collision using the track estimator /// \param hfevselflag flag of the collision associated to utilsEvSelHf.h + /// \param redQVec is the reduced Q vector for EsE void fillThn(const float mass, const float pt, const float eta, @@ -441,7 +460,8 @@ struct HfTaskFlowCharmHadrons { const float sp, const std::vector& outputMl, const float occupancy, - const o2::hf_evsel::HfCollisionRejectionMask hfevselflag) + const o2::hf_evsel::HfCollisionRejectionMask hfevselflag, + const float redQVec) { auto hSparse = registry.get(HIST("hSparseFlowCharm")); const int ndim = hSparse->GetNdimensions(); @@ -480,6 +500,9 @@ struct HfTaskFlowCharmHadrons { values.push_back(evtSelFlags[3]); values.push_back(evtSelFlags[4]); } + if (storeRedQVec) { + values.push_back(redQVec); + } if (static_cast(values.size()) != ndim) { LOGF(fatal, @@ -487,7 +510,6 @@ struct HfTaskFlowCharmHadrons { "does not match THnSparse dimensionality (%d).", static_cast(values.size()), ndim); } - hSparse->Fill(values.data()); } @@ -501,13 +523,16 @@ struct HfTaskFlowCharmHadrons { aod::BCsWithTimestamps const&, float& centrality) { - const auto occupancy = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); + float occupancy{-999.f}; + if (occEstimator != 0) { + occupancy = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); + registry.fill(HIST("trackOccVsFT0COcc"), collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange()); + } const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); centrality = o2::hf_centrality::getCentralityColl(collision, CentEstimator); /// monitor the satisfied event selections hfEvSel.fillHistograms(collision, rejectionMask, centrality, occupancy); - registry.fill(HIST("trackOccVsFT0COcc"), collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange()); return rejectionMask == 0; } @@ -543,6 +568,16 @@ struct HfTaskFlowCharmHadrons { float yQVec = qVecs[1]; float const amplQVec = qVecs[2]; float const evtPl = epHelper.GetEventPlane(xQVec, yQVec, harmonic); + + float redQVec{-999.f}; + std::array qVecRedComps{-999.f, -999.f, -999.f}; + if (storeRedQVec) { + qVecRedComps = getEseQvec(collision, qVecRedDetector.value); + } + float xRedQVec = qVecRedComps[0]; + float yRedQVec = qVecRedComps[1]; + float const amplRedQVec = qVecRedComps[2]; + for (const auto& candidate : candidates) { float massCand = 0.; float signCand = 0.; @@ -681,10 +716,16 @@ struct HfTaskFlowCharmHadrons { etaCand = candidate.eta(); } + float const cosNPhi = std::cos(harmonic * phiCand); + float const sinNPhi = std::sin(harmonic * phiCand); + float const cosDeltaPhi = std::cos(harmonic * (phiCand - evtPl)); + // If TPC is used for the SP estimation, the tracks of the hadron candidate must be removed from the corresponding TPC Q vector to avoid self-correlations - if (qVecDetector == QvecEstimator::TPCNeg || - qVecDetector == QvecEstimator::TPCPos || - qVecDetector == QvecEstimator::TPCTot) { + bool subtractDaugsFromQVec = (qVecDetector == QvecEstimator::TPCNeg || + qVecDetector == QvecEstimator::TPCPos || + qVecDetector == QvecEstimator::TPCTot); + float scalprodCand{-999.f}; + if (subtractDaugsFromQVec) { std::vector tracksQx; std::vector tracksQy; @@ -697,17 +738,13 @@ struct HfTaskFlowCharmHadrons { } // subtract daughters' contribution from the (normalized) Q-vector - for (std::size_t iTrack = 0; iTrack < tracksQx.size(); ++iTrack) { - xQVec -= tracksQx[iTrack]; - yQVec -= tracksQy[iTrack]; - } + float xQVecDaugSubtr = xQVec - std::accumulate(tracksQx.begin(), tracksQx.end(), 0.0); + float yQVecDaugSubtr = yQVec - std::accumulate(tracksQy.begin(), tracksQy.end(), 0.0); + scalprodCand = cosNPhi * xQVecDaugSubtr + sinNPhi * yQVecDaugSubtr; + } else { + scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec; } - float const cosNPhi = std::cos(harmonic * phiCand); - float const sinNPhi = std::sin(harmonic * phiCand); - float const scalprodCand = cosNPhi * xQVec + sinNPhi * yQVec; - float const cosDeltaPhi = std::cos(harmonic * (phiCand - evtPl)); - if (fillMassPtMlTree || fillMassPtMlSpCentTree) { if (downSampleFactor < 1.) { float const pseudoRndm = ptCand * 1000. - static_cast(ptCand * 1000); @@ -722,8 +759,37 @@ struct HfTaskFlowCharmHadrons { rowCandMassPtMlSpCent(massCand, ptCand, outputMl[0], outputMl[1], scalprodCand, cent); } } + + bool subtractDaugsFromRedQVec = storeRedQVec && (qVecRedDetector == QvecEstimator::TPCNeg || + qVecRedDetector == QvecEstimator::TPCPos || + qVecRedDetector == QvecEstimator::TPCTot); + if (subtractDaugsFromRedQVec) { + std::vector tracksRedQx; + std::vector tracksRedQy; + + // IMPORTANT: use the amplitude of the reduced Q-vector to build this Q-vector + if constexpr (std::is_same_v || std::is_same_v) { + getQvecXic0Tracks(candidate, tracksRedQx, tracksRedQy, std::sqrt(amplRedQVec), static_cast(qVecRedDetector.value)); + } else { + getQvecDtracks(candidate, tracksRedQx, tracksRedQy, std::sqrt(amplRedQVec), static_cast(qVecRedDetector.value)); + } + + // subtract daughters' contribution from the (normalized) Q-vector + const float redQVecXDaugSubtr = xRedQVec - std::accumulate(tracksRedQx.begin(), tracksRedQx.end(), 0.0); + const float redQVecYDaugSubtr = yRedQVec - std::accumulate(tracksRedQy.begin(), tracksRedQy.end(), 0.0); + if (qVecRedDetector.value == QvecEstimator::TPCTot || qVecRedDetector.value == QvecEstimator::TPCPos || qVecRedDetector.value == QvecEstimator::TPCNeg) { + // Correct for track multiplicity + redQVec = std::hypot(redQVecXDaugSubtr, redQVecYDaugSubtr) * std::sqrt(amplRedQVec) / std::sqrt(amplRedQVec - tracksRedQx.size()); + } else { + redQVec = std::hypot(xRedQVec, yRedQVec); + } + } + if (storeRedQVec) { + rowRedQVecEsE(redQVec); + } if (fillSparse) { - fillThn(massCand, ptCand, etaCand, signCand, cent, cosNPhi, sinNPhi, cosDeltaPhi, scalprodCand, outputMl, occupancy, hfevflag); + fillThn(massCand, ptCand, etaCand, signCand, cent, cosNPhi, sinNPhi, + cosDeltaPhi, scalprodCand, outputMl, occupancy, hfevflag, redQVec); } } if (hasCandInMassWin) { @@ -878,22 +944,25 @@ struct HfTaskFlowCharmHadrons { float const yQVecFT0m = collision.qvecFT0MIm(); float const xQVecFV0a = collision.qvecFV0ARe(); float const yQVecFV0a = collision.qvecFV0AIm(); - float const xQVecBPos = collision.qvecBPosRe(); - float const yQVecBPos = collision.qvecBPosIm(); - float const xQVecBNeg = collision.qvecBNegRe(); - float const yQVecBNeg = collision.qvecBNegIm(); - float const xQVecBTot = collision.qvecBTotRe(); - float const yQVecBTot = collision.qvecBTotIm(); + float const xQVecTPCPos = collision.qvecTPCposRe(); + float const yQVecTPCPos = collision.qvecTPCposIm(); + float const xQVecTPCNeg = collision.qvecTPCnegRe(); + float const yQVecTPCNeg = collision.qvecTPCnegIm(); + float const xQVecTPCAll = collision.qvecTPCallRe(); + float const yQVecTPCAll = collision.qvecTPCallIm(); centrality = o2::hf_centrality::getCentralityColl(collision, centEstimator); if (storeResoOccu) { - const auto occupancy = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); - registry.fill(HIST("trackOccVsFT0COcc"), collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange()); + float occupancy{-999.f}; + if (occEstimator != 0) { + occupancy = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); + registry.fill(HIST("trackOccVsFT0COcc"), collision.trackOccupancyInTimeRange(), collision.ft0cOccupancyInTimeRange()); + } const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); std::vector evtSelFlags = getEventSelectionFlags(rejectionMask); registry.fill(HIST("spReso/hSparseReso"), centrality, xQVecFT0c * xQVecFV0a + yQVecFT0c * yQVecFV0a, - xQVecFT0c * xQVecBTot + yQVecFT0c * yQVecBTot, - xQVecFV0a * xQVecBTot + yQVecFV0a * yQVecBTot, + xQVecFT0c * xQVecTPCAll + yQVecFT0c * yQVecTPCAll, + xQVecFV0a * xQVecTPCAll + yQVecFV0a * yQVecTPCAll, occupancy, evtSelFlags[0], evtSelFlags[1], evtSelFlags[2], evtSelFlags[3], evtSelFlags[4]); } @@ -901,6 +970,13 @@ struct HfTaskFlowCharmHadrons { registry.fill(HIST("hSparseCentEstimators"), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(0)), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(1)), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(2)), o2::hf_centrality::getCentralityColl(collision, centEstimatorsForSparse->at(3))); } + if (storeRedQVec) { + registry.fill(HIST("redQVecs/hSparseRedQVecs"), centrality, + std::hypot(collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm()), + std::hypot(collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm()), + std::hypot(collision.eseQvecTPCnegRe(), collision.eseQvecTPCnegIm()), + std::hypot(collision.eseQvecTPCallRe(), collision.eseQvecTPCallIm())); + } if (!isCollSelected(collision, bcs, centrality)) { // no selection on the centrality is applied, but on event selection flags @@ -909,30 +985,35 @@ struct HfTaskFlowCharmHadrons { registry.fill(HIST("spReso/hSpResoFT0cFT0a"), centrality, xQVecFT0c * xQVecFT0a + yQVecFT0c * yQVecFT0a); registry.fill(HIST("spReso/hSpResoFT0cFV0a"), centrality, xQVecFT0c * xQVecFV0a + yQVecFT0c * yQVecFV0a); - registry.fill(HIST("spReso/hSpResoFT0cTPCpos"), centrality, xQVecFT0c * xQVecBPos + yQVecFT0c * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFT0cTPCneg"), centrality, xQVecFT0c * xQVecBNeg + yQVecFT0c * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFT0cTPCtot"), centrality, xQVecFT0c * xQVecBTot + yQVecFT0c * yQVecBTot); + registry.fill(HIST("spReso/hSpResoFT0cTPCpos"), centrality, xQVecFT0c * xQVecTPCPos + yQVecFT0c * yQVecTPCPos); + registry.fill(HIST("spReso/hSpResoFT0cTPCneg"), centrality, xQVecFT0c * xQVecTPCNeg + yQVecFT0c * yQVecTPCNeg); + registry.fill(HIST("spReso/hSpResoFT0cTPCtot"), centrality, xQVecFT0c * xQVecTPCAll + yQVecFT0c * yQVecTPCAll); registry.fill(HIST("spReso/hSpResoFT0aFV0a"), centrality, xQVecFT0a * xQVecFV0a + yQVecFT0a * yQVecFV0a); - registry.fill(HIST("spReso/hSpResoFT0aTPCpos"), centrality, xQVecFT0a * xQVecBPos + yQVecFT0a * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFT0aTPCneg"), centrality, xQVecFT0a * xQVecBNeg + yQVecFT0a * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFT0aTPCtot"), centrality, xQVecFT0a * xQVecBTot + yQVecFT0a * yQVecBTot); + registry.fill(HIST("spReso/hSpResoFT0aTPCpos"), centrality, xQVecFT0a * xQVecTPCPos + yQVecFT0a * yQVecTPCPos); + registry.fill(HIST("spReso/hSpResoFT0aTPCneg"), centrality, xQVecFT0a * xQVecTPCNeg + yQVecFT0a * yQVecTPCNeg); + registry.fill(HIST("spReso/hSpResoFT0aTPCtot"), centrality, xQVecFT0a * xQVecTPCAll + yQVecFT0a * yQVecTPCAll); registry.fill(HIST("spReso/hSpResoFT0mFV0a"), centrality, xQVecFT0m * xQVecFV0a + yQVecFT0m * yQVecFV0a); - registry.fill(HIST("spReso/hSpResoFT0mTPCpos"), centrality, xQVecFT0m * xQVecBPos + yQVecFT0m * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFT0mTPCneg"), centrality, xQVecFT0m * xQVecBNeg + yQVecFT0m * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFT0mTPCtot"), centrality, xQVecFT0m * xQVecBTot + yQVecFT0m * yQVecBTot); - registry.fill(HIST("spReso/hSpResoFV0aTPCpos"), centrality, xQVecFV0a * xQVecBPos + yQVecFV0a * yQVecBPos); - registry.fill(HIST("spReso/hSpResoFV0aTPCneg"), centrality, xQVecFV0a * xQVecBNeg + yQVecFV0a * yQVecBNeg); - registry.fill(HIST("spReso/hSpResoFV0aTPCtot"), centrality, xQVecFV0a * xQVecBTot + yQVecFV0a * yQVecBTot); - registry.fill(HIST("spReso/hSpResoTPCposTPCneg"), centrality, xQVecBPos * xQVecBNeg + yQVecBPos * yQVecBNeg); + registry.fill(HIST("spReso/hSpResoFT0mTPCpos"), centrality, xQVecFT0m * xQVecTPCPos + yQVecFT0m * yQVecTPCPos); + registry.fill(HIST("spReso/hSpResoFT0mTPCneg"), centrality, xQVecFT0m * xQVecTPCNeg + yQVecFT0m * yQVecTPCNeg); + registry.fill(HIST("spReso/hSpResoFT0mTPCtot"), centrality, xQVecFT0m * xQVecTPCAll + yQVecFT0m * yQVecTPCAll); + registry.fill(HIST("spReso/hSpResoFV0aTPCpos"), centrality, xQVecFV0a * xQVecTPCPos + yQVecFV0a * yQVecTPCPos); + registry.fill(HIST("spReso/hSpResoFV0aTPCneg"), centrality, xQVecFV0a * xQVecTPCNeg + yQVecFV0a * yQVecTPCNeg); + registry.fill(HIST("spReso/hSpResoFV0aTPCtot"), centrality, xQVecFV0a * xQVecTPCAll + yQVecFV0a * yQVecTPCAll); + registry.fill(HIST("spReso/hSpResoTPCposTPCneg"), centrality, xQVecTPCPos * xQVecTPCNeg + yQVecTPCPos * yQVecTPCNeg); + + registry.fill(HIST("redQVecs/hRedQVecFT0C"), centrality, std::hypot(collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm())); + registry.fill(HIST("redQVecs/hRedQVecTpcPos"), centrality, std::hypot(collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm())); + registry.fill(HIST("redQVecs/hRedQVecTpcNeg"), centrality, std::hypot(collision.eseQvecTPCnegRe(), collision.eseQvecTPCnegIm())); + registry.fill(HIST("redQVecs/hRedQVecTpcAll"), centrality, std::hypot(collision.eseQvecTPCallRe(), collision.eseQvecTPCallIm())); if (saveEpResoHisto) { float const epFT0a = epHelper.GetEventPlane(xQVecFT0a, yQVecFT0a, harmonic); float const epFT0c = epHelper.GetEventPlane(xQVecFT0c, yQVecFT0c, harmonic); float const epFT0m = epHelper.GetEventPlane(xQVecFT0m, yQVecFT0m, harmonic); float const epFV0a = epHelper.GetEventPlane(xQVecFV0a, yQVecFV0a, harmonic); - float const epBPoss = epHelper.GetEventPlane(xQVecBPos, yQVecBPos, harmonic); - float const epBNegs = epHelper.GetEventPlane(xQVecBNeg, yQVecBNeg, harmonic); - float const epBTots = epHelper.GetEventPlane(xQVecBTot, yQVecBTot, harmonic); + float const epBPoss = epHelper.GetEventPlane(xQVecTPCPos, yQVecTPCPos, harmonic); + float const epBNegs = epHelper.GetEventPlane(xQVecTPCNeg, yQVecTPCNeg, harmonic); + float const epBTots = epHelper.GetEventPlane(xQVecTPCAll, yQVecTPCAll, harmonic); registry.fill(HIST("epReso/hEpResoFT0cFT0a"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epFT0a, harmonic))); registry.fill(HIST("epReso/hEpResoFT0cFV0a"), centrality, std::cos(harmonic * getDeltaPsiInRange(epFT0c, epFV0a, harmonic))); diff --git a/PWGHF/D2H/Utils/utilsFlow.h b/PWGHF/D2H/Utils/utilsFlow.h index 7efaedcdda3..911fc60dd2c 100644 --- a/PWGHF/D2H/Utils/utilsFlow.h +++ b/PWGHF/D2H/Utils/utilsFlow.h @@ -201,6 +201,56 @@ std::array getQvec(TCollision const& collision, const int qvecEst) } return std::array{-999.f, -999.f, -999.f}; } + +/// Get the Ese Q vector choosing your favourite estimator +/// \param collision is the collision with the Q vector information +/// \param qvecEst is the chosen Q-vector estimator +template +std::array getEseQvec(TCollision const& collision, const int qvecEst) +{ + switch (qvecEst) { + case QvecEstimator::FV0A: + if constexpr (HasQvecFV0A) { + return std::array{collision.eseQvecFV0ARe(), collision.eseQvecFV0AIm(), collision.sumAmplFV0A()}; + } + break; + case QvecEstimator::FT0A: + if constexpr (HasQvecFT0A) { + return std::array{collision.eseQvecFT0ARe(), collision.eseQvecFT0AIm(), collision.sumAmplFT0A()}; + } + break; + case QvecEstimator::FT0C: + if constexpr (HasQvecFT0C) { + return std::array{collision.eseQvecFT0CRe(), collision.eseQvecFT0CIm(), collision.sumAmplFT0C()}; + } + break; + case QvecEstimator::FT0M: + if constexpr (HasQvecFT0M) { + return std::array{collision.eseQvecFT0MRe(), collision.eseQvecFT0MIm(), collision.sumAmplFT0M()}; + } + break; + case QvecEstimator::TPCPos: + if constexpr (HasQvecTPCpos) { + return std::array{collision.eseQvecTPCposRe(), collision.eseQvecTPCposIm(), static_cast(collision.nTrkTPCpos())}; + } + break; + case QvecEstimator::TPCNeg: + if constexpr (HasQvecTPCneg) { + return std::array{collision.eseQvecTPCnegRe(), collision.eseQvecTPCnegIm(), static_cast(collision.nTrkTPCneg())}; + } + break; + case QvecEstimator::TPCTot: + if constexpr (HasQvecTPCtot) { + return std::array{collision.eseQvecTPCallRe(), collision.eseQvecTPCallIm(), static_cast(collision.nTrkTPCall())}; + } + break; + default: + LOGP(fatal, "Q-vector estimator not valid. Please choose between FV0A, FT0M, FT0A, FT0C, TPCPos, TPCNeg, TPCTot"); + break; + } + return std::array{-999.f, -999.f, -999.f}; +} + } // namespace hf_flow_utils } // namespace o2::analysis