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   fLeadingHadronType(0),
  37   fNLeadingJets(1),
  38   fJetBitMap(0),
  39   fJetTrigger(0),
  40   fParticleContainer(0),
  41   fClusterContainer(0),
  42   fRho(0),
  43   fGeom(0),
  44   fRunNumber(0)
  45 {
  46   // Default constructor.
  47
  48 }
  49
  50 //________________________________________________________________________
  51 AliJetContainer::AliJetContainer(const char *name):
  52   AliEmcalContainer(name),
  53   fJetAcceptanceType(kUser),
  54   fJetRadius(0),
  55   fRhoName(),
  56   fPtBiasJetTrack(0),
  57   fPtBiasJetClus(0),
  58   fJetPtCut(1),
  59   fJetAreaCut(-1),
  60   fAreaEmcCut(0),
  61   fJetMinEta(-0.9),
  62   fJetMaxEta(0.9),
  63   fJetMinPhi(-10),
  64   fJetMaxPhi(10),
  65   fMaxClusterPt(1000),
  66   fMaxTrackPt(100),
  67   fLeadingHadronType(0),
  68   fNLeadingJets(1),
  69   fJetBitMap(0),
  70   fJetTrigger(0),
  71   fParticleContainer(0),
  72   fClusterContainer(0),
  73   fRho(0),
  74   fGeom(0),
  75   fRunNumber(0)
  76 {
  77   // Standard constructor.
  78
  79 }
  80
  81 //________________________________________________________________________
  82 void AliJetContainer::SetJetArray(AliVEvent *event)
  83 {
  84   // Set jet array
  85
  86   SetArray(event, "AliEmcalJet");
  87
  88   if(fJetAcceptanceType==kTPC) {
  89     AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
  90     SetJetEtaPhiTPC();
  91   }
  92   else if(fJetAcceptanceType==kEMCAL) {
  93     AliDebug(2,Form("%s: set EMCAL acceptance cuts",GetName()));
  94     SetJetEtaPhiEMCAL();
  95  }
  96 }
  97
  98
  99 //________________________________________________________________________
 100 void AliJetContainer::SetEMCALGeometry() {
 101   fGeom = AliEMCALGeometry::GetInstance();
 102   if (!fGeom) {
 103     AliError(Form("%s: Can not create geometry", GetName()));
 104     return;
 105   }
 106 }
 107
 108 //________________________________________________________________________
 109 void AliJetContainer::LoadRho(AliVEvent *event)
 110 {
 111   // Load rho
 112
 113   if (!fRhoName.IsNull() && !fRho) {
 114     fRho = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoName));
 115     if (!fRho) {
 116       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
 117       return;
 118     }
 119   }
 120 }
 121
 122 //________________________________________________________________________
 123 AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt) const
 124 {
 125   // Get the leading jet; if opt contains "rho" the sorting is according to pt-A*rho
 126
 127   TString option(opt);
 128   option.ToLower();
 129
 130   AliEmcalJet *jetMax = GetNextAcceptJet(0);
 131   AliEmcalJet *jet = 0;
 132
 133   if (option.Contains("rho")) {
 134     while ((jet = GetNextAcceptJet())) {
 135       if (jet->Pt()-jet->Area()*fRho->GetVal() > jetMax->Pt()-jetMax->Area()*fRho->GetVal()) jetMax = jet;
 136     }
 137   }
 138   else {
 139     while ((jet = GetNextAcceptJet())) {
 140       if (jet->Pt() > jetMax->Pt()) jetMax = jet;
 141     }
 142   }
 143
 144   return jetMax;
 145 }
 146
 147 //________________________________________________________________________
 148 AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
 149
 150   //Get i^th jet in array
 151
 152   if(i<0 || i>fClArray->GetEntriesFast()) return 0;
 153   AliEmcalJet *jet = static_cast<AliEmcalJet*>(fClArray->At(i));
 154   return jet;
 155
 156 }
 157
 158
 159 //________________________________________________________________________
 160 Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
 161   AliEmcalJet *jet = GetJet(i);
 162
 163   return jet->Pt() - fRho->GetVal()*jet->Area();
 164 }
 165
 166 //________________________________________________________________________
 167 AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
 168
 169   //Only return jet if is accepted
 170
 171   AliEmcalJet *jet = GetJet(i);
 172   if(!AcceptJet(jet)) return 0;
 173
 174   return jet;
 175 }
 176
 177 //________________________________________________________________________
 178 AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) const {
 179
 180   //Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found
 181
 182   static Int_t counter = -1;
 183   if (i>=0) counter = i;
 184
 185   const Int_t njets = GetNEntries();
 186   AliEmcalJet *jet = 0;
 187   while (counter < njets && !jet) {
 188     jet = GetAcceptJet(counter);
 189     counter++;
 190   }
 191
 192   return jet;
 193 }
 194
 195 //________________________________________________________________________
 196 void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
 197 {
 198   //Get momentum of the i^th jet in array
 199
 200   AliEmcalJet *jet = GetJet(i);
 201   if(jet) jet->GetMom(mom);
 202 }
 203
 204 //________________________________________________________________________
 205 Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
 206 {
 207   // Accept jet with a bias.
 208
 209   if (fLeadingHadronType == 0) {
 210     if (jet->MaxTrackPt() < fPtBiasJetTrack) return kFALSE;
 211   }
 212   else if (fLeadingHadronType == 1) {
 213     if (jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
 214   }
 215   else {
 216     if (jet->MaxTrackPt() < fPtBiasJetTrack && jet->MaxClusterPt() < fPtBiasJetClus) return kFALSE;
 217   }
 218
 219   return kTRUE;
 220
 221
 222 }
 223
 224 //________________________________________________________________________
 225 Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
 226 {
 227   // Return true if jet is accepted.
 228   if (!jet)
 229     return kFALSE;
 230
 231   if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
 232     return kFALSE;
 233
 234   if (jet->Pt() <= fJetPtCut)
 235     return kFALSE;
 236
 237   if (jet->Area() <= fJetAreaCut)
 238     return kFALSE;
 239
 240   if (jet->AreaEmc()<fAreaEmcCut)
 241     return kFALSE;
 242
 243   if (!AcceptBiasJet(jet))
 244     return kFALSE;
 245
 246   if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
 247     return kFALSE;
 248
 249   Double_t jetPhi = jet->Phi();
 250   Double_t jetEta = jet->Eta();
 251
 252   if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
 253     jetPhi -= TMath::Pi() * 2;
 254
 255   return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
 256 }
 257
 258 //________________________________________________________________________
 259 Double_t AliJetContainer::GetLeadingHadronPt(AliEmcalJet *jet) const
 260 {
 261   if (fLeadingHadronType == 0)       // charged leading hadron
 262     return jet->MaxTrackPt();
 263   else if (fLeadingHadronType == 1)  // neutral leading hadron
 264     return jet->MaxClusterPt();
 265   else                               // charged or neutral
 266     return jet->MaxPartPt();
 267 }
 268
 269 //________________________________________________________________________
 270 void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet *jet) const
 271 {
 272   Double_t maxClusterPt = 0;
 273   Double_t maxClusterEta = 0;
 274   Double_t maxClusterPhi = 0;
 275
 276   Double_t maxTrackPt = 0;
 277   Double_t maxTrackEta = 0;
 278   Double_t maxTrackPhi = 0;
 279
 280   if (fClusterContainer && fClusterContainer->GetArray() && (fLeadingHadronType == 1 || fLeadingHadronType == 2)) {
 281     AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
 282     if (cluster) {
 283       TLorentzVector nPart;
 284       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 285
 286       maxClusterEta = nPart.Eta();
 287       maxClusterPhi = nPart.Phi();
 288       maxClusterPt = nPart.Pt();
 289     }
 290   }
 291
 292   if (fParticleContainer && fParticleContainer->GetArray() && (fLeadingHadronType == 0 || fLeadingHadronType == 2)) {
 293     AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
 294     if (track) {
 295       maxTrackEta = track->Eta();
 296       maxTrackPhi = track->Phi();
 297       maxTrackPt = track->Pt();
 298     }
 299   }
 300
 301   if (maxTrackPt > maxClusterPt)
 302     mom.SetPtEtaPhiM(maxTrackPt,maxTrackEta,maxTrackPhi,0.139);
 303   else
 304     mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
 305 }
 306
 307 //________________________________________________________________________
 308 void AliJetContainer::SetJetEtaPhiEMCAL()
 309 {
 310   //Set default cuts for full jets
 311
 312   if(!fGeom) SetEMCALGeometry();
 313   if(fGeom) {
 314     SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
 315
 316     if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
 317       SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
 318     else
 319       SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
 320
 321   }
 322   else {
 323     AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
 324     SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
 325     SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
 326   }
 327 }
 328
 329 //________________________________________________________________________
 330 void AliJetContainer::SetJetEtaPhiTPC()
 331 {
 332   //Set default cuts for charged jets
 333
 334   SetJetEtaLimits(-0.9+fJetRadius, 0.9-fJetRadius);
 335   SetJetPhiLimits(-10, 10);
 336 }
 337
 338 //________________________________________________________________________
 339 void AliJetContainer::ResetCuts()
 340 {
 341   // Reset cuts to default values
 342
 343   fPtBiasJetTrack = 0;
 344   fPtBiasJetClus = 0;
 345   fJetPtCut = 1;
 346   fJetAreaCut = -1;
 347   fAreaEmcCut = 0;
 348   fJetMinEta = -0.9;
 349   fJetMaxEta = 0.9;
 350   fJetMinPhi = -10;
 351   fJetMaxPhi = 10;
 352   fMaxClusterPt = 1000;
 353   fMaxTrackPt = 100;
 354   fLeadingHadronType = 0;
 355
 356 }