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