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 < 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
  61 AliEventplane::AliEventplane(const AliEventplane& ep) :
  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 {
  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     TNamed::operator=(ep);
  84     ((AliEventplane &) ep).CopyEP(*this);
  85   }
  86   return *this;
  87 }
  88
  89 void AliEventplane::CopyEP(AliEventplane& ep) const
  90 { // copy function
  91
  92   AliEventplane& target = (AliEventplane &) ep;
  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
 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
 221 AliEventplane::~AliEventplane()
 222 {
 223   /// destructor
 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;
 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   }
 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
 262 TVector2* AliEventplane::GetQVector()
 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
 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   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
 444 TVector2* AliEventplane::GetQsub1()
 445 {
 446   return fQsub1;
 447 }
 448
 449 TVector2* AliEventplane::GetQsub2()
 450 {
 451   return fQsub2;
 452 }
 453
 454 Double_t AliEventplane::GetQsubRes()
 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
 466 Double_t AliEventplane::GetQContributionX(AliVTrack* track)
 467 {
 468   return fQContributionX->GetAt(track->GetID());
 469 }
 470
 471 Double_t AliEventplane::GetQContributionY(AliVTrack* track)
 472 {
 473   return fQContributionY->GetAt(track->GetID());
 474 }
 475
 476 Double_t AliEventplane::GetQContributionXsub1(AliVTrack* track)
 477 {
 478   return fQContributionXsub1->GetAt(track->GetID());
 479 }
 480
 481 Double_t AliEventplane::GetQContributionYsub1(AliVTrack* track)
 482 {
 483   return fQContributionYsub1->GetAt(track->GetID());
 484 }
 485
 486 Double_t AliEventplane::GetQContributionXsub2(AliVTrack* track)
 487 {
 488   return fQContributionXsub2->GetAt(track->GetID());
 489 }
 490
 491 Double_t AliEventplane::GetQContributionYsub2(AliVTrack* track)
 492 {
 493   return fQContributionYsub2->GetAt(track->GetID());
 494 }
 495
 496 void AliEventplane::Reset()
 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 }