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 = dJ01/(dR0+dBinsize);
 335         Double_t dVmin = dJ01/(dR0-dBinsize);
 336         dv = dVtheta;
 337         Double_t dvplus = dVplus;
 338         Double_t dvmin = dVmin;
 339         if (fDebug) cout<<"dv = "<<dv<<" and dvplus = "<<dvplus<< " and dvmin = "<<dvmin<<endl;
 340
 341         //fill the histograms
 342         fHistProR0theta->Fill(theta,dR0);
 343         fHistProVtheta->Fill(theta,dVtheta);
 344         //get average value of fVtheta = fV
 345         dV += dVtheta;
 346       } //end of loop over theta
 347
 348     //get average value of fVtheta = fV
 349     dV /=iNtheta;
 350
 351     //sigma2 and chi
 352     Double_t  dSigma2 = 0;
 353     Double_t  dChi= 0;
 354     if (fEventNumber!=0) {
 355       *fQsum /= fEventNumber;
 356       fQ2sum /= fEventNumber;
 357       dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
 358       if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
 359       else dChi = -1.;
 360       fCommonHistsRes->FillChiRP(dChi);
 361
 362       cout<<"*************************************"<<endl;
 363       cout<<"*************************************"<<endl;
 364       cout<<"      Integrated flow from           "<<endl;
 365       cout<<"        Lee-Yang Zeroes              "<<endl;
 366       cout<<endl;
 367       cout<<"Chi = "<<dChi<<endl;
 368       cout<<endl;
 369     }
 370
 371     // recalculate statistical errors on integrated flow
 372     //combining 5 theta angles to 1 relative error BP eq. 89
 373     Double_t dRelErr2comb = 0.;
 374     Int_t iEvts = fEventNumber;
 375     if (iEvts!=0) {
 376       for (Int_t theta=0;theta<iNtheta;theta++){
 377         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 378         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 379                                        TMath::Cos(dTheta));
 380         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 381                                       TMath::Cos(dTheta));
 382         dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
 383                           TMath::BesselJ1(dJ01)))*
 384           (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
 385            dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
 386       }
 387       dRelErr2comb /= iNtheta;
 388     }
 389     Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
 390
 391     //copy content of profile into TH1D and add error
 392     Double_t dv2pro = dV;   //in the case that fv is equal to fV
 393     Double_t dv2Err = dv2pro*dRelErrcomb ;
 394     cout<<"dV = "<<dv2pro<<" +- "<<dv2Err<<endl;
 395     cout<<endl;
 396     cout<<"*************************************"<<endl;
 397     cout<<"*************************************"<<endl;
 398     fCommonHistsRes->FillIntegratedFlow(dv2pro, dv2Err);
 399
 400
 401     if (fDebug) cout<<"****histograms filled****"<<endl;
 402
 403     return kTRUE;
 404     fEventNumber =0; //set to zero for second round over events
 405   }  //firstrun
 406
 407
 408   else {  //second run
 409
 410     //calculate differential flow
 411     //declare variables
 412     TComplex cDenom, cNumerRP, cNumerPOI, cDtheta;
 413     Int_t m = 1;
 414     TComplex i = TComplex::I();
 415     Double_t dBesselRatio[3] = {1., 1.202, 2.69};
 416     Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 417     Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
 418
 419     Double_t dEtaRP, dPtRP, dReRatioRP, dVetaRP, dVPtRP, dEtaPOI, dPtPOI, dReRatioPOI, dVetaPOI, dVPtPOI;
 420
 421     Double_t dR0 = 0.;
 422     Double_t dVtheta = 0.;
 423     Double_t dV = 0.;
 424     Double_t dReDenom = 0.;
 425     Double_t dImDenom = 0.;
 426     for (Int_t theta=0;theta<iNtheta;theta++)  { //loop over theta
 427       if (!fHistProR0theta) {
 428         cout << "Hist pointer R0theta in file does not exist" <<endl;
 429       } else {
 430         dR0 = fHistProR0theta->GetBinContent(theta+1); //histogram starts at bin 1
 431         if (fDebug) cerr<<"dR0 = "<<dR0<<endl;
 432         if (dR0!=0) dVtheta = dJ01/dR0;
 433         dV += dVtheta;
 434         // BP Eq. 9  -> Vn^theta = j01/r0^theta
 435         if (!fHistProReDenom && !fHistProImDenom) {
 436           cout << "Hist pointer fDenom in file does not exist" <<endl;
 437         } else {
 438           dReDenom = fHistProReDenom->GetBinContent(theta+1);
 439           dImDenom = fHistProImDenom->GetBinContent(theta+1);
 440           cDenom(dReDenom,dImDenom);
 441
 442           //for new method and use by others (only with the sum generating function):
 443           if (fUseSum) {
 444             cDtheta = dR0*cDenom/dJ01;
 445             fHistProReDtheta->Fill(theta,cDtheta.Re());
 446             fHistProImDtheta->Fill(theta,cDtheta.Im());
 447           }
 448
 449           cDenom *= TComplex::Power(i, m-1);
 450           //cerr<<"TComplex::Power(i, m-1) = "<<TComplex::Power(i, m-1).Rho()<<endl; //checked ok
 451
 452           //v as a function of eta for RP selection
 453           for (Int_t be=1;be<=iNbinsEta;be++)  {
 454             dEtaRP = fHist2RP[theta]->GetBinCenter(be);
 455             cNumerRP = fHist2RP[theta]->GetNumerEta(be);
 456             if (cNumerRP.Rho()==0) {
 457               if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero in Finish()"<<endl;
 458               dReRatioRP = 0;
 459             }
 460             else if (cDenom.Rho()==0) {
 461               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 462               dReRatioRP = 0;
 463             }
 464             else {
 465               dReRatioRP = (cNumerRP/cDenom).Re();
 466             }
 467             dVetaRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
 468             fHistProVetaRP->Fill(dEtaRP,dVetaRP);
 469           } //loop over bins be
 470
 471
 472           //v as a function of eta for POI selection
 473           for (Int_t be=1;be<=iNbinsEta;be++)  {
 474             dEtaPOI = fHist2POI[theta]->GetBinCenter(be);
 475             cNumerPOI = fHist2POI[theta]->GetNumerEta(be);
 476             if (cNumerPOI.Rho()==0) {
 477               if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero in Finish()"<<endl;
 478               dReRatioPOI = 0;
 479             }
 480             else if (cDenom.Rho()==0) {
 481               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 482               dReRatioPOI = 0;
 483             }
 484             else {
 485               dReRatioPOI = (cNumerPOI/cDenom).Re();
 486             }
 487             dVetaPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
 488             fHistProVetaPOI->Fill(dEtaPOI,dVetaPOI);
 489           } //loop over bins be
 490
 491
 492           //v as a function of Pt for RP selection
 493           for (Int_t bp=1;bp<=iNbinsPt;bp++)  {
 494             dPtRP = fHist2RP[theta]->GetBinCenterPt(bp);
 495             cNumerRP = fHist2RP[theta]->GetNumerPt(bp);
 496             if (cNumerRP.Rho()==0) {
 497               if (fDebug) cerr<<"WARNING: modulus of cNumerRP is zero"<<endl;
 498               dReRatioRP = 0;
 499             }
 500             else if (cDenom.Rho()==0) {
 501               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 502               dReRatioRP = 0;
 503             }
 504             else {
 505               dReRatioRP = (cNumerRP/cDenom).Re();
 506             }
 507             dVPtRP = dBesselRatio[m-1]*dReRatioRP*dVtheta; //BP eq. 12
 508             fHistProVPtRP->Fill(dPtRP,dVPtRP);
 509           } //loop over bins bp
 510
 511
 512
 513           //v as a function of Pt for POI selection
 514           for (Int_t bp=1;bp<=iNbinsPt;bp++)  {
 515             dPtPOI = fHist2POI[theta]->GetBinCenterPt(bp);
 516             cNumerPOI = fHist2POI[theta]->GetNumerPt(bp);
 517             if (cNumerPOI.Rho()==0) {
 518               if (fDebug) cerr<<"WARNING: modulus of cNumerPOI is zero"<<endl;
 519               dReRatioPOI = 0;
 520             }
 521             else if (cDenom.Rho()==0) {
 522               if (fDebug) cerr<<"WARNING: modulus of cDenom is zero"<<endl;
 523               dReRatioPOI = 0;
 524             }
 525             else {
 526               dReRatioPOI = (cNumerPOI/cDenom).Re();
 527             }
 528             dVPtPOI = dBesselRatio[m-1]*dReRatioPOI*dVtheta; //BP eq. 12
 529             fHistProVPtPOI->Fill(dPtPOI,dVPtPOI);
 530           } //loop over bins bp
 531
 532         }
 533       }
 534
 535     }//end of loop over theta
 536
 537     //calculate the average of fVtheta = fV
 538     dV /= iNtheta;
 539     if (dV==0.) { cout<<"dV = 0! Leaving LYZ analysis."<<endl; return kFALSE; }
 540
 541     //sigma2 and chi (for statistical error calculations)
 542     Double_t  dSigma2 = 0;
 543     Double_t  dChi= 0;
 544     if (fEventNumber!=0) {
 545       *fQsum /= fEventNumber;
 546       //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
 547       //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
 548       fQ2sum /= fEventNumber;
 549       //cout<<"fQ2sum = "<<fQ2sum<<endl;
 550       dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
 551       if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
 552       else dChi = -1.;
 553       fCommonHistsRes->FillChiRP(dChi);
 554
 555       // recalculate statistical errors on integrated flow
 556       //combining 5 theta angles to 1 relative error BP eq. 89
 557       Double_t dRelErr2comb = 0.;
 558       Int_t iEvts = fEventNumber;
 559       if (iEvts!=0) {
 560         for (Int_t theta=0;theta<iNtheta;theta++){
 561         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 562         Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 563                                        TMath::Cos(dTheta));
 564         Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 565                                       TMath::Cos(dTheta));
 566         dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
 567                           TMath::BesselJ1(dJ01)))*
 568           (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
 569            dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
 570         }
 571         dRelErr2comb /= iNtheta;
 572       }
 573       Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
 574       Double_t dVErr = dV*dRelErrcomb ;
 575       fCommonHistsRes->FillIntegratedFlow(dV,dVErr);
 576
 577       cout<<"*************************************"<<endl;
 578       cout<<"*************************************"<<endl;
 579       cout<<"      Integrated flow from           "<<endl;
 580       cout<<"        Lee-Yang Zeroes              "<<endl;
 581       cout<<endl;
 582       cout<<"Chi = "<<dChi<<endl;
 583       cout<<"dV = "<<dV<<" +- "<<dVErr<<endl;
 584       cout<<endl;
 585     }
 586
 587     //copy content of profile into TH1D and add error and fill the AliFlowCommonHistResults
 588
 589     //v as a function of eta for RP selection
 590     for(Int_t b=0;b<iNbinsEta;b++) {
 591       Double_t dv2pro  = fHistProVetaRP->GetBinContent(b);
 592       Double_t dNprime = fCommonHists->GetEntriesInEtaBinRP(b);
 593       Double_t dErrdifcomb = 0.;  //set error to zero
 594       Double_t dErr2difcomb = 0.; //set error to zero
 595       //calculate error
 596       if (dNprime!=0.) {
 597         for (Int_t theta=0;theta<iNtheta;theta++) {
 598           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 599           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 600                                            TMath::Cos(dTheta));
 601           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 602                                           TMath::Cos(dTheta));
 603           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 604                                                  TMath::BesselJ1(dJ01)))*
 605             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 606              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 607         } //loop over theta
 608       }
 609
 610       if (dErr2difcomb!=0.) {
 611         dErr2difcomb /= iNtheta;
 612         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 613       }
 614       else {dErrdifcomb = 0.;}
 615       //fill TH1D
 616       fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dErrdifcomb);
 617     } //loop over bins b
 618
 619
 620     //v as a function of eta for POI selection
 621     for(Int_t b=0;b<iNbinsEta;b++) {
 622       Double_t dv2pro  = fHistProVetaPOI->GetBinContent(b);
 623       Double_t dNprime = fCommonHists->GetEntriesInEtaBinPOI(b);
 624       Double_t dErrdifcomb = 0.;  //set error to zero
 625       Double_t dErr2difcomb = 0.; //set error to zero
 626       //calculate error
 627       if (dNprime!=0.) {
 628         for (Int_t theta=0;theta<iNtheta;theta++) {
 629           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 630           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 631                                            TMath::Cos(dTheta));
 632           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 633                                           TMath::Cos(dTheta));
 634           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 635                                                TMath::BesselJ1(dJ01)))*
 636             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 637              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 638         } //loop over theta
 639       }
 640
 641       if (dErr2difcomb!=0.) {
 642         dErr2difcomb /= iNtheta;
 643         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 644       }
 645       else {dErrdifcomb = 0.;}
 646       //fill TH1D
 647       fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dErrdifcomb);
 648     } //loop over bins b
 649
 650
 651
 652     //v as a function of Pt for RP selection
 653     TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
 654     Double_t dVRP = 0.;
 655     Double_t dSum = 0.;
 656     Double_t dErrV =0.;
 657     for(Int_t b=0;b<iNbinsPt;b++) {
 658       Double_t dv2pro  = fHistProVPtRP->GetBinContent(b);
 659       Double_t dNprime = fCommonHists->GetEntriesInPtBinRP(b);
 660       Double_t dErrdifcomb = 0.;  //set error to zero
 661       Double_t dErr2difcomb = 0.; //set error to zero
 662       //calculate error
 663       if (dNprime!=0.) {
 664         for (Int_t theta=0;theta<iNtheta;theta++) {
 665           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 666           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 667                                            TMath::Cos(dTheta));
 668           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 669                                           TMath::Cos(dTheta));
 670           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 671                                                TMath::BesselJ1(dJ01)))*
 672             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 673              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 674         } //loop over theta
 675       }
 676
 677       if (dErr2difcomb!=0.) {
 678         dErr2difcomb /= iNtheta;
 679         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 680         //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
 681       }
 682       else {dErrdifcomb = 0.;}
 683
 684       //fill TH1D
 685       fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dErrdifcomb);
 686       //calculate integrated flow for RP selection
 687       if (fHistPtRP){
 688         Double_t dYieldPt = fHistPtRP->GetBinContent(b);
 689         dVRP += dv2pro*dYieldPt;
 690         dSum +=dYieldPt;
 691         dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
 692       } else { cout<<"fHistPtRP is NULL"<<endl; }
 693     } //loop over bins b
 694
 695     if (dSum != 0.) {
 696       dVRP /= dSum; //the pt distribution should be normalised
 697       dErrV /= (dSum*dSum);
 698       dErrV = TMath::Sqrt(dErrV);
 699     }
 700     cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
 701     //cout<<endl;
 702     fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
 703
 704
 705     //v as a function of Pt for POI selection
 706     TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
 707     Double_t dVPOI = 0.;
 708     dSum = 0.;
 709     dErrV =0.;
 710
 711     for(Int_t b=0;b<iNbinsPt;b++) {
 712       Double_t dv2pro = fHistProVPtPOI->GetBinContent(b);
 713       Double_t dNprime = fCommonHists->GetEntriesInPtBinPOI(b);
 714       Double_t dErrdifcomb = 0.;  //set error to zero
 715       Double_t dErr2difcomb = 0.; //set error to zero
 716       //calculate error
 717       if (dNprime!=0.) {
 718         for (Int_t theta=0;theta<iNtheta;theta++) {
 719           Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 720           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 721                                            TMath::Cos(dTheta));
 722           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 723                                           TMath::Cos(dTheta));
 724           dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
 725                                                  TMath::BesselJ1(dJ01)))*
 726             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 727              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 728         } //loop over theta
 729       }
 730
 731       if (dErr2difcomb!=0.) {
 732         dErr2difcomb /= iNtheta;
 733         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 734         //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
 735       }
 736       else {dErrdifcomb = 0.;}
 737
 738       //fill TH1D
 739       fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dErrdifcomb);
 740       //calculate integrated flow for POI selection
 741       if (fHistPtPOI){
 742         Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
 743         dVPOI += dv2pro*dYieldPt;
 744         dSum +=dYieldPt;
 745         dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
 746       } else { cout<<"fHistPtPOI is NULL"<<endl; }
 747     } //loop over bins b
 748
 749     if (dSum != 0.) {
 750       dVPOI /= dSum; //the pt distribution should be normalised
 751       dErrV /= (dSum*dSum);
 752       dErrV = TMath::Sqrt(dErrV);
 753     }
 754     cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
 755     cout<<endl;
 756     cout<<"*************************************"<<endl;
 757     cout<<"*************************************"<<endl;
 758     fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
 759
 760   } //secondrun
 761
 762   //cout<<"----LYZ analysis finished....----"<<endl<<endl;
 763
 764   return kTRUE;
 765 }
 766
 767
 768 //-----------------------------------------------------------------------
 769
 770  Bool_t AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent(AliFlowEventSimple* anEvent)
 771 {
 772   // Get event quantities from AliFlowEvent for all particles
 773
 774   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl;
 775
 776   if (!anEvent){
 777     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
 778     return kFALSE;
 779   }
 780
 781   //define variables
 782   TComplex cExpo, cGtheta, cGthetaNew, cZ;
 783   Int_t iNtheta = AliFlowLYZConstants::kTheta;
 784   Int_t iNbins = AliFlowLYZConstants::kNbins;
 785
 786
 787   //calculate flow
 788   Double_t dOrder = 2.;
 789
 790   //get the Q vector
 791   AliFlowVector vQ = anEvent->GetQ();
 792   //weight by the multiplicity
 793   Double_t dQX = 0;
 794   Double_t dQY = 0;
 795   if (vQ.GetMult() != 0) {
 796     dQX = vQ.X()/vQ.GetMult();
 797     dQY = vQ.Y()/vQ.GetMult();
 798   }
 799   vQ.Set(dQX,dQY);
 800
 801   //for chi calculation:
 802   *fQsum += vQ;
 803   fHistQsumforChi->SetBinContent(1,fQsum->X());
 804   fHistQsumforChi->SetBinContent(2,fQsum->Y());
 805   fQ2sum += vQ.Mod2();
 806   fHistQsumforChi->SetBinContent(3,fQ2sum);
 807   //cerr<<"fQ2sum = "<<fQ2sum<<endl;
 808
 809   for (Int_t theta=0;theta<iNtheta;theta++)
 810     {
 811       Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
 812
 813       //calculate dQtheta = cos(dOrder*(fPhi-dTheta);the projection of the Q vector on the reference direction dTheta
 814       Double_t dQtheta = GetQtheta(vQ, dTheta);
 815
 816       for (Int_t bin=1;bin<=iNbins;bin++)
 817         {
 818           Double_t dR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram  //FIXED???
 819           //if (theta == 0) cerr<<"cerr::dR = "<<dR<<endl;
 820           if (fUseSum)
 821             {
 822               //calculate the sum generating function
 823               cExpo(0.,dR*dQtheta);                       //Re=0 ; Im=dR*dQtheta
 824               cGtheta = TComplex::Exp(cExpo);
 825             }
 826           else
 827             {
 828               //calculate the product generating function
 829               cGtheta = GetGrtheta(anEvent, dR, dTheta);
 830               if (cGtheta.Rho2() > 100.) break;
 831             }
 832
 833           fHist1[theta]->Fill(dR,cGtheta);              //fill real and imaginary part of cGtheta
 834
 835         } //loop over bins
 836     } //loop over theta
 837
 838
 839   return kTRUE;
 840
 841
 842 }
 843
 844  //-----------------------------------------------------------------------
 845  Bool_t AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent(AliFlowEventSimple* anEvent)
 846 {
 847   //for differential flow
 848
 849   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl;
 850
 851   if (!anEvent){
 852     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
 853     return kFALSE;
 854   }
 855
 856   //define variables
 857   TComplex cExpo, cDenom, cNumerRP, cNumerPOI, cCosTermComplex;
 858   Double_t dPhi, dEta, dPt;
 859   Double_t dR0 = 0.;
 860   Double_t dCosTermRP = 0.;
 861   Double_t dCosTermPOI = 0.;
 862   Double_t m = 1.;
 863   Double_t dOrder = 2.;
 864   Int_t iNtheta = AliFlowLYZConstants::kTheta;
 865
 866
 867   //get the Q vector
 868   AliFlowVector vQ = anEvent->GetQ();
 869   //weight by the multiplicity
 870   Double_t dQX = 0.;
 871   Double_t dQY = 0.;
 872   if (vQ.GetMult() != 0) {
 873     dQX = vQ.X()/vQ.GetMult();
 874     dQY = vQ.Y()/vQ.GetMult();
 875   }
 876   vQ.Set(dQX,dQY);
 877
 878   //for chi calculation:
 879   *fQsum += vQ;
 880   fHistQsumforChi->SetBinContent(1,fQsum->X());
 881   fHistQsumforChi->SetBinContent(2,fQsum->Y());
 882   fQ2sum += vQ.Mod2();
 883   fHistQsumforChi->SetBinContent(3,fQ2sum);
 884
 885   for (Int_t theta=0;theta<iNtheta;theta++)
 886     {
 887       Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
 888
 889       //calculate dQtheta = cos(dOrder*(dPhi-dTheta);the projection of the Q vector on the reference direction dTheta
 890       Double_t dQtheta = GetQtheta(vQ, dTheta);
 891       //cerr<<"dQtheta for fdenom = "<<dQtheta<<endl;
 892
 893       //denominator for differential v
 894       if (fHistProR0theta) {
 895         dR0 = fHistProR0theta->GetBinContent(theta+1);
 896       }
 897       else {
 898         cout <<"pointer fHistProR0theta does not exist" << endl;
 899       }
 900       //cerr<<"dR0 = "<<dR0 <<endl;
 901
 902       if (fUseSum) //sum generating function
 903         {
 904           cExpo(0.,dR0*dQtheta);
 905           cDenom = dQtheta*(TComplex::Exp(cExpo)); //BP eq 12
 906           //loop over tracks in event
 907           Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 908           for (Int_t i=0;i<iNumberOfTracks;i++)  {
 909             AliFlowTrackSimple*  pTrack = anEvent->GetTrack(i);
 910             if (pTrack) {
 911               dEta = pTrack->Eta();
 912               dPt = pTrack->Pt();
 913               dPhi = pTrack->Phi();
 914               if (pTrack->InRPSelection()) { // RP selection
 915                 dCosTermRP = cos(m*dOrder*(dPhi-dTheta));
 916                 cNumerRP = dCosTermRP*(TComplex::Exp(cExpo));
 917                 if (cNumerRP.Rho()==0) {cerr<<"WARNING: modulus of cNumerRP is zero in SecondFillFromFlowEvent"<<endl;}
 918                 if (fDebug) cerr<<"modulus of cNumerRP is "<<cNumerRP.Rho()<<endl;
 919                 if (fHist2RP[theta]) {
 920                   fHist2RP[theta]->Fill(dEta,dPt,cNumerRP); }
 921               }
 922               if (pTrack->InPOISelection()) { //POI selection
 923                 dCosTermPOI = cos(m*dOrder*(dPhi-dTheta));
 924                 cNumerPOI = dCosTermPOI*(TComplex::Exp(cExpo));
 925                 if (cNumerPOI.Rho()==0) {cerr<<"WARNING: modulus of cNumerPOI is zero in SecondFillFromFlowEvent"<<endl;}
 926                 if (fDebug) cerr<<"modulus of cNumerPOI is "<<cNumerPOI.Rho()<<endl;
 927                 if (fHist2POI[theta]) {
 928                   fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI); }
 929               }
 930               //              else {
 931               //                cout << "fHist2 pointer mising" <<endl;
 932               //              }
 933             } //if track
 934             else {cerr << "no particle!!!"<<endl;}
 935           } //loop over tracks
 936
 937         } //sum
 938       else        //product generating function
 939         {
 940           cDenom = GetDiffFlow(anEvent, dR0, theta);
 941
 942         }//product
 943       if (fHistProReDenom && fHistProImDenom) {
 944         fHistProReDenom->Fill(theta,cDenom.Re());               //fill the real part of fDenom
 945         fHistProImDenom->Fill(theta,cDenom.Im());               //fill the imaginary part of fDenom
 946       }
 947       else {
 948         cout << "Pointers to cDenom  mising" << endl;
 949       }
 950
 951     }//end of loop over theta
 952
 953   return kTRUE;
 954
 955
 956 }
 957  //-----------------------------------------------------------------------
 958  Double_t AliFlowAnalysisWithLeeYangZeros::GetQtheta(AliFlowVector aQ, Double_t aTheta)
 959 {
 960   //calculate Qtheta. Qtheta is the sum over all particles of cos(dOrder*(dPhi-dTheta)) BP eq. 3
 961   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl;
 962
 963   Double_t dQtheta = 0.;
 964   Double_t dOrder = 2.;
 965
 966   dQtheta = aQ.X()*cos(dOrder*aTheta)+aQ.Y()*sin(dOrder*aTheta);
 967
 968   return dQtheta;
 969
 970 }
 971
 972
 973 //-----------------------------------------------------------------------
 974 TComplex AliFlowAnalysisWithLeeYangZeros::GetGrtheta(AliFlowEventSimple* const anEvent, Double_t aR, Double_t aTheta)
 975 {
 976   // Product Generating Function for LeeYangZeros method
 977   // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
 978
 979   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
 980
 981
 982   TComplex cG = TComplex::One();
 983   Double_t dOrder =  2.;
 984   Double_t dWgt = 1./anEvent->GetEventNSelTracksRP(); //weight with the multiplicity
 985
 986   Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 987
 988   for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
 989     {
 990       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
 991       if (pTrack){
 992         if (pTrack->InRPSelection()) {
 993           Double_t dPhi = pTrack->Phi();
 994           Double_t dGIm = aR * dWgt*cos(dOrder*(dPhi - aTheta));
 995           TComplex cGi(1., dGIm);
 996           cG *= cGi;     //product over all tracks
 997         }
 998       }
 999       else {cerr << "no particle pointer !!!"<<endl;}
1000     }//loop over tracks
1001
1002   return cG;
1003
1004 }
1005
1006
1007 //-----------------------------------------------------------------------
1008 TComplex AliFlowAnalysisWithLeeYangZeros::GetDiffFlow(AliFlowEventSimple* const anEvent, Double_t aR0, Int_t theta)
1009 {
1010   // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
1011   // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1012   // Also for v1 mixed harmonics: DF Eq. 5
1013   // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
1014
1015   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
1016
1017   TComplex cG = TComplex::One();
1018   TComplex cdGr0(0.,0.);
1019   Double_t dOrder =  2.;
1020   Double_t dWgt = 1.;
1021
1022   Int_t iNumberOfTracks = anEvent->NumberOfTracks();
1023
1024   Int_t iNtheta = AliFlowLYZConstants::kTheta;
1025   Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi()/dOrder;
1026
1027   //for the denominator (use all RP selected particles)
1028   for (Int_t i=0;i<iNumberOfTracks;i++) //loop over tracks in event
1029     {
1030       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1031       if (pTrack){
1032         if (pTrack->InRPSelection()) {
1033           Double_t dPhi = pTrack->Phi();
1034           Double_t dCosTerm = dWgt*cos(dOrder*(dPhi - dTheta));
1035           //GetGr0theta
1036           Double_t dGIm = aR0 * dCosTerm;
1037           TComplex cGi(1., dGIm);
1038           TComplex cCosTermComplex(1., aR0*dCosTerm);
1039           cG *= cGi;     //product over all tracks
1040           //GetdGr0theta
1041           cdGr0 +=(dCosTerm / cCosTermComplex);  //sum over all tracks
1042         }
1043       } //if particle
1044       else {cerr << "no particle!!!"<<endl;}
1045     }//loop over tracks
1046
1047   //for the numerator
1048   for (Int_t i=0;i<iNumberOfTracks;i++)
1049     {
1050       AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
1051       if (pTrack){
1052         Double_t dEta = pTrack->Eta();
1053         Double_t dPt = pTrack->Pt();
1054         Double_t dPhi = pTrack->Phi();
1055         Double_t dCosTerm = cos(dOrder*(dPhi-dTheta));
1056         TComplex cCosTermComplex(1.,aR0*dCosTerm);
1057         //RP selection
1058         if (pTrack->InRPSelection()) {
1059           TComplex cNumerRP = cG*dCosTerm/cCosTermComplex;  //PG Eq. 9
1060           fHist2RP[theta]->Fill(dEta,dPt,cNumerRP);
1061         }
1062         //POI selection
1063         if (pTrack->InPOISelection()) {
1064           TComplex cNumerPOI = cG*dCosTerm/cCosTermComplex;  //PG Eq. 9
1065           fHist2POI[theta]->Fill(dEta,dPt,cNumerPOI);
1066         }
1067       } //if particle
1068       else {cerr << "no particle pointer!!!"<<endl;}
1069     }//loop over tracks
1070
1071   TComplex cDenom = cG*cdGr0;
1072   return cDenom;
1073
1074 }
1075
1076 //-----------------------------------------------------------------------
1077