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