99// granted to it by virtue of its status as an Intergovernmental Organization
1010// or submit itself to any jurisdiction.
1111
12+ // / \file GFW.cxx
13+ // / \brief Initialize regions and calculate n-particle correlations
14+ // / \author Emil Gorm Nielsen, NBI, emil.gorm.nielsen@cern.ch
15+
1216#include " GFW.h"
1317
1418#include " PWGCF/GenericFramework/Core/GFWPowerArray.h"
1519
20+ #include < Framework/Logger.h>
21+
1622#include < complex>
17- #include < cstdio >
23+ #include < cstddef >
1824#include < string>
1925#include < utility>
2026#include < vector>
2127
2228using std::complex ;
2329using std::pair;
2430using std::string;
25- using std::vector;
31+
32+ namespace
33+ {
34+ constexpr std::size_t TwoParticleOrder = 2 ;
35+ constexpr std::size_t ThreeParticleOrder = 3 ;
36+ } // namespace
2637
2738GFW ::GFW () : fInitialized (false ) {}
2839
40+ void GFW::Region::PrintStructure ()
41+ {
42+ LOGF (info, " %s: eta [%f.. %f]." , rName.c_str (), EtaMin, EtaMax);
43+ }
44+
2945GFW ::~GFW ()
3046{
3147 for (auto pItr = fCumulants .begin (); pItr != fCumulants .end (); ++pItr)
@@ -34,15 +50,15 @@ GFW::~GFW()
3450void GFW::AddRegion (const string& refName, double lEtaMin, double lEtaMax, int lNpT, int BitMask)
3551{
3652 if (lNpT < 1 ) {
37- printf ( " Number of pT bins cannot be less than 1! Not adding anything.\n " );
53+ LOGF (warning, " Number of pT bins cannot be less than 1! Not adding anything." );
3854 return ;
3955 }
4056 if (lEtaMin >= lEtaMax) {
41- printf ( " Eta min. cannot be more than eta max! Not adding...\n " );
57+ LOGF (warning, " Eta min. cannot be more than eta max! Not adding..." );
4258 return ;
4359 }
4460 if (refName == " " ) {
45- printf ( " Region must have a name!\n " );
61+ LOGF (warning, " Region must have a name!" );
4662 return ;
4763 }
4864 Region lOneRegion;
@@ -56,23 +72,26 @@ void GFW::AddRegion(const string& refName, double lEtaMin, double lEtaMax, int l
5672 lOneRegion.BitMask = BitMask; // Bit mask
5773 AddRegion (lOneRegion);
5874};
59- void GFW::AddRegion (const string& refName, const vector<int >& lNparVec, double lEtaMin, double lEtaMax, int lNpT, int BitMask)
75+ void GFW::AddRegion (string refName, std:: vector<int > lNparVec, double lEtaMin, double lEtaMax, int lNpT, int BitMask)
6076{
61- AddRegion (std::move (refName), lEtaMin, lEtaMax, lNpT, BitMask);
77+ const auto oldSize = fRegions .size ();
78+ AddRegion (refName, lEtaMin, lEtaMax, lNpT, BitMask);
79+ if (fRegions .size () == oldSize)
80+ return ;
6281 (fRegions .end () - 1 )->Nhar = static_cast <int >(lNparVec.size ());
6382 (fRegions .end () - 1 )->NparVec = lNparVec;
6483 (fRegions .end () - 1 )->powsDefined = true ;
6584};
6685void GFW::AddRegion (string refName, int lNhar, int lNpar, double lEtaMin, double lEtaMax, int lNpT, int BitMask)
6786{
68- vector<int > tVec = {};
87+ std:: vector<int > tVec = {};
6988 for (int i = 0 ; i < lNhar; i++)
7089 tVec.push_back (lNpar);
7190 AddRegion (std::move (refName), tVec, lEtaMin, lEtaMax, lNpT, BitMask);
7291};
7392void GFW::AddRegion (string refName, int lNhar, int * lNparVec, double lEtaMin, double lEtaMax, int lNpT, int BitMask)
7493{
75- vector<int > tVec = {};
94+ std:: vector<int > tVec = {};
7695 for (int i = 0 ; i < lNhar; i++)
7796 tVec.push_back (lNparVec[i]);
7897 AddRegion (std::move (refName), tVec, lEtaMin, lEtaMax, lNpT, BitMask);
@@ -84,14 +103,13 @@ int GFW::CreateRegions()
84103 fCumulants .clear ();
85104 InitializePowerArrays ();
86105 if (fRegions .size () < 1 ) {
87- printf ( " No regions set. Skipping...\n " );
106+ LOGF (warning, " No regions set. Skipping..." );
88107 return 0 ;
89108 }
90109 int nRegions = 0 ;
91110 for (auto pItr = fRegions .begin (); pItr != fRegions .end (); pItr++) {
92- GFWCumulant* lCumulant = new GFWCumulant ();
93- lCumulant->CreateComplexVectorArrayVarPower (pItr->Nhar , pItr->NparVec , pItr->NpT );
94- fCumulants .push_back (*lCumulant);
111+ fCumulants .emplace_back ();
112+ fCumulants .back ().CreateComplexVectorArrayVarPower (pItr->Nhar , pItr->NparVec , pItr->NpT );
95113 ++nRegions;
96114 }
97115 if (nRegions)
@@ -114,22 +132,24 @@ complex<double> GFW::TwoRec(int n1, int n2, int p1, int p2, int ptbin, GFWCumula
114132 complex <double > formula = part1 * part2 - part3;
115133 return formula;
116134};
117- complex <double > GFW::RecursiveCorr (GFWCumulant* qpoi, GFWCumulant* qref, GFWCumulant* qol, int ptbin, vector<int >& hars)
135+ complex <double > GFW::RecursiveCorr (GFWCumulant* qpoi, GFWCumulant* qref, GFWCumulant* qol, int ptbin, std:: vector<int >& hars)
118136{
119- vector<int > pows;
137+ std:: vector<int > pows;
120138 for (int i = 0 ; i < static_cast <int >(hars.size ()); i++)
121139 pows.push_back (1 );
122140 return RecursiveCorr (qpoi, qref, qol, ptbin, hars, pows);
123141};
124142
125- complex <double > GFW::RecursiveCorr (GFWCumulant* qpoi, GFWCumulant* qref, GFWCumulant* qol, int ptbin, vector<int >& hars, vector<int >& pows)
143+ complex <double > GFW::RecursiveCorr (GFWCumulant* qpoi, GFWCumulant* qref, GFWCumulant* qol, int ptbin, std:: vector<int >& hars, std:: vector<int >& pows)
126144{
145+ if (hars.empty () || pows.empty ())
146+ return complex <double >(0 , 0 );
127147 if ((pows.at (0 ) != 1 ) && qol)
128148 qpoi = qol; // if the power of POI is not unity, then always use overlap (if defined).
129149 // Only valid for 1 particle of interest though!
130- if (hars.size () < 2 )
150+ if (hars.size () < TwoParticleOrder )
131151 return qpoi->Vec (hars.at (0 ), pows.at (0 ), ptbin);
132- if (hars.size () < 3 )
152+ if (hars.size () < ThreeParticleOrder )
133153 return TwoRec (hars.at (0 ), hars.at (1 ), pows.at (0 ), pows.at (1 ), ptbin, qpoi, qref, qol);
134154 int harlast = hars.at (hars.size () - 1 );
135155 int powlast = pows.at (pows.size () - 1 );
@@ -141,7 +161,7 @@ complex<double> GFW::RecursiveCorr(GFWCumulant* qpoi, GFWCumulant* qref, GFWCumu
141161 for (int i = harSize - 1 ; i >= 0 ; i--) {
142162 // checking if current configuration is a permutation of the next one.
143163 // Need to have more than 2 harmonics though, otherwise it doesn't make sense.
144- if (i > 2 ) { // only makes sense when we have more than two harmonics remaining
164+ if (i > static_cast < int >(TwoParticleOrder)) { // only makes sense when we have more than two harmonics remaining
145165 if (hars.at (i) == hars.at (i - 1 ) && pows.at (i) == pows.at (i - 1 )) { // if it is a permutation, then increase degeneracy and continue;
146166 lDegeneracy++;
147167 continue ;
@@ -181,30 +201,30 @@ GFW::CorrConfig GFW::GetCorrelatorConfig(string config, string head, bool ptdif)
181201 // Then make sure we don't have any double-spaces:
182202 while (s_index (config, " " ) > -1 )
183203 s_replace_all (config, " " , " " );
184- vector<int > regs;
185- vector<int > hars;
204+ std:: vector<int > regs;
205+ std:: vector<int > hars;
186206 int sz1 = 0 ;
187207 int szend = 0 ;
188208 string ts, ts2;
189209 CorrConfig ReturnConfig;
190210 // Fetch region descriptor
191211 if (!s_tokenize (config, ts, szend, " {" )) {
192- printf ( " Could not find any harmonics!\n " );
212+ LOGF (warning, " Could not find any harmonics!" );
193213 return ReturnConfig;
194214 }
195215 szend = 0 ;
196216 int counter = 0 ;
197217 while (s_tokenize (config, ts, szend, " {" )) {
198218 counter++;
199- ReturnConfig.Regs .push_back (vector<int >{});
200- ReturnConfig.Hars .push_back (vector<int >{});
219+ ReturnConfig.Regs .push_back (std:: vector<int >{});
220+ ReturnConfig.Hars .push_back (std:: vector<int >{});
201221 ReturnConfig.Overlap .push_back (-1 ); // initially, assume no overlap
202222 // Check if there's a particular pT bin I should be using here. If so, store it (otherwise, it's bin 0)
203223 int ptbin = -1 ;
204224 int sz2 = 0 ;
205225 if (s_contains (ts, " (" )) {
206226 if (!s_contains (ts, " )" )) {
207- printf ( " Missing \" )\" in the configurator. Returning...\n " );
227+ LOGF (warning, " Missing \" )\" in the configurator. Returning..." );
208228 return ReturnConfig;
209229 }
210230 sz2 = s_index (ts, " (" );
@@ -231,7 +251,7 @@ GFW::CorrConfig GFW::GetCorrelatorConfig(string config, string head, bool ptdif)
231251 if (ts2 == " " || ts2 == " " )
232252 continue ;
233253 if (ind < 0 ) {
234- printf ( " Could not find region named %s!\n " , ts2.c_str ());
254+ LOGF (warning, " Could not find region named %s!" , ts2.c_str ());
235255 break ;
236256 }
237257 if (!isOverlap)
@@ -253,7 +273,7 @@ GFW::CorrConfig GFW::GetCorrelatorConfig(string config, string head, bool ptdif)
253273 return ReturnConfig;
254274};
255275
256- complex <double > GFW::Calculate (int poi, int ref, vector<int > hars, int ptbin)
276+ complex <double > GFW::Calculate (int poi, int ref, std:: vector<int > hars, int ptbin)
257277{
258278 GFWCumulant* qref = &fCumulants .at (ref);
259279 GFWCumulant* qpoi = &fCumulants .at (poi);
@@ -305,13 +325,13 @@ complex<double> GFW::Calculate(CorrConfig corconf, int ptbin, bool SetHarmsToZer
305325 }
306326 return retval;
307327};
308- vector<pair<int , vector<int >>> GFW::GetHarmonicsSingleConfig (const CorrConfig& incfg)
328+ std:: vector<pair<int , std:: vector<int >>> GFW::GetHarmonicsSingleConfig (const CorrConfig& incfg)
309329{
310- vector<pair<int , vector<int >>> retPair;
330+ std:: vector<pair<int , std:: vector<int >>> retPair;
311331 for (int iR = 0 ; iR < static_cast <int >(incfg.Regs .size ()); iR++) {
312332 if (static_cast <int >(incfg.Regs [iR].size ()) > 1 ) {
313- retPair.push_back (make_pair (incfg.Regs [iR][0 ], vector<int >{incfg.Hars [iR][0 ]})); // If we have a PoI, then it comes with the first harmonic
314- retPair.push_back (make_pair (incfg.Regs [iR][1 ], incfg.Hars [iR])); // Then the second is ref. with full harmonics
333+ retPair.push_back (make_pair (incfg.Regs [iR][0 ], std:: vector<int >{incfg.Hars [iR][0 ]})); // If we have a PoI, then it comes with the first harmonic
334+ retPair.push_back (make_pair (incfg.Regs [iR][1 ], incfg.Hars [iR])); // Then the second is ref. with full harmonics
315335 } else {
316336 retPair.push_back (make_pair (incfg.Regs [iR][0 ], incfg.Hars [iR])); // Otherwise, it's only ref with all harmonics
317337 }
@@ -322,7 +342,7 @@ vector<pair<int, vector<int>>> GFW::GetHarmonicsSingleConfig(const CorrConfig& i
322342};
323343void GFW::InitializePowerArrays ()
324344{
325- vector<vector<vector<int >>> harSets (static_cast <int >(fRegions .size ()));
345+ std:: vector<std:: vector<std:: vector<int >>> harSets (static_cast <int >(fRegions .size ()));
326346 for (const CorrConfig& lConf : fListOfCFGs ) {
327347 auto HarPerReg = GetHarmonicsSingleConfig (lConf);
328348 for (const auto & oneHar : HarPerReg)
@@ -332,13 +352,13 @@ void GFW::InitializePowerArrays()
332352 for (int i = 0 ; i < static_cast <int >(harSets.size ()); i++) {
333353 if (fRegions [i].powsDefined )
334354 continue ; // Only do if powers have not been externally defined
335- vector<int > powerArray = GFWPowerArray::GetPowerArray (harSets[i]);
355+ std:: vector<int > powerArray = GFWPowerArray::GetPowerArray (harSets[i]);
336356 fRegions [i].Nhar = static_cast <int >(powerArray.size ());
337357 fRegions [i].NparVec = powerArray;
338358 fRegions [i].powsDefined = true ;
339359 }
340360};
341- complex <double > GFW::Calculate (int poi, vector<int > hars)
361+ complex <double > GFW::Calculate (int poi, std:: vector<int > hars)
342362{
343363 GFWCumulant* qpoi = &fCumulants .at (poi);
344364 return RecursiveCorr (qpoi, qpoi, qpoi, 0 , hars);
0 commit comments