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   fQsum(NULL),
  84   fQ2sum(0)
  85 {
  86
  87   // Constructor.
  88   fQsum = new TVector2();        // flow vector sum
  89 }
  90
  91
  92
  93  //-----------------------------------------------------------------------
  94
  95
  96  AliFlowAnalysisWithLYZEventPlane::~AliFlowAnalysisWithLYZEventPlane()
  97  {
  98    //destructor
  99    delete fQsum;
 100  }
 101
 102
 103 //-----------------------------------------------------------------------
 104 void AliFlowAnalysisWithLYZEventPlane::Init() {
 105
 106   //Initialise all histograms
 107   cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
 108
 109   //input histograms
 110   if (fSecondRunFile->IsZombie()){ //check if file exists
 111     cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
 112     exit(-1);
 113   } else if (fSecondRunFile->IsOpen()){
 114     cout<<"----secondRunFile is open----"<<endl;
 115     fSecondVPt = (TProfile*)fSecondRunFile->Get("Flow_Differential_Pt_LYZ"); //to compare to
 116   }
 117
 118   if (fFirstRunFile->IsZombie()){ //check if file exists
 119     cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
 120     exit(-1);
 121   } else if (fFirstRunFile->IsOpen()){
 122     cout<<"----firstRunFile is open----"<<endl<<endl;
 123     fFirstr0theta = (TProfile*)fFirstRunFile->Get("First_FlowPro_r0theta_LYZ");
 124   }
 125
 126
 127   //output file
 128   // ********make output file
 129   // analysis file (output)
 130   fOutFile = new TFile(fOutFileName.Data(),"RECREATE") ;
 131
 132   fCommonHists = new AliFlowCommonHist("LYZEP");
 133   fCommonHistsRes = new AliFlowCommonHistResults("LYZEP");
 134
 135   // output histograms
 136   Int_t iNtheta = AliFlowLYZConstants::kTheta;
 137   Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 138   Double_t  dPtMin = AliFlowCommonConstants::GetPtMin();
 139   Double_t  dPtMax = AliFlowCommonConstants::GetPtMax();
 140
 141
 142   fHistProFlow = new TProfile("FlowPro_VPt_LYZEP","FlowPro_VPt_LYZEP",iNbinsPt,dPtMin,dPtMax);
 143   fHistProFlow->SetXTitle("Pt");
 144   fHistProFlow->SetYTitle("v2 (%)");
 145
 146   fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
 147   fHistProWr->SetXTitle("Q");
 148   fHistProWr->SetYTitle("Wr");
 149
 150   fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
 151   fHistDeltaPhi->SetXTitle("DeltaPhi");
 152   fHistDeltaPhi->SetYTitle("Counts");
 153
 154   fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
 155   fHistPhiLYZ->SetXTitle("Phi from LYZ");
 156   fHistPhiLYZ->SetYTitle("Counts");
 157
 158   fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
 159   fHistPhiEP->SetXTitle("Phi from EP");
 160   fHistPhiEP->SetYTitle("Counts");
 161
 162
 163   fHistProR0theta  = new TProfile("FlowPro_r0theta_LYZEP","FlowPro_r0theta_LYZEP",iNtheta,-0.5,iNtheta-0.5);
 164   fHistProR0theta->SetXTitle("#theta");
 165   fHistProR0theta->SetYTitle("r_{0}^{#theta}");
 166
 167   fHistProReDtheta = new TProfile("FlowPro_ReDtheta_LYZEP","FlowPro_ReDtheta_LYZEP",iNtheta, -0.5, iNtheta-0.5);
 168   fHistProReDtheta->SetXTitle("#theta");
 169   fHistProReDtheta->SetYTitle("Re(D^{#theta})");
 170
 171   fHistProImDtheta = new TProfile("FlowPro_ImDtheta_LYZEP","FlowPro_ImDtheta_LYZEP",iNtheta, -0.5, iNtheta-0.5);
 172   fHistProImDtheta->SetXTitle("#theta");
 173   fHistProImDtheta->SetYTitle("Im(D^{#theta})");
 174
 175   fEventNumber = 0;  //set number of events to zero
 176
 177
 178 }
 179
 180 //-----------------------------------------------------------------------
 181
 182 void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* aLYZEP) {
 183
 184   //Get the event plane and weight for each event
 185   if (anEvent) {
 186
 187     //fill control histograms
 188     fCommonHists->FillControlHistograms(anEvent);
 189
 190     //get the Q vector from the FlowEvent
 191     AliFlowVector vQ = anEvent->GetQ();
 192     //Weight with the multiplicity
 193     Double_t dQX = 0.;
 194     Double_t dQY = 0.;
 195     if (vQ.GetMult()!=0.) {
 196       dQX = vQ.X()/vQ.GetMult();
 197       dQY = vQ.Y()/vQ.GetMult();
 198     } else {cerr<<"vQ.GetMult() is zero!"<<endl; }
 199     vQ.Set(dQX,dQY);
 200     //for chi calculation:
 201     *fQsum += vQ;
 202     fQ2sum += vQ.Mod2();
 203     cout<<"fQ2sum = "<<fQ2sum<<endl;
 204
 205     //call AliFlowLYZEventPlane::CalculateRPandW() here!
 206     aLYZEP->CalculateRPandW(vQ);
 207
 208     Double_t dWR = aLYZEP->GetWR();
 209     Double_t dRP = aLYZEP->GetPsi();
 210
 211     //fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
 212     fHistPhiLYZ->Fill(dRP);
 213
 214     //plot difference between event plane from EP-method and LYZ-method
 215     Double_t dRPEP = vQ.Phi()/2;                              //gives distribution from (0 to pi)
 216     //Double_t dRPEP = 0.5*TMath::ATan2(vQ.Y(),vQ.X());       //gives distribution from (-pi/2 to pi/2)
 217     //cout<<"dRPEP = "<< dRPEP <<endl;
 218     fHistPhiEP->Fill(dRPEP);
 219
 220     Double_t dDeltaPhi = dRPEP - dRP;
 221     if (dDeltaPhi < 0.) { dDeltaPhi += TMath::Pi(); }        //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
 222     //cout<<"dDeltaPhi = "<<dDeltaPhi<<endl;
 223     fHistDeltaPhi->Fill(dDeltaPhi);
 224
 225     //Flip sign of WR
 226     Double_t dLow = TMath::Pi()/4.;
 227     Double_t dHigh = 3.*(TMath::Pi()/4.);
 228     if ((dDeltaPhi > dLow) && (dDeltaPhi < dHigh)){
 229       dRP -= (TMath::Pi()/2);
 230       dWR = -dWR;
 231       cerr<<"*** dRP modified ***"<<endl;
 232     }
 233     fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
 234
 235     //calculate flow
 236     //loop over the tracks of the event
 237     Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 238     for (Int_t i=0;i<iNumberOfTracks;i++)
 239       {
 240         AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
 241         if (pTrack){
 242           if (pTrack->UseForDifferentialFlow()) {
 243             Double_t dPhi = pTrack->Phi();
 244             //if (dPhi<0.) fPhi+=2*TMath::Pi();
 245             //calculate flow v2:
 246             Double_t dv2 = dWR * TMath::Cos(2*(dPhi-dRP));
 247             Double_t dPt = pTrack->Pt();
 248             //fill histogram
 249             fHistProFlow->Fill(dPt,100*dv2);
 250           }
 251         }//track selected
 252       }//loop over tracks
 253
 254     fEventNumber++;
 255     cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 256   }
 257 }
 258
 259   //--------------------------------------------------------------------
 260 void AliFlowAnalysisWithLYZEventPlane::Finish() {
 261
 262   //plot histograms etc.
 263   cout<<"AliFlowAnalysisWithLYZEventPlane::Terminate()"<<endl;
 264   //constands:
 265   Double_t  dJ01 = 2.405;
 266   Int_t iNtheta = AliFlowLYZConstants::kTheta;
 267   Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 268
 269   //calculate dV the mean of dVtheta
 270   Double_t  dVtheta = 0;
 271   Double_t  dV = 0;
 272   for (Int_t theta=0;theta<iNtheta;theta++)     {
 273     Double_t dR0 = fFirstr0theta->GetBinContent(theta+1);
 274     if (dR0!=0.) { dVtheta = dJ01/dR0 ;}
 275     dV += dVtheta;
 276   }
 277   dV /= iNtheta;
 278
 279   //calculate dChi
 280   Double_t  dSigma2 = 0;
 281   Double_t  dChi= 0;
 282   if (fEventNumber!=0) {
 283     *fQsum /= fEventNumber;
 284     //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
 285     //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
 286     fQ2sum /= fEventNumber;
 287     cerr<<"fEventNumber = "<<fEventNumber<<endl;
 288     cerr<<"fQ2sum = "<<fQ2sum<<endl;
 289     dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
 290     if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
 291     else dChi = -1.;
 292     fCommonHistsRes->FillChi(dChi);
 293     cerr<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
 294   }
 295
 296   for(Int_t b=0;b<iNbinsPt;b++){
 297     Double_t dv2pro = 0.;
 298     Double_t dErr2difcomb = 0.;
 299     Double_t dErrdifcomb = 0.;
 300     if(fHistProFlow) {
 301       dv2pro = fHistProFlow->GetBinContent(b);
 302       //calculate error
 303       for (Int_t theta=0;theta<iNtheta;theta++) {
 304         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 305         Int_t iNprime = TMath::Nint(fHistProFlow->GetBinEntries(b));
 306         //cerr<<"iNprime = "<<iNprime<<endl;
 307         if (iNprime!=0) {
 308           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 309                                            TMath::Cos(dTheta));
 310           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 311                                           TMath::Cos(dTheta));
 312           dErr2difcomb += (TMath::Cos(dTheta)/(4*iNprime*TMath::BesselJ1(dJ01)*
 313                                                  TMath::BesselJ1(dJ01)))*
 314             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 315              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 316         } //if !=0
 317         //else { cout<<"iNprime = 0"<<endl; }
 318       } //loop over theta
 319
 320       if (dErr2difcomb!=0.) {
 321         dErr2difcomb /= iNtheta;
 322         dErrdifcomb = TMath::Sqrt(dErr2difcomb)*100;
 323         //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
 324       }
 325       else {dErrdifcomb = 0.; }
 326
 327       //fill TH1D
 328       fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb);
 329
 330     } //if fHistProFLow
 331     else  {
 332       cout << "Profile Hist missing" << endl;
 333       break;
 334     }
 335
 336   } //loop over b
 337
 338   // write to file
 339   fOutFile->Write();
 340
 341   cout<<"Making some plots to check the results:"<<endl<<endl;
 342
 343   TCanvas *canvas = new TCanvas("canvas","compare v2 vs pt",800,800);
 344   canvas->cd();
 345   fSecondVPt->SetLineColor(3);
 346   fSecondVPt->SetLineWidth(2);
 347   fSecondVPt->Draw();
 348   fHistFlow = fCommonHistsRes->GetHistDiffFlow();
 349   fHistFlow->Draw("SAME");
 350   // draw the legend
 351   TLegend *legend2 = new TLegend(0.6,0.65,0.88,0.85);
 352   legend2->SetTextFont(72);
 353   legend2->SetTextSize(0.04);
 354   legend2->AddEntry(fSecondVPt,"stand. LYZ","lpe");
 355   legend2->AddEntry(fHistFlow,"new LYZ with calculated errors","lpe");
 356   legend2->Draw();
 357
 358   TCanvas *canvas5 = new TCanvas("canvas5","phi and Delta Phi",800,600);
 359   canvas5->Divide(3,1);
 360   canvas5->cd(1);
 361   fHistDeltaPhi->Draw();
 362   canvas5->cd(2);
 363   fHistPhiLYZ->Draw();
 364   canvas5->cd(3);
 365   fHistPhiEP->Draw();
 366
 367   cout<<".....finished"<<endl;
 368  }
 369