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