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