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