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