Adding FJ3 functionality:

    - AliFJWrapper: interface to generic and constituent subtractor
    - AliEmcalJetTask: flags to execute generic and constituent subtractor if requested + writing info to event and/or AliEmcalJet
    - AliEmcalJet: variables for generic subtraction method (derivatives and subtracted jet shape) -> only implemented for jet mass
    - AliJetShape: jet shape object inheriting from FunctionFromPseudoJet -> only implemented jet mass (AliJetShapeMass)

add FJ3 subtractors from contrib functionality

diff --git a/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg b/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg
index 31780b936b2b53f83559fae58ffd8fd38a4d0cee..2b72fd0b140a32a39189471a6ef13d4ed2160755 100644 (file)
--- a/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg
+++ b/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg
@@ -70,6 +70,8 @@ set ( SRCS
  EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTriggerQA.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalTriggerInfoQA.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskHJetEmbed.cxx
+ EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeDeriv.cxx
+ EMCALJetTasks/UserTasks/AliAnalysisTaskJetShapeConst.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskJetMatching.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskJetV2.cxx
  EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMass.cxx

diff --git a/PWGJE/EMCALJetTasks/AliEmcalJet.cxx b/PWGJE/EMCALJetTasks/AliEmcalJet.cxx
index f99fa83c3e9e38d3b609efe08dd1c676bbca6b36..75d97c2dc18cc601148361f3bc880a24dea4ee72 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliEmcalJet.cxx
+++ b/PWGJE/EMCALJetTasks/AliEmcalJet.cxx
@@ -39,7 +39,12 @@ AliEmcalJet::AliEmcalJet() :
   fTagStatus(-1),
   fPtSub(0),
   fPtVectSub(0),
-  fTriggers(0)
+  fTriggers(0),
+  fJetShapeMassFirstDer(0),
+  fJetShapeMassSecondDer(0),
+  fJetShapeMassFirstSub(0),
+  fJetShapeMassSecondSub(0),
+  fLabel(-1)
 {
   // Constructor.
 
@@ -78,7 +83,12 @@ AliEmcalJet::AliEmcalJet(Double_t px, Double_t py, Double_t pz) :
   fTagStatus(-1),
   fPtSub(0),
   fPtVectSub(0),
-  fTriggers(0)
+  fTriggers(0),
+  fJetShapeMassFirstDer(0),
+  fJetShapeMassSecondDer(0),
+  fJetShapeMassFirstSub(0),
+  fJetShapeMassSecondSub(0),
+  fLabel(-1)
 {    
   // Constructor.
 
@@ -123,7 +133,12 @@ AliEmcalJet::AliEmcalJet(Double_t pt, Double_t eta, Double_t phi, Double_t m) :
   fTagStatus(-1),
   fPtSub(0),
   fPtVectSub(0),
-  fTriggers(0)
+  fTriggers(0),
+  fJetShapeMassFirstDer(0),
+  fJetShapeMassSecondDer(0),
+  fJetShapeMassFirstSub(0),
+  fJetShapeMassSecondSub(0),
+  fLabel(-1)
 {
   // Constructor.
 
@@ -165,7 +180,12 @@ AliEmcalJet::AliEmcalJet(const AliEmcalJet &jet) :
   fTagStatus(jet.fTagStatus),
   fPtSub(jet.fPtSub),
   fPtVectSub(jet.fPtVectSub),
-  fTriggers(jet.fTriggers)
+  fTriggers(jet.fTriggers),
+  fJetShapeMassFirstDer(jet.fJetShapeMassFirstDer),
+  fJetShapeMassSecondDer(jet.fJetShapeMassSecondDer),
+  fJetShapeMassFirstSub(jet.fJetShapeMassFirstSub),
+  fJetShapeMassSecondSub(jet.fJetShapeMassSecondSub),
+  fLabel(jet.fLabel)
 {
   // Copy constructor.
 
@@ -212,6 +232,11 @@ AliEmcalJet &AliEmcalJet::operator=(const AliEmcalJet &jet)
     fPtSub              = jet.fPtSub;
     fPtVectSub          = jet.fPtVectSub;
     fTriggers           = jet.fTriggers;
+    fJetShapeMassFirstDer  = jet.fJetShapeMassFirstDer;
+    fJetShapeMassSecondDer = jet.fJetShapeMassSecondDer;
+    fJetShapeMassFirstSub  = jet.fJetShapeMassFirstSub;
+    fJetShapeMassSecondSub = jet.fJetShapeMassSecondSub;
+    fLabel              = jet.fLabel;
   }
 
   return *this;

diff --git a/PWGJE/EMCALJetTasks/AliEmcalJet.h b/PWGJE/EMCALJetTasks/AliEmcalJet.h
index 2a9672c48be06d2b43b27c26b2c2804a9b42e1c2..25f37daa4fbafb6031146c3aa4e3b05e78e5bccd 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliEmcalJet.h
+++ b/PWGJE/EMCALJetTasks/AliEmcalJet.h
@@ -47,7 +47,7 @@ class AliEmcalJet : public AliVParticle
   Double_t          Eta()                        const { return fEta;    }
   Double_t          Y()                          const { return 0.5*TMath::Log((E()+Pz())/(E()-Pz()));    }
   Short_t           Charge()                     const { return 0;       }
-  Int_t             GetLabel()                   const { return -1;      }
+  Int_t             GetLabel()                   const { return fLabel;  }
   Int_t             PdgCode()                    const { return 0;       }
   const Double_t   *PID()                        const { return 0;       }
   void              GetMom(TLorentzVector &vec)  const;
@@ -97,6 +97,8 @@ class AliEmcalJet : public AliVParticle
   void              Clear(Option_t */*option*/="")     { fClusterIDs.Set(0); fTrackIDs.Set(0); fClosestJets[0] = 0; fClosestJets[1] = 0; 
                                                          fClosestJetsDist[0] = 0; fClosestJetsDist[1] = 0; fMatched = 0; fPtSub = 0; }
   Double_t          DeltaR(const AliVParticle* part) const;
+
+  void              SetLabel(Int_t l)                  { fLabel = l;                       }
   void              SetArea(Double_t a)                { fArea    = a;                     }
   void              SetAreaEta(Double_t a)             { fAreaEta = a;                     }
   void              SetAreaPhi(Double_t a)             { fAreaPhi = a;                     }
@@ -142,6 +144,16 @@ class AliEmcalJet : public AliVParticle
   AliEmcalJet*      GetTaggedJet()                            const { return fTaggedJet                               ; }
   Int_t             GetTagStatus()                            const { return fTagStatus                               ; }
 
+  //jet shape derivatives
+  void              SetFirstDerivative(Double_t d)                  { fJetShapeMassFirstDer = d                       ; }
+  void              SetSecondDerivative(Double_t d)                 { fJetShapeMassSecondDer = d                      ; }
+  void              SetFirstOrderSubtracted(Double_t d)             { fJetShapeMassFirstSub = d                       ; }
+  void              SetSecondOrderSubtracted(Double_t d)            { fJetShapeMassSecondSub = d                      ; }
+  Double_t          GetFirstDerivative()                      const { return fJetShapeMassFirstDer                    ; }
+  Double_t          GetSecondDerivative()                     const { return fJetShapeMassSecondDer                   ; }
+  Double_t          GetFirstOrderSubtracted()                 const { return fJetShapeMassFirstSub                    ; }
+  Double_t          GetSecondOrderSubtracted()                const { return fJetShapeMassSecondSub                   ; }
+
  protected:
   Double32_t        fPt;                  //[0,0,12]   pt 
   Double32_t        fEta;                 //[-1,1,12]  eta
@@ -172,6 +184,11 @@ class AliEmcalJet : public AliVParticle
   Double_t          fPtSub;               //!          background subtracted pt (not stored set from outside) 
   Double_t          fPtVectSub;           //!          background vector subtracted pt (not stored set from outside) 
   UInt_t            fTriggers;            //!          triggers that the jet might have fired (AliVEvent::EOfflineTriggerTypes)
+  Double_t          fJetShapeMassFirstDer;  //!        result from shape derivatives for jet mass: 1st derivative
+  Double_t          fJetShapeMassSecondDer; //!        result from shape derivatives for jet mass: 2nd derivative
+  Double_t          fJetShapeMassFirstSub;  //!        result from shape derivatives for jet mass: 1st order subtracted
+  Double_t          fJetShapeMassSecondSub; //!        result from shape derivatives for jet mass: 2nd order subtracted
+  Int_t             fLabel;               //           label to inclusive jet for constituent subtracted jet    
 
   private:
     struct sort_descend
@@ -180,6 +197,6 @@ class AliEmcalJet : public AliVParticle
         bool operator () (const std::pair<Double_t, Int_t>& p1, const std::pair<Double_t, Int_t>& p2)  { return p1.first > p2.first ; }
         };
 
-  ClassDef(AliEmcalJet,13) // Emcal jet class in cylindrical coordinates
+  ClassDef(AliEmcalJet,14) // Emcal jet class in cylindrical coordinates
 };
 #endif

diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx b/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
index 98ef9b9327c2734535ae90197041e8bcad55cef0..a5df58bfae67dab5ce1f753ffc761a29ae079cd2 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
+++ b/PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
@@ -1,7 +1,7 @@
 //
 // Emcal jet finder task.
 //
-// Authors: C.Loizides, S.Aiola
+// Authors: C.Loizides, S.Aiola, M. Verweij
 
 #include <vector>
 #include "AliEmcalJetTask.h"
@@ -24,6 +24,7 @@
 #include "AliVCluster.h"
 #include "AliVEvent.h"
 #include "AliVParticle.h"
+#include "AliRhoParameter.h"
 using std::cout;
 using std::endl;
 using std::cerr;
@@ -36,6 +37,7 @@ AliEmcalJetTask::AliEmcalJetTask() :
   fTracksName("Tracks"),
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
+  fJetsSubName(""),
   fJetType(kNone),
   fConstSel(0),
   fMCConstSel(0),
@@ -63,10 +65,20 @@ AliEmcalJetTask::AliEmcalJetTask() :
   fIsEmcPart(0),
   fLegacyMode(kFALSE),
   fCodeDebug(kFALSE),
+  fDoGenericSubtraction(kFALSE),
+  fDoConstituentSubtraction(kFALSE),
+  fUseExternalBkg(kFALSE),
+  fRhoName(""),
+  fRhomName(""),
+  fRho(0),
+  fRhom(0),
   fJets(0),
+  fJetsSub(0),
   fEvent(0),
   fTracks(0),
-  fClus(0)
+  fClus(0),
+  fRhoParam(0),
+  fRhomParam(0)
 {
   // Default constructor.
 }
@@ -77,6 +89,7 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) :
   fTracksName("Tracks"),
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
+  fJetsSubName(""),
   fJetType(kAKT|kFullJet|kRX1Jet),
   fConstSel(0),
   fMCConstSel(0),
@@ -104,10 +117,20 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) :
   fIsEmcPart(0),
   fLegacyMode(kFALSE),
   fCodeDebug(kFALSE),
