diff --git a/Common/Core/FFitWeights.cxx b/Common/Core/FFitWeights.cxx index 9c98479627d..3a92114f48c 100644 --- a/Common/Core/FFitWeights.cxx +++ b/Common/Core/FFitWeights.cxx @@ -16,6 +16,10 @@ #include "FFitWeights.h" +#include +#include +#include + #include ClassImp(FFitWeights) @@ -43,20 +47,20 @@ FFitWeights::~FFitWeights() delete qAxis; }; -void FFitWeights::Init() +void FFitWeights::init() { fW_data = new TObjArray(); fW_data->SetName("FFitWeights_Data"); fW_data->SetOwner(kTRUE); if (!qAxis) - this->SetBinAxis(500, 0, 25); + this->setBinAxis(500, 0, 25); for (const auto& qn : qnTYPE) { - fW_data->Add(new TH2D(this->GetQName(qn.first, qn.second.c_str()), this->GetAxisName(qn.first, qn.second.c_str()), CentBin, 0, CentBin, qAxis->GetNbins(), qAxis->GetXmin(), qAxis->GetXmax())); + fW_data->Add(new TH2D(this->getQName(qn.first, qn.second.c_str()), this->getAxisName(qn.first, qn.second.c_str()), CentBin, 0, CentBin, qAxis->GetNbins(), qAxis->GetXmin(), qAxis->GetXmax())); } }; -void FFitWeights::Fill(float centrality, float qn, int nh, const char* pf) +void FFitWeights::fillWeights(float centrality, float qn, int nh, const char* pf) { TObjArray* tar{nullptr}; @@ -64,9 +68,9 @@ void FFitWeights::Fill(float centrality, float qn, int nh, const char* pf) if (!tar) return; - TH2D* th2 = reinterpret_cast(tar->FindObject(this->GetQName(nh, pf))); + TH2D* th2 = reinterpret_cast(tar->FindObject(this->getQName(nh, pf))); if (!th2) { - tar->Add(new TH2D(this->GetQName(nh, pf), this->GetAxisName(nh, pf), CentBin, 0, CentBin, qAxis->GetNbins(), qAxis->GetXmin(), qAxis->GetXmax())); + tar->Add(new TH2D(this->getQName(nh, pf), this->getAxisName(nh, pf), CentBin, 0, CentBin, qAxis->GetNbins(), qAxis->GetXmin(), qAxis->GetXmax())); th2 = reinterpret_cast(tar->At(tar->GetEntries() - 1)); } th2->Fill(centrality, qn); @@ -83,12 +87,12 @@ Long64_t FFitWeights::Merge(TCollection* collist) FFitWeights* l_w = 0; TIter all_w(collist); while ((l_w = (reinterpret_cast(all_w())))) { - AddArray(fW_data, l_w->GetDataArray()); + addArray(fW_data, l_w->getDataArray()); nmerged++; } return nmerged; }; -void FFitWeights::AddArray(TObjArray* targ, TObjArray* sour) +void FFitWeights::addArray(TObjArray* targ, TObjArray* sour) { if (!sour) { printf("Source array does not exist!\n"); @@ -107,7 +111,7 @@ void FFitWeights::AddArray(TObjArray* targ, TObjArray* sour) } }; -void FFitWeights::qSelectionSpline(std::vector nhv, std::vector stv) /* only execute OFFLINE */ +void FFitWeights::qSelection(std::vector nhv, std::vector stv) /* only execute OFFLINE */ { TObjArray* tar{nullptr}; @@ -117,15 +121,15 @@ void FFitWeights::qSelectionSpline(std::vector nhv, std::vector(tar->FindObject(this->GetQName(nh, pf.c_str()))); + TH2D* th2{reinterpret_cast(tar->FindObject(this->getQName(nh, pf.c_str())))}; if (!th2) { printf("qh not found!\n"); return; } - TH1D* tmp = nullptr; - TGraph* tmpgr = nullptr; - TSpline3* spline = nullptr; + TH1D* tmp{nullptr}; + TGraph* tmpgr{nullptr}; + // TSpline3* spline = nullptr; for (int iSP{0}; iSP < 90; iSP++) { tmp = th2->ProjectionY(Form("q%i_%i_%i", nh, iSP, iSP + 1), iSP + 1, iSP + 1); std::vector xq(nResolution); @@ -134,15 +138,16 @@ void FFitWeights::qSelectionSpline(std::vector nhv, std::vector(i + 1) / static_cast(nResolution); tmp->GetQuantiles(nResolution, yq.data(), xq.data()); tmpgr = new TGraph(nResolution, yq.data(), xq.data()); - spline = new TSpline3(Form("sp_q%i%s_%i", nh, pf.c_str(), iSP), tmpgr); - spline->SetName(Form("sp_q%i%s_%i", nh, pf.c_str(), iSP)); - fW_data->Add(spline); + tmpgr->SetName(Form("sp_q%i%s_%i", nh, pf.c_str(), iSP)); + // spline = new TSpline3(Form("sp_q%i%s_%i", nh, pf.c_str(), iSP), tmpgr); + // spline->SetName(Form("sp_q%i%s_%i", nh, pf.c_str(), iSP)); + fW_data->Add(tmpgr); } } } }; -float FFitWeights::EvalSplines(float centr, const float& dqn, const int nh, const char* pf) +float FFitWeights::eval(float centr, const float& dqn, const int nh, const char* pf) { TObjArray* tar{nullptr}; @@ -150,19 +155,19 @@ float FFitWeights::EvalSplines(float centr, const float& dqn, const int nh, cons if (!tar) return -1; - int isp = static_cast(centr); + int isp{static_cast(centr)}; if (isp < 0 || isp > 90) { return -1; } - TSpline3* spline = nullptr; - spline = reinterpret_cast(tar->FindObject(Form("sp_q%i%s_%i", nh, pf, isp))); + TGraph* spline{nullptr}; + spline = reinterpret_cast(tar->FindObject(Form("sp_q%i%s_%i", nh, pf, isp))); if (!spline) { return -1; } - float qn_val = 100. * spline->Eval(dqn); - if (qn_val < 0) { + float qn_val{static_cast(100. * spline->Eval(dqn))}; + if (qn_val < 0 || qn_val > 100.05) { return -1; } diff --git a/Common/Core/FFitWeights.h b/Common/Core/FFitWeights.h index 0a3d285176e..c80165730f7 100644 --- a/Common/Core/FFitWeights.h +++ b/Common/Core/FFitWeights.h @@ -41,22 +41,23 @@ class FFitWeights : public TNamed explicit FFitWeights(const char* name); ~FFitWeights(); - void Init(); - void Fill(float centrality, float qn, int nh, const char* pf = ""); - TObjArray* GetDataArray() { return fW_data; } + void init(); + void fillWeights(float centrality, float qn, int nh, const char* pf = ""); + TObjArray* getDataArray() { return fW_data; } - void SetCentBin(int bin) { CentBin = bin; } - void SetBinAxis(int bin, float min, float max) + void setCentBin(int bin) { CentBin = bin; } + void setBinAxis(int bin, float min, float max) { qAxis = new TAxis(bin, min, max); } - TAxis* GetqVecAx() { return qAxis; } + TAxis* getqVecAx() { return qAxis; } Long64_t Merge(TCollection* collist); - void qSelectionSpline(std::vector nhv, std::vector stv); - float EvalSplines(float centr, const float& dqn, const int nh, const char* pf = ""); - void SetResolution(int res) { nResolution = res; } - void SetQnType(std::vector> qninp) { qnTYPE = qninp; } + void qSelection(std::vector nhv, std::vector stv); + float eval(float centr, const float& dqn, const int nh, const char* pf = ""); + void setResolution(int res) { nResolution = res; } + int getResolution() const { return nResolution; } + void setQnType(std::vector> qninp) { qnTYPE = qninp; } private: TObjArray* fW_data; @@ -67,15 +68,15 @@ class FFitWeights : public TNamed std::vector> qnTYPE; - const char* GetQName(const int nh, const char* pf = "") + const char* getQName(const int nh, const char* pf = "") { return Form("q%i%s", nh, pf); }; - const char* GetAxisName(const int nh, const char* pf = "") + const char* getAxisName(const int nh, const char* pf = "") { return Form(";Centrality;q_{%i}^{%s}", nh, pf); }; - void AddArray(TObjArray* targ, TObjArray* sour); + void addArray(TObjArray* targ, TObjArray* sour); ClassDef(FFitWeights, 1); // calibration class }; diff --git a/Common/DataModel/EseTable.h b/Common/DataModel/EseTable.h index ca9bd41357e..68ffdff450a 100644 --- a/Common/DataModel/EseTable.h +++ b/Common/DataModel/EseTable.h @@ -29,20 +29,23 @@ namespace q_vector DECLARE_SOA_COLUMN(QPERCFT0C, qPERCFT0C, std::vector); DECLARE_SOA_COLUMN(QPERCFT0A, qPERCFT0A, std::vector); DECLARE_SOA_COLUMN(QPERCFV0A, qPERCFV0A, std::vector); -DECLARE_SOA_COLUMN(QPERCTPC, qPERCTPC, std::vector); -DECLARE_SOA_COLUMN(FESECOL, fESECOL, std::vector); +DECLARE_SOA_COLUMN(QPERCTPCall, qPERCTPCall, std::vector); +DECLARE_SOA_COLUMN(QPERCTPCneg, qPERCTPCneg, std::vector); +DECLARE_SOA_COLUMN(QPERCTPCpos, qPERCTPCpos, std::vector); } // namespace q_vector DECLARE_SOA_TABLE(QPercentileFT0Cs, "AOD", "QPERCENTILEFT0C", q_vector::QPERCFT0C); DECLARE_SOA_TABLE(QPercentileFT0As, "AOD", "QPERCENTILEFT0A", q_vector::QPERCFT0A); DECLARE_SOA_TABLE(QPercentileFV0As, "AOD", "QPERCENTILEFV0A", q_vector::QPERCFV0A); -DECLARE_SOA_TABLE(QPercentileTPCs, "AOD", "QPERCENTILETPC", q_vector::QPERCTPC); -DECLARE_SOA_TABLE(FEseCols, "AOD", "FEVENTSHAPE", q_vector::FESECOL); +DECLARE_SOA_TABLE(QPercentileTPCalls, "AOD", "QPERCENTILETPCall", q_vector::QPERCTPCall); +DECLARE_SOA_TABLE(QPercentileTPCnegs, "AOD", "QPERCENTILETPCneg", q_vector::QPERCTPCneg); +DECLARE_SOA_TABLE(QPercentileTPCposs, "AOD", "QPERCENTILETPCpos", q_vector::QPERCTPCpos); using QPercentileFT0C = QPercentileFT0Cs::iterator; using QPercentileFT0A = QPercentileFT0As::iterator; using QPercentileFV0A = QPercentileFV0As::iterator; -using QPercentileTPC = QPercentileTPCs::iterator; -using FEseCol = FEseCols::iterator; +using QPercentileTPCall = QPercentileTPCalls::iterator; +using QPercentileTPCneg = QPercentileTPCnegs::iterator; +using QPercentileTPCpos = QPercentileTPCposs::iterator; } // namespace o2::aod diff --git a/Common/TableProducer/eseTableProducer.cxx b/Common/TableProducer/eseTableProducer.cxx index 34f534dc700..098bfbe400f 100644 --- a/Common/TableProducer/eseTableProducer.cxx +++ b/Common/TableProducer/eseTableProducer.cxx @@ -21,6 +21,8 @@ #include #include #include +#include +#include #include "Framework/ASoA.h" #include "Framework/AnalysisDataModel.h" @@ -40,33 +42,50 @@ using namespace o2; using namespace o2::framework; -using CollWithMults = soa::Join; +using CollWithMults = soa::Join; struct EseTableProducer { Produces qPercsFT0C; Produces qPercsFT0A; Produces qPercsFV0A; - Produces qPercsTPC; - Produces fEseCol; + Produces qPercsTPCall; + Produces qPercsTPCneg; + Produces qPercsTPCpos; OutputObj FFitObj{FFitWeights("weights")}; HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; - Configurable cfgESE{"cfgESE", 1, "ese activation step: false = no ese, true = evaluate splines and fill table"}; + Configurable cfgESE{"cfgESE", 1, "ese activation step: false = no ese, true = evaluate qSelection and fill table"}; Configurable cfgEsePath{"cfgEsePath", "Users/j/joachiha/ESE/local/ffitsplines", "CCDB path for ese splines"}; - Configurable cfgFT0C{"cfgFT0C", 1, "FT0C flag"}; - Configurable cfgFT0A{"cfgFT0A", 0, "FT0A flag"}; - Configurable cfgFV0A{"cfgFV0A", 0, "FV0A flag"}; - Configurable cfgTPC{"cfgTPC", 0, "TPC flag"}; - Configurable> cfgLoopHarmonics{"cfgLoopHarmonics", {2, 3}, "Harmonics to loop over when filling and evaluating splines"}; - Configurable cfgCentEst{"cfgCentEst", "FT0C", "centrality estimator"}; + Configurable> cfgDetectors{"cfgDetectors", {"FT0C"}, "detectors to loop over: ['FT0C', 'FT0A', 'FV0A', 'TPCall', 'TPCneg', 'TPCpos']"}; + Configurable> cfgLoopHarmonics{"cfgLoopHarmonics", {2, 3}, "Harmonics to loop over when filling and evaluating q-Selection"}; Configurable> cfgaxisqn{"cfgaxisqn", {500, 0, 25}, "q_n amplitude range"}; - Configurable cfgnResolution{"cfgnResolution", 3000, "resolution of splines"}; + Configurable cfgnResolution{"cfgnResolution", 3000, "resolution of q-Selection"}; + + Configurable cfgnTotalSystem{"cfgnTotalSystem", 7, "total qvector number // look in Qvector table for this number"}; + Configurable cfgnCorrLevel{"cfgnCorrLevel", 3, "QVector step: 0 = no corr, 1 = rect, 2 = twist, 3 = full"}; int runNumber{-1}; - FFitWeights* splines{nullptr}; + enum class DetID { FT0C, + FT0A, + FT0M, + FV0A, + TPCpos, + TPCneg, + TPCall }; + + std::unordered_map detMap = { + {"FT0C", DetID::FT0C}, + {"FT0A", DetID::FT0A}, + {"FT0M", DetID::FT0M}, + {"FV0A", DetID::FV0A}, + {"TPCpos", DetID::TPCpos}, + {"TPCneg", DetID::TPCneg}, + {"TPCall", DetID::TPCall}}; + + FFitWeights* qSelection{nullptr}; Service ccdb; @@ -85,26 +104,16 @@ struct EseTableProducer { int64_t now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); ccdb->setCreatedNotAfter(now); - std::vector vecStr{}; - if (cfgFT0C) - vecStr.push_back("FT0C"); - if (cfgFT0A) - vecStr.push_back("FT0A"); - if (cfgFV0A) - vecStr.push_back("FV0A"); - if (cfgTPC) - vecStr.push_back("TPC"); - std::vector> veccfg; for (std::size_t i{0}; i < cfgLoopHarmonics->size(); i++) { - for (const auto& j : vecStr) { - veccfg.push_back({cfgLoopHarmonics->at(i), j}); + for (std::size_t j{0}; j < cfgDetectors->size(); j++) { + veccfg.push_back({cfgLoopHarmonics->at(i), cfgDetectors->at(j)}); } } - FFitObj->SetBinAxis(cfgaxisqn->at(0), cfgaxisqn->at(1), cfgaxisqn->at(2)); - FFitObj->SetResolution(cfgnResolution); - FFitObj->SetQnType(veccfg); - FFitObj->Init(); + FFitObj->setBinAxis(cfgaxisqn->at(0), cfgaxisqn->at(1), cfgaxisqn->at(2)); + FFitObj->setResolution(cfgnResolution); + FFitObj->setQnType(veccfg); + FFitObj->init(); } void initCCDB(aod::BCsWithTimestamps::iterator const& bc) @@ -113,10 +122,10 @@ struct EseTableProducer { auto timestamp = bc.timestamp(); if (cfgESE) { - splines = ccdb->getForTimeStamp(cfgEsePath, timestamp); - if (!splines) - LOGF(fatal, "failed loading splines with ese flag"); - LOGF(info, "successfully loaded splines"); + qSelection = ccdb->getForTimeStamp(cfgEsePath, timestamp); + if (!qSelection) + LOGF(fatal, "failed loading qSelection with ese flag"); + LOGF(info, "successfully loaded qSelection"); } } @@ -130,76 +139,88 @@ struct EseTableProducer { return qn; } - bool validQvec(const float& qVec) + constexpr int detIDN(const DetID id) { - if (qVec == 999. || qVec == -999) { - return false; - } else { - return true; + switch (id) { + case DetID::FT0C: + return 0; + case DetID::FT0A: + return 1; + case DetID::FT0M: + return 2; + case DetID::FV0A: + return 3; + case DetID::TPCpos: + return 4; + case DetID::TPCneg: + return 5; + case DetID::TPCall: + return 6; } - }; + return -1; + } - void doSpline(float& splineVal, int& eseval, const float& centr, const float& nHarm, const char* pf, const auto& QX, const auto& QY, const auto& Ampl) + void doSpline(float& splineVal, const float& centr, const float& nHarm, const char* pf, const auto& QX, const auto& QY, const auto& sumAmpl) { - if (validQvec(QX[nHarm - 2]) && validQvec(QY[nHarm - 2]) && Ampl > 1e-8) { - float qnval = Calcqn(QX[nHarm - 2] * Ampl, QY[nHarm - 2] * Ampl, Ampl); - FFitObj->Fill(centr, qnval, nHarm, pf); + if (sumAmpl > 1e-8) { + float qnval = Calcqn(QX * sumAmpl, QY * sumAmpl, sumAmpl); + FFitObj->fillWeights(centr, qnval, nHarm, pf); if (cfgESE) { - splineVal = splines->EvalSplines(centr, qnval, nHarm, pf); - eseval = cfgFT0C ? 1 : 0; + splineVal = qSelection->eval(centr, qnval, nHarm, pf); } } } + template + std::tuple getVectors(const C& col, const int& nHarm, const DetID& id) + { + const int detId = detIDN(id); + const int detInd{detId * 4 + cfgnTotalSystem * 4 * (nHarm - 2)}; + const auto Qx{col.qvecRe()[detInd + cfgnCorrLevel]}; + const auto Qy{col.qvecIm()[detInd + cfgnCorrLevel]}; + const auto sumAmpl{col.qvecAmp()[detId]}; + return {Qx, Qy, sumAmpl}; + } + template void calculateESE(T const& collision, std::vector& qnpFT0C, std::vector& qnpFT0A, std::vector& qnpFV0A, - std::vector& qnpTPC, - std::vector& fIsEseAvailable) + std::vector& qnpTPCall, + std::vector& qnpTPCneg, + std::vector& qnpTPCpos) { - const float centrality = collision.centFT0C(); - registry.fill(HIST("hESEstat"), 0.5); - for (std::size_t i{0}; i < cfgLoopHarmonics->size(); i++) { - float splineValFT0C{-1.0}; - float splineValFT0A{-1.0}; - float splineValFV0A{-1.0}; - qnpTPC.push_back(-1.0); /* not implemented yet */ - int eseAvailable{0}; - - int nHarm = cfgLoopHarmonics->at(i); - if (cfgFT0C) { - const auto QxFT0C_Qvec = collision.qvecFT0CReVec(); - const auto QyFT0C_Qvec = collision.qvecFT0CImVec(); - const auto SumAmplFT0C = collision.sumAmplFT0C(); - doSpline(splineValFT0C, eseAvailable, centrality, nHarm, "FT0C", QxFT0C_Qvec, QyFT0C_Qvec, SumAmplFT0C); - if (i == 0) - registry.fill(HIST("hESEstat"), 1.5); - } - qnpFT0C.push_back(splineValFT0C); - fIsEseAvailable.push_back(eseAvailable); - - if (cfgFT0A) { - const auto QxFT0A_Qvec = collision.qvecFT0AReVec(); - const auto QyFT0A_Qvec = collision.qvecFT0AImVec(); - const auto SumAmplFT0A = collision.sumAmplFT0A(); - doSpline(splineValFT0A, eseAvailable, centrality, nHarm, "FT0A", QxFT0A_Qvec, QyFT0A_Qvec, SumAmplFT0A); - if (i == 0) - registry.fill(HIST("hESEstat"), 2.5); - } - qnpFT0A.push_back(splineValFT0A); - - if (cfgFV0A) { - const auto QxFV0A_Qvec = collision.qvecFV0AReVec(); - const auto QyFV0A_Qvec = collision.qvecFV0AImVec(); - const auto SumAmplFV0A = collision.sumAmplFV0A(); - doSpline(splineValFV0A, eseAvailable, centrality, nHarm, "FV0A", QxFV0A_Qvec, QyFV0A_Qvec, SumAmplFV0A); - if (i == 0) - registry.fill(HIST("hESEstat"), 3.5); + float counter{0.5}; + registry.fill(HIST("hESEstat"), counter++); + + std::unordered_map*> vMap{ + {"FT0C", &qnpFT0C}, + {"FT0A", &qnpFT0A}, + {"FV0A", &qnpFV0A}, + {"TPCall", &qnpTPCall}, + {"TPCneg", &qnpTPCneg}, + {"TPCpos", &qnpTPCpos}}; + + for (std::size_t j{0}; j < cfgDetectors->size(); j++) { + const auto det{cfgDetectors->at(j)}; + const auto iter{detMap.find(det)}; + float splineVal{-1.0}; + + if (iter != detMap.end()) { + for (std::size_t i{0}; i < cfgLoopHarmonics->size(); i++) { + const int nHarm{cfgLoopHarmonics->at(i)}; + const auto [qxt, qyt, st] = getVectors(collision, nHarm, iter->second); + doSpline(splineVal, centrality, nHarm, det.c_str(), qxt, qyt, st); + if (i == 0) + registry.fill(HIST("hESEstat"), counter++); + + if (vMap.find(det) != vMap.end()) { + vMap[det]->push_back(splineVal); + } + } } - qnpFV0A.push_back(splineValFV0A); } }; @@ -211,24 +232,25 @@ struct EseTableProducer { std::vector qnpFT0C{}; std::vector qnpFT0A{}; std::vector qnpFV0A{}; - std::vector qnpTPC{}; - - std::vector fIsEseAvailable{}; + std::vector qnpTPCall{}; + std::vector qnpTPCneg{}; + std::vector qnpTPCpos{}; - auto bc = collision.bc_as(); - int currentRun = bc.runNumber(); + auto bc{collision.bc_as()}; + int currentRun{bc.runNumber()}; if (runNumber != currentRun) { runNumber = currentRun; initCCDB(bc); } registry.fill(HIST("hEventCounter"), counter++); - calculateESE(collision, qnpFT0C, qnpFT0A, qnpFV0A, qnpTPC, fIsEseAvailable); + calculateESE(collision, qnpFT0C, qnpFT0A, qnpFV0A, qnpTPCall, qnpTPCneg, qnpTPCpos); qPercsFT0C(qnpFT0C); qPercsFT0A(qnpFT0A); qPercsFV0A(qnpFV0A); - qPercsTPC(qnpTPC); - fEseCol(fIsEseAvailable); + qPercsTPCall(qnpTPCall); + qPercsTPCneg(qnpTPCneg); + qPercsTPCpos(qnpTPCpos); registry.fill(HIST("hEventCounter"), counter++); } PROCESS_SWITCH(EseTableProducer, processESE, "proccess q vectors to calculate reduced q-vector", true);