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..c801d4c4ab3 100644 --- a/Common/TableProducer/qVectorsTable.cxx +++ b/Common/TableProducer/qVectorsTable.cxx @@ -50,6 +50,7 @@ #include #include +#include #include #include #include @@ -72,7 +73,21 @@ struct qVectorsTable { kFV0A, kTPCpos, kTPCneg, - kTPCall + kTPCall, + kNDetectors + }; + enum Corrections { + kNoCorr = 0, + kRecenter, + kTwist, + kRescale, + kNCorrections + }; + enum MultNorms { + kNoNorm = 0, + kScalarProd, + kEsE, + kMultNormTypes }; // Configurables. @@ -98,7 +113,8 @@ 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 cfgQvecEseCalibPath{"cfgQvecEseCalibPath", "Analysis/EventPlane/QVecEseCorrections", "CCDB path for EsE Q-vector calibration constants"}; Configurable cfgShiftCorr{"cfgShiftCorr", false, "configurable flag for shift correction"}; Configurable cfgShiftPath{"cfgShiftPath", "", "CCDB path for shift correction"}; @@ -112,6 +128,7 @@ struct qVectorsTable { 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; @@ -131,6 +148,24 @@ struct qVectorsTable { 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{}; @@ -150,7 +185,8 @@ struct qVectorsTable { int runNumber{-1}; float cent; - std::vector objQvec{}; + std::vector corrsQvecSp{}; + std::vector corrsQvecEse{}; std::vector shiftprofile{}; // Deprecated, will be removed in future after transition time // @@ -257,19 +293,34 @@ 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 = cfgQvecEseCalibPath; + fullPath += "/eseq"; + fullPath += std::to_string(ind); + auto modeCorrQvecEse = getForTsOrRun(fullPath, timestamp, runnumber); + if (!modeCorrQvecEse) { + fullPath = cfgQvecCalibPath; // cfgQvecEseCalibPath; + fullPath += "/v2"; // "/eseq2"; + modeCorrQvecEse = getForTsOrRun(fullPath, timestamp, runnumber); + } + corrsQvecEse.push_back(modeCorrQvecEse); } if (cfgShiftCorr) { @@ -347,16 +398,175 @@ 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 QvecReRaw, + std::vector QvecImRaw, + std::vector& QvecAmp, + MultNorms normType) + { + for (std::size_t i = 0; i < kNDetectors; i++) { + float qVecDetReNorm{999.}, qVecDetImNorm{999.}; + if (QvecAmp[i] > 1e-8) { + 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; + } + std::cout << "[NORMALIZED] " << i << " Re: " << qVecDetReNorm << ", Im: " << qVecDetImNorm << ", amp: " << QvecAmp[i] << std::endl; + } + for (int iCorr = 0; iCorr < Corrections::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, int nMode) + { + int nCorrections = static_cast(Corrections::kNCorrections); + if (cent < cfgMaxCentrality) { + for (auto i{0u}; i < kTPCall + 1; i++) { + int idxDet = i * Corrections::kNCorrections; + helperEP.DoRecenter(qvecRe[idxDet + Corrections::kRecenter], qvecIm[idxDet + Corrections::kRecenter], + histsCorrs->GetBinContent(static_cast(cent) + 1, 1, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 2, i + 1)); + + helperEP.DoRecenter(qvecRe[idxDet + Corrections::kTwist], qvecIm[idxDet + Corrections::kTwist], + histsCorrs->GetBinContent(static_cast(cent) + 1, 1, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 2, i + 1)); + helperEP.DoTwist(qvecRe[idxDet + Corrections::kTwist], qvecIm[idxDet + Corrections::kTwist], + histsCorrs->GetBinContent(static_cast(cent) + 1, 3, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 4, i + 1)); + + helperEP.DoRecenter(qvecRe[idxDet + Corrections::kRescale], qvecIm[idxDet + Corrections::kRescale], + histsCorrs->GetBinContent(static_cast(cent) + 1, 1, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 2, i + 1)); + helperEP.DoTwist(qvecRe[idxDet + Corrections::kRescale], qvecIm[idxDet + Corrections::kRescale], + histsCorrs->GetBinContent(static_cast(cent) + 1, 3, i + 1), histsCorrs->GetBinContent(static_cast(cent) + 1, 4, i + 1)); + helperEP.DoRescale(qvecRe[idxDet + Corrections::kRescale], qvecIm[idxDet + Corrections::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 = TMath::ATan2(qvecIm[kFT0C * nCorrections + Corrections::kRescale], qvecRe[kFT0C * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefFT0A = TMath::ATan2(qvecIm[kFT0A * nCorrections + Corrections::kRescale], qvecRe[kFT0A * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefFT0M = TMath::ATan2(qvecIm[kFT0M * nCorrections + Corrections::kRescale], qvecRe[kFT0M * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefFV0A = TMath::ATan2(qvecIm[kFV0A * nCorrections + Corrections::kRescale], qvecRe[kFV0A * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefTPCpos = TMath::ATan2(qvecIm[kTPCpos * nCorrections + Corrections::kRescale], qvecRe[kTPCpos * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefTPCneg = TMath::ATan2(qvecIm[kTPCneg * nCorrections + Corrections::kRescale], qvecRe[kTPCneg * nCorrections + Corrections::kRescale]) / static_cast(nMode); + auto psidefTPCall = TMath::ATan2(qvecIm[kTPCall * nCorrections + Corrections::kRescale], qvecRe[kTPCall * nCorrections + Corrections::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 * 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[kFT0C * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0C) - qvecIm[kFT0C * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0C); + float qvecImShiftedFT0C = qvecRe[kFT0C * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0C) + qvecIm[kFT0C * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0C); + float qvecReShiftedFT0A = qvecRe[kFT0A * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0A) - qvecIm[kFT0A * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0A); + float qvecImShiftedFT0A = qvecRe[kFT0A * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0A) + qvecIm[kFT0A * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0A); + float qvecReShiftedFT0M = qvecRe[kFT0M * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0M) - qvecIm[kFT0M * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0M); + float qvecImShiftedFT0M = qvecRe[kFT0M * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFT0M) + qvecIm[kFT0M * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFT0M); + float qvecReShiftedFV0A = qvecRe[kFV0A * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFV0A) - qvecIm[kFV0A * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFV0A); + float qvecImShiftedFV0A = qvecRe[kFV0A * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiFV0A) + qvecIm[kFV0A * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiFV0A); + float qvecReShiftedTPCpos = qvecRe[kTPCpos * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCpos) - qvecIm[kTPCpos * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCpos); + float qvecImShiftedTPCpos = qvecRe[kTPCpos * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCpos) + qvecIm[kTPCpos * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCpos); + float qvecReShiftedTPCneg = qvecRe[kTPCneg * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCneg) - qvecIm[kTPCneg * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCneg); + float qvecImShiftedTPCneg = qvecRe[kTPCneg * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCneg) + qvecIm[kTPCneg * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCneg); + float qvecReShiftedTPCall = qvecRe[kTPCall * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCall) - qvecIm[kTPCall * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCall); + float qvecImShiftedTPCall = qvecRe[kTPCall * nCorrections + Corrections::kRescale] * TMath::Sin(deltapsiTPCall) + qvecIm[kTPCall * nCorrections + Corrections::kRescale] * TMath::Cos(deltapsiTPCall); + + qvecRe[kFT0C * nCorrections + Corrections::kRescale] = qvecReShiftedFT0C; + qvecIm[kFT0C * nCorrections + Corrections::kRescale] = qvecImShiftedFT0C; + qvecRe[kFT0A * nCorrections + Corrections::kRescale] = qvecReShiftedFT0A; + qvecIm[kFT0A * nCorrections + Corrections::kRescale] = qvecImShiftedFT0A; + qvecRe[kFT0M * nCorrections + Corrections::kRescale] = qvecReShiftedFT0M; + qvecIm[kFT0M * nCorrections + Corrections::kRescale] = qvecImShiftedFT0M; + qvecRe[kFV0A * nCorrections + Corrections::kRescale] = qvecReShiftedFV0A; + qvecIm[kFV0A * nCorrections + Corrections::kRescale] = qvecImShiftedFV0A; + qvecRe[kTPCpos * nCorrections + Corrections::kRescale] = qvecReShiftedTPCpos; + qvecIm[kTPCpos * nCorrections + Corrections::kRescale] = qvecImShiftedTPCpos; + qvecRe[kTPCneg * nCorrections + Corrections::kRescale] = qvecReShiftedTPCneg; + qvecIm[kTPCneg * nCorrections + Corrections::kRescale] = qvecImShiftedTPCneg; + qvecRe[kTPCall * nCorrections + Corrections::kRescale] = qvecReShiftedTPCall; + qvecIm[kTPCall * nCorrections + Corrections::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 track 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 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) { - 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.}; + 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); @@ -376,17 +586,13 @@ 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(); } - } else { - qVectFT0A[0] = 999.; - qVectFT0A[1] = 999.; } if (useDetector["QvectorFT0Cs"]) { @@ -398,65 +604,39 @@ 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.; } - } 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.; + if (sumAmplFT0M > 1e-8 && useDetector["QvectorFT0Ms"]) { + qVectFT0M[0] = QvecFT0M.Re(); + qVectFT0M[1] = QvecFT0M.Im(); + } } - } 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 > 1e-8) { + qVectFV0A[0] = QvecDet.Re(); + qVectFV0A[1] = QvecDet.Im(); + } } - } else { - qVectFV0A[0] = -999.; - qVectFV0A[1] = -999.; } int nTrkTPCpos = 0; @@ -474,77 +654,40 @@ 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); + 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); + 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); + 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]); - } + 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); @@ -553,31 +696,57 @@ struct qVectorsTable { QvecAmp.push_back(static_cast(nTrkTPCpos)); QvecAmp.push_back(static_cast(nTrkTPCneg)); QvecAmp.push_back(static_cast(nTrkTPCall)); + + LOG(info) << "[RAW] qVectFT0A: " << qVectFT0A[0] << ", " << qVectFT0A[1] << ", ampl: " << sumAmplFT0A; + LOG(info) << "[RAW] qVectFT0C: " << qVectFT0C[0] << ", " << qVectFT0C[1] << ", ampl: " << sumAmplFT0C; + LOG(info) << "[RAW] qVectFT0M: " << qVectFT0M[0] << ", " << qVectFT0M[1] << ", ampl: " << sumAmplFT0M; + LOG(info) << "[RAW] qVectFV0A: " << qVectFV0A[0] << ", " << qVectFV0A[1] << ", ampl: " << sumAmplFV0A; + LOG(info) << "[RAW] qVectTPCpos: " << qVectTPCpos[0] << ", " << qVectTPCpos[1] << ", nTrk: " << nTrkTPCpos; + LOG(info) << "[RAW] qVectTPCneg: " << qVectTPCneg[0] << ", " << qVectTPCneg[1] << ", nTrk: " << nTrkTPCneg; + LOG(info) << "[RAW] qVectTPCall: " << qVectTPCall[0] << ", " << qVectTPCall[1] << ", nTrk: " << nTrkTPCall; } void process(MyCollisions::iterator const& coll, aod::BCsWithTimestamps const&, aod::FT0s const&, aod::FV0As const&, MyTracks const& tracks) { + LOG(info) << "---------------------------- Processing Event ---------------------------"; 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 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(); @@ -595,162 +764,139 @@ struct qVectorsTable { cent = 110.; 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{}; + CalQvec(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], nMode); + // Add to summary vector + qvecReSp.insert(qvecReSp.end(), nModeQvecReSp.begin(), nModeQvecReSp.end()); + qvecImSp.insert(qvecImSp.end(), nModeQvecImSp.begin(), nModeQvecImSp.end()); + + // 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], nMode); + // Add to summary vector + qvecReEse.insert(qvecReEse.end(), nModeQvecReEse.begin(), nModeQvecReEse.end()); + qvecImEse.insert(qvecImEse.end(), nModeQvecImEse.begin(), nModeQvecImEse.end()); + + // Pick the desired correction level for the Q-vectors to be stored in the analysis table 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]); + int nCorrections = static_cast(Corrections::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]); + + 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]); } // 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); - - 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); + qVectorTPCall(IsCalibrated, qvecReTPCallSp.at(0), qvecImTPCallSp.at(0), qvecAmp[kTPCall], TrkTPCallLabel); + + // Debug prints of values after corrections + std::cout << "[CORRECTED] FT0C, Re: " << qvecReFT0CSp.at(0) << ", Im: " << qvecImFT0CSp.at(0) << std::endl; + std::cout << "[CORRECTED] FT0A, Re: " << qvecReFT0ASp.at(0) << ", Im: " << qvecImFT0ASp.at(0) << std::endl; + std::cout << "[CORRECTED] FT0M, Re: " << qvecReFT0MSp.at(0) << ", Im: " << qvecImFT0MSp.at(0) << std::endl; + std::cout << "[CORRECTED] FV0A, Re: " << qvecReFV0ASp.at(0) << ", Im: " << qvecImFV0ASp.at(0) << std::endl; + std::cout << "[CORRECTED] TPCpos, Re: " << qvecReTPCposSp.at(0) << ", Im: " << qvecImTPCposSp.at(0) << std::endl; + std::cout << "[CORRECTED] TPCneg, Re: " << qvecReTPCnegSp.at(0) << ", Im: " << qvecImTPCnegSp.at(0) << std::endl; + std::cout << "[CORRECTED] TPCall, Re: " << qvecReTPCallSp.at(0) << ", Im: " << qvecImTPCallSp.at(0) << std::endl; + + 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