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