]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/EMCALJetTasks/AliJetContainer.cxx
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetContainer.cxx
index b6b3d0da12785735f83b5807a575933900b2d722..a1a0b60325c4cca9a7b5975e2b7c9ea92d8fb314 100644 (file)
@@ -1,5 +1,6 @@
+// $Id$
 //
-// Jet Container
+// Container with name, TClonesArray and cuts for jets
 //
 // Author: M. Verweij
 
@@ -11,6 +12,7 @@
 #include "AliEMCALGeometry.h"
 #include "AliParticleContainer.h"
 #include "AliClusterContainer.h"
+#include "AliLocalRhoParameter.h"
 
 #include "AliJetContainer.h"
 
@@ -22,6 +24,8 @@ AliJetContainer::AliJetContainer():
   fJetAcceptanceType(kUser),
   fJetRadius(0),
   fRhoName(),
+  fLocalRhoName(),
+  fFlavourSelection(0),
   fPtBiasJetTrack(0),
   fPtBiasJetClus(0),
   fJetPtCut(1),
@@ -33,6 +37,10 @@ AliJetContainer::AliJetContainer():
   fJetMaxPhi(10),
   fMaxClusterPt(1000),
   fMaxTrackPt(100),
+  fZLeadingEmcCut(10.),
+  fZLeadingChCut(10.),
+  fNEFMinCut(-10.),
+  fNEFMaxCut(10.),
   fLeadingHadronType(0),
   fNLeadingJets(1),
   fJetBitMap(0),
@@ -40,6 +48,7 @@ AliJetContainer::AliJetContainer():
   fParticleContainer(0),
   fClusterContainer(0),
   fRho(0),
+  fLocalRho(0),
   fGeom(0),
   fRunNumber(0)
 {
@@ -54,6 +63,8 @@ AliJetContainer::AliJetContainer(const char *name):
   fJetAcceptanceType(kUser),
   fJetRadius(0),
   fRhoName(),
+  fLocalRhoName(),
+  fFlavourSelection(0),
   fPtBiasJetTrack(0),
   fPtBiasJetClus(0),
   fJetPtCut(1),
@@ -65,6 +76,10 @@ AliJetContainer::AliJetContainer(const char *name):
   fJetMaxPhi(10),
   fMaxClusterPt(1000),
   fMaxTrackPt(100),
+  fZLeadingEmcCut(10.),
+  fZLeadingChCut(10.),
+  fNEFMinCut(-10.),
+  fNEFMaxCut(10.),
   fLeadingHadronType(0),
   fNLeadingJets(1),
   fJetBitMap(0),
@@ -72,6 +87,7 @@ AliJetContainer::AliJetContainer(const char *name):
   fParticleContainer(0),
   fClusterContainer(0),
   fRho(0),
+  fLocalRho(0),
   fGeom(0),
   fRunNumber(0)
 {
@@ -121,6 +137,20 @@ void AliJetContainer::LoadRho(AliVEvent *event)
   }
 }
 
+//________________________________________________________________________
+void AliJetContainer::LoadLocalRho(AliVEvent *event)
+{
+  // Load local rho
+
+  if (!fLocalRhoName.IsNull() && !fLocalRho) {
+    fLocalRho = dynamic_cast<AliLocalRhoParameter*>(event->FindListObject(fLocalRhoName));
+    if (!fLocalRho) {
+      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fLocalRhoName.Data()));
+      return;
+    }
+  }
+}
+
 //________________________________________________________________________
 AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt)
 {
@@ -136,7 +166,8 @@ AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt)
 
   if (option.Contains("rho")) {
     while ((jet = GetNextAcceptJet())) {
-      if (jet->Pt()-jet->Area()*fRho->GetVal() > jetMax->Pt()-jetMax->Area()*fRho->GetVal()) jetMax = jet;
+      if ( (jet->Pt()-jet->Area()*GetRhoVal()) > (jetMax->Pt()-jetMax->Area()*GetRhoVal()) )
+       jetMax = jet;
     }
   }
   else {
@@ -161,14 +192,6 @@ AliEmcalJet* AliJetContainer::GetJet(Int_t i) const {
 
 }
 
-
-//________________________________________________________________________
-Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
-  AliEmcalJet *jet = GetJet(i);
-
-  return jet->Pt() - fRho->GetVal()*jet->Area();
-}
-
 //________________________________________________________________________
 AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
 
@@ -180,6 +203,24 @@ AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
   return jet;
 }
 
+//________________________________________________________________________
+AliEmcalJet* AliJetContainer::GetJetWithLabel(Int_t lab) const {
+
+  //Get particle with label lab in array
+  
+  Int_t i = GetIndexFromLabel(lab);
+  return GetJet(i);
+}
+
+//________________________________________________________________________
+AliEmcalJet* AliJetContainer::GetAcceptJetWithLabel(Int_t lab) const {
+
+  //Get particle with label lab in array
+  
+  Int_t i = GetIndexFromLabel(lab);
+  return GetAcceptJet(i);
+}
+
 //________________________________________________________________________
 AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) {
 
@@ -197,6 +238,37 @@ AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) {
   return jet;
 }
 
