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