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 #include "Riostream.h"                 //needed as include
  18 #include "TObject.h"                   //needed as include
  19 #include "AliFlowCommonConstants.h"    //needed as include
  20 #include "AliFlowLYZConstants.h"       //needed as include
  21 #include "AliFlowAnalysisWithLeeYangZeros.h"
  22 #include "AliFlowLYZHist1.h"           //needed as include
  23 #include "AliFlowLYZHist2.h"           //needed as include
  24 #include "AliFlowCommonHist.h"         //needed as include
  25 #include "AliFlowCommonHistResults.h"  //needed as include
  26 #include "AliFlowEventSimple.h"        //needed as include
  27 #include "AliFlowTrackSimple.h"        //needed as include
  28
  29 class AliFlowVector;
  30
  31 #include "TMath.h" //needed as include
  32 #include "TFile.h" //needed as include
  33 #include "TList.h"
  34
  35 class TComplex;
  36 class TProfile;
  37 class TH1F;
  38 class TH1D;
  39
  40
  41 //Description: Maker to analyze Flow using the LeeYangZeros method
  42 //             Equation numbers are from Big Paper (BP): Nucl. Phys. A 727, 373 (2003)
  43 //             Practical Guide (PG):    J. Phys. G: Nucl. Part. Phys. 30, S1213 (2004)
  44 //             Adapted from StFlowLeeYangZerosMaker.cxx
  45 //             by Markus Oldenberg and Art Poskanzer, LBNL
  46 //             with advice from Jean-Yves Ollitrault and Nicolas Borghini
  47 //
  48 //Author: Naomi van der Kolk (kolk@nikhef.nl)
  49
  50
  51
  52 ClassImp(AliFlowAnalysisWithLeeYangZeros)
  53
  54   //-----------------------------------------------------------------------
  55
  56   AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros():
  57     fQsum(NULL),
  58     fQ2sum(0),
  59     fEventNumber(0),
  60     fFirstRun(kTRUE),
  61     fUseSum(kTRUE),
  62     fDoubleLoop(kFALSE),
  63     fDebug(kFALSE),
  64     fHistList(NULL),
  65     fFirstRunList(NULL),
  66     fHistProVtheta(NULL),
  67     fHistProVetaRP(NULL),
  68     fHistProVetaPOI(NULL),
  69     fHistProVPtRP(NULL),
  70     fHistProVPtPOI(NULL),
  71     fHistProR0theta(NULL),
  72     fHistProReDenom(NULL),
  73     fHistProImDenom(NULL),
  74     fHistProReDtheta(NULL),
  75     fHistProImDtheta(NULL),
  76     fHistQsumforChi(NULL),
  77     fCommonHists(NULL),
  78     fCommonHistsRes(NULL)
  79
  80 {
  81   //default constructor
  82   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros default constructor****"<<endl;
  83
  84   fHistList = new TList();
  85   fFirstRunList = new TList();
  86
  87   for(Int_t i = 0;i<5;i++)
  88     {
  89       fHist1[i]=0;
  90       fHist2RP[i]=0;
  91       fHist2POI[i]=0;
  92     }
  93
  94   fQsum = new TVector2();
  95
  96 }
  97
  98 //-----------------------------------------------------------------------
  99
 100
 101  AliFlowAnalysisWithLeeYangZeros::~AliFlowAnalysisWithLeeYangZeros()
 102  {
 103    //default destructor
 104    if (fDebug) cout<<"****~AliFlowAnalysisWithLeeYangZeros****"<<endl;
 105    delete fQsum;
 106    delete fHistList;
 107    delete fFirstRunList;
 108
 109  }
 110
 111 //-----------------------------------------------------------------------
 112
 113 void AliFlowAnalysisWithLeeYangZeros::WriteHistograms(TString* outputFileName)
 114 {
 115  //store the final results in output .root file
 116
 117   TFile *output = new TFile(outputFileName->Data(),"RECREATE");
 118   if (GetFirstRun()) {
 119     output->WriteObject(fHistList, "cobjLYZ1","SingleKey");
 120   }
 121   else {
 122     output->WriteObject(fHistList, "cobjLYZ2","SingleKey");
 123   }
 124   delete output;
 125 }
 126
 127 //-----------------------------------------------------------------------
 128
 129 Bool_t AliFlowAnalysisWithLeeYangZeros::Init()
 130 {
 131   //init method
 132   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Init()****"<<endl;
 133
 134   // Book histograms
 135   Int_t iNtheta = AliFlowLYZConstants::kTheta;
 136   Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 137   Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
 138
 139   Double_t  dPtMin = AliFlowCommonConstants::GetPtMin();
 140   Double_t  dPtMax = AliFlowCommonConstants::GetPtMax();
 141   Double_t  dEtaMin = AliFlowCommonConstants::GetEtaMin();
 142   Double_t  dEtaMax = AliFlowCommonConstants::GetEtaMax();
 143
 144   //for control histograms
 145   if (fFirstRun){ fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ1");
 146   fHistList->Add(fCommonHists);
 147   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ1");
 148   fHistList->Add(fCommonHistsRes);
 149   }
 150   else { fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZ2");
 151   fHistList->Add(fCommonHists);
 152   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZ2");
 153   fHistList->Add(fCommonHistsRes);
 154   }
 155
 156   fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZ","Flow_QsumforChi_LYZ",3,-1.,2.);
 157   fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
 158   fHistQsumforChi->SetYTitle("value");
 159   fHistList->Add(fHistQsumforChi);
 160
 161   //for first loop over events
 162   if (fFirstRun){
 163     fHistProR0theta  = new TProfile("First_FlowPro_r0theta_LYZ","First_FlowPro_r0theta_LYZ",iNtheta,-0.5,iNtheta-0.5);
 164     fHistProR0theta->SetXTitle("#theta");
 165     fHistProR0theta->SetYTitle("r_{0}^{#theta}");
 166     fHistList->Add(fHistProR0theta);
 167
 168     fHistProVtheta  = new TProfile("First_FlowPro_Vtheta_LYZ","First_FlowPro_Vtheta_LYZ",iNtheta,-0.5,iNtheta-0.5);
 169     fHistProVtheta->SetXTitle("#theta");
 170     fHistProVtheta->SetYTitle("V_{n}^{#theta}");
 171     fHistList->Add(fHistProVtheta);
 172
 173     //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta, fHistProR0theta
 174     for (Int_t theta=0;theta<iNtheta;theta++) {
 175       TString name = "AliFlowLYZHist1_";
 176       name += theta;
 177       fHist1[theta]=new AliFlowLYZHist1(theta, name);
 178       fHistList->Add(fHist1[theta]);
 179     }
 180
 181   }
 182   //for second loop over events
 183   else {
 184     fHistProReDenom = new TProfile("Second_FlowPro_ReDenom_LYZ","Second_FlowPro_ReDenom_LYZ" , iNtheta, -0.5, iNtheta-0.5);
 185     fHistProReDenom->SetXTitle("#theta");
 186     fHistProReDenom->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
 187     fHistList->Add(fHistProReDenom);
 188
 189     fHistProImDenom = new TProfile("Second_FlowPro_ImDenom_LYZ","Second_FlowPro_ImDenom_LYZ" , iNtheta, -0.5, iNtheta-0.5);
 190     fHistProImDenom->SetXTitle("#theta");
 191     fHistProImDenom->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
 192     fHistList->Add(fHistProImDenom);
 193
 194     fHistProVetaRP = new TProfile("Second_FlowPro_VetaRP_LYZ","Second_FlowPro_VetaRP_LYZ",iNbinsEta,dEtaMin,dEtaMax);
 195     fHistProVetaRP->SetXTitle("rapidity");
 196     fHistProVetaRP->SetYTitle("v_{2}(#eta) for RP selection");
 197     fHistList->Add(fHistProVetaRP);
 198
 199     fHistProVetaPOI = new TProfile("Second_FlowPro_VetaPOI_LYZ","Second_FlowPro_VetaPOI_LYZ",iNbinsEta,dEtaMin,dEtaMax);
 200     fHistProVetaPOI->SetXTitle("rapidity");
 201     fHistProVetaPOI->SetYTitle("v_{2}(#eta) for POI selection");
 202     fHistList->Add(fHistProVetaPOI);
 203
 204     fHistProVPtRP = new TProfile("Second_FlowPro_VPtRP_LYZ","Second_FlowPro_VPtRP_LYZ",iNbinsPt,dPtMin,dPtMax);
 205     fHistProVPtRP->SetXTitle("Pt");
 206     fHistProVPtRP->SetYTitle("v_{2}(p_{T}) for RP selection");
 207     fHistList->Add(fHistProVPtRP);
 208
 209     fHistProVPtPOI = new TProfile("Second_FlowPro_VPtPOI_LYZ","Second_FlowPro_VPtPOI_LYZ",iNbinsPt,dPtMin,dPtMax);
 210     fHistProVPtPOI->SetXTitle("p_{T}");
 211     fHistProVPtPOI->SetYTitle("v_{2}(p_{T}) for POI selection");
 212     fHistList->Add(fHistProVPtPOI);
 213
 214     fHistProReDtheta = new TProfile("Second_FlowPro_ReDtheta_LYZ","Second_FlowPro_ReDtheta_LYZ",iNtheta, -0.5, iNtheta-0.5);
 215     fHistProReDtheta->SetXTitle("#theta");
 216     fHistProReDtheta->SetYTitle("Re(D^{#theta})");
 217     fHistList->Add(fHistProReDtheta);
 218
 219     fHistProImDtheta = new TProfile("Second_FlowPro_ImDtheta_LYZ","Second_FlowPro_ImDtheta_LYZ",iNtheta, -0.5, iNtheta-0.5);
 220     fHistProImDtheta->SetXTitle("#theta");
 221     fHistProImDtheta->SetYTitle("Im(D^{#theta})");
 222     fHistList->Add(fHistProImDtheta);
 223
 224     //class AliFlowLYZHist2 defines the histograms:
 225     for (Int_t theta=0;theta<iNtheta;theta++)  {
 226       TString nameRP = "AliFlowLYZHist2RP_";
 227       nameRP += theta;
 228       fHist2RP[theta]=new AliFlowLYZHist2(theta,nameRP);
 229       fHistList->Add(fHist2RP[theta]);
 230
 231       TString namePOI = "AliFlowLYZHist2POI_";
 232       namePOI += theta;
 233       fHist2POI[theta]=new AliFlowLYZHist2(theta,namePOI);
 234       fHistList->Add(fHist2POI[theta]);
 235     }
 236
 237     //read histogram fHistProR0theta from the first run list
 238     if (fFirstRunList) {
 239       fHistProR0theta  = (TProfile*)fFirstRunList->FindObject("First_FlowPro_r0theta_LYZ");
 240       if (!fHistProR0theta) {cout<<"fHistProR0theta has a NULL pointer!"<<endl;}
 241       fHistList->Add(fHistProR0theta);
 242     } else { cout<<"list is NULL pointer!"<<endl; }
 243
 244   }
 245
 246
 247   if (fDebug) cout<<"****Histograms initialised****"<<endl;
 248
 249   fEventNumber = 0; //set event counter to zero
 250
 251   return kTRUE;
 252 }
 253
 254  //-----------------------------------------------------------------------
 255
 256 Bool_t AliFlowAnalysisWithLeeYangZeros::Make(AliFlowEventSimple* anEvent)
 257 {
 258   //make method
 259   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Make()****"<<endl;
 260
 261   //get tracks from event
 262   if (anEvent) {
 263     if (fFirstRun){
 264       fCommonHists->FillControlHistograms(anEvent);
 265       FillFromFlowEvent(anEvent);
 266     }
 267     else {
 268       fCommonHists->FillControlHistograms(anEvent);
 269       SecondFillFromFlowEvent(anEvent);
 270     }
 271   }
 272
 273   else {
 274     cout<<"##### FlowLeeYangZero: Stack pointer null"<<endl;
 275     return kFALSE;
 276   }
 277   cout<<"^^^^read event "<<fEventNumber<<endl;
 278   fEventNumber++;
 279
 280
 281   return kTRUE;
 282 }
 283
 284
 285   //-----------------------------------------------------------------------
 286  Bool_t AliFlowAnalysisWithLeeYangZeros::Finish()
 287 {
 288   //finish method
 289   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Finish()****"<<endl;
 290
 291   //define variables for both runs
 292   Double_t  dJ01 = 2.405;
 293   Int_t iNtheta = AliFlowLYZConstants::kTheta;
 294   //set the event number
 295   SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries());
 296   //cout<<"number of events processed is "<<fEventNumber<<endl;
 297
 298   //set the sum of Q vectors
 299   fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
 300   SetQ2sum(fHistQsumforChi->GetBinContent(3));
 301
 302   if (fFirstRun){
 303     Double_t  dR0 = 0;
 304     Double_t  dVtheta = 0;
 305     Double_t  dv = 0;
 306     Double_t  dV = 0;
 307     for (Int_t theta=0;theta<iNtheta;theta++)
 308       {
 309         //get the first minimum r0
 310         fHist1[theta]->FillGtheta();
 311         dR0 = fHist1[theta]->GetR0();
 312
 313         //calculate integrated flow
 314         if (dR0!=0) dVtheta = dJ01/dR0;
 315         //for estimating systematic error resulting from d0
 316         Double_t dBinsize = (AliFlowLYZConstants::fgMax)/(AliFlowLYZConstants::kNbins);
 317         Double_t dVplus = dJ01/(dR0+dBinsize);
 318         Double_t dVmin = dJ01/(dR0-dBinsize);
 319         dv = dVtheta;
 320         Double_t dvplus = dVplus;
 321         Double_t dvmin = dVmin;
 322         if (fDebug) cout<<"dv = "<<dv<<" and dvplus = "<<dvplus<< " and dvmin = "<<dvmin<<endl;
 323
 324         //fill the histograms
 325         fHistProR0theta->Fill(theta,dR0);
 326         fHistProVtheta->Fill(theta,dVtheta);
 327         //get average value of fVtheta = fV
 328         dV += dVtheta;
 329       } //end of loop over theta
 330
 331     //get average value of fVtheta = fV
 332     dV /=iNtheta;
 333
 334     //sigma2 and chi
 335     Double_t  dSigma2 = 0;
 336     Double_t  dChi= 0;
 337     if (fEventNumber!=0) {
 338       *fQsum /= fEventNumber;
 339       fQ2sum /= fEventNumber;
 340       dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
 341       if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
 342       else dChi = -1.;
 343       fCommonHistsRes->FillChiRP(dChi);
 344
 345       cout<<"*************************************"<<endl;
 346       cout<<"*************************************"<<endl;
 347       cout<<"      Integrated flow from           "<<endl;
 348       cout<<"        Lee-Yang Zeroes              "<<endl;
 349       cout<<endl;
 350       cout<<"dV(RP) = "<<dV<<" and chi = "<<dChi<<endl;
 351       cout<<endl;
 352     }
 353
 354     // recalculate statistical errors on integrated flow
 355     //combining 5 theta angles to 1 relative error BP eq. 89
 356     Double_t dRelErr2comb = 0.;
 357     Int_t iEvts = fEventNumber;
 358     if (iEvts!=0) {
 359       for (Int_t theta=0;theta<iNtheta;theta++){
 360         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 361         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 362                                        TMath::Cos(dTheta));
 363         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 364                                       TMath::Cos(dTheta));
 365         dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
 366                           TMath::BesselJ1(dJ01)))*
 367           (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
 368            dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
 369       }
 370       dRelErr2comb /= iNtheta;
 371     }
 372     Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
 373
 374     //copy content of profile into TH1D and add error
 375     Double_t dv2pro = dV;   //in the case that fv is equal to fV
 376     Double_t dv2Err = dv2pro*dRelErrcomb ;
 377     cout<<"dV(RP) = "<<dv2pro<<" +- "<<dv2Err<<endl;
 378     cout<<endl;
 379     cout<<"*************************************"<<endl;
 380     cout<<"*************************************"<<endl;
 381     fCommonHistsRes->FillIntegratedFlowRP(dv2pro, dv2Err);
 382
 383
 384     if (fDebug) cout<<"****histograms filled****"<<endl;
 385
 386     return kTRUE;
 387     fEventNumber =0; //set to zero for second round over events
 388   }  //firstrun
 389
 390
 391   else {  //second run
 392
 393     //calculate differential flow
 394     //declare variables
 395     TComplex cDenom, cNumerRP, cNumerPOI, cDtheta;
 396     Int_t m = 1;
 397     TComplex i = TComplex::I();
 398     Double_t dBesselRatio[3] = {1., 1.202, 2.69};
 399     Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 400     Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
 401
 402     Double_t dEtaRP, dPtRP, dReRatioRP, dVetaRP, dVPtRP, dEtaPOI, dPtPOI, dReRatioPOI, dVetaPOI, dVPtPOI;
 403
 404     Double_t dR0 = 0.;
 405     Double_t dVtheta = 0.;
 406     Double_t dV = 0.;
 407     Double_t dReDenom = 0.;
 408     Double_t dImDenom = 0.;
 409     for (Int_t theta=0;theta<iNtheta;theta++)  { //loop over theta
 410       if (!fHistProR0theta) {
 411         cout << "Hist pointer R0theta in file does not exist" <<endl;
 412       } else {
 413         dR0 = fHistProR0theta->GetBinContent(theta+1); //histogram starts at bin 1
 414         if (fDebug) cerr<<"dR0 = "<<dR0<<endl;
 415         if (dR0!=0) dVtheta = dJ01/dR0;
 416         dV += dVtheta;
 417         // BP Eq. 9  -> Vn^theta = j01/r0^theta
 418         if (!fHistProReDenom && !fHistProImDenom) {
 419           cout << "Hist pointer fDenom in file does not exist" <<endl;
 420         } else {
 421           dReDenom = fHistProReDenom->GetBinContent(theta+1);
 422           dImDenom = fHistProImDenom->GetBinContent(theta+1);
 423           cDenom(dReDenom,dImDenom);
 424
 425           //for new method and use by others (only with the sum generating function):
 426           if (fUseSum) {
 427             cDtheta = dR0*cDenom/dJ01;
 428             fHistProReDtheta->Fill(theta,cDtheta.Re());
 429             fHistProImDtheta->Fill(theta,cDtheta.Im());
 430           }
 431
 432           cDenom *= TComplex::Power(i, m-1);
 433           //cerr<<"TComplex::Power(i, m-1) = "<<TComplex::Power(i, m-1).Rho()<<endl; //checked ok
 434
 435           //v as a function of eta for RP selection
 436           for (Int_t be=1;be<=iNbinsEta;be++)  {
 437             dEtaRP = fHist2RP[theta]->GetBinCenter(be);
 438             cNumerRP = fHist2RP[theta]->GetNumerEta(be);
 439             if (cNumerRP.Rho()==0) {
 440               if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero in Finish()"<<endl;
 441               dReRatioRP = 0;
 442             }
 443             else if (cDenom.Rho()==0) {
 444               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 445               dReRatioRP = 0;
 446             }
 447             else {
 448               dReRatioRP = (cNumerRP/cDenom).Re();
 449             }
 450             dVetaRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
 451             fHistProVetaRP->Fill(dEtaRP,dVetaRP);
 452           } //loop over bins be
 453
 454
 455           //v as a function of eta for POI selection
 456           for (Int_t be=1;be<=iNbinsEta;be++)  {
 457             dEtaPOI = fHist2POI[theta]->GetBinCenter(be);
 458             cNumerPOI = fHist2POI[theta]->GetNumerEta(be);
 459             if (cNumerPOI.Rho()==0) {
 460               if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero in Finish()"<<endl;
 461               dReRatioPOI = 0;
 462             }
 463             else if (cDenom.Rho()==0) {
 464               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 465               dReRatioPOI = 0;
 466             }
 467             else {
 468               dReRatioPOI = (cNumerPOI/cDenom).Re();
 469             }
 470             dVetaPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
 471             fHistProVetaPOI->Fill(dEtaPOI,dVetaPOI);
 472           } //loop over bins be
 473
 474
 475           //v as a function of Pt for RP selection
 476           for (Int_t bp=1;bp<=iNbinsPt;bp++)  {
 477             dPtRP = fHist2RP[theta]->GetBinCenterPt(bp);
 478             cNumerRP = fHist2RP[theta]->GetNumerPt(bp);
 479             if (cNumerRP.Rho()==0) {
 480               if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero"<<endl;
 481               dReRatioRP = 0;
 482             }
 483             else if (cDenom.Rho()==0) {
 484               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 485               dReRatioRP = 0;
 486             }
 487             else {
 488               dReRatioRP = (cNumerRP/cDenom).Re();
 489             }
 490             dVPtRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
 491             fHistProVPtRP->Fill(dPtRP,dVPtRP);
 492           } //loop over bins bp
 493
 494
 495
 496           //v as a function of Pt for POI selection
 497           for (Int_t bp=1;bp<=iNbinsPt;bp++)  {
 498             dPtPOI = fHist2POI[theta]->GetBinCenterPt(bp);
 499             cNumerPOI = fHist2POI[theta]->GetNumerPt(bp);
 500             if (cNumerPOI.Rho()==0) {
 501               if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero"<<endl;
 502               dReRatioPOI = 0;
 503             }
 504             else if (cDenom.Rho()==0) {
 505               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 506               dReRatioPOI = 0;
 507             }
 508             else {
 509               dReRatioPOI = (cNumerPOI/cDenom).Re();
 510             }
 511             dVPtPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
 512             fHistProVPtPOI->Fill(dPtPOI,dVPtPOI);
 513           } //loop over bins bp
 514
 515         }
 516       }
 517
 518     }//end of loop over theta
 519
 520     //calculate the average of fVtheta = fV
 521     dV /= iNtheta;
 522
 523     //sigma2 and chi (for statistical error calculations)
 524     Double_t  dSigma2 = 0;
 525     Double_t  dChi= 0;
 526     if (fEventNumber!=0) {
 527       *fQsum /= fEventNumber;
 528       //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
 529       //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
 530       fQ2sum /= fEventNumber;
 531       //cout<<"fQ2sum = "<<fQ2sum<<endl;
 532       dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
 533       if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
 534       else dChi = -1.;
 535       fCommonHistsRes->FillChiRP(dChi);
 536
 537       // recalculate statistical errors on integrated flow
 538       //combining 5 theta angles to 1 relative error BP eq. 89
 539       Double_t dRelErr2comb = 0.;
 540       Int_t iEvts = fEventNumber;
 541       if (iEvts!=0) {
 542         for (Int_t theta=0;theta<iNtheta;theta++){
 543         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 544         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 545                                        TMath::Cos(dTheta));
 546         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 547                                       TMath::Cos(dTheta));
 548         dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
 549                           TMath::BesselJ1(dJ01)))*
 550           (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
 551            dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
 552         }
 553         dRelErr2comb /= iNtheta;
 554       }
 555       Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
 556       Double_t dVErr = dV*dRelErrcomb ;
 557       fCommonHistsRes->FillIntegratedFlowRP(dV,dVErr);
 558
 559       cout<<"*************************************"<<endl;
 560       cout<<"*************************************"<<endl;
 561       cout<<"      Integrated flow from           "<<endl;
 562       cout<<"        Lee-Yang Zeroes              "<<endl;
 563       cout<<endl;
 564       cout<<"dV(RP) = "<<dV<<" +- "<<dVErr<<" and chi = "<<dChi<<endl;
 565     }
 566
 567     //copy content of profile into TH1D and add error and fill the AliFlowCommonHistResults
 568
 569     //v as a function of eta for RP selection
 570     for(Int_t b=0;b<iNbinsEta;b++) {
 571       Double_t dv2pro  = fHistProVetaRP->GetBinContent(b);
 572       Double_t dNprime = fCommonHists->GetEntriesInEtaBinRP(b);
 573       Double_t dErrdifcomb = 0.;  //set error to zero
 574       Double_t dErr2difcomb = 0.; //set error to zero
 575       //calculate error
 576       if (dNprime!=0.) {
 577         for (Int_t theta=0;theta<iNtheta;theta++) {
 578           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 579           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 580                                            TMath::Cos(dTheta));
 581           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 582                                           TMath::Cos(dTheta));
 583           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 584                                                  TMath::BesselJ1(dJ01)))*
 585             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 586              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 587         } //loop over theta
 588       }
 589
 590       if (dErr2difcomb!=0.) {
 591         dErr2difcomb /= iNtheta;
 592         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 593       }
 594       else {dErrdifcomb = 0.;}
 595       //fill TH1D
 596       fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dErrdifcomb);
 597     } //loop over bins b
 598
 599
 600     //v as a function of eta for POI selection
 601     for(Int_t b=0;b<iNbinsEta;b++) {
 602       Double_t dv2pro  = fHistProVetaPOI->GetBinContent(b);
 603       Double_t dNprime = fCommonHists->GetEntriesInEtaBinPOI(b);
 604       Double_t dErrdifcomb = 0.;  //set error to zero
 605       Double_t dErr2difcomb = 0.; //set error to zero
 606       //calculate error
 607       if (dNprime!=0.) {
 608         for (Int_t theta=0;theta<iNtheta;theta++) {
 609           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 610           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 611                                            TMath::Cos(dTheta));
 612           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 613                                           TMath::Cos(dTheta));
 614           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 615                                                TMath::BesselJ1(dJ01)))*
 616             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 617              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 618         } //loop over theta
 619       }
 620
 621       if (dErr2difcomb!=0.) {
 622         dErr2difcomb /= iNtheta;
 623         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 624       }
 625       else {dErrdifcomb = 0.;}
 626       //fill TH1D
 627       fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dErrdifcomb);
 628     } //loop over bins b
 629
 630
 631
 632     //v as a function of Pt for RP selection
 633     for(Int_t b=0;b<iNbinsPt;b++) {
 634       Double_t dv2pro  = fHistProVPtRP->GetBinContent(b);
 635       Double_t dNprime = fCommonHists->GetEntriesInPtBinRP(b);
 636       Double_t dErrdifcomb = 0.;  //set error to zero
 637       Double_t dErr2difcomb = 0.; //set error to zero
 638       //calculate error
 639       if (dNprime!=0.) {
 640         for (Int_t theta=0;theta<iNtheta;theta++) {
 641           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 642           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 643                                            TMath::Cos(dTheta));
 644           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 645                                           TMath::Cos(dTheta));
 646           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 647                                                TMath::BesselJ1(dJ01)))*
 648             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 649              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 650         } //loop over theta
 651       }
 652
 653       if (dErr2difcomb!=0.) {
 654         dErr2difcomb /= iNtheta;
 655         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 656         //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
 657       }
 658       else {dErrdifcomb = 0.;}
 659
 660       //fill TH1D
 661       fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dErrdifcomb);
 662     } //loop over bins b
 663
 664     //v as a function of Pt for POI selection
 665     TH1F* fHistPtDiff = fCommonHists->GetHistPtDiff(); //for calculating integrated flow
 666     Double_t dVPOI = 0.;
 667     Double_t dSum = 0.;
 668     Double_t dErrV =0.;
 669
 670     for(Int_t b=0;b<iNbinsPt;b++) {
 671       Double_t dv2pro = fHistProVPtPOI->GetBinContent(b);
 672       Double_t dNprime = fCommonHists->GetEntriesInPtBinPOI(b);
 673       Double_t dErrdifcomb = 0.;  //set error to zero
 674       Double_t dErr2difcomb = 0.; //set error to zero
 675       //calculate error
 676       if (dNprime!=0.) {
 677         for (Int_t theta=0;theta<iNtheta;theta++) {
 678           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 679           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 680                                            TMath::Cos(dTheta));
 681           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 682                                           TMath::Cos(dTheta));
 683           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 684                                                  TMath::BesselJ1(dJ01)))*
 685             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 686              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 687         } //loop over theta
 688       }
 689
 690       if (dErr2difcomb!=0.) {
 691         dErr2difcomb /= iNtheta;
 692         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 693         //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
 694       }
 695       else {dErrdifcomb = 0.;}
 696
 697       //fill TH1D
 698       fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dErrdifcomb);
 699       //calculate integrated flow for POI selection
 700       if (fHistPtDiff){
 701         Double_t dYieldPt = fHistPtDiff->GetBinContent(b);
 702         dVPOI += dv2pro*dYieldPt;
 703         dSum +=dYieldPt;
 704         dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
 705       } else { cout<<"fHistPtDiff is NULL"<<endl; }
 706     } //loop over bins b
 707
 708     if (dSum != 0.) {
 709       dVPOI /= dSum; //the pt distribution should be normalised
 710       dErrV /= (dSum*dSum);
 711       dErrV = TMath::Sqrt(dErrV);
 712     }
 713     cout<<"dV(POI) is "<<dVPOI<<" +- "<<dErrV<<endl;
 714     cout<<endl;
 715     cout<<"*************************************"<<endl;
 716     cout<<"*************************************"<<endl;
 717     fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
 718
 719   } //secondrun
 720
 721   cout<<"----LYZ analysis finished....----"<<endl<<endl;
 722
 723   return kTRUE;
 724 }
 725
 726
 727 //-----------------------------------------------------------------------
 728
 729  Bool_t AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent(AliFlowEventSimple* anEvent)
 730 {
 731   // Get event quantities from AliFlowEvent for all particles
 732
 733   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl;
 734
 735   if (!anEvent){
 736     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
 737     return kFALSE;
 738   }
 739
 740   //define variables
 741   TComplex cExpo, cGtheta, cGthetaNew, cZ;
 742   Int_t iNtheta = AliFlowLYZConstants::kTheta;
 743   Int_t iNbins = AliFlowLYZConstants::kNbins;
 744
 745
 746   //calculate flow
 747   Double_t dOrder = 2.;
 748
 749   //get the Q vector
 750   AliFlowVector vQ = anEvent->GetQ();
 751   //weight by the multiplicity
 752   Double_t dQX = 0;
 753   Double_t dQY = 0;
 754   if (vQ.GetMult() != 0) {
 755     dQX = vQ.X()/vQ.GetMult();
 756     dQY = vQ.Y()/vQ.GetMult();
 757   }
 758   vQ.Set(dQX,dQY);
 759
 760   //for chi calculation:
 761   *fQsum += vQ;
 762   fHistQsumforChi->SetBinContent(1,fQsum->X());
 763   fHistQsumforChi->SetBinContent(2,fQsum->Y());
 764   fQ2sum += vQ.Mod2();
 765   fHistQsumforChi->SetBinContent(3,fQ2sum);
 766   //cerr<<"fQ2sum = "<<fQ2sum<<endl;
 767
 768   for (Int_t theta=0;theta<iNtheta;theta++)
 769     {
 770       Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
 771
 772       //calculate dQtheta = cos(dOrder*(fPhi-dTheta);the projection of the Q vector on the reference direction dTheta
 773       Double_t dQtheta = GetQtheta(vQ, dTheta);
 774
 775       for (Int_t bin=1;bin<=iNbins;bin++)
 776         {
 777           Double_t dR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram  //FIXED???
 778           //if (theta == 0) cerr<<"cerr::dR = "<<dR<<endl;
 779           if (fUseSum)
 780             {
 781               //calculate the sum generating function
 782               cExpo(0.,dR*dQtheta);                       //Re=0 ; Im=dR*dQtheta
 783               cGtheta = TComplex::Exp(cExpo);
 784             }
 785           else
 786             {
 787               //calculate the product generating function
 788               cGtheta = GetGrtheta(anEvent, dR, dTheta);
 789               if (cGtheta.Rho2() > 100.) break;
 790             }
 791
 792           fHist1[theta]->Fill(dR,cGtheta);              //fill real and imaginary part of cGtheta
 793
 794         } //loop over bins
 795     } //loop over theta
 796
 797
 798   return kTRUE;
 799
 800
 801 }
 802
 803  //-----------------------------------------------------------------------
 804  Bool_t AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent(AliFlowEventSimple* anEvent)
 805 {
 806   //for differential flow
 807
 808   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl;
 809
 810   if (!anEvent){
 811     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
 812     return kFALSE;
 813   }
 814
 815   //define variables
 816   TComplex cExpo, cDenom, cNumerRP, cNumerPOI, cCosTermComplex;
 817   Double_t dPhi, dEta, dPt;
 818   Double_t dR0 = 0.;
 819   Double_t dCosTermRP = 0.;
 820   Double_t dCosTermPOI = 0.;
 821   Double_t m = 1.;
 822   Double_t dOrder = 2.;
 823   Int_t iNtheta = AliFlowLYZConstants::kTheta;
 824
 825
 826   //get the Q vector
 827   AliFlowVector vQ = anEvent->GetQ();
 828   //weight by the multiplicity
 829   Double_t dQX = 0.;
 830   Double_t dQY = 0.;
 831   if (vQ.GetMult() != 0) {
 832     dQX = vQ.X()/vQ.GetMult();
 833     dQY = vQ.Y()/vQ.GetMult();
 834   }
 835   vQ.Set(dQX,dQY);
 836
 837   //for chi calculation:
 838   *fQsum += vQ;
 839   fHistQsumforChi->SetBinContent(1,fQsum->X());
 840   fHistQsumforChi->SetBinContent(2,fQsum->Y());
 841   fQ2sum += vQ.Mod2();
 842   fHistQsumforChi->SetBinContent(3,fQ2sum);
 843
 844   for (Int_t theta=0;theta<iNtheta;theta++)
 845     {
 846       Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
 847
 848       //calculate dQtheta = cos(dOrder*(dPhi-dTheta);the projection of the Q vector on the reference direction dTheta
 849       Double_t dQtheta = GetQtheta(vQ, dTheta);
 850       //cerr<<"dQtheta for fdenom = "<<dQtheta<<endl;
 851
 852       //denominator for differential v
 853       if (fHistProR0theta) {
 854         dR0 = fHistProR0theta->GetBinContent(theta+1);
 855       }
 856       else {
 857         cout <<"pointer fHistProR0theta does not exist" << endl;
 858       }
 859       //cerr<<"dR0 = "<<dR0 <<endl;
 860
 861       if (fUseSum) //sum generating function
 862         {
 863           cExpo(0.,dR0*dQtheta);
 864           cDenom = dQtheta*(TComplex::Exp(cExpo)); //BP eq 12
 865           //loop over tracks in event
 866           Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 867           for (Int_t i=0;i<iNumberOfTracks;i++)  {
 868             AliFlowTrackSimple*  pTrack = anEvent->GetTrack(i);
 869             if (pTrack) {
 870               dEta = pTrack->Eta();
 871               dPt = pTrack->Pt();
 872               dPhi = pTrack->Phi();
 873               if (pTrack->UseForIntegratedFlow()) { // RP selection
 874                 dCosTermRP = cos(m*dOrder*(dPhi-dTheta));
 875                 cNumerRP = dCosTermRP*(TComplex::Exp(cExpo));
 876                 if (cNumerRP.Rho()==0) {cerr<<"WARNING: modulus of cNumerRP is zero in SecondFillFromFlowEvent"<<endl;}
 877                 if (fDebug) cerr<<"modulus of cNumerRP is "<<cNumerRP.Rho()<<endl;
 878                 if (fHist2RP[theta]) {
 879                   fHist2RP[theta]->Fill(dEta,dPt,cNumerRP); }
 880               }
 881               if (pTrack->UseForDifferentialFlow()) { //POI selection
 882                 dCosTermPOI = cos(m*dOrder*(dPhi-dTheta));
 883                 cNumerPOI = dCosTermPOI*(TComplex::Exp(cExpo));
 884                 if (cNumerPOI.Rho()==0) {cerr<<"WARNING: modulus of cNumerPOI is zero in SecondFillFromFlowEvent"<<endl;}
 885                 if (fDebug) cerr<<"modulus of cNumerPOI is "<<cNumerPOI.Rho()<<endl;
 886                 if (fHist2POI[theta]) {
 887                   fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI); }
 888               }
 889               else {
 890                 cout << "fHist2 pointer mising" <<endl;
 891               }
 892             } //if track
 893             else {cerr << "no particle!!!"<<endl;}
 894           } //loop over tracks
 895
 896         } //sum
 897       else        //product generating function
 898         {
 899           cDenom = GetDiffFlow(anEvent, dR0, theta);
 900
 901         }//product
 902       if (fHistProReDenom && fHistProImDenom) {
 903         fHistProReDenom->Fill(theta,cDenom.Re());               //fill the real part of fDenom
 904         fHistProImDenom->Fill(theta,cDenom.Im());               //fill the imaginary part of fDenom
 905       }
 906       else {
 907         cout << "Pointers to cDenom  mising" << endl;
 908       }
 909
 910     }//end of loop over theta
 911
 912   return kTRUE;
 913
 914
 915 }
 916  //-----------------------------------------------------------------------
 917  Double_t AliFlowAnalysisWithLeeYangZeros::GetQtheta(AliFlowVector aQ, Double_t aTheta)
 918 {
 919   //calculate Qtheta. Qtheta is the sum over all particles of cos(dOrder*(dPhi-dTheta)) BP eq. 3
 920   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl;
 921
 922   Double_t dQtheta = 0.;
 923   Double_t dOrder = 2.;
 924
 925   dQtheta = aQ.X()*cos(dOrder*aTheta)+aQ.Y()*sin(dOrder*aTheta);
 926
 927   return dQtheta;
 928
 929 }
 930
 931
 932 //-----------------------------------------------------------------------
 933 TComplex AliFlowAnalysisWithLeeYangZeros::GetGrtheta(AliFlowEventSimple* anEvent, Double_t aR, Double_t aTheta)
 934 {
 935   // Product Generating Function for LeeYangZeros method
 936   // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
 937
 938   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
 939
 940
 941   TComplex cG = TComplex::One();
 942   Double_t dOrder =  2.;
 943   Double_t dWgt = 1.;
 944
 945   Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 946
 947   for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
 948     {
 949       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
 950       if (pTrack){
 951         if (pTrack->UseForIntegratedFlow()) {
 952           Double_t dPhi = pTrack->Phi();
 953           Double_t dGIm = aR * dWgt*cos(dOrder*(dPhi - aTheta));
 954           TComplex cGi(1., dGIm);
 955           cG *= cGi;     //product over all tracks
 956         }
 957       }
 958       else {cerr << "no particle pointer !!!"<<endl;}
 959     }//loop over tracks
 960
 961   return cG;
 962
 963 }
 964
 965
 966 //-----------------------------------------------------------------------
 967 TComplex AliFlowAnalysisWithLeeYangZeros::GetDiffFlow(AliFlowEventSimple* anEvent, Double_t aR0, Int_t theta)
 968 {
 969   // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
 970   // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
 971   // Also for v1 mixed harmonics: DF Eq. 5
 972   // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
 973
 974   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
 975
 976   TComplex cG = TComplex::One();
 977   TComplex cdGr0(0.,0.);
 978   Double_t dOrder =  2.;
 979   Double_t dWgt = 1.;
 980
 981   Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 982
 983   Int_t iNtheta = AliFlowLYZConstants::kTheta;
 984   Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
 985
 986   //for the denominator (use all RP selected particles)
 987   for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
 988     {
 989       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
 990       if (pTrack){
 991         if (pTrack->UseForIntegratedFlow()) {
 992           Double_t dPhi = pTrack->Phi();
 993           Double_t dCosTerm = dWgt*cos(dOrder*(dPhi - dTheta));
 994           //GetGr0theta
 995           Double_t dGIm = aR0 * dCosTerm;
 996           TComplex cGi(1., dGIm);
 997           TComplex cCosTermComplex(1., aR0*dCosTerm);
 998           cG *= cGi;     //product over all tracks
 999           //GetdGr0theta
1000           cdGr0 +=(dCosTerm / cCosTermComplex);  //sum over all tracks
1001         }
1002       } //if particle
1003       else {cerr << "no particle!!!"<<endl;}
1004     }//loop over tracks
1005
1006   //for the numerator
1007   for (Int_t i=0;i<iNumberOfTracks;i++)
1008     {
1009       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1010       if (pTrack){
1011         Double_t dEta = pTrack->Eta();
1012         Double_t dPt = pTrack->Pt();
1013         Double_t dPhi = pTrack->Phi();
1014         Double_t dCosTerm = cos(dOrder*(dPhi-dTheta));
1015         TComplex cCosTermComplex(1.,aR0*dCosTerm);
1016         //RP selection
1017         if (pTrack->UseForIntegratedFlow()) {
1018           TComplex cNumerRP = cG*dCosTerm/cCosTermComplex;  //PG Eq. 9
1019           fHist2RP[theta]->Fill(dEta,dPt,cNumerRP);
1020         }
1021         //POI selection
1022         if (pTrack->UseForDifferentialFlow()) {
1023           TComplex cNumerPOI = cG*dCosTerm/cCosTermComplex;  //PG Eq. 9
1024           fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI);
1025         }
1026       } //if particle
1027       else {cerr << "no particle pointer!!!"<<endl;}
1028     }//loop over tracks
1029
1030   TComplex cDenom = cG*cdGr0;
1031   return cDenom;
1032
1033 }
1034
1035 //-----------------------------------------------------------------------
1036