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