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