Update master to aliroot
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliEventplane.cxx
CommitLineData
ce7adfe9 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// Class AliEventplane
18// author: Alberica Toia, Johanna Gramling
19//*****************************************************
20/// A container for the event plane stored in AOD in ESD
21
fd5c821a 22#include "AliLog.h"
ce7adfe9 23#include "AliEventplane.h"
24#include "TVector2.h"
aa260dd1 25#include "AliVTrack.h"
ce7adfe9 26#include "TObjArray.h"
27#include "TArrayF.h"
fd5c821a 28#include "AliVEvent.h"
29#include "AliVVZERO.h"
ce7adfe9 30
31ClassImp(AliEventplane)
32
33AliEventplane::AliEventplane() : TNamed("Eventplane", "Eventplane"),
34 fQVector(0),
35 fQContributionX(0),
36 fQContributionY(0),
e51055a0 37 fQContributionXsub1(0),
38 fQContributionYsub1(0),
39 fQContributionXsub2(0),
40 fQContributionYsub2(0),
3e64f60c 41 fEventplaneQ(-1),
ce7adfe9 42 fQsub1(0),
43 fQsub2(0),
44 fQsubRes(0)
45{
46 /// constructor
3e64f60c 47 fQContributionX = new TArrayF(0);
48 fQContributionY = new TArrayF(0);
e51055a0 49 fQContributionXsub1 = new TArrayF(0);
50 fQContributionYsub1 = new TArrayF(0);
51 fQContributionXsub2 = new TArrayF(0);
52 fQContributionYsub2 = new TArrayF(0);
cc4d0e9c 53 for(Int_t i = 0; i < 11; ++i) {
0643f5a7 54 fMeanX2[i]=fMeanY2[i]=0.;
55 fAPlus[i]=fAMinus[i]=1.;
56 fLambdaPlus[i]=fLambdaMinus[i]=0.;
57 fCos8Psi[i]=0.;
58 }
ce7adfe9 59}
60
61AliEventplane::AliEventplane(const AliEventplane& ep) :
fe536aed 62 TNamed(ep),
ce7adfe9 63 fQVector(0),
64 fQContributionX(0),
65 fQContributionY(0),
e51055a0 66 fQContributionXsub1(0),
67 fQContributionYsub1(0),
68 fQContributionXsub2(0),
69 fQContributionYsub2(0),
ce7adfe9 70 fEventplaneQ(0),
71 fQsub1(0),
72 fQsub2(0),
73 fQsubRes(0)
74{
75 /// Copy constructor
76 ((AliEventplane &) ep).CopyEP(*this);
77}
78
79AliEventplane& AliEventplane::operator=(const AliEventplane& ep)
80{
81 /// Assignment operator
fe536aed 82 if (this!=&ep){
83 TNamed::operator=(ep);
ce7adfe9 84 ((AliEventplane &) ep).CopyEP(*this);
fe536aed 85 }
ce7adfe9 86 return *this;
87}
88
fe536aed 89void AliEventplane::CopyEP(AliEventplane& ep) const
ce7adfe9 90{ // copy function
91
92 AliEventplane& target = (AliEventplane &) ep;
fe536aed 93
94 if(fQContributionX)
95 {
96 if(ep.fQContributionX)
97 {
98 *(ep.fQContributionX) = *fQContributionX;
99 }
100 else
101 {
102 ep.fQContributionX = new TArrayF(*fQContributionX);
103 }
104 }
105 else
106 {
107 if(ep.fQContributionX) delete ep.fQContributionX;
108 ep.fQContributionX = 0;
109 }
110
111 if(fQContributionY)
112 {
113 if(ep.fQContributionY)
114 {
115 *(ep.fQContributionY) = *fQContributionY;
116 }
117 else
118 {
119 ep.fQContributionY = new TArrayF(*fQContributionY);
120 }
121 }
122 else
123 {
124 if(ep.fQContributionY) delete ep.fQContributionY;
125 ep.fQContributionY = 0;
126 }
127
128
129 if(fQContributionXsub1)
130 {
131 if(ep.fQContributionXsub1)
132 {
133 *(ep.fQContributionXsub1) = *fQContributionXsub1;
134 }
135 else
136 {
137 ep.fQContributionXsub1 = new TArrayF(*fQContributionXsub1);
138 }
139 }
140 else
141 {
142 if(ep.fQContributionXsub1) delete ep.fQContributionXsub1;
143 ep.fQContributionXsub1 = 0;
144 }
145
146 if(fQContributionYsub1)
147 {
148 if(ep.fQContributionYsub1)
149 {
150 *(ep.fQContributionYsub1) = *fQContributionYsub1;
151 }
152 else
153 {
154 ep.fQContributionYsub1 = new TArrayF(*fQContributionYsub1);
155 }
156 }
157 else
158 {
159 if(ep.fQContributionYsub1) delete ep.fQContributionYsub1;
160 ep.fQContributionYsub1 = 0;
161 }
162
163
164 if(fQContributionXsub2)
165 {
166 if(ep.fQContributionXsub2)
167 {
168 *(ep.fQContributionXsub2) = *fQContributionXsub2;
169 }
170 else
171 {
172 ep.fQContributionXsub2 = new TArrayF(*fQContributionXsub2);
173 }
174 }
175 else
176 {
177 if(ep.fQContributionXsub2) delete ep.fQContributionXsub2;
178 ep.fQContributionXsub2 = 0;
179 }
180
181 if(fQContributionYsub2)
182 {
183 if(ep.fQContributionYsub2)
184 {
185 *(ep.fQContributionYsub2) = *fQContributionYsub2;
186 }
187 else
188 {
189 ep.fQContributionYsub2 = new TArrayF(*fQContributionYsub2);
190 }
191 }
192 else
193 {
194 if(ep.fQContributionYsub2) delete ep.fQContributionYsub2;
195 ep.fQContributionYsub2 = 0;
196 }
197
ce7adfe9 198 if (fEventplaneQ)
199 target.fEventplaneQ = fEventplaneQ;
fe536aed 200
ce7adfe9 201 if (fQVector)
202 target.fQVector = dynamic_cast<TVector2*> (fQVector->Clone());
203 if (fQsub1)
204 target.fQsub1 = dynamic_cast<TVector2*> (fQsub1->Clone());
205 if (fQsub2)
206 target.fQsub2 = dynamic_cast<TVector2*> (fQsub2->Clone());
207 if (fQsubRes)
208 target.fQsubRes = fQsubRes;
0643f5a7 209
cc4d0e9c 210 for(Int_t i = 0; i < 11; ++i) {
0643f5a7 211 target.fMeanX2[i]=fMeanX2[i];
212 target.fMeanY2[i]=fMeanY2[i];
213 target.fAPlus[i]=fAPlus[i];
214 target.fAMinus[i]=fAMinus[i];
215 target.fLambdaPlus[i]=fLambdaPlus[i];
216 target.fLambdaMinus[i]=fLambdaMinus[i];
217 target.fCos8Psi[i]=fCos8Psi[i];
218 }
ce7adfe9 219}
220
221AliEventplane::~AliEventplane()
222{
223 /// destructor
224 if (fQContributionX){
225 delete fQContributionX;
226 fQContributionX = 0;
227 }
228 if (fQContributionY){
229 delete fQContributionY;
230 fQContributionY = 0;
231 }
e51055a0 232 if (fQContributionXsub1){
233 delete fQContributionXsub1;
234 fQContributionXsub1 = 0;
235 }
236 if (fQContributionYsub1){
237 delete fQContributionYsub1;
238 fQContributionYsub1 = 0;
239 }
240 if (fQContributionXsub2){
241 delete fQContributionXsub2;
242 fQContributionXsub2 = 0;
243 }
244 if (fQContributionYsub2){
245 delete fQContributionYsub2;
246 fQContributionYsub2 = 0;
247 }
ce7adfe9 248 if (fQVector){
249 delete fQVector;
250 fQVector = 0;
251 }
252 if (fQsub1){
253 delete fQsub1;
254 fQsub1 = 0;
255 }
256 if (fQsub2){
257 delete fQsub2;
258 fQsub2 = 0;
259 }
260}
261
262TVector2* AliEventplane::GetQVector()
263{
264 return fQVector;
265}
266
fd5c821a 267Double_t AliEventplane::GetEventplane(const char *x, const AliVEvent *event, Int_t harmonic) const
ce7adfe9 268{
269 TString method = x;
fe509094 270 Double_t qx = 0, qy = 0;
ce7adfe9 271 if(method.CompareTo("Q")==0) return fEventplaneQ;
cc4d0e9c 272 else if(method.CompareTo("V0A")==0) return CalculateVZEROEventPlane(event, 8, harmonic, qx, qy);
273 else if(method.CompareTo("V0C")==0) return CalculateVZEROEventPlane(event, 9, harmonic, qx, qy);
274 else if(method.CompareTo("V0")==0) return CalculateVZEROEventPlane(event,10, harmonic, qx, qy);
fd5c821a 275
276 return -1000.;
277}
278
fe509094 279Double_t AliEventplane::CalculateVZEROEventPlane(const AliVEvent * event, Int_t firstRing, Int_t lastRing, Int_t harmonic, Double_t &qxTot, Double_t &qyTot) const
fd5c821a 280{
cc4d0e9c 281 // Calculate the VZERO event-plane by summing rings
282 // The flatenning is done on every ring separately
fd5c821a 283 if(!event) {
284 AliError("No Event received");
285 return -1000.;
286 }
287 AliVVZERO *vzeroData = event->GetVZEROData();
288 if(!vzeroData) {
289 AliError("Enable to get VZERO Data");
290 return -1000.;
291 }
292 if(harmonic <= 0) {
293 AliError("Required harmonic is less or equal to 0");
294 return -1000.;
295 }
296
fe509094 297 qxTot=qyTot=0.;
0643f5a7 298 Double_t qx, qy;
299 Double_t totMult = 0.;
300 for(Int_t ring = firstRing; ring <=lastRing; ++ring) {
301 qx=qy=0.;
0643f5a7 302 Double_t multRing = 0.;
303 for(Int_t iCh = ring*8; iCh < (ring+1)*8; ++iCh) {
304 Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
305 Double_t mult = event->GetVZEROEqMultiplicity(iCh);
306 qx += mult*TMath::Cos(harmonic*phi);
307 qy += mult*TMath::Sin(harmonic*phi);
308 multRing += mult;
309 }
310
311 if (multRing < 1e-6) continue;
312 totMult += multRing;
313
314 if (harmonic == 2) {
315 // Recentering
316 Double_t qxPrime = qx - fMeanX2[ring];
317 Double_t qyPrime = qy - fMeanY2[ring];
318 // Twist
cc4d0e9c 319 Double_t trans[2][2];
320 trans[0][0] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
321 trans[0][1] = -fLambdaMinus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
322 trans[1][0] = -fLambdaPlus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
323 trans[1][1] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
0643f5a7 324 Double_t qxSeconde = trans[0][0]*qxPrime + trans[0][1]*qyPrime;
325 Double_t qySeconde = trans[1][0]*qxPrime + trans[1][1]*qyPrime;
326 // Rescaling
327 Double_t qxTierce = qxSeconde/fAPlus[ring];
328 Double_t qyTierce = qySeconde/fAMinus[ring];
93736b53 329 // 4th Fourier momenta in order to flatten the EP within a sector
330 Double_t psi = TMath::ATan2(qyTierce,qxTierce)/harmonic;
331 Double_t deltaPsi = (fCos8Psi[ring]*TMath::Sin(2.*4.*psi))/2.;
332 Double_t qxQuarte = qxTierce*TMath::Cos(2.*deltaPsi) - qyTierce*TMath::Sin(2.*deltaPsi);
333 Double_t qyQuarte = qxTierce*TMath::Sin(2.*deltaPsi) + qyTierce*TMath::Cos(2.*deltaPsi);
334 qxTot += qxQuarte;
335 qyTot += qyQuarte;
0643f5a7 336 }
337 else {
338 qxTot += qx;
339 qyTot += qy;
340 }
341 }
fd5c821a 342
0643f5a7 343 // In case of no hits return an invalid event-plane angle
344 if (totMult<1e-6) return -999;
fd5c821a 345
93736b53 346 return (TMath::ATan2(qyTot,qxTot)/harmonic);
ce7adfe9 347}
348
cc4d0e9c 349Double_t AliEventplane::CalculateVZEROEventPlane(const AliVEvent * event, Int_t ring, Int_t harmonic, Double_t &qx, Double_t &qy) const
350{
351 // Calculate the VZERO event-plane from an individual ring
352 // Ring 8 - total V0A (rings 4-7), Ring 9 - total V0C (rings 0-3)
353 // Ring 10 - total V0 (rings 0-7)
354 if(!event) {
355 AliError("No Event received");
356 return -1000.;
357 }
358 AliVVZERO *vzeroData = event->GetVZEROData();
359 if(!vzeroData) {
360 AliError("Enable to get VZERO Data");
361 return -1000.;
362 }
363 if(harmonic <= 0) {
364 AliError("Required harmonic is less or equal to 0");
365 return -1000.;
366 }
367
6e4f0e6b 368 Int_t firstCh=-1,lastCh=-1;
369 if (ring < 8) {
370 firstCh = ring*8;
371 lastCh = (ring+1)*8;
372 }
373 else if (ring == 8) {
374 firstCh = 32;
375 lastCh = 64;
376 }
377 else if (ring == 9) {
378 firstCh = 0;
379 lastCh = 32;
380 }
381 else if (ring == 10) {
382 firstCh = 0;
383 lastCh = 64;
384 }
cc4d0e9c 385 qx=qy=0.;
386 Double_t multRing = 0.;
6e4f0e6b 387 for(Int_t iCh = firstCh; iCh < lastCh; ++iCh) {
cc4d0e9c 388 Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
389 Double_t mult = event->GetVZEROEqMultiplicity(iCh);
390 qx += mult*TMath::Cos(harmonic*phi);
391 qy += mult*TMath::Sin(harmonic*phi);
392 multRing += mult;
393 }
394
395 // In case of no hits return an invalid event-plane angle
396 if (multRing < 1e-6) return -999;
397
398 if (harmonic == 2) {
399 // Recentering
400 Double_t qxPrime = qx - fMeanX2[ring];
401 Double_t qyPrime = qy - fMeanY2[ring];
402 // Twist
403 Double_t trans[2][2];
404 trans[0][0] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
405 trans[0][1] = -fLambdaMinus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
406 trans[1][0] = -fLambdaPlus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
407 trans[1][1] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
408 Double_t qxSeconde = trans[0][0]*qxPrime + trans[0][1]*qyPrime;
409 Double_t qySeconde = trans[1][0]*qxPrime + trans[1][1]*qyPrime;
410 // Rescaling
411 Double_t qxTierce = qxSeconde/fAPlus[ring];
412 Double_t qyTierce = qySeconde/fAMinus[ring];
413 // 4th Fourier momenta in order to flatten the EP within a sector
414 Double_t psi = TMath::ATan2(qyTierce,qxTierce)/harmonic;
415 Double_t deltaPsi = (fCos8Psi[ring]*TMath::Sin(2.*4.*psi))/2.;
416 Double_t qxQuarte = qxTierce*TMath::Cos(2.*deltaPsi) - qyTierce*TMath::Sin(2.*deltaPsi);
417 Double_t qyQuarte = qxTierce*TMath::Sin(2.*deltaPsi) + qyTierce*TMath::Cos(2.*deltaPsi);
418 qx = qxQuarte;
419 qy = qyQuarte;
420 }
421
422 return (TMath::ATan2(qy,qx)/harmonic);
423}
424
0643f5a7 425void AliEventplane::SetVZEROEPParams(Int_t ring,
426 Double_t meanX2, Double_t meanY2,
427 Double_t aPlus, Double_t aMinus,
428 Double_t lambdaPlus, Double_t lambdaMinus,
429 Double_t cos8Psi)
430{
431 // Set the VZERO event-plane
432 // flatenning parameters
cc4d0e9c 433 if (ring < 0 || ring >= 11) AliFatal(Form("Bad ring index (%d)",ring));
0643f5a7 434
435 fMeanX2[ring] = meanX2;
436 fMeanY2[ring] = meanY2;
437 fAPlus[ring] = aPlus;
438 fAMinus[ring] = aMinus;
439 fLambdaPlus[ring] = lambdaPlus;
440 fLambdaMinus[ring] = lambdaMinus;
441 fCos8Psi[ring] = cos8Psi;
442}
fd5c821a 443
ce7adfe9 444TVector2* AliEventplane::GetQsub1()
445{
446 return fQsub1;
447}
448
449TVector2* AliEventplane::GetQsub2()
450{
451 return fQsub2;
452}
453
454Double_t AliEventplane::GetQsubRes()
455{
456 return fQsubRes;
457}
458
459Bool_t AliEventplane::IsEventInEventplaneClass(Double_t a, Double_t b, const char *x)
460{
461 TString method = x;
462 if ((method.CompareTo("Q")==0) && (fEventplaneQ >=a && fEventplaneQ < b)) return kTRUE;
463 else return kFALSE;
464}
465
aa260dd1 466Double_t AliEventplane::GetQContributionX(AliVTrack* track)
ce7adfe9 467{
468 return fQContributionX->GetAt(track->GetID());
469}
470
aa260dd1 471Double_t AliEventplane::GetQContributionY(AliVTrack* track)
ce7adfe9 472{
473 return fQContributionY->GetAt(track->GetID());
474}
3e64f60c 475
e51055a0 476Double_t AliEventplane::GetQContributionXsub1(AliVTrack* track)
477{
478 return fQContributionXsub1->GetAt(track->GetID());
479}
480
481Double_t AliEventplane::GetQContributionYsub1(AliVTrack* track)
482{
483 return fQContributionYsub1->GetAt(track->GetID());
484}
485
486Double_t AliEventplane::GetQContributionXsub2(AliVTrack* track)
487{
488 return fQContributionXsub2->GetAt(track->GetID());
489}
490
491Double_t AliEventplane::GetQContributionYsub2(AliVTrack* track)
492{
493 return fQContributionYsub2->GetAt(track->GetID());
494}
495
3e64f60c 496void AliEventplane::Reset()
497{
498 delete fQVector; fQVector=0;
499 fQContributionX->Reset();
500 fQContributionY->Reset();
e51055a0 501 fQContributionXsub1->Reset();
502 fQContributionYsub1->Reset();
503 fQContributionXsub2->Reset();
504 fQContributionYsub2->Reset();
3e64f60c 505 fEventplaneQ = -1;
506 delete fQsub1; fQsub1=0;
507 delete fQsub2; fQsub2=0;
508 fQsubRes = 0;
509}