PWGJE/EMCALJetTasks/AliJetContainer.cxx

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