+  fDoGenericSubtraction(kFALSE),
+  fDoConstituentSubtraction(kFALSE),
+  fUseExternalBkg(kFALSE),
+  fRhoName(""),
+  fRhomName(""),
+  fRho(0),
+  fRhom(0),
   fJets(0),
+  fJetsSub(0),
   fEvent(0),
   fTracks(0),
-  fClus(0)
+  fClus(0),
+  fRhoParam(0),
+  fRhomParam(0)
 {
   // Standard constructor.
 
@@ -127,6 +150,11 @@ void AliEmcalJetTask::UserCreateOutputObjects()
 
   fJets = new TClonesArray("AliEmcalJet");
   fJets->SetName(fJetsName);
+
+  if(!fJetsSubName.IsNull()) {
+    fJetsSub = new TClonesArray("AliEmcalJet");
+    fJetsSub->SetName(fJetsSubName);
+  }
 }
 
 //________________________________________________________________________
@@ -141,6 +169,7 @@ void AliEmcalJetTask::UserExec(Option_t *)
 
   // clear the jet array (normally a null operation)
   fJets->Delete();
+  if(fJetsSub) fJetsSub->Delete();
 
   FindJets();
 }
@@ -160,6 +189,11 @@ void AliEmcalJetTask::FindJets()
     return;
   }
 
