AliRoot Core  3dc7879 (3dc7879)
AliEventplane.cxx
Go to the documentation of this file.
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 //*****************************************************
21 
22 #include "AliLog.h"
23 #include "AliEventplane.h"
24 #include "TVector2.h"
25 #include "AliVTrack.h"
26 #include "TObjArray.h"
27 #include "TArrayF.h"
28 #include "AliVEvent.h"
29 #include "AliVVZERO.h"
30 
31 ClassImp(AliEventplane)
32 
33 AliEventplane::AliEventplane() : TNamed("Eventplane", "Eventplane"),
34  fQVector(0),
35  fQContributionX(0),
36  fQContributionY(0),
37  fQContributionXsub1(0),
38  fQContributionYsub1(0),
39  fQContributionXsub2(0),
40  fQContributionYsub2(0),
41  fEventplaneQ(-1),
42  fQsub1(0),
43  fQsub2(0),
44  fQsubRes(0)
45 {
47  fQContributionX = new TArrayF(0);
48  fQContributionY = new TArrayF(0);
49  fQContributionXsub1 = new TArrayF(0);
50  fQContributionYsub1 = new TArrayF(0);
51  fQContributionXsub2 = new TArrayF(0);
52  fQContributionYsub2 = new TArrayF(0);
53  for(Int_t i = 0; i < 11; ++i) {
54  fMeanX2[i]=fMeanY2[i]=0.;
55  fAPlus[i]=fAMinus[i]=1.;
56  fLambdaPlus[i]=fLambdaMinus[i]=0.;
57  fCos8Psi[i]=0.;
58  }
59 }
60 
62  TNamed(ep),
63  fQVector(0),
64  fQContributionX(0),
65  fQContributionY(0),
66  fQContributionXsub1(0),
67  fQContributionYsub1(0),
68  fQContributionXsub2(0),
69  fQContributionYsub2(0),
70  fEventplaneQ(0),
71  fQsub1(0),
72  fQsub2(0),
73  fQsubRes(0)
74 {
76  ((AliEventplane &) ep).CopyEP(*this);
77 }
78 
80 {
82  if (this!=&ep){
83  TNamed::operator=(ep);
84  ((AliEventplane &) ep).CopyEP(*this);
85  }
86  return *this;
87 }
88 
90 { // copy function
91 
92  AliEventplane& target = (AliEventplane &) ep;
93 
94  if(fQContributionX)
95  {
96  if(ep.fQContributionX)
97  {
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  {
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 
130  {
131  if(ep.fQContributionXsub1)
132  {
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 
147  {
148  if(ep.fQContributionYsub1)
149  {
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 
165  {
166  if(ep.fQContributionXsub2)
167  {
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 
182  {
183  if(ep.fQContributionYsub2)
184  {
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 
198  if (fEventplaneQ)
199  target.fEventplaneQ = fEventplaneQ;
200 
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;
209 
210  for(Int_t i = 0; i < 11; ++i) {
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  }
219 }
220 
222 {
224  if (fQContributionX){
225  delete fQContributionX;
226  fQContributionX = 0;
227  }
228  if (fQContributionY){
229  delete fQContributionY;
230  fQContributionY = 0;
231  }
232  if (fQContributionXsub1){
233  delete fQContributionXsub1;
235  }
236  if (fQContributionYsub1){
237  delete fQContributionYsub1;
239  }
240  if (fQContributionXsub2){
241  delete fQContributionXsub2;
243  }
244  if (fQContributionYsub2){
245  delete fQContributionYsub2;
247  }
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 
263 {
264  return fQVector;
265 }
266 
267 Double_t AliEventplane::GetEventplane(const char *x, const AliVEvent *event, Int_t harmonic) const
268 {
269  TString method = x;
270  Double_t qx = 0, qy = 0;
271  if(method.CompareTo("Q")==0) return fEventplaneQ;
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);
275 
276  return -1000.;
277 }
278 
279 Double_t AliEventplane::CalculateVZEROEventPlane(const AliVEvent * event, Int_t firstRing, Int_t lastRing, Int_t harmonic, Double_t &qxTot, Double_t &qyTot) const
280 {
281  // Calculate the VZERO event-plane by summing rings
282  // The flatenning is done on every ring separately
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 
297  qxTot=qyTot=0.;
298  Double_t qx, qy;
299  Double_t totMult = 0.;
300  for(Int_t ring = firstRing; ring <=lastRing; ++ring) {
301  qx=qy=0.;
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
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]);
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];
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;
336  }
337  else {
338  qxTot += qx;
339  qyTot += qy;
340  }
341  }
342 
343  // In case of no hits return an invalid event-plane angle
344  if (totMult<1e-6) return -999;
345 
346  return (TMath::ATan2(qyTot,qxTot)/harmonic);
347 }
348 
349 Double_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 
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  }
385  qx=qy=0.;
386  Double_t multRing = 0.;
387  for(Int_t iCh = firstCh; iCh < lastCh; ++iCh) {
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 
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
433  if (ring < 0 || ring >= 11) AliFatal(Form("Bad ring index (%d)",ring));
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 }
443 
445 {
446  return fQsub1;
447 }
448 
450 {
451  return fQsub2;
452 }
453 
455 {
456  return fQsubRes;
457 }
458 
459 Bool_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 
467 {
468  return fQContributionX->GetAt(track->GetID());
469 }
470 
472 {
473  return fQContributionY->GetAt(track->GetID());
474 }
475 
477 {
478  return fQContributionXsub1->GetAt(track->GetID());
479 }
480 
482 {
483  return fQContributionYsub1->GetAt(track->GetID());
484 }
485 
487 {
488  return fQContributionXsub2->GetAt(track->GetID());
489 }
490 
492 {
493  return fQContributionYsub2->GetAt(track->GetID());
494 }
495 
497 {
498  delete fQVector; fQVector=0;
499  fQContributionX->Reset();
500  fQContributionY->Reset();
501  fQContributionXsub1->Reset();
502  fQContributionYsub1->Reset();
503  fQContributionXsub2->Reset();
504  fQContributionYsub2->Reset();
505  fEventplaneQ = -1;
506  delete fQsub1; fQsub1=0;
507  delete fQsub2; fQsub2=0;
508  fQsubRes = 0;
509 }
TVector2 * GetQsub2()
TBrowser b
Definition: RunAnaESD.C:12
Double_t fMeanX2[11]
Definition: AliEventplane.h:77
Double_t GetQContributionYsub1(AliVTrack *track)
TArrayF * fQContributionYsub2
Definition: AliEventplane.h:72
AliEventplane()
A container for the event plane stored in AOD in ESD.
AliEventplane & operator=(const AliEventplane &ep)
copy constructor
TVector2 * fQsub1
Definition: AliEventplane.h:74
void SetVZEROEPParams(Int_t ring, Double_t meanX2, Double_t meanY2, Double_t aPlus, Double_t aMinus, Double_t lambdaPlus, Double_t lambdaMinus, Double_t cos8Psi)
TArrayF * fQContributionXsub2
Definition: AliEventplane.h:71
Double_t fLambdaPlus[11]
Definition: AliEventplane.h:81
Bool_t IsEventInEventplaneClass(Double_t a, Double_t b, const char *method)
Double_t fEventplaneQ
Definition: AliEventplane.h:73
Double_t GetQContributionX(AliVTrack *track)
TArrayF * fQContributionX
Definition: AliEventplane.h:67
AliTPCfastTrack * track
Double_t CalculateVZEROEventPlane(const AliVEvent *event, Int_t firstRing, Int_t lastRing, Int_t harmonic, Double_t &qxTot, Double_t &qyTot) const
Double_t fCos8Psi[11]
Definition: AliEventplane.h:83
bool trans(const AliFMDIndex &x, const AliFMDIndex &y, const AliFMDIndex &z)
Definition: TestIndex.C:94
virtual void CopyEP(AliEventplane &ep) const
assignment operator
Double_t fLambdaMinus[11]
Definition: AliEventplane.h:82
TVector2 * GetQVector()
get event plane result
Double_t GetQContributionXsub1(AliVTrack *track)
TArrayF * fQContributionXsub1
Definition: AliEventplane.h:69
Double_t fAMinus[11]
Definition: AliEventplane.h:80
Double_t fQsubRes
Definition: AliEventplane.h:76
virtual Int_t GetID() const =0
Double_t GetQContributionXsub2(AliVTrack *track)
Double_t GetQContributionY(AliVTrack *track)
~AliEventplane()
constructor
TVector2 * GetQsub1()
#define AliFatal(message)
Definition: AliLog.h:640
TArrayF * fQContributionY
Definition: AliEventplane.h:68
TVector2 * fQVector
Definition: AliEventplane.h:66
#define AliError(message)
Definition: AliLog.h:591
Double_t GetEventplane(const char *x, const AliVEvent *event=NULL, Int_t harmonic=2) const
Double_t fAPlus[11]
Definition: AliEventplane.h:79
Double_t GetQContributionYsub2(AliVTrack *track)
Double_t GetQsubRes()
TArrayF * fQContributionYsub1
Definition: AliEventplane.h:70
TVector2 * fQsub2
Definition: AliEventplane.h:75
Double_t fMeanY2[11]
Definition: AliEventplane.h:78