PWGDQ/dielectron/AliReducedEvent.cxx

   1 /*
   2 ***********************************************************
   3   Implementation of reduced ESD information classes for
   4   quick analysis.
   5   Contact: i.c.arsene@gsi.de, i.c.arsene@cern.ch
   6   2012/06/21
   7   *********************************************************
   8 */
   9
  10 #ifndef ALIREDUCEDEVENT_H
  11 #include "AliReducedEvent.h"
  12 #endif
  13
  14 #include <TMath.h>
  15 #include <TClonesArray.h>
  16
  17 ClassImp(AliReducedTrack)
  18 ClassImp(AliReducedPair)
  19 ClassImp(AliReducedEvent)
  20 ClassImp(AliReducedEventFriend)
  21 ClassImp(AliReducedCaloCluster)
  22
  23 TClonesArray* AliReducedEvent::fgTracks = 0;
  24 TClonesArray* AliReducedEvent::fgCandidates = 0;
  25 TClonesArray* AliReducedEvent::fgCaloClusters = 0;
  26
  27 //_______________________________________________________________________________
  28 AliReducedTrack::AliReducedTrack() :
  29   fTrackId(-1),
  30   fStatus(0),
  31   fGlobalPhi(0.0),
  32   fGlobalPt(0.0),
  33   fGlobalEta(0.0),
  34   fTPCPhi(0.0),
  35   fTPCPt(0.0),
  36   fTPCEta(0.0),
  37   fMomentumInner(0.0),
  38   fDCA(),
  39   fITSclusterMap(0),
  40   fITSsignal(0.0),
  41   fITSnSig(),
  42   fTPCNcls(0),
  43   fTPCCrossedRows(0),
  44   fTPCNclsF(0),
  45   fTPCNclsIter1(0),
  46   fTPCClusterMap(0),
  47   fTPCsignal(0),
  48   fTPCnSig(),
  49   fTOFbeta(0.0),
  50   fTOFnSig(),
  51   fTRDntracklets(),
  52   fTRDpid(),
  53   fCaloClusterId(-999),
  54   fBayesPID(),
  55   fFlags(0),
  56   fMoreFlags(0)
  57 {
  58   //
  59   // Constructor
  60   //
  61   fDCA[0] = 0.0; fDCA[1]=0.0;
  62   for(Int_t i=0; i<4; ++i) {fTPCnSig[i]=-999.; fTOFnSig[i]=-999.; fITSnSig[i]=-999.;}
  63   for(Int_t i=0; i<3; ++i) {fBayesPID[i]=-999.;}
  64   fTRDpid[0]=-999.; fTRDpid[1]=-999.;
  65 }
  66
  67
  68 //_______________________________________________________________________________
  69 AliReducedTrack::~AliReducedTrack()
  70 {
  71   //
  72   // De-Constructor
  73   //
  74 }
  75
  76
  77 //_______________________________________________________________________________
  78 AliReducedPair::AliReducedPair() :
  79   fCandidateId(-1),
  80   fPairType(-1),
  81   fLegIds(),
  82   fMass(),
  83   fPhi(0.0),
  84   fPt(0.0),
  85   fEta(0.0),
  86   fLxy(0.0),
  87   fLxyErr(0.0),
  88   fPointingAngle(0.0),
  89   fMCid(0)
  90 {
  91   //
  92   // Constructor
  93   //
  94   fLegIds[0] = -1; fLegIds[1] = -1;
  95   fMass[0]=-999.; fMass[1]=-999.; fMass[2]=-999.; fMass[3]=-999.;
  96 }
  97
  98
  99 //_______________________________________________________________________________
 100 AliReducedPair::AliReducedPair(const AliReducedPair &c) :
 101   TObject(c),
 102   fCandidateId(c.CandidateId()),
 103   fPairType(c.PairType()),
 104   fLegIds(),
 105   fMass(),
 106   fPhi(c.Phi()),
 107   fPt(c.Pt()),
 108   fEta(c.Eta()),
 109   fLxy(c.Lxy()),
 110   fLxyErr(c.LxyErr()),
 111   fPointingAngle(c.PointingAngle()),
 112   fMCid(c.MCid())
 113 {
 114   //
 115   // copy constructor
 116   //
 117   fLegIds[0] = c.LegId(0);
 118   fLegIds[1] = c.LegId(1);
 119   fMass[0] = c.Mass(0); fMass[1] = c.Mass(1); fMass[2] = c.Mass(2); fMass[3] = c.Mass(3);
 120 }
 121
 122
 123 //_______________________________________________________________________________
 124 AliReducedPair::~AliReducedPair() {
 125   //
 126   // destructor
 127   //
 128 }
 129
 130 //____________________________________________________________________________
 131 AliReducedEvent::AliReducedEvent() :
 132   TObject(),
 133   fRunNo(0),
 134   fBC(0),
 135   fTriggerMask(0),
 136   fIsPhysicsSelection(kTRUE),
 137   fVtx(),
 138   fNVtxContributors(0),
 139   fVtxTPC(),
 140   fNVtxTPCContributors(0),
 141   fCentrality(),
 142   fCentQuality(0),
 143   fNV0candidates(),
 144   fNDielectronCandidates(0),
 145   fNtracks(),
 146   fSPDntracklets(0),
 147   fVZEROMult(),
 148   fZDCnEnergy(),
 149   fTracks(0x0),
 150   fCandidates(0x0),
 151   fNCaloClusters(0),
 152   fCaloClusters(0x0)
 153 //fFMDMult()
 154 {
 155   //
 156   // Constructor
 157   //
 158   for(Int_t i=0; i<3; ++i) {fVtx[i]=-999.; fVtxTPC[i]=-999.;}
 159   for(Int_t i=0; i<4; ++i) fCentrality[i]=-1.;
 160   fNV0candidates[0]=0; fNV0candidates[1]=0;
 161   fNtracks[0]=0; fNtracks[1]=0;
 162   for(Int_t i=0; i<64; ++i) fVZEROMult[i] = 0.0;
 163   for(Int_t i=0; i<8; ++i) fZDCnEnergy[i]=0.0;
 164 }
 165
 166
 167 //____________________________________________________________________________
 168 AliReducedEvent::AliReducedEvent(const Char_t* /*name*/) :
 169   TObject(),
 170   fRunNo(0),
 171   fBC(0),
 172   fTriggerMask(0),
 173   fIsPhysicsSelection(kTRUE),
 174   fVtx(),
 175   fNVtxContributors(0),
 176   fVtxTPC(),
 177   fNVtxTPCContributors(0),
 178   fCentrality(),
 179   fCentQuality(0),
 180   fNV0candidates(),
 181   fNDielectronCandidates(0),
 182   fNtracks(),
 183   fSPDntracklets(0),
 184   fVZEROMult(),
 185   fZDCnEnergy(),
 186   fTracks(0x0),
 187   fCandidates(0x0),
 188   fNCaloClusters(0),
 189   fCaloClusters(0x0)
 190 //fFMDMult()
 191 {
 192   //
 193   // Constructor
 194   //
 195   for(Int_t i=0; i<3; ++i) {fVtx[i]=-999.; fVtxTPC[i]=-999.;}
 196   for(Int_t i=0; i<4; ++i) fCentrality[i]=-1.;
 197   fNV0candidates[0]=0; fNV0candidates[1]=0;
 198   fNtracks[0]=0; fNtracks[1]=0;
 199   for(Int_t i=0; i<64; ++i) fVZEROMult[i] = 0.0;
 200   for(Int_t i=0; i<8; ++i) fZDCnEnergy[i]=0.0;
 201
 202   if(!fgTracks) fgTracks = new TClonesArray("AliReducedTrack", 100000);
 203   fTracks = fgTracks;
 204   if(!fgCandidates) fgCandidates = new TClonesArray("AliReducedPair", 100000);
 205   fCandidates = fgCandidates;
 206   if(!fgCaloClusters) fgCaloClusters = new TClonesArray("AliReducedCaloCluster", 50000);
 207   fCaloClusters = fgCaloClusters;
 208 }
 209
 210
 211 //____________________________________________________________________________
 212 AliReducedEvent::~AliReducedEvent()
 213 {
 214   //
 215   // De-Constructor
 216   //
 217   //ClearEvent();
 218 }
 219
 220
 221 //____________________________________________________________________________
 222 Float_t AliReducedEvent::MultVZEROA() const
 223 {
 224   //
 225   // Total VZERO multiplicity in A side
 226   //
 227   Float_t mult=0.0;
 228   for(Int_t i=32;i<64;++i)
 229     mult += fVZEROMult[i];
 230   return mult;
 231 }
 232
 233
 234 //____________________________________________________________________________
 235 Float_t AliReducedEvent::MultVZEROC() const
 236 {
 237   //
 238   // Total VZERO multiplicity in C side
 239   //
 240   Float_t mult=0.0;
 241   for(Int_t i=0;i<32;++i)
 242     mult += fVZEROMult[i];
 243   return mult;
 244 }
 245
 246
 247 //____________________________________________________________________________
 248 Float_t AliReducedEvent::MultVZERO() const
 249 {
 250   //
 251   // Total VZERO multiplicity
 252   //
 253   return MultVZEROA()+MultVZEROC();
 254 }
 255
 256
 257 //____________________________________________________________________________
 258 Float_t AliReducedEvent::MultRingVZEROA(Int_t ring) const
 259 {
 260   //
 261   //  VZERO multiplicity in a given ring on A side
 262   //
 263   if(ring<0 || ring>3) return -1.0;
 264
 265   Float_t mult=0.0;
 266   for(Int_t i=32+8*ring; i<32+8*(ring+1); ++i)
 267     mult += fVZEROMult[i];
 268   return mult;
 269 }
 270
 271
 272 //____________________________________________________________________________
 273 Float_t AliReducedEvent::MultRingVZEROC(Int_t ring) const
 274 {
 275   //
 276   //  VZERO multiplicity in a given ring on C side
 277   //
 278   if(ring<0 || ring>3) return -1.0;
 279
 280   Float_t mult=0.0;
 281   for(Int_t i=8*ring; i<8*(ring+1); ++i)
 282     mult += fVZEROMult[i];
 283   return mult;
 284 }
 285
 286 //_____________________________________________________________________________
 287 void AliReducedEvent::ClearEvent() {
 288   //
 289   // clear the event
 290   //
 291   if(fTracks) fTracks->Clear("C");
 292   if(fCandidates) fCandidates->Clear("C");
 293   if(fCaloClusters) fCaloClusters->Clear("C");
 294   fRunNo = 0;
 295   fBC = 0;
 296   fTriggerMask = 0;
 297   fIsPhysicsSelection = kTRUE;
 298   fCentQuality = 0;
 299   fNV0candidates[0] = 0; fNV0candidates[1] = 0;
 300   fNDielectronCandidates = 0;
 301   fNtracks[0] = 0; fNtracks[1] = 0;
 302   for(Int_t i=0; i<3; ++i) {fVtx[i]=-999.; fVtxTPC[i]=-999.; fCentrality[i]=-1.;}
 303   for(Int_t i=0; i<64; ++i) fVZEROMult[i] = 0.0;
 304   for(Int_t i=0; i<8; ++i) fZDCnEnergy[i]=0.0;
 305 }
 306
 307
 308 //_______________________________________________________________________________
 309 AliReducedEventFriend::AliReducedEventFriend() :
 310  fQvector(),
 311  fEventPlaneStatus()
 312 {
 313   //
 314   // default constructor
 315   //
 316   for(Int_t idet=0; idet<kNdetectors; ++idet) {
 317     for(Int_t ih=0; ih<fgkNMaxHarmonics; ++ih) {
 318       fEventPlaneStatus[idet][ih] = 0;
 319       for(Int_t ic=0; ic<2; ++ic)
 320         fQvector[idet][ih][ic] = 0.0;
 321     }
 322   }
 323 }
 324
 325
 326 //____________________________________________________________________________
 327 AliReducedEventFriend::~AliReducedEventFriend()
 328 {
 329   //
 330   // De-Constructor
 331   //
 332   ClearEvent();
 333 }
 334
 335
 336 //_____________________________________________________________________________
 337 void AliReducedEventFriend::ClearEvent() {
 338   //
 339   // clear the event
 340   //
 341   for(Int_t idet=0; idet<kNdetectors; ++idet) {
 342     for(Int_t ih=0; ih<fgkNMaxHarmonics; ++ih) {
 343       fEventPlaneStatus[idet][ih] = 0;
 344       for(Int_t ic=0; ic<2; ++ic)
 345         fQvector[idet][ih][ic] = 0.0;
 346     }
 347   }
 348 }
 349
 350
 351 //____________________________________________________________________________
 352 void AliReducedEventFriend::CopyEvent(const AliReducedEventFriend* event) {
 353   //
 354   // copy information from another event to this one
 355   //
 356   for(Int_t idet=0; idet<kNdetectors; ++idet) {
 357     for(Int_t ih=0; ih<fgkNMaxHarmonics; ++ih) {
 358       fQvector[idet][ih][0] = event->Qx(idet, ih+1);
 359       fQvector[idet][ih][1] = event->Qy(idet, ih+1);
 360       fEventPlaneStatus[idet][ih] = event->GetEventPlaneStatus(idet, ih+1);
 361     }
 362   }
 363 }
 364
 365
 366 //_____________________________________________________________________________
 367 AliReducedCaloCluster::AliReducedCaloCluster() :
 368  fType(kUndefined),
 369  fEnergy(-999.),
 370  fTrackDx(-999.),
 371  fTrackDz(-999.),
 372  fM20(-999.),
 373  fM02(-999.),
 374  fDispersion(-999.)
 375 {
 376   //
 377   // default constructor
 378   //
 379 }
 380
 381
 382 //_____________________________________________________________________________
 383 AliReducedCaloCluster::~AliReducedCaloCluster()
 384 {
 385   //
 386   // destructor
 387   //
 388 }
 389
 390
 391 //_______________________________________________________________________________
 392 void AliReducedEvent::GetQvector(Double_t Qvec[][2], Int_t det,
 393                                  Float_t etaMin/*=-0.8*/, Float_t etaMax/*=+0.8*/,
 394                                  Bool_t (*IsTrackSelected)(AliReducedTrack*)/*=NULL*/) {
 395   //
 396   // Get the event plane for a specified detector
 397   //
 398   if(det==AliReducedEventFriend::kTPC ||
 399      det==AliReducedEventFriend::kTPCptWeights ||
 400      det==AliReducedEventFriend::kTPCpos ||
 401      det==AliReducedEventFriend::kTPCneg) {
 402     GetTPCQvector(Qvec, det, etaMin, etaMax, IsTrackSelected);
 403     return;
 404   }
 405   if(det==AliReducedEventFriend::kVZEROA ||
 406      det==AliReducedEventFriend::kVZEROC) {
 407     GetVZEROQvector(Qvec, det);
 408     return;
 409   }
 410   if(det==AliReducedEventFriend::kZDCA ||
 411      det==AliReducedEventFriend::kZDCC) {
 412     GetZDCQvector(Qvec, det);
 413     return;
 414   }
 415   if(det==AliReducedEventFriend::kFMD) {
 416     //TODO implementation
 417     return;
 418   }
 419   return;
 420 }
 421
 422
 423 //_________________________________________________________________
 424 Int_t AliReducedEvent::GetTPCQvector(Double_t Qvec[][2], Int_t det,
 425                                      Float_t etaMin/*=-0.8*/, Float_t etaMax/*=+0.8*/,
 426                                      Bool_t (*IsTrackSelected)(AliReducedTrack*)/*=NULL*/) {
 427   //
 428   // Construct the event plane using tracks in the barrel
 429   //
 430   if(!(det==AliReducedEventFriend::kTPC ||
 431        det==AliReducedEventFriend::kTPCpos ||
 432        det==AliReducedEventFriend::kTPCneg))
 433     return 0;
 434   Int_t nUsedTracks = 0;
 435   Short_t charge = 0;
 436   AliReducedTrack* track=0x0;
 437   Double_t weight=0.0; Double_t absWeight = 0.0; Double_t x=0.0; Double_t y=0.0;
 438   TIter nextTrack(fTracks);
 439   while((track=static_cast<AliReducedTrack*>(nextTrack()))) {
 440     if(track->Eta()<etaMin) continue;
 441     if(track->Eta()>etaMax) continue;
 442     charge = track->Charge();
 443     if(det==AliReducedEventFriend::kTPCpos && charge<0) continue;
 444     if(det==AliReducedEventFriend::kTPCneg && charge>0) continue;
 445
 446     if(IsTrackSelected && !IsTrackSelected(track)) continue;
 447     absWeight = 1.0;
 448     if(det==AliReducedEventFriend::kTPCptWeights) {
 449       absWeight = track->Pt();
 450       if(absWeight>2.0) absWeight = 2.0;    // pt is the weight used for the event plane
 451     }
 452     weight = absWeight;
 453
 454     ++nUsedTracks;
 455     x = TMath::Cos(track->Phi());
 456     y = TMath::Sin(track->Phi());
 457     //  1st harmonic
 458     Qvec[0][0] += weight*x;
 459     Qvec[0][1] += weight*y;
 460     //  2nd harmonic
 461     Qvec[1][0] += absWeight*(2.0*TMath::Power(x,2.0)-1);
 462     Qvec[1][1] += absWeight*(2.0*x*y);
 463     //  3rd harmonic
 464     Qvec[2][0] += weight*(4.0*TMath::Power(x,3.0)-3.0*x);
 465     Qvec[2][1] += weight*(3.0*y-4.0*TMath::Power(y,3.0));
 466     //  4th harmonic
 467     Qvec[3][0] += absWeight*(1.0-8.0*TMath::Power(x*y,2.0));
 468     Qvec[3][1] += absWeight*(4.0*x*y-8.0*x*TMath::Power(y,3.0));
 469     //  5th harmonic
 470     Qvec[4][0] += weight*(16.0*TMath::Power(x,5.0)-20.0*TMath::Power(x, 3.0)+5.0*x);
 471     Qvec[4][1] += weight*(16.0*TMath::Power(y,5.0)-20.0*TMath::Power(y, 3.0)+5.0*y);
 472     //  6th harmonic
 473     Qvec[5][0] += absWeight*(32.0*TMath::Power(x,6.0)-48.0*TMath::Power(x, 4.0)+18.0*TMath::Power(x, 2.0)-1.0);
 474     Qvec[5][1] += absWeight*(x*y*(32.0*TMath::Power(y,4.0)-32.0*TMath::Power(y, 2.0)+6.0));
 475   }
 476   return nUsedTracks;
 477 }
 478
 479
 480 //____________________________________________________________________________
 481 void AliReducedEvent::SubtractParticleFromQvector(
 482         AliReducedTrack* particle, Double_t Qvec[][2], Int_t det,
 483         Float_t etaMin/*=-0.8*/, Float_t etaMax/*=+0.8*/,
 484         Bool_t (*IsTrackSelected)(AliReducedTrack*)/*=NULL*/) {
 485   //
 486   // subtract a particle from the event Q-vector
 487   //
 488   Float_t eta = particle->Eta();
 489   if(eta<etaMin) return;
 490   if(eta>etaMax) return;
 491
 492   Float_t charge = particle->Charge();
 493   if(det==AliReducedEventFriend::kTPCpos && charge<0) return;
 494   if(det==AliReducedEventFriend::kTPCneg && charge>0) return;
 495
 496   if(IsTrackSelected && !IsTrackSelected(particle)) return;
 497
 498   Double_t weight=0.0; Double_t absWeight = 0.0;
 499   if(det==AliReducedEventFriend::kTPCptWeights) {
 500     absWeight = particle->Pt();
 501     if(absWeight>2.0) absWeight = 2.0;
 502   }
 503   weight = absWeight;
 504   //  if(eta<0.0) weight *= -1.0;
 505
 506   Float_t x = TMath::Cos(particle->Phi());
 507   Float_t y = TMath::Sin(particle->Phi());
 508
 509   //  1st harmonic
 510   Qvec[0][0] -= weight*x;
 511   Qvec[0][1] -= weight*y;
 512   //  2nd harmonic
 513   Qvec[1][0] -= absWeight*(2.0*TMath::Power(x,2.0)-1);
 514   Qvec[1][1] -= absWeight*(2.0*x*y);
 515   //  3rd harmonic
 516   Qvec[2][0] -= weight*(4.0*TMath::Power(x,3.0)-3.0*x);
 517   Qvec[2][1] -= weight*(3.0*y-4.0*TMath::Power(y,3.0));
 518   //  4th harmonic
 519   Qvec[3][0] -= absWeight*(1.0-8.0*TMath::Power(x*y,2.0));
 520   Qvec[3][1] -= absWeight*(4.0*x*y-8.0*x*TMath::Power(y,3.0));
 521   //  5th harmonic
 522   Qvec[4][0] -= weight*(16.0*TMath::Power(x,5.0)-20.0*TMath::Power(x, 3.0)+5.0*x);
 523   Qvec[4][1] -= weight*(16.0*TMath::Power(y,5.0)-20.0*TMath::Power(y, 3.0)+5.0*y);
 524   //  6th harmonic
 525   Qvec[5][0] -= absWeight*(32.0*TMath::Power(x,6.0)-48.0*TMath::Power(x, 4.0)+18.0*TMath::Power(x, 2.0)-1.0);
 526   Qvec[5][1] -= absWeight*(x*y*(32.0*TMath::Power(y,4.0)-32.0*TMath::Power(y, 2.0)+6.0));
 527
 528   return;
 529 }
 530
 531
 532 //_________________________________________________________________
 533 void AliReducedEvent::GetVZEROQvector(Double_t Qvec[][2], Int_t det) {
 534   //
 535   // Get the reaction plane from the VZERO detector for a given harmonic
 536   //
 537   GetVZEROQvector(Qvec, det, fVZEROMult);
 538 }
 539
 540
 541 //_________________________________________________________________
 542 void AliReducedEvent::GetVZEROQvector(Double_t Qvec[][2], Int_t det, Float_t* vzeroMult) {
 543   //
 544   // Get the reaction plane from the VZERO detector for a given harmonic
 545   //
 546   //  Q{x,y} = SUM_i mult(i) * {cos(n*phi_i), sin(n*phi_i)}
 547   //  phi_i - phi angle of the VZERO sector i
 548   //          Each sector covers 45 degrees(8 sectors per ring). Middle of sector 0 is at 45/2
 549   //        channel 0: 22.5
 550   //                1: 22.5+45
 551   //                2: 22.5+45*2
 552   //               ...
 553   //        at the next ring continues the same
 554   //        channel 8: 22.5
 555   //        channel 9: 22.5 + 45
 556   //
 557   if(!(det==AliReducedEventFriend::kVZEROA ||
 558        det==AliReducedEventFriend::kVZEROC))
 559     return;
 560
 561   const Double_t kX[8] = {0.92388, 0.38268, -0.38268, -0.92388, -0.92388, -0.38268, 0.38268, 0.92388};    // cosines of the angles of the VZERO sectors (8 per ring)
 562   const Double_t kY[8] = {0.38268, 0.92388, 0.92388, 0.38268, -0.38268, -0.92388, -0.92388, -0.38268};    // sines     -- " --
 563   Int_t phi;
 564
 565   for(Int_t iChannel=0; iChannel<64; ++iChannel) {
 566     if(iChannel<32 && det==AliReducedEventFriend::kVZEROA) continue;
 567     if(iChannel>=32 && det==AliReducedEventFriend::kVZEROC) continue;
 568     phi=iChannel%8;
 569     //  1st harmonic
 570     Qvec[0][0] += vzeroMult[iChannel]*kX[phi];
 571     Qvec[0][1] += vzeroMult[iChannel]*kY[phi];
 572     //  2nd harmonic
 573     Qvec[1][0] += vzeroMult[iChannel]*(2.0*TMath::Power(kX[phi],2.0)-1);
 574     Qvec[1][1] += vzeroMult[iChannel]*(2.0*kX[phi]*kY[phi]);
 575     //  3rd harmonic
 576     Qvec[2][0] += vzeroMult[iChannel]*(4.0*TMath::Power(kX[phi],3.0)-3.0*kX[phi]);
 577     Qvec[2][1] += vzeroMult[iChannel]*(3.0*kY[phi]-4.0*TMath::Power(kY[phi],3.0));
 578     //  4th harmonic
 579     Qvec[3][0] += vzeroMult[iChannel]*(1.0-8.0*TMath::Power(kX[phi]*kY[phi],2.0));
 580     Qvec[3][1] += vzeroMult[iChannel]*(4.0*kX[phi]*kY[phi]-8.0*kX[phi]*TMath::Power(kY[phi],3.0));
 581     //  5th harmonic
 582     Qvec[4][0] += vzeroMult[iChannel]*(16.0*TMath::Power(kX[phi],5.0)-20.0*TMath::Power(kX[phi], 3.0)+5.0*kX[phi]);
 583     Qvec[4][1] += vzeroMult[iChannel]*(16.0*TMath::Power(kY[phi],5.0)-20.0*TMath::Power(kY[phi], 3.0)+5.0*kY[phi]);
 584     //  6th harmonic
 585     Qvec[5][0] += vzeroMult[iChannel]*(32.0*TMath::Power(kX[phi],6.0)-48.0*TMath::Power(kX[phi], 4.0)+18.0*TMath::Power(kX[phi], 2.0)-1.0);
 586     Qvec[5][1] += vzeroMult[iChannel]*(kX[phi]*kY[phi]*(32.0*TMath::Power(kY[phi],4.0)-32.0*TMath::Power(kY[phi], 2.0)+6.0));
 587   }    // end loop over channels
 588 }
 589
 590
 591 //_________________________________________________________________
 592 void AliReducedEvent::GetZDCQvector(Double_t Qvec[][2], Int_t det) const {
 593   //
 594   // Construct the event plane using the ZDC
 595   // ZDC has 2 side (A and C) with 4 calorimeters on each side
 596   // The XY position of each calorimeter is specified by the
 597   // zdcNTowerCenters_x and zdcNTowerCenters_y arrays
 598   if(!(det==AliReducedEventFriend::kZDCA ||
 599        det==AliReducedEventFriend::kZDCC)) return;       // bad detector option
 600   const Float_t zdcTowerCenter = 1.75;
 601   const Float_t zdcNTowerCentersX[4] = {-zdcTowerCenter,  zdcTowerCenter, -zdcTowerCenter, zdcTowerCenter};
 602   const Float_t zdcNTowerCentersY[4] = {-zdcTowerCenter, -zdcTowerCenter,  zdcTowerCenter, zdcTowerCenter};
 603
 604   Qvec[0][0] = 0.0; Qvec[0][1] = 0.0;   // first harmonic Q-vector
 605   Float_t zdcNCentroidSum = 0;
 606   Float_t zdcNalpha = 0.5;
 607
 608   for(Int_t i=0; i<4; ++i) {
 609     if(fZDCnEnergy[i+(det==AliReducedEventFriend::kZDCA ? 4 : 0)]>0.0) {
 610       Float_t zdcNenergyAlpha = TMath::Power(fZDCnEnergy[i+(det==AliReducedEventFriend::kZDCA ? 4 : 0)], zdcNalpha);
 611       Qvec[0][0] += zdcNTowerCentersX[i-1]*zdcNenergyAlpha;
 612       Qvec[0][1] += zdcNTowerCentersY[i-1]*zdcNenergyAlpha;
 613       zdcNCentroidSum += zdcNenergyAlpha;
 614     }
 615   }   // end loop over channels
 616
 617   if(zdcNCentroidSum>0.0) {
 618     Qvec[0][0] /= zdcNCentroidSum;
 619     Qvec[0][1] /= zdcNCentroidSum;
 620   }
 621 }