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