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   fTPCNcls(0),
  42   fTPCCrossedRows(0),
  43   fTPCNclsF(0),
  44   fTPCNclsIter1(0),
  45   fTPCClusterMap(0),
  46   fTPCsignal(0),
  47   fTPCnSig(),
  48   fTOFbeta(0.0),
  49   fTOFnSig(),
  50   fTRDntracklets(),
  51   fTRDpid(),
  52   fCaloClusterId(-999),
  53   fBayesPID(),
  54   fFlags(0)
  55 {
  56   //
  57   // Constructor
  58   //
  59   fDCA[0] = 0.0; fDCA[1]=0.0;
  60   for(Int_t i=0; i<4; ++i) {fTPCnSig[i]=-999.; fTOFnSig[i]=-999.;}
  61   for(Int_t i=0; i<3; ++i) {fBayesPID[i]=-999.;}
  62   fTRDpid[0]=-999.; fTRDpid[1]=-999.;
  63 }
  64
  65
  66 //_______________________________________________________________________________
  67 AliReducedTrack::~AliReducedTrack()
  68 {
  69   //
  70   // De-Constructor
  71   //
  72 }
  73
  74
  75 //_______________________________________________________________________________
  76 AliReducedPair::AliReducedPair() :
  77   fCandidateId(-1),
  78   fPairType(-1),
  79   fLegIds(),
  80   fMass(),
  81   fPhi(0.0),
  82   fPt(0.0),
  83   fEta(0.0),
  84   fLxy(0.0),
  85   fLxyErr(0.0),
  86   fOpeningAngle(-1.0),
  87   //fOnTheFly(kFALSE),
  88   fMCid(0)
  89 {
  90   //
  91   // Constructor
  92   //
  93   fLegIds[0] = -1; fLegIds[1] = -1;
  94   fMass[0]=-999.; fMass[1]=-999.; fMass[2]=-999.;
  95 }
  96
  97
  98 //_______________________________________________________________________________
  99 AliReducedPair::AliReducedPair(const AliReducedPair &c) :
 100   TObject(c),
 101   fCandidateId(c.CandidateId()),
 102   fPairType(c.PairType()),
 103   fLegIds(),
 104   fMass(),
 105   fPhi(c.Phi()),
 106   fPt(c.Pt()),
 107   fEta(c.Eta()),
 108   fLxy(c.Lxy()),
 109   fLxyErr(c.LxyErr()),
 110   fOpeningAngle(c.OpeningAngle()),
 111   //fOnTheFly(c.IsOnTheFly()),
 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);
 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 {
 373   //
 374   // default constructor
 375   //
 376 }
 377
 378
 379 //_____________________________________________________________________________
 380 AliReducedCaloCluster::~AliReducedCaloCluster()
 381 {
 382   //
 383   // destructor
 384   //
 385 }
 386
 387
 388 //_______________________________________________________________________________
 389 void AliReducedEvent::GetQvector(Double_t Qvec[][2], Int_t det,
 390                                             Float_t etaMin/*=-0.8*/, Float_t etaMax/*=+0.8*/,
 391                                             Bool_t (*IsTrackSelected)(AliReducedTrack*)/*=NULL*/) {
 392   //
 393   // Get the event plane for a specified detector
 394   //
 395   if(det==AliReducedEventFriend::kTPC ||
 396      det==AliReducedEventFriend::kTPCptWeights ||
 397      det==AliReducedEventFriend::kTPCpos ||
 398      det==AliReducedEventFriend::kTPCneg) {
 399     GetTPCQvector(Qvec, det, etaMin, etaMax, IsTrackSelected);
 400     return;
 401   }
 402   if(det==AliReducedEventFriend::kVZEROA ||
 403      det==AliReducedEventFriend::kVZEROC) {
 404     GetVZEROQvector(Qvec, det);
 405     return;
 406   }
 407   if(det==AliReducedEventFriend::kZDCA ||
 408      det==AliReducedEventFriend::kZDCC) {
 409     GetZDCQvector(Qvec, det);
 410     return;
 411   }
 412   if(det==AliReducedEventFriend::kFMD) {
 413     //TODO implementation
 414     return;
 415   }
 416   return;
 417 }
 418
 419
 420 //_________________________________________________________________
 421 Int_t AliReducedEvent::GetTPCQvector(Double_t Qvec[][2], Int_t det,
 422                                                 Float_t etaMin/*=-0.8*/, Float_t etaMax/*=+0.8*/,
 423                                                 Bool_t (*IsTrackSelected)(AliReducedTrack*)/*=NULL*/) {
 424   //
 425   // Construct the event plane using tracks in the barrel
 426   //
 427   if(!(det==AliReducedEventFriend::kTPC ||
 428        det==AliReducedEventFriend::kTPCpos ||
 429        det==AliReducedEventFriend::kTPCneg))
 430     return 0;
 431   Int_t nUsedTracks = 0;
 432   Short_t charge = 0;
 433   AliReducedTrack* track=0x0;
 434   Double_t weight=0.0; Double_t absWeight = 0.0; Double_t x=0.0; Double_t y=0.0;
 435   TIter nextTrack(fTracks);
 436   while((track=static_cast<AliReducedTrack*>(nextTrack()))) {
 437     if(track->Eta()<etaMin) continue;
 438     if(track->Eta()>etaMax) continue;
 439     charge = track->Charge();
 440     if(det==AliReducedEventFriend::kTPCpos && charge<0) continue;
 441     if(det==AliReducedEventFriend::kTPCneg && charge>0) continue;
 442
 443     if(IsTrackSelected && !IsTrackSelected(track)) continue;
 444     absWeight = 1.0;
 445     if(det==AliReducedEventFriend::kTPCptWeights) {
 446       absWeight = track->Pt();
 447       if(absWeight>2.0) absWeight = 2.0;    // pt is the weight used for the event plane
 448     }
 449     weight = absWeight;
 450     if(track->Eta()<0.0) weight *= -1.0;
 451
 452     ++nUsedTracks;
 453     x = TMath::Cos(track->Phi());
 454     y = TMath::Sin(track->Phi());
 455     //  1st harmonic
 456     Qvec[0][0] += weight*x;
 457     Qvec[0][1] += weight*y;
 458     //  2nd harmonic
 459     Qvec[1][0] += absWeight*(2.0*TMath::Power(x,2.0)-1);
 460     Qvec[1][1] += absWeight*(2.0*x*y);
 461     //  3rd harmonic
 462     Qvec[2][0] += weight*(4.0*TMath::Power(x,3.0)-3.0*x);
 463     Qvec[2][1] += weight*(3.0*y-4.0*TMath::Power(y,3.0));
 464     //  4th harmonic
 465     Qvec[3][0] += absWeight*(1.0-8.0*TMath::Power(x*y,2.0));
 466     Qvec[3][1] += absWeight*(4.0*x*y-8.0*x*TMath::Power(y,3.0));
 467     //  5th harmonic
 468     Qvec[4][0] += weight*(16.0*TMath::Power(x,5.0)-20.0*TMath::Power(x, 3.0)+5.0*x);
 469     Qvec[4][1] += weight*(16.0*TMath::Power(y,5.0)-20.0*TMath::Power(y, 3.0)+5.0*y);
 470     //  6th harmonic
 471     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);
 472     Qvec[5][1] += absWeight*(x*y*(32.0*TMath::Power(y,4.0)-32.0*TMath::Power(y, 2.0)+6.0));
 473   }
 474   return nUsedTracks;
 475 }
 476
 477
 478 //____________________________________________________________________________
 479 void AliReducedEvent::SubtractParticleFromQvector(
 480         AliReducedTrack* particle, Double_t Qvec[][2], Int_t det,
 481         Float_t etaMin/*=-0.8*/, Float_t etaMax/*=+0.8*/,
 482         Bool_t (*IsTrackSelected)(AliReducedTrack*)/*=NULL*/) {
 483   //
 484   // subtract a particle from the event Q-vector
 485   //
 486   Float_t eta = particle->Eta();
 487   if(eta<etaMin) return;
 488   if(eta>etaMax) return;
 489
 490   Float_t charge = particle->Charge();
 491   if(det==AliReducedEventFriend::kTPCpos && charge<0) return;
 492   if(det==AliReducedEventFriend::kTPCneg && charge>0) return;
 493
 494   if(IsTrackSelected && !IsTrackSelected(particle)) return;
 495
 496   Double_t weight=0.0; Double_t absWeight = 0.0;
 497   if(det==AliReducedEventFriend::kTPCptWeights) {
 498     absWeight = particle->Pt();
 499     if(absWeight>2.0) absWeight = 2.0;
 500   }
 501   weight = absWeight;
 502   if(eta<0.0) weight *= -1.0;
 503
 504   Float_t x = TMath::Cos(particle->Phi());
 505   Float_t y = TMath::Sin(particle->Phi());
 506
 507   //  1st harmonic
 508   Qvec[0][0] -= weight*x;
 509   Qvec[0][1] -= weight*y;
 510   //  2nd harmonic
 511   Qvec[1][0] -= absWeight*(2.0*TMath::Power(x,2.0)-1);
 512   Qvec[1][1] -= absWeight*(2.0*x*y);
 513   //  3rd harmonic
 514   Qvec[2][0] -= weight*(4.0*TMath::Power(x,3.0)-3.0*x);
 515   Qvec[2][1] -= weight*(3.0*y-4.0*TMath::Power(y,3.0));
 516   //  4th harmonic
 517   Qvec[3][0] -= absWeight*(1.0-8.0*TMath::Power(x*y,2.0));
 518   Qvec[3][1] -= absWeight*(4.0*x*y-8.0*x*TMath::Power(y,3.0));
 519   //  5th harmonic
 520   Qvec[4][0] -= weight*(16.0*TMath::Power(x,5.0)-20.0*TMath::Power(x, 3.0)+5.0*x);
 521   Qvec[4][1] -= weight*(16.0*TMath::Power(y,5.0)-20.0*TMath::Power(y, 3.0)+5.0*y);
 522   //  6th harmonic
 523   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);
 524   Qvec[5][1] -= absWeight*(x*y*(32.0*TMath::Power(y,4.0)-32.0*TMath::Power(y, 2.0)+6.0));
 525
 526   return;
 527 }
 528
 529
 530 //_________________________________________________________________
 531 void AliReducedEvent::GetVZEROQvector(Double_t Qvec[][2], Int_t det) {
 532   //
 533   // Get the reaction plane from the VZERO detector for a given harmonic
 534   //
 535   GetVZEROQvector(Qvec, det, fVZEROMult);
 536 }
 537
 538
 539 //_________________________________________________________________
 540 void AliReducedEvent::GetVZEROQvector(Double_t Qvec[][2], Int_t det, Float_t* vzeroMult) {
 541   //
 542   // Get the reaction plane from the VZERO detector for a given harmonic
 543   //
 544   //  Q{x,y} = SUM_i mult(i) * {cos(n*phi_i), sin(n*phi_i)}
 545   //  phi_i - phi angle of the VZERO sector i
 546   //          Each sector covers 45 degrees(8 sectors per ring). Middle of sector 0 is at 45/2
 547   //        channel 0: 22.5
 548   //                1: 22.5+45
 549   //                2: 22.5+45*2
 550   //               ...
 551   //        at the next ring continues the same
 552   //        channel 8: 22.5
 553   //        channel 9: 22.5 + 45
 554   //
 555   if(!(det==AliReducedEventFriend::kVZEROA ||
 556        det==AliReducedEventFriend::kVZEROC))
 557     return;
 558
 559   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)
 560   const Double_t kY[8] = {0.38268, 0.92388, 0.92388, 0.38268, -0.38268, -0.92388, -0.92388, -0.38268};    // sines     -- " --
 561   Int_t phi;
 562
 563   for(Int_t iChannel=0; iChannel<64; ++iChannel) {
 564     if(iChannel<32 && det==AliReducedEventFriend::kVZEROA) continue;
 565     if(iChannel>=32 && det==AliReducedEventFriend::kVZEROC) continue;
 566     phi=iChannel%8;
 567     //  1st harmonic
 568     Qvec[0][0] += vzeroMult[iChannel]*kX[phi];
 569     Qvec[0][1] += vzeroMult[iChannel]*kY[phi];
 570     //  2nd harmonic
 571     Qvec[1][0] += vzeroMult[iChannel]*(2.0*TMath::Power(kX[phi],2.0)-1);
 572     Qvec[1][1] += vzeroMult[iChannel]*(2.0*kX[phi]*kY[phi]);
 573     //  3rd harmonic
 574     Qvec[2][0] += vzeroMult[iChannel]*(4.0*TMath::Power(kX[phi],3.0)-3.0*kX[phi]);
 575     Qvec[2][1] += vzeroMult[iChannel]*(3.0*kY[phi]-4.0*TMath::Power(kY[phi],3.0));
 576     //  4th harmonic
 577     Qvec[3][0] += vzeroMult[iChannel]*(1.0-8.0*TMath::Power(kX[phi]*kY[phi],2.0));
 578     Qvec[3][1] += vzeroMult[iChannel]*(4.0*kX[phi]*kY[phi]-8.0*kX[phi]*TMath::Power(kY[phi],3.0));
 579     //  5th harmonic
 580     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]);
 581     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]);
 582     //  6th harmonic
 583     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);
 584     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));
 585   }    // end loop over channels
 586 }
 587
 588
 589 //_________________________________________________________________
 590 void AliReducedEvent::GetZDCQvector(Double_t Qvec[][2], Int_t det) const {
 591   //
 592   // Construct the event plane using the ZDC
 593   // ZDC has 2 side (A and C) with 4 calorimeters on each side
 594   // The XY position of each calorimeter is specified by the
 595   // zdcNTowerCenters_x and zdcNTowerCenters_y arrays
 596   if(!(det==AliReducedEventFriend::kZDCA ||
 597        det==AliReducedEventFriend::kZDCC)) return;       // bad detector option
 598   const Float_t zdcTowerCenter = 1.75;
 599   const Float_t zdcNTowerCentersX[4] = {-zdcTowerCenter,  zdcTowerCenter, -zdcTowerCenter, zdcTowerCenter};
 600   const Float_t zdcNTowerCentersY[4] = {-zdcTowerCenter, -zdcTowerCenter,  zdcTowerCenter, zdcTowerCenter};
 601
 602   Qvec[0][0] = 0.0; Qvec[0][1] = 0.0;   // first harmonic Q-vector
 603   Float_t zdcNCentroidSum = 0;
 604   Float_t zdcNalpha = 0.5;
 605
 606   for(Int_t i=0; i<4; ++i) {
 607     if(fZDCnEnergy[i+(det==AliReducedEventFriend::kZDCA ? 4 : 0)]>0.0) {
 608       Float_t zdcNenergyAlpha = TMath::Power(fZDCnEnergy[i+(det==AliReducedEventFriend::kZDCA ? 4 : 0)], zdcNalpha);
 609       Qvec[0][0] += zdcNTowerCentersX[i-1]*zdcNenergyAlpha;
 610       Qvec[0][1] += zdcNTowerCentersY[i-1]*zdcNenergyAlpha;
 611       zdcNCentroidSum += zdcNenergyAlpha;
 612     }
 613   }   // end loop over channels
 614
 615   if(zdcNCentroidSum>0.0) {
 616     Qvec[0][0] /= zdcNCentroidSum;
 617     Qvec[0][1] /= zdcNCentroidSum;
 618   }
 619 }