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