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()) ) jetMax = jet;
 170     }
 171   }
 172   else {
 173     while ((jet = GetNextAcceptJet())) {
 174       if (jet->Pt() > jetMax->Pt()) jetMax = jet;
 175     }
 176   }
 177
 178   fCurrentID = tempID;
 179
 180   return jetMax;
 181 }
 182
 183 //________________________________________________________________________
 184 AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
 185
 186   //Get i^th jet in array
 187
 188   if(i<0 || i>fClArray->GetEntriesFast()) return 0;
 189   AliEmcalJet *jet = static_cast<AliEmcalJet*>(fClArray->At(i));
 190   return jet;
 191
 192 }
 193
 194 //________________________________________________________________________
 195 AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
 196
 197   //Only return jet if is accepted
 198
 199   AliEmcalJet *jet = GetJet(i);
 200   if(!AcceptJet(jet)) return 0;
 201
 202   return jet;
 203 }
 204
 205 //________________________________________________________________________
 206 AliEmcalJet* AliJetContainer::GetJetWithLabel(Int_t lab) const {
 207
 208   //Get particle with label lab in array
 209
 210   Int_t i = GetIndexFromLabel(lab);
 211   return GetJet(i);
 212 }
 213
 214 //________________________________________________________________________
 215 AliEmcalJet* AliJetContainer::GetAcceptJetWithLabel(Int_t lab) const {
 216
 217   //Get particle with label lab in array
 218
 219   Int_t i = GetIndexFromLabel(lab);
 220   return GetAcceptJet(i);
 221 }
 222
 223 //________________________________________________________________________
 224 AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) {
 225
 226   //Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found
 227
 228   if (i>=0) fCurrentID = i;
 229
 230   const Int_t njets = GetNEntries();
 231   AliEmcalJet *jet = 0;
 232   while (fCurrentID < njets && !jet) {
 233     jet = GetAcceptJet(fCurrentID);
 234     fCurrentID++;
 235   }
 236
 237   return jet;
 238 }
 239
 240 //________________________________________________________________________
 241 AliEmcalJet* AliJetContainer::GetNextJet(Int_t i) {
 242
 243   //Get next jet; if i >= 0 (re)start counter from i; return 0 if no jet could be found
 244
 245   if (i>=0) fCurrentID = i;
 246
 247   const Int_t njets = GetNEntries();
 248   AliEmcalJet *jet = 0;
 249   while (fCurrentID < njets && !jet) {
 250     jet = GetJet(fCurrentID);
 251     fCurrentID++;
 252   }
 253
 254   return jet;
 255 }
 256
 257 //________________________________________________________________________
 258 Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
 259   AliEmcalJet *jet = GetJet(i);
 260
 261   return jet->Pt() - fRho->GetVal()*jet->Area();
 262 }
 263
 264 //________________________________________________________________________
 265 Double_t AliJetContainer::GetJetPtCorrLocal(Int_t i) const {
 266   AliEmcalJet *jet = GetJet(i);
 267
 268   return jet->Pt() - fLocalRho->GetLocalVal(jet->Phi(), fJetRadius)*jet->Area();
 269 }
 270
 271 //________________________________________________________________________
 272 void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
 273 {
 274   //Get momentum of the i^th jet in array
 275
 276   AliEmcalJet *jet = GetJet(i);
 277   if(jet) jet->GetMom(mom);
 278 }
 279
 280 //________________________________________________________________________
 281 Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
 282 {
 283   // Accept jet with a bias.
 284
 285   if (fLeadingHadronType == 0) {
 286     if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
 287   }
 288   else if (fLeadingHadronType == 1) {
 289     if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
 290   }
 291   else {
 292     if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
 293   }
 294
 295   return kTRUE;
 296
 297
 298 }
 299
 300 //________________________________________________________________________
 301 Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
 302 {
 303
 304    // Return true if jet is accepted.
 305
 306    if (!jet)
 307       return kFALSE;
 308
 309    if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
 310       return kFALSE;
 311
 312    if (jet->Pt() <= fJetPtCut)
 313       return kFALSE;
 314
 315    if (jet->Area() <= fJetAreaCut)
 316       return kFALSE;
 317
 318    if (jet->AreaEmc() < fAreaEmcCut)
 319       return kFALSE;
 320
 321    if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut)
 322       return kFALSE;
 323
 324    if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut)
 325       return kFALSE;
 326
 327    if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut)
 328       return kFALSE;
 329
 330    if (!AcceptBiasJet(jet))
 331       return kFALSE;
 332
 333    if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
 334       return kFALSE;
 335
 336    if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection))
 337       return kFALSE;
 338
 339    Double_t jetPhi = jet->Phi();
 340    Double_t jetEta = jet->Eta();
 341
 342    if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
 343       jetPhi -= TMath::Pi() * 2;
 344
 345    return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
 346 }
 347
 348 //________________________________________________________________________
 349 Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
 350 {
 351   if (fLeadingHadronType == 0)       // charged leading hadron
 352     return jet->MaxTrackPt();
 353   else if (fLeadingHadronType == 1)  // neutral leading hadron
 354     return jet->MaxClusterPt();
 355   else                               // charged or neutral
 356     return jet->MaxPartPt();
 357 }
 358
 359 //________________________________________________________________________
 360 void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet *jet) const
 361 {
 362   Double_t maxClusterPt = 0;
 363   Double_t maxClusterEta = 0;
 364   Double_t maxClusterPhi = 0;
 365
 366   Double_t maxTrackPt = 0;
 367   Double_t maxTrackEta = 0;
 368   Double_t maxTrackPhi = 0;
 369
 370   if (fClusterContainer && fClusterContainer->GetArray() && (fLeadingHadronType == 1 || fLeadingHadronType == 2)) {
 371     AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
 372     if (cluster) {
 373       TLorentzVector nPart;
 374       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 375
 376       maxClusterEta = nPart.Eta();
 377       maxClusterPhi = nPart.Phi();
 378       maxClusterPt = nPart.Pt();
 379     }
 380   }
 381
 382   if (fParticleContainer && fParticleContainer->GetArray() && (fLeadingHadronType == 0 || fLeadingHadronType == 2)) {
 383     AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
 384     if (track) {
 385       maxTrackEta = track->Eta();
 386       maxTrackPhi = track->Phi();
 387       maxTrackPt = track->Pt();
 388     }
 389   }
 390
 391   if (maxTrackPt > maxClusterPt)
 392     mom.SetPtEtaPhiM(maxTrackPt,maxTrackEta,maxTrackPhi,0.139);
 393   else
 394     mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
 395 }
 396
 397 //________________________________________________________________________
 398 Double_t AliJetContainer::GetZLeadingEmc(AliEmcalJet *jet) const
 399 {
 400
 401   if (fClusterContainer && fClusterContainer->GetArray()) {
 402     TLorentzVector mom;
 403
 404     AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
 405     if (cluster) {
 406       cluster->GetMomentum(mom, const_cast<Double_t*>(fVertex));
 407
 408       return GetZ(jet,mom);
 409     }
 410     else
 411       return -1;
 412   }
 413   else
 414     return -1;
 415 }
 416
 417 //________________________________________________________________________
 418 Double_t AliJetContainer::GetZLeadingCharged(AliEmcalJet *jet) const
 419 {
 420
 421   if (fParticleContainer && fParticleContainer->GetArray() ) {
 422     TLorentzVector mom;
 423
 424     AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
 425     if (track) {
 426       mom.SetPtEtaPhiM(track->Pt(),track->Eta(),track->Phi(),0.139);
 427
 428       return GetZ(jet,mom);
 429     }
 430     else
 431       return -1;
 432   }
 433   else
 434     return -1;
 435 }
 436
 437 //________________________________________________________________________
 438 Double_t AliJetContainer::GetZ(AliEmcalJet *jet, TLorentzVector mom) const
 439 {
 440
 441   Double_t pJetSq = jet->Px()*jet->Px() + jet->Py()*jet->Py() + jet->Pz()*jet->Pz();
 442
 443   if(pJetSq>1e-6)
 444     return (mom.Px()*jet->Px() + mom.Py()*jet->Py() + mom.Pz()*jet->Pz())/pJetSq;
 445   else {
 446     AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
 447     return -1;
 448   }
 449
 450 }
 451
 452 //________________________________________________________________________
 453 void AliJetContainer::SetJetEtaPhiEMCAL()
 454 {
 455   //Set default cuts for full jets
 456
 457   if(!fGeom) SetEMCALGeometry();
 458   if(fGeom) {
 459     SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
 460
 461     if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
 462       SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
 463     else
 464       SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
 465
 466   }
 467   else {
 468     AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
 469     SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
 470     SetJetPhiLimits(1.405+fJetRadius,3.135-fJetRadius);
 471   }
 472 }
 473
 474 //________________________________________________________________________
 475 void AliJetContainer::SetJetEtaPhiTPC()
 476 {
 477   //Set default cuts for charged jets
 478
 479   SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
 480   SetJetPhiLimits(-10, 10);
 481 }
 482
 483 //________________________________________________________________________
 484 void AliJetContainer::ResetCuts()
 485 {
 486   // Reset cuts to default values
 487
 488   fPtBiasJetTrack = 0;
 489   fPtBiasJetClus = 0;
 490   fJetPtCut = 1;
 491   fJetAreaCut = -1;
 492   fAreaEmcCut = 0;
 493   fJetMinEta = -0.9;
 494   fJetMaxEta = 0.9;
 495   fJetMinPhi = -10;
 496   fJetMaxPhi = 10;
 497   fMaxClusterPt = 1000;
 498   fMaxTrackPt = 100;
 499   fLeadingHadronType = 0;
 500   fZLeadingEmcCut = 10.;
 501   fZLeadingChCut = 10.;
 502 }
 503
 504 //________________________________________________________________________
 505 void AliJetContainer::SetClassName(const char *clname)
 506 {
 507   // Set the class name
 508
 509   TClass cls(clname);
 510   if (cls.InheritsFrom("AliEmcalJet")) fClassName = clname;
 511   else AliError(Form("Unable to set class name %s for a AliJetContainer, it must inherits from AliEmcalJet!",clname));
 512 }