PWGJE/EMCALJetTasks/AliJetContainer.cxx

   1 // $Id$
   2 //
   3 // Container with name, TClonesArray and cuts for jets
   4 //
   5 // Author: M. Verweij
   6
   7 #include <TClonesArray.h>
   8
   9 #include "AliEmcalJet.h"
  10 #include "AliVEvent.h"
  11 #include "AliLog.h"
  12 #include "AliEMCALGeometry.h"
  13 #include "AliParticleContainer.h"
  14 #include "AliClusterContainer.h"
  15 #include "AliLocalRhoParameter.h"
  16
  17 #include "AliJetContainer.h"
  18
  19 ClassImp(AliJetContainer)
  20
  21 //________________________________________________________________________
  22 AliJetContainer::AliJetContainer():
  23   AliEmcalContainer("AliJetContainer"),
  24   fJetAcceptanceType(kUser),
  25   fJetRadius(0),
  26   fRhoName(),
  27   fLocalRhoName(),
  28   fFlavourSelection(0),
  29   fPtBiasJetTrack(0),
  30   fPtBiasJetClus(0),
  31   fJetPtCut(1),
  32   fJetAreaCut(-1),
  33   fAreaEmcCut(0),
  34   fJetMinEta(-0.9),
  35   fJetMaxEta(0.9),
  36   fJetMinPhi(-10),
  37   fJetMaxPhi(10),
  38   fMaxClusterPt(1000),
  39   fMaxTrackPt(100),
  40   fZLeadingEmcCut(10.),
  41   fZLeadingChCut(10.),
  42   fNEFMinCut(-10.),
  43   fNEFMaxCut(10.),
  44   fLeadingHadronType(0),
  45   fNLeadingJets(1),
  46   fJetBitMap(0),
  47   fJetTrigger(0),
  48   fParticleContainer(0),
  49   fClusterContainer(0),
  50   fRho(0),
  51   fLocalRho(0),
  52   fGeom(0),
  53   fRunNumber(0)
  54 {
  55   // Default constructor.
  56
  57   fClassName = "AliEmcalJet";
  58 }
  59
  60 //________________________________________________________________________
  61 AliJetContainer::AliJetContainer(const char *name):
  62   AliEmcalContainer(name),
  63   fJetAcceptanceType(kUser),
  64   fJetRadius(0),
  65   fRhoName(),
  66   fLocalRhoName(),
  67   fFlavourSelection(0),
  68   fPtBiasJetTrack(0),
  69   fPtBiasJetClus(0),
  70   fJetPtCut(1),
  71   fJetAreaCut(-1),
  72   fAreaEmcCut(0),
  73   fJetMinEta(-0.9),
  74   fJetMaxEta(0.9),
  75   fJetMinPhi(-10),
  76   fJetMaxPhi(10),
  77   fMaxClusterPt(1000),
  78   fMaxTrackPt(100),
  79   fZLeadingEmcCut(10.),
  80   fZLeadingChCut(10.),
  81   fNEFMinCut(-10.),
  82   fNEFMaxCut(10.),
  83   fLeadingHadronType(0),
  84   fNLeadingJets(1),
  85   fJetBitMap(0),
  86   fJetTrigger(0),
  87   fParticleContainer(0),
  88   fClusterContainer(0),
  89   fRho(0),
  90   fLocalRho(0),
  91   fGeom(0),
  92   fRunNumber(0)
  93 {
  94   // Standard constructor.
  95
  96   fClassName = "AliEmcalJet";
  97 }
  98
  99 //________________________________________________________________________
 100 void AliJetContainer::SetArray(AliVEvent *event)
 101 {
 102   // Set jet array
 103
 104   AliEmcalContainer::SetArray(event);
 105
 106   if(fJetAcceptanceType==kTPC) {
 107     AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
 108     SetJetEtaPhiTPC();
 109   }
 110   else if(fJetAcceptanceType==kEMCAL) {
 111     AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
 112     SetJetEtaPhiEMCAL();
 113  }
 114 }
 115
 116
 117 //________________________________________________________________________
 118 void AliJetContainer::SetEMCALGeometry() {
 119   fGeom = AliEMCALGeometry::GetInstance();
 120   if (!fGeom) {
 121     AliError(Form("%s: Can not create geometry", GetName()));
 122     return;
 123   }
 124 }
 125
 126 //________________________________________________________________________
 127 void AliJetContainer::LoadRho(AliVEvent *event)
 128 {
 129   // Load rho
 130
 131   if (!fRhoName.IsNull() && !fRho) {
 132     fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
 133     if (!fRho) {
 134       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
 135       return;
 136     }
 137   }
 138 }
 139
 140 //________________________________________________________________________
 141 void AliJetContainer::LoadLocalRho(AliVEvent *event)
 142 {
 143   // Load local rho
 144
 145   if (!fLocalRhoName.IsNull() && !fLocalRho) {
 146     fLocalRho = dynamic_cast<AliLocalRhoParameter*>(event->FindListObject(fLocalRhoName));
 147     if (!fLocalRho) {
 148       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fLocalRhoName.Data()));
 149       return;
 150     }
 151   }
 152 }
 153
 154 //________________________________________________________________________
 155 AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt)
 156 {
 157   // Get the leading jet; if opt contains "rho" the sorting is according to pt-A*rho
 158
 159   TString option(opt);
 160   option.ToLower();
 161
 162   Int_t tempID = fCurrentID;
 163
 164   AliEmcalJet *jetMax = GetNextAcceptJet(0);
 165   AliEmcalJet *jet = 0;
 166
 167   if (option.Contains("rho")) {
 168     while ((jet = GetNextAcceptJet())) {
 169       if ( (jet->Pt()-jet->Area()*GetRhoVal()) > (jetMax->Pt()-jetMax->Area()*GetRhoVal()) )
 170         jetMax = jet;
 171     }
 172   }
 173   else {
 174     while ((jet = GetNextAcceptJet())) {
 175       if (jet->Pt() > jetMax->Pt()) jetMax = jet;
 176     }
 177   }
 178
 179   fCurrentID = tempID;
 180
 181   return jetMax;
 182 }
 183
 184 //________________________________________________________________________
 185 AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
 186
 187   //Get i^th jet in array
 188
 189   if(i<0 || i>fClArray->GetEntriesFast()) return 0;
 190   AliEmcalJet *jet = static_cast<AliEmcalJet*>(fClArray->At(i));
 191   return jet;
 192
 193 }
 194
 195 //________________________________________________________________________
 196 AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
 197
 198   //Only return jet if is accepted
 199
 200   AliEmcalJet *jet = GetJet(i);
 201   if(!AcceptJet(jet)) return 0;
 202
 203   return jet;
 204 }
 205
 206 //________________________________________________________________________
 207 AliEmcalJet* AliJetContainer::GetJetWithLabel(Int_t lab) const {
 208
 209   //Get particle with label lab in array
 210
 211   Int_t i = GetIndexFromLabel(lab);
 212   return GetJet(i);
 213 }
 214
 215 //________________________________________________________________________
 216 AliEmcalJet* AliJetContainer::GetAcceptJetWithLabel(Int_t lab) const {
 217
 218   //Get particle with label lab in array
 219
 220   Int_t i = GetIndexFromLabel(lab);
 221   return GetAcceptJet(i);
 222 }
 223
 224 //________________________________________________________________________
 225 AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) {
 226
 227   //Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found
 228
 229   if (i>=0) fCurrentID = i;
 230
 231   const Int_t njets = GetNEntries();
 232   AliEmcalJet *jet = 0;
 233   while (fCurrentID < njets && !jet) {
 234     jet = GetAcceptJet(fCurrentID);
 235     fCurrentID++;
 236   }
 237
 238   return jet;
 239 }
 240
 241 //________________________________________________________________________
 242 AliEmcalJet* AliJetContainer::GetNextJet(Int_t i) {
 243
 244   //Get next jet; if i >= 0 (re)start counter from i; return 0 if no jet could be found
 245
 246   if (i>=0) fCurrentID = i;
 247
 248   const Int_t njets = GetNEntries();
 249   AliEmcalJet *jet = 0;
 250   while (fCurrentID < njets && !jet) {
 251     jet = GetJet(fCurrentID);
 252     fCurrentID++;
 253   }
 254
 255   return jet;
 256 }
 257
 258 //________________________________________________________________________
 259 Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
 260   AliEmcalJet *jet = GetJet(i);
 261
 262   return jet->Pt() - fRho->GetVal()*jet->Area();
 263 }
 264
 265 //________________________________________________________________________
 266 Double_t AliJetContainer::GetJetPtCorrLocal(Int_t i) const {
 267   AliEmcalJet *jet = GetJet(i);
 268
 269   return jet->Pt() - fLocalRho->GetLocalVal(jet->Phi(), fJetRadius)*jet->Area();
 270 }
 271
 272 //________________________________________________________________________
 273 void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
 274 {
 275   //Get momentum of the i^th jet in array
 276
 277   AliEmcalJet *jet = GetJet(i);
 278   if(jet) jet->GetMom(mom);
 279 }
 280
 281 //________________________________________________________________________
 282 Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
 283 {
 284   // Accept jet with a bias.
 285
 286   if (fLeadingHadronType == 0) {
 287     if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
 288   }
 289   else if (fLeadingHadronType == 1) {
 290     if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
 291   }
 292   else {
 293     if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
 294   }
 295
 296   return kTRUE;
 297
 298
 299 }
 300
 301 //________________________________________________________________________
 302 Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
 303 {
 304
 305    // Return true if jet is accepted.
 306
 307    if (!jet)
 308       return kFALSE;
 309
 310    if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
 311       return kFALSE;
 312
 313    if (jet->Pt() <= fJetPtCut)
 314       return kFALSE;
 315
 316    if (jet->Area() <= fJetAreaCut)
 317       return kFALSE;
 318
 319    if (jet->AreaEmc() < fAreaEmcCut)
 320       return kFALSE;
 321
 322    if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut)
 323       return kFALSE;
 324
 325    if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut)
 326       return kFALSE;
 327
 328    if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut)
 329       return kFALSE;
 330
 331    if (!AcceptBiasJet(jet))
 332       return kFALSE;
 333
 334    if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
 335       return kFALSE;
 336
 337    if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection))
 338       return kFALSE;
 339
 340    Double_t jetPhi = jet->Phi();
 341    Double_t jetEta = jet->Eta();
 342
 343    if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
 344       jetPhi -= TMath::Pi() * 2;
 345
 346    return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
 347 }
 348
 349 //________________________________________________________________________
 350 Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
 351 {
 352   if (fLeadingHadronType == 0)       // charged leading hadron
 353     return jet->MaxTrackPt();
 354   else if (fLeadingHadronType == 1)  // neutral leading hadron
 355     return jet->MaxClusterPt();
 356   else                               // charged or neutral
 357     return jet->MaxPartPt();
 358 }
 359
 360 //________________________________________________________________________
 361 void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet *jet) const
 362 {
 363   Double_t maxClusterPt = 0;
 364   Double_t maxClusterEta = 0;
 365   Double_t maxClusterPhi = 0;
 366
 367   Double_t maxTrackPt = 0;
 368   Double_t maxTrackEta = 0;
 369   Double_t maxTrackPhi = 0;
 370
 371   if (fClusterContainer && fClusterContainer->GetArray() && (fLeadingHadronType == 1 || fLeadingHadronType == 2)) {
 372     AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
 373     if (cluster) {
 374       TLorentzVector nPart;
 375       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 376
 377       maxClusterEta = nPart.Eta();
 378       maxClusterPhi = nPart.Phi();
 379       maxClusterPt = nPart.Pt();
 380     }
 381   }
 382
 383   if (fParticleContainer && fParticleContainer->GetArray() && (fLeadingHadronType == 0 || fLeadingHadronType == 2)) {
 384     AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
 385     if (track) {
 386       maxTrackEta = track->Eta();
 387       maxTrackPhi = track->Phi();
 388       maxTrackPt = track->Pt();
 389     }
 390   }
 391
 392   if (maxTrackPt > maxClusterPt)
 393     mom.SetPtEtaPhiM(maxTrackPt,maxTrackEta,maxTrackPhi,0.139);
 394   else
 395     mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
 396 }
 397
 398 //________________________________________________________________________
 399 Double_t AliJetContainer::GetZLeadingEmc(AliEmcalJet *jet) const
 400 {
 401
 402   if (fClusterContainer && fClusterContainer->GetArray()) {
 403     TLorentzVector mom;
 404
 405     AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
 406     if (cluster) {
 407       cluster->GetMomentum(mom, const_cast<Double_t*>(fVertex));
 408
 409       return GetZ(jet,mom);
 410     }
 411     else
 412       return -1;
 413   }
 414   else
 415     return -1;
 416 }
 417
 418 //________________________________________________________________________
 419 Double_t AliJetContainer::GetZLeadingCharged(AliEmcalJet *jet) const
 420 {
 421
 422   if (fParticleContainer && fParticleContainer->GetArray() ) {
 423     TLorentzVector mom;
 424
 425     AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
 426     if (track) {
 427       mom.SetPtEtaPhiM(track->Pt(),track->Eta(),track->Phi(),0.139);
 428
 429       return GetZ(jet,mom);
 430     }
 431     else
 432       return -1;
 433   }
 434   else
 435     return -1;
 436 }
 437
 438 //________________________________________________________________________
 439 Double_t AliJetContainer::GetZ(AliEmcalJet *jet, TLorentzVector mom) const
 440 {
 441
 442   Double_t pJetSq = jet->Px()*jet->Px() + jet->Py()*jet->Py() + jet->Pz()*jet->Pz();
 443
 444   if(pJetSq>1e-6)
 445     return (mom.Px()*jet->Px() + mom.Py()*jet->Py() + mom.Pz()*jet->Pz())/pJetSq;
 446   else {
 447     AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
 448     return -1;
 449   }
 450
 451 }
 452
 453 //________________________________________________________________________
 454 void AliJetContainer::SetJetEtaPhiEMCAL()
 455 {
 456   //Set default cuts for full jets
 457
 458   if(!fGeom) SetEMCALGeometry();
 459   if(fGeom) {
 460     SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
 461
 462     if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
 463       SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
 464     else
 465       SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
 466
 467   }
 468   else {
 469     AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
 470     SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
 471     SetJetPhiLimits(1.405+fJetRadius,3.135-fJetRadius);
 472   }
 473 }
 474
 475 //________________________________________________________________________
 476 void AliJetContainer::SetJetEtaPhiTPC()
 477 {
 478   //Set default cuts for charged jets
 479
 480   SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
 481   SetJetPhiLimits(-10, 10);
 482 }
 483
 484 //________________________________________________________________________
 485 void AliJetContainer::ResetCuts()
 486 {
 487   // Reset cuts to default values
 488
 489   fPtBiasJetTrack = 0;
 490   fPtBiasJetClus = 0;
 491   fJetPtCut = 1;
 492   fJetAreaCut = -1;
 493   fAreaEmcCut = 0;
 494   fJetMinEta = -0.9;
 495   fJetMaxEta = 0.9;
 496   fJetMinPhi = -10;
 497   fJetMaxPhi = 10;
 498   fMaxClusterPt = 1000;
 499   fMaxTrackPt = 100;
 500   fLeadingHadronType = 0;
 501   fZLeadingEmcCut = 10.;
 502   fZLeadingChCut = 10.;
 503 }
 504
 505 //________________________________________________________________________
 506 void AliJetContainer::SetClassName(const char *clname)
 507 {
 508   // Set the class name
 509
 510   TClass cls(clname);
 511   if (cls.InheritsFrom("AliEmcalJet")) fClassName = clname;
 512   else AliError(Form("Unable to set class name %s for a AliJetContainer, it must inherits from AliEmcalJet!",clname));
 513 }
 514
 515 //________________________________________________________________________
 516 Double_t AliJetContainer::GetFractionSharedPt(AliEmcalJet *jet1) const
 517 {
 518   //
 519   // Get fraction of shared pT between matched full and charged jet
 520   // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
 521   // Only works if tracks array of both jets is the same
 522   //
 523
 524   AliEmcalJet *jet2 = jet1->ClosestJet();
 525   if(!jet2) return -1;
 526
 527   Double_t fraction = 0.;
 528   Double_t jetPt2 = jet2->Pt();
 529
 530   if(jetPt2>0) {
 531     Double_t sumPt = 0.;
 532     AliVParticle *vpf = 0x0;
 533     Int_t iFound = 0;
 534     for(Int_t icc=0; icc<jet2->GetNumberOfTracks(); icc++) {
 535       Int_t idx = (Int_t)jet2->TrackAt(icc);
 536       iFound = 0;
 537       for(Int_t icf=0; icf<jet1->GetNumberOfTracks(); icf++) {
 538         if(idx == jet1->TrackAt(icf) && iFound==0 ) {
 539           iFound=1;
 540           vpf = static_cast<AliVParticle*>(jet1->TrackAt(icf, fParticleContainer->GetArray()));
 541           if(vpf) sumPt += vpf->Pt();
 542           continue;
 543         }
 544       }
 545     }
 546     fraction = sumPt/jetPt2;
 547   } else
 548     fraction = -1;
 549
 550   return fraction;
 551 }
 552