setter to assume pion mass for clusters
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetContainer.cxx
index 4b3e381..c753fd0 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,9 @@ AliJetContainer::AliJetContainer():
   fJetAcceptanceType(kUser),
   fJetRadius(0),
   fRhoName(),
+  fLocalRhoName(),
+  fRhoMassName(),
+  fFlavourSelection(0),
   fPtBiasJetTrack(0),
   fPtBiasJetClus(0),
   fJetPtCut(1),
@@ -33,18 +38,26 @@ AliJetContainer::AliJetContainer():
   fJetMaxPhi(10),
   fMaxClusterPt(1000),
   fMaxTrackPt(100),
+  fZLeadingEmcCut(10.),
+  fZLeadingChCut(10.),
+  fNEFMinCut(-10.),
+  fNEFMaxCut(10.),
   fLeadingHadronType(0),
   fNLeadingJets(1),
   fJetBitMap(0),
   fJetTrigger(0),
+  fTagStatus(-1),
   fParticleContainer(0),
   fClusterContainer(0),
   fRho(0),
+  fLocalRho(0),
+  fRhoMass(0),
   fGeom(0),
   fRunNumber(0)
 {
   // Default constructor.
 
+  fClassName = "AliEmcalJet";
 }
 
 //________________________________________________________________________
@@ -53,6 +66,9 @@ AliJetContainer::AliJetContainer(const char *name):
   fJetAcceptanceType(kUser),
   fJetRadius(0),
   fRhoName(),
+  fLocalRhoName(),
+  fRhoMassName(),
+  fFlavourSelection(0),
   fPtBiasJetTrack(0),
   fPtBiasJetClus(0),
   fJetPtCut(1),
@@ -64,26 +80,34 @@ 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),
   fJetTrigger(0),
+  fTagStatus(-1),
   fParticleContainer(0),
   fClusterContainer(0),
   fRho(0),
+  fLocalRho(0),
+  fRhoMass(0),
   fGeom(0),
   fRunNumber(0)
 {
   // Standard constructor.
 
+  fClassName = "AliEmcalJet";
 }
 
 //________________________________________________________________________
-void AliJetContainer::SetJetArray(AliVEvent *event) 
+void AliJetContainer::SetArray(AliVEvent *event) 
 {
   // Set jet array
 
-  SetArray(event, "AliEmcalJet");
+  AliEmcalContainer::SetArray(event);
 
   if(fJetAcceptanceType==kTPC) {
     AliDebug(2,Form("%s: set TPC acceptance cuts",GetName()));
@@ -120,19 +144,51 @@ void AliJetContainer::LoadRho(AliVEvent *event)
 }
 
 //________________________________________________________________________
-AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt) const
+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;
+    }
+  }
+}
+
+//________________________________________________________________________
+void AliJetContainer::LoadRhoMass(AliVEvent *event)
+{
+  // Load rho
+
+  if (!fRhoMassName.IsNull() && !fRhoMass) {
+    fRhoMass = dynamic_cast<AliRhoParameter*>(event->FindListObject(fRhoMassName));
+    if (!fRhoMass) {
+      AliError(Form("%s: Could not retrieve rho_mass %s!", GetName(), fRhoMassName.Data()));
+      return;
+    }
+  }
+}
+
+
+//________________________________________________________________________
+AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt)
 {
   // Get the leading jet; if opt contains "rho" the sorting is according to pt-A*rho
 
   TString option(opt);
   option.ToLower();
 
+  Int_t tempID = fCurrentID;
+
   AliEmcalJet *jetMax = GetNextAcceptJet(0);
   AliEmcalJet *jet = 0;
 
   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 {
@@ -141,6 +197,8 @@ AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt) const
     }
   }
 
+  fCurrentID = tempID;
+
   return jetMax;
 }
 
@@ -155,14 +213,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 {
 
@@ -175,24 +225,72 @@ AliEmcalJet* AliJetContainer::GetAcceptJet(Int_t i) const {
 }
 
 //________________________________________________________________________
-AliEmcalJet* AliJetContainer::GetNextAcceptJet(Int_t i) const {
+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) {
 
   //Get next accepted jet; if i >= 0 (re)start counter from i; return 0 if no accepted jet could be found
 
-  static Int_t counter = -1;
-  if (i>=0) counter = i;
+  if (i>=0) fCurrentID = i;
 
   const Int_t njets = GetNEntries();
   AliEmcalJet *jet = 0;
-  while (counter < njets && !jet) { 
-    jet = GetAcceptJet(counter);
-    counter++;
+  while (fCurrentID < njets && !jet) { 
+    jet = GetAcceptJet(fCurrentID);
+    fCurrentID++;
   }
 
   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
 {
   //Get momentum of the i^th jet in array
@@ -224,35 +322,76 @@ 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 true if jet is accepted.
+
+  if (!jet) {
+    AliDebug(11,"No jet found");
     return kFALSE;
+  }
 
-  if (jet->Pt() <= fJetPtCut) 
+  if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap) {
+    AliDebug(11,"Cut rejecting jet: Bit map");
     return kFALSE;
-  if (jet->Area() <= fJetAreaCut) 
+  }
+
+  if (jet->Pt() <= fJetPtCut) {
+    AliDebug(11,"Cut rejecting jet: JetPtCut");
     return kFALSE;
+  }
 
-  if (jet->AreaEmc()<fAreaEmcCut)
+  if (jet->Area() <= fJetAreaCut)  {
+    AliDebug(11,"Cut rejecting jet: Area");
     return kFALSE;
+  }
 
-  if (!AcceptBiasJet(jet))
+  if (jet->AreaEmc() < fAreaEmcCut) {
+    AliDebug(11,"Cut rejecting jet: AreaEmc");
     return kFALSE;
-  
-  if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
+  }
+   
+  if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut) {
+    AliDebug(11,"Cut rejecting jet: ZLeading");
     return kFALSE;
+  }
+   
+  if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut) {
+    AliDebug(11,"Cut rejecting jet: ZLeadEmc");
+    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;
+  if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut) {
+    AliDebug(11,"Cut rejecting jet: NEF");
+    return kFALSE;
+  }
+   
+  if (!AcceptBiasJet(jet)) {
+    AliDebug(11,"Cut rejecting jet: Bias");
+    return kFALSE;
+  }
 
-  return (Bool_t)(jetEta > fJetMinEta && jetEta < fJetMaxEta && jetPhi > fJetMinPhi && jetPhi < fJetMaxPhi);
+  if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt) {
+    AliDebug(11,"Cut rejecting jet: MaxTrClPt");
+    return kFALSE;
+  }
+   
+  if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection)) {
+    AliDebug(11,"Cut rejecting jet: Flavour");
+    return kFALSE;
+  }
+   
+  if(fTagStatus>-1 && jet->GetTagStatus()!=fTagStatus) {
+    AliDebug(11,"Cut rejecting jet: tag status");
+    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);
 }
 
 //________________________________________________________________________