+ if (fRhoParam)
+    fRho = fRhoParam->GetVal();
+  if (fRhomParam)
+    fRhom = fRhomParam->GetVal();
+
   TString name("kt");
   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
   if ((fJetType & kAKT) != 0) {
@@ -345,6 +379,18 @@ void AliEmcalJetTask::FindJets()
   // run jet finder
   fjw.Run();
 
+  //run generic subtractor
+  if(fDoGenericSubtraction) {
+    fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+    fjw.DoGenericSubtractionJetMass();
+  }
+
+  //run constituent subtractor
+  if(fDoConstituentSubtraction) {
+    fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+    fjw.DoConstituentSubtraction();
+  }
+
   // get geometry
   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
   if (!geom) {
@@ -373,6 +419,21 @@ void AliEmcalJetTask::FindJets()
 
     AliEmcalJet *jet = new ((*fJets)[jetCount]) 
       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
+    jet->SetLabel(ij);
+
+    //do generic subtraction if requested
+    if(fDoGenericSubtraction) {
+#ifdef FASTJET_VERSION
+      std::vector<fastjet::contrib::GenericSubtractorInfo> jetMassInfo = fjw.GetGenSubtractorInfoJetMass();
+      UInt_t n = (UInt_t)jetMassInfo.size();
+      if(n>ij && n>0) {
+       jet->SetFirstDerivative(jetMassInfo[ij].first_derivative());
+       jet->SetSecondDerivative(jetMassInfo[ij].second_derivative());
+       jet->SetFirstOrderSubtracted(jetMassInfo[ij].first_order_subtracted());
+       jet->SetSecondOrderSubtracted(jetMassInfo[ij].second_order_subtracted());
+      }
+#endif
+    }
 
     // loop over constituents
     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
@@ -530,6 +591,30 @@ void AliEmcalJetTask::FindJets()
     jetCount++;
   }
   //fJets->Sort();
+
+  //run constituent subtractor if requested
+  if(fDoConstituentSubtraction) {
+#ifdef FASTJET_VERSION
+    if(!fJetsSub) AliWarning(Form("No jet branch to write to for subtracted jets. fJetsSubName: %s",fJetsSubName.Data()));
+    else {
+      std::vector<fastjet::PseudoJet> jets_sub;
+      jets_sub = fjw.GetConstituentSubtrJets();
+      AliDebug(1,Form("%d constituent subtracted jets found", (Int_t)jets_sub.size()));
+      for (UInt_t ijet=0, jetCount=0; ijet<jets_sub.size(); ++ijet) {
+       //Only storing 4-vector and jet area of unsubtracted jet        
+       AliEmcalJet *jet_sub = new ((*fJetsSub)[jetCount]) 
+         AliEmcalJet(jets_sub[ijet].perp(), jets_sub[ijet].eta(), jets_sub[ijet].phi(), jets_sub[ijet].m());
+       jet_sub->SetLabel(ijet);
+       fastjet::PseudoJet area(fjw.GetJetAreaVector(ijet));
+       jet_sub->SetArea(area.perp());
+       jet_sub->SetAreaEta(area.eta());
+       jet_sub->SetAreaPhi(area.phi());
+       jet_sub->SetAreaEmc(area.perp());
+       jetCount++;
+      }
+    }
+#endif
+  } //constituent subtraction
 }
 
 //________________________________________________________________________
