99// granted to it by virtue of its status as an Intergovernmental Organization
1010// or submit itself to any jurisdiction.
1111
12+ // / \file FlowContainer.cxx
13+ // / \brief Container class to store and calculate multi-particle azimuthal correlations and cumulants
14+ // / \author Emil Gorm Dahlbæk Nielsen <emil.gorm.nielsen@cern.ch>
15+
1216#include " FlowContainer.h"
1317
1418#include " PWGCF/GenericFramework/Core/ProfileSubset.h"
3034#include < Rtypes.h>
3135#include < RtypesCore.h>
3236
33- #include < cstdio>
37+ #include < fstream>
38+ #include < string>
3439#include < vector>
3540
3641ClassImp (FlowContainer);
@@ -73,15 +78,15 @@ void FlowContainer::Initialize(TObjArray* inputList, const o2::framework::AxisSp
7378 if (nMultiBins <= 0 )
7479 nMultiBins = multiBins.size () - 1 ;
7580 if (nMultiBins <= 0 ) {
76- printf ( " Multiplicity axis does not exist" );
81+ LOGF (error, " Multiplicity axis does not exist" );
7782 return ;
7883 }
7984 if (!inputList) {
80- printf ( " Input list not specified\n " );
85+ LOGF (warning, " Input list not specified" );
8186 return ;
8287 }
8388 if (inputList->GetEntries () < 1 ) {
84- printf ( " Input list empty!\n " );
89+ LOGF (warning, " Input list empty!" );
8590 return ;
8691 }
8792 fProf = new TProfile2D (Form (" %s_CorrProfile" , this ->GetName ()), " CorrProfile" , nMultiBins, &multiBins[0 ], inputList->GetEntries (), 0.5 , inputList->GetEntries () + 0.5 );
@@ -101,11 +106,11 @@ void FlowContainer::Initialize(TObjArray* inputList, const o2::framework::AxisSp
101106void FlowContainer::Initialize (TObjArray* inputList, int nMultiBins, double MultiMin, double MultiMax, int nRandom)
102107{
103108 if (!inputList) {
104- printf ( " Input list not specified\n " );
109+ LOGF (warning, " Input list not specified" );
105110 return ;
106111 }
107112 if (inputList->GetEntries () < 1 ) {
108- printf ( " Input list empty!\n " );
113+ LOGF (warning, " Input list empty!" );
109114 return ;
110115 }
111116 fProf = new TProfile2D (Form (" %s_CorrProfile" , this ->GetName ()), " CorrProfile" , nMultiBins, MultiMin, MultiMax, inputList->GetEntries (), 0.5 , inputList->GetEntries () + 0.5 );
@@ -138,7 +143,7 @@ void FlowContainer::SetXAxis(TAxis* inax)
138143 fXAxis = dynamic_cast <TAxis*>(inax->Clone (" pTAxis" ));
139144 bool success = CreateBinsFromAxis (fXAxis );
140145 if (!success)
141- printf ( " Something went wrong setting the x axis!\n " );
146+ LOGF (warning, " Something went wrong setting the x axis!" );
142147}
143148void FlowContainer::SetXAxis ()
144149{
@@ -161,7 +166,7 @@ int FlowContainer::FillProfile(const char* hname, double multi, double corr, dou
161166 return -1 ;
162167 int yin = fProf ->GetYaxis ()->FindBin (hname);
163168 if (!yin) {
164- printf ( " Could not find bin %s\n " , hname);
169+ LOGF (info, " Could not find bin %s\n " , hname);
165170 return -1 ;
166171 }
167172 fProf ->Fill (multi, yin, corr, w);
@@ -176,23 +181,23 @@ void FlowContainer::OverrideProfileErrors(TProfile2D* inpf)
176181 int nBinsX = fProf ->GetNbinsX ();
177182 int nBinsY = fProf ->GetNbinsY ();
178183 if ((inpf->GetNbinsX () != nBinsX) || (inpf->GetNbinsY () != nBinsY)) {
179- printf ( " Number of bins in two profiles do not match, not doing anything\n " );
184+ LOGF (info, " Number of bins in two profiles do not match, not doing anything\n " );
180185 return ;
181186 }
182187 if (!inpf->GetBinSumw2 ()->fArray ) {
183- printf ( " Input profile has no BinSumw2()! Returning\n " );
188+ LOGF (info, " Input profile has no BinSumw2()! Returning\n " );
184189 return ;
185190 }
186191 if (!fProf ->GetBinSumw2 ()->fArray )
187192 fProf ->Sumw2 ();
188193 double * sumw2Prof = fProf ->GetSumw2 ()->fArray ;
189- double * sumw2Targ = inpf->GetSumw2 ()->fArray ;
194+ const double * sumw2Targ = inpf->GetSumw2 ()->fArray ;
190195 double * binsw2Prof = fProf ->GetBinSumw2 ()->fArray ;
191- double * binsw2Targ = inpf->GetBinSumw2 ()->fArray ;
196+ const double * binsw2Targ = inpf->GetBinSumw2 ()->fArray ;
192197 double * farrProf = fProf ->fArray ;
193198 for (int ix = 1 ; ix <= nBinsX; ix++) {
194199 double xval = fProf ->GetXaxis ()->GetBinCenter (ix);
195- printf ( " Processing x-bin %i\n " , ix);
200+ LOGF (info, " Processing x-bin %i\n " , ix);
196201 for (int iy = 1 ; iy <= nBinsY; iy++) {
197202 double yval = fProf ->GetYaxis ()->GetBinCenter (iy);
198203 int binno = fProf ->FindBin (xval, yval);
@@ -244,34 +249,38 @@ Long64_t FlowContainer::Merge(TCollection* collist)
244249
245250void FlowContainer::ReadAndMerge (const char * filelist)
246251{
247- FILE * flist = fopen (filelist, " r" );
248- char str[150 ];
249- int nFiles = 0 ;
250- while (fscanf (flist, " %s\n " , str) == 1 )
251- nFiles++;
252- rewind (flist);
253- if (nFiles == 0 ) {
254- printf (" No files to read!\n " );
252+ if (!filelist) {
253+ LOGF (error, " File list path is null!" );
254+ return ;
255+ }
256+ std::ifstream input (filelist);
257+ if (!input) {
258+ LOGF (error, " Could not open file list %s!" , filelist);
255259 return ;
256260 }
257- for (int i = 0 ; i < nFiles; i++) {
258- auto retVal = fscanf (flist, " %s\n " , str);
259- (void )retVal;
260- TFile* tf = new TFile (str, " READ" );
261- if (tf->IsZombie ()) {
262- printf (" Could not open file %s!\n " , str);
263- tf->Close ();
261+
262+ std::string filename;
263+ bool hasFiles = false ;
264+ while (input >> filename) {
265+ hasFiles = true ;
266+ TFile tf (filename.c_str (), " READ" );
267+ if (tf.IsZombie ()) {
268+ LOGF (info, " Could not open file %s!" , filename.c_str ());
264269 continue ;
265270 }
266- PickAndMerge (tf);
267- tf->Close ();
271+ PickAndMerge (&tf);
272+ }
273+ if (input.bad ()) {
274+ LOGF (error, " Error reading file list %s!" , filelist);
275+ } else if (!hasFiles) {
276+ LOGF (info, " No files to read!" );
268277 }
269278}
270279void FlowContainer::PickAndMerge (TFile* tfi)
271280{
272281 FlowContainer* lfc = dynamic_cast <FlowContainer*>(tfi->Get (this ->GetName ()));
273282 if (!lfc) {
274- printf ( " Could not pick up the %s from %s\n " , this ->GetName (), tfi->GetName ());
283+ LOGF (info, " Could not pick up the %s from %s" , this ->GetName (), tfi->GetName ());
275284 return ;
276285 }
277286 TProfile2D* spro = lfc->GetProfile ();
@@ -313,13 +322,13 @@ bool FlowContainer::OverrideBinsWithZero(int xb1, int yb1, int xb2, int yb2)
313322bool FlowContainer::OverrideMainWithSub (int ind, bool ExcludeChosen)
314323{
315324 if (!fProfRand ) {
316- printf ( " Cannot override main profile with a randomized one. Random profile array does not exist.\n " );
325+ LOGF (info, " Cannot override main profile with a randomized one. Random profile array does not exist." );
317326 return kFALSE ;
318327 }
319328 if (!ExcludeChosen) {
320329 TProfile2D* tarprof = dynamic_cast <TProfile2D*>(fProfRand ->At (ind));
321330 if (!tarprof) {
322- printf ( " Target random histogram does not exist.\n " );
331+ LOGF (info, " Target random histogram does not exist." );
323332 return kFALSE ;
324333 }
325334 TString ts (fProf ->GetName ());
@@ -345,7 +354,7 @@ bool FlowContainer::OverrideMainWithSub(int ind, bool ExcludeChosen)
345354bool FlowContainer::RandomizeProfile (int nSubsets)
346355{
347356 if (!fProfRand ) {
348- printf ( " Cannot randomize profile, random array does not exist.\n " );
357+ LOGF (info, " Cannot randomize profile, random array does not exist." );
349358 return kFALSE ;
350359 }
351360 int l_Subsets = nSubsets ? nSubsets : fProfRand ->GetEntries ();
@@ -393,7 +402,7 @@ TProfile* FlowContainer::GetCorrXXVsMulti(const char* order, int l_pti)
393402 const char * ybinlab = Form (" %s%s%s" , l_name.Data (), order, ptpf);
394403 int ybinno = fProf ->GetYaxis ()->FindBin (ybinlab);
395404 if (ybinno < 0 ) {
396- printf ( " Could not find %s!\n " , ybinlab);
405+ LOGF (info, " Could not find %s!" , ybinlab);
397406 return 0 ;
398407 }
399408 TProfile* rethist = dynamic_cast <TProfile*>(fProf ->ProfileX (" temp_prof" , ybinno, ybinno));
@@ -434,7 +443,7 @@ TH1D* FlowContainer::GetCorrXXVsPt(const char* order, double lminmulti, double l
434443 int ybn1 = fProf ->GetYaxis ()->FindBin (ybl1.Data ());
435444 int ybn2 = fProf ->GetYaxis ()->FindBin (ybl2.Data ());
436445 if (fNbinsPt != (ybn2 - ybn1 + 1 )) {
437- printf ( " fNbinsPt is not matching the num of found histograms" );
446+ LOGF (info, " fNbinsPt is not matching the num of found histograms" );
438447 return nullptr ;
439448 }
440449 const TString temporaryTag = Form (" %s_%s_%.3f_%.3f" , fIDName .Data (), order, lminmulti, lmaxmulti);
@@ -477,7 +486,7 @@ TH1D* FlowContainer::GetHistCorrXXVsPt(const char* order, double lminmulti, doub
477486{
478487 TH1D * rethist = GetCorrXXVsPt (order, lminmulti, lmaxmulti);
479488 if (!rethist) {
480- printf ( " GetCorrXXVsPt return nullptr!" );
489+ LOGF (info, " GetCorrXXVsPt return nullptr!" );
481490 return nullptr ;
482491 }
483492 TProfile* refflow = GetRefFlowProfile (order, lminmulti, lmaxmulti);
@@ -888,23 +897,24 @@ TH1D* FlowContainer::GetVN8VsX(int n, bool onPt, double arg1, double arg2)
888897 }
889898 return rethist;
890899}
900+
891901TH1D * FlowContainer::GetCNN (int n, int c, bool onPt, double arg1, double arg2)
892902{
893- if (c == 8 )
903+ if (c == kEightParticleOrder )
894904 return GetCN8VsX (n, onPt, arg1, arg2);
895- if (c == 6 )
905+ if (c == kSixParticleOrder )
896906 return GetCN6VsX (n, onPt, arg1, arg2);
897- if (c == 4 )
907+ if (c == kFourParticleOrder )
898908 return GetCN4VsX (n, onPt, arg1, arg2);
899909 return GetCN2VsX (n, onPt, arg1, arg2);
900910};
901911TH1D * FlowContainer::GetVNN (int n, int c, bool onPt, double arg1, double arg2)
902912{
903- if (c == 8 )
913+ if (c == kEightParticleOrder )
904914 return GetVN8VsX (n, onPt, arg1, arg2);
905- if (c == 6 )
915+ if (c == kSixParticleOrder )
906916 return GetVN6VsX (n, onPt, arg1, arg2);
907- if (c == 4 )
917+ if (c == kFourParticleOrder )
908918 return GetVN4VsX (n, onPt, arg1, arg2);
909919 return GetVN2VsX (n, onPt, arg1, arg2);
910920};
@@ -1049,12 +1059,14 @@ double FlowContainer::CN6Error(double cor6e, double cor4, double cor4e, double c
10491059{
10501060 if (!fPropagateErrors )
10511061 return 0 ;
1052- double inters[3 ];
1062+
1063+ constexpr int kCN6Terms = 3 ;
1064+ double inters[kCN6Terms ];
10531065 inters[0 ] = cor6e;
10541066 inters[1 ] = -9 * cor2 * cor4e;
10551067 inters[2 ] = (-9 * cor4 + 36 * cor2 * cor2) * cor2e;
10561068 double sum = 0 ;
1057- for (int i = 0 ; i < 3 ; i++)
1069+ for (int i = 0 ; i < kCN6Terms ; i++)
10581070 sum += (inters[i] * inters[i]);
10591071 return TMath::Sqrt (sum);
10601072};
@@ -1068,14 +1080,15 @@ double FlowContainer::DN6Error(double d6e, double d4, double d4e, double d2,
10681080{
10691081 if (!fPropagateErrors )
10701082 return 0 ;
1071- double inters[5 ];
1083+ constexpr int kDN6Terms = 5 ;
1084+ double inters[kDN6Terms ];
10721085 inters[0 ] = d6e;
10731086 inters[1 ] = -6 * c2 * d4e;
10741087 inters[2 ] = (-3 * c4 + 12 * c2 * c2) * d2e;
10751088 inters[3 ] = -3 * d2 * c4e;
10761089 inters[4 ] = (-6 * d4 + 24 * d2 * c2) * c2e;
10771090 double sum = 0 ;
1078- for (int i = 0 ; i < 5 ; i++)
1091+ for (int i = 0 ; i < kDN6Terms ; i++)
10791092 sum += (inters[i] * inters[i]);
10801093 return TMath::Sqrt (sum);
10811094};
@@ -1124,13 +1137,14 @@ double FlowContainer::CN8Error(double cor8e, double cor6, double cor6e,
11241137{
11251138 if (!fPropagateErrors )
11261139 return 0 ;
1127- double parts[4 ];
1140+ constexpr int kCN8Terms = 4 ;
1141+ double parts[kCN8Terms ];
11281142 parts[0 ] = cor8e;
11291143 parts[1 ] = -16 * cor2 * cor6e;
11301144 parts[2 ] = (-36 * cor4 + 144 * cor2 * cor2) * cor4e;
11311145 parts[3 ] = (-16 * cor6 + 288 * cor4 * cor2 + 576 * cor2 * cor2 * cor2) * cor2e;
11321146 double retval = 0 ;
1133- for (int i = 0 ; i < 4 ; i++)
1147+ for (int i = 0 ; i < kCN8Terms ; i++)
11341148 retval += TMath::Power (parts[i], 2 );
11351149 return TMath::Sqrt (retval);
11361150};
@@ -1145,7 +1159,8 @@ double FlowContainer::DN8Error(double d8e, double d6, double d6e, double d4,
11451159{
11461160 if (!fPropagateErrors )
11471161 return 0 ;
1148- double parts[7 ];
1162+ constexpr int kDN8Terms = 7 ;
1163+ double parts[kDN8Terms ];
11491164 parts[0 ] = d8e; // d/d8'
11501165 parts[1 ] = -12 * c2 * d6e; // d/d6'
11511166 parts[2 ] = -4 * d2 * c6e; // d/d6
@@ -1154,7 +1169,7 @@ double FlowContainer::DN8Error(double d8e, double d6, double d6e, double d4,
11541169 parts[5 ] = (-4 * c6 + 72 * c4 * c2 - 144 * c2 * c2 * c2) * d2e;
11551170 parts[6 ] = (-12 * d6 + 144 * d4 * c2 + 72 * c4 * d2 - 432 * d2 * c2 * c2) * c2e;
11561171 double retval = 0 ;
1157- for (int i = 0 ; i < 7 ; i++)
1172+ for (int i = 0 ; i < kDN8Terms ; i++)
11581173 retval += TMath::Power (parts[i], 2 );
11591174 return TMath::Sqrt (retval);
11601175};
@@ -1195,26 +1210,6 @@ void FlowContainer::SetPtRebin(int nbins, double* binedges)
11951210{
11961211 fPtRebin = nbins;
11971212 fPtRebinEdges = binedges;
1198- return ;
1199- int fPtRebin = 0 ;
1200- // double *lPtRebinEdges=binedges;
1201- if (!fbinsPt)
1202- SetXAxis ();
1203- for (int i = 0 ; i < nbins; i++)
1204- if (binedges[i] < fbinsPt[0 ] || binedges[i] > fbinsPt[fNbinsPt - 1 ])
1205- continue ;
1206- else
1207- fPtRebin ++;
1208- if (fPtRebinEdges )
1209- delete[] fPtRebinEdges ;
1210- fPtRebinEdges = new double [fPtRebin ];
1211- fPtRebin = 0 ;
1212- for (int i = 0 ; i < nbins; i++)
1213- if (binedges[i] < fbinsPt[0 ] || binedges[i] > fbinsPt[fNbinsPt ])
1214- continue ;
1215- else
1216- fPtRebinEdges [fPtRebin ++] = binedges[i];
1217- // fPtRebin--;
12181213}
12191214void FlowContainer::SetMultiRebin (int nbins, double * binedges)
12201215{
0 commit comments