Reorganization of the EMCal jet finder task: the manin task only performs the basic...
authorsaiola <salvatore.aiola@cern.ch>
Fri, 12 Dec 2014 21:03:25 +0000 (16:03 -0500)
committersaiola <salvatore.aiola@cern.ch>
Fri, 12 Dec 2014 21:03:52 +0000 (16:03 -0500)
13 files changed:
PWG/EMCAL/AliEmcalParticle.cxx
PWG/EMCAL/AliEmcalParticle.h
PWGJE/CMakelibPWGJEEMCALJetTasks.pkg
PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx
PWGJE/EMCALJetTasks/AliEmcalJetTask.h
PWGJE/EMCALJetTasks/AliEmcalJetUtility.cxx [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliEmcalJetUtility.h [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.cxx [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.h [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.cxx [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.h [new file with mode: 0644]
PWGJE/EMCALJetTasks/AliFJWrapper.h
PWGJE/PWGJEEMCALJetTasksLinkDef.h

index f6c70a9..80a7939 100644 (file)
@@ -1,4 +1,3 @@
-// $Id$
 //
 // Emcal particle class, which can contain either an AliVTrack or an AliVCluster
 //
index 3b2803f..689bed8 100644 (file)
@@ -1,8 +1,6 @@
 #ifndef ALIEMCALPARTICLE_H
 #define ALIEMCALPARTICLE_H
 
-// $Id$
-
 #include <TLorentzVector.h>
 #include <TMath.h>
 #include <TObjArray.h>
@@ -69,6 +67,10 @@ class AliEmcalParticle: public AliVParticle {
   void              SetMatchedObj(Int_t id, Double_t d)  { ResetMatchedObjects(); fMatchedIds[0] = id; fMatchedDist[0] = d; fNMatched = 1; }
   void              SetMatchedPtr(TObjArray *arr)        { fMatchedPtr = arr; }
 
+  void              SetPt(Double_t pt)           { fPt  = pt ; }
+  void              SetPhi(Double_t phi)         { fPhi = phi; }
+  void              SetEta(Double_t eta)         { fEta = eta; }
+
  protected:
   TLorentzVector   &GetLorentzVector(const Double_t *vertex = 0)  const;
 
index 8c85768..a891508 100644 (file)
@@ -126,6 +126,9 @@ set ( SRCS
 if (FASTJET_FOUND)
   LIST(APPEND SRCS
    EMCALJetTasks/AliFJWrapper.cxx
+   EMCALJetTasks/AliEmcalJetUtility.cxx
+   EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.cxx
+   EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.cxx
    EMCALJetTasks/AliEmcalJetTask.cxx
    EMCALJetTasks/AliEmcalJetFinder.cxx
    EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx
index 856d5ff..620ab9d 100644 (file)
@@ -12,7 +12,6 @@
 #include <TLorentzVector.h>
 #include <TMath.h>
 #include <TRandom3.h>
-#include <TParticle.h>
 
 #include "AliAnalysisManager.h"
 #include "AliCentrality.h"
 #include "AliFJWrapper.h"
 #include "AliMCEvent.h"
 #include "AliVCluster.h"
-#include "AliAODCaloCluster.h"
-#include "AliESDCaloCluster.h"
 #include "AliVEvent.h"
 #include "AliVParticle.h"
-#include "AliRhoParameter.h"
-#include "AliPicoTrack.h"
+#include "AliEmcalJetUtility.h"
 
 using std::cout;
 using std::endl;
@@ -42,13 +38,12 @@ AliEmcalJetTask::AliEmcalJetTask() :
   fTracksName("Tracks"),
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
-  fJetsSubName(""),
-  fTracksSubName(""),
-  fClusSubName(""),
-  fJetType(kNone),
-  fConstSel(0),
-  fMCConstSel(0),
-  fMarkConst(kFALSE),
+  fJetType(kAKT|kFullJet|kRX1Jet),
+  fMinLabelTracks(-kMaxInt),
+  fMaxLabelTracks(kMaxInt),
+  fMinLabelClusters(-kMaxInt),
+  fMaxLabelClusters(kMaxInt),
+  fMinMCLabel(0),
   fRadius(0.4),
   fMinJetTrackPt(0.15),
   fMinJetClusPt(0.15),
@@ -56,51 +51,35 @@ AliEmcalJetTask::AliEmcalJetTask() :
   fPhiMax(TMath::TwoPi()),
   fEtaMin(-0.9),
   fEtaMax(+0.9),
-  fMinJetArea(0.005),
+  fMinJetArea(0.001),
   fMinJetPt(1.0),
   fJetPhiMin(-10),
   fJetPhiMax(+10),
   fJetEtaMin(-1),
   fJetEtaMax(+1),
   fGhostArea(0.005),
-  fMinMCLabel(0),
   fRecombScheme(fastjet::pt_scheme),
   fTrackEfficiency(1.),
   fMCFlag(0),
   fGeneratorIndex(-1),
-  fIsInit(0),
+  fUtilities(0),
   fLocked(0),
+  fIsInit(0),
   fIsPSelSet(0),
-  fIsMcPart(0),
   fIsEmcPart(0),
-  fIsPicoTrack(0),
-  fIsEsdClus(0),
   fLegacyMode(kFALSE),
-  fCodeDebug(kFALSE),
-  fPionMassClusters(kFALSE),
-  fDoGenericSubtractionJetMass(kFALSE),
-  fDoGenericSubtractionGR(kFALSE),
-  fDoGenericSubtractionExtraJetShapes(kFALSE),
-  fDoConstituentSubtraction(kFALSE),
-  fUseExternalBkg(kFALSE),
-  fRhoName(""),
-  fRhomName(""),
-  fRho(0),
-  fRhom(0),
-  fRMax(0.4),
-  fDRStep(0.04),
-  fPtMinGR(40.),
+  fGeom(0),
   fJets(0),
-  fJetsSub(0),
-  fTracksSub(0),
-  fClusSub(0),
   fEvent(0),
   fTracks(0),
   fClus(0),
-  fRhoParam(0),
-  fRhomParam(0)
+  fFastJetWrapper("AliEmcalJetTask","AliEmcalJetTask")
 {
   // Default constructor.
+
+  fVertex[0] = 0.;
+  fVertex[1] = 0.;
+  fVertex[2] = 0.;
 }
 
 //________________________________________________________________________
@@ -109,13 +88,12 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) :
   fTracksName("Tracks"),
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
-  fJetsSubName(""),
-  fTracksSubName(""),
-  fClusSubName(""),
   fJetType(kAKT|kFullJet|kRX1Jet),
-  fConstSel(0),
-  fMCConstSel(0),
-  fMarkConst(kFALSE),
+  fMinLabelTracks(-kMaxInt),
+  fMaxLabelTracks(kMaxInt),
+  fMinLabelClusters(-kMaxInt),
+  fMaxLabelClusters(kMaxInt),
+  fMinMCLabel(0),
   fRadius(0.4),
   fMinJetTrackPt(0.15),
   fMinJetClusPt(0.15),
@@ -130,46 +108,30 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) :
   fJetEtaMin(-1),
   fJetEtaMax(+1),
   fGhostArea(0.005),
-  fMinMCLabel(0),
   fRecombScheme(fastjet::pt_scheme),
   fTrackEfficiency(1.),
   fMCFlag(0),
   fGeneratorIndex(-1),
-  fIsInit(0),
+  fUtilities(0),
   fLocked(0),
+  fIsInit(0),
   fIsPSelSet(0),
-  fIsMcPart(0),
   fIsEmcPart(0),
-  fIsPicoTrack(0),
-  fIsEsdClus(0),
   fLegacyMode(kFALSE),
-  fCodeDebug(kFALSE),
-  fPionMassClusters(kFALSE),
-  fDoGenericSubtractionJetMass(kFALSE),
-  fDoGenericSubtractionGR(kFALSE),
-  fDoGenericSubtractionExtraJetShapes(kFALSE),
-  fDoConstituentSubtraction(kFALSE),
-  fUseExternalBkg(kFALSE),
-  fRhoName(""),
-  fRhomName(""),
-  fRho(0),
-  fRhom(0),
-  fRMax(0.4),
-  fDRStep(0.04),
-  fPtMinGR(40.),
+  fGeom(0),
   fJets(0),
-  fJetsSub(0),
-  fTracksSub(0),
-  fClusSub(0),
   fEvent(0),
   fTracks(0),
   fClus(0),
-  fRhoParam(0),
-  fRhomParam(0)
+  fFastJetWrapper(name,name)
 {
   // Standard constructor.
 
   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
+
+  fVertex[0] = 0.;
+  fVertex[1] = 0.;
+  fVertex[2] = 0.;
 }
 
 //________________________________________________________________________
@@ -179,28 +141,57 @@ AliEmcalJetTask::~AliEmcalJetTask()
 }
 
 //________________________________________________________________________
-void AliEmcalJetTask::UserCreateOutputObjects()
+AliEmcalJetUtility* AliEmcalJetTask::AddUtility(AliEmcalJetUtility* utility)
 {
-  // Create user objects.
+  // Addition of utilities.
 
-  fJets = new TClonesArray("AliEmcalJet");
-  fJets->SetName(fJetsName);
+  if (!fUtilities) fUtilities = new TObjArray();
+  if (fUtilities->FindObject(utility)) {
+     Error("AddSupply", "Jet utility %s already connected.", utility->GetName());
+     return utility;
+  }   
+  fUtilities->Add(utility);
+  utility->SetJetTask(this);
 
-  if(!fJetsSubName.IsNull()) {
-    fJetsSub = new TClonesArray("AliEmcalJet");
-    fJetsSub->SetName(fJetsSubName);
-  }
+  return utility;
+}
 
-  if (fDoConstituentSubtraction && !fTracksSubName.IsNull()) {
-    fTracksSub = new TClonesArray();
-    fTracksSub->SetName(fTracksSubName);
-  }
+//________________________________________________________________________
+void AliEmcalJetTask::InitUtilities()
+{
+  TIter next(fUtilities);
+  AliEmcalJetUtility *utility = 0;
+  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Init();
+}
 
-  if (fDoConstituentSubtraction && !fClusSubName.IsNull()) {
-    fClusSub = new TClonesArray();
-    fClusSub->SetName(fClusSubName);
-  }
-  
+//________________________________________________________________________
+void AliEmcalJetTask::PrepareUtilities()
+{
+  TIter next(fUtilities);
+  AliEmcalJetUtility *utility = 0;
+  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Prepare(fFastJetWrapper);
+}
+
+//________________________________________________________________________
+void AliEmcalJetTask::ExecuteUtilities(AliEmcalJet* jet, Int_t ij)
+{
+  TIter next(fUtilities);
+  AliEmcalJetUtility *utility = 0;
+  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->ProcessJet(jet, ij, fFastJetWrapper);
+}
+
+//________________________________________________________________________
+void AliEmcalJetTask::TerminateUtilities()
+{
+  TIter next(fUtilities);
+  AliEmcalJetUtility *utility = 0;
+  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Terminate(fFastJetWrapper);
+}
+
+//________________________________________________________________________
+void AliEmcalJetTask::UserCreateOutputObjects()
+{
+  // Create user objects.
 }
 
 //________________________________________________________________________
@@ -215,117 +206,68 @@ void AliEmcalJetTask::UserExec(Option_t *)
 
   // clear the jet array (normally a null operation)
   fJets->Delete();
-  if(fJetsSub)   fJetsSub->Delete();
-  if(fTracksSub) fTracksSub->Delete();
-  if(fClusSub)   fClusSub->Delete();
 
-  //Run jet finding
+  // get primary vertex
+  if(fEvent->GetPrimaryVertex()) fEvent->GetPrimaryVertex()->GetXYZ(fVertex);
+
   FindJets();
-}
 
-//________________________________________________________________________
-void AliEmcalJetTask::Terminate(Option_t *)
-{
-  // Called once at the end of the analysis.
+  FillJetBranch();
 }
 
 //________________________________________________________________________
 void AliEmcalJetTask::FindJets()
 {
   // Find jets.
+
   if (!fTracks && !fClus){
-    AliInfo("WARNING NO TRACKS OR CLUSTERS");
+    AliError("No tracks or clusters, returning.");
     return;
   }
 
- if (fRhoParam)
-    fRho = fRhoParam->GetVal();
-  if (fRhomParam)
-    fRhom = fRhomParam->GetVal();
-
-  TString name("kt");
-  fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
-  if ((fJetType & kAKT) != 0) {
-    name  = "antikt";
-    jalgo = fastjet::antikt_algorithm;
-    AliDebug(1,"Using AKT algorithm");
-  }
-  else {
-    AliDebug(1,"Using KT algorithm");
-  }
-
-  if ((fJetType & kR020Jet) != 0)
-    fRadius = 0.2;
-  else if ((fJetType & kR030Jet) != 0)
-    fRadius = 0.3;
-  else if ((fJetType & kR040Jet) != 0)
-    fRadius = 0.4;
-
-  // setup fj wrapper
-  AliFJWrapper fjw(name, name);
-  fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
-  fjw.SetGhostArea(fGhostArea);
-  fjw.SetR(fRadius);
-  fjw.SetAlgorithm(jalgo);
-  fjw.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme));
-  fjw.SetMaxRap(fEtaMax);
-  fjw.Clear();
-
-  // get primary vertex
-  Double_t vertex[3] = {0, 0, 0};
-  if(fEvent->GetPrimaryVertex()) fEvent->GetPrimaryVertex()->GetXYZ(vertex);
+  fFastJetWrapper.Clear();
 
   AliDebug(2,Form("Jet type = %d", fJetType));
 
-  if ((fIsMcPart || ((fJetType & kFullJet) != 0) || ((fJetType & kChargedJet) != 0)) && fTracks) {
+  if (fTracks) {
     const Int_t Ntracks = fTracks->GetEntries();
     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
       AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
-      if (!t)
+
+      if (!t) continue;
+
+      if (t->Pt() < fMinJetTrackPt) continue;
+      Double_t eta = t->Eta();
+      Double_t phi = t->Phi();
+      if ((eta<fEtaMin) || (eta>fEtaMax) ||
+          (phi<fPhiMin) || (phi>fPhiMax))
+        continue;
+
+      if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) {
+        AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks));
         continue;
-      if (fIsMcPart) {
-       if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) {
-         AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks));
-         continue;
-       }
-       if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) {
-         AliDebug(2,Form("Skipping track %d because it is charged.", iTracks));
-         continue;
-       }
       }