@@ -580,6 +665,9 @@ Bool_t AliEmcalJetTask::DoInit()
     return 0;
   }
 
+  if (!(fEvent->FindListObject(fJetsSubName)) && fJetsSub)
+    fEvent->AddObject(fJetsSub);
+
   if (fTracksName == "Tracks")
     am->LoadBranch("Tracks");
   if (!fTracks && !fTracksName.IsNull()) {
@@ -610,5 +698,20 @@ Bool_t AliEmcalJetTask::DoInit()
       fIsEmcPart = 1;
   }
 
+  if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event
+    fRhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
+    if (!fRhoParam) {
+      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
+      return 0;
+    }
+  }
+  if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event
+    fRhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhomName));
+    if (!fRhomParam) {
+      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhomName.Data()));
+      return 0;
+    }
+  }
+  
   return 1;
 }

diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetTask.h b/PWGJE/EMCALJetTasks/AliEmcalJetTask.h
index 81fee038cf61d1b0903fe4496dc8624375ec93ef..3e2231b22feff5c1ef88db1245a829e1d5762f63 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliEmcalJetTask.h
+++ b/PWGJE/EMCALJetTasks/AliEmcalJetTask.h
@@ -3,6 +3,7 @@
 
 class TClonesArray;
 class AliVEvent;
+class AliRhoParameter;
 
 namespace fastjet {
   class PseudoJet;
@@ -10,6 +11,7 @@ namespace fastjet {
 
 #include "AliLog.h"
 #include "AliAnalysisTaskSE.h"
+#include "AliFJWrapper.h"
 
 class AliEmcalJetTask : public AliAnalysisTaskSE {
  public:
@@ -42,6 +44,7 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
   void                   SetAlgo(Int_t a)                 { if (a==0) fJetType |= kKT; else fJetType |= kAKT; }  // for backward compatibility only
   void                   SetClusName(const char *n)       { fCaloName      = n     ; }
   void                   SetJetsName(const char *n)       { fJetsName      = n     ; }
+  void                   SetJetsSubName(const char *n)    { fJetsSubName   = n     ; }
   void                   SetJetType(UInt_t t)             { fJetType       = t     ; }
   void                   SetMarkConstituents(UInt_t m)    { fMarkConst     = m     ; }
   void                   SetMinJetArea(Double_t a)        { fMinJetArea    = a     ; }
@@ -74,11 +77,18 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
       fOfflineTriggerMask = fOfflineTriggerMask | offlineTriggerMask;
     }
   }
-  void                   SetLegacyMode(Bool_t mode)       { fLegacyMode ^= mode; }
+  void                   SetLegacyMode(Bool_t mode)       { fLegacyMode = mode; }
   void                   SetCodeDebug(Bool_t val)         { fCodeDebug = val; }
 
+  void                   SetRhoName(const char *n)              { fRhoName      = n            ; }
+  void                   SetRhomName(const char *n)             { fRhomName     = n            ; }
+  void                   SetGenericSubtraction(Bool_t b)        { fDoGenericSubtraction     = b; }
+  void                   SetConstituentSubtraction(Bool_t b)    { fDoConstituentSubtraction = b; }
+  void                   SetUseExternalBkg(Bool_t b, Double_t rho, Double_t rhom) { fUseExternalBkg = b; fRho = rho; fRhom = rhom;}
+
   UInt_t                 GetJetType()                     { return fJetType; }
   Bool_t                 GetLegacyMode()                  { return fLegacyMode; }
+  TString                GetJetsSubName()                 { return fJetsSubName; }
 
   const char*            GetJetsName()                    { return fJetsName.Data(); }
   Double_t               GetRadius()                      { return fRadius; }
@@ -115,6 +125,7 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
   TString                fTracksName;             // name of track collection
   TString                fCaloName;               // name of calo cluster collection
   TString                fJetsName;               // name of jet collection
+  TString                fJetsSubName;            // name of subtracted jet collection
   UInt_t                 fJetType;                // jet type (algorithm, radius, constituents)
   UInt_t                 fConstSel;               // select constituents from a previous jet finding
   UInt_t                 fMCConstSel;             // select MC constituents (label!=0) from a previous jet finding
@@ -142,15 +153,27 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
   Bool_t                 fIsEmcPart;              //!=true if emcal particles are given as input (for clusters)
   Bool_t                 fLegacyMode;             //! if true, enable FJ 2.x behavior
   Bool_t                 fCodeDebug;              // use nontested code changes 
+
+  Bool_t                 fDoGenericSubtraction;   // calculate generic subtraction
+  Bool_t                 fDoConstituentSubtraction; // calculate constituent subtraction
+  Bool_t                 fUseExternalBkg;         // use external background for generic subtractor
+  TString                fRhoName;                // name of rho
+  TString                fRhomName;               // name of rhom
+  Double_t               fRho;                    // pT background density
+  Double_t               fRhom;                   // mT background density
+
   TClonesArray          *fJets;                   //!jet collection
+  TClonesArray          *fJetsSub;                //!subtracted jet collection
   AliVEvent             *fEvent;                  //!current event
   TClonesArray          *fTracks;                 //!tracks collection
   TClonesArray          *fClus;                   //!cluster collection
+  AliRhoParameter       *fRhoParam;               //!event rho
+  AliRhoParameter       *fRhomParam;              //!event rhom
 
  private:
   AliEmcalJetTask(const AliEmcalJetTask&);            // not implemented
   AliEmcalJetTask &operator=(const AliEmcalJetTask&); // not implemented
 
-  ClassDef(AliEmcalJetTask, 11) // Jet producing task
+  ClassDef(AliEmcalJetTask, 12) // Jet producing task
 };
 #endif

