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