setter to assume pion mass for clusters

[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetContainer.cxx

diff --git a/PWGJE/EMCALJetTasks/AliJetContainer.cxx b/PWGJE/EMCALJetTasks/AliJetContainer.cxx
index 010bed53eb39dad84bcc681cfda20006d60b2f73..c753fd09ee144121a731580960f5fe5b0c1a1971 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliJetContainer.cxx
+++ b/PWGJE/EMCALJetTasks/AliJetContainer.cxx
@@ -25,6 +25,7 @@ AliJetContainer::AliJetContainer():
   fJetRadius(0),
   fRhoName(),
   fLocalRhoName(),
+  fRhoMassName(),
   fFlavourSelection(0),
   fPtBiasJetTrack(0),
   fPtBiasJetClus(0),
@@ -45,10 +46,12 @@ AliJetContainer::AliJetContainer():
   fNLeadingJets(1),
   fJetBitMap(0),
   fJetTrigger(0),
+  fTagStatus(-1),
   fParticleContainer(0),
   fClusterContainer(0),
   fRho(0),
   fLocalRho(0),
+  fRhoMass(0),
   fGeom(0),
   fRunNumber(0)
 {
@@ -64,6 +67,7 @@ AliJetContainer::AliJetContainer(const char *name):
   fJetRadius(0),
   fRhoName(),
   fLocalRhoName(),
+  fRhoMassName(),
   fFlavourSelection(0),
   fPtBiasJetTrack(0),
   fPtBiasJetClus(0),
@@ -84,10 +88,12 @@ AliJetContainer::AliJetContainer(const char *name):
   fNLeadingJets(1),
   fJetBitMap(0),
   fJetTrigger(0),
+  fTagStatus(-1),
   fParticleContainer(0),
   fClusterContainer(0),
   fRho(0),
   fLocalRho(0),
+  fRhoMass(0),
   fGeom(0),
   fRunNumber(0)
 {
@@ -151,6 +157,21 @@ void AliJetContainer::LoadLocalRho(AliVEvent *event)
   }
 }
 
+//________________________________________________________________________
+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)
 {
@@ -166,7 +187,8 @@ AliEmcalJet* AliJetContainer::GetLeadingJet(const char* opt)
 
   if (option.Contains("rho")) {
     while ((jet = GetNextAcceptJet())) {
-      if ( (jet->Pt()-jet->Area()*GetRhoVal()) > (jetMax->Pt()-jetMax->Area()*GetRhoVal()) ) jetMax = jet;
+      if ( (jet->Pt()-jet->Area()*GetRhoVal()) > (jetMax->Pt()-jetMax->Area()*GetRhoVal()) )
+       jetMax = jet;
     }
   }
   else {
@@ -303,38 +325,65 @@ Bool_t AliJetContainer::AcceptJet(AliEmcalJet *jet) const
 
    // Return true if jet is accepted.
 
-   if (!jet)
-      return kFALSE;
+  if (!jet) {
+    AliDebug(11,"No jet found");
+    return kFALSE;
+  }
 
-   if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap)
-      return kFALSE;
+  if (jet->TestBits(fJetBitMap) != (Int_t)fJetBitMap) {
+    AliDebug(11,"Cut rejecting jet: Bit map");
+    return kFALSE;
+  }
 
-   if (jet->Pt() <= fJetPtCut) 
-      return kFALSE;
+  if (jet->Pt() <= fJetPtCut) {
+    AliDebug(11,"Cut rejecting jet: JetPtCut");
+    return kFALSE;
+  }
 
-   if (jet->Area() <= fJetAreaCut) 
-      return kFALSE;
+  if (jet->Area() <= fJetAreaCut)  {
+    AliDebug(11,"Cut rejecting jet: Area");
+    return kFALSE;
+  }
 
-   if (jet->AreaEmc() < fAreaEmcCut)
-      return kFALSE;
+  if (jet->AreaEmc() < fAreaEmcCut) {
+    AliDebug(11,"Cut rejecting jet: AreaEmc");
+    return kFALSE;
+  }
    
-   if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut)
-      return kFALSE;
+  if (fZLeadingChCut < 1 && GetZLeadingCharged(jet) > fZLeadingChCut) {
+    AliDebug(11,"Cut rejecting jet: ZLeading");
+    return kFALSE;
+  }
    
-   if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut)
-      return kFALSE;
+  if (fZLeadingEmcCut < 1 && GetZLeadingEmc(jet) > fZLeadingEmcCut) {
+    AliDebug(11,"Cut rejecting jet: ZLeadEmc");
+    return kFALSE;
+  }
 
-   if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut)
-      return kFALSE;
+  if (jet->NEF() < fNEFMinCut || jet->NEF() > fNEFMaxCut) {
+    AliDebug(11,"Cut rejecting jet: NEF");
+    return kFALSE;
+  }
    
-   if (!AcceptBiasJet(jet))
-      return kFALSE;
+  if (!AcceptBiasJet(jet)) {
+    AliDebug(11,"Cut rejecting jet: Bias");
+    return kFALSE;
+  }
+
+  if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt) {
+    AliDebug(11,"Cut rejecting jet: MaxTrClPt");
+    return kFALSE;
+  }
    
-   if (jet->MaxTrackPt() > fMaxTrackPt || jet->MaxClusterPt() > fMaxClusterPt)
-      return kFALSE;
+  if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection)) {
+    AliDebug(11,"Cut rejecting jet: Flavour");
+    return kFALSE;
+  }
    
-   if (fFlavourSelection != 0 && !jet->TestFlavourTag(fFlavourSelection))
-      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();
@@ -459,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);
 
@@ -480,6 +529,28 @@ void AliJetContainer::SetJetEtaPhiTPC()
   SetJetPhiLimits(-10, 10);
 }
 
+//________________________________________________________________________
+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() 
 {
@@ -510,3 +581,42 @@ void AliJetContainer::SetClassName(const char *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;
+}
+