PWG2/FLOW/AliFlowAnalysisWithLYZEventPlane.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, 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 //#define AliFlowAnalysisWithLYZEventPlane_cxx
  21
  22 #include "Riostream.h"  //needed as include
  23
  24 #include "TFile.h"
  25 #include "TComplex.h"   //needed as include
  26 #include "TCanvas.h"   //needed as include
  27 #include "TLegend.h"   //needed as include
  28 #include "TProfile.h"  //needed as include
  29
  30 class TH1F;
  31
  32 #include "AliFlowLYZConstants.h"    //needed as include
  33 #include "AliFlowCommonConstants.h" //needed as include
  34 #include "AliFlowEventSimple.h"
  35 #include "AliFlowTrackSimple.h"
  36 #include "AliFlowCommonHist.h"
  37 #include "AliFlowCommonHistResults.h"
  38 #include "AliFlowLYZEventPlane.h"
  39 #include "AliFlowAnalysisWithLYZEventPlane.h"
  40
  41 class AliFlowVector;
  42
  43 // AliFlowAnalysisWithLYZEventPlane:
  44 //
  45 // Class to do flow analysis with the event plane from the LYZ method
  46 //
  47 // author: N. van der Kolk (kolk@nikhef.nl)
  48
  49
  50 ClassImp(AliFlowAnalysisWithLYZEventPlane)
  51
  52   //-----------------------------------------------------------------------
  53
  54  AliFlowAnalysisWithLYZEventPlane::AliFlowAnalysisWithLYZEventPlane():
  55   fOutFile(0),
  56   fFirstRunFile(0),
  57   fSecondRunFile(0),
  58   fFirstRunFileName(0),
  59   fSecondRunFileName(0),
  60   fOutFileName(0),
  61   fSecondReDtheta(0),
  62   fSecondImDtheta(0),
  63   fFirstr0theta(0),
  64   fSecondVPt(0),
  65   fHistProFlow(0),
  66   fHistProFlow2(0),
  67   fHistProWr(0),
  68   fHistProWrCorr(0),
  69   fHistFlow(0),
  70   fHistDeltaPhi(0),
  71   fHistDeltaPhi2(0),
  72   fHistDeltaPhihere(0),
  73   fHistPhiEP(0),
  74   fHistPhiEPhere(0),
  75   fHistPhiLYZ(0),
  76   fHistPhiLYZ2(0),
  77   fHistProR0theta(0),
  78   fHistProReDtheta(0),
  79   fHistProImDtheta(0),
  80   fCommonHists(0),
  81   fCommonHistsRes(0),
  82   fEventNumber(0),
  83   fQ(NULL),
  84   fQsum(NULL),
  85   fQ2sum(0),
  86   fQtheta(0),
  87   fTrack(NULL)
  88 //  fLYZEP(NULL)
  89 {
  90
  91   // Constructor.
  92   //  fQ.Set(0.,0.);           // flow vector
  93   //  fQsum.Set(0.,0.);        // flow vector sum
  94   fQ = new AliFlowVector();           // flow vector
  95   fQsum = new TVector2();        // flow vector sum
  96   //  fLYZEP = new AliFlowLYZEventPlane();
  97 }
  98
  99
 100
 101  //-----------------------------------------------------------------------
 102
 103
 104  AliFlowAnalysisWithLYZEventPlane::~AliFlowAnalysisWithLYZEventPlane()
 105  {
 106    //destructor
 107    delete fQ;
 108    delete fQsum;
 109    //   delete fLYZEP;
 110  }
 111
 112
 113 //-----------------------------------------------------------------------
 114 void AliFlowAnalysisWithLYZEventPlane::Init() {
 115
 116   //Initialise all histograms
 117   cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
 118
 119   //input histograms
 120   if (fSecondRunFile->IsZombie()){ //check if file exists
 121     cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
 122     exit(-1);
 123   } else if (fSecondRunFile->IsOpen()){
 124     cout<<"----secondRunFile is open----"<<endl;
 125     fSecondVPt = (TProfile*)fSecondRunFile->Get("Flow_Differential_Pt_LYZ"); //to compare to
 126   }
 127
 128   if (fFirstRunFile->IsZombie()){ //check if file exists
 129     cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
 130     exit(-1);
 131   } else if (fFirstRunFile->IsOpen()){
 132     cout<<"----firstRunFile is open----"<<endl<<endl;
 133     fFirstr0theta = (TProfile*)fFirstRunFile->Get("First_FlowPro_r0theta_LYZ");
 134   }
 135
 136
 137   //output file
 138   // ********make output file
 139   // analysis file (output)
 140   fOutFile = new TFile(fOutFileName.Data(),"RECREATE") ;
 141
 142   fCommonHists = new AliFlowCommonHist("LYZEP");
 143   fCommonHistsRes = new AliFlowCommonHistResults("LYZEP");
 144
 145   // output histograms
 146   Int_t fNtheta = AliFlowLYZConstants::kTheta;
 147   Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 148   Double_t  fPtMin = AliFlowCommonConstants::GetPtMin();
 149   Double_t  fPtMax = AliFlowCommonConstants::GetPtMax();
 150
 151
 152   fHistProFlow = new TProfile("FlowPro_VPt_LYZEP","FlowPro_VPt_LYZEP",fNbinsPt,fPtMin,fPtMax);
 153   fHistProFlow->SetXTitle("Pt");
 154   fHistProFlow->SetYTitle("v2 (%)");
 155
 156   fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
 157   fHistProWr->SetXTitle("Q");
 158   fHistProWr->SetYTitle("Wr");
 159
 160   fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
 161   fHistDeltaPhi->SetXTitle("DeltaPhi");
 162   fHistDeltaPhi->SetYTitle("Counts");
 163
 164   fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
 165   fHistPhiLYZ->SetXTitle("Phi from LYZ");
 166   fHistPhiLYZ->SetYTitle("Counts");
 167
 168   fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
 169   fHistPhiEP->SetXTitle("Phi from EP");
 170   fHistPhiEP->SetYTitle("Counts");
 171
 172
 173   fHistProR0theta  = new TProfile("FlowPro_r0theta_LYZEP","FlowPro_r0theta_LYZEP",fNtheta,-0.5,fNtheta-0.5);
 174   fHistProR0theta->SetXTitle("#theta");
 175   fHistProR0theta->SetYTitle("r_{0}^{#theta}");
 176
 177   fHistProReDtheta = new TProfile("FlowPro_ReDtheta_LYZEP","FlowPro_ReDtheta_LYZEP",fNtheta, -0.5, fNtheta-0.5);
 178   fHistProReDtheta->SetXTitle("#theta");
 179   fHistProReDtheta->SetYTitle("Re(D^{#theta})");
 180
 181   fHistProImDtheta = new TProfile("FlowPro_ImDtheta_LYZEP","FlowPro_ImDtheta_LYZEP",fNtheta, -0.5, fNtheta-0.5);
 182   fHistProImDtheta->SetXTitle("#theta");
 183   fHistProImDtheta->SetYTitle("Im(D^{#theta})");
 184
 185   fEventNumber = 0;  //set number of events to zero
 186
 187
 188 }
 189
 190 //-----------------------------------------------------------------------
 191
 192 void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* fLYZEP) {
 193 //void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent) {
 194
 195   //Get the event plane and weight for each event
 196   if (anEvent) {
 197
 198     //fill control histograms
 199     fCommonHists->FillControlHistograms(anEvent);
 200
 201     //get the Q vector from the FlowEvent
 202     *fQ = anEvent->GetQ();
 203     //Weight with the multiplicity
 204     Double_t fQX = 0.;
 205     Double_t fQY = 0.;
 206     if (fQ->GetMult()!=0.) {
 207       fQX = fQ->X()/fQ->GetMult();
 208       fQY = fQ->Y()/fQ->GetMult();
 209     } else {cerr<<"fQ->GetMult() is zero!"<<endl; }
 210     fQ->Set(fQX,fQY);
 211     //cout<<"fQ->Mod() = " << fQ->Mod() << endl;
 212     //for chi calculation:
 213     *fQsum += *fQ;
 214     //cout<<"fQsum.Mod() = "<<fQsum.Mod()<<endl;
 215     fQ2sum += fQ->Mod2();
 216     cout<<"fQ2sum = "<<fQ2sum<<endl;
 217
 218     //call AliFlowLYZEventPlane::CalculateRPandW() here!
 219     fLYZEP->CalculateRPandW(*fQ);
 220
 221     Double_t fWR = fLYZEP->GetWR();
 222     Double_t fRP = fLYZEP->GetPsi();
 223
 224     //fHistProWr->Fill(fQ.Mod(),fWR); //this weight is always positive
 225     fHistPhiLYZ->Fill(fRP);
 226
 227     //plot difference between event plane from EP-method and LYZ-method
 228     Double_t fRPEP = fQ->Phi()/2;                              //gives distribution from (0 to pi)
 229     //Float_t fRPEP = 0.5*TMath::ATan2(fQ.Y(),fQ.X());       //gives distribution from (-pi/2 to pi/2)
 230     //cout<<"fRPEP = "<< fRPEP <<endl;
 231     fHistPhiEP->Fill(fRPEP);
 232
 233     Double_t fDeltaPhi = fRPEP - fRP;
 234     if (fDeltaPhi < 0.) { fDeltaPhi += TMath::Pi(); }        //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
 235     //cout<<"fDeltaPhi = "<<fDeltaPhi<<endl;
 236     fHistDeltaPhi->Fill(fDeltaPhi);
 237
 238     //Flip sign of WR
 239     Double_t low = TMath::Pi()/4.;
 240     Double_t high = 3.*(TMath::Pi()/4.);
 241     if ((fDeltaPhi > low) && (fDeltaPhi < high)){
 242       fRP -= (TMath::Pi()/2);
 243       fWR = -fWR;
 244       cerr<<"*** fRP modified ***"<<endl;
 245     }
 246     fHistProWr->Fill(fQ->Mod(),fWR); //corrected weight
 247
 248     //calculate flow
 249     //loop over the tracks of the event
 250     Int_t fNumberOfTracks = anEvent->NumberOfTracks();
 251     for (Int_t i=0;i<fNumberOfTracks;i++)
 252       {
 253         fTrack = anEvent->GetTrack(i) ;
 254         if (fTrack){
 255           if (fTrack->UseForDifferentialFlow()) {
 256             Double_t fPhi = fTrack->Phi();
 257             //if (fPhi<0.) fPhi+=2*TMath::Pi();
 258             //calculate flow v2:
 259             Double_t fv2 = fWR * TMath::Cos(2*(fPhi-fRP));
 260             Double_t fPt = fTrack->Pt();
 261             //fill histogram
 262             fHistProFlow->Fill(fPt,100*fv2);
 263           }
 264         }//track selected
 265       }//loop over tracks
 266
 267     fEventNumber++;
 268     cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 269   }
 270 }
 271
 272   //--------------------------------------------------------------------
 273 void AliFlowAnalysisWithLYZEventPlane::Finish() {
 274
 275   //plot histograms etc.
 276   cout<<"AliFlowAnalysisWithLYZEventPlane::Terminate()"<<endl;
 277   //constands:
 278   Double_t  fJ01 = 2.405;
 279   Int_t fNtheta = AliFlowLYZConstants::kTheta;
 280   Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 281   //Double_t fErr2difcomb = 0.;
 282   //Double_t fErrdifcomb = 0.;
 283
 284   //calculate fV the mean of fVtheta
 285   Double_t  fVtheta = 0;
 286   Double_t  fV = 0;
 287   for (Int_t theta=0;theta<fNtheta;theta++)     {
 288     Double_t fR0 = fFirstr0theta->GetBinContent(theta+1);
 289     if (fR0!=0.) { fVtheta = fJ01/fR0 ;}
 290     fV += fVtheta;
 291   }
 292   fV /= fNtheta;
 293
 294   //calculate fChi
 295   Double_t  fSigma2 = 0;
 296   Double_t  fChi= 0;
 297   if (fEventNumber!=0) {
 298     *fQsum /= fEventNumber;
 299     //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
 300     //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
 301     fQ2sum /= fEventNumber;
 302     cerr<<"fEventNumber = "<<fEventNumber<<endl;
 303     cerr<<"fQ2sum = "<<fQ2sum<<endl;
 304     fSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(fV,2.);  //BP eq. 62
 305     if (fSigma2>0) fChi = fV/TMath::Sqrt(fSigma2);
 306     else fChi = -1.;
 307     fCommonHistsRes->FillChi(fChi);
 308     cerr<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
 309   }
 310
 311   for(Int_t b=0;b<fNbinsPt;b++){
 312     Double_t fv2pro = 0.;
 313     Double_t fErr2difcomb = 0.;
 314     Double_t fErrdifcomb = 0.;
 315     if(fHistProFlow) {
 316       fv2pro = fHistProFlow->GetBinContent(b);
 317       //calculate error
 318       for (Int_t theta=0;theta<fNtheta;theta++) {
 319         Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi();
 320         Int_t fNprime = TMath::Nint(fHistProFlow->GetBinEntries(b));
 321         //cerr<<"fNprime = "<<fNprime<<endl;
 322         if (fNprime!=0.) {
 323           Double_t fApluscomb = TMath::Exp((fJ01*fJ01)/(2*fChi*fChi)*
 324                                            TMath::Cos(fTheta));
 325           Double_t fAmincomb = TMath::Exp(-(fJ01*fJ01)/(2*fChi*fChi)*
 326                                           TMath::Cos(fTheta));
 327           fErr2difcomb += (TMath::Cos(fTheta)/(4*fNprime*TMath::BesselJ1(fJ01)*
 328                                                  TMath::BesselJ1(fJ01)))*
 329             ((fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2))) -
 330              (fAmincomb*TMath::BesselJ0(2*fJ01*TMath::Cos(fTheta/2))));
 331         } //if !=0
 332         //else { cout<<"fNprime = 0."<<endl; }
 333       } //loop over theta
 334
 335       if (fErr2difcomb!=0.) {
 336         fErr2difcomb /= fNtheta;
 337         fErrdifcomb = TMath::Sqrt(fErr2difcomb)*100;
 338         //cerr<<"fErrdifcomb = "<<fErrdifcomb<<endl;
 339       }
 340       else {fErrdifcomb = 0.; }
 341
 342       //fill TH1D
 343       fCommonHistsRes->FillDifferentialFlow(b, fv2pro, fErrdifcomb);
 344
 345     } //if fHistProFLow
 346     else  {
 347       cout << "Profile Hist missing" << endl;
 348       break;
 349     }
 350
 351   } //loop over b
 352
 353   // write to file
 354   fOutFile->Write();
 355
 356   cout<<"Making some plots to check the results:"<<endl<<endl;
 357
 358   TCanvas *canvas = new TCanvas("canvas","compare v2 vs pt",800,800);
 359   canvas->cd();
 360   fSecondVPt->SetLineColor(3);
 361   fSecondVPt->SetLineWidth(2);
 362   fSecondVPt->Draw();
 363   fHistFlow = fCommonHistsRes->GetfHistDiffFlow();
 364   fHistFlow->Draw("SAME");
 365   // draw the legend
 366   TLegend *legend2 = new TLegend(0.6,0.65,0.88,0.85);
 367   legend2->SetTextFont(72);
 368   legend2->SetTextSize(0.04);
 369   legend2->AddEntry(fSecondVPt,"stand. LYZ","lpe");
 370   legend2->AddEntry(fHistFlow,"new LYZ with calculated errors","lpe");
 371   legend2->Draw();
 372
 373   TCanvas *canvas5 = new TCanvas("canvas5","phi and Delta Phi",800,600);
 374   canvas5->Divide(3,1);
 375   canvas5->cd(1);
 376   fHistDeltaPhi->Draw();
 377   canvas5->cd(2);
 378   fHistPhiLYZ->Draw();
 379   canvas5->cd(3);
 380   fHistPhiEP->Draw();
 381
 382   cout<<".....finished"<<endl;
 383  }