diff --git a/PWGJE/EMCALJetTasks/AliFJWrapper.h b/PWGJE/EMCALJetTasks/AliFJWrapper.h
index 586093b34253c576d2460a761f2e881a31fc8ea7..7141f35e394833e272cd63d1b80e739e38941c36 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliFJWrapper.h
+++ b/PWGJE/EMCALJetTasks/AliFJWrapper.h
@@ -9,6 +9,7 @@
 #include <TString.h>
 #include "AliLog.h"
 #include "FJ_includes.h"
+#include "AliJetShape.h"
 
 class AliFJWrapper
 {
@@ -23,20 +24,26 @@ class AliFJWrapper
   virtual void  Clear(const Option_t* /*opt*/ = "");
   virtual void  CopySettingsFrom (const AliFJWrapper& wrapper);
   virtual void  GetMedianAndSigma(Double_t& median, Double_t& sigma, Int_t remove = 0) const;
-  fastjet::ClusterSequenceArea*         GetClusterSequence() const { return fClustSeq;                   }
-  const std::vector<fastjet::PseudoJet> GetInputVectors()    const { return fInputVectors;               }
-  const std::vector<fastjet::PseudoJet> GetInclusiveJets()   const { return fInclusiveJets;              }
-  std::vector<fastjet::PseudoJet>       GetJetConstituents(UInt_t idx) const;
-  Double_t                              GetMedianUsedForBgSubtraction() const { return fMedUsedForBgSub; }
-  const char*                           GetName()            const { return fName;                       }
-  const char*                           GetTitle()           const { return fTitle;                      }
-  Double_t                              GetJetArea         (UInt_t idx) const;
-  fastjet::PseudoJet                    GetJetAreaVector   (UInt_t idx) const;
-  Double_t                              GetJetSubtractedPt (UInt_t idx) const;
-  virtual std::vector<double>           GetSubtractedJetsPts(Double_t median_pt = -1, Bool_t sorted = kFALSE);
-  Bool_t                                GetLegacyMode()            { return fLegacyMode; }
+  fastjet::ClusterSequenceArea*           GetClusterSequence() const { return fClustSeq;                   }
+  const std::vector<fastjet::PseudoJet>   GetInputVectors()    const { return fInputVectors;               }
+  const std::vector<fastjet::PseudoJet>   GetInclusiveJets()   const { return fInclusiveJets;              }
+  std::vector<fastjet::PseudoJet>         GetJetConstituents(UInt_t idx) const;
+  Double_t                                GetMedianUsedForBgSubtraction() const { return fMedUsedForBgSub; }
+  const char*                             GetName()            const { return fName;                       }
+  const char*                             GetTitle()           const { return fTitle;                      }
+  Double_t                                GetJetArea         (UInt_t idx) const;
+  fastjet::PseudoJet                      GetJetAreaVector   (UInt_t idx) const;
+  Double_t                                GetJetSubtractedPt (UInt_t idx) const;
+  virtual std::vector<double>             GetSubtractedJetsPts(Double_t median_pt = -1, Bool_t sorted = kFALSE);
+  Bool_t                                  GetLegacyMode()            { return fLegacyMode; }
+#ifdef FASTJET_VERSION
+  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetMass() const {return fGenSubtractorInfoJetMass;}
+  const std::vector<fastjet::PseudoJet>   GetConstituentSubtrJets() const { return fConstituentSubtrJets;  }
+#endif
 
   virtual Int_t Run();
+  virtual Int_t DoGenericSubtractionJetMass();
+  virtual Int_t DoConstituentSubtraction();
 
   void SetStrategy(const fastjet::Strategy &strat)                 { fStrategy = strat;  }
   void SetAlgorithm(const fastjet::JetAlgorithm &algor)            { fAlgor    = algor;  }
@@ -57,40 +64,47 @@ class AliFJWrapper
   void SetupStrategyfromOpt(const char *option);
   void SetLegacyMode (Bool_t mode)      { fLegacyMode ^= mode; }
   void SetLegacyFJ();
+  void SetUseExternalBkg(Bool_t b, Double_t rho, Double_t rhom) { fUseExternalBkg = b; fRho = rho; fRhom = rhom;}
 
  protected:
-  TString                                fName;             //!
-  TString                                fTitle;            //!
-  std::vector<fastjet::PseudoJet>        fInputVectors;     //!
-  std::vector<fastjet::PseudoJet>        fInclusiveJets;    //!
-  std::vector<double>                    fSubtractedJetsPt; //!
-  fastjet::AreaDefinition               *fAreaDef;          //!
-  fastjet::VoronoiAreaSpec              *fVorAreaSpec;      //!
-  fastjet::GhostedAreaSpec              *fGhostedAreaSpec;  //!
-  fastjet::JetDefinition                *fJetDef;           //!
-  fastjet::JetDefinition::Plugin        *fPlugin;           //!
-  fastjet::RangeDefinition              *fRange;            //!
-  fastjet::ClusterSequenceArea          *fClustSeq;         //!
-  fastjet::Strategy                      fStrategy;         //!
-  fastjet::JetAlgorithm                  fAlgor;            //!
-  fastjet::RecombinationScheme           fScheme;           //!
-  fastjet::AreaType                      fAreaType;         //!
-  Int_t                                  fNGhostRepeats;    //!
-  Double_t                               fGhostArea;       //!
-  Double_t                               fMaxRap;          //!
-  Double_t                               fR;                //!
+  TString                                fName;               //!
+  TString                                fTitle;              //!
+  std::vector<fastjet::PseudoJet>        fInputVectors;       //!
+  std::vector<fastjet::PseudoJet>        fInclusiveJets;      //!
+  std::vector<double>                    fSubtractedJetsPt;   //!
+  std::vector<fastjet::PseudoJet>        fConstituentSubtrJets; //!
+  fastjet::AreaDefinition               *fAreaDef;            //!
+  fastjet::VoronoiAreaSpec              *fVorAreaSpec;        //!
+  fastjet::GhostedAreaSpec              *fGhostedAreaSpec;    //!
+  fastjet::JetDefinition                *fJetDef;             //!
+  fastjet::JetDefinition::Plugin        *fPlugin;             //!
+  fastjet::RangeDefinition              *fRange;              //!
+  fastjet::ClusterSequenceArea          *fClustSeq;           //!
+  fastjet::Strategy                      fStrategy;           //!
+  fastjet::JetAlgorithm                  fAlgor;              //!
+  fastjet::RecombinationScheme           fScheme;             //!
+  fastjet::AreaType                      fAreaType;           //!
+  Int_t                                  fNGhostRepeats;      //!
+  Double_t                               fGhostArea;         //!
+  Double_t                               fMaxRap;            //!
+  Double_t                               fR;                  //!
   // no setters for the moment - used default values in the constructor
-  Double_t                               fGridScatter;      //!
-  Double_t                               fKtScatter;       //!
-  Double_t                               fMeanGhostKt;      //!
-  Int_t                                  fPluginAlgor;      //!
+  Double_t                               fGridScatter;        //!
+  Double_t                               fKtScatter;         //!
+  Double_t                               fMeanGhostKt;        //!
+  Int_t                                  fPluginAlgor;        //!
   // extra parameters
-  Double_t                               fMedUsedForBgSub;  //!
-  Bool_t                                 fUseArea4Vector;   //!
+  Double_t                               fMedUsedForBgSub;    //!
+  Bool_t                                 fUseArea4Vector;     //!
 #ifdef FASTJET_VERSION
-  fastjet::JetMedianBackgroundEstimator *fBkrdEstimator;    //!
+  fastjet::JetMedianBackgroundEstimator   *fBkrdEstimator;        //!
+  fastjet::contrib::GenericSubtractor     *fGenSubtractor;        //!
+  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetMass;    //!
 #endif
-  Bool_t                                 fLegacyMode;       //!
+  Bool_t                                   fLegacyMode;           //!
+  Bool_t                                   fUseExternalBkg;       //!
+  Double_t                                 fRho;                  //  pT background density
+  Double_t                                 fRhom;                 //  mT background density
 
   virtual void   SubtractBackground(const Double_t median_pt = -1);
 
@@ -119,6 +133,7 @@ AliFJWrapper::AliFJWrapper(const char *name, const char *title)
   , fInputVectors      ( )
   , fInclusiveJets     ( )
   , fSubtractedJetsPt  ( )
+  , fConstituentSubtrJets ( )
   , fAreaDef           (0)
   , fVorAreaSpec       (0)
   , fGhostedAreaSpec   (0)
@@ -142,8 +157,13 @@ AliFJWrapper::AliFJWrapper(const char *name, const char *title)
   , fUseArea4Vector    (kFALSE)
 #ifdef FASTJET_VERSION
   , fBkrdEstimator     (0)
+  , fGenSubtractor     (0)
+  , fGenSubtractorInfoJetMass ( )
 #endif
   , fLegacyMode        (false)
+  , fUseExternalBkg    (false)
+  , fRho               (0)
+  , fRhom              (0)
 {
   // Constructor.
 }
@@ -161,7 +181,8 @@ AliFJWrapper::~AliFJWrapper()
   delete fRange;
   delete fClustSeq;
 #ifdef FASTJET_VERSION
-  if (fBkrdEstimator) delete fBkrdEstimator;
+  if (fBkrdEstimator)     delete fBkrdEstimator;
+  if (fGenSubtractor)     delete fGenSubtractor;
 #endif
 }
 
