PWG2/FLOW/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 #include "Riostream.h"                 //needed as include
  18 #include "AliFlowCommonConstants.h"    //needed as include
  19 #include "AliFlowLYZConstants.h"       //needed as include
  20 #include "AliFlowAnalysisWithLeeYangZeros.h"
  21 #include "AliFlowLYZHist1.h"           //needed as include
  22 #include "AliFlowLYZHist2.h"           //needed as include
  23 #include "AliFlowCommonHist.h"         //needed as include
  24 #include "AliFlowCommonHistResults.h"  //needed as include
  25 #include "AliFlowEventSimple.h"        //needed as include
  26 #include "AliFlowTrackSimple.h"        //needed as include
  27
  28 class AliFlowVector;
  29
  30 #include "TMath.h" //needed as include
  31 #include "TFile.h" //needed as include
  32
  33 class TComplex;
  34 class TProfile;
  35 class TH1F;
  36 class TH1D;
  37
  38
  39 //Description: Maker to analyze Flow using the LeeYangZeros method
  40 //             Equation numbers are from Big Paper (BP): Nucl. Phys. A 727, 373 (2003)
  41 //             Practical Guide (PG):    J. Phys. G: Nucl. Part. Phys. 30, S1213 (2004)
  42 //             Adapted from StFlowLeeYangZerosMaker.cxx
  43 //             by Markus Oldenberg and Art Poskanzer, LBNL
  44 //             with advice from Jean-Yves Ollitrault and Nicolas Borghini
  45 //
  46 //Author: Naomi van der Kolk (kolk@nikhef.nl)
  47
  48
  49
  50 ClassImp(AliFlowAnalysisWithLeeYangZeros)
  51
  52   //-----------------------------------------------------------------------
  53
  54   AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros():
  55     fQ2sum(0),
  56     fQtheta(0),
  57     fEventNumber(0),
  58     fMult(0),
  59     fNbins(0),
  60     fTheta(0),
  61     fEvent(0),
  62     fTrack(0),
  63     fFirstRun(kTRUE),
  64     fUseSum(kTRUE),
  65     fDoubleLoop(kFALSE),
  66     fDebug(kFALSE),
  67     fHistFileName(0),
  68     fHistFile(0),
  69     fSummaryFile(0),
  70     firstRunFileName(0),
  71     firstRunFile(0),
  72     fHistProVtheta(0),
  73     fHistProVeta(0),
  74     fHistProVPt(0),
  75     fHistProR0theta(0),
  76     fHistProReDenom(0),
  77     fHistProImDenom(0),
  78     fHistProReDtheta(0),
  79     fHistProImDtheta(0),
  80     fCommonHists(0),
  81     fCommonHistsRes(0)
  82
  83 {
  84   //default constructor
  85   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::AliFlowAnalysisWithLeeYangZeros default constructor****"<<endl;
  86
  87   // output file (histograms)
  88   TString fHistFileName = "outputFromLYZAnalysis.root" ;
  89
  90   for(Int_t i = 0;i<5;i++)
  91     {
  92       fHist1[i]=0;
  93       fHist2[i]=0;
  94     }
  95
  96   fQ.Set(0.,0.);
  97   fQ.SetMult(0);
  98   fQsum.Set(0.,0.);
  99
 100 }
 101
 102 //-----------------------------------------------------------------------
 103
 104
 105  AliFlowAnalysisWithLeeYangZeros::~AliFlowAnalysisWithLeeYangZeros()
 106  {
 107    //default destructor
 108    if (fDebug) cout<<"****~AliFlowAnalysisWithLeeYangZeros****"<<endl;
 109    delete fHistFile;
 110  }
 111
 112  //-----------------------------------------------------------------------
 113
 114 Bool_t AliFlowAnalysisWithLeeYangZeros::Init()
 115 {
 116   //init method
 117   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Init()****"<<endl;
 118
 119   // Open output files (->plots)
 120   fHistFile = new TFile(fHistFileName.Data(), "RECREATE");
 121
 122
 123   // Book histograms
 124   Int_t fNtheta = AliFlowLYZConstants::kTheta;
 125   Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 126   Int_t fNbinsEta = AliFlowCommonConstants::GetNbinsEta();
 127
 128   Double_t  fPtMin = AliFlowCommonConstants::GetPtMin();
 129   Double_t  fPtMax = AliFlowCommonConstants::GetPtMax();
 130   Double_t  fEtaMin = AliFlowCommonConstants::GetEtaMin();
 131   Double_t  fEtaMax = AliFlowCommonConstants::GetEtaMax();
 132
 133
 134   //for control histograms
 135   fCommonHists = new AliFlowCommonHist("LYZ");
 136   fCommonHistsRes = new AliFlowCommonHistResults("LYZ");
 137
 138   //for first loop over events
 139   if (fFirstRun){
 140     fHistProR0theta  = new TProfile("First_FlowPro_r0theta_LYZ","First_FlowPro_r0theta_LYZ",fNtheta,-0.5,fNtheta-0.5);
 141     fHistProR0theta->SetXTitle("#theta");
 142     fHistProR0theta->SetYTitle("r_{0}^{#theta}");
 143
 144     fHistProVtheta  = new TProfile("First_FlowPro_Vtheta_LYZ","First_FlowPro_Vtheta_LYZ",fNtheta,-0.5,fNtheta-0.5);
 145     fHistProVtheta->SetXTitle("#theta");
 146     fHistProVtheta->SetYTitle("V_{n}^{#theta}");
 147
 148     //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta, fHistProR0theta
 149     for (Int_t theta=0;theta<fNtheta;theta++) {
 150       fHist1[theta]=new AliFlowLYZHist1(theta);
 151     }
 152
 153   }
 154   //for second loop over events
 155   else {
 156     fHistProReDenom = new TProfile("Second_FlowPro_ReDenom_LYZ","Second_FlowPro_ReDenom_LYZ" , fNtheta, -0.5, fNtheta-0.5);
 157     fHistProReDenom->SetXTitle("#theta");
 158     fHistProReDenom->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
 159
 160     fHistProImDenom = new TProfile("Second_FlowPro_ImDenom_LYZ","Second_FlowPro_ImDenom_LYZ" , fNtheta, -0.5, fNtheta-0.5);
 161     fHistProImDenom->SetXTitle("#theta");
 162     fHistProImDenom->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
 163
 164     fHistProVeta = new TProfile("Second_FlowPro_Veta_LYZ","Second_FlowPro_Veta_LYZ",fNbinsEta,fEtaMin,fEtaMax);
 165     fHistProVeta->SetXTitle("rapidity");
 166     fHistProVeta->SetYTitle("v (%)");
 167
 168     fHistProVPt = new TProfile("Second_FlowPro_VPt_LYZ","Second_FlowPro_VPt_LYZ",fNbinsPt,fPtMin,fPtMax);
 169     fHistProVPt->SetXTitle("Pt");
 170     fHistProVPt->SetYTitle("v (%)");
 171
 172     fHistProReDtheta = new TProfile("Second_FlowPro_ReDtheta_LYZ","Second_FlowPro_ReDtheta_LYZ",fNtheta, -0.5, fNtheta-0.5);
 173     fHistProReDtheta->SetXTitle("#theta");
 174     fHistProReDtheta->SetYTitle("Re(D^{#theta})");
 175
 176     fHistProImDtheta = new TProfile("Second_FlowPro_ImDtheta_LYZ","Second_FlowPro_ImDtheta_LYZ",fNtheta, -0.5, fNtheta-0.5);
 177     fHistProImDtheta->SetXTitle("#theta");
 178     fHistProImDtheta->SetYTitle("Im(D^{#theta})");
 179
 180     //class AliFlowLYZHist2 defines the histograms:
 181     for (Int_t theta=0;theta<fNtheta;theta++)  {
 182       fHist2[theta]=new AliFlowLYZHist2(theta);
 183     }
 184
 185     //read hists from first run file
 186     //firstRunFile = new TFile("fof_flowLYZAnal_firstrun.root","READ");  //default is read
 187     if (firstRunFile->IsZombie()){ //check if file exists
 188       cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
 189       exit(-1);
 190     } else if (firstRunFile->IsOpen()){
 191       cout<<"----firstRunFile is open----"<<endl<<endl;
 192       fHistProR0theta  = (TProfile*)firstRunFile->Get("First_FlowPro_r0theta_LYZ");
 193     }
 194   }
 195
 196
 197   if (fDebug) cout<<"****Histograms initialised****"<<endl;
 198
 199   fEventNumber = 0; //set event counter to zero
 200
 201   return kTRUE;
 202 }
 203
 204  //-----------------------------------------------------------------------
 205
 206 Bool_t AliFlowAnalysisWithLeeYangZeros::Make(AliFlowEventSimple* fEvent)
 207 {
 208   //make method
 209   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Make()****"<<endl;
 210
 211   //get tracks from event
 212   if (fEvent) {
 213     if (fFirstRun){
 214       fCommonHists->FillControlHistograms(fEvent);
 215       FillFromFlowEvent(fEvent);
 216     }
 217     else {
 218       fCommonHists->FillControlHistograms(fEvent);
 219       SecondFillFromFlowEvent(fEvent);
 220     }
 221   }
 222
 223   else {
 224     cout<<"##### FlowLeeYangZero: Stack pointer null"<<endl;
 225     return kFALSE;
 226   }
 227   cout<<"^^^^read event "<<fEventNumber<<endl;
 228   fEventNumber++;
 229
 230
 231   return kTRUE;
 232 }
 233
 234
 235   //-----------------------------------------------------------------------
 236  Bool_t AliFlowAnalysisWithLeeYangZeros::Finish()
 237 {
 238   //finish method
 239   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::Finish()****"<<endl;
 240
 241   //define variables for both runs
 242   Double_t  fJ01 = 2.405;
 243   Int_t fNtheta = AliFlowLYZConstants::kTheta;
 244
 245   if (fFirstRun){
 246     Double_t  fR0 = 0;
 247     Double_t  fVtheta = 0;
 248     Double_t  fv = 0;
 249     Double_t  fV = 0;
 250     for (Int_t theta=0;theta<fNtheta;theta++)
 251       {
 252         //get the first minimum r0
 253         fHist1[theta]->FillGtheta();
 254         fR0 = fHist1[theta]->GetR0();
 255
 256         //calculate integrated flow
 257         if (fR0!=0) fVtheta = fJ01/fR0;
 258         //for estimating systematic error resulting from r0
 259         Double_t binsize = (AliFlowLYZConstants::fgMax)/(AliFlowLYZConstants::kNbins);
 260         Double_t fVplus = fJ01/(fR0+binsize);
 261         Double_t fVmin = fJ01/(fR0-binsize);
 262         fv = fVtheta;
 263         Double_t fvplus = fVplus;
 264         Double_t fvmin = fVmin;
 265         cout<<"fv = "<<fv<<" and fvplus = "<<fvplus<< " and fvmin = "<<fvmin<<endl;
 266
 267         //fill the histograms
 268         fHistProR0theta->Fill(theta,fR0);
 269         fHistProVtheta->Fill(theta,fVtheta);
 270         //get average value of fVtheta = fV
 271         fV += fVtheta;
 272       } //end of loop over theta
 273
 274     //get average value of fVtheta = fV
 275     fV /=fNtheta;
 276
 277     //sigma2 and chi
 278     Double_t  fSigma2 = 0;
 279     Double_t  fChi= 0;
 280     if (fEventNumber!=0) {
 281       fQsum /= fEventNumber;
 282       fQ2sum /= fEventNumber;
 283       fSigma2 = fQ2sum - TMath::Power(fQsum.X(),2.) - TMath::Power(fQsum.Y(),2.) - TMath::Power(fV,2.);  //BP eq. 62
 284       if (fSigma2>0) fChi = fV/TMath::Sqrt(fSigma2);
 285       else fChi = -1.;
 286       fCommonHistsRes->FillChi(fChi);
 287       cout<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
 288     }
 289
 290     // recalculate statistical errors on integrated flow
 291     //combining 5 theta angles to 1 relative error BP eq. 89
 292     Double_t fRelErr2comb = 0.;
 293     Int_t nEvts = fEventNumber;
 294     for (Int_t theta=0;theta<fNtheta;theta++){
 295       Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi();
 296       Double_t fApluscomb = TMath::Exp((fJ01*fJ01)/(2*fChi*fChi)*
 297                                        TMath::Cos(fTheta));
 298       Double_t fAmincomb = TMath::Exp(-(fJ01*fJ01)/(2*fChi*fChi)*
 299                                       TMath::Cos(fTheta));
 300       fRelErr2comb += (1/(2*nEvts*(fJ01*fJ01)*TMath::BesselJ1(fJ01)*
 301                           TMath::BesselJ1(fJ01)))*
 302         (fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2)) +
 303          fAmincomb*TMath::BesselJ0(2*fJ01*TMath::Cos(fTheta/2)));
 304     }
 305     fRelErr2comb /= fNtheta;
 306     Double_t fRelErrcomb = TMath::Sqrt(fRelErr2comb);
 307     cout<<"fRelErrcomb = "<<fRelErrcomb<<endl;
 308
 309     //copy content of profile into TH1D and add error
 310     Double_t fv2pro = fV;   //in the case that fv is equal to fV
 311     Double_t fv2Err = fv2pro*fRelErrcomb ;
 312     cout<<"fv2pro +- fv2Err = "<<fv2pro<<" +- "<<fv2Err<<endl;
 313     fCommonHistsRes->FillIntegratedFlow(fv2pro, fv2Err);
 314
 315
 316     if (fDebug) cout<<"****histograms filled****"<<endl;
 317
 318     //save histograms in file //temp for testing selector
 319     fHistFile->cd();
 320     fHistFile->Write();
 321
 322     return kTRUE;
 323     fEventNumber =0; //set to zero for second round over events
 324   }  //firstrun
 325
 326
 327   else {  //second run
 328
 329     //calculate differential flow
 330     //declare variables
 331     TComplex fDenom, fNumer, fDtheta;
 332     Int_t m = 1;
 333     TComplex i = TComplex::I();
 334     Double_t fBesselRatio[3] = {1., 1.202, 2.69};
 335     Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 336     Int_t fNbinsEta = AliFlowCommonConstants::GetNbinsEta();
 337
 338     Double_t fEta, fPt, fReRatio, fVeta, fVPt;
 339
 340
 341     Double_t fR0 = 0.;
 342     Double_t fVtheta = 0.;
 343     Double_t fV = 0.;
 344     Double_t fReDenom = 0.;
 345     Double_t fImDenom = 0.;
 346     for (Int_t theta=0;theta<fNtheta;theta++)  { //loop over theta
 347       if (fHistProR0theta) {
 348         fR0 = fHistProR0theta->GetBinContent(theta+1);
 349         if (fDebug) cerr<<"fR0 = "<<fR0<<endl;
 350         if (fR0!=0) fVtheta = fJ01/fR0;
 351         fV += fVtheta;
 352         // BP Eq. 9  -> Vn^theta = j01/r0^theta
 353         if (fHistProReDenom && fHistProImDenom) {
 354           fReDenom = fHistProReDenom->GetBinContent(theta+1);
 355           fImDenom = fHistProImDenom->GetBinContent(theta+1);
 356         }
 357         else {
 358           cout << "Hist pointer fDenom in file does not exist" <<endl;
 359         }
 360
 361       }
 362       else {
 363         cout << "Hist pointer R0theta in file does not exist" <<endl;
 364       }
 365       //} //loop over theta
 366
 367       fDenom(fReDenom,fImDenom);
 368
 369
 370       //for new method and use by others (only with the sum generating function):
 371       if (fUseSum) {
 372         fR0 = fHistProR0theta->GetBinContent(theta+1);
 373         fDtheta = fR0*fDenom/fJ01;
 374         fHistProReDtheta->Fill(theta,fDtheta.Re());
 375         fHistProImDtheta->Fill(theta,fDtheta.Im());
 376       }
 377
 378       fDenom *= TComplex::Power(i, m-1);
 379       //cerr<<"TComplex::Power(i, m-1) = "<<TComplex::Power(i, m-1).Rho()<<endl; //checked ok
 380
 381       //v as a function of eta
 382       for (Int_t be=1;be<=fNbinsEta;be++)  {
 383         fEta = fHist2[theta]->GetBinCenter(be);
 384         fNumer = fHist2[theta]->GetfNumer(be);
 385         if (fNumer.Rho()==0) {
 386           if (fDebug) cerr<<"WARNING: modulus of fNumer is zero in Finish()"<<endl;
 387           fReRatio = 0;
 388         }
 389         else if (fDenom.Rho()==0) {
 390           if (fDebug) cerr<<"WARNING: modulus of fDenom is zero"<<endl;
 391           fReRatio = 0;
 392         }
 393         else {
 394           //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
 395           fReRatio = (fNumer/fDenom).Re();
 396         }
 397
 398         fVeta = fBesselRatio[m-1]*fReRatio*fVtheta*100.; //BP eq. 12
 399         //if ( j==1 && theta==0) cerr<<"eta = "<<fEta<<" cerr::fReRatio for eta = "<<fReRatio<<" cerr::fVeta for eta = "<<fVeta<<endl;
 400
 401         fHistProVeta->Fill(fEta,fVeta);
 402       } //loop over bins be
 403
 404       //v as a function of Pt
 405       Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 406       for (Int_t bp=1;bp<=fNbinsPt;bp++)  {
 407         fPt = fHist2[theta]->GetBinCenterPt(bp);
 408         fNumer = fHist2[theta]->GetfNumerPt(bp);
 409         if (fNumer.Rho()==0) {
 410           if (fDebug) cerr<<"modulus of fNumer is zero"<<endl;
 411           fReRatio = 0;
 412         }
 413         else if (fDenom.Rho()==0) {
 414           if (fDebug) cerr<<"modulus of fDenom is zero"<<endl;
 415           fReRatio = 0;
 416         }
 417         else {
 418           //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
 419           fReRatio = (fNumer/fDenom).Re();
 420         }
 421
 422         fVPt = fBesselRatio[m-1]*fReRatio*fVtheta*100.; //BP eq. 12
 423
 424         fHistProVPt->Fill(fPt,fVPt);
 425       } //loop over bins bp
 426
 427     }//end of loop over theta
 428
 429     //calculate the average of fVtheta = fV
 430     fV /= fNtheta;
 431
 432     //sigma2 and chi (for statistical error calculations)
 433     Double_t  fSigma2 = 0;
 434     Double_t  fChi= 0;
 435     if (fEventNumber!=0) {
 436       fQsum /= fEventNumber;
 437       //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
 438       //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
 439       fQ2sum /= fEventNumber;
 440       //cout<<"fQ2sum = "<<fQ2sum<<endl;
 441       fSigma2 = fQ2sum - TMath::Power(fQsum.X(),2.) - TMath::Power(fQsum.Y(),2.) - TMath::Power(fV,2.);  //BP eq. 62
 442       if (fSigma2>0) fChi = fV/TMath::Sqrt(fSigma2);
 443       else fChi = -1.;
 444       fCommonHistsRes->FillChi(fChi);
 445       cout<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
 446     }
 447
 448     //copy content of profile into TH1D and add error
 449     for(Int_t b=0;b<fNbinsPt;b++) {
 450       Double_t fv2pro = fHistProVPt->GetBinContent(b);
 451       Double_t fNprime = fCommonHists->GetEntriesInPtBin(b);
 452       Double_t fErrdifcomb = 0.;  //set error to zero
 453       Double_t fErr2difcomb = 0.; //set error to zero
 454       //calculate error
 455       if (fNprime!=0.) {
 456         for (Int_t theta=0;theta<fNtheta;theta++) {
 457           Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi();
 458           Double_t fApluscomb = TMath::Exp((fJ01*fJ01)/(2*fChi*fChi)*
 459                                            TMath::Cos(fTheta));
 460           Double_t fAmincomb = TMath::Exp(-(fJ01*fJ01)/(2*fChi*fChi)*
 461                                           TMath::Cos(fTheta));
 462           fErr2difcomb += (TMath::Cos(fTheta)/(4*fNprime*TMath::BesselJ1(fJ01)*
 463                                                  TMath::BesselJ1(fJ01)))*
 464             ((fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2))) -
 465              (fAmincomb*TMath::BesselJ0(2*fJ01*TMath::Cos(fTheta/2))));
 466         } //loop over theta
 467       }
 468
 469       if (fErr2difcomb!=0.) {
 470         fErr2difcomb /= fNtheta;
 471         fErrdifcomb = TMath::Sqrt(fErr2difcomb)*100;
 472         //cerr<<"fErrdifcomb = "<<fErrdifcomb<<endl;
 473       }
 474       else {fErrdifcomb = 0.;}
 475
 476       //fill TH1D
 477       fCommonHistsRes->FillDifferentialFlow(b, fv2pro, fErrdifcomb);
 478     } //loop over bins b
 479
 480
 481     //save histograms in file
 482     fHistFile->cd();
 483     fHistFile->Write();
 484     fHistFile->Close();
 485     //Note that when the file is closed, all histograms and Trees in memory associated with this file are deleted
 486     if (fDebug) cout<<"****Histograms saved and fHistFile closed, all histograms deleted****"<<endl;
 487
 488     //close the first run file
 489     firstRunFile->Close();
 490
 491
 492   } //secondrun
 493
 494   cout<<"----LYZ analysis finished....----"<<endl<<endl;
 495
 496   return kTRUE;
 497 }
 498
 499
 500 //-----------------------------------------------------------------------
 501
 502  Bool_t AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent(AliFlowEventSimple* fEvent)
 503 {
 504   // Get event quantities from AliFlowEvent for all particles
 505
 506   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::FillFromFlowEvent()****"<<endl;
 507
 508   if (!fEvent){
 509     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
 510     return kFALSE;
 511   }
 512
 513   //define variables
 514   TComplex fExpo, fGtheta, fGthetaNew, fZ;
 515   Int_t fNtheta = AliFlowLYZConstants::kTheta;
 516   Int_t fNbins = AliFlowLYZConstants::kNbins;
 517
 518
 519   //calculate flow
 520   Double_t fOrder = 2.;
 521
 522   //get the Q vector
 523   fQ = fEvent->GetQ();
 524   //weight by the multiplicity
 525   Double_t fQX = 0;
 526   Double_t fQY = 0;
 527   if (fQ.GetMult() != 0) {
 528     fQX = fQ.X()/fQ.GetMult();
 529     fQY = fQ.Y()/fQ.GetMult();
 530   }
 531   fQ.Set(fQX,fQY);
 532   //for chi calculation:
 533   fQsum += fQ;
 534   fQ2sum += fQ.Mod2();
 535   //cerr<<"fQ2sum = "<<fQ2sum<<endl;
 536
 537   for (Int_t theta=0;theta<fNtheta;theta++)
 538     {
 539       fTheta = ((double)theta/fNtheta)*TMath::Pi()/fOrder;
 540
 541       //calculate fQtheta = cos(fOrder*(fPhi-fTheta);the projection of the Q vector on the reference direction fTheta
 542       fQtheta = GetQtheta(fQ, fTheta);
 543
 544       for (Int_t bin=1;bin<=fNbins;bin++)
 545         {
 546           Double_t fR = fHist1[theta]->GetBinCenter(bin); //bincentre of bins in histogram  //FIXED???
 547           //if (theta == 0) cerr<<"cerr::fR = "<<fR<<endl;
 548           if (fUseSum)
 549             {
 550               //calculate the sum generating function
 551               fExpo(0.,fR*fQtheta);                           //Re=0 ; Im=fR*fQtheta
 552               fGtheta = TComplex::Exp(fExpo);
 553             }
 554           else
 555             {
 556               //calculate the product generating function
 557               fGtheta = GetGrtheta(fEvent, fR, fTheta);  //make this function
 558               if (fGtheta.Rho2() > 100.) break;
 559             }
 560
 561           fHist1[theta]->Fill(fR,fGtheta);              //fill real and imaginary part of fGtheta
 562
 563         } //loop over bins
 564     } //loop over theta
 565
 566
 567   return kTRUE;
 568
 569
 570 }
 571
 572  //-----------------------------------------------------------------------
 573  Bool_t AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent(AliFlowEventSimple* fEvent)
 574 {
 575   //for differential flow
 576
 577   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::SecondFillFromFlowEvent()****"<<endl;
 578
 579   if (!fEvent){
 580     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
 581     return kFALSE;
 582   }
 583
 584   //define variables
 585   //TVector2 fQ;
 586   TComplex fExpo, fDenom, fNumer,fCosTermComplex;
 587   Double_t  fOrder, fPhi, fEta, fPt;
 588   Double_t fR0 = 0.;
 589   Double_t fCosTerm;
 590   Double_t m = 1.;
 591   Int_t fNtheta = AliFlowLYZConstants::kTheta;
 592
 593   //calculate flow
 594   fOrder = 2.;
 595
 596   //get the Q vector
 597   fQ = fEvent->GetQ();
 598   //weight by the multiplicity
 599   Double_t fQX = 0;
 600   Double_t fQY = 0;
 601   if (fQ.GetMult() != 0) {
 602     fQX = fQ.X()/fQ.GetMult();
 603     fQY = fQ.Y()/fQ.GetMult();
 604   }
 605   fQ.Set(fQX,fQY);
 606   //for chi calculation:
 607   fQsum += fQ;
 608   fQ2sum += fQ.Mod2();
 609
 610   for (Int_t theta=0;theta<fNtheta;theta++)
 611     {
 612       fTheta = ((double)theta/fNtheta)*TMath::Pi()/fOrder;
 613
 614       //calculate fQtheta = cos(fOrder*(fPhi-fTheta);the projection of the Q vector on the reference direction fTheta
 615       fQtheta = GetQtheta(fQ, fTheta);
 616       //cerr<<"fQtheta for fdenom = "<<fQtheta<<endl;
 617
 618       //denominator for differential v
 619       if (fHistProR0theta) {
 620         fR0 = fHistProR0theta->GetBinContent(theta+1);
 621       }
 622       else {
 623         cout <<"pointer fHistProR0Theta does not exist" << endl;
 624       }
 625       //cerr<<"fR0 = "<<fR0 <<endl;
 626
 627       if (fUseSum)                                                    //sum generating function
 628         {
 629           fExpo(0.,fR0*fQtheta);
 630           fDenom = fQtheta*(TComplex::Exp(fExpo)); //BP eq 12
 631           //loop over tracks in event
 632           Int_t fNumberOfTracks = fEvent->NumberOfTracks();
 633           for (Int_t i=0;i<fNumberOfTracks;i++)  {
 634             fTrack = fEvent->GetTrack(i);
 635             if (fTrack) {
 636               if (fTrack->UseForDifferentialFlow()) {
 637                 fEta = fTrack->Eta();
 638                 fPt = fTrack->Pt();
 639                 fPhi = fTrack->Phi();
 640                 fCosTerm = cos(m*fOrder*(fPhi-fTheta));
 641                 //cerr<<"fCosTerm = "<<fCosTerm <<endl;
 642                 fNumer = fCosTerm*(TComplex::Exp(fExpo));
 643                 if (fNumer.Rho()==0) {cerr<<"WARNING: modulus of fNumer is zero in SecondFillFromFlowEvent"<<endl;}
 644                 if (fDebug) cerr<<"modulus of fNumer is "<<fNumer.Rho()<<endl;
 645                 if (fHist2[theta]) {
 646                   fHist2[theta]->Fill(fEta,fPt,fNumer);
 647                 }
 648                 else {
 649                   cout << "fHist2 pointer mising" <<endl;
 650                 }
 651               }
 652             } //if particle
 653             else {cerr << "no particle!!!"<<endl;}
 654           } //loop over tracks
 655
 656         } //sum
 657       else                                                        //product generating function
 658         {
 659           fDenom = GetDiffFlow(fEvent, fR0, theta);
 660
 661         }//product
 662       if (fHistProReDenom && fHistProImDenom) {
 663         fHistProReDenom->Fill(theta,fDenom.Re());               //fill the real part of fDenom
 664         fHistProImDenom->Fill(theta,fDenom.Im());               //fill the imaginary part of fDenom
 665       }
 666       else {
 667         cout << "Pointers to fDenom  mising" << endl;
 668       }
 669
 670     }//end of loop over theta
 671
 672   return kTRUE;
 673
 674
 675 }
 676  //-----------------------------------------------------------------------
 677  Double_t AliFlowAnalysisWithLeeYangZeros::GetQtheta(TVector2 fQ, Double_t fTheta)
 678 {
 679   //calculate Qtheta. Qtheta is the sum over all particles of cos(fOrder*(fPhi-fTheta)) BP eq. 3
 680   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetQtheta()****"<<endl;
 681
 682   Double_t fQtheta = 0.;
 683   Double_t fOrder = 2.;
 684
 685   fQtheta = fQ.X()*cos(fOrder*fTheta)+fQ.Y()*sin(fOrder*fTheta);
 686
 687   return fQtheta;
 688
 689 }
 690
 691
 692 //-----------------------------------------------------------------------
 693 TComplex AliFlowAnalysisWithLeeYangZeros::GetGrtheta(AliFlowEventSimple* fEvent, Double_t fR, Double_t fTheta)
 694 {
 695   // Product Generating Function for LeeYangZeros method
 696   // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
 697
 698   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
 699
 700
 701   TComplex fG = TComplex::One();
 702   Double_t fOrder =  2.;
 703   Double_t fWgt = 1.;
 704
 705   Int_t fNumberOfTracks = fEvent->NumberOfTracks();
 706
 707   for (Int_t i=0;i<fNumberOfTracks;i++) //loop over tracks in event
 708     {
 709       fTrack = fEvent->GetTrack(i) ;
 710       if (fTrack){
 711         if (fTrack->UseForIntegratedFlow()) {
 712           Double_t fPhi = fTrack->Phi();
 713           Double_t fGIm = fR * fWgt*cos(fOrder*(fPhi - fTheta));
 714           TComplex fGi(1., fGIm);
 715           fG *= fGi;     //product over all tracks
 716         }
 717       }
 718       else {cerr << "no particle pointer !!!"<<endl;}
 719     }//loop over tracks
 720
 721   return fG;
 722
 723 }
 724
 725
 726 //-----------------------------------------------------------------------
 727 TComplex AliFlowAnalysisWithLeeYangZeros::GetDiffFlow(AliFlowEventSimple* fEvent, Double_t fR0, Int_t theta)
 728 {
 729   // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
 730   // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
 731   // Also for v1 mixed harmonics: DF Eq. 5
 732   // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
 733
 734   if (fDebug) cout<<"****AliFlowAnalysisWithLeeYangZeros::GetGrtheta()****"<<endl;
 735
 736   TComplex fG = TComplex::One();
 737   TComplex fdGr0(0.,0.);
 738   Double_t fOrder =  2.;
 739   Double_t fWgt = 1.;
 740
 741   Int_t fNumberOfTracks = fEvent->NumberOfTracks();
 742
 743   Int_t fNtheta = AliFlowLYZConstants::kTheta;
 744   Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi()/fOrder;
 745
 746   for (Int_t i=0;i<fNumberOfTracks;i++) //loop over tracks in event
 747     {
 748       fTrack = fEvent->GetTrack(i) ;
 749       if (fTrack){
 750         if (fTrack->UseForDifferentialFlow()) {
 751           Double_t fPhi = fTrack->Phi();
 752           Double_t fCosTerm = fWgt*cos(fOrder*(fPhi - fTheta));
 753           //GetGr0theta
 754           Double_t fGIm = fR0 * fCosTerm;
 755           TComplex fGi(1., fGIm);
 756           fG *= fGi;     //product over all tracks
 757           //GetdGr0theta
 758           TComplex fCosTermComplex(1., fR0*fCosTerm);
 759           fdGr0 +=(fCosTerm / fCosTermComplex);  //sum over all tracks
 760         }
 761       } //if particle
 762       else {cerr << "no particle!!!"<<endl;}
 763     }//loop over tracks
 764
 765
 766   for (Int_t i=0;i<fNumberOfTracks;i++)
 767     {
 768       fTrack = fEvent->GetTrack(i) ;
 769       if (fTrack){
 770         if (fTrack->UseForDifferentialFlow()) {
 771           Double_t fEta = fTrack->Eta();
 772           Double_t fPt = fTrack->Pt();
 773           Double_t fPhi = fTrack->Phi();
 774           Double_t fCosTerm = cos(fOrder*(fPhi-fTheta));
 775           TComplex fCosTermComplex(1.,fR0*fCosTerm);
 776           TComplex fNumer = fG*fCosTerm/fCosTermComplex;  //PG Eq. 9
 777           fHist2[theta]->Fill(fEta,fPt,fNumer);
 778         }
 779       } //if particle
 780       else {cerr << "no particle pointer!!!"<<endl;}
 781     }//loop over tracks
 782
 783   TComplex fDenom = fG*fdGr0;
 784   return fDenom;
 785
 786 }
 787
 788 //-----------------------------------------------------------------------
 789