PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliCorrelationReducedEvent.cxx

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