PWG/FLOW/Base/AliFlowAnalysisWithLeeYangZeros.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 //Description: Maker to analyze Flow using the LeeYangZeros method
  18 //             Equation numbers are from Big Paper (BP): Nucl. Phys. A 727, 373 (2003)
  19 //             Practical Guide (PG):    J. Phys. G: Nucl. Part. Phys. 30, S1213 (2004)
  20 //             Adapted from StFlowLeeYangZerosMaker.cxx
  21 //             by Markus Oldenberg and Art Poskanzer, LBNL
  22 //             with advice from Jean-Yves Ollitrault and Nicolas Borghini
  23 //
  24 //Author: Naomi van der Kolk (kolk@nikhef.nl)
  25 /////////////////////////////////////////////////////////////////////////////////////////
  26
  27 #include "Riostream.h"
  28 #include "TObject.h"
  29 #include "TMath.h"
  30 #include "TFile.h"
  31 #include "TList.h"
  32 #include "TVector2.h"
  33 #include "TH1F.h"
  34 #include "TComplex.h"
  35 #include "TProfile.h"
  36 #include "TDirectoryFile.h"
  37
  38 #include "AliFlowCommonConstants.h"
  39 #include "AliFlowLYZConstants.h"
  40 #include "AliFlowAnalysisWithLeeYangZeros.h"
  41 #include "AliFlowLYZHist1.h"
  42 #include "AliFlowLYZHist2.h"
  43 #include "AliFlowCommonHist.h"
  44 #include "AliFlowCommonHistResults.h"
  45 #include "AliFlowEventSimple.h"
  46 #include "AliFlowTrackSimple.h"
  47 #include "AliFlowVector.h"
  48
  49 using std::endl;
  50 using std::cout;
  51 using std::cerr;
  52 ClassImp(AliFlowAnalysisWithLeeYangZeros)
  53
  54   //-----------------------------------------------------------------------
  55
  56   AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros():
  57     fQsum(NULL),
  58     fQ2sum(0),
  59     fEventNumber(0),
  60     fFirstRun(kTRUE),
  61     fUseSum(kTRUE),
  62     fDoubleLoop(kFALSE),
  63     fDebug(kFALSE),
  64     fHistList(NULL),
  65     fFirstRunList(NULL),
  66     fHistVtheta(NULL),
  67     fHistProVetaRP(NULL),
  68     fHistProVetaPOI(NULL),
  69     fHistProVPtRP(NULL),
  70     fHistProVPtPOI(NULL),
  71     fHistR0theta(NULL),
  72     fHistProReDenom(NULL),
  73     fHistProImDenom(NULL),
  74     fHistReDtheta(NULL),
  75     fHistImDtheta(NULL),
  76     fHistQsumforChi(NULL),
  77     fCommonHists(NULL),
  78     fCommonHistsRes(NULL)
  79
  80 {
  81   //default constructor
  82   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros default constructor****"<<endl;
  83
  84   fHistList = new TList();
  85   fFirstRunList = new TList();
  86
  87   for(Int_t i = 0;i<5;i++)
  88     {
  89       fHist1[i]=0;
  90       fHist2RP[i]=0;
  91       fHist2POI[i]=0;
  92     }
  93
  94   fQsum = new TVector2();
  95
  96 }
  97
  98 //-----------------------------------------------------------------------
  99
 100
 101  AliFlowAnalysisWithLeeYangZeros::~AliFlowAnalysisWithLeeYangZeros()
 102  {
 103    //default destructor
 104    if (fDebug) cout<<"****~AliFlowAnalysisWithLeeYangZeros****"<<endl;
 105    delete fQsum;
 106    delete fHistList;
 107    delete fFirstRunList;
 108
 109  }
 110
 111 //-----------------------------------------------------------------------
 112
 113 void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TString* outputFileName)
 114 {
 115  //store the final results in output .root file
 116
 117   TFile *output = new TFile(outputFileName->Data(),"RECREATE");
 118   if (GetFirstRun()) {
 119     //output->WriteObject(fHistList, "cobjLYZ1","SingleKey");
 120     if (fUseSum) { fHistList->SetName("cobjLYZ1SUM");}
 121     else {fHistList->SetName("cobjLYZ1PROD");}
 122     fHistList->SetOwner(kTRUE);
 123     fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 124   }
 125   else {
 126     //output->WriteObject(fHistList, "cobjLYZ2","SingleKey");
 127     if (fUseSum) { fHistList->SetName("cobjLYZ2SUM"); }
 128     else { fHistList->SetName("cobjLYZ2PROD"); }
 129     fHistList->SetOwner(kTRUE);
 130     fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 131   }
 132   delete output;
 133 }
 134
 135 //-----------------------------------------------------------------------
 136
 137 void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TString outputFileName)
 138 {
 139  //store the final results in output .root file
 140
 141   TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 142   if (GetFirstRun()) {
 143     //output->WriteObject(fHistList, "cobjLYZ1","SingleKey");
 144     if (fUseSum) { fHistList->SetName("cobjLYZ1SUM");}
 145     else {fHistList->SetName("cobjLYZ1PROD");}
 146     fHistList->SetOwner(kTRUE);
 147     fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 148   }
 149   else {
 150     //output->WriteObject(fHistList, "cobjLYZ2","SingleKey");
 151     if (fUseSum) { fHistList->SetName("cobjLYZ2SUM"); }
 152     else { fHistList->SetName("cobjLYZ2PROD"); }
 153     fHistList->SetOwner(kTRUE);
 154     fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
 155   }
 156   delete output;
 157 }
 158
 159 //-----------------------------------------------------------------------
 160
 161 void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TDirectoryFile *outputFileName)
 162 {
 163  //store the final results in output .root file
 164  if (GetFirstRun()) {
 165    if (fUseSum) { fHistList->SetName("cobjLYZ1SUM");}
 166    else {fHistList->SetName("cobjLYZ1PROD");}
 167    fHistList->SetOwner(kTRUE);
 168    outputFileName->Add(fHistList);
 169    outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
 170  }
 171  else {
 172   if (fUseSum) { fHistList->SetName("cobjLYZ2SUM"); }
 173   else { fHistList->SetName("cobjLYZ2PROD"); }
 174   fHistList->SetOwner(kTRUE);
 175   outputFileName->Add(fHistList);
 176   outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
 177  }
 178 }
 179
 180 //-----------------------------------------------------------------------
 181
 182 Bool_t AliFlowAnalysisWithLeeYangZeros::Init()
 183 {
 184   //init method
 185   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Init()****"<<endl;
 186
 187   //save old value and prevent histograms from being added to directory
 188   //to avoid name clashes in case multiple analaysis objects are used
 189   //in an analysis
 190   Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
 191   TH1::AddDirectory(kFALSE);
 192
 193   // Book histograms
 194   Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
 195   Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 196   Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 197
 198   Double_t  dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
 199   Double_t  dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
 200   Double_t  dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
 201   Double_t  dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
 202
 203   //for control histograms
 204   if (fFirstRun){
 205     if (fUseSum) {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ1SUM");}
 206     else {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ1PROD");}
 207     fHistList->Add(fCommonHists);
 208     if (fUseSum) {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ1SUM");}
 209     else {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ1PROD");}
 210     fHistList->Add(fCommonHistsRes);
 211   }
 212   else {
 213     if (fUseSum) {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ2SUM");}
 214     else {fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ2PROD");}
 215     fHistList->Add(fCommonHists);
 216     if (fUseSum) {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ2SUM");}
 217     else {fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ2PROD");}
 218     fHistList->Add(fCommonHistsRes);
 219   }
 220
 221   TString nameChiHist;
 222   if (fUseSum) {nameChiHist = "Flow_QsumforChi_LYZSUM";}
 223   else {nameChiHist = "Flow_QsumforChi_LYZPROD";}
 224   fHistQsumforChi = new TH1F(nameChiHist.Data(),nameChiHist.Data(),3,-1.,2.);
 225   fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
 226   fHistQsumforChi->SetYTitle("value");
 227   fHistList->Add(fHistQsumforChi);
 228
 229   //for first loop over events
 230   if (fFirstRun){
 231     TString nameR0Hist;
 232     if (fUseSum) {nameR0Hist = "First_Flow_r0theta_LYZSUM";}
 233     else {nameR0Hist = "First_Flow_r0theta_LYZPROD";}
 234     fHistR0theta  = new TH1D(nameR0Hist.Data(),nameR0Hist.Data(),iNtheta,-0.5,iNtheta-0.5);
 235     fHistR0theta->SetXTitle("#theta");
 236     fHistR0theta->SetYTitle("r_{0}^{#theta}");
 237     fHistList->Add(fHistR0theta);
 238
 239     TString nameVHist;
 240     if (fUseSum) {nameVHist = "First_Flow_Vtheta_LYZSUM";}
 241     else {nameVHist = "First_Flow_Vtheta_LYZPROD";}
 242     fHistVtheta  = new TH1D(nameVHist.Data(),nameVHist.Data(),iNtheta,-0.5,iNtheta-0.5);
 243     fHistVtheta->SetXTitle("#theta");
 244     fHistVtheta->SetYTitle("V_{n}^{#theta}");
 245     fHistList->Add(fHistVtheta);
 246
 247     //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta
 248     for (Int_t theta=0;theta<iNtheta;theta++) {
 249       TString nameHist1;
 250       if (fUseSum) {nameHist1 = "AliFlowLYZHist1_";}
 251       else {nameHist1 = "AliFlowLYZHist1_";}
 252       nameHist1 += theta;
 253       fHist1[theta]=new AliFlowLYZHist1(theta, nameHist1, fUseSum);
 254       fHistList->Add(fHist1[theta]);
 255     }
 256
 257   }
 258   //for second loop over events
 259   else {
 260     TString nameReDenomHist;
 261     if (fUseSum) {nameReDenomHist = "Second_FlowPro_ReDenom_LYZSUM";}
 262     else {nameReDenomHist = "Second_FlowPro_ReDenom_LYZPROD";}
 263     fHistProReDenom = new TProfile(nameReDenomHist.Data(),nameReDenomHist.Data(), iNtheta, -0.5, iNtheta-0.5);
 264     fHistProReDenom->SetXTitle("#theta");
 265     fHistProReDenom->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
 266     fHistList->Add(fHistProReDenom);
 267
 268     TString nameImDenomHist;
 269     if (fUseSum) {nameImDenomHist = "Second_FlowPro_ImDenom_LYZSUM";}
 270     else {nameImDenomHist = "Second_FlowPro_ImDenom_LYZPROD";}
 271     fHistProImDenom = new TProfile(nameImDenomHist.Data(),nameImDenomHist.Data(), iNtheta, -0.5, iNtheta-0.5);
 272     fHistProImDenom->SetXTitle("#theta");
 273     fHistProImDenom->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
 274     fHistList->Add(fHistProImDenom);
 275
 276     TString nameVetaRPHist;
 277     if (fUseSum) {nameVetaRPHist = "Second_FlowPro_VetaRP_LYZSUM";}
 278     else {nameVetaRPHist = "Second_FlowPro_VetaRP_LYZPROD";}
 279     fHistProVetaRP = new TProfile(nameVetaRPHist.Data(),nameVetaRPHist.Data(),iNbinsEta,dEtaMin,dEtaMax);
 280     fHistProVetaRP->SetXTitle("rapidity");
 281     fHistProVetaRP->SetYTitle("v_{2}(#eta) for RP selection");
 282     fHistList->Add(fHistProVetaRP);
 283
 284     TString nameVetaPOIHist;
 285     if (fUseSum) {nameVetaPOIHist = "Second_FlowPro_VetaPOI_LYZSUM";}
 286     else {nameVetaPOIHist = "Second_FlowPro_VetaPOI_LYZPROD";}
 287     fHistProVetaPOI = new TProfile(nameVetaPOIHist.Data(),nameVetaPOIHist.Data(),iNbinsEta,dEtaMin,dEtaMax);
 288     fHistProVetaPOI->SetXTitle("rapidity");
 289     fHistProVetaPOI->SetYTitle("v_{2}(#eta) for POI selection");
 290     fHistList->Add(fHistProVetaPOI);
 291
 292     TString nameVPtRPHist;
 293     if (fUseSum) {nameVPtRPHist = "Second_FlowPro_VPtRP_LYZSUM";}
 294     else {nameVPtRPHist = "Second_FlowPro_VPtRP_LYZPROD";}
 295     fHistProVPtRP = new TProfile(nameVPtRPHist.Data(),nameVPtRPHist.Data(),iNbinsPt,dPtMin,dPtMax);
 296     fHistProVPtRP->SetXTitle("Pt");
 297     fHistProVPtRP->SetYTitle("v_{2}(p_{T}) for RP selection");
 298     fHistList->Add(fHistProVPtRP);
 299
 300     TString nameVPtPOIHist;
 301     if (fUseSum) {nameVPtPOIHist = "Second_FlowPro_VPtPOI_LYZSUM";}
 302     else {nameVPtPOIHist = "Second_FlowPro_VPtPOI_LYZPROD";}
 303     fHistProVPtPOI = new TProfile(nameVPtPOIHist.Data(),nameVPtPOIHist.Data(),iNbinsPt,dPtMin,dPtMax);
 304     fHistProVPtPOI->SetXTitle("p_{T}");
 305     fHistProVPtPOI->SetYTitle("v_{2}(p_{T}) for POI selection");
 306     fHistList->Add(fHistProVPtPOI);
 307
 308     if (fUseSum){
 309       fHistReDtheta = new TH1D("Second_Flow_ReDtheta_LYZSUM","Second_Flow_ReDtheta_LYZSUM",iNtheta, -0.5, (double)iNtheta-0.5);
 310       fHistReDtheta->SetXTitle("#theta");
 311       fHistReDtheta->SetYTitle("Re(D^{#theta})");
 312       fHistList->Add(fHistReDtheta);
 313
 314       fHistImDtheta = new TH1D("Second_Flow_ImDtheta_LYZSUM","Second_Flow_ImDtheta_LYZSUM",iNtheta, -0.5, (double)iNtheta-0.5);
 315       fHistImDtheta->SetXTitle("#theta");
 316       fHistImDtheta->SetYTitle("Im(D^{#theta})");
 317       fHistList->Add(fHistImDtheta);
 318     }
 319
 320     //class AliFlowLYZHist2 defines the histograms:
 321     for (Int_t theta=0;theta<iNtheta;theta++)  {
 322       TString nameRP = "AliFlowLYZHist2RP_";
 323       nameRP += theta;
 324       fHist2RP[theta]=new AliFlowLYZHist2(theta, "RP", nameRP, fUseSum);
 325       fHistList->Add(fHist2RP[theta]);
 326
 327       TString namePOI = "AliFlowLYZHist2POI_";
 328       namePOI += theta;
 329       fHist2POI[theta]=new AliFlowLYZHist2(theta, "POI", namePOI, fUseSum);
 330       fHistList->Add(fHist2POI[theta]);
 331     }
 332
 333     //read histogram fHistR0theta from the first run list
 334     if (fFirstRunList) {
 335       if (fUseSum) { fHistR0theta  = (TH1D*)fFirstRunList->FindObject("First_Flow_r0theta_LYZSUM");}
 336       else{ fHistR0theta  = (TH1D*)fFirstRunList->FindObject("First_Flow_r0theta_LYZPROD");}
 337       if (!fHistR0theta) {cout<<"fHistR0theta has a NULL pointer!"<<endl;}
 338       fHistList->Add(fHistR0theta);
 339     } else { cout<<"list is NULL pointer!"<<endl; }
 340
 341   }
 342
 343
 344   if (fDebug) cout<<"****Histograms initialised****"<<endl;
 345
 346   fEventNumber = 0; //set event counter to zero
 347
 348   //resore old stuff
 349   TH1::AddDirectory(oldHistAddStatus);
 350
 351   return kTRUE;
 352 }
 353
 354  //-----------------------------------------------------------------------
 355
 356 Bool_t AliFlowAnalysisWithLeeYangZeros::Make(AliFlowEventSimple* anEvent)
 357 {
 358   //make method
 359   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Make()****"<<endl;
 360
 361   //get tracks from event
 362   if (anEvent) {
 363     if (fFirstRun){
 364       fCommonHists->FillControlHistograms(anEvent);
 365       FillFromFlowEvent(anEvent);
 366     }
 367     else {
 368       fCommonHists->FillControlHistograms(anEvent);
 369       SecondFillFromFlowEvent(anEvent);
 370     }
 371   }
 372
 373   else {
 374     cout<<"##### FlowLeeYangZero: Stack pointer null"<<endl;
 375     return kFALSE;
 376   }
 377   //  cout<<"^^^^read event "<<fEventNumber<<endl;
 378   fEventNumber++;
 379
 380
 381   return kTRUE;
 382 }
 383
 384    //-----------------------------------------------------------------------
 385 void AliFlowAnalysisWithLeeYangZeros::GetOutputHistograms(TList *outputListHistos) {
 386  // get the pointers to all output histograms before calling Finish()
 387
 388   const Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
 389
 390   if (outputListHistos) {
 391
 392     //define histograms for first and second run
 393     AliFlowCommonHist *pCommonHist = NULL;
 394     AliFlowCommonHistResults *pCommonHistResults = NULL;
 395     TH1D*     pHistR0theta = NULL;
 396     TH1D*     pHistVtheta  = NULL;
 397     TProfile* pHistProReDenom = NULL;
 398     TProfile* pHistProImDenom = NULL;
 399     TProfile* pHistProVetaRP  = NULL;
 400     TProfile* pHistProVetaPOI = NULL;
 401     TProfile* pHistProVPtRP   = NULL;
 402     TProfile* pHistProVPtPOI  = NULL;
 403     TH1F* pHistQsumforChi = NULL;
 404     //AliFlowLYZHist1 *pLYZHist1[iNtheta] = {NULL};      //array of pointers to AliFlowLYZHist1
 405     //AliFlowLYZHist2 *pLYZHist2RP[iNtheta] = {NULL};    //array of pointers to AliFlowLYZHist2
 406     //AliFlowLYZHist2 *pLYZHist2POI[iNtheta] = {NULL};   //array of pointers to AliFlowLYZHist2
 407     AliFlowLYZHist1 **pLYZHist1 = new AliFlowLYZHist1*[iNtheta];      //array of pointers to AliFlowLYZHist1
 408     AliFlowLYZHist2 **pLYZHist2RP = new AliFlowLYZHist2*[iNtheta];    //array of pointers to AliFlowLYZHist2
 409     AliFlowLYZHist2 **pLYZHist2POI = new AliFlowLYZHist2*[iNtheta];   //array of pointers to AliFlowLYZHist2
 410     for (Int_t i=0; i<iNtheta; i++)
 411     {
 412       pLYZHist1[i] = NULL;
 413       pLYZHist2RP[i] = NULL;
 414       pLYZHist2POI[i] = NULL;
 415     }
 416
 417     if (GetFirstRun()) { //first run
 418       //Get the common histograms from the output list
 419       if (GetUseSum()){
 420         pCommonHist = dynamic_cast<AliFlowCommonHist*>
 421           (outputListHistos->FindObject("AliFlowCommonHistLYZ1SUM"));
 422         pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 423           (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ1SUM"));
 424         //Get the histograms from the output list
 425         for(Int_t theta = 0;theta<iNtheta;theta++){
 426           TString name = "AliFlowLYZHist1_";
 427           name += theta;
 428           pLYZHist1[theta] = dynamic_cast<AliFlowLYZHist1*>
 429             (outputListHistos->FindObject(name));
 430         }
 431         pHistVtheta = dynamic_cast<TH1D*>
 432           (outputListHistos->FindObject("First_Flow_Vtheta_LYZSUM"));
 433
 434         pHistR0theta = dynamic_cast<TH1D*>
 435           (outputListHistos->FindObject("First_Flow_r0theta_LYZSUM"));
 436
 437         pHistQsumforChi = dynamic_cast<TH1F*>
 438           (outputListHistos->FindObject("Flow_QsumforChi_LYZSUM"));
 439
 440         //Set the histogram pointers and call Finish()
 441         if (pCommonHist && pCommonHistResults && pLYZHist1[0] &&
 442             pHistVtheta && pHistR0theta && pHistQsumforChi ) {
 443           this->SetCommonHists(pCommonHist);
 444           this->SetCommonHistsRes(pCommonHistResults);
 445           this->SetHist1(pLYZHist1);
 446           this->SetHistVtheta(pHistVtheta);
 447           this->SetHistR0theta(pHistR0theta);
 448           this->SetHistQsumforChi(pHistQsumforChi);
 449         }
 450         else {
 451           cout<<"WARNING: Histograms needed to run Finish() firstrun (SUM) are not accessable!"<<endl;
 452         }
 453       }
 454       else {
 455         pCommonHist = dynamic_cast<AliFlowCommonHist*>
 456           (outputListHistos->FindObject("AliFlowCommonHistLYZ1PROD"));
 457         pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 458           (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ1PROD"));
 459         //Get the histograms from the output list
 460         for(Int_t theta = 0;theta<iNtheta;theta++){
 461           TString name = "AliFlowLYZHist1_";
 462           name += theta;
 463           pLYZHist1[theta] = dynamic_cast<AliFlowLYZHist1*>
 464             (outputListHistos->FindObject(name));
 465         }
 466         pHistVtheta = dynamic_cast<TH1D*>
 467           (outputListHistos->FindObject("First_Flow_Vtheta_LYZPROD"));
 468
 469         pHistR0theta = dynamic_cast<TH1D*>
 470           (outputListHistos->FindObject("First_Flow_r0theta_LYZPROD"));
 471
 472         pHistQsumforChi = dynamic_cast<TH1F*>
 473           (outputListHistos->FindObject("Flow_QsumforChi_LYZPROD"));
 474
 475         //Set the histogram pointers and call Finish()
 476         if (pCommonHist && pCommonHistResults && pLYZHist1[0] &&
 477             pHistVtheta && pHistR0theta && pHistQsumforChi ) {
 478           this->SetCommonHists(pCommonHist);
 479           this->SetCommonHistsRes(pCommonHistResults);
 480           this->SetHist1(pLYZHist1);
 481           this->SetHistVtheta(pHistVtheta);
 482           this->SetHistR0theta(pHistR0theta);
 483           this->SetHistQsumforChi(pHistQsumforChi);
 484         } else {
 485           cout<<"WARNING: Histograms needed to run Finish() firstrun (PROD) are not accessable!"<<endl;
 486         }
 487       }
 488     }
 489     else { //second run
 490       //Get the common histograms from the output list
 491       if (GetUseSum()){
 492         pCommonHist = dynamic_cast<AliFlowCommonHist*>
 493           (outputListHistos->FindObject("AliFlowCommonHistLYZ2SUM"));
 494         pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 495           (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ2SUM"));
 496
 497         pHistR0theta = dynamic_cast<TH1D*>
 498           (outputListHistos->FindObject("First_Flow_r0theta_LYZSUM"));
 499
 500         pHistQsumforChi = dynamic_cast<TH1F*>
 501           (outputListHistos->FindObject("Flow_QsumforChi_LYZSUM"));
 502
 503         //Get the histograms from the output list
 504         for(Int_t theta = 0;theta<iNtheta;theta++){
 505           TString nameRP = "AliFlowLYZHist2RP_";
 506           nameRP += theta;
 507           pLYZHist2RP[theta] = dynamic_cast<AliFlowLYZHist2*>
 508             (outputListHistos->FindObject(nameRP));
 509           TString namePOI = "AliFlowLYZHist2POI_";
 510           namePOI += theta;
 511           pLYZHist2POI[theta] = dynamic_cast<AliFlowLYZHist2*>
 512             (outputListHistos->FindObject(namePOI));
 513         }
 514         pHistProReDenom = dynamic_cast<TProfile*>
 515           (outputListHistos->FindObject("Second_FlowPro_ReDenom_LYZSUM"));
 516         pHistProImDenom = dynamic_cast<TProfile*>
 517           (outputListHistos->FindObject("Second_FlowPro_ImDenom_LYZSUM"));
 518
 519         TH1D* pHistReDtheta = dynamic_cast<TH1D*>
 520           (outputListHistos->FindObject("Second_Flow_ReDtheta_LYZSUM"));
 521         TH1D* pHistImDtheta = dynamic_cast<TH1D*>
 522           (outputListHistos->FindObject("Second_Flow_ImDtheta_LYZSUM"));
 523
 524         pHistProVetaRP = dynamic_cast<TProfile*>
 525           (outputListHistos->FindObject("Second_FlowPro_VetaRP_LYZSUM"));
 526         pHistProVetaPOI = dynamic_cast<TProfile*>
 527           (outputListHistos->FindObject("Second_FlowPro_VetaPOI_LYZSUM"));
 528         pHistProVPtRP = dynamic_cast<TProfile*>
 529           (outputListHistos->FindObject("Second_FlowPro_VPtRP_LYZSUM"));
 530         pHistProVPtPOI = dynamic_cast<TProfile*>
 531           (outputListHistos->FindObject("Second_FlowPro_VPtPOI_LYZSUM"));
 532
 533
 534         //Set the histogram pointers and call Finish()
 535         if (pCommonHist && pCommonHistResults &&
 536             pLYZHist2RP[0] && pLYZHist2POI[0] &&
 537             pHistR0theta &&
 538             pHistProReDenom && pHistProImDenom &&
 539             pHistReDtheta && pHistImDtheta &&
 540             pHistProVetaRP && pHistProVetaPOI &&
 541             pHistProVPtRP && pHistProVPtPOI &&
 542             pHistQsumforChi) {
 543           this->SetCommonHists(pCommonHist);
 544           this->SetCommonHistsRes(pCommonHistResults);
 545           this->SetHist2RP(pLYZHist2RP);
 546           this->SetHist2POI(pLYZHist2POI);
 547           this->SetHistR0theta(pHistR0theta);
 548           this->SetHistProReDenom(pHistProReDenom);
 549           this->SetHistProImDenom(pHistProImDenom);
 550           this->SetHistReDtheta(pHistReDtheta);
 551           this->SetHistImDtheta(pHistImDtheta);
 552           this->SetHistProVetaRP(pHistProVetaRP);
 553           this->SetHistProVetaPOI(pHistProVetaPOI);
 554           this->SetHistProVPtRP(pHistProVPtRP);
 555           this->SetHistProVPtPOI(pHistProVPtPOI);
 556           this->SetHistQsumforChi(pHistQsumforChi);
 557         }
 558         else {
 559           cout<<"WARNING: Histograms needed to run Finish() secondrun (SUM) are not accessable!"<<endl;
 560         }
 561       }
 562       else {
 563         pCommonHist = dynamic_cast<AliFlowCommonHist*>
 564           (outputListHistos->FindObject("AliFlowCommonHistLYZ2PROD"));
 565         pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
 566           (outputListHistos->FindObject("AliFlowCommonHistResultsLYZ2PROD"));
 567
 568
 569         pHistR0theta = dynamic_cast<TH1D*>
 570           (outputListHistos->FindObject("First_Flow_r0theta_LYZPROD"));
 571
 572         pHistQsumforChi = dynamic_cast<TH1F*>
 573           (outputListHistos->FindObject("Flow_QsumforChi_LYZPROD"));
 574
 575         //Get the histograms from the output list
 576         for(Int_t theta = 0;theta<iNtheta;theta++){
 577           TString nameRP = "AliFlowLYZHist2RP_";
 578           nameRP += theta;
 579           pLYZHist2RP[theta] = dynamic_cast<AliFlowLYZHist2*>
 580             (outputListHistos->FindObject(nameRP));
 581           TString namePOI = "AliFlowLYZHist2POI_";
 582           namePOI += theta;
 583           pLYZHist2POI[theta] = dynamic_cast<AliFlowLYZHist2*>
 584             (outputListHistos->FindObject(namePOI));
 585         }
 586         pHistProReDenom = dynamic_cast<TProfile*>
 587           (outputListHistos->FindObject("Second_FlowPro_ReDenom_LYZPROD"));
 588         pHistProImDenom = dynamic_cast<TProfile*>
 589           (outputListHistos->FindObject("Second_FlowPro_ImDenom_LYZPROD"));
 590
 591         pHistProVetaRP = dynamic_cast<TProfile*>
 592           (outputListHistos->FindObject("Second_FlowPro_VetaRP_LYZPROD"));
 593         pHistProVetaPOI = dynamic_cast<TProfile*>
 594           (outputListHistos->FindObject("Second_FlowPro_VetaPOI_LYZPROD"));
 595         pHistProVPtRP = dynamic_cast<TProfile*>
 596           (outputListHistos->FindObject("Second_FlowPro_VPtRP_LYZPROD"));
 597         pHistProVPtPOI = dynamic_cast<TProfile*>
 598           (outputListHistos->FindObject("Second_FlowPro_VPtPOI_LYZPROD"));
 599
 600         //Set the histogram pointers and call Finish()
 601         if (pCommonHist && pCommonHistResults && pLYZHist2RP[0] && pLYZHist2POI[0] &&
 602             pHistR0theta && pHistProReDenom && pHistProImDenom && pHistProVetaRP &&
 603             pHistProVetaPOI && pHistProVPtRP && pHistProVPtPOI && pHistQsumforChi) {
 604           this->SetCommonHists(pCommonHist);
 605           this->SetCommonHistsRes(pCommonHistResults);
 606           this->SetHist2RP(pLYZHist2RP);
 607           this->SetHist2POI(pLYZHist2POI);
 608           this->SetHistR0theta(pHistR0theta);
 609           this->SetHistProReDenom(pHistProReDenom);
 610           this->SetHistProImDenom(pHistProImDenom);
 611           this->SetHistProVetaRP(pHistProVetaRP);
 612           this->SetHistProVetaPOI(pHistProVetaPOI);
 613           this->SetHistProVPtRP(pHistProVPtRP);
 614           this->SetHistProVPtPOI(pHistProVPtPOI);
 615           this->SetHistQsumforChi(pHistQsumforChi);
 616         }
 617         else {
 618           cout<<"WARNING: Histograms needed to run Finish() secondrun (PROD) are not accessable!"<<endl;
 619         }
 620       }
 621     } //secondrun
 622     //outputListHistos->Print();
 623     delete [] pLYZHist1;
 624     delete [] pLYZHist2RP;
 625     delete [] pLYZHist2POI;
 626   } //listhistos
 627   else {
 628     cout << "histogram list pointer is empty in method AliFlowAnalysisWithLeeYangZeros::GetOutputHistograms() " << endl;}
 629 }
 630
 631   //-----------------------------------------------------------------------
 632  Bool_t AliFlowAnalysisWithLeeYangZeros::Finish()
 633 {
 634   //finish method
 635   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Finish()****"<<endl;
 636
 637   //define variables for both runs
 638   Double_t  dJ01 = 2.405;
 639   Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
 640
 641   //set the event number
 642   SetEventNumber((int)fCommonHists->GetHistQ()->GetEntries());
 643
 644   //Get multiplicity for RP selection
 645   Double_t dMultRP = fCommonHists->GetHistMultRP()->GetMean();
 646   if (fDebug) cout<<"The average multiplicity is "<<dMultRP<<endl;
 647
 648   //set the sum of Q vectors
 649   fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
 650   SetQ2sum(fHistQsumforChi->GetBinContent(3));
 651
 652   if (fFirstRun){
 653
 654     //define variables for the first run
 655     Double_t  dR0 = 0;
 656     Double_t  dVtheta = 0;
 657     Double_t  dv = 0;
 658     Double_t  dV = 0;
 659
 660     //reset histograms in case of merged output files
 661     fHistR0theta->Reset();
 662     fHistVtheta->Reset();
 663
 664     for (Int_t theta=0;theta<iNtheta;theta++)
 665       {
 666         //get the first minimum r0
 667         fHist1[theta]->FillGtheta();
 668         dR0 = fHist1[theta]->GetR0();
 669
 670         //calculate integrated flow
 671         if (!TMath::AreEqualAbs(dR0, 0., 1e-100)) { dVtheta = dJ01/dR0; }
 672         else { cout<<"r0 is not found! Leaving LYZ analysis."<<endl; return kFALSE; }
 673
 674         //for estimating systematic error resulting from d0
 675         Double_t dBinsize =0.;
 676         if (fUseSum){ dBinsize = (AliFlowLYZConstants::GetMaster()->GetMaxSUM())/(AliFlowLYZConstants::GetMaster()->GetNbins());}
 677         else { dBinsize = (AliFlowLYZConstants::GetMaster()->GetMaxPROD())/(AliFlowLYZConstants::GetMaster()->GetNbins());}
 678         Double_t dVplus = -1.;
 679         Double_t dVmin  = -1.;
 680         if (!TMath::AreEqualAbs(dR0+dBinsize, 0., 1e-100)) {dVplus = dJ01/(dR0+dBinsize);}
 681         if (!TMath::AreEqualAbs(dR0-dBinsize, 0., 1e-100)) {dVmin = dJ01/(dR0-dBinsize);}
 682         //convert V to v
 683         Double_t dvplus = -1.;
 684         Double_t dvmin= -1.;
 685         if (fUseSum){
 686           //for SUM: V=v because the Q-vector is scaled by 1/M
 687           dv = dVtheta;
 688           dvplus = dVplus;
 689           dvmin = dVmin; }
 690         else {
 691           //for PRODUCT: v=V/M
 692           if (!TMath::AreEqualAbs(dMultRP, 0., 1e-100)){
 693             dv = dVtheta/dMultRP;
 694             dvplus = dVplus/dMultRP;
 695             dvmin = dVmin/dMultRP; }}
 696
 697         if (fDebug) cout<<"dv = "<<dv<<" and dvplus = "<<dvplus<< " and dvmin = "<<dvmin<<endl;
 698
 699         //fill the histograms
 700         fHistR0theta->SetBinContent(theta+1,dR0);
 701         fHistR0theta->SetBinError(theta+1,0.0);
 702         fHistVtheta->SetBinContent(theta+1,dVtheta);
 703         fHistVtheta->SetBinError(theta+1,0.0);
 704         //get average value of fVtheta = fV
 705         dV += dVtheta;
 706       } //end of loop over theta
 707
 708     //get average value of fVtheta = fV
 709     dV /=iNtheta;
 710     if (!fUseSum) { if (dMultRP!=0.){dV /=dMultRP;}} //scale with multiplicity for PRODUCT
 711
 712     //sigma2 and chi
 713     Double_t  dSigma2 = 0;
 714     Double_t  dChi= 0;
 715     if (fEventNumber!=0) {
 716       *fQsum /= fEventNumber;
 717       //cout<<"fQsum is "<<fQsum->X()<<" "<<fQsum->Y()<<endl;
 718       fQ2sum /= fEventNumber;
 719       //cout<<"fQ2sum is "<<fQ2sum<<endl;
 720       dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
 721       //cout<<"dSigma2 is "<<dSigma2<<endl;
 722       if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
 723       else dChi = -1.;
 724       fCommonHistsRes->FillChi(dChi);
 725
 726       cout<<"*************************************"<<endl;
 727       cout<<"*************************************"<<endl;
 728       cout<<"      Integrated flow from           "<<endl;
 729       if (fUseSum) {
 730         cout<<"       Lee-Yang Zeroes SUM          "<<endl;}
 731       else {
 732         cout<<"     Lee-Yang Zeroes PRODUCT      "<<endl;}
 733       cout<<endl;
 734       cout<<"Chi = "<<dChi<<endl;
 735       //cout<<endl;
 736     }
 737
 738     // recalculate statistical errors on integrated flow
 739     //combining 5 theta angles to 1 relative error BP eq. 89
 740     Double_t dRelErr2comb = 0.;
 741     Int_t iEvts = fEventNumber;
 742     if (iEvts!=0) {
 743       for (Int_t theta=0;theta<iNtheta;theta++){
 744         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 745         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 746                                        TMath::Cos(dTheta));
 747         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 748                                       TMath::Cos(dTheta));
 749         dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
 750                           TMath::BesselJ1(dJ01)))*
 751           (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
 752            dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
 753       }
 754       dRelErr2comb /= iNtheta;
 755     }
 756     Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
 757
 758     //copy content of profile into TH1D and add error
 759     Double_t dv2pro = dV;   //in the case that fv is equal to fV
 760     Double_t dv2Err = dv2pro*dRelErrcomb ;
 761     cout<<"dV = "<<dv2pro<<" +- "<<dv2Err<<endl;
 762     cout<<endl;
 763     cout<<"*************************************"<<endl;
 764     cout<<"*************************************"<<endl;
 765     fCommonHistsRes->FillIntegratedFlow(dv2pro, dv2Err);
 766
 767
 768     if (fDebug) cout<<"****histograms filled****"<<endl;
 769
 770     return kTRUE;
 771   }  //firstrun
 772
 773
 774   else {  //second run to calculate differential flow
 775
 776     //declare variables for the second run
 777     TComplex cDenom, cNumerRP, cNumerPOI, cDtheta;
 778     Int_t m = 1;
 779     TComplex i = TComplex::I();
 780     Double_t dBesselRatio[3] = {1., 1.202, 2.69};
 781     Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 782     Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 783
 784     Double_t dR0 = 0.;
 785     Double_t dVtheta = 0.;
 786     Double_t dV = 0.;
 787     Double_t dReDenom = 0.;
 788     Double_t dImDenom = 0.;
 789
 790     //reset histograms in case of merged output files
 791     if (fUseSum) {
 792       fHistReDtheta->Reset();
 793       fHistImDtheta->Reset();
 794     }
 795     fHistProVetaRP ->Reset();
 796     fHistProVetaPOI->Reset();
 797     fHistProVPtRP  ->Reset();
 798     fHistProVPtPOI ->Reset();
 799
 800     //scale fHistR0theta by the number of merged files to undo the merging
 801     if (!fHistR0theta) { cout<<"Hist pointer R0theta in file does not exist"<<endl; }
 802     else {
 803       Int_t iEntries = (int)fHistR0theta->GetEntries();
 804       if (iEntries > iNtheta){
 805         //for each individual file fHistR0theta has iNtheta entries
 806         Int_t iFiles = iEntries/iNtheta;
 807         cout<<iFiles<<" files were merged!"<<endl;
 808         fHistR0theta->Scale(1./iFiles);
 809       }
 810
 811       //loop over theta
 812       for (Int_t theta=0;theta<iNtheta;theta++)  {
 813         dR0 = fHistR0theta->GetBinContent(theta+1); //histogram starts at bin 1
 814         if (fDebug) cerr<<"dR0 = "<<dR0<<endl;
 815         if (!TMath::AreEqualAbs(dR0, 0., 1e-100)) dVtheta = dJ01/dR0;
 816         dV += dVtheta;
 817         // BP Eq. 9  -> Vn^theta = j01/r0^theta
 818         if (!fHistProReDenom || !fHistProImDenom) {
 819           cout << "Hist pointer fDenom in file does not exist" <<endl;
 820           cout<< "Leaving LYZ second pass analysis!" <<endl;
 821           return kFALSE;
 822         } else {
 823           dReDenom = fHistProReDenom->GetBinContent(theta+1);
 824           dImDenom = fHistProImDenom->GetBinContent(theta+1);
 825           cDenom(dReDenom,dImDenom);
 826
 827           //for new method and use by others (only with the sum generating function):
 828           if (fUseSum) {
 829             cDtheta = dR0*cDenom/dJ01;
 830             fHistReDtheta->SetBinContent(theta+1,cDtheta.Re());
 831             fHistReDtheta->SetBinError(theta+1,0.0);
 832             fHistImDtheta->SetBinContent(theta+1,cDtheta.Im());
 833             fHistImDtheta->SetBinError(theta+1,0.0);
 834           }
 835
 836           cDenom *= TComplex::Power(i, m-1);
 837
 838           //v as a function of eta for RP selection
 839           for (Int_t be=1;be<=iNbinsEta;be++)  {
 840             Double_t dReRatioRP = 0.0;
 841             Double_t dEtaRP = fHist2RP[theta]->GetBinCenter(be);
 842             cNumerRP = fHist2RP[theta]->GetNumerEta(be);
 843             if (cNumerRP.Rho()==0) {
 844               if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero in Finish()"<<endl;
 845             }
 846             else if (TMath::AreEqualAbs(cDenom.Rho(), 0, 1e-100)) {
 847               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 848             }
 849             else {
 850               dReRatioRP = (cNumerRP/cDenom).Re();
 851             }
 852             Double_t dVetaRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
 853             fHistProVetaRP->Fill(dEtaRP,dVetaRP);
 854           } //loop over bins be
 855
 856
 857           //v as a function of eta for POI selection
 858           for (Int_t be=1;be<=iNbinsEta;be++)  {
 859             Double_t dReRatioPOI = 0.0;
 860             Double_t dEtaPOI = fHist2POI[theta]->GetBinCenter(be);
 861             cNumerPOI = fHist2POI[theta]->GetNumerEta(be);
 862             if (cNumerPOI.Rho()==0) {
 863               if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero in Finish()"<<endl;
 864             }
 865             else if (cDenom.Rho()==0) {
 866               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 867             }
 868             else {
 869               dReRatioPOI = (cNumerPOI/cDenom).Re();
 870             }
 871             Double_t dVetaPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
 872             fHistProVetaPOI->Fill(dEtaPOI,dVetaPOI);
 873           } //loop over bins be
 874
 875
 876           //v as a function of Pt for RP selection
 877           for (Int_t bp=1;bp<=iNbinsPt;bp++)  {
 878             Double_t dReRatioRP = 0.0;
 879             Double_t dPtRP = fHist2RP[theta]->GetBinCenterPt(bp);
 880             cNumerRP = fHist2RP[theta]->GetNumerPt(bp);
 881             if (cNumerRP.Rho()==0) {
 882               if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero"<<endl;
 883             }
 884             else if (cDenom.Rho()==0) {
 885               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 886             }
 887             else {
 888               dReRatioRP = (cNumerRP/cDenom).Re();
 889             }
 890             Double_t dVPtRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
 891             fHistProVPtRP->Fill(dPtRP,dVPtRP);
 892           } //loop over bins bp
 893
 894
 895           //v as a function of Pt for POI selection
 896           for (Int_t bp=1;bp<=iNbinsPt;bp++)  {
 897             Double_t dReRatioPOI = 0.0;
 898             Double_t dPtPOI = fHist2POI[theta]->GetBinCenterPt(bp);
 899             cNumerPOI = fHist2POI[theta]->GetNumerPt(bp);
 900             if (cNumerPOI.Rho()==0) {
 901               if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero"<<endl;
 902             }
 903             else if (cDenom.Rho()==0) {
 904               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 905             }
 906             else {
 907               dReRatioPOI = (cNumerPOI/cDenom).Re();
 908             }
 909             Double_t dVPtPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
 910             fHistProVPtPOI->Fill(dPtPOI,dVPtPOI);
 911           } //loop over bins bp
 912
 913         }
 914       }
 915
 916     }//end of loop over theta
 917
 918     //calculate the average of fVtheta = fV
 919     dV /= iNtheta;
 920     if (!fUseSum) { if (dMultRP!=0.) { dV /=dMultRP; }} //scale by the multiplicity for PRODUCT
 921     if (TMath::AreEqualAbs(dV, 0., 1e-100)) { cout<<"dV = 0! Leaving LYZ analysis."<<endl; return kFALSE; }
 922
 923     //sigma2 and chi (for statistical error calculations)
 924     Double_t  dSigma2 = 0;
 925     Double_t  dChi= 0;
 926     if (fEventNumber!=0) {
 927       *fQsum /= fEventNumber;
 928       //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
 929       //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
 930       fQ2sum /= fEventNumber;
 931       //cout<<"fQ2sum = "<<fQ2sum<<endl;
 932       dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
 933       if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
 934       else dChi = -1.;
 935       fCommonHistsRes->FillChi(dChi);
 936
 937       // recalculate statistical errors on integrated flow
 938       //combining 5 theta angles to 1 relative error BP eq. 89
 939       Double_t dRelErr2comb = 0.;
 940       Int_t iEvts = fEventNumber;
 941       if (iEvts!=0) {
 942         for (Int_t theta=0;theta<iNtheta;theta++){
 943         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 944         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 945                                        TMath::Cos(dTheta));
 946         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 947                                       TMath::Cos(dTheta));
 948         dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
 949                           TMath::BesselJ1(dJ01)))*
 950           (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
 951            dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
 952         }
 953         dRelErr2comb /= iNtheta;
 954       }
 955       Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
 956       Double_t dVErr = dV*dRelErrcomb ;
 957       fCommonHistsRes->FillIntegratedFlow(dV,dVErr);
 958
 959       cout<<"*************************************"<<endl;
 960       cout<<"*************************************"<<endl;
 961       cout<<"      Integrated flow from           "<<endl;
 962       if (fUseSum) {
 963         cout<<"       Lee-Yang Zeroes SUM          "<<endl;}
 964       else {
 965         cout<<"     Lee-Yang Zeroes PRODUCT      "<<endl;}
 966       cout<<endl;
 967       cout<<"Chi = "<<dChi<<endl;
 968       cout<<"dV = "<<dV<<" +- "<<dVErr<<endl;
 969       //cout<<endl;
 970     }
 971
 972     //copy content of profile into TH1D and add error and fill the AliFlowCommonHistResults
 973
 974     //v as a function of eta for RP selection
 975     for(Int_t b=0;b<iNbinsEta;b++) {
 976       Double_t dv2pro  = fHistProVetaRP->GetBinContent(b);
 977       Double_t dNprime = fCommonHists->GetEntriesInEtaBinRP(b);
 978       Double_t dErrdifcomb = 0.;  //set error to zero
 979       Double_t dErr2difcomb = 0.; //set error to zero
 980       //calculate error
 981       if (!TMath::AreEqualAbs(dNprime, 0., 1e-100)) {
 982         for (Int_t theta=0;theta<iNtheta;theta++) {
 983           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 984           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 985                                            TMath::Cos(dTheta));
 986           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 987                                           TMath::Cos(dTheta));
 988           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 989                                                  TMath::BesselJ1(dJ01)))*
 990             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 991              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 992         } //loop over theta
 993       }
 994
 995       if (!TMath::AreEqualAbs(dErr2difcomb, 0., 1e-100)) {
 996         dErr2difcomb /= iNtheta;
 997         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 998       }
 999       else {dErrdifcomb = 0.;}
1000       //fill TH1D
1001       fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dErrdifcomb);
1002     } //loop over bins b
1003
1004
1005     //v as a function of eta for POI selection
1006     for(Int_t b=0;b<iNbinsEta;b++) {
1007       Double_t dv2pro  = fHistProVetaPOI->GetBinContent(b);
1008       Double_t dNprime = fCommonHists->GetEntriesInEtaBinPOI(b);
1009       Double_t dErrdifcomb = 0.;  //set error to zero
1010       Double_t dErr2difcomb = 0.; //set error to zero
1011       //calculate error
1012       if (dNprime!=0.) {
1013         for (Int_t theta=0;theta<iNtheta;theta++) {
1014           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1015           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1016                                            TMath::Cos(dTheta));
1017           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1018                                           TMath::Cos(dTheta));
1019           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1020                                                TMath::BesselJ1(dJ01)))*
1021             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1022              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
1023         } //loop over theta
1024       }
1025
1026       if (dErr2difcomb!=0.) {
1027         dErr2difcomb /= iNtheta;
1028         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
1029       }
1030       else {dErrdifcomb = 0.;}
1031       //fill TH1D
1032       fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dErrdifcomb);
1033     } //loop over bins b
1034
1035
1036
1037     //v as a function of Pt for RP selection
1038     TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
1039     Double_t dVRP = 0.;
1040     Double_t dSum = 0.;
1041     Double_t dErrV =0.;
1042     for(Int_t b=0;b<iNbinsPt;b++) {
1043       Double_t dv2pro  = fHistProVPtRP->GetBinContent(b);
1044       Double_t dNprime = fCommonHists->GetEntriesInPtBinRP(b);
1045       Double_t dErrdifcomb = 0.;  //set error to zero
1046       Double_t dErr2difcomb = 0.; //set error to zero
1047       //calculate error
1048       if (dNprime!=0.) {
1049         for (Int_t theta=0;theta<iNtheta;theta++) {
1050           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1051           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1052                                            TMath::Cos(dTheta));
1053           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1054                                           TMath::Cos(dTheta));
1055           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1056                                                TMath::BesselJ1(dJ01)))*
1057             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1058              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
1059         } //loop over theta
1060       }
1061
1062       if (dErr2difcomb!=0.) {
1063         dErr2difcomb /= iNtheta;
1064         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
1065         //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
1066       }
1067       else {dErrdifcomb = 0.;}
1068
1069       //fill TH1D
1070       fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dErrdifcomb);
1071       //calculate integrated flow for RP selection
1072       if (fHistPtRP){
1073         Double_t dYieldPt = fHistPtRP->GetBinContent(b);
1074         dVRP += dv2pro*dYieldPt;
1075         dSum +=dYieldPt;
1076         dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
1077       } else { cout<<"fHistPtRP is NULL"<<endl; }
1078     } //loop over bins b
1079
1080     if (!TMath::AreEqualAbs(dSum, 0., 1e-100)) {
1081       dVRP /= dSum; //the pt distribution should be normalised
1082       dErrV /= (dSum*dSum);
1083       dErrV = TMath::Sqrt(dErrV);
1084     }
1085     cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
1086     //cout<<endl;
1087     fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
1088
1089
1090     //v as a function of Pt for POI selection
1091     TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
1092     Double_t dVPOI = 0.;
1093     dSum = 0.;
1094     dErrV =0.;
1095
1096     for(Int_t b=0;b<iNbinsPt;b++) {
1097       Double_t dv2pro = fHistProVPtPOI->GetBinContent(b);
1098       Double_t dNprime = fCommonHists->GetEntriesInPtBinPOI(b);
1099       Double_t dErrdifcomb = 0.;  //set error to zero
1100       Double_t dErr2difcomb = 0.; //set error to zero
1101       //calculate error
1102       if (dNprime!=0.) {
1103         for (Int_t theta=0;theta<iNtheta;theta++) {
1104           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1105           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1106                                            TMath::Cos(dTheta));
1107           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1108                                           TMath::Cos(dTheta));
1109           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1110                                                  TMath::BesselJ1(dJ01)))*
1111             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1112              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
1113         } //loop over theta
1114       }
1115
1116       if (dErr2difcomb!=0.) {
1117         dErr2difcomb /= iNtheta;
1118         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
1119         //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
1120       }
1121       else {dErrdifcomb = 0.;}
1122
1123       //fill TH1D
1124       fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dErrdifcomb);
1125       //calculate integrated flow for POI selection
1126       if (fHistPtPOI){
1127         Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
1128         dVPOI += dv2pro*dYieldPt;
1129         dSum +=dYieldPt;
1130         dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
1131       } else { cout<<"fHistPtPOI is NULL"<<endl; }
1132     } //loop over bins b
1133
1134     if (dSum != 0.) {
1135       dVPOI /= dSum; //the pt distribution should be normalised
1136       dErrV /= (dSum*dSum);
1137       dErrV = TMath::Sqrt(dErrV);
1138     }
1139     cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
1140     cout<<endl;
1141     cout<<"*************************************"<<endl;
1142     cout<<"*************************************"<<endl;
1143     fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
1144
1145   } //secondrun
1146
1147   //cout<<"----LYZ analysis finished....----"<<endl<<endl;
1148
1149   return kTRUE;
1150 }
1151
1152
1153 //-----------------------------------------------------------------------
1154
1155  Bool_t AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent(AliFlowEventSimple* anEvent)
1156 {
1157   // Get event quantities from AliFlowEvent for all particles
1158
1159   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl;
1160
1161   if (!anEvent){
1162     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
1163     return kFALSE;
1164   }
1165
1166   //define variables
1167   TComplex cExpo, cGtheta, cGthetaNew, cZ;
1168   Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
1169   Int_t iNbins = AliFlowLYZConstants::GetMaster()->GetNbins();
1170
1171   //calculate flow
1172   Double_t dOrder = 2.;
1173
1174   //get the Q vector
1175   AliFlowVector vQ = anEvent->GetQ();
1176   //weight by the multiplicity
1177   Double_t dQX = 0;
1178   Double_t dQY = 0;
1179   if (!TMath::AreEqualAbs(vQ.GetMult(), 0., 1e-100)) {
1180     dQX = vQ.X()/vQ.GetMult();
1181     dQY = vQ.Y()/vQ.GetMult();
1182   }
1183   vQ.Set(dQX,dQY);
1184
1185   //for chi calculation:
1186   *fQsum += vQ;
1187   fHistQsumforChi->SetBinContent(1,fQsum->X());
1188   fHistQsumforChi->SetBinContent(2,fQsum->Y());
1189   fQ2sum += vQ.Mod2();
1190   fHistQsumforChi->SetBinContent(3,fQ2sum);
1191
1192   for (Int_t theta=0;theta<iNtheta;theta++) {
1193     Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
1194
1195     //calculate dQtheta = cos(dOrder*(fPhi-dTheta);the projection of the Q vector on the reference direction dTheta
1196     Double_t dQtheta = GetQtheta(vQ, dTheta);
1197
1198     for (Int_t bin=1;bin<=iNbins;bin++) {
1199       Double_t dR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram  //FIXED???
1200       if (fUseSum) {
1201         //calculate the sum generating function
1202         cExpo(0.,dR*dQtheta); //Re=0 ; Im=dR*dQtheta
1203         cGtheta = TComplex::Exp(cExpo);
1204       }
1205       else {
1206         //calculate the product generating function
1207         cGtheta = GetGrtheta(anEvent, dR, dTheta);
1208         if (cGtheta.Rho2() > 100.) break;
1209       }
1210       //fill real and imaginary part of cGtheta
1211       fHist1[theta]->Fill(dR,cGtheta);
1212     } //loop over bins
1213   } //loop over theta
1214
1215   return kTRUE;
1216
1217 }
1218
1219  //-----------------------------------------------------------------------
1220  Bool_t AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent(AliFlowEventSimple* anEvent)
1221 {
1222   //for differential flow
1223
1224   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl;
1225
1226   if (!anEvent){
1227     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
1228     return kFALSE;
1229   }
1230
1231   //define variables
1232   TComplex cExpo, cDenom, cNumerRP, cNumerPOI, cCosTermComplex;
1233   Double_t dR0 = 0.;
1234   Double_t dCosTermRP = 0.;
1235   Double_t dCosTermPOI = 0.;
1236   Double_t m = 1.;
1237   Double_t dOrder = 2.;
1238   Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
1239
1240   //get the Q vector
1241   AliFlowVector vQ = anEvent->GetQ();
1242   //weight by the multiplicity
1243   Double_t dQX = 0.;
1244   Double_t dQY = 0.;
1245   if (vQ.GetMult() != 0) {
1246     dQX = vQ.X()/vQ.GetMult();
1247     dQY = vQ.Y()/vQ.GetMult();
1248   }
1249   vQ.Set(dQX,dQY);
1250
1251   //for chi calculation:
1252   *fQsum += vQ;
1253   fHistQsumforChi->SetBinContent(1,fQsum->X());
1254   fHistQsumforChi->SetBinContent(2,fQsum->Y());
1255   fQ2sum += vQ.Mod2();
1256   fHistQsumforChi->SetBinContent(3,fQ2sum);
1257
1258   for (Int_t theta=0;theta<iNtheta;theta++)
1259     {
1260       Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
1261
1262       //calculate dQtheta = cos(dOrder*(dPhi-dTheta);the projection of the Q vector on the reference direction dTheta
1263       Double_t dQtheta = GetQtheta(vQ, dTheta);
1264
1265       //denominator for differential v
1266       if (fHistR0theta) {
1267         dR0 = fHistR0theta->GetBinContent(theta+1);
1268       }
1269       else { cout <<"pointer fHistR0theta does not exist" << endl;
1270       }
1271
1272       if (fUseSum) //sum generating function
1273         {
1274           cExpo(0.,dR0*dQtheta);
1275           cDenom = dQtheta*(TComplex::Exp(cExpo)); //BP eq 12
1276           //loop over tracks in event
1277           Int_t iNumberOfTracks = anEvent->NumberOfTracks();
1278           for (Int_t i=0;i<iNumberOfTracks;i++)  {
1279             AliFlowTrackSimple*  pTrack = anEvent->GetTrack(i);
1280             if (pTrack) {
1281               Double_t dEta = pTrack->Eta();
1282               Double_t dPt = pTrack->Pt();
1283               Double_t dPhi = pTrack->Phi();
1284               if (pTrack->InRPSelection()) { // RP selection
1285                 dCosTermRP = cos(m*dOrder*(dPhi-dTheta));
1286                 cNumerRP = dCosTermRP*(TComplex::Exp(cExpo));
1287                 if (cNumerRP.Rho()==0) { cerr<<"WARNING: modulus of cNumerRP is zero in SecondFillFromFlowEvent"<<endl;}
1288                 if (fDebug) { cerr<<"modulus of cNumerRP is "<<cNumerRP.Rho()<<endl;}
1289                 if (fHist2RP[theta]) {
1290                   fHist2RP[theta]->Fill(dEta,dPt,cNumerRP);
1291                 }
1292               }
1293               if (pTrack->InPOISelection()) { //POI selection
1294                 dCosTermPOI = cos(m*dOrder*(dPhi-dTheta));
1295                 cNumerPOI = dCosTermPOI*(TComplex::Exp(cExpo));
1296                 if (cNumerPOI.Rho()==0) { cerr<<"WARNING: modulus of cNumerPOI is zero in SecondFillFromFlowEvent"<<endl;}
1297                 if (fDebug) { cerr<<"modulus of cNumerPOI is "<<cNumerPOI.Rho()<<endl;}
1298                 if (fHist2POI[theta]) {
1299                   fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI);
1300                 }
1301               }
1302             } //if track
1303             else {cerr << "no particle!!!"<<endl;}
1304           } //loop over tracks
1305         } //sum
1306       else {    //product generating function
1307         cDenom = GetDiffFlow(anEvent, dR0, theta);
1308       }//product
1309
1310       if (fHistProReDenom && fHistProImDenom) {
1311         fHistProReDenom->Fill(theta,cDenom.Re());               //fill the real part of fDenom
1312         fHistProImDenom->Fill(theta,cDenom.Im());               //fill the imaginary part of fDenom
1313       }
1314       else { cout << "Pointers to cDenom  mising" << endl;}
1315
1316     }//end of loop over theta
1317
1318   return kTRUE;
1319
1320 }
1321  //-----------------------------------------------------------------------
1322  Double_t AliFlowAnalysisWithLeeYangZeros::GetQtheta(AliFlowVector aQ, Double_t aTheta)
1323 {
1324   //calculate Qtheta. Qtheta is the sum over all particles of cos(dOrder*(dPhi-dTheta)) BP eq. 3
1325   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl;
1326
1327   Double_t dQtheta = 0.;
1328   Double_t dOrder = 2.;
1329
1330   dQtheta = aQ.X()*cos(dOrder*aTheta)+aQ.Y()*sin(dOrder*aTheta);
1331
1332   return dQtheta;
1333
1334 }
1335
1336
1337 //-----------------------------------------------------------------------
1338 TComplex AliFlowAnalysisWithLeeYangZeros::GetGrtheta(AliFlowEventSimple* const anEvent, Double_t aR, Double_t aTheta)
1339 {
1340   // Product Generating Function for LeeYangZeros method
1341   // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1342
1343   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
1344
1345
1346   TComplex cG = TComplex::One();
1347   Double_t dOrder =  2.;
1348   Double_t dWgt = 1.;
1349   //Double_t dWgt = 1./anEvent->GetEventNSelTracksRP(); //weight with the multiplicity
1350
1351   Int_t iNumberOfTracks = anEvent->NumberOfTracks();
1352
1353   for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
1354     {
1355       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1356       if (pTrack){
1357         if (pTrack->InRPSelection()) {
1358           Double_t dPhi = pTrack->Phi();
1359           Double_t dGIm = aR * dWgt*cos(dOrder*(dPhi - aTheta));
1360           TComplex cGi(1., dGIm);
1361           cG *= cGi;     //product over all tracks
1362         }
1363       }
1364       else {cerr << "no particle pointer !!!"<<endl;}
1365     }//loop over tracks
1366
1367   return cG;
1368
1369 }
1370
1371
1372 //-----------------------------------------------------------------------
1373 TComplex AliFlowAnalysisWithLeeYangZeros::GetDiffFlow(AliFlowEventSimple* const anEvent, Double_t aR0, Int_t theta)
1374 {
1375   // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
1376   // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1377   // Also for v1 mixed harmonics: DF Eq. 5
1378   // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
1379
1380   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
1381
1382   TComplex cG = TComplex::One();
1383   TComplex cdGr0(0.,0.);
1384   Double_t dOrder =  2.;
1385   Double_t dWgt = 1.;
1386   //Double_t dWgt = 1./anEvent->GetEventNSelTracksRP(); //weight with the multiplicity
1387
1388   Int_t iNumberOfTracks = anEvent->NumberOfTracks();
1389
1390   Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
1391   Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
1392
1393   //for the denominator (use all RP selected particles)
1394   for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
1395     {
1396       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1397       if (pTrack){
1398         if (pTrack->InRPSelection()) {
1399           Double_t dPhi = pTrack->Phi();
1400           Double_t dCosTerm = dWgt*cos(dOrder*(dPhi - dTheta));
1401           //GetGr0theta
1402           Double_t dGIm = aR0 * dCosTerm;
1403           TComplex cGi(1., dGIm);
1404           TComplex cCosTermComplex(1., aR0*dCosTerm);
1405           cG *= cGi;     //product over all tracks
1406           //GetdGr0theta
1407           cdGr0 +=(dCosTerm / cCosTermComplex);  //sum over all tracks
1408         }
1409       } //if particle
1410       else {cerr << "no particle!!!"<<endl;}
1411     }//loop over tracks
1412
1413   //for the numerator
1414   for (Int_t i=0;i<iNumberOfTracks;i++)
1415     {
1416       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1417       if (pTrack){
1418         Double_t dEta = pTrack->Eta();
1419         Double_t dPt = pTrack->Pt();
1420         Double_t dPhi = pTrack->Phi();
1421         Double_t dCosTerm = cos(dOrder*(dPhi-dTheta));
1422         TComplex cCosTermComplex(1.,aR0*dCosTerm);
1423         //RP selection
1424         if (pTrack->InRPSelection()) {
1425           TComplex cNumerRP = cG*dCosTerm/cCosTermComplex;  //PG Eq. 9
1426           fHist2RP[theta]->Fill(dEta,dPt,cNumerRP);
1427         }
1428         //POI selection
1429         if (pTrack->InPOISelection()) {
1430           TComplex cNumerPOI = cG*dCosTerm/cCosTermComplex;  //PG Eq. 9
1431           fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI);
1432         }
1433       } //if particle
1434       else {cerr << "no particle pointer!!!"<<endl;}
1435     }//loop over tracks
1436
1437   TComplex cDenom = cG*cdGr0;
1438   return cDenom;
1439
1440 }
1441
1442 //-----------------------------------------------------------------------
1443