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