869164495a37e52ef7a9130ce3596f9b90ad829a
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowAnalysisWithQCumulants.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  * 
14 **************************************************************************/
15
16 /********************************** 
17  * flow analysis with Q-cumulants * 
18  *                                * 
19  * author:  Ante Bilandzic        * 
20  *           (anteb@nikhef.nl)    *
21  *********************************/ 
22
23 #define AliFlowAnalysisWithQCumulants_cxx
24
25 #include "Riostream.h"
26 #include "AliFlowCommonConstants.h"
27 #include "AliFlowCommonHist.h"
28 #include "AliFlowCommonHistResults.h"
29 #include "TChain.h"
30 #include "TFile.h"
31 #include "TList.h"
32 #include "TGraph.h"
33 #include "TParticle.h"
34 #include "TRandom3.h"
35 #include "TStyle.h"
36 #include "TProfile.h"
37 #include "TProfile2D.h" 
38 #include "TProfile3D.h"
39 #include "TMath.h"
40 #include "TArrow.h"
41 #include "TPaveLabel.h"
42 #include "TCanvas.h"
43 #include "AliFlowEventSimple.h"
44 #include "AliFlowTrackSimple.h"
45 #include "AliFlowAnalysisWithQCumulants.h"
46 #include "TArrayD.h"
47 #include "TRandom.h"
48 #include "TF1.h"
49
50 class TH1;
51 class TH2;
52 class TGraph;
53 class TPave;
54 class TLatex;
55 class TMarker;
56 class TRandom3;
57 class TObjArray;
58 class TList;
59 class TCanvas;
60 class TSystem;
61 class TROOT;
62 class AliFlowVector;
63 class TVector;
64
65
66 //================================================================================================================
67
68
69 ClassImp(AliFlowAnalysisWithQCumulants)
70
71 AliFlowAnalysisWithQCumulants::AliFlowAnalysisWithQCumulants(): 
72  // 0.) base:
73  fHistList(NULL),
74  // 1.) common:
75  fCommonHists(NULL),
76  fCommonHists2nd(NULL), 
77  fCommonHists4th(NULL),
78  fCommonHists6th(NULL),
79  fCommonHists8th(NULL),
80  fCommonHistsResults2nd(NULL),
81  fCommonHistsResults4th(NULL),
82  fCommonHistsResults6th(NULL),
83  fCommonHistsResults8th(NULL),
84  fnBinsPhi(0),
85  fPhiMin(0),
86  fPhiMax(0),
87  fPhiBinWidth(0),
88  fnBinsPt(0),
89  fPtMin(0),
90  fPtMax(0),
91  fPtBinWidth(0),
92  fnBinsEta(0),
93  fEtaMin(0),
94  fEtaMax(0),
95  fEtaBinWidth(0),
96  fHarmonic(2),
97  fAnalysisLabel(NULL),
98  // 2.) weights:
99  fWeightsList(NULL),
100  fUsePhiWeights(kFALSE),
101  fUsePtWeights(kFALSE),
102  fUseEtaWeights(kFALSE),
103  fUseParticleWeights(NULL),
104  fPhiWeights(NULL),
105  fPtWeights(NULL),
106  fEtaWeights(NULL),
107  // 3.) integrated flow:
108  fIntFlowList(NULL), 
109  fIntFlowProfiles(NULL),
110  fIntFlowResults(NULL),
111  fReQ(NULL),
112  fImQ(NULL),
113  fSMpk(NULL),
114  fAvMultiplicity(NULL),
115  // 4.) differential flow:
116  fDiffFlowList(NULL),
117  fDiffFlowProfiles(NULL),
118  fDiffFlowResults(NULL),
119  // 5.) distributions:
120  fDistributionsList(NULL),
121  // x.) debugging and cross-checking:
122  fNestedLoopsList(NULL),
123  fEvaluateNestedLoopsForIntFlow(kFALSE),
124  fEvaluateNestedLoopsForDiffFlow(kFALSE),  
125  fEvaluateNestedLoops(NULL),
126  fDirectCorrelations(NULL),
127  fDirectCorrectionsCos(NULL),
128  fDirectCorrectionsSin(NULL),
129  fDirectCorrelationsDiffFlow(NULL),
130  fDirectCorrectionsDiffFlowCos(NULL),
131  fDirectCorrectionsDiffFlowSin(NULL),
132  fDirectCorrelationsW(NULL),
133  fDirectCorrectionsCosW(NULL),
134  fDirectCorrectionsSinW(NULL),
135  fDirectCorrelationsDiffFlowW(NULL),
136  fDirectCorrectionsDiffFlowCosW(NULL),
137  fDirectCorrectionsDiffFlowSinW(NULL)
138  {
139   // constructor  
140   
141   // base list to hold all output objects:
142   fHistList = new TList();
143   fHistList->SetName("cobjQC");
144   fHistList->SetOwner(kTRUE);
145   
146   // list to hold histograms with phi, pt and eta weights:      
147   fWeightsList = new TList();
148     
149   // analysis label;
150   fAnalysisLabel = new TString();
151       
152   // initialize all arrays:  
153   this->InitializeArraysForIntFlow();
154   this->InitializeArraysForDiffFlow();
155   this->InitializeArraysForDistributions();
156   
157  } // end of constructor
158  
159
160 //================================================================================================================  
161
162
163 AliFlowAnalysisWithQCumulants::~AliFlowAnalysisWithQCumulants()
164 {
165  // destructor
166  
167  delete fHistList;
168
169 } // end of AliFlowAnalysisWithQCumulants::~AliFlowAnalysisWithQCumulants()
170
171
172 //================================================================================================================
173
174
175 void AliFlowAnalysisWithQCumulants::Init()
176 {
177  // initialize all constants and book everything
178  
179  // access constants:
180  this->AccessConstants();
181  
182  // booking:
183  this->BookAndFillWeightsHistograms();
184  this->BookAndNestAllLists();
185  this->BookCommonHistograms();
186  this->BookEverythingForIntegratedFlow(); 
187  this->BookEverythingForDifferentialFlow(); 
188  this->BookEverythingForDistributions();
189  this->BookEverythingForNestedLoops();
190  
191  // set harmonic in common control histograms (to be improved (should I do this somewhere else?)):
192  (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic);
193  (fCommonHists2nd->GetHarmonic())->Fill(0.5,fHarmonic);
194  (fCommonHists4th->GetHarmonic())->Fill(0.5,fHarmonic);
195  (fCommonHists6th->GetHarmonic())->Fill(0.5,fHarmonic);
196  (fCommonHists8th->GetHarmonic())->Fill(0.5,fHarmonic);
197  
198 } // end of void AliFlowAnalysisWithQCumulants::Init()
199
200
201 //================================================================================================================
202
203
204 void AliFlowAnalysisWithQCumulants::Make(AliFlowEventSimple* anEvent)
205 {
206  // running over data only in this method
207  
208  // a) fill the common control histograms and call method to fill fAvMultiplicity;
209  // b) loop over data to calculate e-b-e quantities;
210  // c) call the methods;
211  // d) debugging and cross-checking (evaluate nested loops);
212  // e) reset e-b-e quantities.
213  
214  Double_t dPhi = 0.; // azimuthal angle in the laboratory frame
215  Double_t dPt  = 0.; // transverse momentum
216  Double_t dEta = 0.; // pseudorapidity
217
218  Double_t wPhi = 1.; // phi weight
219  Double_t wPt  = 1.; // pt weight
220  Double_t wEta = 1.; // eta weight
221                                                                                                                                 
222  // ********************************************
223  // **** FILL THE COMMON CONTROL HISTOGRAMS ****
224  // ********************************************
225                                          
226  Int_t nRP = anEvent->GetEventNSelTracksRP(); // number of RPs (i.e. number of particles used to determine the reaction plane)
227
228  fCommonHists->FillControlHistograms(anEvent); 
229  
230  if(nRP>1)
231  {
232   fCommonHists2nd->FillControlHistograms(anEvent);                                        
233   if(nRP>3)
234   {
235    fCommonHists4th->FillControlHistograms(anEvent);                                        
236    if(nRP>5)
237    {
238     fCommonHists6th->FillControlHistograms(anEvent);                                        
239     if(nRP>7)
240     {
241      fCommonHists8th->FillControlHistograms(anEvent);                                        
242     } // end of if(nRP>7)  
243    } // end of if(nRP>5) 
244   } // end of if(nRP>3)                                                                                                                      
245  } // end of if(nRP>1) 
246                              
247  this->FillAverageMultiplicities(nRP);                                                                  
248                                                                                                                                             
249  // *******************************************************
250  // **** LOOP OVER DATA AND CALCULATE E-B-E QUANTITIES ****
251  // *******************************************************
252  
253  Int_t nPrim = anEvent->NumberOfTracks();  // nPrim = total number of primary tracks, i.e. nPrim = nRP + nPOI + rest, where:
254                                            // nRP   = # of particles used to determine the reaction plane;
255                                            // nPOI  = # of particles of interest for a detailed flow analysis;
256                                            // rest  = # of particles which are not niether RPs nor POIs.  
257  
258  AliFlowTrackSimple *aftsTrack = NULL;
259  
260  for(Int_t i=0;i<nPrim;i++) 
261  { 
262   aftsTrack=anEvent->GetTrack(i);
263   if(aftsTrack)
264   {
265    if(!(aftsTrack->InRPSelection() || aftsTrack->InPOISelection())) continue; // consider only tracks which are RPs or POIs
266    
267    Int_t n = fHarmonic; // shortcut for the harmonic
268  
269    if(aftsTrack->InRPSelection()) // RP condition:
270    {    
271     dPhi = aftsTrack->Phi();
272     dPt  = aftsTrack->Pt();
273     dEta = aftsTrack->Eta();
274   
275     if(fUsePhiWeights && fPhiWeights && fnBinsPhi) // determine phi weight for this particle:
276     {
277      wPhi = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*fnBinsPhi/TMath::TwoPi())));
278     }
279     if(fUsePtWeights && fPtWeights && fnBinsPt) // determine pt weight for this particle:
280     {
281      wPt = fPtWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-fPtMin)/fPtBinWidth))); 
282     }            
283     if(fUseEtaWeights && fEtaWeights && fEtaBinWidth) // determine eta weight for this particle: 
284     {
285      wEta = fEtaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-fEtaMin)/fEtaBinWidth))); 
286     } 
287     
288     // integrated flow: 
289     // calculate Re[Q_{m*n,k}] and Im[Q_{m*n,k}], m = 1,2,3,4, for this event:
290     for(Int_t m=0;m<4;m++)
291     {
292      for(Int_t k=0;k<9;k++)
293      {
294       (*fReQ)(m,k)+=pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1)*n*dPhi); 
295       (*fImQ)(m,k)+=pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1)*n*dPhi); 
296      } 
297     }
298     // calculate S^{M}_{p,k} for this event 
299     // Remark: final calculation of S^{M}_{p,k} follows after the loop over data bellow:
300     for(Int_t p=0;p<8;p++)
301     {
302      for(Int_t k=0;k<9;k++)
303      {     
304       (*fSMpk)(p,k)+=pow(wPhi*wPt*wEta,k);
305      }
306     } 
307
308     // differential flow:
309     // (r_{m*m,k}(pt,eta)): 
310     for(Int_t m=0;m<4;m++)
311     {
312      for(Int_t k=0;k<9;k++)
313      {
314       fReEBE[0][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1.)*n*dPhi),1.);
315       fImEBE[0][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1.)*n*dPhi),1.);
316      }
317     } 
318     // s_{k}(pt,eta) for RPs // to be improved (clarified)
319     // Remark: final calculation of s_{p,k}(pt,eta) follows after the loop over data bellow:
320     for(Int_t k=0;k<9;k++)
321     {
322      fs[0][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k),1.);
323     } 
324      
325     if(aftsTrack->InPOISelection())
326     {
327      // (q_{m*m,k}(pt,eta)): 
328      for(Int_t m=0;m<4;m++)
329      {
330       for(Int_t k=0;k<9;k++)
331       {
332        fReEBE[2][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k)*TMath::Cos((m+1.)*n*dPhi),1.);
333        fImEBE[2][m][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k)*TMath::Sin((m+1.)*n*dPhi),1.);
334       }
335      } 
336      // s_{k}(pt,eta) for RP&&POIs // to be improved (clarified)
337      // Remark: final calculation of s_{p,k}(pt,eta) follows after the loop over data bellow:
338      for(Int_t k=0;k<9;k++)
339      {
340       fs[2][k]->Fill(dPt,dEta,pow(wPhi*wPt*wEta,k),1.);
341      }
342       
343     } // end of if(aftsTrack->InPOISelection())
344    } // end of if(pTrack->InRPSelection())
345
346    if(aftsTrack->InPOISelection())
347    {
348     dPhi = aftsTrack->Phi();
349     dPt  = aftsTrack->Pt();
350     dEta = aftsTrack->Eta();
351        
352     // p_n(m*n,0):   
353     for(Int_t m=0;m<4;m++)
354     {
355      fReEBE[1][m][0]->Fill(dPt,dEta,TMath::Cos((m+1.)*n*dPhi),1.);
356      fImEBE[1][m][0]->Fill(dPt,dEta,TMath::Sin((m+1.)*n*dPhi),1.);
357     } 
358     
359    } // end of if(pTrack->InPOISelection() )   
360  
361   } else // to if(aftsTrack)
362     {
363      cout<<endl;
364      cout<<" WARNING: no particle! (i.e. aftsTrack is a NULL pointer in AFAWQC::Make().)"<<endl;
365      cout<<endl;       
366     }
367  } // end of for(Int_t i=0;i<nPrim;i++) 
368
369  // calculate the final expressions for S^{M}_{p,k}:
370  for(Int_t p=0;p<8;p++)
371  {
372   for(Int_t k=0;k<9;k++)
373   {
374    (*fSMpk)(p,k)=pow((*fSMpk)(p,k),p+1);
375   }  
376  } 
377  
378  // *****************************
379  // **** CALL THE METHODS *******
380  // *****************************
381
382  if(!fEvaluateNestedLoopsForIntFlow)
383  {
384   // without weights:
385   if(nRP>1) this->CalculateCorrelationsForIntegratedFlow();
386   if(nRP>0) this->CalculateCorrectionsForNonUniformAcceptanceForIntFlowCosTerms();
387   if(nRP>0) this->CalculateCorrectionsForNonUniformAcceptanceForIntFlowSinTerms();
388   if(nRP>3) this->CalculateQProductsForIntFlow();
389   if(nRP>1) this->CalculateSumAndProductOfEventWeights();
390   // with weights:
391   if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
392   {
393    if(nRP>1) this->CalculateWeightedCorrelationsForIntegratedFlow();
394    if(nRP>3) this->CalculateWeightedQProductsForIntFlow();
395   } 
396  }
397
398  if(!fEvaluateNestedLoopsForDiffFlow)
399  {
400   // without weights:
401   if(nRP>1) this->CalculateCorrelationsForDifferentialFlow("RP");
402   if(nRP>1) this->CalculateCorrelationsForDifferentialFlow("POI");
403   
404   // with weights:
405   if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
406   {
407    if(nRP>1) this->CalculateWeightedCorrelationsForDifferentialFlow("RP");
408    if(nRP>1) this->CalculateWeightedCorrelationsForDifferentialFlow("POI");
409   } 
410  }
411  
412  // **************************************************************
413  // **** DEBUGGING AND CROSS-CHECKING (EVALUATE NESTED LOOPS) ****
414  // **************************************************************
415
416  if(fEvaluateNestedLoopsForIntFlow)
417  {
418   if(nPrim>0 && nPrim<15) // only for these multiplicities it is feasible to evaluate 8 nested loops in short time 
419   {
420    // without weights:
421    if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights))
422    {
423     this->CalculateCorrelationsForIntegratedFlow();
424     this->CalculateCorrectionsForNonUniformAcceptanceForIntFlowCosTerms();
425     this->CalculateCorrectionsForNonUniformAcceptanceForIntFlowSinTerms();
426    }
427    // with weights:
428    if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
429    {
430     this->CalculateWeightedCorrelationsForIntegratedFlow();
431    }
432     
433    this->EvaluateNestedLoopsForIntegratedFlow(anEvent);  
434   }
435  } 
436  
437  if(fEvaluateNestedLoopsForDiffFlow)
438  {
439   if(nPrim>0 && nPrim<15) // only for these multiplicities it is feasible to evaluate 8 nested loops in short time 
440   {
441    // without weights:
442    if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights))
443    {
444     this->CalculateCorrelationsForDifferentialFlow("RP");
445     this->CalculateCorrelationsForDifferentialFlow("POI");
446    }
447    // with weights:
448    if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
449    {
450     this->CalculateWeightedCorrelationsForDifferentialFlow("RP");
451     this->CalculateWeightedCorrelationsForDifferentialFlow("POI");
452    }
453     
454    this->EvaluateNestedLoopsForDifferentialFlow(anEvent);  
455   }
456  } 
457  
458  // ********************************
459  // **** RESET E-B-E QUANTITIES ****
460  // ********************************
461  
462  // integrated flow:
463  fReQ->Zero();
464  fImQ->Zero();
465  fSMpk->Zero();
466  for(Int_t pW=0;pW<1+(Int_t)(fUsePhiWeights||fUsePtWeights||fUseEtaWeights);pW++)
467  {
468   fQCorrelationsEBE[pW]->Reset();
469   for(Int_t sc=0;sc<2;sc++)
470   {
471     fQCorrectionsEBE[pW][sc]->Reset();
472   } 
473  }  
474  
475  // differential flow:
476  for(Int_t t=0;t<3;t++) // type (RP, POI, POI&&RP)
477  {
478   for(Int_t m=0;m<4;m++) // multiple of harmonic
479   {
480    for(Int_t k=0;k<9;k++) // power of weight
481    {
482     if(fReEBE[t][m][k]) fReEBE[t][m][k]->Reset();
483     if(fImEBE[t][m][k]) fImEBE[t][m][k]->Reset();
484    }   
485   }
486  }
487  
488  for(Int_t t=0;t<3;t++) // type (0 = RP, 1 = POI, 2 = RP&&POI )
489  { 
490   for(Int_t k=0;k<9;k++)
491   {
492    if(fs[t][k]) fs[t][k]->Reset();
493   }
494  }
495  
496 } // end of AliFlowAnalysisWithQCumulants::Make(AliFlowEventSimple* anEvent)
497
498
499 //================================================================================================================================
500
501
502 void AliFlowAnalysisWithQCumulants::Finish()
503 {
504  // calculate the final results
505
506  // a) acces the constants;
507  // b) access the flags;
508  // c) calculate the final results for integrated flow (without and with weights);
509  // d) store in AliFlowCommonHistResults and print the final results for integrated flow;
510  // e) calculate the final results for differential flow (without and with weights);
511  // f) print the final results for integrated flow obtained from differential flow (to be improved (terminology));
512  // g) COMPARE RESULTS FROM NESTED LOOPS vs RESULTS FROM Q-VECTORS FOR INTEGRATED FLOW
513  
514  // ******************************
515  // **** ACCESS THE CONSTANTS ****
516  // ******************************
517  
518  this->AccessConstants();          
519  
520  // **************************
521  // **** ACCESS THE FLAGS ****
522  // **************************    
523
524  fUsePhiWeights = (Int_t)fUseParticleWeights->GetBinContent(1); 
525  fUsePtWeights = (Int_t)fUseParticleWeights->GetBinContent(2); 
526  fUseEtaWeights = (Int_t)fUseParticleWeights->GetBinContent(3); 
527  fEvaluateNestedLoopsForIntFlow = (Int_t)fEvaluateNestedLoops->GetBinContent(1);
528  fEvaluateNestedLoopsForDiffFlow = (Int_t)fEvaluateNestedLoops->GetBinContent(2); 
529     
530  // *********************************************************
531  // **** CALCULATE THE FINAL RESULTS FOR INTEGRATED FLOW ****
532  // *********************************************************    
533  
534  // without weights:
535  this->FinalizeCorrelationsForIntFlow(kFALSE,"exact");
536  this->CalculateFinalCorrectionsForNonUniformAcceptanceForCumulantsForIntFlow(kFALSE,"exact");
537  this->CalculateCovariancesForIntFlow(kFALSE,"exact");
538  this->CalculateCumulantsForIntFlow(kFALSE,"exact");
539  this->ApplyCorrectionForNonUniformAcceptanceToCumulantsForIntFlow(kFALSE,"exact");
540  this->CalculateIntFlow(kFALSE,"exact",kFALSE); // pW = 0, eW = 0, not corrected for non-uniform acceptance
541  this->CalculateIntFlow(kFALSE,"exact",kTRUE); // pW = 0, eW = 0, corrected for non-uniform acceptance
542  
543  // with weights:
544  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
545  {
546   this->FinalizeCorrelationsForIntFlow(kTRUE,"exact"); 
547   // this->CalculateFinalCorrectionsForNonUniformAcceptanceForCumulantsForIntFlow(kTRUE,"exact");
548   this->CalculateCovariancesForIntFlow(kTRUE,"exact");
549   this->CalculateCumulantsForIntFlow(kTRUE,"exact");
550   // this->ApplyCorrectionForNonUniformAcceptanceToCumulantsForIntFlow(kTRUE,"exact");
551   this->CalculateIntFlow(kTRUE,"exact",kFALSE); // weighted and not corrected for non-uniform acceptance
552   // this->CalculateIntFlow(kTRUE,"exact",kTRUE); // weighted and corrected for non-uniform acceptance
553  }
554  
555  // ***************************************************************
556  // **** STORE AND PRINT THE FINAL RESULTS FOR INTEGRATED FLOW ****
557  // ***************************************************************
558  
559  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
560  {        
561   this->FillCommonHistResultsIntFlow(kTRUE,"exact",kFALSE); // weighted and not corrected for non-uniform acceptance       
562   // this->FillCommonHistResultsIntFlow(kTRUE,kTRUE); // weighted and corrected for non-uniform acceptance (to be improved (enabled))    
563   // this->PrintQuantifyingCorrectionsForNonUniformAcceptance(kTRUE,"exact"); // (to be improved (enabled))
564  } else 
565    {
566     this->FillCommonHistResultsIntFlow(kFALSE,"exact",kTRUE); // non-weighted and corrected for non-uniform acceptance 
567     this->PrintQuantifyingCorrectionsForNonUniformAcceptance(kFALSE,"exact"); 
568    }
569  
570  this->PrintFinalResultsForIntegratedFlow("NONAME"); // to be improved (name)
571  
572  // ***********************************************************
573  // **** CALCULATE THE FINAL RESULTS FOR DIFFERENTIAL FLOW ****
574  // ***********************************************************    
575  
576  // without weights:
577  this->FinalizeCorrelationsForDiffFlow("RP",kFALSE,"exact");
578  this->FinalizeCorrelationsForDiffFlow("POI",kFALSE,"exact");
579  this->CalculateCumulantsForDiffFlow("RP",kFALSE,"exact");
580  this->CalculateCumulantsForDiffFlow("POI",kFALSE,"exact");
581  this->CalculateDiffFlow("RP",kFALSE,"exact");
582  this->CalculateDiffFlow("POI",kFALSE,"exact");
583  this->CalculateFinalResultsForRPandPOIIntegratedFlow("RP",kFALSE,"exact");
584  this->CalculateFinalResultsForRPandPOIIntegratedFlow("POI",kFALSE,"exact");
585  
586  // with weights:
587  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
588  {
589   this->FinalizeCorrelationsForDiffFlow("RP",kTRUE,"exact");
590   this->FinalizeCorrelationsForDiffFlow("POI",kTRUE,"exact");
591   this->CalculateCumulantsForDiffFlow("RP",kTRUE,"exact");
592   this->CalculateCumulantsForDiffFlow("POI",kTRUE,"exact");
593   this->CalculateDiffFlow("RP",kTRUE,"exact");
594   this->CalculateDiffFlow("POI",kTRUE,"exact");
595   this->CalculateFinalResultsForRPandPOIIntegratedFlow("RP",kTRUE,"exact");
596   this->CalculateFinalResultsForRPandPOIIntegratedFlow("POI",kTRUE,"exact");
597  }
598  
599    
600   //this->CalculateFinalCorrectionsForNonUniformAcceptanceForDifferentialFlow(kFALSE,"POI"); // to be improved (to calculate also when weights are used) 
601   //this->CalculateFinalCorrectionsForNonUniformAcceptanceForDifferentialFlow(kFALSE,"RP"); // to be improved (to calculate also when weights are used)
602
603  
604  // *****************************************************************
605  // **** STORE AND PRINT THE FINAL RESULTS FOR DIFFERENTIAL FLOW ****
606  // *****************************************************************
607  if(fUsePhiWeights||fUsePtWeights||fUseEtaWeights)
608  {
609   this->FillCommonHistResultsDiffFlow("RP",kTRUE,"exact",kFALSE);
610   this->FillCommonHistResultsDiffFlow("POI",kTRUE,"exact",kFALSE);
611  } else
612    {
613     this->FillCommonHistResultsDiffFlow("RP",kFALSE,"exact",kFALSE);
614     this->FillCommonHistResultsDiffFlow("POI",kFALSE,"exact",kFALSE);
615    }
616  
617  this->PrintFinalResultsForIntegratedFlow("RP"); 
618  this->PrintFinalResultsForIntegratedFlow("POI"); 
619   
620  // *****************************************************************************************
621  // **** COMPARE RESULTS FROM NESTED LOOPS vs RESULTS FROM Q-VECTORS FOR INTEGRATED FLOW ****
622  // *****************************************************************************************    
623  
624  if(fEvaluateNestedLoopsForIntFlow) 
625  {
626   this->CompareResultsFromNestedLoopsAndFromQVectorsForIntFlow(fUsePhiWeights||fUsePtWeights||fUseEtaWeights);
627  } 
628  
629  if(fEvaluateNestedLoopsForDiffFlow) 
630  {
631   this->CompareResultsFromNestedLoopsAndFromQVectorsForDiffFlow(fUsePhiWeights||fUsePtWeights||fUseEtaWeights);
632  } 
633                                                                                                                                                                                                                                                                                                                                    
634 } // end of AliFlowAnalysisWithQCumulants::Finish()
635
636
637 //================================================================================================================================
638
639
640 void AliFlowAnalysisWithQCumulants::CalculateCorrectionsForNonUniformAcceptanceForIntFlowCosTerms()
641 {
642  // calculate corrections for non-uniform acceptance of the detector for no-name integrated flow (cos terms)
643  
644  // multiplicity:
645  Double_t dMult = (*fSMpk)(0,0);
646  
647  // real and imaginary parts of non-weighted Q-vectors evaluated in harmonics n, 2n, 3n and 4n: 
648  Double_t dReQ1n = (*fReQ)(0,0);
649  Double_t dReQ2n = (*fReQ)(1,0);
650  //Double_t dReQ3n = (*fReQ)(2,0);
651  //Double_t dReQ4n = (*fReQ)(3,0);
652  Double_t dImQ1n = (*fImQ)(0,0);
653  Double_t dImQ2n = (*fImQ)(1,0);
654  //Double_t dImQ3n = (*fImQ)(2,0);
655  //Double_t dImQ4n = (*fImQ)(3,0);
656         
657  //                                  *************************************************************
658  //                                  **** corrections for non-uniform acceptance (cos terms): ****
659  //                                  *************************************************************
660  //
661  // Remark 1: corrections for non-uniform acceptance (cos terms) calculated with non-weighted Q-vectors 
662  //           are stored in 1D profile fQCorrectionsCos.
663  // Remark 2: binning of fQCorrectionsCos is organized as follows:
664  // --------------------------------------------------------------------------------------------------------------------
665  // 1st bin: <<cos(n*(phi1))>> = cosP1n
666  // 2nd bin: <<cos(n*(phi1+phi2))>> = cosP1nP1n
667  // 3rd bin: <<cos(n*(phi1-phi2-phi3))>> = cosP1nM1nM1n
668  // ...
669  // --------------------------------------------------------------------------------------------------------------------
670   
671  // 1-particle:
672  Double_t cosP1n = 0.; // <<cos(n*(phi1))>>
673    
674  if(dMult>0)
675  {
676   cosP1n = dReQ1n/dMult; 
677   
678   // average non-weighted 1-particle correction (cos terms) for non-uniform acceptance for single event:
679   fQCorrectionsEBE[0][1]->SetBinContent(1,cosP1n);
680   
681   // final average non-weighted 1-particle correction (cos terms) for non-uniform acceptance for all events:
682   fQCorrections[0][0][1]->Fill(0.5,cosP1n,dMult);  
683  } 
684  
685  // 2-particle:
686  Double_t cosP1nP1n = 0.; // <<cos(n*(phi1+phi2))>>
687  
688  if(dMult>1)
689  {
690   cosP1nP1n = (pow(dReQ1n,2)-pow(dImQ1n,2)-dReQ2n)/(dMult*(dMult-1)); 
691   
692   // average non-weighted 2-particle correction (cos terms) for non-uniform acceptance for single event:
693   fQCorrectionsEBE[0][1]->SetBinContent(2,cosP1nP1n);
694   
695   // final average non-weighted 2-particle correction (cos terms) for non-uniform acceptance for all events:
696   fQCorrections[0][0][1]->Fill(1.5,cosP1nP1n,dMult*(dMult-1));  
697  } 
698  
699  // 3-particle:
700  Double_t cosP1nM1nM1n = 0.; // <<cos(n*(phi1-phi2-phi3))>>
701  
702  if(dMult>2)
703  {
704   cosP1nM1nM1n = (dReQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))-dReQ1n*dReQ2n-dImQ1n*dImQ2n-2.*(dMult-1)*dReQ1n)
705                / (dMult*(dMult-1)*(dMult-2)); 
706   
707   // average non-weighted 3-particle correction (cos terms) for non-uniform acceptance for single event:
708   fQCorrectionsEBE[0][1]->SetBinContent(3,cosP1nM1nM1n);
709   
710   // final average non-weighted 3-particle correction (cos terms) for non-uniform acceptance for all events:
711   fQCorrections[0][0][1]->Fill(2.5,cosP1nM1nM1n,dMult*(dMult-1)*(dMult-2));  
712  } 
713  
714 } // end of AliFlowAnalysisWithQCumulants::CalculateCorrectionsForNonUniformAcceptanceForIntFlowCosTerms()
715
716
717 //================================================================================================================================
718
719
720 void AliFlowAnalysisWithQCumulants::CalculateCorrectionsForNonUniformAcceptanceForIntFlowSinTerms()
721 {
722  // calculate corrections for non-uniform acceptance of the detector for no-name integrated flow (sin terms)
723  
724  // multiplicity:
725  Double_t dMult = (*fSMpk)(0,0);
726  
727  // real and imaginary parts of non-weighted Q-vectors evaluated in harmonics n, 2n, 3n and 4n: 
728  Double_t dReQ1n = (*fReQ)(0,0);
729  Double_t dReQ2n = (*fReQ)(1,0);
730  //Double_t dReQ3n = (*fReQ)(2,0);
731  //Double_t dReQ4n = (*fReQ)(3,0);
732  Double_t dImQ1n = (*fImQ)(0,0);
733  Double_t dImQ2n = (*fImQ)(1,0);
734  //Double_t dImQ3n = (*fImQ)(2,0);
735  //Double_t dImQ4n = (*fImQ)(3,0);
736         
737  //                                  *************************************************************
738  //                                  **** corrections for non-uniform acceptance (sin terms): ****
739  //                                  *************************************************************
740  //
741  // Remark 1: corrections for non-uniform acceptance (sin terms) calculated with non-weighted Q-vectors 
742  //           are stored in 1D profile fQCorrectionsSin.
743  // Remark 2: binning of fQCorrectionsSin is organized as follows:
744  // --------------------------------------------------------------------------------------------------------------------
745  // 1st bin: <<sin(n*(phi1))>> = sinP1n
746  // 2nd bin: <<sin(n*(phi1+phi2))>> = sinP1nP1n
747  // 3rd bin: <<sin(n*(phi1-phi2-phi3))>> = sinP1nM1nM1n
748  // ...
749  // --------------------------------------------------------------------------------------------------------------------
750  
751  // 1-particle:
752  Double_t sinP1n = 0.; // <sin(n*(phi1))>
753  
754  if(dMult>0)
755  {
756   sinP1n = dImQ1n/dMult; 
757      
758   // average non-weighted 1-particle correction (sin terms) for non-uniform acceptance for single event:
759   fQCorrectionsEBE[0][0]->SetBinContent(1,sinP1n);
760   
761   // final average non-weighted 1-particle correction (sin terms) for non-uniform acceptance for all events:   
762   fQCorrections[0][0][0]->Fill(0.5,sinP1n,dMult);  
763  } 
764  
765  // 2-particle:
766  Double_t sinP1nP1n = 0.; // <<sin(n*(phi1+phi2))>>
767  
768  if(dMult>1)
769  {
770   sinP1nP1n = (2.*dReQ1n*dImQ1n-dImQ2n)/(dMult*(dMult-1)); 
771      
772   // average non-weighted 2-particle correction (sin terms) for non-uniform acceptance for single event:
773   fQCorrectionsEBE[0][0]->SetBinContent(2,sinP1nP1n);
774   
775   // final average non-weighted 1-particle correction (sin terms) for non-uniform acceptance for all events:      
776   fQCorrections[0][0][0]->Fill(1.5,sinP1nP1n,dMult*(dMult-1));  
777  } 
778  
779  // 3-particle:
780  Double_t sinP1nM1nM1n = 0.; // <<sin(n*(phi1-phi2-phi3))>>
781  
782  if(dMult>2)
783  {
784   sinP1nM1nM1n = (-dImQ1n*(pow(dReQ1n,2)+pow(dImQ1n,2))+dReQ1n*dImQ2n-dImQ1n*dReQ2n+2.*(dMult-1)*dImQ1n)
785                / (dMult*(dMult-1)*(dMult-2)); 
786   
787   // average non-weighted 3-particle correction (sin terms) for non-uniform acceptance for single event:
788   fQCorrectionsEBE[0][0]->SetBinContent(3,sinP1nM1nM1n);
789   
790   // final average non-weighted 3-particle correction (sin terms) for non-uniform acceptance for all events:  
791   fQCorrections[0][0][0]->Fill(2.5,sinP1nM1nM1n,dMult*(dMult-1)*(dMult-2));  
792  } 
793  
794 } // end of AliFlowAnalysisWithQCumulants::CalculateCorrectionsForNonUniformAcceptanceForIntFlowSinTerms()
795
796
797 //================================================================================================================================
798
799
800 void AliFlowAnalysisWithQCumulants::CalculateCorrectionsForNonUniformAcceptanceForDifferentialFlowCosTerms(TString type)
801 {
802  // calculate corrections for non-uniform acceptance of the detector for differential flow (cos terms)
803  
804  // multiplicity:
805  //Double_t dMult = (*fSMpk)(0,0);
806  
807  // real and imaginary parts of non-weighted Q-vectors evaluated in harmonics n, 2n, 3n and 4n: 
808  //Double_t dReQ1n = (*fReQ)(0,0);
809  //Double_t dReQ2n = (*fReQ)(1,0);
810  //Double_t dReQ3n = (*fReQ)(2,0);
811  //Double_t dReQ4n = (*fReQ)(3,0);
812  //Double_t dImQ1n = (*fImQ)(0,0);
813  //Double_t dImQ2n = (*fImQ)(1,0);
814  //Double_t dImQ3n = (*fImQ)(2,0);
815  //Double_t dImQ4n = (*fImQ)(3,0);
816
817  // looping over all (pt,eta) bins and calculating correlations needed for differential flow: 
818  for(Int_t p=1;p<=fnBinsPt;p++)
819  {
820   for(Int_t e=1;e<=fnBinsEta;e++)
821   {
822    // real and imaginary parts of q_n (non-weighted Q-vector evaluated only for POIs in harmonic n for each (pt,eta) bin): 
823    //Double_t dReqnPtEta = 0.;
824    //Double_t dImqnPtEta = 0.;
825
826    // number of POIs in each (pt,eta) bin:
827    Double_t dmPtEta = 0.;
828
829    // real and imaginary parts of q''_{n}, q''_{2n}, ... 
830    // (non-weighted Q-vectors evaluated only for particles which are both RPs and POIs in harmonic n, 2n, ... for each (pt,eta) bin): 
831    //Double_t dReqPrimePrime1nPtEta = 0.;
832    //Double_t dImqPrimePrime1nPtEta = 0.;
833    //Double_t dReqPrimePrime2nPtEta = 0.;
834    //Double_t dImqPrimePrime2nPtEta = 0.;
835
836    // number of particles which are both RPs and POIs in each (pt,eta) bin:
837    //Double_t dmPrimePrimePtEta = 0.;
838    
839    if(type == "POI")
840    {
841     // q''_{n}, q''_{2n}:
842     //...............................................................................................
843     //dReqPrimePrime1nPtEta = fReqPrimePrime1nPtEta->GetBinContent(fReqPrimePrime1nPtEta->GetBin(p,e));
844     //dImqPrimePrime1nPtEta = fImqPrimePrime1nPtEta->GetBinContent(fImqPrimePrime1nPtEta->GetBin(p,e));
845     //dReqPrimePrime2nPtEta = fReqPrimePrime2nPtEta->GetBinContent(fReqPrimePrime2nPtEta->GetBin(p,e));
846     //dImqPrimePrime2nPtEta = fImqPrimePrime2nPtEta->GetBinContent(fImqPrimePrime2nPtEta->GetBin(p,e));
847     //...............................................................................................
848    
849     // m'':
850     //dmPrimePrimePtEta = fmPrimePrimePtEta->GetBinContent(fmPrimePrimePtEta->GetBin(p,e));
851    
852     // q'_{n}: 
853     //dReqnPtEta = fReqnPtEta->GetBinContent(fReqnPtEta->GetBin(p,e));
854     //dImqnPtEta = fImqnPtEta->GetBinContent(fImqnPtEta->GetBin(p,e));
855     //dmPtEta    = fmPtEta->GetBinContent(fmPtEta->GetBin(p,e));
856    }
857    else if(type == "RP")
858    {
859     // q_RP{n}, q_RP{2n}:
860     //...............................................................................................
861     //dReqPrimePrime1nPtEta = fReqRP1nPtEta->GetBinContent(fReqRP1nPtEta->GetBin(p,e));
862     //dImqPrimePrime1nPtEta = fImqRP1nPtEta->GetBinContent(fImqRP1nPtEta->GetBin(p,e));
863     //dReqPrimePrime2nPtEta = fReqRP2nPtEta->GetBinContent(fReqRP2nPtEta->GetBin(p,e));
864     //dImqPrimePrime2nPtEta = fImqRP2nPtEta->GetBinContent(fImqRP2nPtEta->GetBin(p,e));
865     //...............................................................................................
866    
867     // m'':
868     //dmPrimePrimePtEta = fmRPPtEta->GetBinContent(fmRPPtEta->GetBin(p,e));
869    
870     //dReqnPtEta = fReqRP1nPtEta->GetBinContent(fReqRP1nPtEta->GetBin(p,e)); // not a bug ;-)
871     //dImqnPtEta = fImqRP1nPtEta->GetBinContent(fImqRP1nPtEta->GetBin(p,e)); // not a bug ;-)
872     //dmPtEta    = fmRPPtEta->GetBinContent(fmRPPtEta->GetBin(p,e));         // not a bug ;-) 
873    }
874    
875    // 1'-p correction:
876    //Double_t oneCosP1nPsiPtEta = 0.;
877    
878    if(dmPtEta)
879    {
880     //oneCosP1nPsiPtEta = dReqnPtEta/dmPtEta;
881    
882     // fill the 2D profile to get the average 1'-p correction for each (pt, eta) bin:
883     if(type == "POI")
884     { 
885      //fCorrectionsCosP1nPsiPtEtaPOI->Fill(fPtMin+(p-1)*fPtBinWidth,fEtaMin+(e-1)*fEtaBinWidth,
886      //                                    oneCosP1nPsiPtEta,dmPtEta);
887     }
888     else if(type == "RP")
889     {
890      //fCorrectionsCosP1nPsiPtEtaRP->Fill(fPtMin+(p-1)*fPtBinWidth,fEtaMin+(e-1)*fEtaBinWidth,
891      //                                    oneCosP1nPsiPtEta,dmPtEta);
892     }
893    } // end of if(dmPtEta*dMult-dmPrimePrimePtEta)
894    
895    /*
896    
897    // 4'-particle correlation:
898    Double_t four1n1n1n1nPtEta = 0.;
899    if((dmPtEta-dmPrimePrimePtEta)*dMult*(dMult-1.)*(dMult-2.)
900        + dmPrimePrimePtEta*(dMult-1.)*(dMult-2.)*(dMult-3.)) // to be improved (introduce a new variable for this expression)
901    {
902     four1n1n1n1nPtEta = ((pow(dReQ1n,2.)+pow(dImQ1n,2.))*(dReqnPtEta*dReQ1n+dImqnPtEta*dImQ1n)
903                       - dReqPrimePrime2nPtEta*(pow(dReQ1n,2.)-pow(dImQ1n,2.))
904                       - 2.*dImqPrimePrime2nPtEta*dReQ1n*dImQ1n
905                       - dReqnPtEta*(dReQ1n*dReQ2n+dImQ1n*dImQ2n)
906                       + dImqnPtEta*(dImQ1n*dReQ2n-dReQ1n*dImQ2n)
907                       - 2.*dMult*(dReqnPtEta*dReQ1n+dImqnPtEta*dImQ1n)
908                       - 2.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))*dmPrimePrimePtEta                      
909                       + 6.*(dReqPrimePrime1nPtEta*dReQ1n+dImqPrimePrime1nPtEta*dImQ1n)                                            
910                       + 1.*(dReqPrimePrime2nPtEta*dReQ2n+dImqPrimePrime2nPtEta*dImQ2n)                      
911                       + 2.*(dReqnPtEta*dReQ1n+dImqnPtEta*dImQ1n)                       
912                       + 2.*dmPrimePrimePtEta*dMult                      
913                       - 6.*dmPrimePrimePtEta)        
914                       / ((dmPtEta-dmPrimePrimePtEta)*dMult*(dMult-1.)*(dMult-2.)
915                           + dmPrimePrimePtEta*(dMult-1.)*(dMult-2.)*(dMult-3.)); 
916     
917     // fill the 2D profile to get the average correlation for each (pt, eta) bin:
918     if(type == "POI")
919     {
920      f4pPtEtaPOI->Fill(fPtMin+(p-1)*fPtBinWidth,fEtaMin+(e-1)*fEtaBinWidth,four1n1n1n1nPtEta,
921                        (dmPtEta-dmPrimePrimePtEta)*dMult*(dMult-1.)*(dMult-2.)
922                         + dmPrimePrimePtEta*(dMult-1.)*(dMult-2.)*(dMult-3.));
923     }
924     else if(type == "RP")
925     {
926      f4pPtEtaRP->Fill(fPtMin+(p-1)*fPtBinWidth,fEtaMin+(e-1)*fEtaBinWidth,four1n1n1n1nPtEta,
927                       (dmPtEta-dmPrimePrimePtEta)*dMult*(dMult-1.)*(dMult-2.)
928                        + dmPrimePrimePtEta*(dMult-1.)*(dMult-2.)*(dMult-3.));   
929     }
930    } // end of if((dmPtEta-dmPrimePrimePtEta)*dMult*(dMult-1.)*(dMult-2.)
931      //            +dmPrimePrimePtEta*(dMult-1.)*(dMult-2.)*(dMult-3.))
932    
933   */
934    
935   } // end of for(Int_t e=1;e<=fnBinsEta;e++)
936  } // end of for(Int_t p=1;p<=fnBinsPt;p++)
937  
938 } // end of AliFlowAnalysisWithQCumulants::CalculateCorrectionsForNonUniformAcceptanceForDifferentialFlowCosTerms(TString type)
939
940
941 //================================================================================================================================
942
943
944 void AliFlowAnalysisWithQCumulants::CalculateCorrectionsForNonUniformAcceptanceForDifferentialFlowSinTerms(TString type)
945 {
946  // calculate corrections for non-uniform acceptance of the detector for differential flow (sin terms)
947  
948  // multiplicity:
949  //Double_t dMult = (*fSMpk)(0,0);
950  
951  // real and imaginary parts of non-weighted Q-vectors evaluated in harmonics n, 2n, 3n and 4n: 
952  //Double_t dReQ1n = (*fReQ)(0,0);
953  //Double_t dReQ2n = (*fReQ)(1,0);
954  //Double_t dReQ3n = (*fReQ)(2,0);
955  //Double_t dReQ4n = (*fReQ)(3,0);
956  //Double_t dImQ1n = (*fImQ)(0,0);
957  //Double_t dImQ2n = (*fImQ)(1,0);
958  //Double_t dImQ3n = (*fImQ)(2,0);
959  //Double_t dImQ4n = (*fImQ)(3,0);
960
961  // looping over all (pt,eta) bins and calculating correlations needed for differential flow: 
962  for(Int_t p=1;p<=fnBinsPt;p++)
963  {
964   for(Int_t e=1;e<=fnBinsEta;e++)
965   {
966    // real and imaginary parts of q_n (non-weighted Q-vector evaluated only for POIs in harmonic n for each (pt,eta) bin): 
967    //Double_t dReqnPtEta = 0.;
968    //Double_t dImqnPtEta = 0.;
969
970    // number of POIs in each (pt,eta) bin:
971    Double_t dmPtEta = 0.;
972
973    // real and imaginary parts of q''_{n}, q''_{2n}, ... 
974    // (non-weighted Q-vectors evaluated only for particles which are both RPs and POIs in harmonic n, 2n, ... for each (pt,eta) bin): 
975    //Double_t dReqPrimePrime1nPtEta = 0.;
976    //Double_t dImqPrimePrime1nPtEta = 0.;
977    //Double_t dReqPrimePrime2nPtEta = 0.;
978    //Double_t dImqPrimePrime2nPtEta = 0.;
979
980    // number of particles which are both RPs and POIs in each (pt,eta) bin:
981    //Double_t dmPrimePrimePtEta = 0.;
982    
983    if(type == "POI")
984    {
985     // q''_{n}, q''_{2n}:
986     //...............................................................................................
987     //dReqPrimePrime1nPtEta = fReqPrimePrime1nPtEta->GetBinContent(fReqPrimePrime1nPtEta->GetBin(p,e));
988     //dImqPrimePrime1nPtEta = fImqPrimePrime1nPtEta->GetBinContent(fImqPrimePrime1nPtEta->GetBin(p,e));
989     //dReqPrimePrime2nPtEta = fReqPrimePrime2nPtEta->GetBinContent(fReqPrimePrime2nPtEta->GetBin(p,e));
990     //dImqPrimePrime2nPtEta = fImqPrimePrime2nPtEta->GetBinContent(fImqPrimePrime2nPtEta->GetBin(p,e));
991     //...............................................................................................
992    
993     // m'':
994     //dmPrimePrimePtEta = fmPrimePrimePtEta->GetBinContent(fmPrimePrimePtEta->GetBin(p,e));
995    
996     // q'_{n}: 
997     //dReqnPtEta = fReqnPtEta->GetBinContent(fReqnPtEta->GetBin(p,e));
998     //dImqnPtEta = fImqnPtEta->GetBinContent(fImqnPtEta->GetBin(p,e));
999     //dmPtEta    = fmPtEta->GetBinContent(fmPtEta->GetBin(p,e));
1000    }
1001    else if(type == "RP")
1002    {
1003     // q_RP{n}, q_RP{2n}:
1004     //...............................................................................................
1005     //dReqPrimePrime1nPtEta = fReqRP1nPtEta->GetBinContent(fReqRP1nPtEta->GetBin(p,e));
1006     //dImqPrimePrime1nPtEta = fImqRP1nPtEta->GetBinContent(fImqRP1nPtEta->GetBin(p,e));
1007     //dReqPrimePrime2nPtEta = fReqRP2nPtEta->GetBinContent(fReqRP2nPtEta->GetBin(p,e));
1008     //dImqPrimePrime2nPtEta = fImqRP2nPtEta->GetBinContent(fImqRP2nPtEta->GetBin(p,e));
1009     //...............................................................................................
1010    
1011     // m'':
1012     //dmPrimePrimePtEta = fmRPPtEta->GetBinContent(fmRPPtEta->GetBin(p,e));
1013    
1014     //dReqnPtEta = fReqRP1nPtEta->GetBinContent(fReqRP1nPtEta->GetBin(p,e)); // not a bug ;-)
1015     //dImqnPtEta = fImqRP1nPtEta->GetBinContent(fImqRP1nPtEta->GetBin(p,e)); // not a bug ;-)
1016     //dmPtEta    = fmRPPtEta->GetBinContent(fmRPPtEta->GetBin(p,e));         // not a bug ;-) 
1017    }
1018    
1019    // 1'-p correction:
1020    //Double_t oneSinP1nPsiPtEta = 0.;
1021    
1022    if(dmPtEta)
1023    {
1024     //oneSinP1nPsiPtEta = dImqnPtEta/dmPtEta;
1025    
1026     // fill the 2D profile to get the average 1'-p correction for each (pt, eta) bin:
1027     if(type == "POI")
1028     { 
1029      //fCorrectionsSinP1nPsiPtEtaPOI->Fill(fPtMin+(p-1)*fPtBinWidth,fEtaMin+(e-1)*fEtaBinWidth,
1030      //                                    oneSinP1nPsiPtEta,dmPtEta);
1031     }
1032     else if(type == "RP")
1033     {
1034      //fCorrectionsSinP1nPsiPtEtaRP->Fill(fPtMin+(p-1)*fPtBinWidth,fEtaMin+(e-1)*fEtaBinWidth,
1035      //                                    oneSinP1nPsiPtEta,dmPtEta);
1036     }
1037    } // end of if(dmPtEta*dMult-dmPrimePrimePtEta)
1038    
1039    /*
1040    
1041    // 4'-particle correlation:
1042    Double_t four1n1n1n1nPtEta = 0.;
1043    if((dmPtEta-dmPrimePrimePtEta)*dMult*(dMult-1.)*(dMult-2.)
1044        + dmPrimePrimePtEta*(dMult-1.)*(dMult-2.)*(dMult-3.)) // to be improved (introduce a new variable for this expression)
1045    {
1046     four1n1n1n1nPtEta = ((pow(dReQ1n,2.)+pow(dImQ1n,2.))*(dReqnPtEta*dReQ1n+dImqnPtEta*dImQ1n)
1047                       - dReqPrimePrime2nPtEta*(pow(dReQ1n,2.)-pow(dImQ1n,2.))
1048                       - 2.*dImqPrimePrime2nPtEta*dReQ1n*dImQ1n
1049                       - dReqnPtEta*(dReQ1n*dReQ2n+dImQ1n*dImQ2n)
1050                       + dImqnPtEta*(dImQ1n*dReQ2n-dReQ1n*dImQ2n)
1051                       - 2.*dMult*(dReqnPtEta*dReQ1n+dImqnPtEta*dImQ1n)
1052                       - 2.*(pow(dReQ1n,2.)+pow(dImQ1n,2.))*dmPrimePrimePtEta                      
1053                       + 6.*(dReqPrimePrime1nPtEta*dReQ1n+dImqPrimePrime1nPtEta*dImQ1n)                                            
1054                       + 1.*(dReqPrimePrime2nPtEta*dReQ2n+dImqPrimePrime2nPtEta*dImQ2n)                      
1055                       + 2.*(dReqnPtEta*dReQ1n+dImqnPtEta*dImQ1n)                       
1056                       + 2.*dmPrimePrimePtEta*dMult                      
1057                       - 6.*dmPrimePrimePtEta)        
1058                       / ((dmPtEta-dmPrimePrimePtEta)*dMult*(dMult-1.)*(dMult-2.)
1059                           + dmPrimePrimePtEta*(dMult-1.)*(dMult-2.)*(dMult-3.)); 
1060     
1061     // fill the 2D profile to get the average correlation for each (pt, eta) bin:
1062     if(type == "POI")
1063     {
1064      f4pPtEtaPOI->Fill(fPtMin+(p-1)*fPtBinWidth,fEtaMin+(e-1)*fEtaBinWidth,four1n1n1n1nPtEta,
1065                        (dmPtEta-dmPrimePrimePtEta)*dMult*(dMult-1.)*(dMult-2.)
1066                         + dmPrimePrimePtEta*(dMult-1.)*(dMult-2.)*(dMult-3.));
1067     }
1068     else if(type == "RP")
1069     {
1070      f4pPtEtaRP->Fill(fPtMin+(p-1)*fPtBinWidth,fEtaMin+(e-1)*fEtaBinWidth,four1n1n1n1nPtEta,
1071                       (dmPtEta-dmPrimePrimePtEta)*dMult*(dMult-1.)*(dMult-2.)
1072                        + dmPrimePrimePtEta*(dMult-1.)*(dMult-2.)*(dMult-3.));   
1073     }
1074    } // end of if((dmPtEta-dmPrimePrimePtEta)*dMult*(dMult-1.)*(dMult-2.)
1075      //            +dmPrimePrimePtEta*(dMult-1.)*(dMult-2.)*(dMult-3.))
1076    
1077   */
1078    
1079   } // end of for(Int_t e=1;e<=fnBinsEta;e++)
1080  } // end of for(Int_t p=1;p<=fnBinsPt;p++)
1081  
1082 } // end of AliFlowAnalysisWithQCumulants::CalculateCorrectionsForNonUniformAcceptanceForDifferentialFlowSinTerms(TString type)
1083
1084
1085 //================================================================================================================================
1086
1087
1088 void AliFlowAnalysisWithQCumulants::EvaluateNestedLoopsForDifferentialFlow(AliFlowEventSimple* anEvent)
1089 {
1090  // evaluate the nested loops relevant for differential flow (needed for cross-checking the results)
1091  
1092  Int_t nPrim = anEvent->NumberOfTracks(); 
1093  AliFlowTrackSimple *aftsTrack = NULL;
1094  
1095  Double_t psi1=0., phi2=0., phi3=0., phi4=0.;// phi5=0., phi6=0., phi7=0., phi8=0.;
1096  Double_t wPhi1=1., wPhi2=1., wPhi3=1., wPhi4=1.;// wPhi5=1., wPhi6=1., wPhi7=1., wPhi8=1.;
1097  
1098  Int_t n = fHarmonic; // to be improved
1099  
1100  //                                          ********************************************
1101  //                                          **** NESTED LOOPS FOR DIFFERENTIAL FLOW ****
1102  //                                          ******************************************** 
1103  
1104  // Remark 1: (pt,eta) bin in which the cross-checking will be performed is given by 1.1 < pt < 1.2 GeV and -0.55 < eta < -0.525 
1105  
1106  // Remark 2: multi-particle correlations needed for diff. flow calculated with nested loops without weights are stored in 1D profile  
1107  //           fDirectCorrelationsDiffFlow
1108  
1109  // Remark 3: multi-particle correlations needed for diff. flow calculated with nested loops with weights are stored in 1D profile  
1110  //           fDirectCorrelationsDiffFlowW;
1111  
1112  // Remark 4: binning of fDirectCorrelationsDiffFlow is organized as follows:
1113  //......................................................................................
1114  //       ---- bins 1-20: 2-particle correlations ----
1115  //  1st bin: <2'>_{1n|1n} = twoPrime1n1n = <cos(n*(psi1-phi2))>
1116  //       ---- bins 21-40: 3-particle correlations ----
1117  //       ---- bins 41-60: 4-particle correlations ----
1118  // 41st bin: <4'>_{1n,1n|1n,1n} = fourPrime1n1n1n1n  = <cos(n*(psi1+phi2-phi3-phi4))>
1119  //......................................................................................
1120  
1121  // Remark 5: binning of fDirectCorrelationsDiffFlow is organized as follows:
1122  //......................................................................................
1123  //       ---- bins 1-20: 2-particle correlations ----
1124  //  1st bin: twoPrime1n1nW0W1 = <w2 cos(n*(psi1-phi2))>
1125  //       ---- bins 21-40: 3-particle correlations ----
1126  //       ---- bins 41-60: 4-particle correlations ----
1127  // 41st bin: fourPrime1n1n1n1nW0W1W1W1 = <w2 w3 w4 cos(n*(psi1+phi2-phi3-phi4))>
1128  //......................................................................................
1129  
1130  // 2'-particle:
1131  for(Int_t i1=0;i1<nPrim;i1++)
1132  {
1133   aftsTrack=anEvent->GetTrack(i1);
1134   // POI condition (first particle in the correlator must be POI): 
1135   if(!((aftsTrack->Pt()>=1.1 && aftsTrack->Pt()<1.2) && (aftsTrack->Eta()>=-0.55 && aftsTrack->Eta()<-0.525) && (aftsTrack->InPOISelection())))continue;
1136   psi1=aftsTrack->Phi(); 
1137   if(fUsePhiWeights && fPhiWeights) wPhi1 = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(psi1*fnBinsPhi/TMath::TwoPi())));
1138   
1139   fDirectCorrectionsDiffFlowCos->Fill(0.,cos(1.*n*(psi1)),1.); // <<cos(n*(psi1))>>
1140   fDirectCorrectionsDiffFlowSin->Fill(0.,sin(1.*n*(psi1)),1.); // <<sin(n*(psi1))>>
1141   
1142   for(Int_t i2=0;i2<nPrim;i2++)
1143   {
1144    if(i2==i1)continue;
1145    aftsTrack=anEvent->GetTrack(i2);
1146    // RP condition (!(first) particle in the correlator must be RP):
1147    if(!(aftsTrack->InRPSelection()))continue;
1148    phi2=aftsTrack->Phi();
1149    if(fUsePhiWeights && fPhiWeights) wPhi2 = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*fnBinsPhi/TMath::TwoPi())));
1150     
1151    // non-weighted: 
1152    //.....................................................................................  
1153    fDirectCorrelationsDiffFlow->Fill(0.,cos(1.*n*(psi1-phi2)),1.); // <cos(n*(psi1-phi2))>
1154    //.....................................................................................  
1155    // weighted:
1156    //.....................................................................................   
1157    if(fUsePhiWeights) fDirectCorrelationsDiffFlowW->Fill(0.,cos(1.*n*(psi1-phi2)),wPhi2); // <w2 cos(n*(psi1-phi2))>
1158    //.....................................................................................  
1159    
1160    //fDirectCorrelations->Fill(103.,cos(1.*n*(phi1-phi2)),pow(wPhi1,2)*wPhi2);//<2'>_{n,n}
1161    //fDirectCorrelations->Fill(104.,cos(2.*n*(phi1-phi2)),wPhi1*pow(wPhi2,2));//<2'>_{n,n}
1162    //fDirectCorrelations->Fill(105.,cos(1.*n*(phi1-phi2)),pow(wPhi2,3));//<2'>_{n,n}  
1163    //fDirectCorrelations->Fill(41.,cos(2.*n*(phi1-phi2)),1);//<2'>_{2n,2n}
1164    //fDirectCorrelations->Fill(42.,cos(3.*n*(phi1-phi2)),1);//<2'>_{3n,3n}
1165    //fDirectCorrelations->Fill(43.,cos(4.*n*(phi1-phi2)),1);//<2'>_{4n,4n}   
1166     
1167   }//end of for(Int_t i2=0;i2<nPrim;i2++)
1168  }//end of for(Int_t i1=0;i1<nPrim;i1++)
1169  
1170  
1171  
1172  /*
1173  
1174  //<3'>_{2n|n,n}
1175  for(Int_t i1=0;i1<nPrim;i1++)
1176  {
1177   aftsTrack=anEvent->GetTrack(i1);
1178   if(!((aftsTrack->Pt()>=0.5&&aftsTrack->Pt()<0.6)&&(aftsTrack->InPOISelection())))continue;//POI condition
1179   psi1=aftsTrack->Phi();
1180   if(fUsePhiWeights && fPhiWeights) wPhi1 = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(psi1*fnBinsPhi/TMath::TwoPi())));
1181   for(Int_t i2=0;i2<nPrim;i2++)
1182   {
1183    if(i2==i1)continue;
1184    aftsTrack=anEvent->GetTrack(i2);
1185    if(!(aftsTrack->InRPSelection()))continue;//RP condition
1186    phi2=aftsTrack->Phi();
1187    if(fUsePhiWeights && fPhiWeights) wPhi2 = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*fnBinsPhi/TMath::TwoPi())));
1188    for(Int_t i3=0;i3<nPrim;i3++)
1189    {
1190     if(i3==i1||i3==i2)continue;
1191     aftsTrack=anEvent->GetTrack(i3);
1192     if(!(aftsTrack->InRPSelection()))continue;//RP condition
1193     phi3=aftsTrack->Phi();
1194     if(fUsePhiWeights && fPhiWeights) wPhi3 = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*fnBinsPhi/TMath::TwoPi())));
1195     //fill the fDirectCorrelations:     
1196     
1197     // 2-p
1198     //fDirectCorrelations->Fill(101.,cos(n*(phi2-phi3)),wPhi1*wPhi2*wPhi3); // <w1 w2 w3 cos(n(phi2-phi3))>
1199     //fDirectCorrelations->Fill(102.,cos(n*(phi1-phi3)),pow(wPhi2,2.)*wPhi3); // <w2^2 w3 cos(n(psi1-phi2))>
1200     
1201     // 3-p            
1202     //fDirectCorrelations->Fill(110.,cos(n*(2.*phi1-phi2-phi3)),wPhi1*wPhi2*wPhi3); // <w1 w2 w3 cos(n(2psi1-phi2-phi3))>
1203     //fDirectCorrelations->Fill(111.,cos(n*(phi1+phi2-2.*phi3)),wPhi2*pow(wPhi3,2.)); // <w2 w3^2 cos(n(psi1+phi2-2.*phi3))>
1204     
1205     
1206     //fDirectCorrelations->Fill(46.,cos(n*(phi1+phi2-2.*phi3)),1);//<3'>_{n,n|2n}    
1207    }//end of for(Int_t i3=0;i3<nPrim;i3++)  
1208   }//end of for(Int_t i2=0;i2<nPrim;i2++)  
1209  }//end of for(Int_t i1=0;i1<nPrim;i1++)
1210  */
1211  
1212  // 4'-particle:
1213  for(Int_t i1=0;i1<nPrim;i1++)
1214  {
1215   aftsTrack=anEvent->GetTrack(i1);
1216   // POI condition (first particle in the correlator must be POI): 
1217   if(!((aftsTrack->Pt()>=1.1 && aftsTrack->Pt()<1.2) && (aftsTrack->Eta()>=-0.55 && aftsTrack->Eta()<-0.525) && (aftsTrack->InPOISelection())))continue;
1218   psi1=aftsTrack->Phi();
1219   if(fUsePhiWeights && fPhiWeights) wPhi1 = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(psi1*fnBinsPhi/TMath::TwoPi())));
1220   for(Int_t i2=0;i2<nPrim;i2++)
1221   {
1222    if(i2==i1)continue;
1223    aftsTrack=anEvent->GetTrack(i2);
1224    // RP condition (!(first) particle in the correlator must be RP): 
1225    if(!(aftsTrack->InRPSelection()))continue;
1226    phi2=aftsTrack->Phi();
1227    if(fUsePhiWeights && fPhiWeights) wPhi2 = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi2*fnBinsPhi/TMath::TwoPi())));
1228    for(Int_t i3=0;i3<nPrim;i3++)
1229    { 
1230     if(i3==i1||i3==i2)continue;
1231     aftsTrack=anEvent->GetTrack(i3);
1232     // RP condition (!(first) particle in the correlator must be RP):
1233     if(!(aftsTrack->InRPSelection()))continue;
1234     phi3=aftsTrack->Phi();
1235     if(fUsePhiWeights && fPhiWeights) wPhi3 = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi3*fnBinsPhi/TMath::TwoPi())));
1236     for(Int_t i4=0;i4<nPrim;i4++)
1237     {
1238      if(i4==i1||i4==i2||i4==i3)continue;
1239      aftsTrack=anEvent->GetTrack(i4);
1240      // RP condition (!(first) particle in the correlator must be RP):
1241      if(!(aftsTrack->InRPSelection()))continue;  
1242      phi4=aftsTrack->Phi();
1243      if(fUsePhiWeights && fPhiWeights) wPhi4 = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(phi4*fnBinsPhi/TMath::TwoPi())));
1244      
1245      // non-weighted:
1246      //.........................................................................................................................
1247      fDirectCorrelationsDiffFlow->Fill(40.,cos(n*(psi1+phi2-phi3-phi4)),1.); // <cos(n(psi1+phi1-phi2-phi3))> 
1248      //.........................................................................................................................     
1249      // weighted:
1250      //...............................................................................................................................
1251      if(fUsePhiWeights) fDirectCorrelationsDiffFlowW->Fill(40.,cos(n*(psi1+phi2-phi3-phi4)),wPhi2*wPhi3*wPhi4); // <w2 w3 w4 cos(n(psi1+phi2-phi3-phi4))> 
1252      //...............................................................................................................................     
1253           
1254     }//end of for(Int_t i4=0;i4<nPrim;i4++)
1255    }//end of for(Int_t i3=0;i3<nPrim;i3++)
1256   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
1257  }//end of for(Int_t i1=0;i1<nPrim;i1++)
1258  
1259  /*                
1260  //<5'>_{2n,n|n,n,n}
1261  for(Int_t i1=0;i1<nPrim;i1++)
1262  {
1263   aftsTrack=anEvent->GetTrack(i1);
1264   if(!((aftsTrack->Pt()>=0.5&&aftsTrack->Pt()<0.6)&&(aftsTrack->InPOISelection())))continue;//POI condition
1265   phi1=aftsTrack->Phi();
1266   for(Int_t i2=0;i2<nPrim;i2++)
1267   {
1268    if(i2==i1)continue;
1269    aftsTrack=anEvent->GetTrack(i2);
1270    if(!(aftsTrack->InRPSelection()))continue;//RP condition   
1271    phi2=aftsTrack->Phi();
1272    for(Int_t i3=0;i3<nPrim;i3++)
1273    { 
1274     if(i3==i1||i3==i2)continue;
1275     aftsTrack=anEvent->GetTrack(i3);
1276     if(!(aftsTrack->InRPSelection()))continue;//RP condition   
1277     phi3=aftsTrack->Phi();
1278     for(Int_t i4=0;i4<nPrim;i4++)
1279     {
1280      if(i4==i1||i4==i2||i4==i3)continue;
1281      aftsTrack=anEvent->GetTrack(i4);
1282      if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1283      phi4=aftsTrack->Phi();//
1284      for(Int_t i5=0;i5<nPrim;i5++)
1285      {
1286       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
1287       aftsTrack=anEvent->GetTrack(i5);
1288       if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1289       phi5=aftsTrack->Phi();    
1290       //fill the fDirectCorrelations:if(bNestedLoops)
1291       //fDirectCorrelations->Fill(55.,cos(2.*n*phi1+n*phi2-n*phi3-n*phi4-n*phi5),1);//<5'>_{2n,n|n,n,n}
1292      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
1293     }//end of for(Int_t i4=0;i4<nPrim;i4++)
1294    }//end of for(Int_t i3=0;i3<nPrim;i3++)
1295   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
1296  }//end of for(Int_t i1=0;i1<nPrim;i1++)
1297  
1298
1299   
1300  */
1301  /*
1302  
1303  
1304  
1305  //<6'>_{n,n,n|n,n,n}
1306  for(Int_t i1=0;i1<nPrim;i1++)
1307  {
1308   aftsTrack=anEvent->GetTrack(i1);
1309   if(!((aftsTrack->Pt()>=0.5&&aftsTrack->Pt()<0.6)&&(aftsTrack->InPOISelection())))continue;//POI condition
1310   phi1=aftsTrack->Phi();
1311   for(Int_t i2=0;i2<nPrim;i2++)
1312   {
1313    if(i2==i1)continue;
1314    aftsTrack=anEvent->GetTrack(i2);
1315    if(!(aftsTrack->InRPSelection()))continue;//RP condition   
1316    phi2=aftsTrack->Phi();
1317    for(Int_t i3=0;i3<nPrim;i3++)
1318    { 
1319     if(i3==i1||i3==i2)continue;
1320     aftsTrack=anEvent->GetTrack(i3);
1321     if(!(aftsTrack->InRPSelection()))continue;//RP condition   
1322     phi3=aftsTrack->Phi();
1323     for(Int_t i4=0;i4<nPrim;i4++)
1324     {
1325      if(i4==i1||i4==i2||i4==i3)continue;
1326      aftsTrack=anEvent->GetTrack(i4);
1327      if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1328      phi4=aftsTrack->Phi();
1329      for(Int_t i5=0;i5<nPrim;i5++)
1330      {
1331       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
1332       aftsTrack=anEvent->GetTrack(i5);
1333       if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1334       phi5=aftsTrack->Phi();    
1335       for(Int_t i6=0;i6<nPrim;i6++)
1336       {
1337        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
1338        aftsTrack=anEvent->GetTrack(i6);
1339        if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1340        phi6=aftsTrack->Phi();  
1341        //fill the fDirectCorrelations:
1342        //fDirectCorrelations->Fill(60.,cos(n*(phi1+phi2+phi3-phi4-phi5-phi6)),1);//<6'>_{n,n,n|n,n,n}
1343       }//end of for(Int_t i6=0;i6<nPrim;i6++)   
1344      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
1345     }//end of for(Int_t i4=0;i4<nPrim;i4++)
1346    }//end of for(Int_t i3=0;i3<nPrim;i3++)
1347   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
1348  }//end of for(Int_t i1=0;i1<nPrim;i1++)
1349
1350  
1351  */
1352  /* 
1353    
1354      
1355  //<7'>_{2n,n,n|n,n,n,n}
1356  for(Int_t i1=0;i1<nPrim;i1++)
1357  {
1358   aftsTrack=anEvent->GetTrack(i1);
1359   if(!((aftsTrack->Pt()>=0.5&&aftsTrack->Pt()<0.6)&&(aftsTrack->InPOISelection())))continue;//POI condition
1360   phi1=aftsTrack->Phi();
1361   for(Int_t i2=0;i2<nPrim;i2++)
1362   {
1363    if(i2==i1)continue;
1364    aftsTrack=anEvent->GetTrack(i2);
1365    if(!(aftsTrack->InRPSelection()))continue;//RP condition   
1366    phi2=aftsTrack->Phi();
1367    for(Int_t i3=0;i3<nPrim;i3++)
1368    { 
1369     if(i3==i1||i3==i2)continue;
1370     aftsTrack=anEvent->GetTrack(i3);
1371     if(!(aftsTrack->InRPSelection()))continue;//RP condition   
1372     phi3=aftsTrack->Phi();
1373     for(Int_t i4=0;i4<nPrim;i4++)
1374     {
1375      if(i4==i1||i4==i2||i4==i3)continue;
1376      aftsTrack=anEvent->GetTrack(i4);
1377      if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1378      phi4=aftsTrack->Phi();
1379      for(Int_t i5=0;i5<nPrim;i5++)
1380      {
1381       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
1382       aftsTrack=anEvent->GetTrack(i5);
1383       if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1384       phi5=aftsTrack->Phi();    
1385       for(Int_t i6=0;i6<nPrim;i6++)
1386       {
1387        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
1388        aftsTrack=anEvent->GetTrack(i6);
1389        if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1390        phi6=aftsTrack->Phi();
1391        for(Int_t i7=0;i7<nPrim;i7++)
1392        {
1393         if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6)continue;
1394         aftsTrack=anEvent->GetTrack(i7);
1395         if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1396         phi7=aftsTrack->Phi();   
1397         //fill the fDirectCorrelations:
1398         //fDirectCorrelations->Fill(65.,cos(2.*n*phi1+n*phi2+n*phi3-n*phi4-n*phi5-n*phi6-n*phi7),1);//<7'>_{2n,n,n|n,n,n,n}
1399        }//end of for(Int_t i7=0;i7<nPrim;i7++)  
1400       }//end of for(Int_t i6=0;i6<nPrim;i6++)   
1401      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
1402     }//end of for(Int_t i4=0;i4<nPrim;i4++)
1403    }//end of for(Int_t i3=0;i3<nPrim;i3++)
1404   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
1405  }//end of for(Int_t i1=0;i1<nPrim;i1++)
1406
1407  
1408   
1409  */
1410  /*  
1411     
1412      
1413        
1414  //<8'>_{n,n,n,n|n,n,n,n}
1415  for(Int_t i1=0;i1<nPrim;i1++)
1416  {
1417   aftsTrack=anEvent->GetTrack(i1);
1418   if(!((aftsTrack->Pt()>=0.5&&aftsTrack->Pt()<0.6)&&(aftsTrack->InPOISelection())))continue;//POI condition
1419   phi1=aftsTrack->Phi();
1420   for(Int_t i2=0;i2<nPrim;i2++)
1421   {
1422    if(i2==i1)continue;
1423    aftsTrack=anEvent->GetTrack(i2);
1424    if(!(aftsTrack->InRPSelection()))continue;//RP condition   
1425    phi2=aftsTrack->Phi();
1426    for(Int_t i3=0;i3<nPrim;i3++)
1427    { 
1428     if(i3==i1||i3==i2)continue;
1429     aftsTrack=anEvent->GetTrack(i3);
1430     if(!(aftsTrack->InRPSelection()))continue;//RP condition   
1431     phi3=aftsTrack->Phi();
1432     for(Int_t i4=0;i4<nPrim;i4++)
1433     {
1434      if(i4==i1||i4==i2||i4==i3)continue;
1435      aftsTrack=anEvent->GetTrack(i4);
1436      if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1437      phi4=aftsTrack->Phi();
1438      for(Int_t i5=0;i5<nPrim;i5++)
1439      {
1440       if(i5==i1||i5==i2||i5==i3||i5==i4)continue;
1441       aftsTrack=anEvent->GetTrack(i5);
1442       if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1443       phi5=aftsTrack->Phi();    
1444       for(Int_t i6=0;i6<nPrim;i6++)
1445       {
1446        if(i6==i1||i6==i2||i6==i3||i6==i4||i6==i5)continue;
1447        aftsTrack=anEvent->GetTrack(i6);
1448        if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1449        phi6=aftsTrack->Phi();
1450        for(Int_t i7=0;i7<nPrim;i7++)
1451        {
1452         if(i7==i1||i7==i2||i7==i3||i7==i4||i7==i5||i7==i6)continue;
1453         aftsTrack=anEvent->GetTrack(i7);
1454         if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1455         phi7=aftsTrack->Phi();
1456         for(Int_t i8=0;i8<nPrim;i8++)
1457         {
1458          if(i8==i1||i8==i2||i8==i3||i8==i4||i8==i5||i8==i6||i8==i7)continue;
1459          aftsTrack=anEvent->GetTrack(i8);
1460          if(!(aftsTrack->InRPSelection()))continue;//RP condition  
1461          phi8=aftsTrack->Phi();           
1462          //fill the fDirectCorrelations:
1463          //fDirectCorrelations->Fill(70.,cos(n*(phi1+phi2+phi3+phi4-phi5-phi6-phi7-phi8)),1);//<8'>_{n,n,n,n|n,n,n,n}
1464         }//end of for(Int_t i8=0;i8<nPrim;i8++) 
1465        }//end of for(Int_t i7=0;i7<nPrim;i7++)  
1466       }//end of for(Int_t i6=0;i6<nPrim;i6++)   
1467      }//end of for(Int_t i5=0;i5<nPrim;i5++)  
1468     }//end of for(Int_t i4=0;i4<nPrim;i4++)
1469    }//end of for(Int_t i3=0;i3<nPrim;i3++)
1470   }//end of for(Int_t i2=0;i2<nPrim;i2++) 
1471  }//end of for(Int_t i1=0;i1<nPrim;i1++)
1472  
1473  
1474  
1475  */ 
1476  
1477  
1478  
1479  
1480 } // end of AliFlowAnalysisWithQCumulants::EvaluateNestedLoopsForDifferentialFlow(AliFlowEventSimple* anEvent)
1481
1482
1483 //================================================================================================================================
1484
1485
1486 void AliFlowAnalysisWithQCumulants::GetOutputHistograms(TList *outputListHistos)
1487 {
1488  // get pointers to all output histograms (called before Finish())
1489  if(outputListHistos)
1490  {      
1491   // 1.) common control histograms and common histograms for final results:
1492   AliFlowCommonHist *commonHist = dynamic_cast<AliFlowCommonHist*>(outputListHistos->FindObject("AliFlowCommonHistQC"));
1493   if(commonHist) this->SetCommonHists(commonHist); 
1494   AliFlowCommonHist *commonHist2nd = dynamic_cast<AliFlowCommonHist*>(outputListHistos->FindObject("AliFlowCommonHist2ndOrderQC"));
1495   if(commonHist2nd) this->SetCommonHists2nd(commonHist2nd); 
1496   AliFlowCommonHist *commonHist4th = dynamic_cast<AliFlowCommonHist*>(outputListHistos->FindObject("AliFlowCommonHist4thOrderQC"));
1497   if(commonHist4th) this->SetCommonHists4th(commonHist4th);
1498   AliFlowCommonHist *commonHist6th = dynamic_cast<AliFlowCommonHist*>(outputListHistos->FindObject("AliFlowCommonHist6thOrderQC"));
1499   if(commonHist6th) this->SetCommonHists6th(commonHist6th);
1500   AliFlowCommonHist *commonHist8th = dynamic_cast<AliFlowCommonHist*>(outputListHistos->FindObject("AliFlowCommonHist8thOrderQC"));
1501   if(commonHist8th) this->SetCommonHists8th(commonHist8th);
1502   AliFlowCommonHistResults *commonHistRes2nd = dynamic_cast<AliFlowCommonHistResults*>
1503                                                (outputListHistos->FindObject("AliFlowCommonHistResults2ndOrderQC"));
1504   if(commonHistRes2nd) this->SetCommonHistsResults2nd(commonHistRes2nd); 
1505   AliFlowCommonHistResults *commonHistRes4th = dynamic_cast<AliFlowCommonHistResults*>
1506                                                (outputListHistos->FindObject("AliFlowCommonHistResults4thOrderQC"));
1507   if(commonHistRes4th) this->SetCommonHistsResults4th(commonHistRes4th);
1508   AliFlowCommonHistResults *commonHistRes6th = dynamic_cast<AliFlowCommonHistResults*>
1509                                                (outputListHistos->FindObject("AliFlowCommonHistResults6thOrderQC"));
1510   if(commonHistRes6th) this->SetCommonHistsResults6th(commonHistRes6th);
1511   AliFlowCommonHistResults *commonHistRes8th = dynamic_cast<AliFlowCommonHistResults*>
1512                                                (outputListHistos->FindObject("AliFlowCommonHistResults8thOrderQC"));  
1513   if(commonHistRes8th) this->SetCommonHistsResults8th(commonHistRes8th);
1514   
1515   // 2.) weights: 
1516   TList *weightsList = dynamic_cast<TList*>(outputListHistos->FindObject("Weights"));
1517   if(weightsList) this->SetWeightsList(weightsList);
1518   Bool_t bUsePhiWeights = kFALSE;
1519   Bool_t bUsePtWeights = kFALSE;
1520   Bool_t bUseEtaWeights = kFALSE;
1521   TString fUseParticleWeightsName = "fUseParticleWeightsQC";
1522   fUseParticleWeightsName += fAnalysisLabel->Data();
1523   TProfile *useParticleWeights = dynamic_cast<TProfile*>(weightsList->FindObject(fUseParticleWeightsName.Data()));
1524   if(useParticleWeights)
1525   {
1526    this->SetUseParticleWeights(useParticleWeights);  
1527    bUsePhiWeights = (Int_t)useParticleWeights->GetBinContent(1);
1528    bUsePtWeights = (Int_t)useParticleWeights->GetBinContent(2);
1529    bUseEtaWeights = (Int_t)useParticleWeights->GetBinContent(3);
1530   }
1531   
1532   // 3.) integrated flow:
1533   TList *intFlowList = NULL;
1534   TList *intFlowProfiles = NULL;
1535   TList *intFlowResults = NULL;
1536    
1537   intFlowList = dynamic_cast<TList*>(outputListHistos->FindObject("Integrated Flow"));
1538   if(intFlowList) 
1539   {
1540    intFlowProfiles = dynamic_cast<TList*>(intFlowList->FindObject("Profiles"));
1541    intFlowResults = dynamic_cast<TList*>(intFlowList->FindObject("Results"));
1542   } else
1543     {
1544      cout<<"WARNING: intFlowList is NULL in AFAWQC::GOH() !!!!"<<endl;
1545     }  
1546    
1547   // profiles:   
1548   if(intFlowProfiles)  
1549   {
1550    // average multiplicities:
1551    TProfile *avMultiplicity = dynamic_cast<TProfile*>(intFlowProfiles->FindObject("fAvMultiplicity"));
1552    if(avMultiplicity) 
1553    {
1554     this->SetAvMultiplicity(avMultiplicity);
1555    } else 
1556      {
1557       cout<<"WARNING: avMultiplicity is NULL in AFAWQC::GOH() !!!!"<<endl;
1558      }
1559    
1560    // flags: (to be improved (united with other flags in this method))
1561    TString pWeightsFlag[2] = {"pWeights not used","pWeights used"};
1562    TString eWeightsFlag[2] = {"exact eWeights","non-exact eWeights"};
1563    //TString nuaFlag[2] = {"not corrected","corrected"};
1564    TString sinCosFlag[2] = {"sin","cos"};
1565     
1566    for(Int_t pW=0;pW<1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights);pW++)
1567    {
1568     for(Int_t eW=0;eW<1;eW++)
1569     {
1570      // correlations (profiles):
1571      TProfile *qCorrelations = dynamic_cast<TProfile*>(intFlowProfiles->FindObject(Form("fQCorrelations: %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1572      if(qCorrelations) 
1573      {
1574       this->SetQCorrelations(qCorrelations,pW,eW);
1575      } else 
1576        {
1577         cout<<"WARNING: qCorrelations is NULL in AFAWQC::GOH() !!!!"<<endl;
1578         cout<<"pW = "<<pW<<endl;
1579         cout<<"eW = "<<eW<<endl;
1580        } 
1581      // products (profiles):  
1582      TProfile *qProducts = dynamic_cast<TProfile*>(intFlowProfiles->FindObject(Form("fQProducts: %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1583      if(qProducts) 
1584      {
1585       this->SetQProducts(qProducts,pW,eW);
1586      } else 
1587        {
1588         cout<<"WARNING: qProducts is NULL in AFAWQC::GOH() !!!!"<<endl;
1589         cout<<"pW = "<<pW<<endl;
1590         cout<<"eW = "<<eW<<endl;
1591        }     
1592      // corrections (profiles):
1593      for(Int_t sc=0;sc<2;sc++)
1594      {
1595       TProfile *qCorrections = dynamic_cast<TProfile*>(intFlowProfiles->FindObject((Form("fQCorrections: %s, %s, %s terms",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data(),sinCosFlag[sc].Data()))));
1596       if(qCorrections) 
1597       {
1598        this->SetQCorrections(qCorrections,pW,eW,sc);
1599       } else 
1600         {
1601          cout<<"WARNING: qCorrections is NULL in AFAWQC::GOH() !!!!"<<endl;
1602          cout<<"pW = "<<pW<<endl;
1603          cout<<"eW = "<<eW<<endl;
1604          cout<<"sc = "<<sc<<endl;
1605         } 
1606      } // end of for(Int_t sc=0;sc<2;sc++)           
1607     } // end of for(Int_t eW=0;eW<1;eW++)
1608    } // end of for(Int_t pW=0;pW<1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights);pW++)
1609   } else // to if(intFlowProfiles)  
1610     {
1611      cout<<"WARNING: intFlowProfiles is NULL in AFAWQC::GOH() !!!!"<<endl;
1612     }
1613    
1614   // results:   
1615   if(intFlowResults)  
1616   {
1617    for(Int_t pW=0;pW<1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights);pW++)
1618    {
1619     for(Int_t eW=0;eW<2;eW++)
1620     {
1621      // flags: (to be improved (united with other flags in this method))
1622      TString pWeightsFlag[2] = {"pWeights not used","pWeights used"};
1623      TString eWeightsFlag[2] = {"exact eWeights","non-exact eWeights"};
1624      TString nuaFlag[2] = {"not corrected","corrected"};
1625      TString powerFlag[2] = {"linear","quadratic"};
1626      //TString sinCosFlag[2] = {"sin","cos"};     
1627      // correlations (results):
1628      TH1D *correlations = dynamic_cast<TH1D*>(intFlowResults->FindObject(Form("fCorrelations: %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1629      if(correlations) 
1630      {
1631       this->SetCorrelations(correlations,pW,eW);
1632      } else 
1633        {
1634         cout<<"WARNING: correlations is NULL in AFAWQC::GOH() !!!!"<<endl;
1635         cout<<"pW = "<<pW<<endl;
1636         cout<<"eW = "<<eW<<endl;
1637        }     
1638      // corrections (results):
1639      TH1D *corrections = dynamic_cast<TH1D*>(intFlowResults->FindObject(Form("fCorrections: %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1640      if(corrections) 
1641      {
1642       this->SetCorrections(corrections,pW,eW);
1643      } else 
1644        {
1645         cout<<"WARNING: corrections is NULL in AFAWQC::GOH() !!!!"<<endl;
1646         cout<<"pW = "<<pW<<endl;
1647         cout<<"eW = "<<eW<<endl;
1648        }
1649      // covariances (results):
1650      TH1D *covariances = dynamic_cast<TH1D*>(intFlowResults->FindObject(Form("fCovariances: %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1651      if(covariances) 
1652      {
1653       this->SetCovariances(covariances,pW,eW);
1654      } else 
1655        {
1656         cout<<"WARNING: covariances is NULL in AFAWQC::GOH() !!!!"<<endl;
1657         cout<<"pW = "<<pW<<endl;
1658         cout<<"eW = "<<eW<<endl;
1659        } 
1660      // sum of linear and quadratic event weights (results):
1661      for(Int_t power=0;power<2;power++)
1662      {
1663       TH1D *sumOfEventWeights = dynamic_cast<TH1D*>(intFlowResults->FindObject(Form("fSumOfEventWeights: %s, %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data(),powerFlag[power].Data())));
1664       if(sumOfEventWeights) 
1665       {
1666        this->SetSumOfEventWeights(sumOfEventWeights,pW,eW,power);
1667       } else 
1668         {
1669          cout<<"WARNING: sumOfEventWeights is NULL in AFAWQC::GOH() !!!!"<<endl;
1670          cout<<"pW    = "<<pW<<endl;
1671          cout<<"eW    = "<<eW<<endl;
1672          cout<<"power = "<<power<<endl;
1673         }                                   
1674      } // end of for(Int_t power=0;power<2;power++)                                                                  
1675      // products of event weights (results):
1676      TH1D *productOfEventWeights = dynamic_cast<TH1D*>(intFlowResults->FindObject(Form("fProductOfEventWeights: %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1677      if(productOfEventWeights) 
1678      {
1679       this->SetProductOfEventWeights(productOfEventWeights,pW,eW);
1680      } else 
1681        {
1682         cout<<"WARNING: productOfEventWeights is NULL in AFAWQC::GOH() !!!!"<<endl;
1683         cout<<"pW = "<<pW<<endl;
1684         cout<<"eW = "<<eW<<endl;
1685        } 
1686        
1687      for(Int_t nua=0;nua<2;nua++)
1688      {
1689       // integrated Q-cumulants:
1690       TH1D *cumulants = dynamic_cast<TH1D*>(intFlowResults->FindObject(Form("fCumulants: %s, %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data(),nuaFlag[nua].Data())));
1691       if(cumulants) 
1692       {
1693        this->SetCumulants(cumulants,pW,eW,nua);
1694       } else 
1695         {
1696          cout<<"WARNING: cumulants is NULL in AFAWQC::GOH() !!!!"<<endl;
1697          cout<<"pW = "<<pW<<endl;
1698          cout<<"eW = "<<eW<<endl;
1699          cout<<"nua = "<<nua<<endl;
1700         }  
1701       // integrated flow estimates from Q-cumulants:
1702       TH1D *intFlow = dynamic_cast<TH1D*>(intFlowResults->FindObject(Form("fIntFlow: %s, %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data(),nuaFlag[nua].Data())));
1703       if(intFlow) 
1704       {
1705        this->SetIntFlow(intFlow,pW,eW,nua);
1706       } else 
1707         {
1708          cout<<"WARNING: intFlow is NULL in AFAWQC::GOH() !!!!"<<endl;
1709          cout<<"pW = "<<pW<<endl;
1710          cout<<"eW = "<<eW<<endl;
1711          cout<<"nua = "<<nua<<endl;
1712         }   
1713      } // end of for(Int_t nua=0;nua<2;nua++)  
1714     } // end of for(Int_t eW=0;eW<1;eW++) 
1715    } // end of for(Int_t pW=0;pW<1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights);pW++)  
1716   } else // to if(intFlowResults)
1717     {
1718      cout<<"WARNING: intFlowResults is NULL in AFAWQC::GOH() !!!!"<<endl;
1719     }
1720           
1721   // differential flow:
1722   TString typeFlag[2] = {"RP","POI"}; 
1723   TString pWeightsFlag[2] = {"not used","used"};
1724   TString eWeightsFlag[2] = {"exact","non-exact"}; 
1725   TString sinCosFlag[2] = {"sin","cos"};
1726   TString nuaFlag[2] = {"not corrected","corrected"}; // nua = non-uniform acceptance
1727   TString ptEtaFlag[2] = {"p_{t}","#eta"};
1728   // base list fDiffFlowList "Differential Flow":
1729   TList *diffFlowList = NULL;
1730   diffFlowList = dynamic_cast<TList*>(outputListHistos->FindObject("Differential Flow"));  
1731   // list holding nested lists containing profiles:
1732   TList *diffFlowListProfiles = NULL;
1733   // list holding nested lists containing 2D and 1D histograms with final results:
1734   TList *diffFlowListResults = NULL;
1735   if(diffFlowList)
1736   {  
1737    diffFlowListProfiles = dynamic_cast<TList*>(diffFlowList->FindObject("Profiles"));
1738    diffFlowListResults = dynamic_cast<TList*>(diffFlowList->FindObject("Results"));
1739   } else
1740     {
1741      cout<<"WARNING: diffFlowList is NULL in AFAWQC::GOH() !!!!"<<endl;
1742     }      
1743   
1744   // nested list in the list of profiles fDiffFlowListProfiles "Profiles":
1745   TList *dfpType[2] = {NULL};
1746   TList *dfpParticleWeights[2][2] = {{NULL}};
1747   TList *dfpEventWeights[2][2][2] = {{{NULL}}};
1748   TList *diffFlowCorrelations[2][2][2] = {{{NULL}}};
1749   TList *diffFlowProductsOfCorrelations[2][2][2] = {{{NULL}}};
1750   TList *diffFlowCorrectionTerms[2][2][2][2] = {{{{NULL}}}};
1751   
1752   if(diffFlowListProfiles)
1753   {
1754    for(Int_t t=0;t<2;t++)
1755    {
1756     dfpType[t] = dynamic_cast<TList*>(diffFlowListProfiles->FindObject(typeFlag[t].Data()));
1757     if(dfpType[t])
1758     {
1759      for(Int_t pW=0;pW<(1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights));pW++)
1760      {
1761       dfpParticleWeights[t][pW] = dynamic_cast<TList*>(dfpType[t]->FindObject(Form("%s, pWeights %s",typeFlag[t].Data(),pWeightsFlag[pW].Data())));
1762       if(dfpParticleWeights[t][pW])
1763       {
1764        for(Int_t eW=0;eW<2;eW++)
1765        {
1766         dfpEventWeights[t][pW][eW] = dynamic_cast<TList*>(dfpParticleWeights[t][pW]->FindObject(Form("%s, pWeights %s, eWeights %s",typeFlag[t].Data(),pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1767         if(dfpEventWeights[t][pW][eW])
1768         {
1769          // correlations: 
1770          diffFlowCorrelations[t][pW][eW] = dynamic_cast<TList*>(dfpEventWeights[t][pW][eW]->FindObject(Form("Correlations (%s, pWeights %s, eWeights %s)",typeFlag[t].Data(),pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1771          // products of correlations:
1772          diffFlowProductsOfCorrelations[t][pW][eW] = dynamic_cast<TList*>(dfpEventWeights[t][pW][eW]->FindObject(Form("Products of correlations (%s, pWeights %s, eWeights %s)",typeFlag[t].Data(),pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1773          // correction terms:
1774          for(Int_t sc=0;sc<2;sc++)
1775          {
1776           diffFlowCorrectionTerms[t][pW][eW][sc] = dynamic_cast<TList*>(dfpEventWeights[t][pW][eW]->FindObject(Form("Corrections for NUA (%s, pWeights %s, eWeights %s, %s terms)",typeFlag[t].Data(),pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data(),sinCosFlag[sc].Data())));
1777           //this->SetDiffFlowCorrectionTerms(diffFlowCorrectionTerms[t][pW][sc],t,pW,sc);   
1778          }
1779         } else // to if(dfpEventWeights[t][pW][eW])
1780           {
1781            cout<<"WARNING: dfpEventWeights[t][pW][eW] is NULL in AFAWQC::GOH() !!!!"<<endl;
1782            cout<<"t = "<<t<<endl;
1783            cout<<"pW = "<<pW<<endl;
1784            cout<<"eW = "<<eW<<endl;
1785           }
1786        } // end of for(Int_t eW=0;eW<2;eW++)   
1787       } else // to if(dfpParticleWeights[t][pW])
1788         {
1789          cout<<"WARNING: dfpParticleWeights[t][pW] is NULL in AFAWQC::GOH() !!!!"<<endl;
1790          cout<<"t = "<<t<<endl;
1791          cout<<"pW = "<<pW<<endl;
1792         }
1793      } // end of for(Int_t pW=0;pW<(1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights));pW++) 
1794      } else // if(dfpType[t]) 
1795       {
1796        cout<<"WARNING: dfpType[t] is NULL in AFAWQC::GOH() !!!!"<<endl;
1797        cout<<"t = "<<t<<endl;
1798       }
1799    } // end of for(Int_t t=0;t<2;t++)
1800   } else // to if(diffFlowListProfiles)
1801     {
1802      cout<<"WARNING: diffFlowListProfiles is NULL in AFAWQC::GOH() !!!!"<<endl;
1803     }  
1804   
1805   TProfile2D *correlationsPro[2][2][2][4] = {{{{NULL}}}};
1806   TProfile2D *productsOfCorrelationsPro[2][2][2][5] = {{{{NULL}}}};
1807   TProfile2D *correctionTermsPro[2][2][2][2][2] = {{{{{NULL}}}}};
1808   
1809   TString correlationName[4] = {"<2'>","<4'>","<6'>","<8'>"};
1810   TString cumulantName[4] = {"QC{2'}","QC{4'}","QC{6'}","QC{8'}"};
1811   TString productOfCorrelationsName[5] = {"<2><2'>","<2><4'>","<4><2'>","<4><4'>","<2'><4'>"};
1812   TString correctionsSinTermsName[2] = {"sin(1)","sin(2)"};
1813   TString correctionsCosTermsName[2] = {"cos(1)","cos(2)"};
1814   TString correctionName[4] = {"corr. to QC{2'}","corr. to QC{4'}","corr. to QC{6'}","corr. to QC{8'}"};
1815   TString covarianceName[4] = {"Cov(2,2')","Cov(2,4')","Cov(2',4')","Cov(4,2')"}; // to be improved (reorganized)
1816   TString flowName[4] = {"v'{2}","v'{4}","v'{6}","v'{8}"}; 
1817   
1818   for(Int_t t=0;t<2;t++)
1819   {
1820    for(Int_t pW=0;pW<(1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights));pW++)
1821    {
1822     for(Int_t eW=0;eW<2;eW++)
1823     {
1824      // correlations:
1825      if(diffFlowCorrelations[t][pW][eW])
1826      {
1827       for(Int_t correlationIndex=0;correlationIndex<4;correlationIndex++)
1828       {
1829        correlationsPro[t][pW][eW][correlationIndex] = dynamic_cast<TProfile2D*>(diffFlowCorrelations[t][pW][eW]->FindObject(correlationName[correlationIndex].Data()));
1830        if(correlationsPro[t][pW][eW][correlationIndex])
1831        {
1832         this->SetCorrelationsPro(correlationsPro[t][pW][eW][correlationIndex],t,pW,eW,correlationIndex);
1833        } else // to if(correlationsPro[t][pW][ew][correlationIndex])
1834          {
1835           cout<<"WARNING: correlationsPro[t][pW][eW][correlationIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
1836           cout<<"t = "<<t<<endl;
1837           cout<<"pW = "<<pW<<endl;
1838           cout<<"eW = "<<eW<<endl;
1839           cout<<"ci = "<<correlationIndex<<endl;
1840          } 
1841       }
1842      } else // to if(diffFlowCorrelations[t][pW][eW])
1843        {
1844         cout<<"WARNING: diffFlowCorrelations[t][pW][eW] is NULL in AFAWQC::GOH() !!!!"<<endl;
1845         cout<<"t = "<<t<<endl;
1846         cout<<"pW = "<<pW<<endl;
1847         cout<<"eW = "<<eW<<endl;
1848        } 
1849      // products of correlations:
1850      if(diffFlowProductsOfCorrelations[t][pW][eW])
1851      {
1852       for(Int_t productOfCorrelationsIndex=0;productOfCorrelationsIndex<5;productOfCorrelationsIndex++)
1853       {
1854        productsOfCorrelationsPro[t][pW][eW][productOfCorrelationsIndex] = dynamic_cast<TProfile2D*>(diffFlowProductsOfCorrelations[t][pW][eW]->FindObject(productOfCorrelationsName[productOfCorrelationsIndex].Data()));
1855        if(productsOfCorrelationsPro[t][pW][eW][productOfCorrelationsIndex])
1856        {
1857         this->SetProductsOfCorrelationsPro(productsOfCorrelationsPro[t][pW][eW][productOfCorrelationsIndex],t,pW,eW,productOfCorrelationsIndex);
1858        } else // to if(productsOfCorrelationsPro[t][pW][eW][productOfCorrelationsIndex])
1859          {
1860           cout<<"WARNING: productsOfCorrelationsPro[t][pW][eW][productOfCorrelationsIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
1861           cout<<"t = "<<t<<endl;
1862           cout<<"pW = "<<pW<<endl;
1863           cout<<"eW = "<<eW<<endl;
1864           cout<<"ci = "<<productOfCorrelationsIndex<<endl;
1865          } 
1866       }
1867      } else // to if(diffFlowProductsOfCorrelations[t][pW][eW])
1868        {
1869         cout<<"WARNING: diffFlowProductsOfCorrelations[t][pW][eW] is NULL in AFAWQC::GOH() !!!!"<<endl;
1870         cout<<"t = "<<t<<endl;
1871         cout<<"pW = "<<pW<<endl;
1872         cout<<"eW = "<<eW<<endl;
1873        }
1874      // correction terms:
1875      for(Int_t sc=0;sc<2;sc++)
1876      {
1877       if(diffFlowCorrectionTerms[t][pW][eW][sc])
1878       {
1879        for(Int_t correctionIndex=0;correctionIndex<2;correctionIndex++)
1880        {
1881         if(sc==0)
1882         {
1883          correctionTermsPro[t][pW][eW][sc][correctionIndex] = dynamic_cast<TProfile2D*>(diffFlowCorrectionTerms[t][pW][eW][sc]->FindObject(correctionsSinTermsName[correctionIndex].Data())); 
1884          if(correctionTermsPro[t][pW][eW][sc][correctionIndex])
1885          {
1886           this->SetCorrectionTermsPro(correctionTermsPro[t][pW][eW][sc][correctionIndex],t,pW,eW,sc,correctionIndex);
1887          } else 
1888            {
1889             cout<<"WARNING: correctionTermsPro[t][pW][eW][sc][correctionIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
1890             cout<<"t = "<<t<<endl;
1891             cout<<"pW = "<<pW<<endl;
1892             cout<<"eW = "<<eW<<endl;
1893             cout<<"sc = "<<sc<<endl;
1894             cout<<"ci = "<<correctionIndex<<endl;
1895            }
1896         } 
1897         if(sc==1)
1898         {
1899          correctionTermsPro[t][pW][eW][sc][correctionIndex] = dynamic_cast<TProfile2D*>(diffFlowCorrectionTerms[t][pW][eW][sc]->FindObject(correctionsCosTermsName[correctionIndex].Data())); 
1900          if(correctionTermsPro[t][pW][eW][sc][correctionIndex])
1901          {
1902           this->SetCorrectionTermsPro(correctionTermsPro[t][pW][eW][sc][correctionIndex],t,pW,eW,sc,correctionIndex);
1903          } else 
1904            {
1905             cout<<"WARNING: correctionTermsPro[t][pW][eW][sc][correctionIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
1906             cout<<"t = "<<t<<endl;
1907             cout<<"pW = "<<pW<<endl;
1908             cout<<"eW = "<<eW<<endl;
1909             cout<<"sc = "<<sc<<endl;
1910             cout<<"ci = "<<correctionIndex<<endl;
1911            }         
1912         }  
1913        } // end of for(Int_t correctionIndex=0;correctionIndex<2;correctionIndex++)
1914       } else // to if(diffFlowCorrectionTerms[t][pW][eW][sc])
1915         {
1916          cout<<"WARNING: diffFlowCorrectionTerms[t][pW][eW][sc] is NULL in AFAWQC::GOH() !!!!"<<endl;
1917          cout<<"t = "<<t<<endl;
1918          cout<<"pW = "<<pW<<endl;
1919          cout<<"eW = "<<eW<<endl;
1920          cout<<"sc = "<<sc<<endl;
1921         }
1922      } // end of for(Int_t sc=0;sc<2;sc++)  
1923     } // end of for(Int_t eW=0;eW<2;eW++)
1924    }  // end of for(Int_t pW=0;pW<(1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights));pW++)
1925   } // end of for(Int_t t=0;t<2;t++)
1926   
1927   // nested list in the list of results fDiffFlowListResults "Results":
1928   TList *dfrType[2] = {NULL};
1929   TList *dfrParticleWeights[2][2] = {{NULL}};
1930   TList *dfrEventWeights[2][2][2] = {{{NULL}}};
1931   TList *dfrCorrections[2][2][2][2] = {{{{NULL}}}}; 
1932   TList *diffFlowFinalCorrelations[2][2][2] = {{{NULL}}}; 
1933   TList *diffFlowFinalCorrections[2][2][2] = {{{NULL}}}; 
1934   TList *diffFlowFinalCovariances[2][2][2] = {{{NULL}}};   
1935   TList *diffFlowFinalCumulants[2][2][2][2] = {{{{NULL}}}};  
1936   TList *diffFlowFinalFlow[2][2][2][2] = {{{{NULL}}}}; 
1937  
1938   if(diffFlowListResults)
1939   {
1940    for(Int_t t=0;t<2;t++)
1941    {
1942     dfrType[t] = dynamic_cast<TList*>(diffFlowListResults->FindObject(typeFlag[t].Data()));
1943     if(dfrType[t])
1944     {
1945      for(Int_t pW=0;pW<(1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights));pW++)
1946      {
1947       dfrParticleWeights[t][pW] = dynamic_cast<TList*>(dfrType[t]->FindObject(Form("%s, pWeights %s",typeFlag[t].Data(),pWeightsFlag[pW].Data())));
1948       if(dfrParticleWeights[t][pW])
1949       {
1950        for(Int_t eW=0;eW<2;eW++)
1951        {
1952         dfrEventWeights[t][pW][eW] = dynamic_cast<TList*>(dfrParticleWeights[t][pW]->FindObject(Form("%s, pWeights %s, eWeights %s",typeFlag[t].Data(),pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1953         if(dfrEventWeights[t][pW][eW])
1954         {
1955          diffFlowFinalCorrelations[t][pW][eW] = dynamic_cast<TList*>(dfrEventWeights[t][pW][eW]->FindObject(Form("Correlations (%s, pWeights %s, eWeights %s)",typeFlag[t].Data(),pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1956          diffFlowFinalCorrections[t][pW][eW] = dynamic_cast<TList*>(dfrEventWeights[t][pW][eW]->FindObject(Form("Corrections (%s, pWeights %s, eWeights %s)",typeFlag[t].Data(),pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));
1957          diffFlowFinalCovariances[t][pW][eW] = dynamic_cast<TList*>(dfrEventWeights[t][pW][eW]->FindObject(Form("Covariances (%s, pWeights %s, eWeights %s)",typeFlag[t].Data(),pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data())));      
1958          for(Int_t nua=0;nua<2;nua++)
1959          {
1960           dfrCorrections[t][pW][eW][nua] = dynamic_cast<TList*>(dfrEventWeights[t][pW][eW]->FindObject(Form("%s, pWeights %s, eWeights %s, %s for NUA",typeFlag[t].Data(),pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data(),nuaFlag[nua].Data())));
1961           if(dfrCorrections[t][pW][eW][nua])
1962           {
1963            diffFlowFinalCumulants[t][pW][eW][nua] = dynamic_cast<TList*>(dfrCorrections[t][pW][eW][nua]->FindObject(Form("Cumulants (%s, pWeights %s, eWeights %s, %s for NUA)",typeFlag[t].Data(),pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data(),nuaFlag[nua].Data())));
1964            diffFlowFinalFlow[t][pW][eW][nua] = dynamic_cast<TList*>(dfrCorrections[t][pW][eW][nua]->FindObject(Form("Differential Flow (%s, pWeights %s, eWeights %s, %s for NUA)",typeFlag[t].Data(),pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data(),nuaFlag[nua].Data())));
1965           } else // to if(dfrCorrections[t][pW][eW][nua])
1966             {
1967              cout<<"WARNING: dfrCorrections[t][pW][eW][nua] is NULL in AFAWQC::GOH() !!!!"<<endl;
1968              cout<<"t = "<<t<<endl;
1969              cout<<"pW = "<<pW<<endl;
1970              cout<<"eW = "<<eW<<endl;
1971              cout<<"nua = "<<nua<<endl;
1972             }
1973          } // end of for(Int_t nua=0;nua<2;nua++)
1974         } else // to if(dfrEventWeights[t][pW][eW])
1975           {
1976            cout<<"WARNING: dfrEventWeights[t][pW][eW] is NULL in AFAWQC::GOH() !!!!"<<endl;
1977            cout<<"t = "<<t<<endl;
1978            cout<<"pW = "<<pW<<endl;
1979            cout<<"eW = "<<eW<<endl;
1980           }
1981        } // end of for(Int_t eW=0;eW<2;eW++)
1982       } else // to if(dfrParticleWeights[t][pW])
1983         {
1984          cout<<"WARNING: dfrParticleWeights[t][pW] is NULL in AFAWQC::GOH() !!!!"<<endl;
1985          cout<<"t = "<<t<<endl;
1986          cout<<"pW = "<<pW<<endl;
1987         }
1988      } // end of for(Int_t pW=0;pW<(1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights));pW++)
1989     } else // to if(dfrType[t])
1990       {
1991        cout<<"WARNING: dfrType[t] is NULL in AFAWQC::GOH() !!!!"<<endl;
1992        cout<<"t = "<<t<<endl;
1993       }
1994    } // end of for(Int_t t=0;t<2;t++)
1995   } else // to if(diffFlowListResults)
1996     {
1997      cout<<"WARNING: diffFlowListResults is NULL in AFAWQC::GOH() !!!!"<<endl;
1998     }
1999  
2000  TH2D *finalCorrelations2D[2][2][2][4] = {{{{NULL}}}};
2001  TH1D *finalCorrelations1D[2][2][2][2][4] = {{{{{NULL}}}}}; 
2002  TH2D *finalCumulants2D[2][2][2][2][4] = {{{{{NULL}}}}};
2003  TH1D *finalCumulantsPt[2][2][2][2][4] = {{{{{NULL}}}}}; 
2004  TH1D *finalCumulantsEta[2][2][2][2][4] = {{{{{NULL}}}}}; 
2005  TH2D *finalCorrections2D[2][2][2][4] = {{{{NULL}}}};
2006  TH1D *finalCorrections1D[2][2][2][2][4] = {{{{{NULL}}}}}; 
2007  TH2D *finalCovariances2D[2][2][2][4] = {{{{NULL}}}};
2008  TH1D *finalCovariances1D[2][2][2][2][4] = {{{{{NULL}}}}}; 
2009  TH2D *finalFlow2D[2][2][2][2][4] = {{{{{NULL}}}}};
2010  TH1D *finalFlowPt[2][2][2][2][4] = {{{{{NULL}}}}}; 
2011  TH1D *finalFlowEta[2][2][2][2][4] = {{{{{NULL}}}}}; 
2012  TH2D *nonEmptyBins2D[2] = {NULL};
2013  TH1D *nonEmptyBins1D[2][2] = {{NULL}};
2014  
2015  for(Int_t t=0;t<2;t++)
2016  {
2017   // 2D:
2018   nonEmptyBins2D[t] = dynamic_cast<TH2D*>(dfrType[t]->FindObject(Form("%s, (p_{t},#eta)",typeFlag[t].Data())));
2019   if(nonEmptyBins2D[t])
2020   {
2021    this->SetNonEmptyBins2D(nonEmptyBins2D[t],t);
2022   } else
2023     { 
2024      cout<<"WARNING: nonEmptyBins2D[t] is NULL in AFAWQC::GOH() !!!!"<<endl;
2025      cout<<"t = "<<t<<endl;
2026     }
2027   // 1D:
2028   for(Int_t pe=0;pe<2;pe++)
2029   {
2030    nonEmptyBins1D[t][pe] = dynamic_cast<TH1D*>(dfrType[t]->FindObject(Form("%s, %s",typeFlag[t].Data(),ptEtaFlag[pe].Data())));
2031    if(nonEmptyBins1D[t][pe])
2032    {
2033     this->SetNonEmptyBins1D(nonEmptyBins1D[t][pe],t,pe);
2034    } else
2035      { 
2036       cout<<"WARNING: nonEmptyBins1D[t][pe] is NULL in AFAWQC::GOH() !!!!"<<endl;
2037       cout<<"t = "<<t<<endl;
2038       cout<<"pe = "<<pe<<endl;
2039      }
2040   }  
2041   
2042   for(Int_t pW=0;pW<(1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights));pW++)
2043   {
2044    for(Int_t eW=0;eW<2;eW++)
2045    {
2046     // correlations:
2047     if(diffFlowFinalCorrelations[t][pW][eW])
2048     {
2049      for(Int_t correlationIndex=0;correlationIndex<4;correlationIndex++)
2050      {
2051       // 2D:
2052       finalCorrelations2D[t][pW][eW][correlationIndex] = dynamic_cast<TH2D*>(diffFlowFinalCorrelations[t][pW][eW]->FindObject(Form("%s, (p_{t},#eta)",correlationName[correlationIndex].Data())));
2053       if(finalCorrelations2D[t][pW][eW][correlationIndex])
2054       {
2055        this->SetFinalCorrelations2D(finalCorrelations2D[t][pW][eW][correlationIndex],t,pW,eW,correlationIndex);
2056       } else
2057         {
2058          cout<<"WARNING: finalCorrelations2D[t][pW][eW][correlationIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2059          cout<<"t = "<<t<<endl;
2060          cout<<"pW = "<<pW<<endl;
2061          cout<<"eW = "<<eW<<endl;
2062          cout<<"ci = "<<correlationIndex<<endl;
2063         }
2064        // 1D
2065        for(Int_t pe=0;pe<2;pe++)
2066        {
2067         if(pe==0)
2068         {
2069          finalCorrelations1D[t][pW][eW][pe][correlationIndex] = dynamic_cast<TH1D*>(diffFlowFinalCorrelations[t][pW][eW]->FindObject(Form("%s, (p_{t})",correlationName[correlationIndex].Data())));
2070         }
2071         if(pe==1)
2072         {
2073          finalCorrelations1D[t][pW][eW][pe][correlationIndex] = dynamic_cast<TH1D*>(diffFlowFinalCorrelations[t][pW][eW]->FindObject(Form("%s, (#eta)",correlationName[correlationIndex].Data())));
2074         }          
2075         if(finalCorrelations1D[t][pW][eW][pe][correlationIndex])
2076         {
2077          this->SetFinalCorrelations1D(finalCorrelations1D[t][pW][eW][pe][correlationIndex],t,pW,eW,pe,correlationIndex);
2078         } else
2079           {
2080            cout<<"WARNING: finalCorrelations1D[t][pW][eW][pe][correlationIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2081            cout<<"t = "<<t<<endl;
2082            cout<<"pW = "<<pW<<endl;
2083            cout<<"eW = "<<eW<<endl;
2084            cout<<"pe = "<<pe<<endl;
2085            cout<<"ci = "<<correlationIndex<<endl;
2086           }   
2087        } // end of for(Int_t pe=0;pe<2;pe++)  
2088      } // end of for(Int_t correlationIndex=0;correlationIndex<4;correlationIndex++)
2089     } else // to if(diffFlowFinalCorrelations[t][pW][eW])
2090       { 
2091        cout<<"WARNING: diffFlowFinalCorrelations[t][pW] is NULL in AFAWQC::GOH() !!!!"<<endl;
2092        cout<<"t = "<<t<<endl;
2093        cout<<"pW = "<<pW<<endl;
2094        cout<<"eW = "<<eW<<endl;
2095       }
2096     // corrections:
2097     if(diffFlowFinalCorrections[t][pW][eW])
2098     {
2099      for(Int_t correctionIndex=0;correctionIndex<4;correctionIndex++)
2100      {
2101       // 2D:
2102       finalCorrections2D[t][pW][eW][correctionIndex] = dynamic_cast<TH2D*>(diffFlowFinalCorrections[t][pW][eW]->FindObject(Form("%s, (p_{t},#eta)",correctionName[correctionIndex].Data())));
2103       if(finalCorrections2D[t][pW][eW][correctionIndex])
2104       {
2105        this->SetFinalCorrections2D(finalCorrections2D[t][pW][eW][correctionIndex],t,pW,eW,correctionIndex);
2106       } else
2107         {
2108          cout<<"WARNING: finalCorrections2D[t][pW][eW][correctionIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2109          cout<<"t = "<<t<<endl;
2110          cout<<"pW = "<<pW<<endl;
2111          cout<<"eW = "<<eW<<endl;
2112          cout<<"ci = "<<correctionIndex<<endl;
2113         }
2114       // 1D
2115       for(Int_t pe=0;pe<2;pe++)
2116       {
2117        if(pe==0)
2118        {
2119         finalCorrections1D[t][pW][eW][pe][correctionIndex] = dynamic_cast<TH1D*>(diffFlowFinalCorrections[t][pW][eW]->FindObject(Form("%s, (p_{t})",correctionName[correctionIndex].Data())));
2120        }
2121        if(pe==1)
2122        {
2123         finalCorrections1D[t][pW][eW][pe][correctionIndex] = dynamic_cast<TH1D*>(diffFlowFinalCorrections[t][pW][eW]->FindObject(Form("%s, (#eta)",correctionName[correctionIndex].Data())));
2124        }          
2125        if(finalCorrections1D[t][pW][eW][pe][correctionIndex])
2126        {
2127         this->SetFinalCorrections1D(finalCorrections1D[t][pW][eW][pe][correctionIndex],t,pW,eW,pe,correctionIndex);
2128        } else
2129          {
2130           cout<<"WARNING: finalCorrections1D[t][pW][eW][pe][correctionIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2131           cout<<"t = "<<t<<endl;
2132           cout<<"pW = "<<pW<<endl;
2133           cout<<"eW = "<<eW<<endl;
2134           cout<<"pe = "<<pe<<endl;
2135           cout<<"ci = "<<correctionIndex<<endl;
2136          }   
2137       } // end of for(Int_t pe=0;pe<2;pe++)  
2138      } // end of for(Int_t correctionIndex=0;correctionIndex<4;correctionIndex++)
2139     } else // to if(diffFlowFinalCorrections[t][pW][eW])
2140       { 
2141        cout<<"WARNING: diffFlowFinalCorrections[t][pW][eW] is NULL in AFAWQC::GOH() !!!!"<<endl;
2142        cout<<"t = "<<t<<endl;
2143        cout<<"pW = "<<pW<<endl;
2144        cout<<"eW = "<<eW<<endl;
2145       }     
2146     // covariances:
2147     if(diffFlowFinalCovariances[t][pW][eW])
2148     {
2149      for(Int_t covarianceIndex=0;covarianceIndex<4;covarianceIndex++)
2150      {
2151       // 2D:
2152       finalCovariances2D[t][pW][eW][covarianceIndex] = dynamic_cast<TH2D*>(diffFlowFinalCovariances[t][pW][eW]->FindObject(Form("%s, (p_{t},#eta)",covarianceName[covarianceIndex].Data())));
2153       if(finalCovariances2D[t][pW][eW][covarianceIndex])
2154       {
2155        this->SetFinalCovariances2D(finalCovariances2D[t][pW][eW][covarianceIndex],t,pW,eW,covarianceIndex);
2156       } else
2157         {
2158          cout<<"WARNING: finalCovariances2D[t][pW][eW][nua][covarianceIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2159          cout<<"t = "<<t<<endl;
2160          cout<<"pW = "<<pW<<endl;
2161          cout<<"eW = "<<eW<<endl;
2162          cout<<"ci = "<<covarianceIndex<<endl;
2163         }
2164       // 1D:
2165       for(Int_t pe=0;pe<2;pe++)
2166       {
2167        if(pe==0)
2168        {
2169         finalCovariances1D[t][pW][eW][pe][covarianceIndex] = dynamic_cast<TH1D*>(diffFlowFinalCovariances[t][pW][eW]->FindObject(Form("%s, (p_{t})",covarianceName[covarianceIndex].Data())));
2170        }
2171        if(pe==1)
2172        {
2173         finalCovariances1D[t][pW][eW][pe][covarianceIndex] = dynamic_cast<TH1D*>(diffFlowFinalCovariances[t][pW][eW]->FindObject(Form("%s, (#eta)",covarianceName[covarianceIndex].Data())));
2174        }          
2175        if(finalCovariances1D[t][pW][eW][pe][covarianceIndex])
2176        {
2177         this->SetFinalCovariances1D(finalCovariances1D[t][pW][eW][pe][covarianceIndex],t,pW,eW,pe,covarianceIndex);
2178        } else
2179          {
2180           cout<<"WARNING: finalCovariances1D[t][pW][eW][pe][covarianceIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2181           cout<<"t = "<<t<<endl;
2182           cout<<"pW = "<<pW<<endl;
2183           cout<<"eW = "<<eW<<endl;
2184           cout<<"pe = "<<pe<<endl;
2185           cout<<"ci = "<<covarianceIndex<<endl;
2186          }   
2187       } // end of for(Int_t pe=0;pe<2;pe++)  
2188      } // end of for(Int_t covarianceIndex=0;covarianceIndex<4;covarianceIndex++)
2189     } else // to if(diffFlowFinalCovariances[t][pW][eW])
2190       { 
2191        cout<<"WARNING: diffFlowFinalCovariances[t][pW][eW] is NULL in AFAWQC::GOH() !!!!"<<endl;
2192        cout<<"t = "<<t<<endl;
2193        cout<<"pW = "<<pW<<endl;
2194        cout<<"eW = "<<eW<<endl;
2195       }        
2196        
2197     for(Int_t nua=0;nua<2;nua++)
2198     {  
2199      // cumulants:
2200      if(diffFlowFinalCumulants[t][pW][eW][nua])
2201      {
2202       for(Int_t cumulantIndex=0;cumulantIndex<4;cumulantIndex++)
2203       {
2204        // 2D:
2205        finalCumulants2D[t][pW][eW][nua][cumulantIndex] = dynamic_cast<TH2D*>(diffFlowFinalCumulants[t][pW][eW][nua]->FindObject(Form("%s, (p_{t},#eta)",cumulantName[cumulantIndex].Data())));
2206        if(finalCumulants2D[t][pW][eW][nua][cumulantIndex])
2207        {
2208         this->SetFinalCumulants2D(finalCumulants2D[t][pW][eW][nua][cumulantIndex],t,pW,eW,nua,cumulantIndex);
2209        } else
2210          {
2211           cout<<"WARNING: finalCumulants2D[t][pW][eW][nua][cumulantIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2212           cout<<"t = "<<t<<endl;
2213           cout<<"pW = "<<pW<<endl;
2214           cout<<"eW = "<<eW<<endl;
2215           cout<<"nua = "<<nua<<endl;
2216           cout<<"ci = "<<cumulantIndex<<endl;
2217          }
2218        // pt:  
2219        finalCumulantsPt[t][pW][eW][nua][cumulantIndex] = dynamic_cast<TH1D*>(diffFlowFinalCumulants[t][pW][eW][nua]->FindObject(Form("%s, (p_{t})",cumulantName[cumulantIndex].Data())));
2220        if(finalCumulantsPt[t][pW][eW][nua][cumulantIndex])
2221        {
2222         this->SetFinalCumulantsPt(finalCumulantsPt[t][pW][eW][nua][cumulantIndex],t,pW,eW,nua,cumulantIndex);
2223        } else
2224          {
2225           cout<<"WARNING: finalCumulantsPt[t][pW][eW][nua][cumulantIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2226           cout<<"t = "<<t<<endl;
2227           cout<<"pW = "<<pW<<endl;
2228           cout<<"eW = "<<eW<<endl;
2229           cout<<"nua = "<<nua<<endl;
2230           cout<<"ci = "<<cumulantIndex<<endl;
2231          } 
2232        // eta:
2233        finalCumulantsEta[t][pW][eW][nua][cumulantIndex] = dynamic_cast<TH1D*>(diffFlowFinalCumulants[t][pW][eW][nua]->FindObject(Form("%s, (#eta)",cumulantName[cumulantIndex].Data())));
2234        if(finalCumulantsEta[t][pW][eW][nua][cumulantIndex])
2235        {
2236         this->SetFinalCumulantsEta(finalCumulantsEta[t][pW][eW][nua][cumulantIndex],t,pW,eW,nua,cumulantIndex);
2237        } else
2238          {
2239           cout<<"WARNING: finalCumulantsEta[t][pW][eW][nua][cumulantIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2240           cout<<"t = "<<t<<endl;
2241           cout<<"pW = "<<pW<<endl;
2242           cout<<"eW = "<<eW<<endl;
2243           cout<<"nua = "<<nua<<endl;
2244           cout<<"ci = "<<cumulantIndex<<endl;
2245          }   
2246       } // end of for(Int_t cumulantIndex=0;cumulantIndex<4;cumulantIndex++)
2247      } else // to if(diffFlowFinalCumulants[t][pW][eW][nua])
2248        { 
2249         cout<<"WARNING: diffFlowFinalCumulants[t][pW][eW][nua] is NULL in AFAWQC::GOH() !!!!"<<endl;
2250         cout<<"t = "<<t<<endl;
2251         cout<<"pW = "<<pW<<endl;
2252         cout<<"eW = "<<eW<<endl;
2253         cout<<"nua = "<<nua<<endl;
2254        }      
2255      // flow:
2256      if(diffFlowFinalFlow[t][pW][eW][nua])
2257      {
2258       for(Int_t flowIndex=0;flowIndex<4;flowIndex++)
2259       {
2260        // 2D:
2261        finalFlow2D[t][pW][eW][nua][flowIndex] = dynamic_cast<TH2D*>(diffFlowFinalFlow[t][pW][eW][nua]->FindObject(Form("%s, (p_{t},#eta)",flowName[flowIndex].Data())));
2262        if(finalFlow2D[t][pW][eW][nua][flowIndex])
2263        {
2264         this->SetFinalFlow2D(finalFlow2D[t][pW][eW][nua][flowIndex],t,pW,eW,nua,flowIndex);
2265        } else
2266          {
2267           cout<<"WARNING: finalFlow2D[t][pW][eW][nua][flowIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2268           cout<<"t = "<<t<<endl;
2269           cout<<"pW = "<<pW<<endl;
2270           cout<<"eW = "<<eW<<endl;
2271           cout<<"nua = "<<nua<<endl;
2272           cout<<"ci = "<<flowIndex<<endl;
2273          }
2274        // pt:
2275        finalFlowPt[t][pW][eW][nua][flowIndex] = dynamic_cast<TH1D*>(diffFlowFinalFlow[t][pW][eW][nua]->FindObject(Form("%s, (p_{t})",flowName[flowIndex].Data())));
2276        if(finalFlowPt[t][pW][eW][nua][flowIndex])
2277        {
2278         this->SetFinalFlowPt(finalFlowPt[t][pW][eW][nua][flowIndex],t,pW,eW,nua,flowIndex);
2279        } else
2280          {
2281           cout<<"WARNING: finalFlow1D[t][pW][nua][flowIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2282           cout<<"t = "<<t<<endl;
2283           cout<<"pW = "<<pW<<endl;
2284           cout<<"eW = "<<eW<<endl;
2285           cout<<"nua = "<<nua<<endl;
2286           cout<<"ci = "<<flowIndex<<endl;
2287          }   
2288        // eta: 
2289        finalFlowEta[t][pW][eW][nua][flowIndex] = dynamic_cast<TH1D*>(diffFlowFinalFlow[t][pW][eW][nua]->FindObject(Form("%s, (#eta)",flowName[flowIndex].Data())));          
2290        if(finalFlowEta[t][pW][eW][nua][flowIndex])
2291        {
2292         this->SetFinalFlowEta(finalFlowEta[t][pW][eW][nua][flowIndex],t,pW,eW,nua,flowIndex);
2293        } else
2294          {
2295           cout<<"WARNING: finalFlow1D[t][pW][nua][flowIndex] is NULL in AFAWQC::GOH() !!!!"<<endl;
2296           cout<<"t = "<<t<<endl;
2297           cout<<"pW = "<<pW<<endl;
2298           cout<<"eW = "<<eW<<endl;
2299           cout<<"nua = "<<nua<<endl;
2300           cout<<"ci = "<<flowIndex<<endl;
2301          }   
2302       } // end of for(Int_t flowIndex=0;flowIndex<4;flowIndex++)
2303      } else // to if(diffFlowFinalFlow[t][pW][eW][nua])
2304        { 
2305         cout<<"WARNING: diffFlowFinalFlow[t][pW][eW][nua] is NULL in AFAWQC::GOH() !!!!"<<endl;
2306         cout<<"t = "<<t<<endl;
2307         cout<<"pW = "<<pW<<endl;
2308         cout<<"eW = "<<eW<<endl;
2309         cout<<"nua = "<<nua<<endl;
2310        }               
2311     } // end of for(Int_t nua=0;nua<2;nua++)
2312    } // end of for(Int_t eW=0;eW<2;eW++) 
2313   } // end of for(Int_t pW=0;pW<(1+(Int_t)(bUsePhiWeights||bUsePtWeights||bUseEtaWeights));pW++)
2314  } // end of for(Int_t t=0;t<2;t++)
2315  
2316   // x.) nested loops:
2317   TList *nestedLoopsList = dynamic_cast<TList*>(outputListHistos->FindObject("Nested Loops"));
2318   if(nestedLoopsList) this->SetNestedLoopsList(nestedLoopsList);
2319   TProfile *evaluateNestedLoops = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fEvaluateNestedLoops"));
2320   Bool_t bEvaluateNestedLoopsForIntFlow = kFALSE;
2321   Bool_t bEvaluateNestedLoopsForDiffFlow = kFALSE;
2322   if(evaluateNestedLoops)
2323   {
2324    this->SetEvaluateNestedLoops(evaluateNestedLoops);
2325    bEvaluateNestedLoopsForIntFlow = (Int_t)evaluateNestedLoops->GetBinContent(1);
2326    bEvaluateNestedLoopsForDiffFlow = (Int_t)evaluateNestedLoops->GetBinContent(2);
2327   }
2328   
2329   if(bEvaluateNestedLoopsForIntFlow)
2330   {
2331    TProfile *directCorrelations = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrelation"));
2332    if(directCorrelations) this->SetDirectCorrelations(directCorrelations);
2333    TProfile *directCorrectionsCos = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrectionsCos"));
2334    if(directCorrectionsCos) this->SetDirectCorrectionsCos(directCorrectionsCos);
2335    TProfile *directCorrectionsSin = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrectionsSin"));
2336    if(directCorrectionsSin) this->SetDirectCorrectionsSin(directCorrectionsSin);
2337    if(bUsePhiWeights||bUsePtWeights||bUseEtaWeights)
2338    {
2339     TProfile *directCorrelationsW = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrelationW"));
2340     if(directCorrelationsW) this->SetDirectCorrelationsW(directCorrelationsW);
2341     TProfile *directCorrectionsCosW = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrectionsCosW"));
2342     if(directCorrectionsCosW) this->SetDirectCorrectionsCosW(directCorrectionsCosW);
2343     TProfile *directCorrectionsSinW = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrectionsSinW")); 
2344     if(directCorrectionsSinW) this->SetDirectCorrectionsSinW(directCorrectionsSinW);
2345    }
2346   }
2347   
2348   if(bEvaluateNestedLoopsForDiffFlow)
2349   {
2350    TProfile *directCorrelationsDiffFlow = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrelationsDiffFlow"));
2351    if(directCorrelationsDiffFlow) this->SetDirectCorrelationsDiffFlow(directCorrelationsDiffFlow);
2352    TProfile *directCorrectionsDiffFlowCos = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrectionsDiffFlowCos"));
2353    if(directCorrectionsDiffFlowCos) this->SetDirectCorrectionsDiffFlowCos(directCorrectionsDiffFlowCos);
2354    TProfile *directCorrectionsDiffFlowSin = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrectionsDiffFlowSin"));
2355    if(directCorrectionsDiffFlowSin) this->SetDirectCorrectionsDiffFlowSin(directCorrectionsDiffFlowSin);
2356    if(bUsePhiWeights||bUsePtWeights||bUseEtaWeights)
2357    {
2358     TProfile *directCorrelationsDiffFlowW = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrelationsDiffFlowW")); 
2359     if(directCorrelationsDiffFlowW) this->SetDirectCorrelationsDiffFlowW(directCorrelationsDiffFlowW);
2360     TProfile *directCorrectionsDiffFlowCosW = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrectionsDiffFlowCosW"));
2361     if(directCorrectionsDiffFlowCosW) this->SetDirectCorrectionsDiffFlowCosW(directCorrectionsDiffFlowCosW);
2362     TProfile *directCorrectionsDiffFlowSinW = dynamic_cast<TProfile*>(nestedLoopsList->FindObject("fDirectCorrectionsDiffFlowSinW"));
2363     if(directCorrectionsDiffFlowSinW) this->SetDirectCorrectionsDiffFlowSinW(directCorrectionsDiffFlowSinW);
2364    }
2365   }
2366   
2367  } 
2368 }
2369
2370
2371 //================================================================================================================================
2372
2373
2374 TProfile* AliFlowAnalysisWithQCumulants::MakePtProjection(TProfile2D *profilePtEta) const
2375 {
2376  // project 2D profile onto pt axis to get 1D profile
2377  
2378  Int_t nBinsPt   = profilePtEta->GetNbinsX();
2379  Double_t dPtMin = (profilePtEta->GetXaxis())->GetXmin();
2380  Double_t dPtMax = (profilePtEta->GetXaxis())->GetXmax();
2381  
2382  Int_t nBinsEta   = profilePtEta->GetNbinsY();
2383  
2384  TProfile *profilePt = new TProfile("","",nBinsPt,dPtMin,dPtMax); 
2385  
2386  for(Int_t p=1;p<=nBinsPt;p++)
2387  {
2388   Double_t contentPt = 0.;
2389   Double_t entryPt = 0.;
2390   Double_t spreadPt = 0.;
2391   Double_t sum1 = 0.;
2392   Double_t sum2 = 0.;
2393   Double_t sum3 = 0.;
2394   for(Int_t e=1;e<=nBinsEta;e++)
2395   {
2396    contentPt += (profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)))
2397               * (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
2398    entryPt   += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
2399    
2400    sum1 += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)))
2401          * (pow(profilePtEta->GetBinError(profilePtEta->GetBin(p,e)),2.)
2402             + pow(profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)),2.)); 
2403    sum2 += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
2404    sum3 += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)))
2405          * (profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)));            
2406   }
2407   if(sum2>0. && sum1/sum2-pow(sum3/sum2,2.) > 0.)
2408   {
2409    spreadPt = pow(sum1/sum2-pow(sum3/sum2,2.),0.5);
2410   }
2411   profilePt->SetBinContent(p,contentPt);
2412   profilePt->SetBinEntries(p,entryPt);
2413   {
2414    profilePt->SetBinError(p,spreadPt);
2415   }
2416   
2417  }
2418  
2419  return profilePt;
2420  
2421 } // end of TProfile* AliFlowAnalysisWithQCumulants::MakePtProjection(TProfile2D *profilePtEta)
2422
2423
2424 //================================================================================================================================
2425
2426
2427 TProfile* AliFlowAnalysisWithQCumulants::MakeEtaProjection(TProfile2D *profilePtEta) const
2428 {
2429  // project 2D profile onto eta axis to get 1D profile
2430  
2431  Int_t nBinsEta   = profilePtEta->GetNbinsY();
2432  Double_t dEtaMin = (profilePtEta->GetYaxis())->GetXmin();
2433  Double_t dEtaMax = (profilePtEta->GetYaxis())->GetXmax();
2434  
2435  Int_t nBinsPt = profilePtEta->GetNbinsX();
2436  
2437  TProfile *profileEta = new TProfile("","",nBinsEta,dEtaMin,dEtaMax); 
2438  
2439  for(Int_t e=1;e<=nBinsEta;e++)
2440  {
2441   Double_t contentEta = 0.;
2442   Double_t entryEta = 0.;
2443   for(Int_t p=1;p<=nBinsPt;p++)
2444   {
2445    contentEta += (profilePtEta->GetBinContent(profilePtEta->GetBin(p,e)))
2446               * (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
2447    entryEta   += (profilePtEta->GetBinEntries(profilePtEta->GetBin(p,e)));
2448   }
2449   profileEta->SetBinContent(e,contentEta);
2450   profileEta->SetBinEntries(e,entryEta);
2451  }
2452  
2453  return profileEta;
2454  
2455 } // end of TProfile* AliFlowAnalysisWithQCumulants::MakeEtaProjection(TProfile2D *profilePtEta)
2456
2457
2458 //================================================================================================================================
2459
2460
2461 void AliFlowAnalysisWithQCumulants::CalculateFinalCorrectionsForNonUniformAcceptanceForCumulantsForIntFlow(Bool_t useParticleWeights, TString eventWeights)
2462 {
2463  // calculate final corrections for non-uniform acceptance for QC{2}, QC{4}, QC{6} and QC{8}
2464  
2465  // corrections for non-uniform acceptance (NUA) are stored in histogram fCorrectionsForNUA,
2466  // binning of fCorrectionsForNUA is organized as follows:
2467  //
2468  // 1st bin: correction to QC{2}
2469  // 2nd bin: correction to QC{4}
2470  // 3rd bin: correction to QC{6}
2471  // 4th bin: correction to QC{8}
2472  
2473  // shortcuts flags:
2474  Int_t pW = (Int_t)(useParticleWeights);
2475  
2476  Int_t eW = -1;
2477  
2478  if(eventWeights == "exact")
2479  {
2480   eW = 0;
2481  }
2482
2483  for(Int_t sc=0;sc<2;sc++) // sin or cos terms flag
2484  {
2485   if(!(fQCorrelations[pW][eW] && fQCorrections[pW][eW][sc] && fCorrections[pW][eW]))
2486   {
2487    cout<<"WARNING: fQCorrelations[pW][eW] && fQCorrections[pW][eW][sc] && fCorrections[pW][eW] is NULL in AFAWQC::CFCFNUAFIF() !!!!"<<endl;
2488    cout<<"pW = "<<pW<<endl;
2489    cout<<"eW = "<<eW<<endl;
2490    cout<<"sc = "<<sc<<endl;
2491    exit(0);
2492   }
2493  }  
2494
2495  // measured 2-, 4-, 6- and 8-particle azimuthal correlations (biased with non-uniform acceptance!):
2496  Double_t two = fQCorrelations[pW][eW]->GetBinContent(1); // <<2>>
2497  //Double_t four = fQCorrelations[pW][eW]->GetBinContent(11); // <<4>>
2498  //Double_t six = fQCorrelations[pW][eW]->GetBinContent(24); // <<6>>
2499  //Double_t eight = fQCorrelations[pW][eW]->GetBinContent(31); // <<8>>
2500  
2501  // correction terms to QC{2}:
2502  // <<cos(n*phi1)>>^2
2503  Double_t two1stTerm = pow(fQCorrections[pW][eW][1]->GetBinContent(1),2); 
2504  // <<sin(n*phi1)>>^2
2505  Double_t two2ndTerm = pow(fQCorrections[pW][eW][0]->GetBinContent(1),2); 
2506  // final corrections for non-uniform acceptance to QC{2}:
2507  Double_t correctionQC2 = -1.*two1stTerm-1.*two2ndTerm;
2508  fCorrections[pW][eW]->SetBinContent(1,correctionQC2); 
2509  
2510  // correction terms to QC{4}:
2511  // <<cos(n*phi1)>> <<cos(n*(phi1-phi2-phi3))>>
2512  Double_t four1stTerm = fQCorrections[pW][eW][1]->GetBinContent(1)*fQCorrections[pW][eW][1]->GetBinContent(3);  
2513  // <<sin(n*phi1)>> <<sin(n*(phi1-phi2-phi3))>>
2514  Double_t four2ndTerm = fQCorrections[pW][eW][0]->GetBinContent(1)*fQCorrections[pW][eW][0]->GetBinContent(3);  
2515  // <<cos(n*(phi1+phi2))>>^2
2516  Double_t four3rdTerm = pow(fQCorrections[pW][eW][1]->GetBinContent(2),2); 
2517  // <<sin(n*(phi1+phi2))>>^2
2518  Double_t four4thTerm = pow(fQCorrections[pW][eW][0]->GetBinContent(2),2); 
2519  // <<cos(n*(phi1+phi2))>> (<<cos(n*phi1)>>^2 - <<sin(n*phi1)>>^2)
2520  Double_t four5thTerm = fQCorrections[pW][eW][1]->GetBinContent(2)
2521                       * (pow(fQCorrections[pW][eW][1]->GetBinContent(1),2)-pow(fQCorrections[pW][eW][0]->GetBinContent(1),2));
2522  // <<sin(n*(phi1+phi2))>> <<cos(n*phi1)>> <<sin(n*phi1)>>
2523  Double_t four6thTerm = fQCorrections[pW][eW][0]->GetBinContent(2)
2524                       * fQCorrections[pW][eW][1]->GetBinContent(1)
2525                       * fQCorrections[pW][eW][0]->GetBinContent(1);         
2526  // <<cos(n*(phi1-phi2))>> (<<cos(n*phi1)>>^2 + <<sin(n*phi1)>>^2)
2527  Double_t four7thTerm = two*(pow(fQCorrections[pW][eW][1]->GetBinContent(1),2)+pow(fQCorrections[pW][eW][0]->GetBinContent(1),2));  
2528  // (<<cos(n*phi1)>>^2 + <<sin(n*phi1)>>^2)^2
2529  Double_t four8thTerm = pow(pow(fQCorrections[pW][eW][1]->GetBinContent(1),2)+pow(fQCorrections[pW][eW][0]->GetBinContent(1),2),2);      
2530  // final correction to QC{4}:
2531  Double_t correctionQC4 = -4.*four1stTerm+4.*four2ndTerm-four3rdTerm-four4thTerm
2532                         + 4.*four5thTerm+8.*four6thTerm+8.*four7thTerm-6.*four8thTerm;                            
2533  fCorrections[pW][eW]->SetBinContent(2,correctionQC4);   
2534
2535  // ... to be improved (continued for 6th and 8th order)                                                    
2536
2537 } // end of AliFlowAnalysisWithQCumulants::CalculateFinalCorrectionsForNonUniformAcceptanceForCumulantsForIntFlow(Bool_t useParticleWeights, TString eventWeights)
2538
2539
2540 //================================================================================================================================
2541
2542
2543 void AliFlowAnalysisWithQCumulants::CalculateFinalCorrectionsForNonUniformAcceptanceForDifferentialFlow(Bool_t useParticleWeights,TString type)
2544 {
2545  
2546  useParticleWeights=kFALSE;
2547  type="ac";
2548   
2549  // calculate final corrections due to non-uniform acceptance of the detector to reduced multi-particle correlations
2550  /*
2551  if(!(useParticleWeights))
2552  {
2553   if(type == "POI")
2554   { 
2555    // **** corrections for non-uniform acceptance for 2nd order QC' for POI's ****
2556    
2557    // 1st term: <<cos(n*psi)>><<cos(n*phi)>>:
2558    if(fCorrectionsCosP1nPsiPtEtaPOI && fQCorrectionsCos)
2559    {
2560     // pt,eta: 
2561     if(f2pFinalCorrectionsForNUAPtEtaPOI) f2pFinalCorrectionsForNUAPtEtaPOI->Reset(); // to be improved
2562     TH2D *correctionPtEta1stTerm = new TH2D(*(fCorrectionsCosP1nPsiPtEtaPOI->ProjectionXY("","e")));
2563     correctionPtEta1stTerm->Scale(fQCorrectionsCos->GetBinContent(1)); // to be improved: are errors propagated correctly here?   
2564     if(f2pFinalCorrectionsForNUAPtEtaPOI) f2pFinalCorrectionsForNUAPtEtaPOI->Add(correctionPtEta1stTerm); // to be improved (if condition goes somewhere else)
2565     delete correctionPtEta1stTerm;
2566     // pt:
2567     if(f2pFinalCorrectionsForNUAPtPOI) f2pFinalCorrectionsForNUAPtPOI->Reset(); // to be improved
2568     TH1D *correctionPt1stTerm = new TH1D(*((this->MakePtProjection(fCorrectionsCosP1nPsiPtEtaPOI))->ProjectionX("","e"))); // to be improved: are errors propagated correctly here? 
2569     correctionPt1stTerm->Scale(fQCorrectionsCos->GetBinContent(1)); // to be improved: are errors propagated correctly here? 
2570     if(f2pFinalCorrectionsForNUAPtPOI) f2pFinalCorrectionsForNUAPtPOI->Add(correctionPt1stTerm); // to be improved (if condition goes somewhere else)
2571     delete correctionPt1stTerm;
2572     // eta:
2573     if(f2pFinalCorrectionsForNUAEtaPOI) f2pFinalCorrectionsForNUAEtaPOI->Reset(); // to be improved    
2574     TH1D *correctionEta1stTerm = new TH1D(*((this->MakeEtaProjection(fCorrectionsCosP1nPsiPtEtaPOI))->ProjectionX("","e"))); // to be improved: are errors propagated correctly here? 
2575     correctionEta1stTerm->Scale(fQCorrectionsCos->GetBinContent(1)); // to be improved: are errors propagated correctly here? 
2576     if(f2pFinalCorrectionsForNUAEtaPOI) f2pFinalCorrectionsForNUAEtaPOI->Add(correctionEta1stTerm); // to be improved (if condition goes somewhere else)
2577     delete correctionEta1stTerm;    
2578    } else
2579      { 
2580       cout<<"WARNING: (fCorrectionsCosP1nPsiPtEtaPOI && fQCorrectionsCos && f2pFinalCorrectionsForNUAPtEtaPOI) is NULL in QC::CFCFNUAFDF() !!!!  "<<endl;
2581       cout<<"         Corrections for non-uniform acceptance for differential flow are not correct."<<endl;
2582      } 
2583      
2584    // 2nd term: <<sin(n*psi)>><<sin(n*phi)>>:  
2585    if(fCorrectionsSinP1nPsiPtEtaPOI && fQCorrectionsSin)
2586    {
2587     // pt,eta:
2588     TH2D *correctionPtEta2ndTerm = new TH2D(*(fCorrectionsSinP1nPsiPtEtaPOI->ProjectionXY("","e")));
2589     correctionPtEta2ndTerm->Scale(fQCorrectionsSin->GetBinContent(1)); // to be improved: are errors propagated correctly here?    
2590     if(f2pFinalCorrectionsForNUAPtEtaPOI) f2pFinalCorrectionsForNUAPtEtaPOI->Add(correctionPtEta2ndTerm); // to be improved (if condition goes somewhere else)
2591     delete correctionPtEta2ndTerm;
2592     // pt:
2593     TH1D *correctionPt2ndTerm = new TH1D(*((this->MakePtProjection(fCorrectionsSinP1nPsiPtEtaPOI))->ProjectionX("","e"))); // to be improved: are errors propagated correctly here? 
2594     correctionPt2ndTerm->Scale(fQCorrectionsSin->GetBinContent(1)); // to be improved: are errors propagated correctly here? 
2595     if(f2pFinalCorrectionsForNUAPtPOI) f2pFinalCorrectionsForNUAPtPOI->Add(correctionPt2ndTerm); // to be improved (if condition goes somewhere else)
2596     delete correctionPt2ndTerm;
2597     // eta:
2598     TH1D *correctionEta2ndTerm = new TH1D(*((this->MakeEtaProjection(fCorrectionsSinP1nPsiPtEtaPOI))->ProjectionX("","e"))); // to be improved: are errors propagated correctly here? 
2599     correctionEta2ndTerm->Scale(fQCorrectionsSin->GetBinContent(1)); // to be improved: are errors propagated correctly here? 
2600     if(f2pFinalCorrectionsForNUAEtaPOI) f2pFinalCorrectionsForNUAEtaPOI->Add(correctionEta2ndTerm); // to be improved (if condition goes somewhere else)
2601     delete correctionEta2ndTerm; 
2602    } else
2603      { 
2604       cout<<"WARNING: (fCorrectionsSinP1nPsiPtEtaPOI && fQCorrectionsSin) is NULL in QC::CFCFNUAFDF() !!!!  "<<endl;
2605       cout<<"         Corrections for non-uniform acceptance for differential flow are not correct."<<endl;
2606      } 
2607   } else if(type == "RP")
2608     {
2609      // **** corrections for non-uniform acceptance for 2nd order QC' for RP's ****
2610    
2611      // 1st term: <<cos(n*psi)>><<cos(n*phi)>>:
2612      if(fCorrectionsCosP1nPsiPtEtaRP && fQCorrectionsCos)
2613      {
2614       // pt,eta: 
2615       if(f2pFinalCorrectionsForNUAPtEtaRP) f2pFinalCorrectionsForNUAPtEtaRP->Reset(); // to be improved
2616       TH2D *correctionPtEta1stTerm = new TH2D(*(fCorrectionsCosP1nPsiPtEtaRP->ProjectionXY("","e")));
2617       correctionPtEta1stTerm->Scale(fQCorrectionsCos->GetBinContent(1)); // to be improved: are errors propagated correctly here?    
2618       if(f2pFinalCorrectionsForNUAPtEtaRP) f2pFinalCorrectionsForNUAPtEtaRP->Add(correctionPtEta1stTerm); // to be improved (if condition goes somewhere else)
2619       delete correctionPtEta1stTerm;
2620       // pt:
2621       if(f2pFinalCorrectionsForNUAPtRP) f2pFinalCorrectionsForNUAPtRP->Reset(); // to be improved
2622       TH1D *correctionPt1stTerm = new TH1D(*((this->MakePtProjection(fCorrectionsCosP1nPsiPtEtaRP))->ProjectionX("","e"))); // to be improved: are errors propagated correctly here? 
2623       correctionPt1stTerm->Scale(fQCorrectionsCos->GetBinContent(1)); // to be improved: are errors propagated correctly here? 
2624       if(f2pFinalCorrectionsForNUAPtRP) f2pFinalCorrectionsForNUAPtRP->Add(correctionPt1stTerm); // to be improved (if condition goes somewhere else)
2625       delete correctionPt1stTerm;
2626       // eta:
2627       if(f2pFinalCorrectionsForNUAEtaRP) f2pFinalCorrectionsForNUAEtaRP->Reset(); // to be improved
2628       TH1D *correctionEta1stTerm = new TH1D(*((this->MakeEtaProjection(fCorrectionsCosP1nPsiPtEtaRP))->ProjectionX("","e"))); // to be improved: are errors propagated correctly here? 
2629       correctionEta1stTerm->Scale(fQCorrectionsCos->GetBinContent(1)); // to be improved: are errors propagated correctly here? 
2630       if(f2pFinalCorrectionsForNUAEtaRP) f2pFinalCorrectionsForNUAEtaRP->Add(correctionEta1stTerm); // to be improved (if condition goes somewhere else)
2631       delete correctionEta1stTerm;    
2632      } else
2633        { 
2634         cout<<"WARNING: (fCorrectionsCosP1nPsiPtEtaRP && fQCorrectionsCos) is NULL in QC::CFCFNUAFDF() !!!!  "<<endl;
2635         cout<<"         Corrections for non-uniform acceptance for differential flow are not correct."<<endl;
2636        } 
2637      // 2nd term: <<sin(n*psi)>><<sin(n*phi)>>:  
2638      if(fCorrectionsSinP1nPsiPtEtaRP && fQCorrectionsSin)
2639      {
2640       // pt,eta: 
2641       TH2D *correctionPtEta2ndTerm = new TH2D(*(fCorrectionsSinP1nPsiPtEtaRP->ProjectionXY("","e")));
2642       correctionPtEta2ndTerm->Scale(fQCorrectionsSin->GetBinContent(1)); // to be improved: are errors propagated correctly here?    
2643       if(f2pFinalCorrectionsForNUAPtEtaRP) f2pFinalCorrectionsForNUAPtEtaRP->Add(correctionPtEta2ndTerm); // to be improved (if condition goes somewhere else)
2644       delete correctionPtEta2ndTerm;
2645       // pt:
2646       TH1D *correctionPt2ndTerm = new TH1D(*((this->MakePtProjection(fCorrectionsSinP1nPsiPtEtaRP))->ProjectionX("","e"))); // to be improved: are errors propagated correctly here? 
2647       correctionPt2ndTerm->Scale(fQCorrectionsSin->GetBinContent(1)); // to be improved: are errors propagated correctly here? 
2648       if(f2pFinalCorrectionsForNUAPtRP) f2pFinalCorrectionsForNUAPtRP->Add(correctionPt2ndTerm); // to be improved (if condition goes somewhere else)
2649       delete correctionPt2ndTerm;
2650       // eta:
2651       TH1D *correctionEta2ndTerm = new TH1D(*((this->MakeEtaProjection(fCorrectionsSinP1nPsiPtEtaRP))->ProjectionX("","e"))); // to be improved: are errors propagated correctly here? 
2652       correctionEta2ndTerm->Scale(fQCorrectionsSin->GetBinContent(1)); // to be improved: are errors propagated correctly here? 
2653       if(f2pFinalCorrectionsForNUAEtaRP) f2pFinalCorrectionsForNUAEtaRP->Add(correctionEta2ndTerm); // to be improved (if condition goes somewhere else)
2654       delete correctionEta2ndTerm; 
2655      } else
2656        { 
2657         cout<<"WARNING: (fCorrectionsSinP1nPsiPtEtaRP && fQCorrectionsSin) is NULL in QC::CFCFNUAFDF() !!!!  "<<endl;
2658         cout<<"         Corrections for non-uniform acceptance for differential flow are not correct."<<endl;
2659        }              
2660     } else // to else if(type == "RP")
2661       {
2662        cout<<"WARNING: Type must be either POI or RP in QC::CFCFNUAFDF() !!!!                           "<<endl;
2663        cout<<"         Corrections for non-uniform acceptance for differential flow were not calculated."<<endl;
2664       }  
2665  } else // to if(!(useParticleWeights))
2666    {
2667     // ...
2668    }
2669  */
2670 } // end of AliFlowAnalysisWithQCumulants::CalculateFinalCorrectionsForNonUniformAcceptanceForDifferentialFlow(Bool_t useParticleWeights,TString type)
2671
2672
2673 //==================================================================================================================================
2674
2675 /*
2676 void AliFlowAnalysisWithQCumulants::CalculateFinalResultsForDifferentialFlow(
2677                                                         TH2D *flowPtEta, TH1D *flowPt, TH1D *flowEta, 
2678                                                         TProfile2D *profile2ndPtEta, TProfile2D *profile4thPtEta, 
2679                                                         TProfile2D *profile6thPtEta, TProfile2D *profile8thPtEta)
2680 {
2681  // calculate and store the final results for integrated flow
2682  
2683  TString *namePtEta = new TString();
2684  TString *type = new TString();
2685  TString *order2nd = new TString();
2686  TString *order4th = new TString();
2687  TString *order6th = new TString();
2688  TString *order8th = new TString(); 
2689  TString *pW = new TString();
2690
2691  if(profile2ndPtEta) *namePtEta = profile2ndPtEta->GetName();
2692  if(namePtEta->Contains("POI")) *type = "POI";
2693  if(namePtEta->Contains("RP")) *type  = "RP";
2694  if(namePtEta->Contains("W")) *pW      = "W";
2695  if(namePtEta->Contains("2")) *order2nd  = "2";
2696  if(profile4thPtEta) *namePtEta = profile4thPtEta->GetName();
2697  if(namePtEta->Contains("4")) *order4th = "4";
2698
2699  if(profile6thPtEta) *namePtEta = profile6thPtEta->GetName();
2700  if(namePtEta->Contains("6")) *order6th = "6";
2701  
2702  if(profile8thPtEta) *namePtEta = profile8thPtEta->GetName();
2703  if(namePtEta->Contains("8")) *order8th = "8";
2704
2705  TProfile *profile2ndPt = NULL;
2706  TProfile *profile4thPt = NULL;
2707  TProfile *profile6thPt = NULL;
2708  TProfile *profile8thPt = NULL;
2709
2710  TProfile *profile2ndEta = NULL;
2711  TProfile *profile4thEta = NULL;
2712  TProfile *profile6thEta = NULL;
2713  TProfile *profile8thEta = NULL;
2714   
2715  if(*order2nd == "2")
2716  {
2717   profile2ndPt  = new TProfile(*(this->MakePtProjection(profile2ndPtEta))); 
2718   profile2ndEta = new TProfile(*(this->MakeEtaProjection(profile2ndPtEta))); 
2719   if(*order4th == "4")
2720   {
2721    profile4thPt  = new TProfile(*(this->MakePtProjection(profile4thPtEta))); 
2722    profile4thEta = new TProfile(*(this->MakeEtaProjection(profile4thPtEta))); 
2723    if(*order6th == "6")
2724    {
2725     profile6thPt  = new TProfile(*(this->MakePtProjection(profile6thPtEta))); 
2726     profile6thEta = new TProfile(*(this->MakeEtaProjection(profile6thPtEta))); 
2727     if(*order8th == "8")
2728     {
2729      profile8thPt  = new TProfile(*(this->MakePtProjection(profile8thPtEta))); 
2730      profile8thEta = new TProfile(*(this->MakeEtaProjection(profile8thPtEta))); 
2731     }     
2732    }    
2733   } 
2734  }
2735  
2736  Int_t nBinsPt  = profile2ndPt->GetNbinsX();
2737  Int_t nBinsEta = profile2ndEta->GetNbinsX();
2738  
2739  Double_t dV2 = 0.;
2740  Double_t dV4 = 0.;
2741  Double_t dV6 = 0.;
2742  Double_t dV8 = 0.; 
2743  
2744  if(!(*pW == "W"))
2745  {
2746   dV2 = fIntFlowResultsQC->GetBinContent(1); 
2747   dV4 = fIntFlowResultsQC->GetBinContent(2); 
2748   dV6 = fIntFlowResultsQC->GetBinContent(3); 
2749   dV8 = fIntFlowResultsQC->GetBinContent(4); 
2750  } 
2751  else if(*pW == "W")
2752  {
2753   dV2 = fIntFlowResultsQCW->GetBinContent(1);  
2754   dV4 = fIntFlowResultsQCW->GetBinContent(2); 
2755   dV6 = fIntFlowResultsQCW->GetBinContent(3); 
2756   dV8 = fIntFlowResultsQCW->GetBinContent(4); 
2757  }    
2758  
2759  // 3D (pt,eta): 
2760  Double_t twoPrimePtEta   = 0.; // <<2'>> (pt,eta) 
2761  Double_t fourPrimePtEta  = 0.; // <<4'>> (pt,eta)  
2762  //Double_t sixPrimePtEta   = 0.; // <<6'>> (pt,eta) 
2763  //Double_t eightPrimePtEta = 0.; // <<8'>> (pt,eta) 
2764  Double_t secondOrderDiffFlowCumulantPtEta = 0.; // d_n{2,Q} (pt,eta)
2765  Double_t fourthOrderDiffFlowCumulantPtEta = 0.; // d_n{4,Q} (pt,eta) 
2766  //Double_t sixthOrderDiffFlowCumulantPtEta = 0.; // d_n{6,Q} (pt,eta)
2767  //Double_t eightOrderDiffFlowCumulantPtEta = 0.; // d_n{8,Q} (pt,eta)2nd
2768  Double_t dv2PtEta = 0.; // v'_n{2} (pt,eta) 
2769  Double_t dv4PtEta = 0.; // v'_n{4} (pt,eta) 
2770  //Double_t dv6PtEta = 0.; // v'_n{6} (pt,eta) 
2771  //Double_t dv8PtEta = 0.; // v'_n{8} (pt,eta)  
2772
2773  // 2D (pt):   
2774  Double_t twoPrimePt   = 0.; // <<2'>> (pt)  
2775  Double_t fourPrimePt  = 0.; // <<4'>> (pt) 
2776  //Double_t sixPrimePt   = 0.; // <<6'>> (pt) 
2777  //Double_t eightPrimePt = 0.; // <<8'>> (pt)          
2778  Double_t secondOrderDiffFlowCumulantPt = 0.; // d_n{2,Q} (pt) 
2779  Double_t fourthOrderDiffFlowCumulantPt = 0.; // d_n{4,Q} (pt)  
2780  //Double_t sixthOrderDiffFlowCumulantPt = 0.; // d_n{6,Q} (pt)
2781  //Double_t eightOrderDiffFlowCumulantPt = 0.; // d_n{8,Q} (pt)
2782  Double_t dv2Pt = 0.; // v'_n{2} (pt)
2783  Double_t dv4Pt = 0.; // v'_n{4} (pt)
2784  //Double_t dv6Pt = 0.; // v'_n{6} (pt) 
2785  //Double_t dv8Pt = 0.; // v'_n{8} (pt)  
2786
2787  // 2D (eta):           
2788  Double_t twoPrimeEta   = 0.; // <<2'>> (eta)  
2789  Double_t fourPrimeEta  = 0.; // <<4>> (eta) 
2790  //Double_t sixPrimeEta   = 0.; // <<6>> (eta) 
2791  //Double_t eightPrimeEta = 0.; // <<8'>> (eta)  
2792  Double_t secondOrderDiffFlowCumulantEta = 0.; // d_n{2,Q} (eta)
2793  Double_t fourthOrderDiffFlowCumulantEta = 0.; // d_n{4,Q} (eta) 
2794  //Double_t sixthOrderDiffFlowCumulantEta = 0.; // d_n{6,Q} (eta) 
2795  //Double_t eightOrderDiffFlowCumulantEta = 0.; // d_n{8,Q} (eta) 
2796  Double_t dv2Eta = 0.; // v'_n{2} (eta)
2797  Double_t dv4Eta = 0.; // v'_n{4} (eta)
2798  //Double_t dv6Eta = 0.; // v'_n{6} (eta) 
2799  //Double_t dv8Eta = 0.; // v'_n{8} (eta)
2800  
2801
2802  // looping over (pt,eta) bins to calculate v'(pt,eta) 
2803  for(Int_t p=1;p<nBinsPt+1;p++)
2804  {
2805   for(Int_t e=1;e<nBinsEta+1;e++)
2806   {
2807   
2808    // 2nd order: 
2809    twoPrimePtEta = profile2ndPtEta->GetBinContent(profile2ndPtEta->GetBin(p,e));
2810    secondOrderDiffFlowCumulantPtEta = twoPrimePtEta;
2811    
2812
2813    //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
2814    // to be improved (applying correction for NUA):
2815    if(namePtEta->Contains("POI"))
2816    {
2817     if(f2pFinalCorrectionsForNUAPtEtaPOI) secondOrderDiffFlowCumulantPtEta = twoPrimePtEta 
2818                                      - f2pFinalCorrectionsForNUAPtEtaPOI->GetBinContent(f2pFinalCorrectionsForNUAPtEtaPOI->GetBin(p,e)) ;
2819    } else if (namePtEta->Contains("RP"))
2820      {  
2821       if(f2pFinalCorrectionsForNUAPtEtaRP) secondOrderDiffFlowCumulantPtEta = twoPrimePtEta 
2822                                        - f2pFinalCorrectionsForNUAPtEtaRP->GetBinContent(f2pFinalCorrectionsForNUAPtEtaRP->GetBin(p,e));
2823      }
2824    //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
2825    
2826    
2827    if(dV2)
2828    {
2829     dv2PtEta = secondOrderDiffFlowCumulantPtEta/dV2;
2830     if(*order2nd == "2") 
2831     {
2832      flowPtEta->SetBinContent(p,e,dv2PtEta);   
2833     } 
2834    }
2835    
2836    // 4th order: 
2837    if(*order4th == "4" || *order6th == "6" || *order8th == "8")
2838    {
2839     fourPrimePtEta = profile4thPtEta->GetBinContent(profile4thPtEta->GetBin(p,e));
2840     fourthOrderDiffFlowCumulantPtEta = fourPrimePtEta - 2.*twoPrimePtEta*pow(dV2,2.); // to be improved (correlations instead of pow(dV2,2.))
2841     if(dV4)
2842     {
2843      dv4PtEta = -fourthOrderDiffFlowCumulantPtEta/pow(dV4,3);
2844      if(*order4th == "4")
2845      {
2846       flowPtEta->SetBinContent(p,e,dv4PtEta);
2847      } 
2848     }
2849    }    
2850    
2851   } // end of for(Int_t e=1;e<nBinsEta+1;e++)
2852  } // end of for(Int_t p=1;p<nBinsPt+1;p++) 
2853    
2854    
2855  // looping over (pt) bins to calcualate v'(pt)
2856  for(Int_t p=1;p<nBinsPt+1;p++)
2857  {
2858  
2859   // 2nd order: 
2860   twoPrimePt = profile2ndPt->GetBinContent(p);
2861   secondOrderDiffFlowCumulantPt = twoPrimePt;
2862   
2863   
2864   //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
2865   // to be improved (applying correction for NUA):
2866   if(namePtEta->Contains("POI"))
2867   {
2868    if(f2pFinalCorrectionsForNUAPtPOI) secondOrderDiffFlowCumulantPt = twoPrimePt
2869                                     - f2pFinalCorrectionsForNUAPtPOI->GetBinContent(p) ;
2870   } else if (namePtEta->Contains("RP"))
2871     {
2872      if(f2pFinalCorrectionsForNUAPtRP) secondOrderDiffFlowCumulantPt = twoPrimePt
2873                                       - f2pFinalCorrectionsForNUAPtRP->GetBinContent(p);
2874     }
2875   //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
2876   
2877   
2878   if(dV2)
2879   {
2880    dv2Pt = secondOrderDiffFlowCumulantPt/dV2;
2881    if(*order2nd == "2") 
2882    {
2883     flowPt->SetBinContent(p,dv2Pt);
2884    }
2885    
2886    // common control histos: (to be improved fill only once. now they are filled first without weights and then with weights):
2887    if(namePtEta->Contains("POI") && *order2nd == "2")
2888    {
2889     fCommonHistsResults2nd->FillDifferentialFlowPtPOI(p,dv2Pt,0.); //to be improved (errors && bb or bb+1 ?)
2890    } 
2891    else if(namePtEta->Contains("RP") && *order2nd == "2")
2892    {
2893     fCommonHistsResults2nd->FillDifferentialFlowPtRP(p,dv2Pt,0.); //to be improved (errors && bb or bb+1 ?)
2894    }
2895    
2896   }
2897   
2898   // 4th order: 
2899   if(*order4th == "4" || *order6th == "6" || *order8th == "8")
2900   {
2901    fourPrimePt = profile4thPt->GetBinContent(profile4thPt->GetBin(p));
2902    fourthOrderDiffFlowCumulantPt = fourPrimePt - 2.*twoPrimePt*pow(dV2,2.); // to be improved (correlations instead of pow(dV2,2.))
2903    if(dV4)
2904    {
2905     dv4Pt = -fourthOrderDiffFlowCumulantPt/pow(dV4,3);
2906     if(*order4th == "4") 
2907     {
2908      flowPt->SetBinContent(p,dv4Pt);
2909     }
2910     
2911     // common control histos: (to be improved):
2912     if(namePtEta->Contains("POI") && *order4th == "4")
2913     {
2914      fCommonHistsResults4th->FillDifferentialFlowPtPOI(p,dv4Pt,0.); //to be improved (errors && bb or bb+1 ?)
2915     } 
2916     else if(namePtEta->Contains("RP") && *order4th == "4" )
2917     {
2918      fCommonHistsResults4th->FillDifferentialFlowPtRP(p,dv4Pt,0.); //to be improved (errors && bb or bb+1 ?)
2919     }
2920         
2921    }
2922   }    
2923   
2924  } // end of for(Int_t p=1;p<nBinsPt+1;p++)  
2925  
2926  
2927  // looping over (eta) bins to calcualate v'(eta)
2928  for(Int_t e=1;e<nBinsEta+1;e++)
2929  {
2930  
2931   // 2nd order: 
2932   twoPrimeEta = profile2ndEta->GetBinContent(e);
2933   secondOrderDiffFlowCumulantEta = twoPrimeEta;
2934   
2935   
2936   //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
2937   // to be improved (applying correction for NUA):
2938   if(namePtEta->Contains("POI"))
2939   {
2940    if(f2pFinalCorrectionsForNUAEtaPOI) secondOrderDiffFlowCumulantEta = twoPrimeEta
2941                                     - f2pFinalCorrectionsForNUAEtaPOI->GetBinContent(e) ;
2942   } else if (namePtEta->Contains("RP"))
2943     {
2944      if(f2pFinalCorrectionsForNUAEtaRP) secondOrderDiffFlowCumulantEta = twoPrimeEta
2945                                       - f2pFinalCorrectionsForNUAEtaRP->GetBinContent(e);
2946     }
2947   //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
2948   
2949   
2950   if(dV2)
2951   {
2952    dv2Eta = secondOrderDiffFlowCumulantEta/dV2;
2953    if(*order2nd == "2") 
2954    {
2955     flowEta->SetBinContent(e,dv2Eta);
2956    }
2957    
2958    // common control histos: (to be improved):
2959    if(namePtEta->Contains("POI") && *order2nd == "2")
2960    {
2961     fCommonHistsResults2nd->FillDifferentialFlowEtaPOI(e,dv2Eta,0.); //to be improved (errors && bb or bb+1 ?)
2962    } 
2963    else if(namePtEta->Contains("RP") && *order2nd == "2")
2964    {
2965     fCommonHistsResults2nd->FillDifferentialFlowEtaRP(e,dv2Eta,0.); //to be improved (errors && bb or bb+1 ?)
2966    }
2967      
2968
2969   }
2970   
2971   // 4th order: 
2972   if(*order4th == "4" || *order6th == "6" || *order8th == "8")
2973   {
2974    fourPrimeEta = profile4thEta->GetBinContent(profile4thEta->GetBin(e));
2975    fourthOrderDiffFlowCumulantEta = fourPrimeEta - 2.*twoPrimeEta*pow(dV2,2.); // to be improved (correlations instead of pow(dV2,2.))
2976    if(dV4)
2977    {
2978     dv4Eta = -fourthOrderDiffFlowCumulantEta/pow(dV4,3);
2979     if(*order4th == "4")
2980     {
2981      flowEta->SetBinContent(e,dv4Eta);
2982     }
2983     
2984     // common control histos: (to be improved):
2985     if(namePtEta->Contains("POI") && *order4th == "4")
2986     {
2987      fCommonHistsResults4th->FillDifferentialFlowEtaPOI(e,dv4Eta,0.); //to be improved (errors && bb or bb+1 ?)
2988     } 
2989     else if(namePtEta->Contains("RP") && *order4th == "4")
2990     {
2991      fCommonHistsResults4th->FillDifferentialFlowEtaRP(e,dv4Eta,0.); //to be improved (errors && bb or bb+1 ?)
2992     }
2993    
2994    }
2995   }    
2996   
2997  } // end of for(Int_t e=1;e<nBinsEta+1;e++)    
2998     
2999  delete namePtEta;
3000  delete type;
3001  delete order2nd;
3002  delete order4th;
3003  delete order6th;
3004  delete order8th;
3005  delete pW;
3006  delete profile2ndPt;
3007  delete profile4thPt;
3008  delete profile6thPt;
3009  delete profile8thPt;
3010  delete profile2ndEta;
3011  delete profile4thEta;
3012  delete profile6thEta;
3013  delete profile8thEta;
3014
3015 } // end of AliFlowAnalysisWithQCumulants::CalculateFinalResultsForDifferentialFlow(Bool_t useParticleWeights, TString type)
3016 */
3017
3018 //================================================================================================================================
3019
3020
3021 void AliFlowAnalysisWithQCumulants::PrintFinalResultsForIntegratedFlow(TString type)
3022 {
3023  // printing on the screen the final results for integrated flow (NONAME, POI and RP) // to be improved (NONAME) 
3024  
3025  Int_t n = fHarmonic; 
3026  
3027  if(type == "NONAME" || type == "RP" || type == "POI")
3028  {
3029   if(!(fCommonHistsResults2nd && fCommonHistsResults4th && fCommonHistsResults6th && fCommonHistsResults8th))
3030   {
3031    cout<<"WARNING: fCommonHistsResults2nd && fCommonHistsResults4th && fCommonHistsResults6th && fCommonHistsResults8th"<<endl;
3032    cout<<"         is NULL in AFAWQC::PFRFIF() !!!!"<<endl;
3033   }
3034  } else
3035    {
3036     cout<<"WARNING: type in not from {NONAME, RP, POI} in AFAWQC::PFRFIF() !!!!"<<endl;
3037     exit(0);
3038    }
3039  
3040  Double_t dVn[4] = {0.}; // array to hold Vn{2}, Vn{4}, Vn{6} and Vn{8}   
3041  Double_t dVnErr[4] = {0.}; // array to hold errors of Vn{2}, Vn{4}, Vn{6} and Vn{8}   
3042  
3043  if(type == "NONAME")
3044  {
3045   dVn[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinContent(1); 
3046   dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlow())->GetBinError(1); 
3047   dVn[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinContent(1); 
3048   dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlow())->GetBinError(1); 
3049   dVn[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinContent(1); 
3050   dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlow())->GetBinError(1); 
3051   dVn[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinContent(1); 
3052   dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlow())->GetBinError(1); 
3053  } else if(type == "RP")
3054    {
3055     dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinContent(1); 
3056     dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowRP())->GetBinError(1); 
3057     dVn[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinContent(1); 
3058     dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowRP())->GetBinError(1); 
3059     dVn[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinContent(1); 
3060     dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowRP())->GetBinError(1); 
3061     dVn[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinContent(1); 
3062     dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlowRP())->GetBinError(1); 
3063    } else if(type == "POI")
3064      {
3065       dVn[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinContent(1); 
3066       dVnErr[0] = (fCommonHistsResults2nd->GetHistIntFlowPOI())->GetBinError(1); 
3067       dVn[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinContent(1); 
3068       dVnErr[1] = (fCommonHistsResults4th->GetHistIntFlowPOI())->GetBinError(1); 
3069       dVn[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinContent(1); 
3070       dVnErr[2] = (fCommonHistsResults6th->GetHistIntFlowPOI())->GetBinError(1); 
3071       dVn[3] = (fCommonHistsResults8th->GetHistIntFlowPOI())->GetBinContent(1); 
3072       dVnErr[3] = (fCommonHistsResults8th->GetHistIntFlowPOI())->GetBinError(1); 
3073      }
3074  
3075  TString title = " flow estimates from Q-cumulants"; 
3076  TString subtitle = "    ("; 
3077  
3078  if(!(fUsePhiWeights||fUsePtWeights||fUseEtaWeights))
3079  {
3080   subtitle.Append(type);
3081   subtitle.Append(", without weights)");
3082  } else  
3083    {
3084     subtitle.Append(type);
3085     subtitle.Append(", with weights)");
3086    }
3087   
3088  cout<<endl;
3089  cout<<"*************************************"<<endl;
3090  cout<<"*************************************"<<endl;
3091  cout<<title.Data()<<endl; 
3092  cout<<subtitle.Data()<<endl; 
3093  cout<<endl;
3094   
3095  for(Int_t i=0;i<4;i++)
3096  {
3097   if(dVn[i]>=0.)
3098   {
3099    cout<<"  v_"<<n<<"{"<<2*(i+1)<<"} = "<<dVn[i]<<" +/- "<<dVnErr[i]<<endl;
3100   }
3101   else
3102   {
3103    cout<<"  v_"<<n<<"{"<<2*(i+1)<<"} = Im"<<endl;
3104   }  
3105  }
3106
3107  cout<<endl;
3108  /*
3109  if(type == "NONAME")
3110  {
3111   cout<<"     nEvts = "<<nEvtsNoName<<", AvM = "<<dMultNoName<<endl; // to be improved
3112  }
3113  else if (type == "RP")
3114  {
3115   cout<<"     nEvts = "<<nEvtsRP<<", AvM = "<<dMultRP<<endl; // to be improved  
3116  } 
3117  else if (type == "POI")
3118  {
3119   cout<<"     nEvts = "<<nEvtsPOI<<", AvM = "<<dMultPOI<<endl; // to be improved  
3120  } 
3121  */
3122  cout<<"*************************************"<<endl;
3123  cout<<"*************************************"<<endl;
3124  cout<<endl; 
3125   
3126 }// end of AliFlowAnalysisWithQCumulants::PrintFinalResultsForIntegratedFlow(TString type="NONAME");
3127
3128
3129 //================================================================================================================================
3130
3131
3132 void AliFlowAnalysisWithQCumulants::CompareResultsFromNestedLoopsAndFromQVectorsForDiffFlow(Bool_t useParticleWeights)
3133 {
3134  // compare correlations needed for diff. flow calculated with nested loops and those calculated from Q-vectors
3135
3136  cout<<endl;
3137  cout<<"   *************************************"<<endl;
3138  cout<<"   **** cross-checking the formulas ****"<<endl;
3139  cout<<"   ****    for differential flow    ****"<<endl;
3140  cout<<"   ****                             ****"<<endl;
3141  cout<<"   ****        (pt,eta) bin:        ****"<<endl; 
3142  cout<<"   ****    1.1  < pt  <  1.2  GeV   ****"<<endl;  
3143  cout<<"   ****   -0.55 < eta < -0.525      ****"<<endl; 
3144  cout<<"   *************************************"<<endl;                             
3145  cout<<endl;
3146  
3147  if(!useParticleWeights)
3148  {                                        
3149   cout<<"<cos(n(psi1-phi2))> from Q-vectors    = "<<fCorrelationsPro[0][0][0][0]->GetBinContent(fCorrelationsPro[0][0][0][0]->GetBin(12,19))<<endl;
3150   cout<<"<cos(n(psi1-phi2))> from nested loops = "<<fDirectCorrelationsDiffFlow->GetBinContent(1)<<endl;
3151   cout<<endl;  
3152   cout<<"<cos(n(psi1+phi2-phi3-phi4))> from Q-vectors    = "<<fCorrelationsPro[0][0][0][1]->GetBinContent(fCorrelationsPro[0][0][0][1]->GetBin(12,19))<<endl;
3153   cout<<"<cos(n(psi1+phi2-phi3-phi4))> from nested loops = "<<fDirectCorrelationsDiffFlow->GetBinContent(41)<<endl;
3154   cout<<endl;   
3155   cout<<"****************************************************"<<endl;
3156   cout<<"****************************************************"<<endl;
3157   cout<<endl;
3158   /*
3159   cout<<"<cos(n(psi1))> from Q-vectors    = "<<fCorrectionsCosP1nPsiPtEtaPOI->GetBinContent(fCorrectionsCosP1nPsiPtEtaPOI->GetBin(12,19))<<endl;
3160   cout<<"<cos(n(psi1))> from nested loops = "<<fDirectCorrectionsDiffFlowCos->GetBinContent(1)<<endl;
3161   cout<<endl;  
3162   cout<<"<sin(n(psi1))> from Q-vectors    = "<<fCorrectionsSinP1nPsiPtEtaPOI->GetBinContent(fCorrectionsSinP1nPsiPtEtaPOI->GetBin(12,19))<<endl;
3163   cout<<"<sin(n(psi1))> from nested loops = "<<fDirectCorrectionsDiffFlowSin->GetBinContent(1)<<endl;
3164   cout<<endl;
3165   */
3166  }
3167  
3168  if(useParticleWeights)
3169  {
3170   cout<<"<w2 cos(n(psi1-phi2))> from Q-vectors (RP)   = "<<fCorrelationsPro[0][1][0][0]->GetBinContent(fCorrelationsPro[0][1][0][0]->GetBin(12,19))<<endl;
3171   cout<<"<w2 cos(n(psi1-phi2))> from Q-vectors (POI)  = "<<fCorrelationsPro[1][1][0][0]->GetBinContent(fCorrelationsPro[1][1][0][0]->GetBin(12,19))<<endl;
3172   cout<<"<w2 cos(n(psi1-phi2))> from nested loops     = "<<fDirectCorrelationsDiffFlowW->GetBinContent(1)<<endl;
3173   cout<<endl;  
3174   cout<<"<w2 w3 w4 cos(n(psi1+phi2-phi3-phi4))> from Q-vectors (RP)  = "<<fCorrelationsPro[0][1][0][1]->GetBinContent(fCorrelationsPro[0][1][0][1]->GetBin(12,19))<<endl;
3175   cout<<"<w2 w3 w4 cos(n(psi1+phi2-phi3-phi4))> from Q-vectors (POI) = "<<fCorrelationsPro[1][1][0][1]->GetBinContent(fCorrelationsPro[1][1][0][1]->GetBin(12,19))<<endl;
3176   cout<<"<w2 w3 w4 cos(n(psi1+phi2-phi3-phi4))> from nested loops    = "<<fDirectCorrelationsDiffFlowW->GetBinContent(41)<<endl;
3177   cout<<endl;  
3178  }
3179  
3180 } // end of void AliFlowAnalysisWithQCumulants::CompareResultsFromNestedLoopsAndFromQVectorsForDiffFlow()
3181
3182
3183 //================================================================================================================================
3184
3185
3186 void AliFlowAnalysisWithQCumulants::WriteHistograms(TString outputFileName)
3187 {
3188  //store the final results in output .root file
3189  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
3190  //output->WriteObject(fHistList, "cobjQC","SingleKey");
3191  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
3192  delete output;
3193 }
3194
3195
3196 //================================================================================================================================
3197
3198
3199 void AliFlowAnalysisWithQCumulants::BookCommonHistograms()
3200 {
3201  // book common control histograms and common histograms for final results
3202  
3203  // common control histogram (ALL events)
3204  fCommonHists = new AliFlowCommonHist("AliFlowCommonHistQC");
3205  fHistList->Add(fCommonHists);  
3206  
3207  // common control histogram (for events with 2 and more particles)
3208  fCommonHists2nd = new AliFlowCommonHist("AliFlowCommonHist2ndOrderQC");
3209  fHistList->Add(fCommonHists2nd);  
3210  
3211  // common control histogram (for events with 4 and more particles)
3212  fCommonHists4th = new AliFlowCommonHist("AliFlowCommonHist4thOrderQC");
3213  fHistList->Add(fCommonHists4th);  
3214  
3215  // common control histogram (for events with 6 and more particles)
3216  fCommonHists6th = new AliFlowCommonHist("AliFlowCommonHist6thOrderQC");
3217  fHistList->Add(fCommonHists6th);  
3218  
3219  // common control histogram (for events with 8 and more particles)
3220  fCommonHists8th = new AliFlowCommonHist("AliFlowCommonHist8thOrderQC");
3221  fHistList->Add(fCommonHists8th);  
3222   
3223  // common histograms for final results (calculated for events with 2 and more particles)
3224  fCommonHistsResults2nd = new AliFlowCommonHistResults("AliFlowCommonHistResults2ndOrderQC");
3225  fHistList->Add(fCommonHistsResults2nd); 
3226  
3227  // common histograms for final results (calculated for events with 4 and more particles)
3228  fCommonHistsResults4th = new AliFlowCommonHistResults("AliFlowCommonHistResults4thOrderQC");
3229  fHistList->Add(fCommonHistsResults4th);
3230  
3231  // common histograms for final results (calculated for events with 6 and more particles)
3232  fCommonHistsResults6th = new AliFlowCommonHistResults("AliFlowCommonHistResults6thOrderQC");
3233  fHistList->Add(fCommonHistsResults6th); 
3234  
3235  // common histograms for final results (calculated for events with 8 and more particles)
3236  fCommonHistsResults8th = new AliFlowCommonHistResults("AliFlowCommonHistResults8thOrderQC");
3237  fHistList->Add(fCommonHistsResults8th); 
3238  
3239 } // end of void AliFlowAnalysisWithQCumulants::BookCommonHistograms()
3240
3241
3242 //================================================================================================================================
3243
3244
3245 void AliFlowAnalysisWithQCumulants::BookAndFillWeightsHistograms()
3246 {
3247  // book and fill histograms which hold phi, pt and eta weights
3248
3249  if(!fWeightsList)
3250  {
3251   cout<<"WARNING: fWeightsList is NULL in AFAWQC::BAFWH() !!!!"<<endl;
3252   exit(0);  
3253  }
3254     
3255  TString fUseParticleWeightsName = "fUseParticleWeightsQC";
3256  fUseParticleWeightsName += fAnalysisLabel->Data();
3257  fUseParticleWeights = new TProfile(fUseParticleWeightsName.Data(),"0 = particle weight not used, 1 = particle weight used ",3,0,3);
3258  fUseParticleWeights->SetLabelSize(0.06);
3259  (fUseParticleWeights->GetXaxis())->SetBinLabel(1,"w_{#phi}");
3260  (fUseParticleWeights->GetXaxis())->SetBinLabel(2,"w_{p_{T}}");
3261  (fUseParticleWeights->GetXaxis())->SetBinLabel(3,"w_{#eta}");
3262  fUseParticleWeights->Fill(0.5,(Int_t)fUsePhiWeights);
3263  fUseParticleWeights->Fill(1.5,(Int_t)fUsePtWeights);
3264  fUseParticleWeights->Fill(2.5,(Int_t)fUseEtaWeights);
3265  fWeightsList->Add(fUseParticleWeights); 
3266   
3267  if(fUsePhiWeights)
3268  {
3269   if(fWeightsList->FindObject("phi_weights"))
3270   {
3271    fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
3272    if(fPhiWeights->GetBinWidth(1) != fPhiBinWidth)
3273    {
3274     cout<<"WARNING: fPhiWeights->GetBinWidth(1) != fPhiBinWidth in AFAWQC::BAFWH() !!!!        "<<endl;
3275     cout<<"         This indicates inconsistent binning in phi histograms throughout the code."<<endl;
3276     exit(0);
3277    }
3278   } else 
3279     {
3280      cout<<"WARNING: fWeightsList->FindObject(\"phi_weights\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
3281      exit(0);
3282     }
3283  } // end of if(fUsePhiWeights)
3284  
3285  if(fUsePtWeights) 
3286  {
3287   if(fWeightsList->FindObject("pt_weights"))
3288   {
3289    fPtWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("pt_weights"));
3290    if(fPtWeights->GetBinWidth(1) != fPtBinWidth)
3291    {
3292     cout<<"WARNING: fPtWeights->GetBinWidth(1) != fPtBinWidth in AFAWQC::BAFWH() !!!!         "<<endl;
3293     cout<<"         This indicates insconsistent binning in pt histograms throughout the code."<<endl;
3294     exit(0);
3295    }
3296   } else 
3297     {
3298      cout<<"WARNING: fWeightsList->FindObject(\"pt_weights\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
3299      exit(0);
3300     }
3301  } // end of if(fUsePtWeights)    
3302
3303  if(fUseEtaWeights) 
3304  {
3305   if(fWeightsList->FindObject("eta_weights"))
3306   {
3307    fEtaWeights = dynamic_cast<TH1D*>(fWeightsList->FindObject("eta_weights"));
3308    if(fEtaWeights->GetBinWidth(1) != fEtaBinWidth)
3309    {
3310     cout<<"WARNING: fEtaWeights->GetBinWidth(1) != fEtaBinWidth in AFAWQC::BAFWH() !!!!        "<<endl;
3311     cout<<"         This indicates insconsistent binning in eta histograms throughout the code."<<endl;
3312     exit(0);
3313    }
3314   } else 
3315     {
3316      cout<<"WARNING: fUseEtaWeights && fWeightsList->FindObject(\"eta_weights\") is NULL in AFAWQC::BAFWH() !!!!"<<endl;
3317      exit(0);
3318     }
3319  } // end of if(fUseEtaWeights)
3320  
3321 } // end of AliFlowAnalysisWithQCumulants::BookAndFillWeightsHistograms()
3322
3323
3324 //================================================================================================================================
3325
3326
3327 void AliFlowAnalysisWithQCumulants::BookEverythingForIntegratedFlow()
3328 {
3329  // book all objects for integrated flow 
3330  
3331  // a) common;
3332  // b) profiles;
3333  // c) results;
3334  
3335  // ****************
3336  // **** COMMON ****
3337  // **************** 
3338
3339  // Re[Q_{m*n,k}], Im[Q_{m*n,k}] and S_{p,k}^M: 
3340  fReQ  = new TMatrixD(4,9);
3341  fImQ  = new TMatrixD(4,9);
3342  fSMpk = new TMatrixD(8,9);
3343  
3344  // profile to hold average multiplicities and number of events for events with nRP>=0, nRP>=1, ... , and nRP>=8
3345  fAvMultiplicity = new TProfile("fAvMultiplicity","Average Multiplicities of RPs",9,0,9);
3346  fAvMultiplicity->SetTickLength(-0.01,"Y");
3347  fAvMultiplicity->SetMarkerStyle(25);
3348  fAvMultiplicity->SetLabelSize(0.05);
3349  fAvMultiplicity->SetLabelOffset(0.02,"Y");
3350  fAvMultiplicity->SetYTitle("Average Multiplicity");
3351  (fAvMultiplicity->GetXaxis())->SetBinLabel(1,"all evts");
3352  (fAvMultiplicity->GetXaxis())->SetBinLabel(2,"n_{RP} #geq 1");
3353  (fAvMultiplicity->GetXaxis())->SetBinLabel(3,"n_{RP} #geq 2");
3354  (fAvMultiplicity->GetXaxis())->SetBinLabel(4,"n_{RP} #geq 3");
3355  (fAvMultiplicity->GetXaxis())->SetBinLabel(5,"n_{RP} #geq 4");
3356  (fAvMultiplicity->GetXaxis())->SetBinLabel(6,"n_{RP} #geq 5");
3357  (fAvMultiplicity->GetXaxis())->SetBinLabel(7,"n_{RP} #geq 6");
3358  (fAvMultiplicity->GetXaxis())->SetBinLabel(8,"n_{RP} #geq 7");
3359  (fAvMultiplicity->GetXaxis())->SetBinLabel(9,"n_{RP} #geq 8");
3360  fIntFlowProfiles->Add(fAvMultiplicity);
3361   
3362  TString pWeightsFlag[2] = {"pWeights not used","pWeights used"};
3363  TString eWeightsFlag[2] = {"exact eWeights","non-exact eWeights"};
3364  TString nuaFlag[2] = {"not corrected","corrected"};
3365  TString sinCosFlag[2] = {"sin","cos"};
3366  TString powerFlag[2] = {"linear","quadratic"};
3367   
3368  for(Int_t pW=0;pW<1+(Int_t)(fUsePhiWeights||fUsePtWeights||fUseEtaWeights);pW++)
3369  {
3370   // ***************
3371   // **** E-B-E ****
3372   // ***************
3373   
3374   // average multiparticle correlations for single event calculated from Q-vectors
3375   // (Remark: binning is organized in the same way as in fQCorrelations[pW][eW] bellow):   
3376   fQCorrelationsEBE[pW] = new TH1D(Form("fQCorrelationsEBE: %s",pWeightsFlag[pW].Data()),Form("Average multi-particle correlations for single event calculated from Q-vectors (%s)",pWeightsFlag[pW].Data()),32,0,32);
3377   
3378   for(Int_t sc=0;sc<2;sc++)
3379   {
3380    // correction terms for non-uniform acceptance for single event calculated from Q-vectors
3381    // (Remark: binning is organized in the same way as in fQCorrections[pW][sc]):
3382    fQCorrectionsEBE[pW][sc] = new TH1D(Form("fQCorrectionsEBE: %s, %s terms",pWeightsFlag[pW].Data(),sinCosFlag[sc].Data()),Form("Correction terms for non-uniform acceptance for single event (%s, %s terms)",pWeightsFlag[pW].Data(),sinCosFlag[sc].Data()),10,0,10);  
3383   }
3384   
3385   for(Int_t eW=0;eW<2;eW++)
3386   {
3387    // ******************
3388    // **** PROFILES ****
3389    // ******************
3390   
3391    // final average multiparticle correlations for all events calculated from Q-vectors:
3392    fQCorrelations[pW][eW] = new TProfile(Form("fQCorrelations: %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data()),Form("Final average multi-particle correlations calculated from Q-vectors (%s, %s)",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data()),32,0,32,"s");
3393    fQCorrelations[pW][eW]->SetTickLength(-0.01,"Y");
3394    fQCorrelations[pW][eW]->SetMarkerStyle(25);
3395    fQCorrelations[pW][eW]->SetLabelSize(0.03);
3396    fQCorrelations[pW][eW]->SetLabelOffset(0.01,"Y");
3397    // 2-p correlations:
3398    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(1,"<<2>>_{n|n}");
3399    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(2,"<<2>>_{2n|2n}");
3400    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(3,"<<2>>_{3n|3n}");
3401    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(4,"<<2>>_{4n|4n}");
3402    // 3-p correlations:
3403    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(6,"<<3>>_{2n|n,n}");
3404    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(7,"<<3>>_{3n|2n,n}");
3405    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(8,"<<3>>_{4n|2n,2n}");
3406    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(9,"<<3>>_{4n|3n,n}");
3407    // 4-p correlations:
3408    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(11,"<<4>>_{n,n|n,n}"); 
3409    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(12,"<<4>>_{2n,n|2n,n}");
3410    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(13,"<<4>>_{2n,2n|2n,2n}");
3411    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(14,"<<4>>_{3n|n,n,n}");
3412    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(15,"<<4>>_{3n,n|3n,n}");
3413    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(16,"<<4>>_{3n,n|2n,2n}"); 
3414    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(17,"<<4>>_{4n|2n,n,n}");
3415    // 5-p correlations:
3416    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(19,"<<5>>_{2n|n,n,n,n}"); 
3417    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(20,"<<5>>_{2n,2n|2n,n,n}");
3418    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(21,"<<5>>_{3n,n|2n,n,n}");
3419    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(22,"<<5>>_{4n|n,n,n,n}");
3420    // 6-p correlations:
3421    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(24,"<<6>>_{n,n,n|n,n,n}");
3422    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(25,"<<6>>_{2n,n,n|2n,n,n}");
3423    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(26,"<<6>>_{2n,2n|n,n,n,n}");
3424    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(27,"<<6>>_{3n,n|n,n,n,n}");
3425    // 7-p correlations:
3426    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(29,"<<7>>_{2n,n,n|n,n,n,n}");
3427    // 8-p correlations:
3428    (fQCorrelations[pW][eW]->GetXaxis())->SetBinLabel(31,"<<8>>_{n,n,n,n|n,n,n,n}");
3429    // add fQCorrelations[0] to the list fIntFlowList:
3430    fIntFlowProfiles->Add(fQCorrelations[pW][eW]);
3431   
3432    // averages <<2><4>>, <<2><6>>, <<4><6>>, etc,  needed to calculate covariances:
3433    fQProducts[pW][eW] = new TProfile(Form("fQProducts: %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data()),Form("Averages of products (%s, %s)",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data()),6,0,6);
3434    fQProducts[pW][eW]->SetTickLength(-0.01,"Y");
3435    fQProducts[pW][eW]->SetMarkerStyle(25);
3436    fQProducts[pW][eW]->SetLabelSize(0.05);
3437    fQProducts[pW][eW]->SetLabelOffset(0.01,"Y");
3438    (fQProducts[pW][eW]->GetXaxis())->SetBinLabel(1,"<<2><4>>");
3439    (fQProducts[pW][eW]->GetXaxis())->SetBinLabel(2,"<<2><6>>");
3440    (fQProducts[pW][eW]->GetXaxis())->SetBinLabel(3,"<<2><8>>");
3441    (fQProducts[pW][eW]->GetXaxis())->SetBinLabel(4,"<<4><6>>");
3442    (fQProducts[pW][eW]->GetXaxis())->SetBinLabel(5,"<<4><8>>");
3443    (fQProducts[pW][eW]->GetXaxis())->SetBinLabel(6,"<<6><8>>");
3444    fIntFlowProfiles->Add(fQProducts[pW][eW]);
3445  
3446    for(Int_t sc=0;sc<2;sc++)
3447    {
3448     // final average correction terms for non-uniform acceptance calculated from Q-vectors: 
3449     fQCorrections[pW][eW][sc] = new TProfile(Form("fQCorrections: %s, %s, %s terms",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data(),sinCosFlag[sc].Data()),Form("Correction terms for non-uniform acceptance (%s, %s, %s terms)",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data(),sinCosFlag[sc].Data()),10,0,10,"s");
3450     fQCorrections[pW][eW][sc]->SetTickLength(-0.01,"Y");
3451     fQCorrections[pW][eW][sc]->SetMarkerStyle(25);
3452     fQCorrections[pW][eW][sc]->SetLabelSize(0.03);
3453     fQCorrections[pW][eW][sc]->SetLabelOffset(0.01,"Y");
3454     // ......................................................................... 
3455     // 1-p terms:
3456     (fQCorrections[pW][eW][sc]->GetXaxis())->SetBinLabel(1,Form("%s(n(#phi_{1}))>",sinCosFlag[sc].Data()));
3457     // 2-p terms:
3458     // 3-p terms:
3459     // ...
3460     // ......................................................................... 
3461     // add fQCorrectionsCos to the list fIntFlowList:
3462     fIntFlowProfiles->Add(fQCorrections[pW][eW][sc]);
3463    } // end of for(Int_t sc=0;sc<2;sc++) 
3464   
3465    // *****************
3466    // **** RESULTS ****
3467    // *****************
3468   
3469    // final results for average multi-particle correlations with correct statistical errors:
3470    fCorrelations[pW][eW] = new TH1D(Form("fCorrelations: %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data()),Form("Average multi-particle correlations (%s, %s)",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data()),32,0,32);
3471    fCorrelations[pW][eW]->SetTickLength(-0.01,"Y");
3472    fCorrelations[pW][eW]->SetMarkerStyle(25);
3473    fCorrelations[pW][eW]->SetLabelSize(0.03);
3474    fCorrelations[pW][eW]->SetLabelOffset(0.01,"Y");
3475    // 2-p correlations:
3476    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(1,"<<2>>_{n|n}");
3477    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(2,"<<2>>_{2n|2n}");
3478    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(3,"<<2>>_{3n|3n}");
3479    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(4,"<<2>>_{4n|4n}");
3480    // 3-p correlations:
3481    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(6,"<<3>>_{2n|n,n}");
3482    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(7,"<<3>>_{3n|2n,n}");
3483    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(8,"<<3>>_{4n|2n,2n}");
3484    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(9,"<<3>>_{4n|3n,n}");
3485    // 4-p correlations:
3486    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(11,"<<4>>_{n,n|n,n}"); 
3487    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(12,"<<4>>_{2n,n|2n,n}");
3488    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(13,"<<4>>_{2n,2n|2n,2n}");
3489    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(14,"<<4>>_{3n|n,n,n}");
3490    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(15,"<<4>>_{3n,n|3n,n}");
3491    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(16,"<<4>>_{3n,n|2n,2n}"); 
3492    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(17,"<<4>>_{4n|2n,n,n}");
3493    // 5-p correlations:
3494    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(19,"<<5>>_{2n|n,n,n,n}"); 
3495    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(20,"<<5>>_{2n,2n|2n,n,n}");
3496    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(21,"<<5>>_{3n,n|2n,n,n}");
3497    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(22,"<<5>>_{4n|n,n,n,n}");
3498    // 6-p correlations:
3499    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(24,"<<6>>_{n,n,n|n,n,n}");
3500    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(25,"<<6>>_{2n,n,n|2n,n,n}");
3501    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(26,"<<6>>_{2n,2n|n,n,n,n}");
3502    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(27,"<<6>>_{3n,n|n,n,n,n}");
3503    // 7-p correlations:
3504    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(29,"<<7>>_{2n,n,n|n,n,n,n}");
3505    // 8-p correlations:
3506    (fCorrelations[pW][eW]->GetXaxis())->SetBinLabel(31,"<<8>>_{n,n,n,n|n,n,n,n}");
3507    // add fCorrelations to the list fIntFlowList:
3508    fIntFlowResults->Add(fCorrelations[pW][eW]);
3509  
3510    // final corrections for non-uniform acceptance for QC{2}, QC{4}, QC{6} and QC{8}:
3511    fCorrections[pW][eW] = new TH1D(Form("fCorrections: %s, %s",pWeightsFlag[pW].Data(),eWeightsFlag[eW].Data()),Form("Corrections for non-uniform acceptance for Q-cumulants (%s, %s)",eWeightsFlag[eW].Data()),4,0,4);
3512    fCorrections[pW][eW]->SetLabelSize(0.05);
3513    fCorrections[pW][eW]->SetMarkerStyle(25);
3514    (fCorrections[pW][eW]->GetXaxis())->SetBinLabel(1,"corr. for QC{2}");
3515    (fCorrections[pW][eW]->GetXaxis())->SetBinLabel(2,"corr. for QC{4}");
3516    (fCorrections[pW][eW]->GetXaxis())->SetBinLabel(3,"corr. for QC{6}");