-      if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) {
-       if (t->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
-         AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
-         continue;
-       }
-       else {
-         AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
-       }
+
+      if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) {
+        AliDebug(2,Form("Skipping track %d because it is charged.", iTracks));
+        continue;
       }
-      else {
-       if (t->TestBits(fConstSel) != (Int_t)fConstSel) {
-         AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
-         continue;
-       }
-       else {
-         AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
-       }
+      
+      Int_t lab = TMath::Abs(t->GetLabel());
+      if (lab < fMinLabelTracks || lab > fMaxLabelTracks) {
+        AliDebug(2,Form("Skipping track %d because label %d is not in range (%d, %d)", iTracks, lab, fMinLabelTracks, fMaxLabelTracks));
+        continue;
       }
+
       if ((t->GetFlag() & fMCFlag) != fMCFlag) {
        AliDebug(2,Form("Skipping track %d because it does not match the MC flags", iTracks));
        continue;
       }
+
       if (fGeneratorIndex >= 0 && t->GetGeneratorIndex() != fGeneratorIndex) {
        AliDebug(2,Form("Skipping track %d because it does not match the MC generator index", iTracks));
        continue;
       }
-      if (t->Pt() < fMinJetTrackPt)
-        continue;
-      Double_t eta = t->Eta();
-      Double_t phi = t->Phi();
-      if ((eta<fEtaMin) || (eta>fEtaMax) ||
-          (phi<fPhiMin) || (phi>fPhiMax))
-        continue;
 
       // artificial inefficiency
       if (fTrackEfficiency < 1.) {
@@ -337,12 +279,12 @@ void AliEmcalJetTask::FindJets()
       }
 
       // offset of 100 for consistency with cluster ids
-      AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, t->GetLabel(), t->Pt()));
-      fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100);
+      AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, lab, t->Pt()));
+      fFastJetWrapper.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100);
     }
   }
 
-  if ((((fJetType & kFullJet) != 0) || ((fJetType & kNeutralJet) != 0)) && fClus) {
+  if (fClus) {
     const Int_t Nclus = fClus->GetEntries();
     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
       AliVCluster *c = 0;
@@ -350,31 +292,9 @@ void AliEmcalJetTask::FindJets()
       Double_t cPx=0,cPy=0,cPz=0;
       if (fIsEmcPart) {
        AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
-       if (!ep)
-         continue;
-
+       if (!ep) continue;
        c = ep->GetCluster();
-       if (!c)
-         continue;
-
-       if (c->GetLabel() > fMinMCLabel) {
-         if (ep->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
-           AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
-           continue;
-         }
-         else {
-           AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
-         }
-       }
-       else {
-         if (ep->TestBits(fConstSel) != (Int_t)fConstSel) {
-           AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
-           continue;
-         }
-         else {
-           AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
-         }
-       }
+       if (!c) continue;
 
        cEta = ep->Eta();
        cPhi = ep->Phi();
@@ -382,32 +302,13 @@ void AliEmcalJetTask::FindJets()
        cPx  = ep->Px();
        cPy  = ep->Py();
        cPz  = ep->Pz();
-      } else {
+      } 
+      else {
        c = static_cast<AliVCluster*>(fClus->At(iClus));
-       if (!c)
-         continue;
-
-       if (c->GetLabel() > fMinMCLabel) {
-         if (c->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
-           AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
-           continue;
-         }
-         else {
-           AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
-         }
-       }
-       else {
-         if (c->TestBits(fConstSel) != (Int_t)fConstSel) {
-           AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
-           continue;
-         }
-         else {
-           AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
-         }
-       }
+       if (!c) continue;
 
        TLorentzVector nP;
-       c->GetMomentum(nP, vertex);
+       c->GetMomentum(nP, fVertex);
        cEta = nP.Eta();
        cPhi = nP.Phi();
        cPt  = nP.Pt();
@@ -415,278 +316,83 @@ void AliEmcalJetTask::FindJets()
        cPy  = nP.Py();
        cPz  = nP.Pz();
       }
-      if (!c->IsEMCAL())
-       continue;
-      if (cPt < fMinJetClusPt)
-       continue;
+
+      if (!c->IsEMCAL()) continue;
+      if (cPt < fMinJetClusPt) continue;
       if ((cEta<fEtaMin) || (cEta>fEtaMax) ||
          (cPhi<fPhiMin) || (cPhi>fPhiMax))
        continue;
+
+      Int_t lab = TMath::Abs(c->GetLabel());
+      if (lab < fMinLabelClusters || lab > fMaxLabelClusters) {
+        AliDebug(2,Form("Skipping cluster %d because label %d is not in range (%d, %d)", iClus, lab, fMinLabelClusters, fMaxLabelClusters));
+        continue;
+      }
+
       // offset of 100 to skip ghost particles uid = -1
       AliDebug(2,Form("Cluster %d accepted (label = %d)", iClus, c->GetLabel()));
       Double_t e = TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz);
-      if(fPionMassClusters) e = TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz + 0.13957*0.13957); //MV: dirty, need better solution
-      fjw.AddInputVector(cPx, cPy, cPz, e, -iClus - 100);
-      //      fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);
+      fFastJetWrapper.AddInputVector(cPx, cPy, cPz, e, -iClus - 100);
     }
   }
 
-  // setting legacy mode
-  if (fLegacyMode) {
-    fjw.SetLegacyMode(kTRUE);
-  }
-
   // run jet finder
-  fjw.Run();
-
-  //run generic subtractor
-  if(fDoGenericSubtractionJetMass) {
-    fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
-    fjw.DoGenericSubtractionJetMass();
-  }
-
-  if(fDoGenericSubtractionExtraJetShapes) {
-    fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
-    fjw.DoGenericSubtractionJetAngularity();
-    fjw.DoGenericSubtractionJetpTD();
-    fjw.DoGenericSubtractionJetCircularity();
-    fjw.DoGenericSubtractionJetSigma2();
-    fjw.DoGenericSubtractionJetConstituent();
-    fjw.DoGenericSubtractionJetLeSub();
-  }
+  fFastJetWrapper.Run();
+}
 
-  //run constituent subtractor
-  if(fDoConstituentSubtraction) {
-    fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
-    fjw.DoConstituentSubtraction();
-  }
+//________________________________________________________________________
+void AliEmcalJetTask::FillJetBranch()
+{
+  // Fill jet branch with fastjet jets.
 
-  // get geometry
-  AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
-  if (!geom) {
-    AliFatal(Form("%s: Can not create geometry", GetName()));
-    return;
-  }
+  PrepareUtilities();
 
   // loop over fastjet jets
-  std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
+  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper.GetInclusiveJets();
   // sort jets according to jet pt
   static Int_t indexes[9999] = {-1};
   GetSortedArray(indexes, jets_incl);
 
   AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
-  for (UInt_t ijet=0, jetCount=0; ijet<jets_incl.size(); ++ijet) {
+  for (UInt_t ijet = 0, jetCount = 0; ijet < jets_incl.size(); ++ijet) {
     Int_t ij = indexes[ijet];
-    AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fjw.GetJetArea(ij)));
+    AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fFastJetWrapper.GetJetArea(ij)));
 
-    if (jets_incl[ij].perp()<fMinJetPt)
-      continue;
-    if (fjw.GetJetArea(ij)<fMinJetArea)
-      continue;
-    if ((jets_incl[ij].eta()<fJetEtaMin) || (jets_incl[ij].eta()>fJetEtaMax) ||
-       (jets_incl[ij].phi()<fJetPhiMin) || (jets_incl[ij].phi()>fJetPhiMax))
+    if (jets_incl[ij].perp() < fMinJetPt) continue;
+    if (fFastJetWrapper.GetJetArea(ij) < fMinJetArea) continue;
+    if ((jets_incl[ij].eta() < fJetEtaMin) || (jets_incl[ij].eta() > fJetEtaMax) ||
+       (jets_incl[ij].phi() < fJetPhiMin) || (jets_incl[ij].phi() > fJetPhiMax))
       continue;
 
     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
