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