@@ -172,20 +193,23 @@ void AliFJWrapper::CopySettingsFrom(const AliFJWrapper& wrapper)
   // You very often want to keep most of the settings 
   // but change only the algorithm or R - do it after call to this function
 
-  fStrategy       = wrapper.fStrategy;
-  fAlgor          = wrapper.fAlgor;
-  fScheme         = wrapper.fScheme;
-  fAreaType       = wrapper.fAreaType;
-  fNGhostRepeats  = wrapper.fNGhostRepeats;
-  fGhostArea      = wrapper.fGhostArea;
-  fMaxRap         = wrapper.fMaxRap;
-  fR              = wrapper.fR;
-  fGridScatter    = wrapper.fGridScatter;
-  fKtScatter      = wrapper.fKtScatter;
-  fMeanGhostKt    = wrapper.fMeanGhostKt;
-  fPluginAlgor    = wrapper.fPluginAlgor;
-  fUseArea4Vector = wrapper.fUseArea4Vector;
-  fLegacyMode     = wrapper.fLegacyMode;
+  fStrategy         = wrapper.fStrategy;
+  fAlgor            = wrapper.fAlgor;
+  fScheme           = wrapper.fScheme;
+  fAreaType         = wrapper.fAreaType;
+  fNGhostRepeats    = wrapper.fNGhostRepeats;
+  fGhostArea        = wrapper.fGhostArea;
+  fMaxRap           = wrapper.fMaxRap;
+  fR                = wrapper.fR;
+  fGridScatter      = wrapper.fGridScatter;
+  fKtScatter        = wrapper.fKtScatter;
+  fMeanGhostKt      = wrapper.fMeanGhostKt;
+  fPluginAlgor      = wrapper.fPluginAlgor;
+  fUseArea4Vector   = wrapper.fUseArea4Vector;
+  fLegacyMode       = wrapper.fLegacyMode;
+  fUseExternalBkg   = wrapper.fUseExternalBkg;
+  fRho              = wrapper.fRho;
+  fRhom             = wrapper.fRhom;
 }
 
 //_________________________________________________________________________________________________