-#ifdef FASTJET_VERSION
-    if(fDoGenericSubtractionJetMass) {
-      std::vector<fastjet::contrib::GenericSubtractorInfo> jetMassInfo = fjw.GetGenSubtractorInfoJetMass();
-      Int_t n = (Int_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());
-      }
-    }
-
-    //here do generic subtraction for angular structure function
-    Double_t ptcorr = jets_incl[ij].perp()-fjw.GetJetArea(ij)*fRho;
-    fRMax = fRadius+0.2;
-    if(fDoGenericSubtractionGR && ptcorr>fPtMinGR) {
-      fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
-      fjw.SetRMaxAndStep(fRMax,fDRStep);
-      fjw.DoGenericSubtractionGR(ij);
-      std::vector<double> num = fjw.GetGRNumerator();
-      std::vector<double> den = fjw.GetGRDenominator();
-      std::vector<double> nums = fjw.GetGRNumeratorSub();
-      std::vector<double> dens = fjw.GetGRDenominatorSub();
-      //pass this to AliEmcalJet
-      jet->SetGRNumSize(num.size());
-      jet->SetGRDenSize(den.size());
-      jet->SetGRNumSubSize(nums.size());
-      jet->SetGRDenSubSize(dens.size());
-      Int_t nsize = (Int_t)num.size();
-      for(Int_t g = 0; g<nsize; ++g) {
-        jet->AddGRNumAt(num[g],g);
-        jet->AddGRNumSubAt(nums[g],g);
-      }
-      Int_t dsize = (Int_t)den.size();
-      for(Int_t g = 0; g<dsize; ++g) {
-        jet->AddGRDenAt(den[g],g);
-        jet->AddGRDenSubAt(dens[g],g);
-      }
-    }
-
-   if(fDoGenericSubtractionExtraJetShapes){
-      std::vector<fastjet::contrib::GenericSubtractorInfo> jetAngularityInfo = fjw.GetGenSubtractorInfoJetAngularity();
-      Int_t na = (Int_t)jetAngularityInfo.size();
-      if(na>ij && na>0) {
-       jet->SetFirstDerivativeAngularity(jetAngularityInfo[ij].first_derivative());
-       jet->SetSecondDerivativeAngularity(jetAngularityInfo[ij].second_derivative());
-       jet->SetFirstOrderSubtractedAngularity(jetAngularityInfo[ij].first_order_subtracted());
-       jet->SetSecondOrderSubtractedAngularity(jetAngularityInfo[ij].second_order_subtracted());
-      }
-
-      std::vector<fastjet::contrib::GenericSubtractorInfo> jetpTDInfo = fjw.GetGenSubtractorInfoJetpTD();
-      Int_t np = (Int_t)jetpTDInfo.size();
-      if(np>ij && np>0) {
-       jet->SetFirstDerivativepTD(jetpTDInfo[ij].first_derivative());
-       jet->SetSecondDerivativepTD(jetpTDInfo[ij].second_derivative());
-       jet->SetFirstOrderSubtractedpTD(jetpTDInfo[ij].first_order_subtracted());
-       jet->SetSecondOrderSubtractedpTD(jetpTDInfo[ij].second_order_subtracted());
-      }
-
-      std::vector<fastjet::contrib::GenericSubtractorInfo> jetCircularityInfo = fjw.GetGenSubtractorInfoJetCircularity();
-      Int_t nc = (Int_t)jetCircularityInfo.size();
-      if(nc>ij && nc>0) {
-       jet->SetFirstDerivativeCircularity(jetCircularityInfo[ij].first_derivative());
-       jet->SetSecondDerivativeCircularity(jetCircularityInfo[ij].second_derivative());
-       jet->SetFirstOrderSubtractedCircularity(jetCircularityInfo[ij].first_order_subtracted());
-       jet->SetSecondOrderSubtractedCircularity(jetCircularityInfo[ij].second_order_subtracted());
-      }
-
-      std::vector<fastjet::contrib::GenericSubtractorInfo> jetSigma2Info = fjw.GetGenSubtractorInfoJetSigma2();
-      Int_t ns = (Int_t)jetSigma2Info.size();
-      if(ns>ij && ns>0) {
-       jet->SetFirstDerivativeSigma2(jetSigma2Info[ij].first_derivative());
-       jet->SetSecondDerivativeSigma2(jetSigma2Info[ij].second_derivative());
-       jet->SetFirstOrderSubtractedSigma2(jetSigma2Info[ij].first_order_subtracted());
-       jet->SetSecondOrderSubtractedSigma2(jetSigma2Info[ij].second_order_subtracted());
-      }
-
-
-      std::vector<fastjet::contrib::GenericSubtractorInfo> jetConstituentInfo = fjw.GetGenSubtractorInfoJetConstituent();
-      Int_t nco= (Int_t)jetConstituentInfo.size();
-      if(nco>ij && nco>0) {
-       jet->SetFirstDerivativeConstituent(jetConstituentInfo[ij].first_derivative());
-       jet->SetSecondDerivativeConstituent(jetConstituentInfo[ij].second_derivative());
-       jet->SetFirstOrderSubtractedConstituent(jetConstituentInfo[ij].first_order_subtracted());
-       jet->SetSecondOrderSubtractedConstituent(jetConstituentInfo[ij].second_order_subtracted());
-      }
-
-      std::vector<fastjet::contrib::GenericSubtractorInfo> jetLeSubInfo = fjw.GetGenSubtractorInfoJetLeSub();
-      Int_t nlsub= (Int_t)jetLeSubInfo.size();
-      if(nlsub>ij && nlsub>0) {
-       jet->SetFirstDerivativeLeSub(jetLeSubInfo[ij].first_derivative());
-       jet->SetSecondDerivativeLeSub(jetLeSubInfo[ij].second_derivative());
-       jet->SetFirstOrderSubtractedLeSub(jetLeSubInfo[ij].first_order_subtracted());
-       jet->SetSecondOrderSubtractedLeSub(jetLeSubInfo[ij].second_order_subtracted());
-      }
-   }
-#endif
-
-    // Fill constituent info
-    std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
-    jet->SetNumberOfTracks(constituents.size());
-    jet->SetNumberOfClusters(constituents.size());
-
-    Int_t nt            = 0;
-    Int_t nc            = 0;
-    Double_t neutralE   = 0;
-    Double_t maxCh      = 0;
-    Double_t maxNe      = 0;
-    Int_t gall          = 0;
-    Int_t gemc          = 0;
-    Int_t cemc          = 0;
-    Int_t ncharged      = 0;
-    Int_t nneutral      = 0;
-    Double_t mcpt       = 0;
-    Double_t emcpt      = 0;
-
-    FillJetConstituents(constituents,jet,vertex,jetCount,nt,nc,maxCh,maxNe,ncharged,nneutral,neutralE,mcpt,cemc,emcpt,gall,gemc,constituents,0);
-    jet->SetNumberOfTracks(nt);
-    jet->SetNumberOfClusters(nc);
-    jet->SortConstituents();
-    jet->SetMaxNeutralPt(maxNe);
-    jet->SetMaxChargedPt(maxCh);
-    jet->SetNEF(neutralE / jet->E());
-    fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
+    fastjet::PseudoJet area(fFastJetWrapper.GetJetAreaVector(ij));
     jet->SetArea(area.perp());
     jet->SetAreaEta(area.eta());
     jet->SetAreaPhi(area.phi());
     jet->SetAreaE(area.E());
-    jet->SetNumberOfCharged(ncharged);
-    jet->SetNumberOfNeutrals(nneutral);
-    jet->SetMCPt(mcpt);
-    jet->SetNEmc(cemc);
-    jet->SetPtEmc(emcpt);
-
-    if (gall > 0)
-      jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
-    else
-      jet->SetAreaEmc(-1);
-    if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
-       (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
-       (jet->Eta() > geom->GetArm1EtaMin()) &&
-       (jet->Eta() < geom->GetArm1EtaMax()))
-      jet->SetAxisInEmcal(kTRUE);
-
-    AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size()));
+
+    // Fill constituent info
+    std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper.GetJetConstituents(ij));
+    FillJetConstituents(jet, constituents, fTracks, fClus, constituents);
+
+    if (fGeom) {
+      if ((jet->Phi() > fGeom->GetArm1PhiMin() * TMath::DegToRad()) &&
+          (jet->Phi() < fGeom->GetArm1PhiMax() * TMath::DegToRad()) &&
+          (jet->Eta() > fGeom->GetArm1EtaMin()) &&
+          (jet->Eta() < fGeom->GetArm1EtaMax()))
+        jet->SetAxisInEmcal(kTRUE);
+    }
+
+    ExecuteUtilities(jet, ij);
+
+    AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), jet->GetNumberOfConstituents()));
     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);
-        // Fill constituent info
-        std::vector<fastjet::PseudoJet> constituents_unsub(fjw.GetJetConstituents(ijet));
-       std::vector<fastjet::PseudoJet> constituents_sub = jets_sub[ijet].constituents();
-        jet_sub->SetNumberOfTracks(constituents_sub.size());
-       jet_sub->SetNumberOfClusters(constituents_sub.size());
-       Int_t nt            = 0;
-       Int_t nc            = 0;
-       Double_t neutralE   = 0;
-       Double_t maxCh      = 0;
-       Double_t maxNe      = 0;
-       Int_t gall          = 0;
-       Int_t gemc          = 0;
-       Int_t cemc          = 0;
-       Int_t ncharged      = 0;
-       Int_t nneutral      = 0;
-       Double_t mcpt       = 0;
-       Double_t emcpt      = 0;
-
-       FillJetConstituents(constituents_sub,jet_sub,vertex,jetCount,nt,nc,maxCh,maxNe,ncharged,nneutral,neutralE,mcpt,cemc,emcpt,gall,gemc,constituents_unsub,1);
-       jet_sub->SetNumberOfTracks(nt);
-       jet_sub->SetNumberOfClusters(nc);
-       jet_sub->SortConstituents();
-
-       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
+  TerminateUtilities();
 }
 
 //________________________________________________________________________
@@ -719,6 +425,12 @@ Bool_t AliEmcalJetTask::DoInit()
     gRandom = new TRandom3(0);
   }
 
+  // get geometry
+  fGeom = AliEMCALGeometry::GetInstance();
+  if (!fGeom) {
+    AliWarning(Form("%s: Can not create EMCal geometry, some features will not work...", GetName()));
+  }
+
   // get input collections
   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 
@@ -738,15 +450,6 @@ Bool_t AliEmcalJetTask::DoInit()
       return 0;
     }
   }
-  if (fTracks) {
-    TClass cls(fTracks->GetClass()->GetName());
-    if (cls.InheritsFrom("AliPicoTrack"))
-      fIsPicoTrack = 1; 
-    else if (cls.InheritsFrom("AliMCParticle") || cls.InheritsFrom("AliAODMCParticle"))
-     fIsMcPart = 1;
-    else if (cls.InheritsFrom("AliEmcalParticle"))
-     fIsEmcPart = 1;
-  }
 
   if (fCaloName == "CaloClusters")
     am->LoadBranch("CaloClusters");
@@ -763,249 +466,225 @@ Bool_t AliEmcalJetTask::DoInit()
       fIsEmcPart = 1;
   }
 
-  //Add constituent subtracted jets to event
-  if(!fJetsSubName.IsNull()) {
-    if (!(fEvent->FindListObject(fJetsSubName)) ) {
-      if(fJetsSub) fEvent->AddObject(fJetsSub);
-    } else {
-      AliError(Form("%s: Object for subtracted jet branch with name %s already in event! Returning", GetName(), fJetsSubName.Data()));
-      return 0;
-    }
+  // add jets to event if not yet there
+  if (!(fEvent->FindListObject(fJetsName))) {
+    fJets = new TClonesArray("AliEmcalJet");
+    fJets->SetName(fJetsName);
+    fEvent->AddObject(fJets);
   }
-
-  //Add tracks from constituent subtracted jets to event
-  if(!fTracksSubName.IsNull()) {
-    if (!(fEvent->FindListObject(fTracksSubName))) {
-      if(fTracksSub) fEvent->AddObject(fTracksSub);
-    } else {
-      AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fTracksSubName.Data()));
-      return 0;
-    }
+  else {
+    AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
+    return 0;
   }
 
-  if(fDoConstituentSubtraction && fTracksSub) {
-    if(fIsEmcPart==1) { 
-      AliFatal("Constituent subtraction method not implemented for AliEmcalParticle as input tracks. Use AliPicoTracks instead. Aborting!");
-      return 0;
-    }
-    if(fIsPicoTrack==1)   fTracksSub->SetClass("AliPicoTrack");
-    else if(fIsMcPart==1) fTracksSub->SetClass("AliMCParticle");
-    else {
-      AliError(Form("%s Object type of subtracted track branch %s not known",GetName(),fTracksSubName.Data()));
-      return 0;
-    }
-    fTracksSub->SetName(fTracksSubName);
+  TString name("kt");
+  fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
+  if ((fJetType & kAKT) != 0) {
+    name  = "antikt";
+    jalgo = fastjet::antikt_algorithm;
+    AliDebug(1,"Using AKT algorithm");
   }
-
-  //Add clusters from constituent subtracted jets to event
-  if(!fClusSubName.IsNull()) {
-    if (!(fEvent->FindListObject(fClusSubName))) {
-      if(fClusSub) fEvent->AddObject(fClusSub);
-    } else {
-      AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fClusSubName.Data()));
-      return 0;
-    }
+  else {
+    AliDebug(1,"Using KT algorithm");
   }
 
-  if(fDoConstituentSubtraction && fClusSub) {
-    TClass cls(fClus->GetClass()->GetName());
-    if (cls.InheritsFrom("AliESDCaloCluster")) {
-      fClusSub->SetClass("AliESDCaloCluster");
-      fIsEsdClus = kTRUE;
-    }
-    else if (cls.InheritsFrom("AliAODCaloCluster"))
-      fClusSub->SetClass("AliAODCaloCluster");
-    else if (cls.InheritsFrom("AliEmcalParticle")) { 
-      AliFatal("Constituent subtraction method not implemented for AliEmcalParticle as input clusters. Use AliESD(AOD)CaloCluster instead. Aborting!");
-      return 0;
-    } else {
-      AliError(Form("%s Object type of subtracted cluster branch %s not known",GetName(),fClusSubName.Data()));
-      return 0;
-    }
-    fClusSub->SetName(fClusSubName);
-  }
+  if ((fJetType & kR020Jet) != 0)
+    fRadius = 0.2;
+  else if ((fJetType & kR030Jet) != 0)
+    fRadius = 0.3;
+  else if ((fJetType & kR040Jet) != 0)
+    fRadius = 0.4;
 
-  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_m %s!", GetName(), fRhomName.Data()));
-      return 0;
-    }
-  }
+  // setup fj wrapper
+  fFastJetWrapper.SetName(name);
+  fFastJetWrapper.SetTitle(name);
+  fFastJetWrapper.SetAreaType(fastjet::active_area_explicit_ghosts);
+  fFastJetWrapper.SetGhostArea(fGhostArea);
+  fFastJetWrapper.SetR(fRadius);
+  fFastJetWrapper.SetAlgorithm(jalgo);
+  fFastJetWrapper.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme));
+  fFastJetWrapper.SetMaxRap(fEtaMax);
 
