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->FillChiRP(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     fEventNumber =0; //set to zero for second round over events
 769   }  //firstrun
 770
 771
 772   else {  //second run to calculate differential flow
 773
 774     //declare variables for the second run
 775     TComplex cDenom, cNumerRP, cNumerPOI, cDtheta;
 776     Int_t m = 1;
 777     TComplex i = TComplex::I();
 778     Double_t dBesselRatio[3] = {1., 1.202, 2.69};
 779     Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
 780     Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
 781
 782     Double_t dR0 = 0.;
 783     Double_t dVtheta = 0.;
 784     Double_t dV = 0.;
 785     Double_t dReDenom = 0.;
 786     Double_t dImDenom = 0.;
 787
 788     //reset histograms in case of merged output files
 789     if (fUseSum) {
 790       fHistReDtheta->Reset();
 791       fHistImDtheta->Reset();
 792     }
 793     fHistProVetaRP ->Reset();
 794     fHistProVetaPOI->Reset();
 795     fHistProVPtRP  ->Reset();
 796     fHistProVPtPOI ->Reset();
 797
 798     //scale fHistR0theta by the number of merged files to undo the merging
 799     if (!fHistR0theta) { cout<<"Hist pointer R0theta in file does not exist"<<endl; }
 800     else {
 801       Int_t iEntries = (int)fHistR0theta->GetEntries();
 802       if (iEntries > iNtheta){
 803         //for each individual file fHistR0theta has iNtheta entries
 804         Int_t iFiles = iEntries/iNtheta;
 805         cout<<iFiles<<" files were merged!"<<endl;
 806         fHistR0theta->Scale(1./iFiles);
 807       }
 808
 809       //loop over theta
 810       for (Int_t theta=0;theta<iNtheta;theta++)  {
 811         dR0 = fHistR0theta->GetBinContent(theta+1); //histogram starts at bin 1
 812         if (fDebug) cerr<<"dR0 = "<<dR0<<endl;
 813         if (!TMath::AreEqualAbs(dR0, 0., 1e-100)) dVtheta = dJ01/dR0;
 814         dV += dVtheta;
 815         // BP Eq. 9  -> Vn^theta = j01/r0^theta
 816         if (!fHistProReDenom && !fHistProImDenom) {
 817           cout << "Hist pointer fDenom in file does not exist" <<endl;
 818         } else {
 819           dReDenom = fHistProReDenom->GetBinContent(theta+1);
 820           dImDenom = fHistProImDenom->GetBinContent(theta+1);
 821           cDenom(dReDenom,dImDenom);
 822
 823           //for new method and use by others (only with the sum generating function):
 824           if (fUseSum) {
 825             cDtheta = dR0*cDenom/dJ01;
 826             fHistReDtheta->SetBinContent(theta+1,cDtheta.Re());
 827             fHistReDtheta->SetBinError(theta+1,0.0);
 828             fHistImDtheta->SetBinContent(theta+1,cDtheta.Im());
 829             fHistImDtheta->SetBinError(theta+1,0.0);
 830           }
 831
 832           cDenom *= TComplex::Power(i, m-1);
 833
 834           //v as a function of eta for RP selection
 835           for (Int_t be=1;be<=iNbinsEta;be++)  {
 836             Double_t dReRatioRP = 0.0;
 837             Double_t dEtaRP = fHist2RP[theta]->GetBinCenter(be);
 838             cNumerRP = fHist2RP[theta]->GetNumerEta(be);
 839             if (cNumerRP.Rho()==0) {
 840               if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero in Finish()"<<endl;
 841             }
 842             else if (TMath::AreEqualAbs(cDenom.Rho(), 0, 1e-100)) {
 843               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 844             }
 845             else {
 846               dReRatioRP = (cNumerRP/cDenom).Re();
 847             }
 848             Double_t dVetaRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
 849             fHistProVetaRP->Fill(dEtaRP,dVetaRP);
 850           } //loop over bins be
 851
 852
 853           //v as a function of eta for POI selection
 854           for (Int_t be=1;be<=iNbinsEta;be++)  {
 855             Double_t dReRatioPOI = 0.0;
 856             Double_t dEtaPOI = fHist2POI[theta]->GetBinCenter(be);
 857             cNumerPOI = fHist2POI[theta]->GetNumerEta(be);
 858             if (cNumerPOI.Rho()==0) {
 859               if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero in Finish()"<<endl;
 860             }
 861             else if (cDenom.Rho()==0) {
 862               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 863             }
 864             else {
 865               dReRatioPOI = (cNumerPOI/cDenom).Re();
 866             }
 867             Double_t dVetaPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
 868             fHistProVetaPOI->Fill(dEtaPOI,dVetaPOI);
 869           } //loop over bins be
 870
 871
 872           //v as a function of Pt for RP selection
 873           for (Int_t bp=1;bp<=iNbinsPt;bp++)  {
 874             Double_t dReRatioRP = 0.0;
 875             Double_t dPtRP = fHist2RP[theta]->GetBinCenterPt(bp);
 876             cNumerRP = fHist2RP[theta]->GetNumerPt(bp);
 877             if (cNumerRP.Rho()==0) {
 878               if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero"<<endl;
 879             }
 880             else if (cDenom.Rho()==0) {
 881               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 882             }
 883             else {
 884               dReRatioRP = (cNumerRP/cDenom).Re();
 885             }
 886             Double_t dVPtRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
 887             fHistProVPtRP->Fill(dPtRP,dVPtRP);
 888           } //loop over bins bp
 889
 890
 891           //v as a function of Pt for POI selection
 892           for (Int_t bp=1;bp<=iNbinsPt;bp++)  {
 893             Double_t dReRatioPOI = 0.0;
 894             Double_t dPtPOI = fHist2POI[theta]->GetBinCenterPt(bp);
 895             cNumerPOI = fHist2POI[theta]->GetNumerPt(bp);
 896             if (cNumerPOI.Rho()==0) {
 897               if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero"<<endl;
 898             }
 899             else if (cDenom.Rho()==0) {
 900               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 901             }
 902             else {
 903               dReRatioPOI = (cNumerPOI/cDenom).Re();
 904             }
 905             Double_t dVPtPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
 906             fHistProVPtPOI->Fill(dPtPOI,dVPtPOI);
 907           } //loop over bins bp
 908
 909         }
 910       }
 911
 912     }//end of loop over theta
 913
 914     //calculate the average of fVtheta = fV
 915     dV /= iNtheta;
 916     if (!fUseSum) { if (dMultRP!=0.) { dV /=dMultRP; }} //scale by the multiplicity for PRODUCT
 917     if (TMath::AreEqualAbs(dV, 0., 1e-100)) { cout<<"dV = 0! Leaving LYZ analysis."<<endl; return kFALSE; }
 918
 919     //sigma2 and chi (for statistical error calculations)
 920     Double_t  dSigma2 = 0;
 921     Double_t  dChi= 0;
 922     if (fEventNumber!=0) {
 923       *fQsum /= fEventNumber;
 924       //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
 925       //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
 926       fQ2sum /= fEventNumber;
 927       //cout<<"fQ2sum = "<<fQ2sum<<endl;
 928       dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
 929       if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
 930       else dChi = -1.;
 931       fCommonHistsRes->FillChiRP(dChi);
 932
 933       // recalculate statistical errors on integrated flow
 934       //combining 5 theta angles to 1 relative error BP eq. 89
 935       Double_t dRelErr2comb = 0.;
 936       Int_t iEvts = fEventNumber;
 937       if (iEvts!=0) {
 938         for (Int_t theta=0;theta<iNtheta;theta++){
 939         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 940         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 941                                        TMath::Cos(dTheta));
 942         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 943                                       TMath::Cos(dTheta));
 944         dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
 945                           TMath::BesselJ1(dJ01)))*
 946           (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
 947            dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
 948         }
 949         dRelErr2comb /= iNtheta;
 950       }
 951       Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
 952       Double_t dVErr = dV*dRelErrcomb ;
 953       fCommonHistsRes->FillIntegratedFlow(dV,dVErr);
 954
 955       cout<<"*************************************"<<endl;
 956       cout<<"*************************************"<<endl;
 957       cout<<"      Integrated flow from           "<<endl;
 958       if (fUseSum) {
 959         cout<<"       Lee-Yang Zeroes SUM          "<<endl;}
 960       else {
 961         cout<<"     Lee-Yang Zeroes PRODUCT      "<<endl;}
 962       cout<<endl;
 963       cout<<"Chi = "<<dChi<<endl;
 964       cout<<"dV = "<<dV<<" +- "<<dVErr<<endl;
 965       //cout<<endl;
 966     }
 967
 968     //copy content of profile into TH1D and add error and fill the AliFlowCommonHistResults
 969
 970     //v as a function of eta for RP selection
 971     for(Int_t b=0;b<iNbinsEta;b++) {
 972       Double_t dv2pro  = fHistProVetaRP->GetBinContent(b);
 973       Double_t dNprime = fCommonHists->GetEntriesInEtaBinRP(b);
 974       Double_t dErrdifcomb = 0.;  //set error to zero
 975       Double_t dErr2difcomb = 0.; //set error to zero
 976       //calculate error
 977       if (!TMath::AreEqualAbs(dNprime, 0., 1e-100)) {
 978         for (Int_t theta=0;theta<iNtheta;theta++) {
 979           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 980           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 981                                            TMath::Cos(dTheta));
 982           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 983                                           TMath::Cos(dTheta));
 984           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 985                                                  TMath::BesselJ1(dJ01)))*
 986             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 987              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 988         } //loop over theta
 989       }
 990
 991       if (!TMath::AreEqualAbs(dErr2difcomb, 0., 1e-100)) {
 992         dErr2difcomb /= iNtheta;
 993         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 994       }
 995       else {dErrdifcomb = 0.;}
 996       //fill TH1D
 997       fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dErrdifcomb);
 998     } //loop over bins b
 999
