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