PWG2/FLOW/AliFlowCommon/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 "TList.h"
  26 #include "TComplex.h"   //needed as include
  27 #include "TCanvas.h"   //needed as include
  28 #include "TLegend.h"   //needed as include
  29 #include "TProfile.h"  //needed as include
  30 #include "TVector2.h"
  31
  32 class TH1F;
  33
  34 #include "AliFlowLYZConstants.h"    //needed as include
  35 #include "AliFlowCommonConstants.h" //needed as include
  36 #include "AliFlowEventSimple.h"
  37 #include "AliFlowTrackSimple.h"
  38 #include "AliFlowCommonHist.h"
  39 #include "AliFlowCommonHistResults.h"
  40 #include "AliFlowLYZEventPlane.h"
  41 #include "AliFlowAnalysisWithLYZEventPlane.h"
  42
  43 class AliFlowVector;
  44
  45 // AliFlowAnalysisWithLYZEventPlane:
  46 //
  47 // Class to do flow analysis with the event plane from the LYZ method
  48 //
  49 // author: N. van der Kolk (kolk@nikhef.nl)
  50
  51
  52 ClassImp(AliFlowAnalysisWithLYZEventPlane)
  53
  54   //-----------------------------------------------------------------------
  55
  56  AliFlowAnalysisWithLYZEventPlane::AliFlowAnalysisWithLYZEventPlane():
  57    fHistList(NULL),
  58    fSecondRunList(NULL),
  59    fSecondReDtheta(NULL),
  60    fSecondImDtheta(NULL),
  61    fFirstr0theta(NULL),
  62    fHistProFlow(NULL),
  63    fHistProFlow2(NULL),
  64    fHistProWr(NULL),
  65    fHistProWrCorr(NULL),
  66    fHistQsumforChi(NULL),
  67    fHistDeltaPhi(NULL),
  68    fHistDeltaPhi2(NULL),
  69    fHistDeltaPhihere(NULL),
  70    fHistPhiEP(NULL),
  71    fHistPhiEPhere(NULL),
  72    fHistPhiLYZ(NULL),
  73    fHistPhiLYZ2(NULL),
  74    fCommonHists(NULL),
  75    fCommonHistsRes(NULL),
  76    fEventNumber(0),
  77    fQsum(NULL),
  78    fQ2sum(0)
  79 {
  80
  81   // Constructor.
  82   fQsum = new TVector2();        // flow vector sum
  83
  84   fHistList = new TList();
  85   fSecondRunList = new TList();
  86 }
  87
  88
  89
  90  //-----------------------------------------------------------------------
  91
  92
  93  AliFlowAnalysisWithLYZEventPlane::~AliFlowAnalysisWithLYZEventPlane()
  94  {
  95    //destructor
  96    delete fQsum;
  97    delete fHistList;
  98    delete fSecondRunList;
  99  }
 100
 101
 102 //-----------------------------------------------------------------------
 103
 104 void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TString* outputFileName)
 105 {
 106  //store the final results in output .root file
 107
 108   TFile *output = new TFile(outputFileName->Data(),"RECREATE");
 109   output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
 110   delete output;
 111 }
 112
 113 //-----------------------------------------------------------------------
 114
 115 void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TString outputFileName)
 116 {
 117  //store the final results in output .root file
 118
 119   TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 120   output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
 121   delete output;
 122 }
 123
 124 //-----------------------------------------------------------------------
 125 void AliFlowAnalysisWithLYZEventPlane::Init() {
 126
 127   //Initialise all histograms
 128   cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
 129
 130   //input histograms
 131   if (fSecondRunList) {
 132     fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZ");
 133     fHistList->Add(fSecondReDtheta);
 134
 135     fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZ");
 136     fHistList->Add(fSecondImDtheta);
 137
 138     fFirstr0theta = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZ");
 139     fHistList->Add(fFirstr0theta);
 140
 141     //warnings
 142     if (!fSecondReDtheta) {cout<<"fSecondReDtheta is NULL!"<<endl; }
 143     if (!fSecondImDtheta) {cout<<"fSecondImDtheta is NULL!"<<endl; }
 144     if (!fFirstr0theta)   {cout<<"fFirstr0theta is NULL!"<<endl; }
 145
 146   }
 147
 148   fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZEP");
 149   fHistList->Add(fCommonHists);
 150
 151   fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZEP");
 152   fHistList->Add(fCommonHistsRes);
 153
 154   Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 155   Double_t  dPtMin = AliFlowCommonConstants::GetPtMin();
 156   Double_t  dPtMax = AliFlowCommonConstants::GetPtMax();
 157
 158   fHistProFlow = new TProfile("FlowPro_VPt_LYZEP","FlowPro_VPt_LYZEP",iNbinsPt,dPtMin,dPtMax);
 159   fHistProFlow->SetXTitle("Pt");
 160   fHistProFlow->SetYTitle("v2 (%)");
 161   fHistList->Add(fHistProFlow);
 162
 163   fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
 164   fHistProWr->SetXTitle("Q");
 165   fHistProWr->SetYTitle("Wr");
 166   fHistList->Add(fHistProWr);
 167
 168   fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZEP","Flow_QsumforChi_LYZEP",3,-1.,2.);
 169   fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
 170   fHistQsumforChi->SetYTitle("value");
 171   fHistList->Add(fHistQsumforChi);
 172
 173   fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
 174   fHistDeltaPhi->SetXTitle("DeltaPhi");
 175   fHistDeltaPhi->SetYTitle("Counts");
 176   fHistList->Add(fHistDeltaPhi);
 177
 178   fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
 179   fHistPhiLYZ->SetXTitle("Phi from LYZ");
 180   fHistPhiLYZ->SetYTitle("Counts");
 181   fHistList->Add(fHistPhiLYZ);
 182
 183   fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
 184   fHistPhiEP->SetXTitle("Phi from EP");
 185   fHistPhiEP->SetYTitle("Counts");
 186   fHistList->Add(fHistPhiEP);
 187
 188   fEventNumber = 0;  //set number of events to zero
 189
 190 }
 191
 192 //-----------------------------------------------------------------------
 193
 194 void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* aLYZEP) {
 195
 196   //Get the event plane and weight for each event
 197   if (anEvent) {
 198
 199     //fill control histograms
 200     fCommonHists->FillControlHistograms(anEvent);
 201
 202     //get the Q vector from the FlowEvent
 203     AliFlowVector vQ = anEvent->GetQ();
 204     if (vQ.X()== 0. && vQ.Y()== 0. ) { cout<<"Q vector is NULL!"<<endl; }
 205     //Weight with the multiplicity
 206     Double_t dQX = 0.;
 207     Double_t dQY = 0.;
 208     if (vQ.GetMult()!=0.) {
 209       dQX = vQ.X()/vQ.GetMult();
 210       dQY = vQ.Y()/vQ.GetMult();
 211     } else {cerr<<"vQ.GetMult() is zero!"<<endl; }
 212     vQ.Set(dQX,dQY);
 213     //cout<<"vQ("<<dQX<<","<<dQY<<")"<<endl;
 214
 215     //for chi calculation:
 216     *fQsum += vQ;
 217     fHistQsumforChi->SetBinContent(1,fQsum->X());
 218     fHistQsumforChi->SetBinContent(2,fQsum->Y());
 219     fQ2sum += vQ.Mod2();
 220     fHistQsumforChi->SetBinContent(3,fQ2sum);
 221     //cout<<"fQ2sum = "<<fQ2sum<<endl;
 222
 223     //call AliFlowLYZEventPlane::CalculateRPandW() here!
 224     aLYZEP->CalculateRPandW(vQ);
 225
 226     Double_t dWR = aLYZEP->GetWR();
 227     Double_t dRP = aLYZEP->GetPsi();
 228
 229     //fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
 230     fHistPhiLYZ->Fill(dRP);
 231
 232     //plot difference between event plane from EP-method and LYZ-method
 233     Double_t dRPEP = vQ.Phi()/2;                              //gives distribution from (0 to pi)
 234     //Double_t dRPEP = 0.5*TMath::ATan2(vQ.Y(),vQ.X());       //gives distribution from (-pi/2 to pi/2)
 235     //cout<<"dRPEP = "<< dRPEP <<endl;
 236     fHistPhiEP->Fill(dRPEP);
 237
 238     Double_t dDeltaPhi = dRPEP - dRP;
 239     if (dDeltaPhi < 0.) { dDeltaPhi += TMath::Pi(); }        //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
 240     //cout<<"dDeltaPhi = "<<dDeltaPhi<<endl;
 241     fHistDeltaPhi->Fill(dDeltaPhi);
 242
 243     //Flip sign of WR
 244     Double_t dLow = TMath::Pi()/4.;
 245     Double_t dHigh = 3.*(TMath::Pi()/4.);
 246     if ((dDeltaPhi > dLow) && (dDeltaPhi < dHigh)){
 247       dRP -= (TMath::Pi()/2);
 248       dWR = -dWR;
 249       cerr<<"*** dRP modified ***"<<endl;
 250     }
 251     fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
 252
 253     //calculate flow
 254     //loop over the tracks of the event
 255     Int_t iNumberOfTracks = anEvent->NumberOfTracks();
 256     for (Int_t i=0;i<iNumberOfTracks;i++)
 257       {
 258         AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
 259         if (pTrack){
 260           if (pTrack->UseForDifferentialFlow()) {
 261             Double_t dPhi = pTrack->Phi();
 262             //if (dPhi<0.) fPhi+=2*TMath::Pi();
 263             //calculate flow v2:
 264             Double_t dv2 = dWR * TMath::Cos(2*(dPhi-dRP));
 265             Double_t dPt = pTrack->Pt();
 266             //fill histogram
 267             fHistProFlow->Fill(dPt,dv2);
 268           }
 269         }//track selected
 270       }//loop over tracks
 271
 272     fEventNumber++;
 273     //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
 274   }
 275 }
 276
 277   //--------------------------------------------------------------------
 278 void AliFlowAnalysisWithLYZEventPlane::Finish() {
 279
 280   //plot histograms etc.
 281   cout<<"AliFlowAnalysisWithLYZEventPlane::Finish()"<<endl;
 282
 283   //constants:
 284   Double_t  dJ01 = 2.405;
 285   Int_t iNtheta = AliFlowLYZConstants::kTheta;
 286   Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
 287   //set the event number
 288   SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries());
 289   //cout<<"number of events processed is "<<fEventNumber<<endl;
 290
 291   //set the sum of Q vectors
 292   fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
 293   SetQ2sum(fHistQsumforChi->GetBinContent(3));
 294
 295   //calculate dV the mean of dVtheta
 296   Double_t  dVtheta = 0;
 297   Double_t  dV = 0;
 298   for (Int_t theta=0;theta<iNtheta;theta++)     {
 299     Double_t dR0 = fFirstr0theta->GetBinContent(theta+1);
 300     if (dR0!=0.) { dVtheta = dJ01/dR0 ;}
 301     dV += dVtheta;
 302   }
 303   dV /= iNtheta;
 304
 305   //calculate dChi
 306   Double_t  dSigma2 = 0;
 307   Double_t  dChi= 0;
 308   if (fEventNumber!=0) {
 309     *fQsum /= fEventNumber;
 310     //cerr<<"fQsum->X() = "<<fQsum->X()<<endl;
 311     //cerr<<"fQsum->Y() = "<<fQsum->Y()<<endl;
 312     fQ2sum /= fEventNumber;
 313     //cerr<<"fEventNumber = "<<fEventNumber<<endl;
 314     //cerr<<"fQ2sum = "<<fQ2sum<<endl;
 315     dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
 316     //cerr<<"dSigma2"<<dSigma2<<endl;
 317     if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
 318     else dChi = -1.;
 319     fCommonHistsRes->FillChi(dChi);
 320     cerr<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
 321   }
 322
 323   for(Int_t b=0;b<iNbinsPt;b++){
 324     Double_t dv2pro = 0.;
 325     Double_t dErr2difcomb = 0.;
 326     Double_t dErrdifcomb = 0.;
 327     if(fHistProFlow) {
 328       dv2pro = fHistProFlow->GetBinContent(b);
 329       //calculate error
 330       for (Int_t theta=0;theta<iNtheta;theta++) {
 331         Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
 332         Int_t iNprime = TMath::Nint(fHistProFlow->GetBinEntries(b));
 333         //cerr<<"iNprime = "<<iNprime<<endl;
 334         if (iNprime!=0) {
 335           Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
 336                                            TMath::Cos(dTheta));
 337           Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
 338                                           TMath::Cos(dTheta));
 339           dErr2difcomb += (TMath::Cos(dTheta)/(4*iNprime*TMath::BesselJ1(dJ01)*
 340                                                  TMath::BesselJ1(dJ01)))*
 341             ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
 342              (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
 343         } //if !=0
 344         //else { cout<<"iNprime = 0"<<endl; }
 345       } //loop over theta
 346
 347       if (dErr2difcomb!=0.) {
 348         dErr2difcomb /= iNtheta;
 349         dErrdifcomb = TMath::Sqrt(dErr2difcomb);
 350         //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
 351       }
 352       else {dErrdifcomb = 0.; }
 353
 354       //fill TH1D
 355       fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb);
 356
 357     } //if fHistProFLow
 358     else  {
 359       cout << "Profile Hist missing" << endl;
 360       break;
 361     }
 362
 363   } //loop over b
 364
 365   cout<<".....finished"<<endl;
 366  }
 367