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 {
  53   // Constructor.
  54
  55   fClosestJets[0] = 0;
  56   fClosestJets[1] = 0;
  57   fClosestJetsDist[0] = 999;
  58   fClosestJetsDist[1] = 999;
  59 }
  60
  61 //__________________________________________________________________________________________________
  62 AliEmcalJet::AliEmcalJet(Double_t px, Double_t py, Double_t pz) :
  63   AliVParticle(),
  64   fPt(TMath::Sqrt(px*px+py*py)),
  65   fEta(TMath::ASinH(pz/fPt)),
  66   fPhi(0),
  67   fM(0),
  68   fNEF(0),
  69   fArea(0),
  70   fAreaEta(0),
  71   fAreaPhi(0),
  72   fAreaEmc(-1),
  73   fAxisInEmcal(0),
  74   fFlavourTagging(0),
  75   fMaxCPt(0),
  76   fMaxNPt(0),
  77   fMCPt(0),
  78   fNn(0),
  79   fNch(0),
  80   fPtEmc(0),
  81   fNEmc(0),
  82   fClusterIDs(),
  83   fTrackIDs(),
  84   fMatched(2),
  85   fMatchingType(0),
  86   fTaggedJet(0x0),
  87   fTagStatus(-1),
  88   fPtSub(0),
  89   fPtVectSub(0),
  90   fTriggers(0),
  91   fJetShapeMassFirstDer(0),
  92   fJetShapeMassSecondDer(0),
  93   fJetShapeMassFirstSub(0),
  94   fJetShapeMassSecondSub(0),
  95   fLabel(-1),
  96   fGRNumerator(0),
  97   fGRDenominator(0),
  98   fGRNumeratorSub(0),
  99   fGRDenominatorSub(0)
 100 {
 101   // Constructor.
 102
 103   if (fPt != 0) {
 104     fPhi = TMath::ATan2(py, px);
 105     if (fPhi<0.)
 106       fPhi += 2. * TMath::Pi();
 107   }
 108
 109   fClosestJets[0] = 0;
 110   fClosestJets[1] = 0;
 111   fClosestJetsDist[0] = 999;
 112   fClosestJetsDist[1] = 999;
 113 }
 114
 115 //_________________________________________________________________________________________________
 116 AliEmcalJet::AliEmcalJet(Double_t pt, Double_t eta, Double_t phi, Double_t m) :
 117   AliVParticle(),
 118   fPt(pt),
 119   fEta(eta),
 120   fPhi(phi),
 121   fM(m),
 122   fNEF(0),
 123   fArea(0),
 124   fAreaEta(0),
 125   fAreaPhi(0),
 126   fAreaEmc(-1),
 127   fAxisInEmcal(0),
 128   fFlavourTagging(0),
 129   fMaxCPt(0),
 130   fMaxNPt(0),
 131   fMCPt(0),
 132   fNn(0),
 133   fNch(0),
 134   fPtEmc(0),
 135   fNEmc(0),
 136   fClusterIDs(),
 137   fTrackIDs(),
 138   fMatched(2),
 139   fMatchingType(0),
 140   fTaggedJet(0x0),
 141   fTagStatus(-1),
 142   fPtSub(0),
 143   fPtVectSub(0),
 144   fTriggers(0),
 145   fJetShapeMassFirstDer(0),
 146   fJetShapeMassSecondDer(0),
 147   fJetShapeMassFirstSub(0),
 148   fJetShapeMassSecondSub(0),
 149   fLabel(-1),
 150   fGRNumerator(0),
 151   fGRDenominator(0),
 152   fGRNumeratorSub(0),
 153   fGRDenominatorSub(0)
 154 {
 155   // Constructor.
 156
 157   if (fPhi<0.)
 158     fPhi += TMath::TwoPi();
 159
 160   fClosestJets[0] = 0;
 161   fClosestJets[1] = 0;
 162   fClosestJetsDist[0] = 999;
 163   fClosestJetsDist[1] = 999;
 164 }
 165
 166 //_________________________________________________________________________________________________
 167 AliEmcalJet::AliEmcalJet(const AliEmcalJet &jet) :
 168   AliVParticle(jet),
 169   fPt(jet.fPt),
 170   fEta(jet.fEta),
 171   fPhi(jet.fPhi),
 172   fM(jet.fM),
 173   fNEF(jet.fNEF),
 174   fArea(jet.fArea),
 175   fAreaEta(jet.fAreaEta),
 176   fAreaPhi(jet.fAreaPhi),
 177   fAreaEmc(jet.fAreaEmc),
 178   fAxisInEmcal(jet.fAxisInEmcal),
 179   fFlavourTagging(jet.fFlavourTagging),
 180   fMaxCPt(jet.fMaxCPt),
 181   fMaxNPt(jet.fMaxNPt),
 182   fMCPt(jet.fMCPt),
 183   fNn(jet.fNn),
 184   fNch(jet.fNch),
 185   fPtEmc(jet.fPtEmc),
 186   fNEmc(jet.fNEmc),
 187   fClusterIDs(jet.fClusterIDs),
 188   fTrackIDs(jet.fTrackIDs),
 189   fMatched(jet.fMatched),
 190   fMatchingType(jet.fMatchingType),
 191   fTaggedJet(jet.fTaggedJet),
 192   fTagStatus(jet.fTagStatus),
 193   fPtSub(jet.fPtSub),
 194   fPtVectSub(jet.fPtVectSub),
 195   fTriggers(jet.fTriggers),
 196   fJetShapeMassFirstDer(jet.fJetShapeMassFirstDer),
 197   fJetShapeMassSecondDer(jet.fJetShapeMassSecondDer),
 198   fJetShapeMassFirstSub(jet.fJetShapeMassFirstSub),
 199   fJetShapeMassSecondSub(jet.fJetShapeMassSecondSub),
 200   fLabel(jet.fLabel),
 201   fGRNumerator(jet.fGRNumerator),
 202   fGRDenominator(jet.fGRDenominator),
 203   fGRNumeratorSub(jet.fGRNumeratorSub),
 204   fGRDenominatorSub(jet.fGRDenominatorSub)
 205 {
 206   // Copy constructor.
 207
 208   fClosestJets[0]     = jet.fClosestJets[0];
 209   fClosestJets[1]     = jet.fClosestJets[1];
 210   fClosestJetsDist[0] = jet.fClosestJetsDist[0];
 211   fClosestJetsDist[1] = jet.fClosestJetsDist[1];
 212 }
 213
 214 //_________________________________________________________________________________________________
 215 AliEmcalJet &AliEmcalJet::operator=(const AliEmcalJet &jet)
 216 {
 217   // Assignment operator.
 218
 219   if (this!=&jet) {
 220     AliVParticle::operator=(jet);
 221     fPt                 = jet.fPt;
 222     fEta                = jet.fEta;
 223     fPhi                = jet.fPhi;
 224     fM                  = jet.fM;
 225     fNEF                = jet.fNEF;
 226     fArea               = jet.fArea;
 227     fAreaEta            = jet.fAreaEta;
 228     fAreaPhi            = jet.fAreaPhi;
 229     fAreaEmc            = jet.fAreaEmc;
 230     fAxisInEmcal        = jet.fAxisInEmcal;
 231     fFlavourTagging     = jet.fFlavourTagging;
 232     fMaxCPt             = jet.fMaxCPt;
 233     fMaxNPt             = jet.fMaxNPt;
 234     fMCPt               = jet.fMCPt;
 235     fNn                 = jet.fNn;
 236     fNch                = jet.fNch;
 237     fPtEmc              = jet.fPtEmc;
 238     fNEmc               = jet.fNEmc;
 239     fClusterIDs         = jet.fClusterIDs;
 240     fTrackIDs           = jet.fTrackIDs;
 241     fClosestJets[0]     = jet.fClosestJets[0];
 242     fClosestJets[1]     = jet.fClosestJets[1];
 243     fClosestJetsDist[0] = jet.fClosestJetsDist[0];
 244     fClosestJetsDist[1] = jet.fClosestJetsDist[1];
 245     fMatched            = jet.fMatched;
 246     fTaggedJet          = jet.fTaggedJet;
 247     fTagStatus          = jet.fTagStatus;
 248     fPtSub              = jet.fPtSub;
 249     fPtVectSub          = jet.fPtVectSub;
 250     fTriggers           = jet.fTriggers;
 251     fJetShapeMassFirstDer  = jet.fJetShapeMassFirstDer;
 252     fJetShapeMassSecondDer = jet.fJetShapeMassSecondDer;
 253     fJetShapeMassFirstSub  = jet.fJetShapeMassFirstSub;
 254     fJetShapeMassSecondSub = jet.fJetShapeMassSecondSub;
 255     fLabel              = jet.fLabel;
 256     fGRNumerator        = jet.fGRNumerator;
 257     fGRDenominator      = jet.fGRDenominator;
 258     fGRNumeratorSub     = jet.fGRNumeratorSub;
 259     fGRDenominatorSub   = jet.fGRDenominatorSub;
 260   }
 261
 262   return *this;
 263 }
 264
 265 //_________________________________________________________________________________________________
 266 Int_t AliEmcalJet::Compare(const TObject* obj) const
 267 {
 268   //Return -1 if this is smaller than obj, 0 if objects are equal and 1 if this is larger than obj.
 269
 270   const AliEmcalJet *jet = static_cast<const AliEmcalJet *>(obj);
 271   if (!obj)
 272     return 0;
 273   if (Pt()>jet->Pt())
 274     return -1;
 275   return 1;
 276 }
 277
 278 //__________________________________________________________________________________________________
 279 void AliEmcalJet::GetMom(TLorentzVector &vec) const
 280 {
 281   // Return momentum as four vector.
 282
 283   Double_t p = fPt *TMath::CosH(fEta);
 284   vec.SetPtEtaPhiE(fPt,fEta,fPhi,TMath::Sqrt(p*p+fM*fM));
 285 }
 286
 287 //__________________________________________________________________________________________________
 288 void AliEmcalJet::Print(Option_t* /*option*/) const
 289 {
 290   // Print jet information.
 291
 292   Printf("Jet pt=%.2f, eta=%.2f, phi=%.2f, area=%.2f, NEF=%.2f", fPt, fEta, fPhi, fArea, fNEF);
 293 }
 294
 295 //__________________________________________________________________________________________________
 296 Double_t AliEmcalJet::PtSubVect(Double_t rho) const
 297 {
 298   // Return vectorial subtracted transverse momentum.
 299
 300   Double_t dx = Px() - rho * fArea * TMath::Cos(fAreaPhi);
 301   Double_t dy = Py() - rho * fArea * TMath::Sin(fAreaPhi);
 302   //Double_t dz = Pz() - rho * fArea * TMath::SinH(fAreaEta);
 303   return TMath::Sqrt(dx*dx+dy*dy);
 304 }
 305
 306 //__________________________________________________________________________________________________
 307 void AliEmcalJet::SortConstituents()
 308 {
 309   // Sort constituent by index (increasing).
 310
 311   std::sort(fClusterIDs.GetArray(), fClusterIDs.GetArray() + fClusterIDs.GetSize());
 312   std::sort(fTrackIDs.GetArray(), fTrackIDs.GetArray() + fTrackIDs.GetSize());
 313 }
 314
 315 //__________________________________________________________________________________________________
 316 Double_t AliEmcalJet::DeltaR (const AliVParticle* part) const
 317 { // Helper function to calculate the distance between two jets or a jet and a particle
 318     Double_t dPhi = this->Phi() - part->Phi();
 319     Double_t dEta = this->Eta() - part->Eta();
 320     dPhi = TVector2::Phi_mpi_pi ( dPhi );
 321
 322     return TMath::Sqrt ( dPhi * dPhi + dEta * dEta );
 323 }
 324
 325
 326 //__________________________________________________________________________________________________
 327 std::vector<int> AliEmcalJet::SortConstituentsPt( TClonesArray *tracks ) const
 328 {   //___________________________________________
 329     // Sorting by p_T (decreasing) jet constituents
 330     //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 331     typedef std::pair<Double_t, Int_t> ptidx_pair;
 332
 333     // Create vector for Pt sorting
 334     std::vector<ptidx_pair> pair_list ;
 335
 336     for ( Int_t i_entry = 0; i_entry < GetNumberOfTracks(); i_entry++ )
 337         {
 338         AliVParticle *track = TrackAt(i_entry, tracks);
 339         if (!track)
 340             {
 341             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()));
 342             continue;
 343             }
 344
 345         pair_list.push_back( std::make_pair ( track->Pt(), i_entry ) );
 346         }
 347
 348     std::stable_sort( pair_list.begin() , pair_list.end() , sort_descend() );
 349
 350     // return an vector of indexes of constituents (sorted descending by pt)
 351     std::vector <int> index_sorted_list;
 352
 353     for ( std::vector< std::pair<Double_t,Int_t> >::iterator it = pair_list.begin(); it != pair_list.end(); ++it)
 354         { index_sorted_list.push_back( (*it).second ); } // populating the return object with indexes of sorted tracks
 355
 356     return index_sorted_list;
 357 }
 358
 359 //__________________________________________________________________________________________________
 360 AliVParticle* AliEmcalJet::GetLeadingTrack(TClonesArray *tracks) const
 361 {
 362   AliVParticle* maxTrack = 0;
 363   for (Int_t i = 0; i < GetNumberOfTracks(); i++) {
 364     AliVParticle *track = TrackAt(i, tracks);
 365     if (!track) {
 366       AliError(Form("Unable to find jet track %d in collection %s (pos in collection %d, max %d)",
 367                     i,tracks->GetName(),TrackAt(i),tracks->GetEntriesFast()));
 368       continue;
 369     }
 370     if (!maxTrack || track->Pt() > maxTrack->Pt())
 371       maxTrack = track;
 372   }
 373
 374   return maxTrack;
 375 }
 376
 377 //__________________________________________________________________________________________________
 378 AliVCluster* AliEmcalJet::GetLeadingCluster(TClonesArray *clusters) const
 379 {
 380   AliVCluster* maxCluster = 0;
 381   for (Int_t i = 0; i < GetNumberOfClusters(); i++) {
 382     AliVCluster *cluster = ClusterAt(i, clusters);
 383     if (!cluster) {
 384       AliError(Form("Unable to find jet cluster %d in collection %s (pos in collection %d, max %d)",
 385                     i,clusters->GetName(),ClusterAt(i),clusters->GetEntriesFast()));
 386       continue;
 387     }
 388     if (!maxCluster || cluster->E() > maxCluster->E())
 389       maxCluster = cluster;
 390   }
 391
 392   return maxCluster;
 393 }
 394
 395 //__________________________________________________________________________________________________
 396 void AliEmcalJet::ResetMatching()
 397 {
 398   fClosestJets[0] = 0;
 399   fClosestJets[1] = 0;
 400   fClosestJetsDist[0] = 999;
 401   fClosestJetsDist[1] = 999;
 402   fMatched = 2;
 403 }
 404
 405 //__________________________________________________________________________________________________
 406 void AliEmcalJet::PrintGR() {
 407   for(Int_t i = 0; i<fGRNumerator.GetSize(); i++) {
 408     Printf("num[%d] = %f",i,fGRNumerator.At(i));
 409   }
 410 }