+//________________________________________________________________________
+AliEmcalJet* AliJetContainer::GetNextJet(Int_t i) {
+
+  //Get next jet; if i >= 0 (re)start counter from i; return 0 if no jet could be found
+
+  if (i>=0) fCurrentID = i;
+
+  const Int_t njets = GetNEntries();
+  AliEmcalJet *jet = 0;
+  while (fCurrentID < njets && !jet) { 
+    jet = GetJet(fCurrentID);
+    fCurrentID++;
+  }
+
+  return jet;
+}
+
+//________________________________________________________________________
+Double_t AliJetContainer::GetJetPtCorr(Int_t i) const {
+  AliEmcalJet *jet = GetJet(i);
+
+  return jet->Pt() - fRho->GetVal()*jet->Area();
+}
+
+//________________________________________________________________________
+Double_t AliJetContainer::GetJetPtCorrLocal(Int_t i) const {
+  AliEmcalJet *jet = GetJet(i);
+
+  return jet->Pt() - fLocalRho->GetLocalVal(jet->Phi(), fJetRadius)*jet->Area();
+}
+
 //________________________________________________________________________
 void AliJetContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
 {
@@ -229,35 +301,49 @@ Bool_t AliJetContainer::AcceptBiasJet(AliEmcalJet *jet) const
 //________________________________________________________________________
 Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
 {   
-  // Return true if jet is accepted.
-  if (!jet)
-    return kFALSE;
 
-  if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
-    return kFALSE;
-
-  if (jet->Pt() <= fJetPtCut) 
-    return kFALSE;
-  if (jet->Area() <= fJetAreaCut) 
-    return kFALSE;
-
-  if (jet->AreaEmc()<fAreaEmcCut)
-    return kFALSE;
-
-  if (!AcceptBiasJet(jet))
-    return kFALSE;
-  
-  if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
-    return kFALSE;
-
-  Double_t jetPhi = jet->Phi();
-  Double_t jetEta = jet->Eta();
-  
-  if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
-    jetPhi -= TMath::Pi() * 2;
-
-  return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
+   // Return true if jet is accepted.
+
+   if (!jet)
+      return kFALSE;
+
+   if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
+      return kFALSE;
+
+   if (jet->Pt() <= fJetPtCut) 
+      return kFALSE;
+
+   if (jet->Area() <= fJetAreaCut) 
+      return kFALSE;
+
+   if (jet->AreaEmc() < fAreaEmcCut)
+      return kFALSE;
+   
+   if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut)
+      return kFALSE;
+   
+   if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut)
+      return kFALSE;
+
+   if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut)
+      return kFALSE;
+   
+   if (!AcceptBiasJet(jet))
+      return kFALSE;
+   
+   if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
+      return kFALSE;
+   
+   if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection))
+      return kFALSE;
+   
+   Double_t jetPhi = jet->Phi();
+   Double_t jetEta = jet->Eta();
+   
+   if (fJetMinPhi < 0) // if limits are given in (-pi, pi) range
+      jetPhi -= TMath::Pi() * 2;
+   
+   return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
 }
 
 //________________________________________________________________________
@@ -309,6 +395,61 @@ void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet
     mom.SetPtEtaPhiM(maxClusterPt,maxClusterEta,maxClusterPhi,0.139);
 }
 
+//________________________________________________________________________
+Double_t AliJetContainer::GetZLeadingEmc(AliEmcalJet *jet) const
+{
+
+  if (fClusterContainer && fClusterContainer->GetArray()) {
+    TLorentzVector mom;
+    
+    AliVCluster *cluster = jet->GetLeadingCluster(fClusterContainer->GetArray());
+    if (cluster) {
+      cluster->GetMomentum(mom, const_cast<Double_t*>(fVertex));
+      
+      return GetZ(jet,mom);
+    }
+    else
+      return -1;
+  }
+  else
+    return -1;
+}
+
+//________________________________________________________________________
+Double_t AliJetContainer::GetZLeadingCharged(AliEmcalJet *jet) const
+{
+
+  if (fParticleContainer && fParticleContainer->GetArray() ) {
+    TLorentzVector mom;
+    
+    AliVParticle *track = jet->GetLeadingTrack(fParticleContainer->GetArray());
+    if (track) {
+      mom.SetPtEtaPhiM(track->Pt(),track->Eta(),track->Phi(),0.139);
+      
+      return GetZ(jet,mom);
+    }
+    else
+      return -1;
+  }
+  else
+    return -1;
+}
+
+//________________________________________________________________________
+Double_t AliJetContainer::GetZ(AliEmcalJet *jet, TLorentzVector mom) const
+{
+
+  Double_t pJetSq = jet->Px()*jet->Px() + jet->Py()*jet->Py() + jet->Pz()*jet->Pz();
+
+  if(pJetSq>1e-6)
+    return (mom.Px()*jet->Px() + mom.Py()*jet->Py() + mom.Pz()*jet->Pz())/pJetSq;
+  else {
+    AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
+    return -1;
+  }
+
+}
+
 //________________________________________________________________________
 void AliJetContainer::SetJetEtaPhiEMCAL()
 {
@@ -327,7 +468,7 @@ void AliJetContainer::SetJetEtaPhiEMCAL()
   else {
     AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
     SetJetEtaLimits(-0.7+fJetRadius,0.7-fJetRadius);
-    SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
+    SetJetPhiLimits(1.405+fJetRadius,3.135-fJetRadius);
   }
 }
 
@@ -357,7 +498,8 @@ void AliJetContainer::ResetCuts()
   fMaxClusterPt = 1000;
   fMaxTrackPt = 100;
   fLeadingHadronType = 0;
-
+  fZLeadingEmcCut = 10.;
+  fZLeadingChCut = 10.;
 }
 
 //________________________________________________________________________