-  // add jets to event if not yet there
-  if (!(fEvent->FindListObject(fJetsName)))
-    fEvent->AddObject(fJets);
-  else {
-    AliError(Form("%s: Object for jet branch with name %s already in event! Returning", GetName(), fJetsName.Data()));
-    return 0;
+  // setting legacy mode
+  if (fLegacyMode) {
+    fFastJetWrapper.SetLegacyMode(kTRUE);
   }
 
-  return 1;
+  InitUtilities();
+
+  return kTRUE;
 }
 
 //___________________________________________________________________________________________________________________
-void  AliEmcalJetTask::FillJetConstituents(std::vector<fastjet::PseudoJet>& constituents,AliEmcalJet *jet,Double_t vertex[3],UInt_t jetCount,Int_t& nt,Int_t& nc,Double_t& maxCh,Double_t& maxNe,Int_t& ncharged,Int_t& nneutral,Double_t& neutralE,Double_t& mcpt,Int_t& cemc,Double_t& emcpt,Int_t& gall,Int_t& gemc,std::vector<fastjet::PseudoJet>& constituentsunsub,Int_t flagsub){
-    nt         = 0;
-    nc         = 0;
-    neutralE   = 0;
-    maxCh      = 0;
-    maxNe      = 0;
-    gall       = 0;
-    gemc       = 0;
-    cemc       = 0;
-    ncharged   = 0;
-    nneutral   = 0;
-    mcpt       = 0;
-    emcpt      = 0;
-    Int_t uid   = -1;
-    AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
-    for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
-      if(flagsub==0) uid = constituents[ic].user_index();
-      else           uid = GetIndexSub(constituents[ic].phi(),constituentsunsub);
-      AliDebug(3,Form("Processing constituent %d", uid));
-      if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
-        ++gall;
+void AliEmcalJetTask::FillJetConstituents(AliEmcalJet *jet, std::vector<fastjet::PseudoJet>& constituents, TClonesArray *tracks, TClonesArray *clusters,
+                                          std::vector<fastjet::PseudoJet>& constituents_unsub, Int_t flag, TClonesArray *particles_sub)
+{
+  Int_t nt            = 0;
+  Int_t nc            = 0;
+  Double_t neutralE   = 0.;
+  Double_t maxCh      = 0.;
+  Double_t maxNe      = 0.;
+  Int_t gall          = 0;
+  Int_t gemc          = 0;
+  Int_t cemc          = 0;
+  Int_t ncharged      = 0;
+  Int_t nneutral      = 0;
+  Double_t mcpt       = 0.;
+  Double_t emcpt      = 0.;
+
+  Int_t uid   = -1;
+
+  jet->SetNumberOfTracks(constituents.size());
+  jet->SetNumberOfClusters(constituents.size());
+
+  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
+
+    if (flag == 0) uid = constituents[ic].user_index();
+    else uid = GetIndexSub(constituents[ic].phi(), constituents_unsub);
+
+    AliDebug(3,Form("Processing constituent %d", uid));
+    if (uid == -1) { //ghost particle
+      ++gall;
+      if (fGeom) {
         Double_t gphi = constituents[ic].phi();
-        if (gphi<0)
-          gphi += TMath::TwoPi();
+        if (gphi < 0) gphi += TMath::TwoPi();
         gphi *= TMath::RadToDeg();
         Double_t geta = constituents[ic].eta();
-        if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
-            (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
+        if ((gphi > fGeom->GetArm1PhiMin()) && (gphi < fGeom->GetArm1PhiMax()) &&
+            (geta > fGeom->GetArm1EtaMin()) && (geta < fGeom->GetArm1EtaMax()))
           ++gemc;
-      }        else if ((uid > 0) && fTracks) { // track constituent
-       Int_t tid = uid - 100;
-        AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
-        if (!t) {
-         AliError(Form("Could not find track %d",tid));
-          continue;
-       }
-       if (jetCount < fMarkConst) {
-         AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
-         t->SetBit(fJetType);
-       }
-        Double_t cEta = t->Eta();
-        Double_t cPhi = t->Phi();
-        Double_t cPt  = t->Pt();
-        Double_t cP   = t->P();
-        if (t->Charge() == 0) {
-          neutralE += cP;
-          ++nneutral;
-          if (cPt > maxNe)
-            maxNe = cPt;
-        } else {
-          ++ncharged;
-          if (cPt > maxCh)
-            maxCh = cPt;
-        }
-        if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
-          mcpt += cPt;
+      }
+    }  
+    else if ((uid >= 100) && tracks) { // track constituent
+      Int_t tid = uid - 100;
+      AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
+      if (!t) {
+        AliError(Form("Could not find track %d",tid));
+        continue;
+      }
 
-        if (cPhi<0)
-          cPhi += TMath::TwoPi();
+      Double_t cEta = t->Eta();
+      Double_t cPhi = t->Phi();
+      Double_t cPt  = t->Pt();
+      Double_t cP   = t->P();
+      if (t->Charge() == 0) {
+        neutralE += cP;
+        ++nneutral;
+        if (cPt > maxNe) maxNe = cPt;
+      } else {
+        ++ncharged;
+        if (cPt > maxCh) maxCh = cPt;
+      }
+
+      // check if MC particle
+      if (TMath::Abs(t->GetLabel()) > fMinMCLabel) mcpt += cPt;
+
+      if (fGeom) {
+        if (cPhi < 0) cPhi += TMath::TwoPi();
         cPhi *= TMath::RadToDeg();
-        if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
-            (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
+        if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
+            (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
           emcpt += cPt;
           ++cemc;
         }
+      }
 
-        if(flagsub!=0 && fTracksSub) {       //fill the subtracted track branch
-         if(fIsPicoTrack) 
-           new ((*fTracksSub)[tid]) AliPicoTrack(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),0,0,0,0,0,0,0,constituents[ic].m());
-         else if(fIsMcPart) { 
-           AliMCParticle *tmc = static_cast<AliMCParticle*>(fTracks->At(tid));
-           TParticle *mcPart = new TParticle(tmc->PdgCode(),tmc->Particle()->GetStatusCode(),tmc->GetMother(),0,tmc->GetFirstDaughter(),tmc->GetLastDaughter(),constituents[ic].px(),constituents[ic].py(),constituents[ic].pz(),constituents[ic].e(),0.,0.,0.,0.);
-           new ((*fTracksSub)[tid]) AliMCParticle(mcPart,0x0,tid);
-         }
-       }
+      if (flag == 0 || particles_sub == 0) {
         jet->AddTrackAt(tid, nt);
-        ++nt;
-      } else if (fClus) { // cluster constituent
-       Int_t cid = -uid - 100;
-       AliVCluster *c = 0;
-        Double_t cEta=0,cPhi=0,cPt=0,cP=0;
-        if (fIsEmcPart) {
-          AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
-          if (!ep)
-            continue;
-          c = ep->GetCluster();
-          if (!c)
-            continue;
-         if (jetCount < fMarkConst)
-           ep->SetBit(fJetType);
-          cEta = ep->Eta();
-          cPhi = ep->Phi();
-          cPt  = ep->Pt();
-          cP   = ep->P();
-        } else {
-          c = static_cast<AliVCluster*>(fClus->At(cid));
-          if (!c)
-            continue;
-         if (jetCount < fMarkConst)
-           c->SetBit(fJetType);
-          TLorentzVector nP;
-          c->GetMomentum(nP, vertex);
-          cEta = nP.Eta();
-          cPhi = nP.Phi();
-          cPt  = nP.Pt();
-          cP   = nP.P();
-        }
+      }
+      else {
+        Int_t part_sub_id = particles_sub->GetEntriesFast();
+        AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(t);
+        part_sub->SetPt(constituents[ic].perp());
+        jet->AddTrackAt(part_sub_id, nt);
+      }
 
-        neutralE += cP;
-        if (cPt > maxNe)
-          maxNe = cPt;
+      ++nt;
+    } 
+    else if ((uid <= -100) && clusters) { // cluster constituent
+      Int_t cid = -uid - 100;
+      AliVCluster *c = 0;
+      Double_t cEta=0,cPhi=0,cPt=0,cP=0;
+      if (fIsEmcPart) {
+        AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
+        if (!ep) continue;
+        c = ep->GetCluster();
+        if (!c) continue;
+
+        cEta = ep->Eta();
+        cPhi = ep->Phi();
+        cPt  = ep->Pt();
+        cP   = ep->P();
+      } 
+      else {
+        c = static_cast<AliVCluster*>(fClus->At(cid));
+        if (!c) continue;
+
+        TLorentzVector nP;
+        c->GetMomentum(nP, fVertex);
+        cEta = nP.Eta();
+        cPhi = nP.Phi();
+        cPt  = nP.Pt();
+        cP   = nP.P();
+      }
+
+      neutralE += cP;
+      if (cPt > maxNe) maxNe = cPt;
 
-        if (c->GetLabel() > fMinMCLabel) // MC particle
-          mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
+      // MC particle
+      if (TMath::Abs(c->GetLabel()) > fMinMCLabel) mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
 
-        if (cPhi<0)
-          cPhi += TMath::TwoPi();
+      if (fGeom) {
+        if (cPhi<0) cPhi += TMath::TwoPi();
         cPhi *= TMath::RadToDeg();
-        if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
-            (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
+        if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
+            (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
           emcpt += cPt;
           ++cemc;
         }
+      }
 
-       //Add constituent subtracted clusters 
-       if(flagsub!=0 && fClusSub && !fIsEmcPart) {       //fill the subtracted track branch
-         AliVCluster *oc;
-         if(fIsEsdClus) {
-           AliESDCaloCluster *ec = dynamic_cast<AliESDCaloCluster*>(c);
-           if (!ec) continue;
-           oc = new ((*fClusSub)[cid]) AliESDCaloCluster(*ec);
-         } else {
-           AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(c);
-           if (!ac) continue;
-           oc = new ((*fClusSub)[cid]) AliAODCaloCluster(*ac);
-         }
-         oc->SetE(constituents[ic].E());
-       }
-
+      if (flag == 0 || particles_sub == 0) {
         jet->AddClusterAt(cid, nc);
-        ++nc;
-        ++nneutral;
-      } else {
-        AliError(Form("%s: No logical way to end up here.", GetName()));
-        continue;
       }
+      else {
+        Int_t part_sub_id = particles_sub->GetEntriesFast();
+        AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(c);
+        part_sub->SetPt(constituents[ic].perp());
+        jet->AddTrackAt(part_sub_id, nt);
+      }
+
+      ++nc;
+      ++nneutral;
+    } 
+    else {
+      AliError(Form("%s: No logical way to end up here.", GetName()));
+      continue;
     }
+  }
+
+  jet->SetNumberOfTracks(nt);
+  jet->SetNumberOfClusters(nc);
+  jet->SetNEF(neutralE / jet->E());
+  jet->SetMaxChargedPt(maxCh);
+  jet->SetMaxNeutralPt(maxNe);
+  if (gall > 0) jet->SetAreaEmc(jet->Area() * gemc / gall);
+  else jet->SetAreaEmc(-1);
+  jet->SetNEmc(cemc);
+  jet->SetNumberOfCharged(ncharged);
+  jet->SetNumberOfNeutrals(nneutral);
+  jet->SetMCPt(mcpt);
+  jet->SetPtEmc(emcpt);
+  jet->SortConstituents();
 }
+
 //______________________________________________________________________________________
-Int_t AliEmcalJetTask::GetIndexSub(Double_t phisub,std::vector<fastjet::PseudoJet>& constituentsunsub) {
-  for(UInt_t ii=0;ii<constituentsunsub.size();ii++) {
-    Double_t phiunsub=constituentsunsub[ii].phi();
-    if(phisub==phiunsub) return constituentsunsub[ii].user_index();
-  } 
+Int_t AliEmcalJetTask::GetIndexSub(Double_t phi_sub, std::vector<fastjet::PseudoJet>& constituents_unsub) 
+{
+  for (UInt_t ii = 0; ii < constituents_unsub.size(); ii++) {
+    Double_t phi_unsub = constituents_unsub[ii].phi();
+    if (phi_sub == phi_unsub) return constituents_unsub[ii].user_index();
+  }
+
   return -1;
 }
-//______________________________________________________________________________________
index 69b2a58..56e5fd0 100644 (file)
@@ -2,8 +2,11 @@
 #define ALIEMCALJETTASK_H
 
 class TClonesArray;
+class TObjArray;
 class AliVEvent;
-class AliRhoParameter;
+class AliEMCALGeometry;
+class AliEmcalJetUtility;
+class AliFJWrapper;
 
 namespace fastjet {
   class PseudoJet;
@@ -37,40 +40,83 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
   AliEmcalJetTask(const char *name);
   virtual ~AliEmcalJetTask();
 
-  void                   UserCreateOutputObjects();
-  void                   UserExec(Option_t *option);
-  void                   Terminate(Option_t *option);
-  Bool_t                 IsLocked()                       { if(fLocked) {AliFatal("Jet finder task is locked! Changing properties is not allowed."); return kTRUE;} else return kFALSE;};
-  void                   SetLocked()                      { fLocked = kTRUE;}
-  void                   SelectConstituents(UInt_t constSel, UInt_t MCconstSel)  { fConstSel = constSel; fMCConstSel = MCconstSel; };
-  void                   SetAlgo(Int_t a)                 { if(IsLocked()) return; if (a==0) fJetType |= kKT; else fJetType |= kAKT; }  // for backward compatibility only
-  void                   SetClusName(const char *n)       { if(IsLocked()) return; fCaloName      = n     ; }
-  void                   SetJetsName(const char *n)       { if(IsLocked()) return; fJetsName      = n     ; }
-  void                   SetJetsSubName(const char *n)    { fJetsSubName   = n     ; }
-  void                   SetTracksSubName(const char *n)  { fTracksSubName = n     ; }
-  void                   SetClusSubName(const char *n)    { fClusSubName   = n     ; }
-  void                   SetJetType(UInt_t t)             { if(IsLocked()) return; fJetType       = t     ; }
-  void                   SetMarkConstituents(UInt_t m)    { if(IsLocked()) return; fMarkConst     = m     ; }
-  void                   SetMinJetArea(Double_t a)        { if(IsLocked()) return; fMinJetArea    = a     ; }
-  void                   SetMinJetClusPt(Double_t min)    { if(IsLocked()) return; fMinJetClusPt  = min   ; }
-  void                   SetMinJetPt(Double_t j)          { if(IsLocked()) return; fMinJetPt      = j     ; }
-  void                   SetMinJetTrackPt(Double_t min)   { if(IsLocked()) return; fMinJetTrackPt = min   ; }
-  void                   SetRadius(Double_t r)            { if(IsLocked()) return; fRadius        = r     ; if ((fJetType & (kRX1Jet|kRX2Jet|kRX3Jet)) == 0) AliWarning("Radius value will be ignored if jet type is not set to a user defined radius (kRX1Jet,kRX2Jet,kRX3Jet)."); }
-  void                   SetTrackEfficiency(Double_t t)   { if(IsLocked()) return; fTrackEfficiency = t   ; }
-  void                   SetTracksName(const char *n)     { if(IsLocked()) return; fTracksName    = n     ; }
-  void                   SetType(Int_t t)                 { if(IsLocked()) return; 
-                                                            if (t==0) fJetType |= kFullJet; 
-                                                            else if (t==1) fJetType |= kChargedJet; 
-                                                            else if (t==2) fJetType |= kNeutralJet; } // for backward compatibility only
-  void                   SetEtaRange(Double_t emi, Double_t ema) {if(IsLocked()) return; fEtaMin = emi; fEtaMax = ema; }
-  void                   SetPhiRange(Double_t pmi, Double_t pma) {if(IsLocked()) return; fPhiMin = pmi; fPhiMax = pma; }
-  void                   SetJetEtaRange(Double_t emi, Double_t ema) {if(IsLocked()) return; fJetEtaMin = emi; fJetEtaMax = ema; }
-  void                   SetJetPhiRange(Double_t pmi, Double_t pma) {if(IsLocked()) return; fJetPhiMin = pmi; fJetPhiMax = pma; }
-  void                   SetGhostArea(Double_t gharea)    { if(IsLocked()) return; fGhostArea      = gharea;  }
-  void                   SetMinMCLabel(Int_t s)           { if(IsLocked()) return; fMinMCLabel     = s     ;  }
-  void                   SetRecombScheme(Int_t scheme)    { if(IsLocked()) return; fRecombScheme   = scheme;  }
-  void                   SelectCollisionCandidates(UInt_t offlineTriggerMask = AliVEvent::kMB)
+  void UserCreateOutputObjects();
+  void UserExec(Option_t* option);
+  void Terminate(Option_t* /*option*/) {;}
+
+  void                   SetAlgo(Int_t a)                           { if (IsLocked()) return; if (a==0) fJetType |= kKT; else fJetType |= kAKT; }  // for backward compatibility only
+  void                   SetClusName(const char *n)                 { if (IsLocked()) return; fCaloName         = n     ; }
+  void                   SetClusLabelRange(Int_t min, Int_t max)    { if (IsLocked()) return; fMinLabelClusters = min   ; fMaxLabelClusters = max; }
+  void                   SetEtaRange(Double_t emi, Double_t ema)    { if (IsLocked()) return; fEtaMin           = emi   ; fEtaMax = ema   ; }
+  void                   SetGhostArea(Double_t gharea)              { if (IsLocked()) return; fGhostArea        = gharea; }
+  void                   SetJetsName(const char *n)                 { if (IsLocked()) return; fJetsName         = n     ; }
+  void                   SetJetEtaRange(Double_t emi, Double_t ema) { if (IsLocked()) return; fJetEtaMin        = emi   ; fJetEtaMax = ema; }
+  void                   SetJetPhiRange(Double_t pmi, Double_t pma) { if (IsLocked()) return; fJetPhiMin        = pmi   ; fJetPhiMax = pma; }
+  void                   SetJetType(UInt_t t)                       { if (IsLocked()) return; fJetType          = t     ; }
+  void                   SetLocked()                                { fLocked = kTRUE;}
+  void                   SetMinJetArea(Double_t a)                  { if (IsLocked()) return; fMinJetArea       = a     ; }
+  void                   SetMinJetClusPt(Double_t min)              { if (IsLocked()) return; fMinJetClusPt     = min   ; }
+  void                   SetMinJetPt(Double_t j)                    { if (IsLocked()) return; fMinJetPt         = j     ; }
+  void                   SetMinJetTrackPt(Double_t min)             { if (IsLocked()) return; fMinJetTrackPt    = min   ; }
+  void                   SetMinMCLabel(Int_t s)                     { if (IsLocked()) return; fMinMCLabel       = s     ; }
+  void                   SetPhiRange(Double_t pmi, Double_t pma)    { if (IsLocked()) return; fPhiMin           = pmi   ; fPhiMax = pma; }
+  void                   SetRecombScheme(Int_t scheme)              { if (IsLocked()) return; fRecombScheme     = scheme; }
+  void                   SetTracksName(const char *n)               { if (IsLocked()) return; fTracksName       = n     ; }
+  void                   SetTrackEfficiency(Double_t t)             { if (IsLocked()) return; fTrackEfficiency  = t     ; }
+  void                   SetTrackLabelRange(Int_t min, Int_t max)   { if (IsLocked()) return; fMinLabelTracks   = min   ; fMaxLabelTracks = max; }
+  void                   SetLegacyMode(Bool_t mode)                 { if (IsLocked()) return; fLegacyMode       = mode  ; }
+  void                   SetMCFlag(UInt_t m)                        { if (IsLocked()) return; fMCFlag           = m     ; }
+  void                   SelectHIJING(Bool_t s)                     { if (IsLocked()) return; if (s) fGeneratorIndex = 0; else fGeneratorIndex = -1; }
+  void                   SetGeneratorIndex(Short_t i)               { if (IsLocked()) return; fGeneratorIndex   = i     ; }
+
+  AliEmcalJetUtility*    AddUtility(AliEmcalJetUtility* utility);
+
+  const char*            GetClusName()                    { return fCaloName.Data()   ; }
+  Double_t               GetEtaMin()                      { return fEtaMin            ; }
+  Double_t               GetEtaMax()                      { return fEtaMax            ; }
+  Double_t               GetGhostArea()                   { return fGhostArea         ; }
+  const char*            GetJetsName()                    { return fJetsName.Data()   ; }
+  Double_t               GetJetEtaMin()                   { return fJetEtaMin         ; }
+  Double_t               GetJetEtaMax()                   { return fJetEtaMax         ; }
+  Double_t               GetJetPhiMin()                   { return fJetPhiMin         ; }
+  Double_t               GetJetPhiMax()                   { return fJetPhiMax         ; }
+  UInt_t                 GetJetType()                     { return fJetType           ; }
+  Bool_t                 GetLegacyMode()                  { return fLegacyMode        ; }
+  Double_t               GetMinJetArea()                  { return fMinJetArea        ; }
+  Double_t               GetMinJetClusPt()                { return fMinJetClusPt      ; }
+  Double_t               GetMinJetPt()                    { return fMinJetPt          ; }
+  Double_t               GetMinJetTrackPt()               { return fMinJetTrackPt     ; }
+  Int_t                  GetMinMCLabel()                  { return fMinMCLabel        ; }
+  Double_t               GetPhiMin()                      { return fPhiMin            ; }
+  Double_t               GetPhiMax()                      { return fPhiMax            ; }
+  Double_t               GetRadius()                      { return fRadius            ; }
+  Int_t                  GetRecombScheme()                { return fRecombScheme      ; }
+  const char*            GetTracksName()                  { return fTracksName.Data() ; }
+  Double_t               GetTrackEfficiency()             { return fTrackEfficiency   ; }
+
+  AliVEvent*             GetEvent()                       { return fEvent             ; }
+  TClonesArray*          GetClusters()                    { return fClus              ; }
+  TClonesArray*          GetTracks()                      { return fTracks            ; }
+  TClonesArray*          GetJets()                        { return fJets              ; }
+  TObjArray*             GetUtilities()                   { return fUtilities         ; }
+
+  void                   FillJetConstituents(AliEmcalJet *jet, std::vector<fastjet::PseudoJet>& constituents, TClonesArray *tracks, TClonesArray *clusters,
+                                             std::vector<fastjet::PseudoJet>& constituents_sub, Int_t flag = 0, TClonesArray *particles_sub = 0);
+
+  Int_t                  GetIndexSub(Double_t phi_sub, std::vector<fastjet::PseudoJet>& constituents_unsub);
+
+  Bool_t IsLocked() 
+  {
+    if (fLocked) {
+      AliFatal("Jet finder task is locked! Changing properties is not allowed."); 
+      return kTRUE;
+    } 
+    else {
+      return kFALSE;
+    }
+  }
+
+  void SelectCollisionCandidates(UInt_t offlineTriggerMask = AliVEvent::kMB)
   {
     if(!fIsPSelSet)
     {
@@ -83,78 +129,52 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
       fOfflineTriggerMask = fOfflineTriggerMask | offlineTriggerMask;
     }
   }
-  void                   SetLegacyMode(Bool_t mode)           { if(IsLocked()) return; fLegacyMode = mode; }
-  void                   SetMCFlag(UInt_t m)                  { if(IsLocked()) return; fMCFlag     = m   ; }
-  void                   SelectHIJING(Bool_t s)               { if(IsLocked()) return; if (s) fGeneratorIndex = 0; else fGeneratorIndex = -1; }
-  void                   SetGeneratorIndex(Short_t i)         { if(IsLocked()) return; fGeneratorIndex = i; }
-  void                   SelectPhysicalPrimaries(Bool_t s)    { if(IsLocked()) return; 
-                                                                if (s) fMCFlag |=  AliAODMCParticle::kPhysicalPrim ; 
-                                                                else   fMCFlag &= ~AliAODMCParticle::kPhysicalPrim ; }
-
-  void                   SetCodeDebug(Bool_t val)               { fCodeDebug    = val          ; }
-  void                   SetForceIsMcPart(Bool_t b)             { fIsMcPart     = b            ; }
-  void                   SetPionMassForClusters(Bool_t b)       { fPionMassClusters = b        ; }
-
-  void                   SetRhoName(const char *n)              { fRhoName      = n            ; }
-  void                   SetRhomName(const char *n)             { fRhomName     = n            ; }
-  void                   SetGenericSubtractionJetMass(Bool_t b) { fDoGenericSubtractionJetMass = b; }
-  void                   SetGenericSubtractionGR(Bool_t b, Double_t rmax = 2., Double_t dr = 0.04, Double_t ptmin = 0.) { fDoGenericSubtractionGR = b; fRMax=rmax; fDRStep=dr; fPtMinGR=ptmin;}
-  void                   SetConstituentSubtraction(Bool_t b)    { fDoConstituentSubtraction = b; }
-  void                   SetGenericSubtractionExtraJetShapes(Bool_t b)            { fDoGenericSubtractionExtraJetShapes =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; }
-  const char*            GetTracksName()                  { return fTracksName.Data(); }
-  const char*            GetClusName()                    { return fCaloName.Data(); }
-  const char*            GetRhoName()                     { return fRhoName.Data(); }
-  UInt_t                 GetMarkConstituents()            { return fMarkConst; }
-  Double_t               GetMinJetArea()                  { return fMinJetArea; }
-  Double_t               GetMinJetClusPt()                { return fMinJetClusPt; }
-  Double_t               GetMinJetPt()                    { return fMinJetPt; }
-  Double_t               GetMinJetTrackPt()               { return fMinJetTrackPt; }
-  Double_t               GetTrackEfficiency()             { return fTrackEfficiency; }
-  Double_t               GetGhostArea()                   { return fGhostArea; }
-  Int_t                  GetMinMCLabel()                  { return fMinMCLabel; }
-  Int_t                  GetRecombScheme()                { return fRecombScheme; }
-  Double_t               GetEtaMin()                      { return fEtaMin; }
-  Double_t               GetEtaMax()                      { return fEtaMax; }
-  Double_t               GetPhiMin()                      { return fPhiMin; }
-  Double_t               GetPhiMax()                      { return fPhiMax; }
-  Double_t               GetJetEtaMin()                   { return fJetEtaMin; }
-  Double_t               GetJetEtaMax()                   { return fJetEtaMax; }
-  Double_t               GetJetPhiMin()                   { return fJetPhiMin; }
-  Double_t               GetJetPhiMax()                   { return fJetPhiMax; }
-
-  AliVEvent*             GetEvent()                       { return fEvent;}
-  TClonesArray*          GetClusters()                    { return fClus;}
-  TClonesArray*          GetTracks()                      { return fTracks;}
-  TClonesArray*          GetJets()                        { return fJets;}
+
+  void SetRadius(Double_t r)            
+  { 
+    if (IsLocked()) return; 
+    fRadius = r; 
+    if ((fJetType & (kRX1Jet|kRX2Jet|kRX3Jet)) == 0) 
+      AliWarning("Radius value will be ignored if jet type is not set to a user defined radius (kRX1Jet,kRX2Jet,kRX3Jet)."); 
+  }
+
+  void SetType(Int_t t)                 
+  { 
+    if(IsLocked()) return; 
+    if (t==0) fJetType |= kFullJet; 
+    else if (t==1) fJetType |= kChargedJet; 
+    else if (t==2) fJetType |= kNeutralJet; 
+  } // for backward compatibility only
+
+  void SelectPhysicalPrimaries(Bool_t s)    
+  { 
+    if(IsLocked()) return; 
+    if (s) fMCFlag |=  AliAODMCParticle::kPhysicalPrim; 
+    else   fMCFlag &= ~AliAODMCParticle::kPhysicalPrim; 
+  }
 
  protected:
+
   void                   FindJets();
+  void                   FillJetBranch();
   Bool_t                 DoInit();
+  void                   InitUtilities();
+  void                   PrepareUtilities();
+  void                   ExecuteUtilities(AliEmcalJet* jet, Int_t ij);
+  void                   TerminateUtilities();
   Bool_t                 GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const;
 
-   void                   FillJetConstituents(std::vector<fastjet::PseudoJet>& constituents,AliEmcalJet *jet,Double_t vertex[3],UInt_t jetCount,Int_t& nt,Int_t& nc,Double_t& maxCh,Double_t& maxNe,Int_t& ncharged,Int_t& nneutral,Double_t& neutralE,Double_t& mcpt,Int_t& cemc,Double_t& emcpt,Int_t& gall,Int_t& gemc,std::vector<fastjet::PseudoJet>& constituentsunsub,Int_t flag); 
-  Int_t                  GetIndexSub(Double_t phisub,std::vector<fastjet::PseudoJet>& constituentsunsub);
-
   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
-  TString                fTracksSubName;          // name of subtracted track collection
-  TString                fClusSubName;            // name of subtracted clusters 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
-  UInt_t                 fMarkConst;              // constituents are marked (via TObject::SetBit) as belonging to the # leading jets
+  Int_t                  fMinLabelTracks;         // select track constituents with a minimum label index (track->GetLabel())
+  Int_t                  fMaxLabelTracks;         // select track constituents with a maximum label index (track->GetLabel())
+  Int_t                  fMinLabelClusters;       // select cluster constituents with a minimum label index (cluster->GetLabel())
+  Int_t                  fMaxLabelClusters;       // select cluster constituents with a maximum label index (cluster->GetLabel())
+  Int_t                  fMinMCLabel;             // minimum MC label value for the tracks/clusters being considered MC particles
   Double_t               fRadius;                 // jet radius
-  Double_t               fMinJetTrackPt;          // min jet track momentum   (applied before clustering)
+  Double_t               fMinJetTrackPt;          // min jet track momentum (applied before clustering)
   Double_t               fMinJetClusPt;           // min jet cluster momentum (applied before clustering)
   Double_t               fPhiMin;                 // minimum phi for constituents (applied before clustering)
   Double_t               fPhiMax;                 // maximum phi for constituents (applied before clustering)
@@ -167,44 +187,25 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
   Double_t               fJetEtaMin;              // minimum eta to keep jet in output
   Double_t               fJetEtaMax;              // maximum eta to keep jet in output
   Double_t               fGhostArea;              // ghost area
-  Int_t                  fMinMCLabel;             // minimum MC label value for the tracks/clusters being considered MC particles
   Int_t                  fRecombScheme;           // recombination scheme used by fastjet
   Double_t               fTrackEfficiency;        // artificial tracking inefficiency (0...1)
-  UInt_t                 fMCFlag;                 // select MC particles with flags
-  Short_t                fGeneratorIndex;         // select MC particles with generator index (default = -1 = switch off selection)
-  Bool_t                 fIsInit;                 //!=true if already initialized
+  UInt_t                 fMCFlag;                 // select MC particles with flags (e.g. to select physical primaries)
+  Short_t                fGeneratorIndex;         // select MC particles with generator index (default = -1 to switch off selection)
+  TObjArray             *fUtilities;              // jet utilities (gen subtractor, constituent subtractor etc.)
   Bool_t                 fLocked;                 // true if lock is set
+  Bool_t                 fIsInit;                 //!=true if already initialized
   Bool_t                 fIsPSelSet;              //!=true if physics selection was set
-  Bool_t                 fIsMcPart;               //!=true if MC particles are given as input
   Bool_t                 fIsEmcPart;              //!=true if emcal particles are given as input (for clusters)
-  Bool_t                 fIsPicoTrack;            //!=true if pico tracks are given as input
-  Bool_t                 fIsEsdClus;              //!=true if AliESDCaloCluster is given as input
-  Bool_t                 fLegacyMode;             //! if true, enable FJ 2.x behavior
-  Bool_t                 fCodeDebug;              // use nontested code changes 
-  Bool_t                 fPionMassClusters;       // assume pion mass for clusters
-  Bool_t                 fDoGenericSubtractionJetMass;        // calculate generic subtraction
-  Bool_t                 fDoGenericSubtractionGR;             // calculate generic subtraction for angular structure function GR
-  Bool_t                 fDoGenericSubtractionExtraJetShapes; // calculate generic subtraction for other jet shapes like radialmoment,pTD etc
-  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
-  Double_t               fRMax;                               // R max for GR calculation
-  Double_t               fDRStep;                             // step width for GR calculation
-  Double_t               fPtMinGR;                            // min pT for GR calculation
+  Bool_t                 fLegacyMode;             //!=true to enable FJ 2.x behavior
 
+  Double_t               fVertex[3];              //!vertex of the current event
+  AliEMCALGeometry      *fGeom;                   //!emcal geometry object
   TClonesArray          *fJets;                   //!jet collection
-  TClonesArray          *fJetsSub;                //!subtracted jet collection
-  TClonesArray          *fTracksSub;              //!subtracted track collection
-  TClonesArray          *fClusSub;                //!subtracted cluster collection
   AliVEvent             *fEvent;                  //!current event
   TClonesArray          *fTracks;                 //!tracks collection
   TClonesArray          *fClus;                   //!cluster collection
-  AliRhoParameter       *fRhoParam;               //!event rho
-  AliRhoParameter       *fRhomParam;              //!event rhom
+  AliFJWrapper           fFastJetWrapper;         //!fastjet wrapper
 
  private:
   AliEmcalJetTask(const AliEmcalJetTask&);            // not implemented
diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtility.cxx b/PWGJE/EMCALJetTasks/AliEmcalJetUtility.cxx
new file mode 100644 (file)
index 0000000..35e53e5
--- /dev/null
@@ -0,0 +1,42 @@
+#include "AliEmcalJetUtility.h"
+
+ClassImp(AliEmcalJetUtility)
+
+//______________________________________________________________________________
+AliEmcalJetUtility::AliEmcalJetUtility() :
+TNamed(),
+  fJetTask(0),
+  fInit(kFALSE)
+{
+  // Dummy constructor.
+
+}
+
+//______________________________________________________________________________
+AliEmcalJetUtility::AliEmcalJetUtility(const char* name) :
+  TNamed(name, name),
+  fJetTask(0),
+  fInit(kFALSE)
+{
+  // Default constructor.
+}
+
+//______________________________________________________________________________
+AliEmcalJetUtility::AliEmcalJetUtility(const AliEmcalJetUtility &other) :
+  TNamed(other),
+  fJetTask(other.fJetTask),
+  fInit(other.fInit)
+{
+  // Copy constructor.
+}
+
+//______________________________________________________________________________
+AliEmcalJetUtility& AliEmcalJetUtility::operator=(const AliEmcalJetUtility &other)
+{
+  // Assignment
+  if (&other == this) return *this;
+  TNamed::operator=(other);
+  fJetTask = other.fJetTask;
+  fInit = other.fInit;
+  return *this;
+}
diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtility.h b/PWGJE/EMCALJetTasks/AliEmcalJetUtility.h
new file mode 100644 (file)
index 0000000..3283afd
--- /dev/null
@@ -0,0 +1,36 @@
+#ifndef ALIEMCALJETUTILITY_H
+#define ALIEMCALJETUTILITY_H
+
+#include <TNamed.h>
+
+#include "AliFJWrapper.h"
+
+class AliEmcalJetTask;
+class AliEmcalJet;
+class AliFJWrapper;
+
+class AliEmcalJetUtility : public TNamed
+{
+ public:
+
+  AliEmcalJetUtility();
+  AliEmcalJetUtility(const char* name);
+  AliEmcalJetUtility(const AliEmcalJetUtility &jet);
+  AliEmcalJetUtility& operator=(const AliEmcalJetUtility &jet);
+  ~AliEmcalJetUtility() {;}
+
+  virtual void Init() = 0;                                                        // Executed only once at the end of AliEmcalJetTask::DoInit()
+  virtual void Prepare(AliFJWrapper& fjw) = 0;                                    // Executed for each event at the beginning of AliEmcalJetTask::FillJetBranch()
+  virtual void ProcessJet(AliEmcalJet* jet, Int_t ij, AliFJWrapper& fjw) = 0;     // Executed for each jet in the loop in AliEmcalJetTask::FillJetBranch()
+  virtual void Terminate(AliFJWrapper& fjw) = 0;                                  // Executed for each event at the end of AliEmcalJetTask::FillJetBranch()
+
+  void SetJetTask(AliEmcalJetTask* jetTask) { fJetTask = jetTask; }
+
+ protected:
+
+  AliEmcalJetTask       *fJetTask     ; // pointer to the main jet task
+  Bool_t                 fInit        ; //! whether or not the utility has been initialized
+
+  ClassDef(AliEmcalJetUtility, 1) // Abstract Emcal jet utility class
+};
+#endif
diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.cxx b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.cxx
new file mode 100644 (file)
index 0000000..e46568f
--- /dev/null
@@ -0,0 +1,195 @@
+#include "AliEmcalJetUtilityConstSubtractor.h"
+#include "AliEmcalJet.h"
+#include "AliRhoParameter.h"
+#include "AliEmcalJetTask.h"
+
+ClassImp(AliEmcalJetUtilityConstSubtractor)
+
+//______________________________________________________________________________
+AliEmcalJetUtilityConstSubtractor::AliEmcalJetUtilityConstSubtractor() :
+  AliEmcalJetUtility(),
+  fJetsSubName(""),
+  fParticlesSubName(""),
+  fUseExternalBkg(kFALSE),
+  fRhoName(""),
+  fRhomName(""),
+  fRho(0),
+  fRhom(0),
+  fJetsSub(0x0),
+  fParticlesSub(0x0),
+  fRhoParam(0),
+  fRhomParam(0)
+{
+  // Dummy constructor.
+
+}
+
+//______________________________________________________________________________
+AliEmcalJetUtilityConstSubtractor::AliEmcalJetUtilityConstSubtractor(const char* name) :
+  AliEmcalJetUtility(name),
+  fJetsSubName(""),
+  fParticlesSubName(""),
+  fUseExternalBkg(kFALSE),
+  fRhoName(""),
+  fRhomName(""),
+  fRho(0),
+  fRhom(0),
+  fJetsSub(0x0),
+  fParticlesSub(0x0),
+  fRhoParam(0),
+  fRhomParam(0)
+{
+  // Default constructor.
+}
+
+//______________________________________________________________________________
+AliEmcalJetUtilityConstSubtractor::AliEmcalJetUtilityConstSubtractor(const AliEmcalJetUtilityConstSubtractor &other) :
+  AliEmcalJetUtility(other),
+  fJetsSubName(other.fJetsSubName),
+  fParticlesSubName(other.fParticlesSubName),
+  fUseExternalBkg(other.fUseExternalBkg),
+  fRhoName(other.fRhoName),
+  fRhomName(other.fRhomName),
+  fRho(other.fRho),
+  fRhom(other.fRhom),
+  fJetsSub(other.fJetsSub),
+  fParticlesSub(other.fParticlesSub),
+  fRhoParam(other.fRhoParam),
+  fRhomParam(other.fRhomParam)
+{
+  // Copy constructor.
+}
+
+//______________________________________________________________________________
+AliEmcalJetUtilityConstSubtractor& AliEmcalJetUtilityConstSubtractor::operator=(const AliEmcalJetUtilityConstSubtractor &other)
+{
+  // Assignment.
+
+  if (&other == this) return *this;
+  AliEmcalJetUtility::operator=(other);
+  fJetsSubName = other.fJetsSubName;
+  fParticlesSubName = other.fParticlesSubName;
+  fUseExternalBkg = other.fUseExternalBkg;
+  fRhoName = other.fRhoName;
+  fRhomName = other.fRhomName;
+  fRho = other.fRho;
+  fRhom = other.fRhom;
+  fJetsSub = other.fJetsSub;
+  fParticlesSub = other.fParticlesSub;
+  fRhoParam = other.fRhoParam;
+  fRhomParam = other.fRhomParam;
+  return *this;
+}
+
+//______________________________________________________________________________
+void AliEmcalJetUtilityConstSubtractor::Init()
+{
+  // Initialize the utility.
+
+  // Add constituent subtracted jets to event
+  if (!fJetsSubName.IsNull()) {
+    if (!(fJetTask->GetEvent()->FindListObject(fJetsSubName)) ) {
+      fJetsSub = new TClonesArray("AliEmcalJet");
+      fJetsSub->SetName(fJetsSubName);
+      fJetTask->GetEvent()->AddObject(fJetsSub);
+    } 
+    else {
+      AliError(Form("%s: Object for subtracted jet branch with name %s already in event! Returning", GetName(), fJetsSubName.Data()));
+      return;
+    }
+  }
+
+  // Add tracks from constituent subtracted jets to event
+  if (!fParticlesSubName.IsNull()) {
+    if (!(fJetTask->GetEvent()->FindListObject(fParticlesSubName))) {
+      fParticlesSub = new TClonesArray("AliEmcalParticle");
+      fParticlesSub->SetName(fParticlesSubName);
+      fJetTask->GetEvent()->AddObject(fParticlesSub);
+    } else {
+      AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fParticlesSubName.Data()));
+      return;
+    }
+  }
+
+  if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event
+    fRhoParam = dynamic_cast<AliRhoParameter*>(fJetTask->GetEvent()->FindListObject(fRhoName));
+    if (!fRhoParam) {
+      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
+      return;
+    }
+  }
+
+  if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event
+    fRhomParam = dynamic_cast<AliRhoParameter*>(fJetTask->GetEvent()->FindListObject(fRhomName));
+    if (!fRhomParam) {
+      AliError(Form("%s: Could not retrieve rho_m %s!", GetName(), fRhomName.Data()));
+      return;
+    }
+  }
+
+  fInit = kTRUE;
+}
+
+//______________________________________________________________________________
+void AliEmcalJetUtilityConstSubtractor::Prepare(AliFJWrapper& fjw)
+{
+  // Prepare the utility.
+
+  if (!fInit) return;
+
+  if (fRhoParam) fRho = fRhoParam->GetVal();
+  if (fRhomParam) fRhom = fRhomParam->GetVal();
+
+  if (fJetsSub) fJetsSub->Delete();
+
+  fjw.SetUseExternalBkg(fUseExternalBkg, fRho, fRhom);
+  fjw.DoConstituentSubtraction();
+}
+
+//______________________________________________________________________________
+void AliEmcalJetUtilityConstSubtractor::ProcessJet(AliEmcalJet* /*jet*/, Int_t /*ij*/, AliFJWrapper& /*fjw*/)
+{
+  // Proceess each jet.
+}
+
+//______________________________________________________________________________
+void AliEmcalJetUtilityConstSubtractor::Terminate(AliFJWrapper& fjw)
+{
+  // Run termination of the utility (after each event).
+
+  if (!fInit) return;
+
+  if (!fJetsSub) {
+    AliWarning(Form("No jet branch to write to for subtracted jets. fJetsSubName: %s", fJetsSubName.Data()));
+    return;
+  }
+
+#ifdef FASTJET_VERSION
+  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());
+
+    // Fill constituent info
+    std::vector<fastjet::PseudoJet> constituents_unsub(fjw.GetJetConstituents(ijet));
+    std::vector<fastjet::PseudoJet> constituents_sub = jets_sub[ijet].constituents();
+    jet_sub->SetNumberOfTracks(constituents_sub.size());
+    fJetTask->FillJetConstituents(jet_sub, constituents_sub, fJetTask->GetTracks(), fJetTask->GetClusters(), constituents_unsub, 1, fParticlesSub);
+
+    jetCount++;
+  }
+
+#endif
+
+}
+
diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.h b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityConstSubtractor.h
new file mode 100644 (file)
index 0000000..2d7dc77
--- /dev/null
@@ -0,0 +1,53 @@
+#ifndef ALIEMCALJETUTILITYCONSTSUBTRACTOR_H
+#define ALIEMCALJETUTILITYCONSTSUBTRACTOR_H
+
+#include <TNamed.h>
+
+#include "AliEmcalJetUtility.h"
+#include "AliFJWrapper.h"
+
+class AliEmcalJetTask;
+class AliEmcalJet;
+class AliFJWrapper;
+class AliRhoParameter;
+
+class AliEmcalJetUtilityConstSubtractor : public AliEmcalJetUtility
+{
+ public:
+
+  AliEmcalJetUtilityConstSubtractor();
+  AliEmcalJetUtilityConstSubtractor(const char* name);
+  AliEmcalJetUtilityConstSubtractor(const AliEmcalJetUtilityConstSubtractor &jet);
+  AliEmcalJetUtilityConstSubtractor& operator=(const AliEmcalJetUtilityConstSubtractor &jet);
+  ~AliEmcalJetUtilityConstSubtractor() {;}
+
+  void                   SetRhoName(const char *n)           { fRhoName      = n         ; }
+  void                   SetRhomName(const char *n)          { fRhomName     = n         ; }
+  void                   SetUseExternalBkg(Bool_t b)         { fUseExternalBkg   = b     ; }
+
+  void                   SetJetsSubName(const char *n)       { fJetsSubName      = n     ; }
+  void                   SetParticlesSubName(const char *n)  { fParticlesSubName = n     ; }
+
+  void Init();
+  void Prepare(AliFJWrapper& fjw);
+  void ProcessJet(AliEmcalJet* jet, Int_t ij, AliFJWrapper& fjw);
+  void Terminate(AliFJWrapper& fjw);
+
+ protected:
+
+  TString                fJetsSubName;                        // name of subtracted jet collection
+  TString                fParticlesSubName;                   // name of subtracted particle collection
+  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          *fJetsSub;                            //!subtracted jet collection
+  TClonesArray          *fParticlesSub;                       //!subtracted particle collection
+  AliRhoParameter       *fRhoParam;                           //!event rho
+  AliRhoParameter       *fRhomParam;                          //!event rhom
+
+  ClassDef(AliEmcalJetUtilityConstSubtractor, 1) // Emcal jet utility that implements the constituent subtractor form the fastjet contrib
+};
+#endif
diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.cxx b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.cxx
new file mode 100644 (file)
index 0000000..68f9adf
--- /dev/null
@@ -0,0 +1,256 @@
+#include "AliEmcalJetUtilityGenSubtractor.h"
+#include "AliEmcalJet.h"
+#include "AliRhoParameter.h"
+#include "AliEmcalJetTask.h"
+
+ClassImp(AliEmcalJetUtilityGenSubtractor)
+
+//______________________________________________________________________________
+AliEmcalJetUtilityGenSubtractor::AliEmcalJetUtilityGenSubtractor() :
+AliEmcalJetUtility(),
+  fDoGenericSubtractionJetMass(kFALSE),
+  fDoGenericSubtractionGR(kFALSE),
+  fDoGenericSubtractionExtraJetShapes(kFALSE),
+  fUseExternalBkg(kFALSE),
+  fRhoName(""),
+  fRhomName(""),
+  fRho(0),
+  fRhom(0),
+  fRMax(0.4),
+  fDRStep(0.04),
+  fPtMinGR(40.),
+  fRhoParam(0),
+  fRhomParam(0)
+{
+  // Dummy constructor.
+
+}
+
+//______________________________________________________________________________
+AliEmcalJetUtilityGenSubtractor::AliEmcalJetUtilityGenSubtractor(const char* name) :
+  AliEmcalJetUtility(name),
+  fDoGenericSubtractionJetMass(kFALSE),
+  fDoGenericSubtractionGR(kFALSE),
+  fDoGenericSubtractionExtraJetShapes(kFALSE),
+  fUseExternalBkg(kFALSE),
+  fRhoName(""),
+  fRhomName(""),
+  fRho(0),
+  fRhom(0),
+  fRMax(0.4),
+  fDRStep(0.04),
+  fPtMinGR(40.),
+  fRhoParam(0),
+  fRhomParam(0)
+{
+  // Default constructor.
+}
+
+//______________________________________________________________________________
+AliEmcalJetUtilityGenSubtractor::AliEmcalJetUtilityGenSubtractor(const AliEmcalJetUtilityGenSubtractor &other) :
+  AliEmcalJetUtility(other),
+  fDoGenericSubtractionJetMass(other.fDoGenericSubtractionJetMass),
+  fDoGenericSubtractionGR(other.fDoGenericSubtractionGR),
+  fDoGenericSubtractionExtraJetShapes(other.fDoGenericSubtractionExtraJetShapes),
+  fUseExternalBkg(other.fUseExternalBkg),
+  fRhoName(other.fRhoName),
+  fRhomName(other.fRhomName),
+  fRho(other.fRho),
+  fRhom(other.fRhom),
+  fRMax(other.fRMax),
+  fDRStep(other.fDRStep),
+  fPtMinGR(other.fPtMinGR),
+  fRhoParam(other.fRhoParam),
+  fRhomParam(other.fRhomParam)
+{
+  // Copy constructor.
+}
+
+//______________________________________________________________________________
+AliEmcalJetUtilityGenSubtractor& AliEmcalJetUtilityGenSubtractor::operator=(const AliEmcalJetUtilityGenSubtractor &other)
+{
+  // Assignment.
+
+  if (&other == this) return *this;
+  AliEmcalJetUtility::operator=(other);
+  fDoGenericSubtractionJetMass = other.fDoGenericSubtractionJetMass;
+  fDoGenericSubtractionGR = other.fDoGenericSubtractionGR;
+  fDoGenericSubtractionExtraJetShapes = other.fDoGenericSubtractionExtraJetShapes;
+  fUseExternalBkg = other.fUseExternalBkg;
+  fRhoName = other.fRhoName;
+  fRhomName = other.fRhomName;
+  fRho = other.fRho;
+  fRhom = other.fRhom;
+  fRMax = other.fRMax;
+  fDRStep = other.fDRStep;
+  fPtMinGR = other.fPtMinGR;
+  fRhoParam = other.fRhoParam;
+  fRhomParam = other.fRhomParam;
+  return *this;
+}
+
+//______________________________________________________________________________
+void AliEmcalJetUtilityGenSubtractor::Init()
+{
+  // Initialize the utility.
+
+  if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event
+    if(!fJetTask) return;
+    fRhoParam = dynamic_cast<AliRhoParameter*>(fJetTask->GetEvent()->FindListObject(fRhoName));
+    if (!fRhoParam) {
+      AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
+      return;
+    }
+  }
+  
+  if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event
+    fRhomParam = dynamic_cast<AliRhoParameter*>(fJetTask->GetEvent()->FindListObject(fRhomName));
+    if (!fRhomParam) {
+      AliError(Form("%s: Could not retrieve rho_m %s!", GetName(), fRhomName.Data()));
+      return;
+    }
+  }
+
+  fInit = kTRUE;
+}
+
+//______________________________________________________________________________
+void AliEmcalJetUtilityGenSubtractor::Prepare(AliFJWrapper& fjw)
+{
+  // Prepare the utility.
+
+  if (!fInit) return;
+  
+  if (fRhoParam) fRho = fRhoParam->GetVal();
+  if (fRhomParam) fRhom = fRhomParam->GetVal();
+
+  //run generic subtractor
+  if (fDoGenericSubtractionJetMass) {
+    fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+    fjw.DoGenericSubtractionJetMass();
+  }
+ if (fDoGenericSubtractionExtraJetShapes) {
+   fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
+   fjw.DoGenericSubtractionJetAngularity();
+   fjw.DoGenericSubtractionJetpTD();
+   fjw.DoGenericSubtractionJetCircularity();
+   fjw.DoGenericSubtractionJetSigma2();
+   fjw.DoGenericSubtractionJetConstituent();
+   fjw.DoGenericSubtractionJetLeSub();
+ }
+}
+
+//______________________________________________________________________________
+void AliEmcalJetUtilityGenSubtractor::ProcessJet(AliEmcalJet* jet, Int_t ij, AliFJWrapper& fjw)
+{
+  // Proceess each jet.
+
+  if (!fInit) return;
+
+#ifdef FASTJET_VERSION
+
+  if (fDoGenericSubtractionJetMass) {
+    std::vector<fastjet::contrib::GenericSubtractorInfo> jetMassInfo = fjw.GetGenSubtractorInfoJetMass();
+    Int_t n = (Int_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());
+    }
+  }
+
+  //here do generic subtraction for angular structure function
+  Double_t ptcorr = jet->Pt()-fjw.GetJetArea(ij)*fRho;
+  if (fDoGenericSubtractionGR && ptcorr>fPtMinGR) {
+    fjw.SetUseExternalBkg(fUseExternalBkg, fRho, fRhom);
+    fRMax = fJetTask->GetRadius()+0.2;
+    fjw.SetRMaxAndStep(fRMax, fDRStep);
+    fjw.DoGenericSubtractionGR(ij);
+    std::vector<double> num = fjw.GetGRNumerator();
+    std::vector<double> den = fjw.GetGRDenominator();
+    std::vector<double> nums = fjw.GetGRNumeratorSub();
+    std::vector<double> dens = fjw.GetGRDenominatorSub();
+    //pass this to AliEmcalJet
+    jet->SetGRNumSize(num.size());
+    jet->SetGRDenSize(den.size());
+    jet->SetGRNumSubSize(nums.size());
+    jet->SetGRDenSubSize(dens.size());
+    Int_t nsize = (Int_t)num.size();
+    for (Int_t g = 0; g < nsize; ++g) {
+      jet->AddGRNumAt(num[g],g);
+      jet->AddGRNumSubAt(nums[g],g);
+    }
+    Int_t dsize = (Int_t)den.size();
+    for (Int_t g = 0; g < dsize; ++g) {
+      jet->AddGRDenAt(den[g], g);
+      jet->AddGRDenSubAt(dens[g], g);
+    }
+  }
+
+  if (fDoGenericSubtractionExtraJetShapes) {
+    std::vector<fastjet::contrib::GenericSubtractorInfo> jetAngularityInfo = fjw.GetGenSubtractorInfoJetAngularity();
+    Int_t na = (Int_t)jetAngularityInfo.size();
+    if(na > ij && na > 0) {
+      jet->SetFirstDerivativeAngularity(jetAngularityInfo[ij].first_derivative());
+      jet->SetSecondDerivativeAngularity(jetAngularityInfo[ij].second_derivative());
+      jet->SetFirstOrderSubtractedAngularity(jetAngularityInfo[ij].first_order_subtracted());
+      jet->SetSecondOrderSubtractedAngularity(jetAngularityInfo[ij].second_order_subtracted());
+    }
+
+    std::vector<fastjet::contrib::GenericSubtractorInfo> jetpTDInfo = fjw.GetGenSubtractorInfoJetpTD();
+    Int_t np = (Int_t)jetpTDInfo.size();
+    if(np > ij && np > 0) {
+      jet->SetFirstDerivativepTD(jetpTDInfo[ij].first_derivative());
+      jet->SetSecondDerivativepTD(jetpTDInfo[ij].second_derivative());
+      jet->SetFirstOrderSubtractedpTD(jetpTDInfo[ij].first_order_subtracted());
+      jet->SetSecondOrderSubtractedpTD(jetpTDInfo[ij].second_order_subtracted());
+    }
+
+    std::vector<fastjet::contrib::GenericSubtractorInfo> jetCircularityInfo = fjw.GetGenSubtractorInfoJetCircularity();
+    Int_t nc = (Int_t)jetCircularityInfo.size();
+    if(nc > ij && nc > 0) {
+      jet->SetFirstDerivativeCircularity(jetCircularityInfo[ij].first_derivative());
+      jet->SetSecondDerivativeCircularity(jetCircularityInfo[ij].second_derivative());
+      jet->SetFirstOrderSubtractedCircularity(jetCircularityInfo[ij].first_order_subtracted());
+      jet->SetSecondOrderSubtractedCircularity(jetCircularityInfo[ij].second_order_subtracted());
+    }
+
+    std::vector<fastjet::contrib::GenericSubtractorInfo> jetSigma2Info = fjw.GetGenSubtractorInfoJetSigma2();
+    Int_t ns = (Int_t)jetSigma2Info.size();
+    if (ns > ij && ns > 0) {
+      jet->SetFirstDerivativeSigma2(jetSigma2Info[ij].first_derivative());
+      jet->SetSecondDerivativeSigma2(jetSigma2Info[ij].second_derivative());
+      jet->SetFirstOrderSubtractedSigma2(jetSigma2Info[ij].first_order_subtracted());
+      jet->SetSecondOrderSubtractedSigma2(jetSigma2Info[ij].second_order_subtracted());
+    }
+
+
+    std::vector<fastjet::contrib::GenericSubtractorInfo> jetConstituentInfo = fjw.GetGenSubtractorInfoJetConstituent();
+    Int_t nco = (Int_t)jetConstituentInfo.size();
+    if(nco > ij && nco > 0) {
+      jet->SetFirstDerivativeConstituent(jetConstituentInfo[ij].first_derivative());
+      jet->SetSecondDerivativeConstituent(jetConstituentInfo[ij].second_derivative());
+      jet->SetFirstOrderSubtractedConstituent(jetConstituentInfo[ij].first_order_subtracted());
+      jet->SetSecondOrderSubtractedConstituent(jetConstituentInfo[ij].second_order_subtracted());
+    }
+
+    std::vector<fastjet::contrib::GenericSubtractorInfo> jetLeSubInfo = fjw.GetGenSubtractorInfoJetLeSub();
+    Int_t nlsub = (Int_t)jetLeSubInfo.size();
+    if(nlsub > ij && nlsub > 0) {
+      jet->SetFirstDerivativeLeSub(jetLeSubInfo[ij].first_derivative());
+      jet->SetSecondDerivativeLeSub(jetLeSubInfo[ij].second_derivative());
+      jet->SetFirstOrderSubtractedLeSub(jetLeSubInfo[ij].first_order_subtracted());
+      jet->SetSecondOrderSubtractedLeSub(jetLeSubInfo[ij].second_order_subtracted());
+    }
+  }
+
+#endif
+}
+
+//______________________________________________________________________________
+void AliEmcalJetUtilityGenSubtractor::Terminate(AliFJWrapper& /*fjw*/)
+{
+  // Run termination of the utility (after each event).
+}
diff --git a/PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.h b/PWGJE/EMCALJetTasks/AliEmcalJetUtilityGenSubtractor.h
new file mode 100644 (file)
index 0000000..f66084b
--- /dev/null
@@ -0,0 +1,57 @@
+#ifndef ALIEMCALJETUTILITYGENSUBTRACTOR_H
+#define ALIEMCALJETUTILITYGENSUBTRACTOR_H
+
+#include <TNamed.h>
+
+#include "AliEmcalJetUtility.h"
+#include "AliFJWrapper.h"
+
+class AliEmcalJetTask;
+class AliEmcalJet;
+class AliFJWrapper;
+class AliRhoParameter;
+
+class AliEmcalJetUtilityGenSubtractor : public AliEmcalJetUtility
+{
+ public:
+
+  AliEmcalJetUtilityGenSubtractor();
+  AliEmcalJetUtilityGenSubtractor(const char* name);
+  AliEmcalJetUtilityGenSubtractor(const AliEmcalJetUtilityGenSubtractor &jet);
+  AliEmcalJetUtilityGenSubtractor& operator=(const AliEmcalJetUtilityGenSubtractor &jet);
+  ~AliEmcalJetUtilityGenSubtractor() {;}
+
+  void Init();
+  void Prepare(AliFJWrapper& fjw);
+  void ProcessJet(AliEmcalJet* jet, Int_t ij, AliFJWrapper& fjw);
+  void Terminate(AliFJWrapper& fjw);
+
+  void                   SetRhoName(const char *n)                                        { fRhoName      = n                      ; }
+  void                   SetRhomName(const char *n)                                       { fRhomName     = n                      ; }
+
+  void                   SetGenericSubtractionJetMass(Bool_t b)                           { fDoGenericSubtractionJetMass        = b; }
+  void                   SetGenericSubtractionGR(Bool_t b, Double_t rmax = 2., 
+                                                 Double_t dr = 0.04, Double_t ptmin = 0.) { fDoGenericSubtractionGR             = b; fRMax = rmax; fDRStep = dr; fPtMinGR = ptmin;}
+  void                   SetGenericSubtractionExtraJetShapes(Bool_t b)                    { fDoGenericSubtractionExtraJetShapes = b; }
+  void                   SetUseExternalBkg(Bool_t b)                                      { fUseExternalBkg                     = b; }
+
+ protected:
+
+  Bool_t                 fDoGenericSubtractionJetMass;        // calculate generic subtraction
+  Bool_t                 fDoGenericSubtractionGR;             // calculate generic subtraction for angular structure function GR
+  Bool_t                 fDoGenericSubtractionExtraJetShapes; // calculate generic subtraction for other jet shapes like radialmoment,pTD etc
+  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
+  Double_t               fRMax;                               // R max for GR calculation
+  Double_t               fDRStep;                             // step width for GR calculation
+  Double_t               fPtMinGR;                            // min pT for GR calculation
+
+  AliRhoParameter       *fRhoParam;                           //!event rho
+  AliRhoParameter       *fRhomParam;                          //!event rhom
+
+  ClassDef(AliEmcalJetUtilityGenSubtractor, 1) // Emcal jet utility that implements generic subtractors form the fastjet contrib
+};
+#endif
index 1d4fd71..c919da4 100644 (file)
@@ -1,8 +1,6 @@
 #ifndef AliFJWrapper_H
 #define AliFJWrapper_H
 
