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