1000
1001     //v as a function of eta for POI selection
1002     for(Int_t b=0;b<iNbinsEta;b++) {
1003       Double_t dv2pro  = fHistProVetaPOI->GetBinContent(b);
1004       Double_t dNprime = fCommonHists->GetEntriesInEtaBinPOI(b);
1005       Double_t dErrdifcomb = 0.;  //set error to zero
1006       Double_t dErr2difcomb = 0.; //set error to zero
1007       //calculate error
1008       if (dNprime!=0.) {
1009         for (Int_t theta=0;theta<iNtheta;theta++) {
1010           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1011           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1012                                            TMath::Cos(dTheta));
1013           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1014                                           TMath::Cos(dTheta));
1015           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1016                                                TMath::BesselJ1(dJ01)))*
1017             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1018              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
1019         } //loop over theta
1020       }
1021
1022       if (dErr2difcomb!=0.) {
1023         dErr2difcomb /= iNtheta;
1024         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
1025       }
1026       else {dErrdifcomb = 0.;}
1027       //fill TH1D
1028       fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dErrdifcomb);
1029     } //loop over bins b
1030
1031
1032
1033     //v as a function of Pt for RP selection
1034     TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
1035     Double_t dVRP = 0.;
1036     Double_t dSum = 0.;
1037     Double_t dErrV =0.;
1038     for(Int_t b=0;b<iNbinsPt;b++) {
1039       Double_t dv2pro  = fHistProVPtRP->GetBinContent(b);
1040       Double_t dNprime = fCommonHists->GetEntriesInPtBinRP(b);
1041       Double_t dErrdifcomb = 0.;  //set error to zero
1042       Double_t dErr2difcomb = 0.; //set error to zero
1043       //calculate error
1044       if (dNprime!=0.) {
1045         for (Int_t theta=0;theta<iNtheta;theta++) {
1046           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1047           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1048                                            TMath::Cos(dTheta));
1049           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1050                                           TMath::Cos(dTheta));
1051           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1052                                                TMath::BesselJ1(dJ01)))*
1053             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1054              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
1055         } //loop over theta
1056       }
1057
1058       if (dErr2difcomb!=0.) {
1059         dErr2difcomb /= iNtheta;
1060         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
1061         //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
1062       }
1063       else {dErrdifcomb = 0.;}
1064
1065       //fill TH1D
1066       fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dErrdifcomb);
1067       //calculate integrated flow for RP selection
1068       if (fHistPtRP){
1069         Double_t dYieldPt = fHistPtRP->GetBinContent(b);
1070         dVRP += dv2pro*dYieldPt;
1071         dSum +=dYieldPt;
1072         dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
1073       } else { cout<<"fHistPtRP is NULL"<<endl; }
1074     } //loop over bins b
1075
1076     if (!TMath::AreEqualAbs(dSum, 0., 1e-100)) {
1077       dVRP /= dSum; //the pt distribution should be normalised
1078       dErrV /= (dSum*dSum);
1079       dErrV = TMath::Sqrt(dErrV);
1080     }
1081     cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
1082     //cout<<endl;
1083     fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
1084
1085
1086     //v as a function of Pt for POI selection
1087     TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
1088     Double_t dVPOI = 0.;
1089     dSum = 0.;
1090     dErrV =0.;
1091
1092     for(Int_t b=0;b<iNbinsPt;b++) {
1093       Double_t dv2pro = fHistProVPtPOI->GetBinContent(b);
1094       Double_t dNprime = fCommonHists->GetEntriesInPtBinPOI(b);
1095       Double_t dErrdifcomb = 0.;  //set error to zero
1096       Double_t dErr2difcomb = 0.; //set error to zero
1097       //calculate error
1098       if (dNprime!=0.) {
1099         for (Int_t theta=0;theta<iNtheta;theta++) {
1100           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
1101           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
1102                                            TMath::Cos(dTheta));
1103           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
1104                                           TMath::Cos(dTheta));
1105           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
1106                                                  TMath::BesselJ1(dJ01)))*
1107             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
1108              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
1109         } //loop over theta
1110       }
1111
1112       if (dErr2difcomb!=0.) {
1113         dErr2difcomb /= iNtheta;
1114         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
1115         //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
1116       }
1117       else {dErrdifcomb = 0.;}
1118
1119       //fill TH1D
1120       fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dErrdifcomb);
1121       //calculate integrated flow for POI selection
1122       if (fHistPtPOI){
1123         Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
1124         dVPOI += dv2pro*dYieldPt;
1125         dSum +=dYieldPt;
1126         dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
1127       } else { cout<<"fHistPtPOI is NULL"<<endl; }
1128     } //loop over bins b
1129
1130     if (dSum != 0.) {
1131       dVPOI /= dSum; //the pt distribution should be normalised
1132       dErrV /= (dSum*dSum);
1133       dErrV = TMath::Sqrt(dErrV);
1134     }
1135     cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
1136     cout<<endl;
1137     cout<<"*************************************"<<endl;
1138     cout<<"*************************************"<<endl;
1139     fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
1140
1141   } //secondrun
1142
1143   //cout<<"----LYZ analysis finished....----"<<endl<<endl;
1144
1145   return kTRUE;
1146 }
1147
1148
1149 //-----------------------------------------------------------------------
1150
1151  Bool_t AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent(AliFlowEventSimple* anEvent)
1152 {
1153   // Get event quantities from AliFlowEvent for all particles
1154
1155   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl;
1156
1157   if (!anEvent){
1158     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
1159     return kFALSE;
1160   }
1161
1162   //define variables
1163   TComplex cExpo, cGtheta, cGthetaNew, cZ;
1164   Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
1165   Int_t iNbins = AliFlowLYZConstants::GetMaster()->GetNbins();
1166
1167   //calculate flow
1168   Double_t dOrder = 2.;
1169
1170   //get the Q vector
1171   AliFlowVector vQ = anEvent->GetQ();
1172   //weight by the multiplicity
1173   Double_t dQX = 0;
1174   Double_t dQY = 0;
1175   if (!TMath::AreEqualAbs(vQ.GetMult(), 0., 1e-100)) {
1176     dQX = vQ.X()/vQ.GetMult();
1177     dQY = vQ.Y()/vQ.GetMult();
1178   }
1179   vQ.Set(dQX,dQY);
1180
1181   //for chi calculation:
1182   *fQsum += vQ;
1183   fHistQsumforChi->SetBinContent(1,fQsum->X());
1184   fHistQsumforChi->SetBinContent(2,fQsum->Y());
1185   fQ2sum += vQ.Mod2();
1186   fHistQsumforChi->SetBinContent(3,fQ2sum);
1187
1188   for (Int_t theta=0;theta<iNtheta;theta++) {
1189     Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
1190
1191     //calculate dQtheta = cos(dOrder*(fPhi-dTheta);the projection of the Q vector on the reference direction dTheta
1192     Double_t dQtheta = GetQtheta(vQ, dTheta);
1193
1194     for (Int_t bin=1;bin<=iNbins;bin++) {
1195       Double_t dR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram  //FIXED???
1196       if (fUseSum) {
1197         //calculate the sum generating function
1198         cExpo(0.,dR*dQtheta); //Re=0 ; Im=dR*dQtheta
1199         cGtheta = TComplex::Exp(cExpo);
1200       }
1201       else {
1202         //calculate the product generating function
1203         cGtheta = GetGrtheta(anEvent, dR, dTheta);
1204         if (cGtheta.Rho2() > 100.) break;
1205       }
1206       //fill real and imaginary part of cGtheta
1207       fHist1[theta]->Fill(dR,cGtheta);
1208     } //loop over bins
1209   } //loop over theta
1210
1211   return kTRUE;
1212
1213 }
1214
1215  //-----------------------------------------------------------------------
1216  Bool_t AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent(AliFlowEventSimple* anEvent)
1217 {
1218   //for differential flow
1219
1220   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl;
1221
1222   if (!anEvent){
1223     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
1224     return kFALSE;
1225   }
1226
1227   //define variables
1228   TComplex cExpo, cDenom, cNumerRP, cNumerPOI, cCosTermComplex;
1229   Double_t dR0 = 0.;
1230   Double_t dCosTermRP = 0.;
1231   Double_t dCosTermPOI = 0.;
1232   Double_t m = 1.;
1233   Double_t dOrder = 2.;
1234   Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
1235
1236   //get the Q vector
1237   AliFlowVector vQ = anEvent->GetQ();
1238   //weight by the multiplicity
1239   Double_t dQX = 0.;
1240   Double_t dQY = 0.;
1241   if (vQ.GetMult() != 0) {
1242     dQX = vQ.X()/vQ.GetMult();
1243     dQY = vQ.Y()/vQ.GetMult();
1244   }
1245   vQ.Set(dQX,dQY);
1246
1247   //for chi calculation:
1248   *fQsum += vQ;
1249   fHistQsumforChi->SetBinContent(1,fQsum->X());
1250   fHistQsumforChi->SetBinContent(2,fQsum->Y());
1251   fQ2sum += vQ.Mod2();
1252   fHistQsumforChi->SetBinContent(3,fQ2sum);
1253
1254   for (Int_t theta=0;theta<iNtheta;theta++)
1255     {
1256       Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
1257
1258       //calculate dQtheta = cos(dOrder*(dPhi-dTheta);the projection of the Q vector on the reference direction dTheta
1259       Double_t dQtheta = GetQtheta(vQ, dTheta);
1260
1261       //denominator for differential v
1262       if (fHistR0theta) {
1263         dR0 = fHistR0theta->GetBinContent(theta+1);
1264       }
1265       else { cout <<"pointer fHistR0theta does not exist" << endl;
1266       }
1267
1268       if (fUseSum) //sum generating function
1269         {
1270           cExpo(0.,dR0*dQtheta);
1271           cDenom = dQtheta*(TComplex::Exp(cExpo)); //BP eq 12
1272           //loop over tracks in event
1273           Int_t iNumberOfTracks = anEvent->NumberOfTracks();
1274           for (Int_t i=0;i<iNumberOfTracks;i++)  {
1275             AliFlowTrackSimple*  pTrack = anEvent->GetTrack(i);
1276             if (pTrack) {
1277               Double_t dEta = pTrack->Eta();
1278               Double_t dPt = pTrack->Pt();
1279               Double_t dPhi = pTrack->Phi();
1280               if (pTrack->InRPSelection()) { // RP selection
1281                 dCosTermRP = cos(m*dOrder*(dPhi-dTheta));
1282                 cNumerRP = dCosTermRP*(TComplex::Exp(cExpo));
1283                 if (cNumerRP.Rho()==0) { cerr<<"WARNING: modulus of cNumerRP is zero in SecondFillFromFlowEvent"<<endl;}
1284                 if (fDebug) { cerr<<"modulus of cNumerRP is "<<cNumerRP.Rho()<<endl;}
1285                 if (fHist2RP[theta]) {
1286                   fHist2RP[theta]->Fill(dEta,dPt,cNumerRP);
1287                 }
1288               }
1289               if (pTrack->InPOISelection()) { //POI selection
1290                 dCosTermPOI = cos(m*dOrder*(dPhi-dTheta));
1291                 cNumerPOI = dCosTermPOI*(TComplex::Exp(cExpo));
1292                 if (cNumerPOI.Rho()==0) { cerr<<"WARNING: modulus of cNumerPOI is zero in SecondFillFromFlowEvent"<<endl;}
1293                 if (fDebug) { cerr<<"modulus of cNumerPOI is "<<cNumerPOI.Rho()<<endl;}
1294                 if (fHist2POI[theta]) {
1295                   fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI);
1296                 }
1297               }
1298             } //if track
1299             else {cerr << "no particle!!!"<<endl;}
1300           } //loop over tracks
1301         } //sum
1302       else {    //product generating function
1303         cDenom = GetDiffFlow(anEvent, dR0, theta);
1304       }//product
1305
1306       if (fHistProReDenom && fHistProImDenom) {
1307         fHistProReDenom->Fill(theta,cDenom.Re());               //fill the real part of fDenom
1308         fHistProImDenom->Fill(theta,cDenom.Im());               //fill the imaginary part of fDenom
1309       }
1310       else { cout << "Pointers to cDenom  mising" << endl;}
1311
1312     }//end of loop over theta
1313
1314   return kTRUE;
1315
1316 }
1317  //-----------------------------------------------------------------------
1318  Double_t AliFlowAnalysisWithLeeYangZeros::GetQtheta(AliFlowVector aQ, Double_t aTheta)
1319 {
1320   //calculate Qtheta. Qtheta is the sum over all particles of cos(dOrder*(dPhi-dTheta)) BP eq. 3
1321   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl;
1322
1323   Double_t dQtheta = 0.;
1324   Double_t dOrder = 2.;
1325
1326   dQtheta = aQ.X()*cos(dOrder*aTheta)+aQ.Y()*sin(dOrder*aTheta);
1327
1328   return dQtheta;
1329
1330 }
1331
1332
1333 //-----------------------------------------------------------------------
1334 TComplex AliFlowAnalysisWithLeeYangZeros::GetGrtheta(AliFlowEventSimple* const anEvent, Double_t aR, Double_t aTheta)
1335 {
1336   // Product Generating Function for LeeYangZeros method
1337   // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1338
1339   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
1340
1341
1342   TComplex cG = TComplex::One();
1343   Double_t dOrder =  2.;
1344   Double_t dWgt = 1.;
1345   //Double_t dWgt = 1./anEvent->GetEventNSelTracksRP(); //weight with the multiplicity
1346
1347   Int_t iNumberOfTracks = anEvent->NumberOfTracks();
1348
1349   for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
1350     {
1351       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1352       if (pTrack){
1353         if (pTrack->InRPSelection()) {
1354           Double_t dPhi = pTrack->Phi();
1355           Double_t dGIm = aR * dWgt*cos(dOrder*(dPhi - aTheta));
1356           TComplex cGi(1., dGIm);
1357           cG *= cGi;     //product over all tracks
1358         }
1359       }
1360       else {cerr << "no particle pointer !!!"<<endl;}
1361     }//loop over tracks
1362
1363   return cG;
1364
1365 }
1366
1367
1368 //-----------------------------------------------------------------------
1369 TComplex AliFlowAnalysisWithLeeYangZeros::GetDiffFlow(AliFlowEventSimple* const anEvent, Double_t aR0, Int_t theta)
1370 {
1371   // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
1372   // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1373   // Also for v1 mixed harmonics: DF Eq. 5
1374   // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
1375
1376   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
1377
1378   TComplex cG = TComplex::One();
1379   TComplex cdGr0(0.,0.);
1380   Double_t dOrder =  2.;
1381   Double_t dWgt = 1.;
1382   //Double_t dWgt = 1./anEvent->GetEventNSelTracksRP(); //weight with the multiplicity
1383
1384   Int_t iNumberOfTracks = anEvent->NumberOfTracks();
1385
1386   Int_t iNtheta = AliFlowLYZConstants::GetMaster()->GetNtheta();
1387   Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
1388
1389   //for the denominator (use all RP selected particles)
1390   for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
1391     {
1392       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1393       if (pTrack){
1394         if (pTrack->InRPSelection()) {
1395           Double_t dPhi = pTrack->Phi();
1396           Double_t dCosTerm = dWgt*cos(dOrder*(dPhi - dTheta));
1397           //GetGr0theta
1398           Double_t dGIm = aR0 * dCosTerm;
1399           TComplex cGi(1., dGIm);
1400           TComplex cCosTermComplex(1., aR0*dCosTerm);
1401           cG *= cGi;     //product over all tracks
1402           //GetdGr0theta
1403           cdGr0 +=(dCosTerm / cCosTermComplex);  //sum over all tracks
1404         }
1405       } //if particle
1406       else {cerr << "no particle!!!"<<endl;}
1407     }//loop over tracks
1408
1409   //for the numerator
1410   for (Int_t i=0;i<iNumberOfTracks;i++)
1411     {
1412       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1413       if (pTrack){
1414         Double_t dEta = pTrack->Eta();
1415         Double_t dPt = pTrack->Pt();
1416         Double_t dPhi = pTrack->Phi();
1417         Double_t dCosTerm = cos(dOrder*(dPhi-dTheta));
1418         TComplex cCosTermComplex(1.,aR0*dCosTerm);
1419         //RP selection
1420         if (pTrack->InRPSelection()) {
1421           TComplex cNumerRP = cG*dCosTerm/cCosTermComplex;  //PG Eq. 9
1422           fHist2RP[theta]->Fill(dEta,dPt,cNumerRP);
1423         }
1424         //POI selection
1425         if (pTrack->InPOISelection()) {
1426           TComplex cNumerPOI = cG*dCosTerm/cCosTermComplex;  //PG Eq. 9
1427           fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI);
1428         }
1429       } //if particle
1430       else {cerr << "no particle pointer!!!"<<endl;}
1431     }//loop over tracks
1432
1433   TComplex cDenom = cG*cdGr0;
1434   return cDenom;
1435
1436 }
1437
1438 //-----------------------------------------------------------------------
1439