@@ -207,7 +231,8 @@ void AliFJWrapper::Clear(const Option_t */*opt*/)
   delete fRange;           fRange           = 0;
   delete fClustSeq;        fClustSeq        = 0;
 #ifdef FASTJET_VERSION
-  if (fBkrdEstimator) delete fBkrdEstimator; fBkrdEstimator = 0;
+  if (fBkrdEstimator)     delete fBkrdEstimator     ;  fBkrdEstimator     = 0;
+  if (fGenSubtractor)     delete fGenSubtractor     ;  fGenSubtractor     = 0;
 #endif
 }
 
@@ -415,7 +440,7 @@ Int_t AliFJWrapper::Run()
 
   // FJ3 :: Define an JetMedianBackgroundEstimator just in case it will be used 
 #ifdef FASTJET_VERSION
-  //fBkrdEstimator = new fj::JetMedianBackgroundEstimator(*fRange, *fJetDef, *fAreaDef) ;
+  fBkrdEstimator     = new fj::JetMedianBackgroundEstimator(fj::SelectorAbsRapMax(fMaxRap));
 #endif
 
   if (fLegacyMode) { SetLegacyFJ(); } // for FJ 2.x even if fLegacyMode is set, SetLegacyFJ is dummy
@@ -499,6 +524,54 @@ void AliFJWrapper::SubtractBackground(Double_t median_pt)
   }
 }
 
