PWGJE/EMCALJetTasks/AliEmcalJet.cxx

   1 //
   2 // Emcal jet class.
   3 //
   4 // Author: C.Loizides
   5
   6 #include "AliEmcalJet.h"
   7
   8 #include "AliLog.h"
   9 #include "Riostream.h"
  10
  11 ClassImp(AliEmcalJet)
  12
  13 //__________________________________________________________________________________________________
  14 AliEmcalJet::AliEmcalJet() :
  15   AliVParticle(),
  16   fPt(0),
  17   fEta(0),
  18   fPhi(0),
  19   fM(0),
  20   fNEF(0),
  21   fArea(0),
  22   fAreaEta(0),
  23   fAreaPhi(0),
  24   fAreaE(0),
  25   fAreaEmc(-1),
  26   fAxisInEmcal(0),
  27   fFlavourTagging(0),
  28   fMaxCPt(0),
  29   fMaxNPt(0),
  30   fMCPt(0),
  31   fNn(0),
  32   fNch(0),
  33   fPtEmc(0),
  34   fNEmc(0),
  35   fClusterIDs(),
  36   fTrackIDs(),
  37   fMatched(2),
  38   fMatchingType(0),
  39   fTaggedJet(0x0),
  40   fTagStatus(-1),
  41   fPtSub(0),
  42   fPtSubVect(0),
  43   fTriggers(0),
  44   fJetShapeMassFirstDer(0),
  45   fJetShapeMassSecondDer(0),
  46   fJetShapeMassFirstSub(0),
  47   fJetShapeMassSecondSub(0),
  48   fLabel(-1),
  49   fGRNumerator(0),
  50   fGRDenominator(0),
  51   fGRNumeratorSub(0),
  52   fGRDenominatorSub(0),
  53   fJetShapeAngularityFirstDer(0),
  54   fJetShapeAngularitySecondDer(0),
  55   fJetShapeAngularityFirstSub(0),
  56   fJetShapeAngularitySecondSub(0),
  57   fJetShapepTDFirstDer(0),
  58   fJetShapepTDSecondDer(0),
  59   fJetShapepTDFirstSub(0),
  60   fJetShapepTDSecondSub(0),
  61   fJetShapeCircularityFirstDer(0),
  62   fJetShapeCircularitySecondDer(0),
  63   fJetShapeCircularityFirstSub(0),
  64   fJetShapeCircularitySecondSub(0),
  65   fJetShapeSigma2FirstDer(0),
  66   fJetShapeSigma2SecondDer(0),
  67   fJetShapeSigma2FirstSub(0),
  68   fJetShapeSigma2SecondSub(0),
  69   fJetShapeConstituentFirstDer(0),
  70   fJetShapeConstituentSecondDer(0),
  71   fJetShapeConstituentFirstSub(0),
  72   fJetShapeConstituentSecondSub(0),
  73   fJetShapeLeSubFirstDer(0),
  74   fJetShapeLeSubSecondDer(0),
  75   fJetShapeLeSubFirstSub(0),
  76   fJetShapeLeSubSecondSub(0)
  77 {
  78   // Constructor.
  79   fClosestJets[0] = 0;
  80   fClosestJets[1] = 0;
  81   fClosestJetsDist[0] = 999;
  82   fClosestJetsDist[1] = 999;
  83 }
  84
  85 //__________________________________________________________________________________________________
  86 AliEmcalJet::AliEmcalJet(Double_t px, Double_t py, Double_t pz) :
  87   AliVParticle(),
  88   fPt(TMath::Sqrt(px * px + py* py)),
  89   fEta(TMath::ASinH(pz / fPt)),
  90   fPhi(0),
  91   fM(0),
  92   fNEF(0),
  93   fArea(0),
  94   fAreaEta(0),
  95   fAreaPhi(0),
  96   fAreaE(0),
  97   fAreaEmc(-1),
  98   fAxisInEmcal(0),
  99   fFlavourTagging(0),
 100   fMaxCPt(0),
 101   fMaxNPt(0),
 102   fMCPt(0),
 103   fNn(0),
 104   fNch(0),
 105   fPtEmc(0),
 106   fNEmc(0),
 107   fClusterIDs(),
 108   fTrackIDs(),
 109   fMatched(2),
 110   fMatchingType(0),
 111   fTaggedJet(0x0),
 112   fTagStatus(-1),
 113   fPtSub(0),
 114   fPtSubVect(0),
 115   fTriggers(0),
 116   fJetShapeMassFirstDer(0),
 117   fJetShapeMassSecondDer(0),
 118   fJetShapeMassFirstSub(0),
 119   fJetShapeMassSecondSub(0),
 120   fLabel(-1),
 121   fGRNumerator(0),
 122   fGRDenominator(0),
 123   fGRNumeratorSub(0),
 124   fGRDenominatorSub(0),
 125   fJetShapeAngularityFirstDer(0),
 126   fJetShapeAngularitySecondDer(0),
 127   fJetShapeAngularityFirstSub(0),
 128   fJetShapeAngularitySecondSub(0),
 129   fJetShapepTDFirstDer(0),
 130   fJetShapepTDSecondDer(0),
 131   fJetShapepTDFirstSub(0),
 132   fJetShapepTDSecondSub(0),
 133   fJetShapeCircularityFirstDer(0),
 134   fJetShapeCircularitySecondDer(0),
 135   fJetShapeCircularityFirstSub(0),
 136   fJetShapeCircularitySecondSub(0),
 137   fJetShapeSigma2FirstDer(0),
 138   fJetShapeSigma2SecondDer(0),
 139   fJetShapeSigma2FirstSub(0),
 140   fJetShapeSigma2SecondSub(0),
 141   fJetShapeConstituentFirstDer(0),
 142   fJetShapeConstituentSecondDer(0),
 143   fJetShapeConstituentFirstSub(0),
 144   fJetShapeConstituentSecondSub(0),
 145   fJetShapeLeSubFirstDer(0),
 146   fJetShapeLeSubSecondDer(0),
 147   fJetShapeLeSubFirstSub(0),
 148   fJetShapeLeSubSecondSub(0)
 149 {
 150   // Constructor.
 151
 152   if(fPt != 0) {
 153     fPhi = TVector2::Phi_0_2pi(TMath::ATan2(py, px));
 154   }
 155
 156   fClosestJets[0] = 0;
 157   fClosestJets[1] = 0;
 158   fClosestJetsDist[0] = 999;
 159   fClosestJetsDist[1] = 999;
 160 }
 161
 162 //_________________________________________________________________________________________________
 163 AliEmcalJet::AliEmcalJet(Double_t pt, Double_t eta, Double_t phi, Double_t m) :
 164   AliVParticle(),
 165   fPt(pt),
 166   fEta(eta),
 167   fPhi(phi),
 168   fM(m),
 169   fNEF(0),
 170   fArea(0),
 171   fAreaEta(0),
 172   fAreaPhi(0),
 173   fAreaE(0),
 174   fAreaEmc(-1),
 175   fAxisInEmcal(0),
 176   fFlavourTagging(0),
 177   fMaxCPt(0),
 178   fMaxNPt(0),
 179   fMCPt(0),
 180   fNn(0),
 181   fNch(0),
 182   fPtEmc(0),
 183   fNEmc(0),
 184   fClusterIDs(),
 185   fTrackIDs(),
 186   fMatched(2),
 187   fMatchingType(0),
 188   fTaggedJet(0x0),
 189   fTagStatus(-1),
 190   fPtSub(0),
 191   fPtSubVect(0),
 192   fTriggers(0),
 193   fJetShapeMassFirstDer(0),
 194   fJetShapeMassSecondDer(0),
 195   fJetShapeMassFirstSub(0),
 196   fJetShapeMassSecondSub(0),
 197   fLabel(-1),
 198   fGRNumerator(0),
 199   fGRDenominator(0),
 200   fGRNumeratorSub(0),
 201   fGRDenominatorSub(0),
 202   fJetShapeAngularityFirstDer(0),
 203   fJetShapeAngularitySecondDer(0),
 204   fJetShapeAngularityFirstSub(0),
 205   fJetShapeAngularitySecondSub(0),
 206   fJetShapepTDFirstDer(0),
 207   fJetShapepTDSecondDer(0),
 208   fJetShapepTDFirstSub(0),
 209   fJetShapepTDSecondSub(0),
 210   fJetShapeCircularityFirstDer(0),
 211   fJetShapeCircularitySecondDer(0),
 212   fJetShapeCircularityFirstSub(0),
 213   fJetShapeCircularitySecondSub(0),
 214   fJetShapeSigma2FirstDer(0),
 215   fJetShapeSigma2SecondDer(0),
 216   fJetShapeSigma2FirstSub(0),
 217   fJetShapeSigma2SecondSub(0),
 218   fJetShapeConstituentFirstDer(0),
 219   fJetShapeConstituentSecondDer(0),
 220   fJetShapeConstituentFirstSub(0),
 221   fJetShapeConstituentSecondSub(0),
 222   fJetShapeLeSubFirstDer(0),
 223   fJetShapeLeSubSecondDer(0),
 224   fJetShapeLeSubFirstSub(0),
 225   fJetShapeLeSubSecondSub(0)
 226
 227 {
 228   // Constructor.
 229
 230   fPhi = TVector2::Phi_0_2pi(fPhi);
 231
 232   fClosestJets[0] = 0;
 233   fClosestJets[1] = 0;
 234   fClosestJetsDist[0] = 999;
 235   fClosestJetsDist[1] = 999;
 236 }
 237
 238 //_________________________________________________________________________________________________
 239 AliEmcalJet::AliEmcalJet(const AliEmcalJet& jet) :
 240   AliVParticle(jet),
 241   fPt(jet.fPt),
 242   fEta(jet.fEta),
 243   fPhi(jet.fPhi),
 244   fM(jet.fM),
 245   fNEF(jet.fNEF),
 246   fArea(jet.fArea),
 247   fAreaEta(jet.fAreaEta),
 248   fAreaPhi(jet.fAreaPhi),
 249   fAreaE(jet.fAreaE),
 250   fAreaEmc(jet.fAreaEmc),
 251   fAxisInEmcal(jet.fAxisInEmcal),
 252   fFlavourTagging(jet.fFlavourTagging),
 253   fMaxCPt(jet.fMaxCPt),
 254   fMaxNPt(jet.fMaxNPt),
 255   fMCPt(jet.fMCPt),
 256   fNn(jet.fNn),
 257   fNch(jet.fNch),
 258   fPtEmc(jet.fPtEmc),
 259   fNEmc(jet.fNEmc),
 260   fClusterIDs(jet.fClusterIDs),
 261   fTrackIDs(jet.fTrackIDs),
 262   fMatched(jet.fMatched),
 263   fMatchingType(jet.fMatchingType),
 264   fTaggedJet(jet.fTaggedJet),
 265   fTagStatus(jet.fTagStatus),
 266   fPtSub(jet.fPtSub),
 267   fPtSubVect(jet.fPtSubVect),
 268   fTriggers(jet.fTriggers),
 269   fJetShapeMassFirstDer(jet.fJetShapeMassFirstDer),
 270   fJetShapeMassSecondDer(jet.fJetShapeMassSecondDer),
 271   fJetShapeMassFirstSub(jet.fJetShapeMassFirstSub),
 272   fJetShapeMassSecondSub(jet.fJetShapeMassSecondSub),
 273   fLabel(jet.fLabel),
 274   fGRNumerator(jet.fGRNumerator),
 275   fGRDenominator(jet.fGRDenominator),
 276   fGRNumeratorSub(jet.fGRNumeratorSub),
 277   fGRDenominatorSub(jet.fGRDenominatorSub),
 278   fJetShapeAngularityFirstDer(jet.fJetShapeAngularityFirstDer),
 279   fJetShapeAngularitySecondDer(jet.fJetShapeAngularitySecondDer),
 280   fJetShapeAngularityFirstSub(jet.fJetShapeAngularityFirstSub),
 281   fJetShapeAngularitySecondSub(jet.fJetShapeAngularitySecondSub),
 282   fJetShapepTDFirstDer(jet.fJetShapepTDFirstDer),
 283   fJetShapepTDSecondDer(jet.fJetShapepTDSecondDer),
 284   fJetShapepTDFirstSub(jet.fJetShapepTDFirstSub),
 285   fJetShapepTDSecondSub(jet.fJetShapepTDSecondSub),
 286   fJetShapeCircularityFirstDer(jet.fJetShapeCircularityFirstDer),
 287   fJetShapeCircularitySecondDer(jet.fJetShapeCircularitySecondDer),
 288   fJetShapeCircularityFirstSub(jet.fJetShapeCircularityFirstSub),
 289   fJetShapeCircularitySecondSub(jet.fJetShapeCircularitySecondSub),
 290   fJetShapeSigma2FirstDer(jet.fJetShapeSigma2FirstDer),
 291   fJetShapeSigma2SecondDer(jet.fJetShapeSigma2SecondDer),
 292   fJetShapeSigma2FirstSub(jet.fJetShapeSigma2FirstSub),
 293   fJetShapeSigma2SecondSub(jet.fJetShapeSigma2SecondSub),
 294   fJetShapeConstituentFirstDer(jet.fJetShapeConstituentFirstDer),
 295   fJetShapeConstituentSecondDer(jet.fJetShapeConstituentSecondDer),
 296   fJetShapeConstituentFirstSub(jet.fJetShapeConstituentFirstSub),
 297   fJetShapeConstituentSecondSub(jet.fJetShapeConstituentSecondSub),
 298   fJetShapeLeSubFirstDer(jet.fJetShapeLeSubFirstDer),
 299   fJetShapeLeSubSecondDer(jet.fJetShapeLeSubSecondDer),
 300   fJetShapeLeSubFirstSub(jet.fJetShapeLeSubFirstSub),
 301   fJetShapeLeSubSecondSub(jet.fJetShapeLeSubSecondSub)
 302 {
 303   // Copy constructor.
 304   fClosestJets[0]     = jet.fClosestJets[0];
 305   fClosestJets[1]     = jet.fClosestJets[1];
 306   fClosestJetsDist[0] = jet.fClosestJetsDist[0];
 307   fClosestJetsDist[1] = jet.fClosestJetsDist[1];
 308 }
 309
 310 //_________________________________________________________________________________________________
 311 AliEmcalJet& AliEmcalJet::operator=(const AliEmcalJet& jet)
 312 {
 313   // Assignment operator.
 314
 315   if(this != &jet) {
 316     AliVParticle::operator=(jet);
 317     fPt                 = jet.fPt;
 318     fEta                = jet.fEta;
 319     fPhi                = jet.fPhi;
 320     fM                  = jet.fM;
 321     fNEF                = jet.fNEF;
 322     fArea               = jet.fArea;
 323     fAreaEta            = jet.fAreaEta;
 324     fAreaPhi            = jet.fAreaPhi;
 325     fAreaE              = jet.fAreaE;
 326     fAreaEmc            = jet.fAreaEmc;
 327     fAxisInEmcal        = jet.fAxisInEmcal;
 328     fFlavourTagging     = jet.fFlavourTagging;
 329     fMaxCPt             = jet.fMaxCPt;
 330     fMaxNPt             = jet.fMaxNPt;
 331     fMCPt               = jet.fMCPt;
 332     fNn                 = jet.fNn;
 333     fNch                = jet.fNch;
 334     fPtEmc              = jet.fPtEmc;
 335     fNEmc               = jet.fNEmc;
 336     fClusterIDs         = jet.fClusterIDs;
 337     fTrackIDs           = jet.fTrackIDs;
 338     fClosestJets[0]     = jet.fClosestJets[0];
 339     fClosestJets[1]     = jet.fClosestJets[1];
 340     fClosestJetsDist[0] = jet.fClosestJetsDist[0];
 341     fClosestJetsDist[1] = jet.fClosestJetsDist[1];
 342     fMatched            = jet.fMatched;
 343     fTaggedJet          = jet.fTaggedJet;
 344     fTagStatus          = jet.fTagStatus;
 345     fPtSub              = jet.fPtSub;
 346     fPtSubVect          = jet.fPtSubVect;
 347     fTriggers           = jet.fTriggers;
 348     fJetShapeMassFirstDer  = jet.fJetShapeMassFirstDer;
 349     fJetShapeMassSecondDer = jet.fJetShapeMassSecondDer;
 350     fJetShapeMassFirstSub  = jet.fJetShapeMassFirstSub;
 351     fJetShapeMassSecondSub = jet.fJetShapeMassSecondSub;
 352     fLabel              = jet.fLabel;
 353     fGRNumerator        = jet.fGRNumerator;
 354     fGRDenominator      = jet.fGRDenominator;
 355     fGRNumeratorSub     = jet.fGRNumeratorSub;
 356     fGRDenominatorSub   = jet.fGRDenominatorSub;
 357     fJetShapeAngularityFirstDer  = jet.fJetShapeAngularityFirstDer;
 358     fJetShapeAngularitySecondDer = jet.fJetShapeAngularitySecondDer;
 359     fJetShapeAngularityFirstSub  = jet.fJetShapeAngularityFirstSub;
 360     fJetShapeAngularitySecondSub = jet.fJetShapeAngularitySecondSub;
 361     fJetShapepTDFirstDer  = jet.fJetShapepTDFirstDer;
 362     fJetShapepTDSecondDer = jet.fJetShapepTDSecondDer;
 363     fJetShapepTDFirstSub  = jet.fJetShapepTDFirstSub;
 364     fJetShapepTDSecondSub = jet.fJetShapepTDSecondSub;
 365     fJetShapeCircularityFirstDer  = jet.fJetShapeCircularityFirstDer;
 366     fJetShapeCircularitySecondDer = jet.fJetShapeCircularitySecondDer;
 367     fJetShapeCircularityFirstSub  = jet.fJetShapeCircularityFirstSub;
 368     fJetShapeCircularitySecondSub = jet.fJetShapeCircularitySecondSub;
 369     fJetShapeSigma2FirstDer  = jet.fJetShapeSigma2FirstDer;
 370     fJetShapeSigma2SecondDer = jet.fJetShapeSigma2SecondDer;
 371     fJetShapeSigma2FirstSub  = jet.fJetShapeSigma2FirstSub;
 372     fJetShapeSigma2SecondSub = jet.fJetShapeSigma2SecondSub;
 373     fJetShapeConstituentFirstDer  = jet.fJetShapeConstituentFirstDer;
 374     fJetShapeConstituentSecondDer = jet.fJetShapeConstituentSecondDer;
 375     fJetShapeConstituentFirstSub  = jet.fJetShapeConstituentFirstSub;
 376     fJetShapeConstituentSecondSub = jet.fJetShapeConstituentSecondSub;
 377     fJetShapeLeSubFirstDer  = jet.fJetShapeLeSubFirstDer;
 378     fJetShapeLeSubSecondDer = jet.fJetShapeLeSubSecondDer;
 379     fJetShapeLeSubFirstSub  = jet.fJetShapeLeSubFirstSub;
 380     fJetShapeLeSubSecondSub = jet.fJetShapeLeSubSecondSub;
 381   }
 382
 383   return *this;
 384 }
 385
 386 //_________________________________________________________________________________________________
 387 Int_t AliEmcalJet::Compare(const TObject* obj) const
 388 {
 389   //Return -1 if this is smaller than obj, 0 if objects are equal and 1 if this is larger than obj.
 390
 391   const AliEmcalJet* jet = static_cast<const AliEmcalJet*>(obj);
 392   if(!obj)
 393     return 0;
 394   if(Pt() > jet->Pt())
 395     return -1;
 396   return 1;
 397 }
 398
 399 //__________________________________________________________________________________________________
 400 void AliEmcalJet::GetMom(TLorentzVector& vec) const
 401 {
 402   // Return momentum as four-vector.
 403
 404   vec.SetPtEtaPhiE(fPt, fEta, fPhi, E());
 405 }
 406
 407 //__________________________________________________________________________________________________
 408 void AliEmcalJet::Print(Option_t* /*option*/) const
 409 {
 410   // Print jet information.
 411
 412   Printf("Jet pt=%.2f, eta=%.2f, phi=%.2f, area=%.2f, NEF=%.2f", fPt, fEta, fPhi, fArea, fNEF);
 413 }
 414
 415 //__________________________________________________________________________________________________
 416 Double_t AliEmcalJet::PtSub(Double_t rho, Bool_t save)
 417 {
 418   // Return transverse momentum after scalar subtraction. Save the result if required.
 419   // Result can be negative.
 420
 421   Double_t ptcorr = fPt - rho * fArea;
 422   if(save)
 423     fPtSub = ptcorr;
 424   return ptcorr;
 425 }
 426
 427 //__________________________________________________________________________________________________
 428 Double_t AliEmcalJet::PtSubVect(Double_t rho, Bool_t save)
 429 {
 430   // Return transverse momentum after vectorial subtraction. Save the result if required.
 431   // Result cannot be negative.
 432
 433   Double_t dx = Px() - rho * fArea * TMath::Cos(fAreaPhi);
 434   Double_t dy = Py() - rho * fArea * TMath::Sin(fAreaPhi);
 435   //Double_t dz = Pz() - rho * fArea * TMath::SinH(fAreaEta);
 436   Double_t ptcorr = TMath::Sqrt(dx * dx + dy * dy);
 437   if(save)
 438     fPtSubVect = ptcorr;
 439   return ptcorr;
 440 }
 441
 442 //__________________________________________________________________________________________________
 443 TLorentzVector AliEmcalJet::SubtractRhoVect(Double_t rho, Bool_t save)
 444 {
 445   // Return four-momentum after vectorial subtraction. Save pt if required.
 446   // Saved value of pt is negative if the corrected momentum is pointing to the opposite half-plane in the x-y plane w.r.t. the raw momentum.
 447
 448   TLorentzVector vecCorr;
 449   GetMom(vecCorr);
 450   TLorentzVector vecBg;
 451   vecBg.SetPtEtaPhiE(fArea, fAreaEta, fAreaPhi, fAreaE);
 452   vecBg *= rho;
 453   vecCorr -= vecBg;
 454   if(save)
 455   {
 456     Double_t dPhi = TMath::Abs(TVector2::Phi_mpi_pi(Phi() - vecCorr.Phi()));
 457     Int_t signum = dPhi <= TMath::PiOver2() ? 1 : -1;
 458     fPtSubVect = signum * vecCorr.Pt();
 459   }
 460   return vecCorr;
 461 }
 462
 463 //__________________________________________________________________________________________________
 464 void AliEmcalJet::SortConstituents()
 465 {
 466   // Sort constituent by index (increasing).
 467
 468   std::sort(fClusterIDs.GetArray(), fClusterIDs.GetArray() + fClusterIDs.GetSize());
 469   std::sort(fTrackIDs.GetArray(), fTrackIDs.GetArray() + fTrackIDs.GetSize());
 470 }
 471
 472 //__________________________________________________________________________________________________
 473 Double_t AliEmcalJet::DeltaR(const AliVParticle* part) const
 474 {
 475   // Helper function to calculate the distance between two jets or a jet and a particle
 476
 477   Double_t dPhi = Phi() - part->Phi();
 478   Double_t dEta = Eta() - part->Eta();
 479   dPhi = TVector2::Phi_mpi_pi(dPhi);
 480   return TMath::Sqrt(dPhi * dPhi + dEta * dEta);
 481 }
 482
 483
 484 //__________________________________________________________________________________________________
 485 std::vector<int> AliEmcalJet::SortConstituentsPt(TClonesArray* tracks) const
 486 {
 487   // Sorting by p_T (decreasing) jet constituents
 488
 489   typedef std::pair<Double_t, Int_t> ptidx_pair;
 490
 491   // Create vector for Pt sorting
 492   std::vector<ptidx_pair> pair_list ;
 493
 494   for(Int_t i_entry = 0; i_entry < GetNumberOfTracks(); i_entry++)
 495   {
 496     AliVParticle* track = TrackAt(i_entry, tracks);
 497     if(!track)
 498     {
 499       AliError(Form("Unable to find jet track %d in collection %s (pos in collection %d, max %d)", i_entry, tracks->GetName(), TrackAt(i_entry), tracks->GetEntriesFast()));
 500       continue;
 501     }
 502
 503     pair_list.push_back(std::make_pair(track->Pt(), i_entry));
 504   }
 505
 506   std::stable_sort(pair_list.begin() , pair_list.end() , sort_descend());
 507
 508   // return an vector of indexes of constituents (sorted descending by pt)
 509   std::vector <int> index_sorted_list;
 510
 511   for(std::vector< std::pair<Double_t, Int_t> >::iterator it = pair_list.begin(); it != pair_list.end(); ++it)
 512   { index_sorted_list.push_back((*it).second); }   // populating the return object with indexes of sorted tracks
 513
 514   return index_sorted_list;
 515 }
 516
 517 //________________________________________________________________________
 518 Double_t AliEmcalJet::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz) const
 519 {
 520   // Get the z of a constituent inside of a jet
 521
 522   Double_t pJetSq = P();
 523   pJetSq *= pJetSq;
 524
 525   if(pJetSq > 1e-6)
 526   { return (trkPx * Px() + trkPy * Py() + trkPz * Pz()) / pJetSq ; }
 527   else
 528   { AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f", GetName(), pJetSq)); return -1; }
 529 }
 530
 531 //________________________________________________________________________
 532 Double_t AliEmcalJet::GetZ(const AliVParticle* trk) const
 533 {
 534   // Get Z of constituent trk
 535
 536   return GetZ(trk->Px(), trk->Py(), trk->Pz());
 537 }
 538
 539 //__________________________________________________________________________________________________
 540 AliVParticle* AliEmcalJet::GetLeadingTrack(TClonesArray* tracks) const
 541 {
 542   AliVParticle* maxTrack = 0;
 543   for(Int_t i = 0; i < GetNumberOfTracks(); i++) {
 544     AliVParticle* track = TrackAt(i, tracks);
 545     if(!track) {
 546       AliError(Form("Unable to find jet track %d in collection %s (pos in collection %d, max %d)",
 547                     i, tracks->GetName(), TrackAt(i), tracks->GetEntriesFast()));
 548       continue;
 549     }
 550     if(!maxTrack || track->Pt() > maxTrack->Pt())
 551       maxTrack = track;
 552   }
 553
 554   return maxTrack;
 555 }
 556
 557 //__________________________________________________________________________________________________
 558 AliVCluster* AliEmcalJet::GetLeadingCluster(TClonesArray* clusters) const
 559 {
 560   AliVCluster* maxCluster = 0;
 561   for(Int_t i = 0; i < GetNumberOfClusters(); i++) {
 562     AliVCluster* cluster = ClusterAt(i, clusters);
 563     if(!cluster) {
 564       AliError(Form("Unable to find jet cluster %d in collection %s (pos in collection %d, max %d)",
 565                     i, clusters->GetName(), ClusterAt(i), clusters->GetEntriesFast()));
 566       continue;
 567     }
 568     if(!maxCluster || cluster->E() > maxCluster->E())
 569       maxCluster = cluster;
 570   }
 571
 572   return maxCluster;
 573 }
 574
 575 //__________________________________________________________________________________________________
 576 void AliEmcalJet::ResetMatching()
 577 {
 578   fClosestJets[0] = 0;
 579   fClosestJets[1] = 0;
 580   fClosestJetsDist[0] = 999;
 581   fClosestJetsDist[1] = 999;
 582   fMatched = 2;
 583 }
 584
 585 //__________________________________________________________________________________________________
 586 void AliEmcalJet::PrintGR()
 587 {
 588   for(Int_t i = 0; i < fGRNumerator.GetSize(); i++) {
 589     Printf("num[%d] = %f", i, fGRNumerator.At(i));
 590   }
 591 }