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