+//_________________________________________________________________________________________________
+Int_t AliFJWrapper::DoGenericSubtractionJetMass() {
+  //Do generic subtraction for jet mass
+#ifdef FASTJET_VERSION
+  if(fUseExternalBkg)   fGenSubtractor     = new fj::contrib::GenericSubtractor(fRho,fRhom);
+  else                  fGenSubtractor     = new fj::contrib::GenericSubtractor(fBkrdEstimator);
+
+  // Define jet shape
+  AliJetShapeMass shapeMass;
+
+  // clear the generic subtractor info vector
+  fGenSubtractorInfoJetMass.clear();
+  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
+    fj::contrib::GenericSubtractorInfo info;
+    if(fInclusiveJets[i].perp()>0.)
+      double subtracted_shape = (*fGenSubtractor)(shapeMass, fInclusiveJets[i], info);
+    fGenSubtractorInfoJetMass.push_back(info);
+  }
+  
+#endif
+  return 0;
+}
+
+//_________________________________________________________________________________________________
+Int_t AliFJWrapper::DoConstituentSubtraction() {
+  //Do constituent subtraction
+#ifdef FASTJET_VERSION
+  fj::contrib::ConstituentSubtractor *subtractor;
+  if(fUseExternalBkg)  {
+    Printf("Using external background rho: %f  rhom: %f",fRho,fRhom);
+    subtractor     = new fj::contrib::ConstituentSubtractor(fRho,fRhom,kFALSE,kTRUE);
+  }
+  else                 subtractor     = new fj::contrib::ConstituentSubtractor(fBkrdEstimator);
+
+  //clear constituent subtracted jets
+  fConstituentSubtrJets.clear();
+  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
+    fj::PseudoJet subtracted_jet(0.,0.,0.,0.);
+    if(fInclusiveJets[i].perp()>0.)
+      subtracted_jet = (*subtractor)(fInclusiveJets[i]);
+    fConstituentSubtrJets.push_back(subtracted_jet);
+  }
+  if(subtractor) delete subtractor;
+
+#endif
+  return 0;
+}
+
 //_________________________________________________________________________________________________
 void AliFJWrapper::SetupAlgorithmfromOpt(const char *option)
 {

diff --git a/PWGJE/EMCALJetTasks/AliJetShape.cxx b/PWGJE/EMCALJetTasks/AliJetShape.cxx
new file mode 100644 (file)
index 0000000..feb6886
--- /dev/null
+++ b/PWGJE/EMCALJetTasks/AliJetShape.cxx
@@ -0,0 +1,6 @@
+// $Id$
+//
+// Author: M. Verweij
+
+#define AliJetShape_CXX
+#include "AliJetShape.h"

diff --git a/PWGJE/EMCALJetTasks/AliJetShape.h b/PWGJE/EMCALJetTasks/AliJetShape.h
new file mode 100644 (file)
index 0000000..18bbfa3
--- /dev/null
+++ b/PWGJE/EMCALJetTasks/AliJetShape.h
@@ -0,0 +1,21 @@
+#ifndef AliJetShape_H
+#define AliJetShape_H
+
+// $Id$
+
+#include <vector>
+#include <TString.h>
+#include "fastjet/PseudoJet.hh"
+#ifdef FASTJET_VERSION
+#include "fastjet/FunctionOfPseudoJet.hh"
+#endif
+
+#ifdef FASTJET_VERSION
+class AliJetShapeMass : public fastjet::FunctionOfPseudoJet<Double32_t>
+{
+ public:
+  virtual std::string description() const{return "jet mass";}
+  Double32_t result(const fastjet::PseudoJet &jet) const{ return jet.m();}
+};
+#endif
+#endif

diff --git a/PWGJE/EMCALJetTasks/FJ_includes.h b/PWGJE/EMCALJetTasks/FJ_includes.h
index 573c8bf0999c403e5e4afd6fe28b5c4cd9eb60ec..fb60b77753abb75a04b67e48e4b9cd9a42935b87 100644 (file)
--- a/PWGJE/EMCALJetTasks/FJ_includes.h
+++ b/PWGJE/EMCALJetTasks/FJ_includes.h
@@ -16,7 +16,16 @@
 #include <fastjet/SISConePlugin.hh>
 #include <fastjet/CDFMidPointPlugin.hh>
 #ifdef FASTJET_VERSION
+//if FJ_VERSION>2
+#include <fastjet/Selector.hh>
+#include <fastjet/FunctionOfPseudoJet.hh>
 #include <fastjet/tools/JetMedianBackgroundEstimator.hh>
+#include <fastjet/tools/BackgroundEstimatorBase.hh>
+#include <fastjet/tools/Subtractor.hh>
+#include <fastjet/contrib/GenericSubtractor.hh>
+#include <fastjet/contrib/ShapeWithComponents.hh>
+#include <fastjet/contrib/ConstituentSubtractor.hh>
 #endif
 #endif
+
 #endif

diff --git a/PWGJE/PWGJEEMCALJetTasksLinkDef.h b/PWGJE/PWGJEEMCALJetTasksLinkDef.h
index e9a2cbacbc6522bd0b681e3694996ad771846773..9e44464c7284bdb19b74144847423451fde7f5cf 100644 (file)
--- a/PWGJE/PWGJEEMCALJetTasksLinkDef.h
+++ b/PWGJE/PWGJEEMCALJetTasksLinkDef.h
@@ -52,6 +52,8 @@
 #pragma link C++ class AliAnalysisTaskEmcalJetTriggerQA+;
 #pragma link C++ class AliAnalysisTaskEmcalTriggerInfoQA+;
 #pragma link C++ class AliAnalysisTaskHJetEmbed+;
+#pragma link C++ class AliAnalysisTaskJetShapeDeriv+;
+#pragma link C++ class AliAnalysisTaskJetShapeConst+;
 #pragma link C++ class AliAnalysisTaskJetMatching+;
 #pragma link C++ class AliAnalysisTaskJetV2+;
 #pragma link C++ class AliAnalysisTaskRhoMass+;