PWG2/FLOW/AliFlowLeeYangZerosMaker.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 $Log$
  18 */
  19
  20 #include "Riostream.h"
  21 #include "AliFlowLeeYangZerosMaker.h"
  22 #include "AliFlowEvent.h"
  23 #include "AliFlowSelection.h"
  24 #include "AliFlowConstants.h"     //??
  25 #include "AliFlowLYZHist1.h"
  26 #include "AliFlowLYZHist2.h"
  27 #include "AliFlowLYZConstants.h"  //??
  28
  29 #include "TMath.h"
  30 #include "TComplex.h"
  31 #include "TProfile.h"
  32 #include "TObjArray.h"
  33 #include "TFile.h"
  34 #include "TVector.h"
  35 #include "TVector2.h"
  36 #include "TGraphErrors.h"
  37 #include "TCanvas.h"
  38
  39 class TTree;
  40 class TH1F;
  41 class TH1D;
  42 class TVector3;
  43 class TProfile2D;
  44 class TObject;
  45
  46 //class Riostream; //does not compile
  47 //class TMath;     //does not compile
  48 //class TVector;   //does not compile
  49
  50 //Description: Maker to analyze Flow using the LeeYangZeros method
  51 //             Equation numbers are from Big Paper (BP): Nucl. Phys. A 727, 373 (2003)
  52 //             Practical Guide (PG):    J. Phys. G: Nucl. Part. Phys. 30, S1213 (2004)
  53 //             Adapted from StFlowLeeYangZerosMaker.cxx
  54 //             by Markus Oldenberg and Art Poskanzer, LBNL
  55 //             with advice from Jean-Yves Ollitrault and Nicolas Borghini
  56 //
  57 //Author: Naomi van der Kolk (kolk@nikhef.nl)
  58
  59
  60
  61 ClassImp(AliFlowLeeYangZerosMaker)
  62
  63  //-----------------------------------------------------------------------
  64
  65   AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker():
  66     fFirstRun(kTRUE),
  67     fUseSum(kTRUE),
  68     fDebug(kFALSE),
  69     fHistFileName(0),
  70     fHistFile(0),
  71     fSummaryFile(0),
  72     firstRunFile(0)
  73
  74 {
  75   //default constructor
  76   if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker default constructor****"<<endl;
  77
  78   fFlowSelect = new AliFlowSelection();
  79   if (fDebug) { cerr<<"****fFlowSelect in constructor AliFlowLeeYangZerosMaker ("<<fFlowSelect<<")****"<<endl;}
  80   // output file (histograms)
  81   TString fHistFileName = "flowLYZAnalysPlot.root" ;
  82 }
  83
  84
  85 AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker(const AliFlowSelection* flowSelect):
  86   fFirstRun(kTRUE),
  87   fUseSum(kTRUE),
  88   fDebug(kFALSE),
  89   fHistFileName(0),
  90   fHistFile(0),
  91   fSummaryFile(0),
  92   firstRunFile(0)
  93 {
  94   //custum constructor
  95   if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker custum constructor****"<<endl;
  96
  97   if(flowSelect) { fFlowSelect = new AliFlowSelection(*flowSelect); }
  98   else    {
  99     fFlowSelect = new AliFlowSelection() ;
 100     if (fDebug) cerr<<"****fFlowSelect in constructor AliFlowLeeYangZerosMaker ("<<fFlowSelect<<")****"<<endl;
 101   }
 102   // output file (histograms)
 103   TString fHistFileName = "flowLYZAnalysPlot.root" ;
 104 }
 105
 106  //-----------------------------------------------------------------------
 107
 108
 109  AliFlowLeeYangZerosMaker::~AliFlowLeeYangZerosMaker()
 110  {
 111    //default destructor
 112    if (fDebug) cout<<"****~AliFlowLeeYangZerosMaker****"<<endl;
 113    delete fHistFile;
 114  }
 115
 116  //-----------------------------------------------------------------------
 117
 118 Bool_t AliFlowLeeYangZerosMaker::Init()
 119 {
 120   //init method
 121   if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::Init()****"<<endl;
 122
 123   // Open output files (->plots)
 124   fHistFile = new TFile(fHistFileName.Data(), "RECREATE");
 125   //fHistFile->cd() ;  //all histograms will be saved in this file
 126
 127   //for each harmonic ???
 128   fQsum.Set(0.,0.);
 129   fQ2sum = 0.;
 130
 131   // Book histograms
 132   Int_t fNtheta = AliFlowLYZConstants::kTheta;
 133
 134
 135   //for control histograms
 136   fHistOrigMult = new TH1F("Control_FlowLYZ_OrigMult", "Control_FlowLYZ_OrigMult",1000, 0., 10000.);
 137   fHistOrigMult->SetXTitle("Original Multiplicity");
 138   fHistOrigMult->SetYTitle("Counts");
 139
 140   fHistMult = new TH1F("Control_FlowLYZ_Mult", "Control_FlowLYZ_Mult",1000, 0., 10000.);
 141   fHistMult->SetXTitle("Multiplicity from selection");
 142   fHistMult->SetYTitle("Counts");
 143
 144   fHistQ = new TH1F("Control_FlowLYZ_Q","Control_FlowLYZ_Q",500, 0., 10.);
 145   fHistQ->SetXTitle("Qvector");
 146   fHistQ->SetYTitle("Counts");
 147
 148   fHistPt = new TH1F("Control_FlowLYZ_Pt","Control_FlowLYZ_Pt",200, 0., 10.);
 149   fHistPt->SetXTitle("Pt (GeV/c)");
 150   fHistPt->SetYTitle("Counts");
 151
 152   fHistPhi = new TH1F("Control_FlowLYZ_Phi","Control_FlowLYZ_Phi",70, 0., 7.);
 153   fHistPhi->SetXTitle("Phi");
 154   fHistPhi->SetYTitle("Counts");
 155
 156   fHistEta = new TH1F("Control_FlowLYZ_Eta","Control_FlowLYZ_Eta",40, 0., 2.);
 157   fHistEta->SetXTitle("Eta");
 158   fHistEta->SetYTitle("Counts");
 159
 160   fHistQtheta = new TH1F("Control_FlowLYZ_Qtheta","Control_FlowLYZ_Qtheta",50,-1000.,1000.);
 161   fHistQtheta->SetXTitle("Qtheta");
 162   fHistQtheta->SetYTitle("Counts");
 163
 164   if (fFirstRun){
 165     //for first loop over events
 166     fHistProR0thetaHar1  = new TProfile("First_FlowProLYZ_r0theta_Har1","First_FlowProLYZ_r0theta_Har1",fNtheta,-0.5,fNtheta-0.5);
 167     fHistProR0thetaHar1->SetXTitle("#theta");
 168     fHistProR0thetaHar1->SetYTitle("r_{0}^{#theta}");
 169
 170     fHistProR0thetaHar2  = new TProfile("First_FlowProLYZ_r0theta_Har2","First_FlowProLYZ_r0theta_Har2",fNtheta,-0.5,fNtheta-0.5);
 171     fHistProR0thetaHar2->SetXTitle("#theta");
 172     fHistProR0thetaHar2->SetYTitle("r_{0}^{#theta}");
 173
 174     fHistProVthetaHar1  = new TProfile("First_FlowProLYZ_Vtheta_Har1","First_FlowProLYZ_Vtheta_Har1",fNtheta,-0.5,fNtheta-0.5);
 175     fHistProVthetaHar1->SetXTitle("#theta");
 176     fHistProVthetaHar1->SetYTitle("V_{n}^{#theta}");
 177
 178     fHistProVthetaHar2  = new TProfile("First_FlowProLYZ_Vtheta_Har2","First_FlowProLYZ_Vtheta_Har2",fNtheta,-0.5,fNtheta-0.5);
 179     fHistProVthetaHar2->SetXTitle("#theta");
 180     fHistProVthetaHar2->SetYTitle("V_{n}^{#theta}");
 181
 182     fHistProVR0 = new TProfile("First_FlowProLYZ_vR0","First_FlowProLYZ_vR0",2,0.5,2.5,-100.,100.);
 183     fHistProVR0->SetXTitle("Harmonic");
 184     fHistProVR0->SetYTitle("v integrated from r0 (%)");
 185
 186     fHistVR0 = new TH1D("First_FlowLYZ_vR0","First_FlowLYZ_vR0",2,0.5,2.5);
 187     fHistVR0->SetXTitle("Harmonic");
 188     fHistVR0->SetYTitle("v integrated from r0 (%)");
 189
 190     fHistProV = new TProfile("First_FlowProLYZ_V","First_FlowProLYZ_V",2,0.5,2.5,-1000.,1000.);
 191     fHistProV->SetXTitle("Harmonic");
 192     fHistProV->SetYTitle("v integrated");
 193
 194     //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta, fHistProR0theta
 195     for (Int_t j=0;j<AliFlowConstants::kHars;j++)
 196       {
 197         for (Int_t theta=0;theta<fNtheta;theta++)
 198           {
 199             fHist1[j][theta]=new AliFlowLYZHist1(theta,j+1);
 200           }
 201       }
 202   }
 203   else {
 204     //for second loop over events
 205     fHistProReDenomHar1 = new TProfile("Second_FlowProLYZ_ReDenom_Har1","Second_FlowProLYZ_ReDenom_Har1" , fNtheta, -0.5, fNtheta-0.5);
 206     fHistProReDenomHar1->SetXTitle("#theta");
 207     fHistProReDenomHar1->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
 208
 209     fHistProReDenomHar2 = new TProfile("Second_FlowProLYZ_ReDenom_Har2","Second_FlowProLYZ_ReDenom_Har2" , fNtheta, -0.5, fNtheta-0.5);
 210     fHistProReDenomHar2->SetXTitle("#theta");
 211     fHistProReDenomHar2->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
 212
 213     fHistProImDenomHar1 = new TProfile("Second_FlowProLYZ_ImDenom_Har1","Second_FlowProLYZ_ImDenom_Har1" , fNtheta, -0.5, fNtheta-0.5);
 214     fHistProImDenomHar1->SetXTitle("#theta");
 215     fHistProImDenomHar1->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
 216
 217     fHistProImDenomHar2 = new TProfile("Second_FlowProLYZ_ImDenom_Har2","Second_FlowProLYZ_ImDenom_Har2" , fNtheta, -0.5, fNtheta-0.5);
 218     fHistProImDenomHar2->SetXTitle("#theta");
 219     fHistProImDenomHar2->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");
 220
 221     fHistProVetaHar1 = new TProfile("Second_FlowProLYZ_Veta_Har1","Second_FlowProLYZ_Veta_Har1",40,-2.,2.);
 222     fHistProVetaHar1->SetXTitle("rapidity");
 223     fHistProVetaHar1->SetYTitle("v (%)");
 224
 225     fHistProVetaHar2 = new TProfile("Second_FlowProLYZ_Veta_Har2","Second_FlowProLYZ_Veta_Har2",40,-2.,2.);
 226     fHistProVetaHar2->SetXTitle("rapidity");
 227     fHistProVetaHar2->SetYTitle("v (%)");
 228
 229     fHistProVPtHar1 = new TProfile("Second_FlowProLYZ_VPt_Har1","Second_FlowProLYZ_VPt_Har1",100,0.,10.);
 230     fHistProVPtHar1->SetXTitle("Pt");
 231     fHistProVPtHar1->SetYTitle("v (%)");
 232
 233     fHistProVPtHar2 = new TProfile("Second_FlowProLYZ_VPt_Har2","Second_FlowProLYZ_VPt_Har2",100,0.,10.);
 234     fHistProVPtHar2->SetXTitle("Pt");
 235     fHistProVPtHar2->SetYTitle("v (%)");
 236
 237     fHistVPtHar2 = new TH1D("Second_FlowLYZ_VPt_Har2","Second_FlowLYZ_VPt_Har2",100,0.,10.);
 238     fHistVPtHar2->SetXTitle("Pt");
 239     fHistVPtHar2->SetYTitle("v (%)");
 240
 241     fHistProReDtheta = new TProfile("Second_FlowProLYZ_ReDtheta_Har2","Second_FlowProLYZ_ReDtheta_Har2",fNtheta, -0.5, fNtheta-0.5);
 242     fHistProReDtheta->SetXTitle("#theta");
 243     fHistProReDtheta->SetYTitle("Re(D^{#theta})");
 244
 245     fHistProImDtheta = new TProfile("Second_FlowProLYZ_ImDtheta_Har2","Second_FlowProLYZ_ImDtheta_Har2",fNtheta, -0.5, fNtheta-0.5);
 246     fHistProImDtheta->SetXTitle("#theta");
 247     fHistProImDtheta->SetYTitle("Im(D^{#theta})");
 248
 249
 250     //class AliFlowLYZHist2 defines the histograms:
 251     for (Int_t j=0;j<AliFlowConstants::kHars;j++)
 252       {
 253         for (Int_t theta=0;theta<fNtheta;theta++)
 254           {
 255             fHist2[j][theta]=new AliFlowLYZHist2(theta,j+1);
 256           }
 257       }
 258
 259     //read hists from first run file
 260     //firstRunFile = new TFile("fof_flowLYZAnal_firstrun.root","READ");  //default is read
 261     if (firstRunFile->IsZombie()){ //check if file exists
 262       cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
 263       exit(-1);
 264     } else if (firstRunFile->IsOpen()){
 265       cout<<"----firstRunFile is open----"<<endl<<endl;
 266       fHistProVthetaHar1  = (TProfile*)firstRunFile->Get("First_FlowProLYZ_Vtheta_Har1");
 267       fHistProVthetaHar2  = (TProfile*)firstRunFile->Get("First_FlowProLYZ_Vtheta_Har2");
 268       fHistProR0thetaHar2  = (TProfile*)firstRunFile->Get("First_FlowProLYZ_r0theta_Har2");
 269       fHistProV = (TProfile*)firstRunFile->Get("First_FlowProLYZ_V");
 270     }
 271   }
 272
 273
 274   if (fDebug) cout<<"****Histograms initialised****"<<endl;
 275   if (fDebug) cout<<"****fFlowSelect in Init() "<<fFlowSelect<<"****"<<endl;
 276
 277   fEventNumber = 0; //set event counter to zero
 278   /*
 279   if (fUseSum)
 280     {
 281       //initialize LYZ summary class
 282       fLYZSummary = new AliFlowLYZSummary();
 283       fSummaryFile = new TFile("fFlowSummary.root","RECREATE","Flow LYZ summary file");
 284       fSummaryFile->SetCompressionLevel(1);
 285       fFlowTree = new TTree("FlowTree", "Flow Summary Tree");
 286       fFlowTree->SetAutoSave(1000000);  // autosave when 1 Mbyte written
 287       fFlowTree->Branch("fLYZSummary","AliFlowLYZSummary",&fLYZSummary,25000,99);
 288     }
 289   */
 290   return kTRUE;
 291 }
 292
 293  //-----------------------------------------------------------------------
 294
 295 Bool_t AliFlowLeeYangZerosMaker::Make(AliFlowEvent* fFlowEvent)
 296 {
 297   //make method
 298   if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::Make()****"<<endl;
 299
 300   //get tracks from event
 301   if (fFlowEvent) {
 302     fFlowTracks = fFlowEvent->TrackCollection();
 303     if (fDebug) cout<<"****fFlowSelect in Make() "<<fFlowSelect<<"****"<<endl;
 304     if (fDebug) fFlowSelect->PrintSelectionList() ;
 305     if (fDebug) fFlowSelect->PrintList() ;
 306
 307     if (fFlowSelect && fFlowSelect->Select(fFlowEvent))   // check if event is selected
 308       {
 309         fFlowTracks = fFlowEvent->TrackCollection();     //get tracks from event
 310         fFlowEvent->SetSelections(fFlowSelect) ;                // does the selection of tracks for r.p. calculation (sets flags in AliFlowTrack)
 311         if (fFirstRun){
 312           MakeControlHistograms(fFlowEvent);
 313           FillFromFlowEvent(fFlowEvent);
 314         }
 315         else {
 316           MakeControlHistograms(fFlowEvent);
 317           SecondFillFromFlowEvent(fFlowEvent);
 318         }
 319       }
 320   }
 321   else {
 322     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
 323     return kFALSE;
 324   }
 325
 326   fEventNumber++;
 327
 328   return kTRUE;
 329 }
 330
 331
 332   //-----------------------------------------------------------------------
 333  Bool_t AliFlowLeeYangZerosMaker::Finish()
 334 {
 335   //finish method
 336   if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::Finish()****"<<endl;
 337
 338  //define variables for both runs
 339   Float_t  fR0 = 0;
 340   Float_t  fv = 0;
 341   Float_t  fVtheta = 0;
 342   Float_t  fSigma2 = 0;
 343   Float_t  fChi= 0;
 344   Float_t  fJ01 = 2.405;
 345   Int_t fNtheta = AliFlowLYZConstants::kTheta;
 346   Float_t  fV = 0;
 347
 348   if (fFirstRun){
 349     for (Int_t j=0;j<AliFlowConstants::kHars;j++)
 350       //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
 351       {
 352         Float_t  fMeanMult = fHistMult->GetMean();
 353         //cerr<<"fMeanMult = "<<fMeanMult<<endl;
 354
 355         for (Int_t theta=0;theta<fNtheta;theta++)
 356           {
 357             //get the first minimum r0
 358             fHist1[j][theta]->FillGtheta();
 359             fR0 = fHist1[j][theta]->GetR0();
 360             //cerr<<"fR0 = "<<fR0<<endl;
 361
 362             //calculate integrated flow
 363             if (fR0!=0) fVtheta = fJ01/fR0;
 364             if (fMeanMult!=0.)
 365               {
 366                 fv = fVtheta/fMeanMult;
 367                 cerr<<"fv = "<<fv<<endl;
 368               }
 369
 370
 371             //fill the histograms
 372             if (j==0)
 373               {
 374                 fHistProR0thetaHar1->Fill(theta,fR0);
 375                 fHistProVthetaHar1->Fill(theta,fVtheta);
 376                 fHistProV->Fill(j+1,fVtheta);    //profile takes care of the averaging over theta.
 377               }
 378             if (j==1)
 379               {
 380                 fHistProR0thetaHar2->Fill(theta,fR0);
 381                 fHistProVthetaHar2->Fill(theta,fVtheta);
 382                 fHistProV->Fill(j+1,fVtheta);   //profile takes care of the averaging over theta.
 383               }
 384
 385             fHistProVR0->Fill(j+1,fv*100);    //*100 to get %     //profile takes care of the averaging over theta.
 386
 387           } //end of loop over theta
 388
 389         //sigma2 and chi
 390         if (j==1)   //second harmonic only temporarily
 391           {
 392             if (fEventNumber!=0) {
 393               fQsum /= fEventNumber;
 394               //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
 395               //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
 396               fQ2sum /= fEventNumber;
 397               //cerr<<"fQ2sum = "<<fQ2sum<<endl;
 398               fV = fHistProV->GetBinContent(j+1);
 399               fSigma2 = fQ2sum - TMath::Power(fQsum.X(),2.) - TMath::Power(fQsum.Y(),2.) - TMath::Power(fV,2.);  //BP eq. 62
 400               if (fSigma2>0) fChi = fV/TMath::Sqrt(fSigma2);
 401               else fChi = -1.;
 402               cerr<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
 403             }
 404
 405           } //j==1
 406
 407       } //end of loop over harmonics
 408
 409     // recalculate statistical errors on integrated flow
 410     //combining 5 theta angles to 1 relative error BP eq. 89
 411     Double_t fRelErr2comb = 0.;
 412     Int_t nEvts = fEventNumber;
 413     for (Int_t theta=0;theta<fNtheta;theta++){
 414       Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi();
 415       Double_t fApluscomb = TMath::Exp((fJ01*fJ01)/(2*fChi*fChi)*
 416                                        TMath::Cos(fTheta));
 417       Double_t fAmincomb = TMath::Exp(-(fJ01*fJ01)/(2*fChi*fChi)*
 418                                       TMath::Cos(fTheta));
 419       fRelErr2comb += (1/(2*nEvts*(fJ01*fJ01)*TMath::BesselJ1(fJ01)*
 420                           TMath::BesselJ1(fJ01)))*
 421         (fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2)) +
 422          fAmincomb*TMath::BesselJ0(2*fJ01)*TMath::Cos(fTheta/2));
 423     }
 424     fRelErr2comb /= fNtheta;
 425     Double_t fRelErrcomb = TMath::Sqrt(fRelErr2comb);
 426     cerr<<"fRelErrcomb = "<<fRelErrcomb<<endl;
 427
 428     //copy content of profile into TH1D and add error
 429     for(Int_t b=0;b<2;b++){
 430       Double_t fv2pro = fHistProVR0->GetBinContent(b+1);
 431       fHistVR0->SetBinContent(b+1,fv2pro);
 432       Double_t fv2Err = fv2pro*fRelErrcomb ;
 433       cerr<<"fv2pro +- fv2Err = "<<fv2pro<<" +- "<<fv2Err<<" for bin "<<b+1<<endl;
 434       fHistVR0->SetBinError(b+1,fv2Err);
 435     }
 436
 437
 438     if (fDebug) cout<<"****histograms filled****"<<endl;
 439     /*
 440       if (fUseSum)
 441       {
 442       if (fSummaryFile->IsOpen())
 443       {
 444       fSummaryFile->Write(0,TObject::kOverwrite);
 445       fSummaryFile->Close();
 446       }
 447       }
 448     */
 449
 450     //save histograms in file //temp for testing selector
 451     fHistFile->cd();
 452     fHistFile->Write();
 453
 454     return kTRUE;
 455     fEventNumber =0; //set to zero for second round over events
 456   }  //firstrun
 457
 458
 459   else {  //second run
 460
 461     //calculate differential flow
 462     //declare variables
 463     Float_t fEta, fPt, fReRatio, fVeta, fVPt;
 464     Float_t fReDenom = 0;
 465     Float_t fImDenom = 0;
 466     Double_t fR0 = 0;
 467     TComplex fDenom, fNumer, fDtheta;
 468     Int_t m = 1;
 469     TComplex i = TComplex::I();
 470     Float_t fBesselRatio[3] = {1., 1.202, 2.69};
 471     Double_t fErrdifcomb = 0.;  //set error to zero
 472     Double_t fErr2difcomb = 0.; //set error to zero
 473
 474     /*
 475     //v1 integrated
 476     for (Int_t j=0;j<AliFlowConstants::kHars;j++)   //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
 477     {
 478     for (Int_t theta=0;theta<fNtheta;theta++)
 479     {
 480     //get the first minimum r0
 481     //fR0 = GetR0(fHist1[j][theta]->FillGtheta());
 482
 483     fReDenom = fHistProReDenomHar2->GetBinContent(theta+1);
 484     fImDenom = fHistProImDenomHar2->GetBinContent(theta+1);
 485     TComplex fDenom(fReDenom,fImDenom);
 486
 487     //complete!!
 488
 489     }//end of loop over theta
 490     }//end of loop over harmonics
 491     */
 492
 493     //differential flow
 494     //for (Int_t j=0;j<AliFlowConstants::kHars;j++)    //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
 495     //{
 496     Int_t j=1; //temp only harm 2
 497     //fFlowSelect->SetHarmonic(j);  //not needed here?
 498
 499     for (Int_t theta=0;theta<fNtheta;theta++)  { //loop over theta
 500       if (j==0) {
 501         fR0 = fHistProR0thetaHar1->GetBinContent(theta+1);
 502         fVtheta = 2.405/fR0;
 503         //fVtheta =  fHistProVthetaHar1->GetBinContent(theta+1);  // BP Eq. 9  -> Vn^theta = j01/r0^theta
 504         fReDenom = fHistProReDenomHar1->GetBinContent(theta+1);
 505         fImDenom = fHistProImDenomHar1->GetBinContent(theta+1);
 506       }
 507       if (j==1) {
 508         fR0 = fHistProR0thetaHar2->GetBinContent(theta+1);
 509         if (fDebug) cerr<<"fR0 = "<<fR0<<endl;
 510         fVtheta = 2.405/fR0;                                    // BP Eq. 9  -> Vn^theta = j01/r0^theta
 511         fReDenom = fHistProReDenomHar2->GetBinContent(theta+1);
 512         fImDenom = fHistProImDenomHar2->GetBinContent(theta+1);
 513       }
 514
 515       fDenom(fReDenom,fImDenom);
 516
 517
 518       //for new method and use by others (only with the sum generating function):
 519       if (fUseSum) {
 520         fR0 = fHistProR0thetaHar2->GetBinContent(theta+1);
 521         fDtheta = fR0*fDenom;
 522         fHistProReDtheta->Fill(theta,fDtheta.Re());
 523         fHistProImDtheta->Fill(theta,fDtheta.Im());
 524       }
 525
 526       fDenom *= TComplex::Power(i, m-1);
 527       //cerr<<"TComplex::Power(i, m-1) = "<<TComplex::Power(i, m-1).Rho()<<endl; //checked ok
 528
 529       //v as a function of eta
 530       Int_t fEtaBins = AliFlowLYZConstants::kEtaBins;
 531       for (Int_t be=1;be<=fEtaBins;be++)  {
 532         fEta = fHist2[j][theta]->GetBinCenter(be);
 533         fNumer = fHist2[j][theta]->GetfNumer(be);
 534         if (fNumer.Rho()==0) {
 535           if (fDebug) cerr<<"WARNING: modulus of fNumer is zero in Finish()"<<endl;
 536           fReRatio = 0;
 537         }
 538         else if (fDenom.Rho()==0) {
 539           if (fDebug) cerr<<"WARNING: modulus of fDenom is zero"<<endl;
 540           fReRatio = 0;
 541         }
 542         else {
 543           //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
 544           fReRatio = (fNumer/fDenom).Re();
 545         }
 546
 547         fVeta = fBesselRatio[m-1]*fReRatio*fVtheta*100.; //BP eq. 12
 548         //if ( j==1 && theta==0) cerr<<"eta = "<<fEta<<" cerr::fReRatio for eta = "<<fReRatio<<" cerr::fVeta for eta = "<<fVeta<<endl;
 549
 550         if (j==0) fHistProVetaHar1->Fill(fEta,fVeta);
 551         if (j==1) fHistProVetaHar2->Fill(fEta,fVeta);
 552       } //loop over bins be
 553
 554         //v as a function of Pt
 555       Int_t fPtBins = AliFlowLYZConstants::kPtBins;
 556       for (Int_t bp=1;bp<=fPtBins;bp++)  {
 557         fPt = fHist2[j][theta]->GetBinCenterPt(bp);
 558         fNumer = fHist2[j][theta]->GetfNumerPt(bp);
 559         if (fNumer.Rho()==0) {
 560           if (fDebug) cerr<<"modulus of fNumer is zero"<<endl;
 561           fReRatio = 0;
 562         }
 563         else if (fDenom.Rho()==0) {
 564           if (fDebug) cerr<<"modulus of fDenom is zero"<<endl;
 565           fReRatio = 0;
 566         }
 567         else {
 568           //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
 569           fReRatio = (fNumer/fDenom).Re();
 570         }
 571
 572         fVPt = fBesselRatio[m-1]*fReRatio*fVtheta*100.; //BP eq. 12
 573         //cerr<<"fBesselRatio[m-1] = "<<fBesselRatio[m-1]<<endl;   //checked ok
 574         //if ( j==1 && theta==0) cerr<<"pt = "<<fPt<<" cerr::fReRatio for pt = "<<fReRatio<<" cerr::fVPt for pt = "<<fVeta<<endl;
 575
 576         if (j==0) fHistProVPtHar1->Fill(fPt,fVPt);
 577         if (j==1) fHistProVPtHar2->Fill(fPt,fVPt);
 578       } //loop over bins bp
 579
 580     }//end of loop over theta
 581
 582
 583     //sigma2 and chi (for statistical error calculations)
 584     if (j==1) {  //second harmonic only temporarily
 585
 586       if (fEventNumber!=0) {
 587         fQsum /= fEventNumber;
 588         //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
 589         //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
 590         fQ2sum /= fEventNumber;
 591         cerr<<"fQ2sum = "<<fQ2sum<<endl;
 592         fV = fHistProV->GetBinContent(j+1);
 593         fSigma2 = fQ2sum - TMath::Power(fQsum.X(),2.) - TMath::Power(fQsum.Y(),2.) - TMath::Power(fV,2.);  //BP eq. 62
 594         if (fSigma2>0) fChi = fV/TMath::Sqrt(fSigma2);
 595         else fChi = -1.;
 596         cerr<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
 597       }
 598
 599     } //j==1
 600
 601
 602     //copy content of profile into TH1D and add error
 603     for(Int_t b=0;b<100;b++){
 604       Double_t fv2pro = fHistProVPtHar2->GetBinContent(b);
 605       //calculate error
 606       for (Int_t theta=0;theta<fNtheta;theta++) {
 607         Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi();
 608         Int_t Nprime = fHist2[j][theta]->GetNprimePt(b);
 609         //cerr<<"Nprime = "<<Nprime<<endl;
 610         if (Nprime!=0.) {
 611           Double_t fApluscomb = TMath::Exp((fJ01*fJ01)/(2*fChi*fChi)*
 612                                            TMath::Cos(fTheta));
 613           Double_t fAmincomb = TMath::Exp(-(fJ01*fJ01)/(2*fChi*fChi)*
 614                                           TMath::Cos(fTheta));
 615           fErr2difcomb += (TMath::Cos(fTheta)/(4*Nprime*TMath::BesselJ1(fJ01)*
 616                                                  TMath::BesselJ1(fJ01)))*
 617             ((fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2))) -
 618              (fAmincomb*TMath::BesselJ0(2*fJ01*TMath::Cos(fTheta/2))));
 619         }
 620       } //loop over theta
 621
 622       if (fErr2difcomb!=0.) {
 623         fErr2difcomb /= fNtheta;
 624         fErrdifcomb = TMath::Sqrt(fErr2difcomb)*100;
 625         //cerr<<"fErrdifcomb = "<<fErrdifcomb<<endl;
 626       }
 627       else {fErrdifcomb = 0.;}
 628
 629       //fill TH1D
 630       if (j==1) {
 631         fHistVPtHar2->SetBinContent(b,fv2pro);
 632         fHistVPtHar2->SetBinError(b,fErrdifcomb);
 633       }
 634       //check that error is set
 635       //if (j==1) { cout<<"difference between calculated error and error in hostogram: "<< fErrdifcomb - fHistVPtHar2->GetBinError(b)<<endl; }
 636
 637     } //loop over bins b
 638
 639
 640     //} //end of loop over harmonics         //temporarily out
 641
 642     //save histograms in file
 643     fHistFile->cd();
 644     fHistFile->Write();
 645     fHistFile->Close();
 646     //Note that when the file is closed, all histograms and Trees in memory associated with this file are deleted
 647     if (fDebug) cout<<"****Histograms saved and fHistFile closed, all histograms deleted****"<<endl;
 648
 649     //close the first run file
 650     firstRunFile->Close();
 651
 652
 653   } //secondrun
 654
 655   delete fFlowSelect;
 656   cout<<"----LYZ analysis finished....----"<<endl<<endl;
 657
 658   return kTRUE;
 659 }
 660
 661
 662 //-----------------------------------------------------------------------
 663   Bool_t AliFlowLeeYangZerosMaker::MakeControlHistograms(AliFlowEvent* fFlowEvent)
 664 {
 665   //contol histograms of pt, eta, phi, Q, mult
 666   if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::MakeControlHistograms()****"<<endl;
 667
 668   //define variables
 669   TVector2 fQ;
 670   Float_t fPt, fPhi, fEta, fQmult;
 671
 672   if (!fFlowEvent){
 673     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
 674     return kFALSE;
 675   }
 676
 677   if (!fFlowSelect){
 678     cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
 679     return kFALSE;
 680   }
 681
 682   //set selection and harmonic
 683   fFlowSelect->SetSelection(1);
 684   fFlowSelect->SetHarmonic(1); //second harmonic
 685
 686   //cerr<<"selection in MakeControlHistograms()"<<endl;
 687   //fFlowSelect->PrintSelectionList() ;
 688   //fFlowSelect->PrintList() ;
 689
 690   fFlowTracks = fFlowEvent->TrackCollection();
 691   Int_t fNumberOfTracks = fFlowTracks->GetEntries();
 692   fHistOrigMult->Fill(fNumberOfTracks);
 693   //cerr<<"fNumberOfTracks = "<<fNumberOfTracks<<endl;
 694   Int_t fMult = fFlowEvent->Mult(fFlowSelect);  // Multiplicity of tracks in the specified Selection
 695   fHistMult->Fill(fMult);
 696   //cerr<<"Mult = "<<fMult<<endl;
 697
 698   fQ = fFlowEvent ->Q(fFlowSelect);
 699   fQmult = fQ.Mod()/TMath::Sqrt(fNumberOfTracks);
 700   fHistQ->Fill(fQmult);
 701
 702   Int_t tempmult = 0; //for testing
 703   for (Int_t i=0;i<fNumberOfTracks;i++)
 704     {
 705       fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;   //get object at array position i
 706       //if (fFlowSelect->SelectPart(fFlowTrack))           //if track is selected
 707       if (fFlowSelect->Select(fFlowTrack))
 708         {
 709           fPt = fFlowTrack->Pt();
 710           fHistPt->Fill(fPt);
 711           tempmult++;
 712           fPhi = fFlowTrack->Phi();
 713           if (fPhi<0.) fPhi+=2*TMath::Pi();
 714           fHistPhi->Fill(fPhi);
 715           fEta = fFlowTrack->Eta();
 716           fHistEta->Fill(fEta);
 717         }
 718     }
 719   if (fMult!=tempmult){cerr<<"ERROR: Mult() is not tempmult! "<<fMult<<" :: "<<tempmult<<endl<<endl;}
 720   //else {cerr<<"Mult()= tempmult "<<fMult<<" :: "<<tempmult<<endl<<endl;}
 721
 722   return kTRUE;
 723
 724
 725 }
 726
 727
 728
 729 //-----------------------------------------------------------------------
 730
 731  Bool_t AliFlowLeeYangZerosMaker::FillFromFlowEvent(AliFlowEvent* fFlowEvent)
 732 {
 733   // Get event quantities from AliFlowEvent for all particles
 734
 735   if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::FillFromFlowEvent()****"<<endl;
 736
 737   if (!fFlowEvent){
 738     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
 739     return kFALSE;
 740   }
 741
 742   if (!fFlowSelect){
 743     cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
 744     return kFALSE;
 745   }
 746
 747
 748   //define variables
 749   TComplex fExpo, fGtheta, fGthetaNew, fZ;
 750   //Int_t m;
 751   Int_t fNtheta = AliFlowLYZConstants::kTheta;
 752   Int_t fNbins = AliFlowLYZConstants::kNbins;
 753
 754
 755   //calculate flow
 756   fFlowSelect->SetSelection(1);
 757
 758   for (Int_t j=0;j<AliFlowConstants::kHars;j++)   //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
 759     {
 760       Int_t m=1;
 761       fFlowSelect->SetHarmonic(j);
 762       Float_t fOrder = (double)((j+1)/m);
 763       //cerr<<"fOrder = "<<fOrder<<endl;
 764
 765       //get the Q vector from the FlowEvent
 766       TVector2 fQ = fFlowEvent->Q(fFlowSelect);
 767       //for chi calculation:
 768       if (j==1)  //second harmonic only temporarily
 769         {
 770           fQsum += fQ;
 771           fQ2sum += fQ.Mod2();
 772           //cerr<<"fQ2sum = "<<fQ2sum<<endl;
 773         }
 774
 775       fFlowTracks = fFlowEvent->TrackCollection();
 776
 777       for (Int_t theta=0;theta<fNtheta;theta++)
 778         {
 779           fTheta = ((float)theta/fNtheta)*TMath::Pi()/fOrder;
 780           //cerr<<"fTheta = "<<fTheta<<endl;
 781
 782           //calculate fQtheta = cos(fOrder*(fPhi-fTheta);the projection of the Q vector on the reference direction fTheta
 783           fQtheta = GetQtheta(fFlowSelect,fFlowTracks,fTheta);
 784           //save fQtheta in AliFlowLYZSummary class
 785
 786           //something
 787           //AliFlowLYZSummary::SetQtheta(theta,fQtheta);
 788
 789           if (j==1 && theta==0) fHistQtheta->Fill(fQtheta);
 790           //cerr<<"fQtheta = "<<fQtheta<<endl;
 791
 792           for (Int_t bin=1;bin<=fNbins;bin++)
 793             {
 794               Float_t fR = fHist1[j][theta]->GetBinCenter(bin); //bincentre of bins in histogram
 795               //if (theta == 0) cerr<<"cerr::fR = "<<fR<<endl;
 796               if (fUseSum)
 797                 {
 798                   //calculate the sum generating function
 799                   fExpo(0.,fR*fQtheta);                           //Re=0 ; Im=fR*fQtheta
 800                   fGtheta = TComplex::Exp(fExpo);
 801                 }
 802               else
 803                 {
 804                   //calculate the product generating function
 805                   fGtheta = GetGrtheta(fFlowSelect,fR,fTheta);  //make this function
 806                   if (fGtheta.Rho2() > 100.) break;
 807                 }
 808
 809               fHist1[j][theta]->Fill(fR,fGtheta);              //fill real and imaginary part of fGtheta
 810
 811             } //loop over bins
 812         } //loop over theta
 813     } //loop over harmonics
 814
 815   return kTRUE;
 816
 817
 818 }
 819
 820  //-----------------------------------------------------------------------
 821  Bool_t AliFlowLeeYangZerosMaker::SecondFillFromFlowEvent(AliFlowEvent* fFlowEvent)
 822 {
 823   //for differential flow
 824
 825   if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::SecondFillFromFlowEvent()****"<<endl;
 826
 827   if (!fFlowEvent){
 828     cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
 829     return kFALSE;
 830   }
 831
 832   if (!fFlowSelect){
 833     cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
 834     return kFALSE;
 835   }
 836
 837   //define variables
 838   //TVector2 fQ;
 839   TComplex fExpo, fDenom, fNumer,fCosTermComplex;
 840   Float_t  fOrder, fR0, fPhi, fEta, fPt;
 841   Double_t fCosTerm;
 842   Double_t m = 1.;
 843   Int_t fNtheta = AliFlowLYZConstants::kTheta;
 844
 845   //calculate flow
 846   fFlowSelect->SetSelection(1);
 847
 848   //cerr<<"selection in SecondFillFromFlowEvent()"<<endl;
 849   //fFlowSelect->PrintSelectionList() ;
 850   //fFlowSelect->PrintList() ;
 851
 852
 853   for (Int_t j=0;j<AliFlowConstants::kHars;j++)   //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
 854     {
 855       m=1;
 856       fFlowSelect->SetHarmonic(j);
 857       fOrder = (double)((j+1)/m);
 858
 859       //get the Q vector from the FlowEvent
 860       TVector2 fQ = fFlowEvent->Q(fFlowSelect);
 861       //for chi calculation:
 862       if (j==1)  //second harmonic only temporarily
 863         {
 864           fQsum += fQ;
 865           fQ2sum += fQ.Mod2();
 866           //cerr<<"fQ2sum = "<<fQ2sum<<endl;
 867         }
 868
 869       fFlowTracks = fFlowEvent->TrackCollection();
 870
 871       for (Int_t theta=0;theta<fNtheta;theta++)
 872         {
 873           fTheta = ((float)theta/fNtheta)*TMath::Pi()/fOrder;
 874
 875           //calculate fQtheta = cos(fOrder*(fPhi-fTheta);the projection of the Q vector on the reference direction fTheta
 876           fQtheta = GetQtheta(fFlowSelect,fFlowTracks,fTheta);
 877           //cerr<<"fQtheta for fdenom = "<<fQtheta<<endl;
 878
 879           /*
 880             if (j==0)  //first harmonic
 881             {
 882             //denominator for differential v
 883             fR0 = fHistProR0thetaHar1->GetBinContent(theta+1);
 884             fExpo(0.,fR0*fQtheta);
 885             fDenom = fQtheta*(TComplex::Exp(fExpo)); //BP eq 12
 886             //cerr<<"fDenom.Re() = "<<fDenom.Re()<<endl;
 887             //cerr<<"fDenom.Im() = "<<fDenom.Im()<<endl;
 888
 889             //denominator for differential v
 890             // ****put in product generating function!!
 891
 892
 893             fHistProReDenomHar1->Fill(theta,fDenom.Re());               //fill the real part of fDenom
 894             fHistProImDenomHar1->Fill(theta,fDenom.Im());               //fill the imaginary part of fDenom
 895             }
 896           */
 897
 898           if (j==1)  //second harmonic
 899             {
 900               //denominator for differential v
 901               fR0 = fHistProR0thetaHar2->GetBinContent(theta+1);
 902               //cerr<<"fR0 = "<<fR0 <<endl;
 903
 904               if (fUseSum)                                                    //sum generating function
 905                 {
 906                   fExpo(0.,fR0*fQtheta);
 907                   fDenom = fQtheta*(TComplex::Exp(fExpo)); //BP eq 12
 908                   //loop over tracks in event
 909                   fFlowTracks = fFlowEvent->TrackCollection();
 910                   Int_t fNumberOfTracks = fFlowTracks->GetEntries();
 911                   for (Int_t i=0;i<fNumberOfTracks;i++)
 912                     {
 913                       fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;                //get object at array position i
 914                       if (fFlowSelect->SelectPart(fFlowTrack))                        //if track is selected
 915                         {
 916                           fPhi = fFlowTrack->Phi();
 917                           fCosTerm = cos(m*fOrder*(fPhi-fTheta));
 918                           //cerr<<"fCosTerm = "<<fCosTerm <<endl;
 919                           fNumer = fCosTerm*(TComplex::Exp(fExpo));
 920                           if (fNumer.Rho()==0) {cerr<<"WARNING: modulus of fNumer is zero in SecondFillFromFlowEvent"<<endl;}
 921                           if (fDebug) cerr<<"modulus of fNumer is "<<fNumer.Rho()<<endl;
 922                           fEta = fFlowTrack->Eta();  //rapidity
 923                           fPt = fFlowTrack->Pt();
 924                           fHist2[j][theta]->Fill(fEta,fPt,fNumer);
 925                         } //if
 926                     } //loop over tracks
 927
 928                 } //sum
 929               else                                                        //product generating function
 930                 {
 931                   fDenom = GetDiffFlow(fFlowSelect,fR0,theta);
 932
 933                 }//product
 934
 935               fHistProReDenomHar2->Fill(theta,fDenom.Re());               //fill the real part of fDenom
 936               fHistProImDenomHar2->Fill(theta,fDenom.Im());               //fill the imaginary part of fDenom
 937             }//j==1
 938
 939
 940         }//end of loop over theta
 941
 942     }//loop over harmonics
 943
 944
 945   return kTRUE;
 946
 947
 948
 949 }
 950  //-----------------------------------------------------------------------
 951  Double_t AliFlowLeeYangZerosMaker::GetQtheta(AliFlowSelection* fFlowSelect, TObjArray* fFlowTracks, Float_t fTheta)
 952 {
 953   //calculate Qtheta. Qtheta is the sum over all particles of cos(fOrder*(fPhi-fTheta)) BP eq. 3
 954   if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::GetQtheta()****"<<endl;
 955
 956   if (!fFlowSelect){
 957      cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
 958      return kFALSE;
 959    }
 960
 961
 962   Double_t fQtheta = 0.;
 963   Int_t fHarmonic = fFlowSelect->Har();
 964   Double_t fOrder = (double)(fHarmonic+1);
 965
 966   Int_t fNumberOfTracks = fFlowTracks->GetEntries();
 967   //cerr<<"GetQtheta::fNumberOfTracks = "<<fNumberOfTracks<<endl;
 968
 969   for (Int_t i=0;i<fNumberOfTracks;i++)                  //loop over tracks in event
 970     {
 971       fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;   //get object at array position i
 972       if(fFlowSelect->Select(fFlowTrack))                //if track is selected  //gives the same number of particles as Mult(fFlowSelect) method
 973         {
 974           Float_t fPhi = fFlowTrack->Phi();
 975           fQtheta += cos(fOrder*(fPhi-fTheta));
 976         }
 977
 978     }//loop over tracks
 979
 980   return fQtheta;
 981
 982 }
 983
 984 //-----------------------------------------------------------------------
 985 TComplex AliFlowLeeYangZerosMaker::GetGrtheta(AliFlowSelection* fFlowSelect, Float_t fR, Float_t fTheta)
 986 {
 987   // Product Generating Function for LeeYangZeros method
 988   // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
 989
 990   if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::GetGrtheta()****"<<endl;
 991
 992   if (!fFlowSelect){
 993     cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
 994     return kFALSE;
 995   }
 996
 997   TComplex fG = TComplex::One();
 998   Int_t fHarmonic = fFlowSelect->Har();
 999   Double_t fOrder = (double)(fHarmonic+1);
1000   Double_t fWgt = 1.;
1001
1002   Int_t fNumberOfTracks = fFlowTracks->GetEntries();
1003   //cerr<<"GetGrtheta::fNumberOfTracks = "<<fNumberOfTracks<<endl;
1004
1005   for (Int_t i=0;i<fNumberOfTracks;i++) //loop over tracks in event
1006     {
1007       fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;                //get object at array position i
1008       if (fFlowSelect->SelectPart(fFlowTrack))                        //if track is selected
1009         {
1010           Float_t fPhi = fFlowTrack->Phi();
1011           Double_t fGIm = fR * fWgt*cos(fOrder*(fPhi - fTheta));
1012           TComplex fGi(1., fGIm);
1013           fG *= fGi;     //product over all tracks
1014         }//if
1015     }//loop over tracks
1016
1017   return fG;
1018
1019   }
1020
1021
1022 //-----------------------------------------------------------------------
1023 TComplex AliFlowLeeYangZerosMaker::GetDiffFlow(AliFlowSelection* fFlowSelect, Float_t fR0, Int_t theta)
1024 {
1025   // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
1026   // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
1027   // Also for v1 mixed harmonics: DF Eq. 5
1028   // It is the deriverative of Grtheta at r0 divided by Grtheta at r0
1029
1030   if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::GetGrtheta()****"<<endl;
1031
1032   if (!fFlowSelect){
1033     cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
1034     return kFALSE;
1035   }
1036
1037
1038   TComplex fG = TComplex::One();
1039   TComplex fdGr0(0.,0.);
1040   Int_t fHarmonic = fFlowSelect->Har();
1041   Double_t fOrder = (double)(fHarmonic+1);
1042   Double_t fWgt = 1.;
1043
1044   Int_t fNumberOfTracks = fFlowTracks->GetEntries();
1045   //cerr<<"GetGrtheta::fNumberOfTracks = "<<fNumberOfTracks<<endl;
1046
1047   Int_t fNtheta = AliFlowLYZConstants::kTheta;
1048   Float_t fTheta = ((float)theta/fNtheta)*TMath::Pi()/fOrder;
1049
1050   for (Int_t i=0;i<fNumberOfTracks;i++) //loop over tracks in event
1051     {
1052       fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;                //get object at array position i
1053       if (fFlowSelect->SelectPart(fFlowTrack))                        //if track is selected
1054         {
1055           Float_t fPhi = fFlowTrack->Phi();
1056           Double_t fCosTerm = fWgt*cos(fOrder*(fPhi - fTheta));
1057           //GetGr0theta
1058           Double_t fGIm = fR0 * fCosTerm;
1059           TComplex fGi(1., fGIm);
1060           fG *= fGi;     //product over all tracks
1061           //GetdGr0theta
1062           TComplex fCosTermComplex(1., fR0*fCosTerm);
1063           fdGr0 +=(fCosTerm / fCosTermComplex);  //sum over all tracks
1064         }//if
1065     }//loop over tracks
1066
1067   for (Int_t i=0;i<fNumberOfTracks;i++)
1068     {
1069       fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;                //get object at array position i
1070       if (fFlowSelect->SelectPart(fFlowTrack))                        //if track is selected
1071         {
1072           Float_t fPhi = fFlowTrack->Phi();
1073           Double_t fCosTerm = cos(fOrder*(fPhi-fTheta));
1074           TComplex fCosTermComplex(1.,fR0*fCosTerm);
1075
1076           TComplex fNumer = fG*fCosTerm/fCosTermComplex;  //PG Eq. 9
1077           Float_t fEta = fFlowTrack->Eta();  //rapidity
1078           Float_t fPt = fFlowTrack->Pt();
1079           fHist2[1][theta]->Fill(fEta,fPt,fNumer);
1080         }//if
1081     }//loop over tracks
1082
1083   TComplex fDenom = fG*fdGr0;
1084   return fDenom;
1085
1086   }
1087
1088 //-------------------------------------------------