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