Introduction of VZERO event-plane selection task that can be used in order to flatten...
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliEventplane.cxx
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  
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 {
46   /// constructor
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 < 8; ++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
61 AliEventplane::AliEventplane(const AliEventplane& ep) : 
62   TNamed(),
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 {
75   /// Copy constructor
76   ((AliEventplane &) ep).CopyEP(*this);
77 }
78
79 AliEventplane& AliEventplane::operator=(const AliEventplane& ep)
80 {
81   /// Assignment operator
82   if (this!=&ep)
83     ((AliEventplane &) ep).CopyEP(*this);
84
85   return *this;
86 }
87
88 void AliEventplane::CopyEP(AliEventplane& ep) const
89 { // copy function
90
91   AliEventplane& target = (AliEventplane &) ep;
92   if (fQContributionX)
93       target.fQContributionX = fQContributionX;
94   if (fQContributionY)
95       target.fQContributionY = fQContributionY;
96   if (fQContributionXsub1)
97       target.fQContributionXsub1 = fQContributionXsub1;
98   if (fQContributionYsub1)
99       target.fQContributionYsub1 = fQContributionYsub1;
100   if (fQContributionXsub2)
101       target.fQContributionXsub2 = fQContributionXsub2;
102   if (fQContributionYsub2)
103       target.fQContributionYsub2 = fQContributionYsub2;
104   if (fEventplaneQ)
105       target.fEventplaneQ = fEventplaneQ;
106   if (fQVector)
107       target.fQVector = dynamic_cast<TVector2*> (fQVector->Clone());
108   if (fQsub1)
109       target.fQsub1 = dynamic_cast<TVector2*> (fQsub1->Clone());
110   if (fQsub2)
111       target.fQsub2 = dynamic_cast<TVector2*> (fQsub2->Clone());
112   if (fQsubRes)
113       target.fQsubRes = fQsubRes;
114
115   for(Int_t i = 0; i < 8; ++i) {
116     target.fMeanX2[i]=fMeanX2[i];
117     target.fMeanY2[i]=fMeanY2[i];
118     target.fAPlus[i]=fAPlus[i];
119     target.fAMinus[i]=fAMinus[i];
120     target.fLambdaPlus[i]=fLambdaPlus[i];
121     target.fLambdaMinus[i]=fLambdaMinus[i];
122     target.fCos8Psi[i]=fCos8Psi[i];
123   }
124 }
125
126 AliEventplane::~AliEventplane()
127 {
128   /// destructor
129   if (fQContributionX){
130       delete fQContributionX;
131       fQContributionX = 0;
132   }
133   if (fQContributionY){
134       delete fQContributionY;
135       fQContributionY = 0;
136   }
137   if (fQContributionXsub1){
138       delete fQContributionXsub1;
139       fQContributionXsub1 = 0;
140   }
141   if (fQContributionYsub1){
142       delete fQContributionYsub1;
143       fQContributionYsub1 = 0;
144   }
145   if (fQContributionXsub2){
146       delete fQContributionXsub2;
147       fQContributionXsub2 = 0;
148   }
149   if (fQContributionYsub2){
150       delete fQContributionYsub2;
151       fQContributionYsub2 = 0;
152   }
153   if (fQVector){
154       delete fQVector;
155       fQVector = 0;
156   }
157   if (fQsub1){
158       delete fQsub1;
159       fQsub1 = 0;
160   }
161     if (fQsub2){
162       delete fQsub2;
163       fQsub2 = 0;
164   }
165 }
166
167 TVector2* AliEventplane::GetQVector()
168 {
169   return fQVector;
170 }
171
172 Double_t AliEventplane::GetEventplane(const char *x, const AliVEvent *event, Int_t harmonic) const
173 {
174   TString method = x;
175   if(method.CompareTo("Q")==0)      return fEventplaneQ;
176   else if(method.CompareTo("V0A")==0) return CalculateVZEROEventPlane(event, 4, 7, harmonic);
177   else if(method.CompareTo("V0C")==0) return CalculateVZEROEventPlane(event, 0, 3, harmonic);
178   else if(method.CompareTo("V0")==0)  return CalculateVZEROEventPlane(event, 0, 7, harmonic);
179
180   return -1000.;
181 }
182
183 Double_t AliEventplane::CalculateVZEROEventPlane(const AliVEvent *  event, Int_t firstRing, Int_t lastRing, Int_t harmonic) const
184 {
185   if(!event) {
186     AliError("No Event received");
187     return -1000.;
188   }
189   AliVVZERO *vzeroData = event->GetVZEROData();
190   if(!vzeroData) {
191     AliError("Enable to get VZERO Data");
192     return -1000.;
193   }
194   if(harmonic <= 0) {
195     AliError("Required harmonic is less or equal to 0");
196     return -1000.;
197   }
198
199   Double_t qxTot=0., qyTot=0.;
200   Double_t qx, qy;
201   Double_t totMult = 0.;
202   for(Int_t ring = firstRing; ring <=lastRing; ++ring) {
203     qx=qy=0.;
204     Double_t trans[2][2];
205     trans[0][0] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
206     trans[0][1] = -fLambdaMinus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
207     trans[1][0] = -fLambdaPlus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
208     trans[1][1] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
209     Double_t multRing = 0.;
210     for(Int_t iCh = ring*8; iCh < (ring+1)*8; ++iCh) {
211       Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
212       Double_t mult = event->GetVZEROEqMultiplicity(iCh);
213       qx += mult*TMath::Cos(harmonic*phi);
214       qy += mult*TMath::Sin(harmonic*phi);
215       multRing += mult;
216     }
217
218     if (multRing < 1e-6) continue;
219     totMult += multRing;
220
221     if (harmonic == 2) {
222       // Recentering
223       Double_t qxPrime = qx - fMeanX2[ring];
224       Double_t qyPrime = qy - fMeanY2[ring];
225       // Twist
226       Double_t qxSeconde = trans[0][0]*qxPrime + trans[0][1]*qyPrime;
227       Double_t qySeconde = trans[1][0]*qxPrime + trans[1][1]*qyPrime;
228       // Rescaling
229       Double_t qxTierce = qxSeconde/fAPlus[ring];
230       Double_t qyTierce = qySeconde/fAMinus[ring];
231       qxTot += qxTierce;
232       qyTot += qyTierce;
233     }
234     else {
235       qxTot += qx;
236       qyTot += qy;
237     }
238   }
239
240   // In case of no hits return an invalid event-plane angle
241   if (totMult<1e-6) return -999;
242
243   // 4th Fourier momenta in order to flatten the EP within a sector
244   Double_t psi = TMath::ATan2(qyTot,qxTot)/harmonic;
245   if ((harmonic == 2) && (firstRing == lastRing)) {
246     psi += (fCos8Psi[firstRing]*TMath::Sin(2.*4.*psi))/2.;
247     if (psi > TMath::Pi()/2) psi -= TMath::Pi();
248     if (psi <-TMath::Pi()/2) psi += TMath::Pi();
249   }
250
251   return psi;
252 }
253
254 void AliEventplane::SetVZEROEPParams(Int_t ring,
255                                      Double_t meanX2, Double_t meanY2,
256                                      Double_t aPlus, Double_t aMinus,
257                                      Double_t lambdaPlus, Double_t lambdaMinus,
258                                      Double_t cos8Psi)
259 {
260   // Set the VZERO event-plane
261   // flatenning parameters
262   if (ring < 0 || ring >=8) AliFatal(Form("Bad ring index (%d)",ring));
263
264   fMeanX2[ring] = meanX2;
265   fMeanY2[ring] = meanY2;
266   fAPlus[ring] = aPlus;
267   fAMinus[ring] = aMinus;
268   fLambdaPlus[ring] = lambdaPlus;
269   fLambdaMinus[ring] = lambdaMinus;
270   fCos8Psi[ring] = cos8Psi;
271 }
272
273 TVector2* AliEventplane::GetQsub1()
274 {
275   return fQsub1;
276 }
277
278 TVector2* AliEventplane::GetQsub2()
279 {
280   return fQsub2;
281 }
282
283 Double_t AliEventplane::GetQsubRes()
284 {
285   return fQsubRes;
286 }
287
288 Bool_t AliEventplane::IsEventInEventplaneClass(Double_t a, Double_t b, const char *x)
289 {
290   TString method = x;
291   if ((method.CompareTo("Q")==0) && (fEventplaneQ >=a && fEventplaneQ < b)) return kTRUE;
292   else return kFALSE;
293 }
294
295 Double_t AliEventplane::GetQContributionX(AliVTrack* track)
296
297   return fQContributionX->GetAt(track->GetID());
298 }
299
300 Double_t AliEventplane::GetQContributionY(AliVTrack* track)
301
302   return fQContributionY->GetAt(track->GetID());
303 }
304
305 Double_t AliEventplane::GetQContributionXsub1(AliVTrack* track)
306
307   return fQContributionXsub1->GetAt(track->GetID());
308 }
309
310 Double_t AliEventplane::GetQContributionYsub1(AliVTrack* track)
311
312   return fQContributionYsub1->GetAt(track->GetID());
313 }
314
315 Double_t AliEventplane::GetQContributionXsub2(AliVTrack* track)
316
317   return fQContributionXsub2->GetAt(track->GetID());
318 }
319
320 Double_t AliEventplane::GetQContributionYsub2(AliVTrack* track)
321
322   return fQContributionYsub2->GetAt(track->GetID());
323 }
324
325 void AliEventplane::Reset()
326
327   delete fQVector; fQVector=0;
328   fQContributionX->Reset();
329   fQContributionY->Reset();
330   fQContributionXsub1->Reset();
331   fQContributionYsub1->Reset();
332   fQContributionXsub2->Reset();
333   fQContributionYsub2->Reset();
334   fEventplaneQ = -1;
335   delete fQsub1; fQsub1=0;
336   delete fQsub2; fQsub2=0;
337   fQsubRes = 0;
338 }