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