-// $Id$
-
 #if !defined(__CINT__) 
 
 #include <vector>
@@ -62,6 +60,8 @@ class AliFJWrapper
   virtual Int_t DoGenericSubtractionJetLeSub();
   virtual Int_t DoConstituentSubtraction();
 
+  void SetName(const char* name)        { fName           = name;    }
+  void SetTitle(const char* title)      { fTitle          = title;   }
   void SetStrategy(const fastjet::Strategy &strat)                 { fStrategy = strat;  }
   void SetAlgorithm(const fastjet::JetAlgorithm &algor)            { fAlgor    = algor;  }
   void SetRecombScheme(const fastjet::RecombinationScheme &scheme) { fScheme   = scheme; }
index 2c4dc3a..baf5143 100644 (file)
 #pragma link C++ class AliAnalysisTaskEmcalTriggerTreeWriter+;
 
 #ifdef HAVE_FASTJET
+#pragma link C++ class AliEmcalJetUtility+;
+#pragma link C++ class AliEmcalJetUtilityGenSubtractor+;
+#pragma link C++ class AliEmcalJetUtilityConstSubtractor+;
 #pragma link C++ class AliEmcalJetTask+;
 #pragma link C++ class AliEmcalJetFinder+;
 #pragma link C++ class AliJetEmbeddingFromAODTask+;