@@ -305,6 +444,61 @@ void AliJetContainer::GetLeadingHadronMomentum(TLorentzVector &mom, AliEmcalJet
 }
 
 //________________________________________________________________________
+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()
 {
   //Set default cuts for full jets
@@ -314,7 +508,7 @@ void AliJetContainer::SetJetEtaPhiEMCAL()
     SetJetEtaLimits(fGeom->GetArm1EtaMin() + fJetRadius, fGeom->GetArm1EtaMax() - fJetRadius);
 
     if(fRunNumber>=177295 && fRunNumber<=197470) //small SM masked in 2012 and 2013
-      SetJetPhiLimits(1.4+fJetRadius,TMath::Pi()-fJetRadius);
+      SetJetPhiLimits(1.405+fJetRadius,3.135-fJetRadius);
     else
       SetJetPhiLimits(fGeom->GetArm1PhiMin() * TMath::DegToRad() + fJetRadius, fGeom->GetArm1PhiMax() * TMath::DegToRad() - fJetRadius);
 
@@ -322,7 +516,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);
   }
 }
 
@@ -336,6 +530,28 @@ void AliJetContainer::SetJetEtaPhiTPC()
 }
 
 //________________________________________________________________________
+void AliJetContainer::PrintCuts() 
+{
+  // Reset cuts to default values
+
+  Printf("PtBiasJetTrack: %f",fPtBiasJetTrack);
+  Printf("PtBiasJetClus: %f",fPtBiasJetClus);
+  Printf("JetPtCut: %f", fJetPtCut);
+  Printf("JetAreaCut: %f",fJetAreaCut);
+  Printf("AreaEmcCut: %f",fAreaEmcCut);
+  Printf("JetMinEta: %f", fJetMinEta);
+  Printf("JetMaxEta: %f", fJetMaxEta);
+  Printf("JetMinPhi: %f", fJetMinPhi);
+  Printf("JetMaxPhi: %f", fJetMaxPhi);
+  Printf("MaxClusterPt: %f",fMaxClusterPt);
+  Printf("MaxTrackPt: %f",fMaxTrackPt);
+  Printf("LeadingHadronType: %d",fLeadingHadronType);
+  Printf("ZLeadingEmcCut: %f",fZLeadingEmcCut);
+  Printf("ZLeadingChCut: %f",fZLeadingChCut);
+
+}
+
+//________________________________________________________________________
 void AliJetContainer::ResetCuts() 
 {
   // Reset cuts to default values
@@ -352,5 +568,55 @@ void AliJetContainer::ResetCuts()
   fMaxClusterPt = 1000;
   fMaxTrackPt = 100;
   fLeadingHadronType = 0;
+  fZLeadingEmcCut = 10.;
+  fZLeadingChCut = 10.;
+}
 
+//________________________________________________________________________
+void AliJetContainer::SetClassName(const char *clname)
+{
+  // Set the class name
+
+  TClass cls(clname);
+  if (cls.InheritsFrom("AliEmcalJet")) fClassName = clname;
+  else AliError(Form("Unable to set class name %s for a AliJetContainer, it must inherits from AliEmcalJet!",clname));
 }
+
+//________________________________________________________________________
+Double_t AliJetContainer::GetFractionSharedPt(AliEmcalJet *jet1) const
+{
+  //
+  // Get fraction of shared pT between matched full and charged jet
+  // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
+  // Only works if tracks array of both jets is the same
+  //
+
+  AliEmcalJet *jet2 = jet1->ClosestJet();
+  if(!jet2) return -1;
+
+  Double_t fraction = 0.;
+  Double_t jetPt2 = jet2->Pt();
+  if(jetPt2>0) {
+    Double_t sumPt = 0.;
+    AliVParticle *vpf = 0x0;
+    Int_t iFound = 0;
+    for(Int_t icc=0; icc<jet2->GetNumberOfTracks(); icc++) {
+      Int_t idx = (Int_t)jet2->TrackAt(icc);
+      iFound = 0;
+      for(Int_t icf=0; icf<jet1->GetNumberOfTracks(); icf++) {
+       if(idx == jet1->TrackAt(icf) && iFound==0 ) {
+         iFound=1;
+         vpf = static_cast<AliVParticle*>(jet1->TrackAt(icf, fParticleContainer->GetArray()));
+         if(vpf) sumPt += vpf->Pt();
+         continue;
+       }
+      }
+    }
+    fraction = sumPt/jetPt2;
+  } else 
+    fraction = -1;
+  
+  return fraction;
+}
+