name change Int/Diff RP/POI
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowEventSimple.cxx
CommitLineData
f1d945a1 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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
f1d945a1 16#include "Riostream.h"
17#include "TObjArray.h"
26c4cbb9 18#include "TFile.h"
03a02aca 19#include "TList.h"
f1d945a1 20#include "TMath.h"
26c4cbb9 21#include "TH1F.h"
22#include "TH1D.h"
23#include "TProfile.h"
c076fda8 24#include "TBrowser.h"
f1d945a1 25#include "AliFlowVector.h"
26#include "AliFlowTrackSimple.h"
27#include "AliFlowEventSimple.h"
28
26c4cbb9 29/**************************************
30 * AliFlowEventSimple: A simple event *
31 * for flow analysis *
32 * *
33 * authors: Naomi van der Kolk *
34 * (kolk@nikhef.nl) *
35 * Ante Bilandzic *
36 * (anteb@nikhef.nl) *
37 * ***********************************/
38
f1d945a1 39ClassImp(AliFlowEventSimple)
40
41//-----------------------------------------------------------------------
42
46bec39c 43AliFlowEventSimple::AliFlowEventSimple():
44 fTrackCollection(NULL),
45 fNumberOfTracks(0),
1918addd 46 fEventNSelTracksRP(0),
7183fe85 47 fMCReactionPlaneAngle(0.),
c076fda8 48 fNumberOfTracksWrap(NULL),
1918addd 49 fEventNSelTracksRPWrap(NULL),
a12990bb 50 fMCReactionPlaneAngleWrap(NULL)
46bec39c 51{
52 cout << "AliFlowEventSimple: Default constructor to be used only by root for io" << endl;
53}
54
55//-----------------------------------------------------------------------
56
57AliFlowEventSimple::AliFlowEventSimple(Int_t aLenght):
58 fTrackCollection(NULL),
59 fNumberOfTracks(0),
1918addd 60 fEventNSelTracksRP(0),
7183fe85 61 fMCReactionPlaneAngle(0.),
c076fda8 62 fNumberOfTracksWrap(NULL),
1918addd 63 fEventNSelTracksRPWrap(NULL),
a12990bb 64 fMCReactionPlaneAngleWrap(NULL)
f1d945a1 65{
66 //constructor
e35ddff0 67 fTrackCollection = new TObjArray(aLenght) ;
f1d945a1 68}
69
70//-----------------------------------------------------------------------
71
e35ddff0 72AliFlowEventSimple::AliFlowEventSimple(const AliFlowEventSimple& anEvent):
bc6b015e 73 TObject(),
e35ddff0 74 fTrackCollection(anEvent.fTrackCollection),
75 fNumberOfTracks(anEvent.fNumberOfTracks),
1918addd 76 fEventNSelTracksRP(anEvent.fEventNSelTracksRP),
a12990bb 77 fMCReactionPlaneAngle(anEvent.fMCReactionPlaneAngle),
c076fda8 78 fNumberOfTracksWrap(anEvent.fNumberOfTracksWrap),
1918addd 79 fEventNSelTracksRPWrap(anEvent.fEventNSelTracksRPWrap),
a12990bb 80 fMCReactionPlaneAngleWrap(anEvent.fMCReactionPlaneAngleWrap)
f1d945a1 81{
82 //copy constructor
f1d945a1 83}
84
85//-----------------------------------------------------------------------
86
e35ddff0 87AliFlowEventSimple& AliFlowEventSimple::operator=(const AliFlowEventSimple& anEvent)
f1d945a1 88{
26c4cbb9 89 *fTrackCollection = *anEvent.fTrackCollection ;
e35ddff0 90 fNumberOfTracks = anEvent.fNumberOfTracks;
1918addd 91 fEventNSelTracksRP = anEvent.fEventNSelTracksRP;
a12990bb 92 fMCReactionPlaneAngle = anEvent.fMCReactionPlaneAngle;
c076fda8 93 fNumberOfTracksWrap = anEvent.fNumberOfTracksWrap;
1918addd 94 fEventNSelTracksRPWrap = anEvent.fEventNSelTracksRPWrap;
a12990bb 95 fMCReactionPlaneAngleWrap=anEvent.fMCReactionPlaneAngleWrap;
c076fda8 96
f1d945a1 97 return *this;
f1d945a1 98}
99
f1d945a1 100//-----------------------------------------------------------------------
101
102AliFlowEventSimple::~AliFlowEventSimple()
103{
104 //destructor
c076fda8 105 if (fTrackCollection) fTrackCollection->Delete(); delete fTrackCollection;
106 if (fNumberOfTracksWrap) delete fNumberOfTracksWrap;
1918addd 107 if (fEventNSelTracksRPWrap) delete fEventNSelTracksRPWrap;
a12990bb 108 if (fMCReactionPlaneAngleWrap) delete fMCReactionPlaneAngleWrap;
f1d945a1 109}
110
111//-----------------------------------------------------------------------
112
113AliFlowTrackSimple* AliFlowEventSimple::GetTrack(Int_t i)
114{
115 //get track i from collection
e35ddff0 116 AliFlowTrackSimple* pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i) ;
117 return pTrack;
f1d945a1 118}
119
120//-----------------------------------------------------------------------
03a02aca 121AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
f1d945a1 122{
ae733b3b 123 // calculate Q-vector in harmonic n without weights (default harmonic n=2)
e35ddff0 124 Double_t dQX = 0.;
125 Double_t dQY = 0.;
126 AliFlowVector vQ;
127 vQ.Set(0.,0.);
9825d4a9 128
129 Int_t iOrder = n;
ae733b3b 130 Double_t iUsedTracks = 0;
26c4cbb9 131 Double_t dPhi=0.;
ae733b3b 132 Double_t dPt=0.;
133 Double_t dEta=0.;
26c4cbb9 134
135 AliFlowTrackSimple* pTrack = NULL;
136
137 Int_t nBinsPhi=0;
138 Double_t dBinWidthPt=0.;
ae733b3b 139 Double_t dPtMin=0.;
26c4cbb9 140 Double_t dBinWidthEta=0.;
ae733b3b 141 Double_t dEtaMin=0.;
26c4cbb9 142
ae733b3b 143 Double_t wPhi=1.; // weight Phi
144 Double_t wPt=1.; // weight Pt
145 Double_t wEta=1.; // weight Eta
26c4cbb9 146
03a02aca 147 TH1F *phiWeights = NULL;
ae733b3b 148 TH1D *ptWeights = NULL;
03a02aca 149 TH1D *etaWeights = NULL;
150
ae733b3b 151 Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 2)
152 Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 3)
153 Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 4)
154 Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 5)
155 Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 6)
156 Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 7)
157 Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 8)
158
03a02aca 159 if(weightsList)
26c4cbb9 160 {
03a02aca 161 if(usePhiWeights)
26c4cbb9 162 {
03a02aca 163 phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
164 if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
165 }
166 if(usePtWeights)
26c4cbb9 167 {
03a02aca 168 ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
169 if(ptWeights)
170 {
171 dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
ae733b3b 172 dPtMin = (ptWeights->GetXaxis())->GetXmin();
03a02aca 173 }
174 }
175 if(useEtaWeights)
26c4cbb9 176 {
03a02aca 177 etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
178 if(etaWeights)
179 {
180 dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
ae733b3b 181 dEtaMin = (etaWeights->GetXaxis())->GetXmin();
03a02aca 182 }
183 }
184 } // end of if(weightsList)
26c4cbb9 185
03a02aca 186 // loop over tracks
26c4cbb9 187 for(Int_t i=0;i<fNumberOfTracks;i++)
188 {
189 pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i);
190 if(pTrack)
191 {
1918addd 192 if(pTrack->InRPSelection())
f1d945a1 193 {
26c4cbb9 194 dPhi = pTrack->Phi();
195 dPt = pTrack->Pt();
196 dEta = pTrack->Eta();
ae733b3b 197
198 // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
199 if(phiWeights && nBinsPhi)
26c4cbb9 200 {
ae733b3b 201 wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
26c4cbb9 202 }
ae733b3b 203 // determine v'(pt) weight:
204 if(ptWeights && dBinWidthPt)
26c4cbb9 205 {
ae733b3b 206 wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
26c4cbb9 207 }
ae733b3b 208 // determine v'(eta) weight:
209 if(etaWeights && dBinWidthEta)
26c4cbb9 210 {
ae733b3b 211 wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
03a02aca 212 }
ae733b3b 213
214 // building up the weighted Q-vector:
26c4cbb9 215 dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
216 dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
ae733b3b 217
218 // weighted multiplicity:
219 iUsedTracks+=wPhi*wPt*wEta;
220
221 // weights raised to various powers are summed up:
222 dSumOfWeightsToPower2+=pow(wPhi*wPt*wEta, 2);
223 dSumOfWeightsToPower3+=pow(wPhi*wPt*wEta, 3);
224 dSumOfWeightsToPower4+=pow(wPhi*wPt*wEta, 4);
225 dSumOfWeightsToPower5+=pow(wPhi*wPt*wEta, 5);
226 dSumOfWeightsToPower6+=pow(wPhi*wPt*wEta, 6);
227 dSumOfWeightsToPower7+=pow(wPhi*wPt*wEta, 7);
228 dSumOfWeightsToPower8+=pow(wPhi*wPt*wEta, 8);
03a02aca 229
1918addd 230 } // end of if (pTrack->InRPSelection())
ae733b3b 231 } // end of if (pTrack)
26c4cbb9 232 else {cerr << "no particle!!!"<<endl;}
ae733b3b 233 } // loop over particles
26c4cbb9 234
e35ddff0 235 vQ.Set(dQX,dQY);
236 vQ.SetMult(iUsedTracks);
ae733b3b 237 vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
238 vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
239 vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
240 vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
241 vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
242 vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
243 vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
26c4cbb9 244
e35ddff0 245 return vQ;
f1d945a1 246
5fef318d 247}
248
395fadba 249//-----------------------------------------------------------------------
250AliFlowVector AliFlowEventSimple::GetQsub(Double_t etaMin, Double_t etaMax, Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights)
251{
252 //for eta subevents
253
254 // calculate Q-vector in harmonic n without weights (default harmonic n=2)
255 Double_t dQX = 0.;
256 Double_t dQY = 0.;
257 AliFlowVector vQ;
258 vQ.Set(0.,0.);
259
260 Int_t iOrder = n;
261 Double_t iUsedTracks = 0;
262 Double_t dPhi = 0.;
263 Double_t dPt = 0.;
264 Double_t dEta = 0.;
265
266 AliFlowTrackSimple* pTrack = NULL;
267
268 Int_t nBinsPhi = 0;
269 Double_t dBinWidthPt = 0.;
270 Double_t dPtMin = 0.;
271 Double_t dBinWidthEta= 0.;
272 Double_t dEtaMin = 0.;
273
274 Double_t wPhi = 1.; // weight Phi
275 Double_t wPt = 1.; // weight Pt
276 Double_t wEta = 1.; // weight Eta
277
278 TH1F *phiWeights = NULL;
279 TH1D *ptWeights = NULL;
280 TH1D *etaWeights = NULL;
281
282 Double_t dSumOfWeightsToPower2 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 2)
283 Double_t dSumOfWeightsToPower3 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 3)
284 Double_t dSumOfWeightsToPower4 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 4)
285 Double_t dSumOfWeightsToPower5 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 5)
286 Double_t dSumOfWeightsToPower6 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 6)
287 Double_t dSumOfWeightsToPower7 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 7)
288 Double_t dSumOfWeightsToPower8 = 0.; // sum_{i=1}^{n} pow((wPhi*wPt*wEta)_i, 8)
289
290 if(weightsList)
291 {
292 if(usePhiWeights)
293 {
294 phiWeights = dynamic_cast<TH1F *>(weightsList->FindObject("phi_weights"));
295 if(phiWeights) nBinsPhi = phiWeights->GetNbinsX();
296 }
297 if(usePtWeights)
298 {
299 ptWeights = dynamic_cast<TH1D *>(weightsList->FindObject("pt_weights"));
300 if(ptWeights)
301 {
302 dBinWidthPt = ptWeights->GetBinWidth(1); // assuming that all bins have the same width
303 dPtMin = (ptWeights->GetXaxis())->GetXmin();
304 }
305 }
306 if(useEtaWeights)
307 {
308 etaWeights = dynamic_cast<TH1D *>(weightsList->FindObject("eta_weights"));
309 if(etaWeights)
310 {
311 dBinWidthEta = etaWeights->GetBinWidth(1); // assuming that all bins have the same width
312 dEtaMin = (etaWeights->GetXaxis())->GetXmin();
313 }
314 }
315 } // end of if(weightsList)
316
317 // loop over tracks
318 for(Int_t i=0;i<fNumberOfTracks;i++)
319 {
320 pTrack = (AliFlowTrackSimple*)TrackCollection()->At(i);
321 if(pTrack)
322 {
1918addd 323 if(pTrack->InRPSelection())
395fadba 324 {
325 dPhi = pTrack->Phi();
326 dPt = pTrack->Pt();
327 dEta = pTrack->Eta();
328 if (dEta>etaMin && dEta<etaMax) {
329 // determine Phi weight: (to be improved, I should here only access it + the treatment of gaps in the if statement)
330 if(phiWeights && nBinsPhi)
331 {
332 wPhi = phiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*nBinsPhi/TMath::TwoPi())));
333 }
334 // determine v'(pt) weight:
335 if(ptWeights && dBinWidthPt)
336 {
337 wPt=ptWeights->GetBinContent(1+(Int_t)(TMath::Floor((dPt-dPtMin)/dBinWidthPt)));
338 }
339 // determine v'(eta) weight:
340 if(etaWeights && dBinWidthEta)
341 {
342 wEta=etaWeights->GetBinContent(1+(Int_t)(TMath::Floor((dEta-dEtaMin)/dBinWidthEta)));
343 }
344
345 // building up the weighted Q-vector:
346 dQX += wPhi*wPt*wEta*TMath::Cos(iOrder*dPhi);
347 dQY += wPhi*wPt*wEta*TMath::Sin(iOrder*dPhi);
348
349 // weighted multiplicity:
350 iUsedTracks+=wPhi*wPt*wEta;
351
352 // weights raised to various powers are summed up:
353 dSumOfWeightsToPower2+=pow(wPhi*wPt*wEta, 2);
354 dSumOfWeightsToPower3+=pow(wPhi*wPt*wEta, 3);
355 dSumOfWeightsToPower4+=pow(wPhi*wPt*wEta, 4);
356 dSumOfWeightsToPower5+=pow(wPhi*wPt*wEta, 5);
357 dSumOfWeightsToPower6+=pow(wPhi*wPt*wEta, 6);
358 dSumOfWeightsToPower7+=pow(wPhi*wPt*wEta, 7);
359 dSumOfWeightsToPower8+=pow(wPhi*wPt*wEta, 8);
360 } // end of if dEta in eta range
1918addd 361 } // end of if (pTrack->InRPSelection())
395fadba 362 } // end of if (pTrack)
363 else {cerr << "no particle!!!"<<endl;}
364 } // loop over particles
365
366 vQ.Set(dQX,dQY);
367 vQ.SetMult(iUsedTracks);
368 vQ.SetSumOfWeightsToPower2(dSumOfWeightsToPower2);
369 vQ.SetSumOfWeightsToPower3(dSumOfWeightsToPower3);
370 vQ.SetSumOfWeightsToPower4(dSumOfWeightsToPower4);
371 vQ.SetSumOfWeightsToPower5(dSumOfWeightsToPower5);
372 vQ.SetSumOfWeightsToPower6(dSumOfWeightsToPower6);
373 vQ.SetSumOfWeightsToPower7(dSumOfWeightsToPower7);
374 vQ.SetSumOfWeightsToPower8(dSumOfWeightsToPower8);
375
376 return vQ;
377
378}
379
380
c076fda8 381//-----------------------------------------------------------------------
382void AliFlowEventSimple::Print(Option_t *option) const
383{
384 // -*-*-*-*-*Print some global quantities for this histogram collection class *-*-*-*-*-*-*-*
385 // ===============================================
386 // printf( "TH1.Print Name = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
7183fe85 387 printf( "Class.Print Name = %s, Total number of tracks= %d, Number of selected tracks= %d, MC EventPlaneAngle= %f",
1918addd 388 GetName(),fNumberOfTracks, fEventNSelTracksRP, fMCReactionPlaneAngle );
c076fda8 389
390 if (fTrackCollection) {
391 fTrackCollection->Print(option);
392 }
393 else {
394 printf( "Empty track collection \n");
395 }
396}
397
398//-----------------------------------------------------------------------
399 void AliFlowEventSimple::Browse(TBrowser *b)
400{
401 if (!b) return;
402 if (!fNumberOfTracksWrap) {
403 fNumberOfTracksWrap = new TParameter<int>("fNumberOfTracks", fNumberOfTracks);
404 b->Add(fNumberOfTracksWrap);
405 }
1918addd 406 if (!fEventNSelTracksRPWrap) {
407 fEventNSelTracksRPWrap = new TParameter<int>("fEventNSelTracksRP", fEventNSelTracksRP);
408 b->Add(fEventNSelTracksRPWrap);
c076fda8 409 }
a12990bb 410 if (!fMCReactionPlaneAngleWrap) {
7183fe85 411 fMCReactionPlaneAngleWrap = new TParameter<double>(" fMCReactionPlaneAngle", fMCReactionPlaneAngle);
a12990bb 412 b->Add( fMCReactionPlaneAngleWrap);
413 }
c076fda8 414 if (fTrackCollection) b->Add(fTrackCollection,"AliFlowTracksSimple");
415}
416
395fadba 417