AliAnalysisTaskDmesonJetsSub.h

Go to the documentation of this file.

#ifndef ALIANALYSISTASKDMESONJETSSUB_H
#define ALIANALYSISTASKDMESONJETSSUB_H

/**************************************************************************
* Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

class TClonesArray;
class AliRDHFCuts;
class AliAODEvent;
class AliAODRecoDecay;
class AliAODRecoDecayHF2Prong;
class AliAODRecoCascadeHF;
class AliVParticle;
class AliAODMCParticle;
class AliHFAODMCParticleContainer;
class AliHFTrackContainer;
class AliParticleContainer;
class AliClusterContainer;
class THnSparse;
class AliFJWrapper;
class THashList;
class TTree;
class AliEMCALGeometry;
class TRandom;
class AliRhoParameter;

//C++
#include <exception>
#include <list>
#include <vector>
#include <map>

#include "AliTLorentzVector.h"
#include "THistManager.h"

#include "AliAnalysisTaskEmcalLight.h"
#include "AliJetContainer.h"

class AliAnalysisTaskDmesonJetsSub : public AliAnalysisTaskEmcalLight
{
 public:

  class AliEventNotFound : public std::exception
  {
  public:
    AliEventNotFound(const std::string& class_name, const std::string& method_name);
#if !(defined(__CINT__) || defined(__MAKECINT__))
    const char* what() const noexcept;
#endif


  private:
    std::string  fClassName        ; 
    std::string  fAccessMethodName ; 
    std::string  fWhat             ; 
  };

  typedef AliJetContainer::EJetType_t EJetType_t;
  typedef AliJetContainer::EJetAlgo_t EJetAlgo_t;
  typedef AliJetContainer::ERecoScheme_t ERecoScheme_t;

  enum EOutputType_t { kNoOutput, kTreeOutput, kTHnOutput, kOnlyQAOutput, kTreeExtendedOutput };
  enum ECandidateType_t  { kD0toKpi, kDstartoKpipi, kD0toKpiLikeSign };
  enum EMCMode_t { kNoMC, kSignalOnly, kBackgroundOnly, kMCTruth, kD0Reflection, kOnlyWrongPIDAccepted };
  enum EMesonOrigin_t {
    kUnknownQuark = BIT(0),
    kFromDown     = BIT(1),
    kFromUp       = BIT(2),
    kFromStrange  = BIT(3),
    kFromCharm    = BIT(4),
    kFromBottom   = BIT(5),
    kFromTop      = BIT(6),
    kFromGluon    = BIT(7),
    kAnyOrigin    = kUnknownQuark | kFromDown | kFromUp | kFromStrange | kFromCharm | kFromBottom | kFromTop | kFromGluon
  };

  enum EMesonDecayChannel_t {
    kAnyDecay            = 0,
    kUnknownDecay        = BIT(0),
    kDecayD0toKpi        = BIT(1),
    kDecayDStartoKpipi   = BIT(2)
  };

  enum EAxis_t {
   kPositionD         = BIT(0) , // Add the D meson eta/phi axis in the THnSparse
   kInvMass           = BIT(1) , // Add the invariant mass axis in the THnSparse
   k2ProngInvMass     = BIT(2) , // Add the 2 prong invariant mass axis in the THnSparse (for D* this is the inv mass of the D0)
   kSoftPionPt        = BIT(3) , // Add the soft pion pt axis in the THnSparse (for D*)
   kDeltaR            = BIT(4) , // Add the delta R axis in the THnSparse
   kDeltaEta          = BIT(5) , // Add the delta eta axis in the THnSparse
   kDeltaPhi          = BIT(6) , // Add the delta phi axis in the THnSparse
   kPositionJet       = BIT(7) , // Add the jet eta/phi axis in the THnSparse
   kJetConstituents   = BIT(8)   // Add the jet constituent axis in the THnSparse
  };

  class EventInfo {
  public:
    EventInfo() : fCent(-1), fEP(-1), fWeight(1), fPtHard(0) {}
    EventInfo(double cent, double ep, double w, double pt) : fCent(cent), fEP(ep), fWeight(w), fPtHard(pt) {}

    double fCent;
    double fEP;
    double fWeight;
    double fPtHard;
  };

  class AliEventInfoSummary {
  public:
    AliEventInfoSummary() : fWeight(1), fPtHard(0) {;}
    AliEventInfoSummary(EventInfo event);

    virtual ~AliEventInfoSummary() {;}

    void Reset();
    void Set(EventInfo event);

    Double32_t  fWeight          ; //[0,0,12]
    Double32_t  fPtHard          ; //[0,512,10]

    ClassDef(AliEventInfoSummary, 1);
  };

  class AliJetInfo {
  public:
    AliJetInfo() : fMomentum(), fNConstituents(0), fNEF(0), fMaxChargedPt(0), fMaxNeutralPt(0), fArea(0), fCorrPt(0) {}
    AliJetInfo(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t nconst, Double_t nef, Double_t cpt, Double_t npt) :
      fMomentum(px, py, pz, E), fNConstituents(nconst), fNEF(nef), fMaxChargedPt(cpt), fMaxNeutralPt(npt), fArea(0), fCorrPt(0) {}

    virtual ~AliJetInfo() {;}

    Double_t Pt()        const { return fMomentum.Pt()       ; }
    Double_t Eta()       const { return fMomentum.Eta()      ; }
    Double_t Phi()       const { return fMomentum.Phi()      ; }
    Double_t Phi_0_2pi() const { return fMomentum.Phi_0_2pi(); }
    Double_t CorrPt()    const { return fCorrPt              ; }
    Int_t GetNConstituents() const { return  fNConstituents; }
    Double_t GetDistance(const AliJetInfo& jet, Double_t& deta, Double_t& dphi) const;
    Double_t GetDistance(const AliJetInfo& jet) const;

    AliTLorentzVector fMomentum             ; 
    Int_t             fNConstituents        ; 
    Double_t          fNEF                  ; 
    Double_t          fMaxChargedPt         ; 
    Double_t          fMaxNeutralPt         ; 
    Double_t          fArea                 ; 
    Double_t          fCorrPt               ; 

    ClassDef(AliJetInfo, 2);
  };

  class AliDmesonJetInfo {
  public:
    AliDmesonJetInfo();
    AliDmesonJetInfo(const AliDmesonJetInfo &source);
    AliDmesonJetInfo& operator=(const AliDmesonJetInfo& source);
    virtual ~AliDmesonJetInfo() {;}

    AliVParticle      *fDmesonParticle          ; 
    AliTLorentzVector  fD                       ; 
    Double_t           fSoftPionPt              ; 
    Double_t           fInvMass2Prong           ; 
    std::map<std::string, AliJetInfo>
                       fJets                    ; 
    Int_t              fMCLabel                 ; 
    Bool_t             fReconstructed           ; 
    AliAODMCParticle  *fParton                  ; 
    Short_t            fPartonType              ; 
    AliAODMCParticle  *fAncestor                ; 
    Bool_t             fD0D0bar                 ; 
    Byte_t             fSelectionType           ; 
    const AliVEvent   *fEvent                   ; 

    Byte_t GetSelectionTypeSummary() const { return fD0D0bar && fSelectionType ? 3 : fSelectionType; }

    const AliJetInfo* GetJet(std::string n) const;
    AliJetInfo* GetJet(std::string n);
    void Reset();
    Double_t GetZ(std::string n) const;
    Double_t GetCorrZ(std::string n) const;
    Double_t GetDistance(std::string n, Double_t& deta, Double_t& dphi) const;
    Double_t GetDistance(std::string n) const;
    Double_t GetDistance(const AliJetInfo& jet, Double_t& deta, Double_t& dphi) const;
    Double_t GetDistance(const AliJetInfo& jet) const;
    void Print() const;

    ClassDef(AliDmesonJetInfo, 2);
  };

  class AliJetInfoSummary {
  public:
    AliJetInfoSummary() : fPt(0), fEta(0), fPhi(0), fR(0), fZ(0), fN(0) {;}
    AliJetInfoSummary(const AliDmesonJetInfo& source, std::string n);
    virtual ~AliJetInfoSummary() {}

    virtual void Reset();
    virtual void Set(const AliDmesonJetInfo& source, std::string n);
    virtual void Set(const AliJetInfo& source);

    Double32_t  fPt        ; //[0,409.6,13]
    Double32_t  fEta       ; //[-2.048,2.048,10]
    Double32_t  fPhi       ; //[0,2*pi,10]
    Double32_t  fR         ; //[0,2.56,7]
    Double32_t  fZ         ; //[0,1.024,10]
    Double32_t  fN         ; //[0, 64, 6]

    ClassDef(AliJetInfoSummary, 4);
  };

  class AliJetInfoPbPbSummary : public AliJetInfoSummary {
  public:
    AliJetInfoPbPbSummary() : AliJetInfoSummary(), fCorrPt(0), fCorrZ(0), fArea(0) {;}
    AliJetInfoPbPbSummary(const AliDmesonJetInfo& source, std::string n);
    virtual ~AliJetInfoPbPbSummary() {}

    virtual void Reset();
    virtual void Set(const AliDmesonJetInfo& source, std::string n);
    virtual void Set(const AliJetInfo& source);

    Double32_t  fCorrPt    ; //[-409.6,409.6,14]

    Double32_t  fCorrZ     ; //[0,1.024,10]

    Double32_t  fArea      ; //[0,2.048,8]

    ClassDef(AliJetInfoPbPbSummary, 2);
  };

  class AliDmesonInfoSummary {
  public:
    AliDmesonInfoSummary() : fPt(0), fEta(0), fPhi(0) {;}
    AliDmesonInfoSummary(const AliDmesonJetInfo& source);
    virtual ~AliDmesonInfoSummary() {}

    virtual void Reset();
    virtual void Set(const AliDmesonJetInfo& source);

    Double32_t   fPt     ; //[0,204.8,12]
    Double32_t   fEta    ; //[-2.048,2.048,10]
    Double32_t   fPhi    ; //[0,2*pi,10]

    ClassDef(AliDmesonInfoSummary, 2);
  };

  class AliDmesonMCInfoSummary : public AliDmesonInfoSummary {
  public:
    AliDmesonMCInfoSummary() : AliDmesonInfoSummary(), fPartonType(0), fPartonPt(0), fAncestorPDG(0) {;}
    AliDmesonMCInfoSummary(const AliDmesonJetInfo& source);
    virtual ~AliDmesonMCInfoSummary() {}

    virtual void Reset();
    virtual void Set(const AliDmesonJetInfo& source);

    Double32_t   fPartonType  ; //[0, 16, 4]
    Double32_t   fPartonPt    ; //[0,819.2,14]

    UShort_t     fAncestorPDG ; 

    ClassDef(AliDmesonMCInfoSummary, 3);
  };

  class AliD0InfoSummary : public AliDmesonInfoSummary {
  public:
    AliD0InfoSummary() : AliDmesonInfoSummary(), fInvMass(0), fSelectionType(0) {}
    AliD0InfoSummary(const AliDmesonJetInfo& source);
    virtual ~AliD0InfoSummary() {}

    virtual void Reset();
    virtual void Set(const AliDmesonJetInfo& source);

    Double32_t          fInvMass        ; //[0,6.5536,16]
    Double32_t          fSelectionType  ; //[0,4,2]

    ClassDef(AliD0InfoSummary, 3);
  };

  class AliD0ExtendedInfoSummary : public AliD0InfoSummary {
  public:
    AliD0ExtendedInfoSummary() : AliD0InfoSummary(), fDCA(0), fCosThetaStar(0), fPtK(0), fPtPi(0), fd0K(0), fd0Pi(0), fd0d0(0), fCosPointing(0), fMaxNormd0(0) {}
    AliD0ExtendedInfoSummary(const AliDmesonJetInfo& source);
    virtual ~AliD0ExtendedInfoSummary() {}

    virtual void Reset();
    virtual void Set(const AliDmesonJetInfo& source);

    Double32_t          fDCA           ; //[0,2.56,9]
    Double32_t          fCosThetaStar  ; //[-1.28,1.28,8]
    Double32_t          fPtK           ; //[0,102.4,11]
    Double32_t          fPtPi          ; //[0,102.4,11]
    Double32_t          fd0K           ; //[-2.56,2.56,8]
    Double32_t          fd0Pi          ; //[-2.56,2.56,8]
    Double32_t          fd0d0          ; //[-0.000524288,0.000524288,20]
    Double32_t          fCosPointing   ; //[0,1.28,7]
    Double32_t          fMaxNormd0     ; //[-102.4,102.4,11]

    ClassDef(AliD0ExtendedInfoSummary, 4);
  };

  class AliDStarInfoSummary : public AliDmesonInfoSummary {
  public:
    AliDStarInfoSummary() : AliDmesonInfoSummary(), f2ProngInvMass(0), fDeltaInvMass(0) {}
    AliDStarInfoSummary(const AliDmesonJetInfo& source);
    virtual ~AliDStarInfoSummary() {}

    virtual void Reset();
    virtual void Set(const AliDmesonJetInfo& source);

    Double32_t   f2ProngInvMass   ; //[0,8.192,14]
    Double32_t   fDeltaInvMass    ; //[0,0.8192,16]

    ClassDef(AliDStarInfoSummary, 2);
  };

  class AliHFJetDefinition : public TObject {
  public:
    AliHFJetDefinition();
    AliHFJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco);
    AliHFJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco, TString rhoName);
    AliHFJetDefinition(const AliHFJetDefinition &source);

    AliHFJetDefinition& operator=(const AliHFJetDefinition& source);

    const char* GetName() const;

    void SetJetPhiRange(Double_t min, Double_t max)       { fMinJetPhi    = min; fMaxJetPhi    = max; }
    void SetJetEtaRange(Double_t min, Double_t max)       { fMinJetEta    = min; fMaxJetEta    = max; }
    void SetJetPtRange(Double_t min, Double_t max)        { fMinJetPt     = min; fMaxJetPt     = max; }
    void SetChargedPtRange(Double_t min, Double_t max)    { fMinChargedPt = min; fMaxChargedPt = max; }
    void SetNeutralPtRange(Double_t min, Double_t max)    { fMinNeutralPt = min; fMaxNeutralPt = max; }
    void SetRhoName(TString n)                            { fRhoName      = n  ; }
    void SetRecombinationScheme(ERecoScheme_t reco)       { fRecoScheme   = reco;}
    Double_t GetRadius() const { return fRadius; }
    TString GetRhoName() const { return fRhoName; }

    Bool_t IsJetInAcceptance(const AliJetInfo& jet) const;
    Bool_t IsJetInAcceptance(const AliDmesonJetInfo& dMesonJet, std::string n) const;

    friend bool        operator< (const AliHFJetDefinition& lhs, const AliHFJetDefinition& rhs);
    friend inline bool operator> (const AliHFJetDefinition& lhs, const AliHFJetDefinition& rhs){ return rhs < lhs    ; }
    friend inline bool operator<=(const AliHFJetDefinition& lhs, const AliHFJetDefinition& rhs){ return !(lhs > rhs) ; }
    friend inline bool operator>=(const AliHFJetDefinition& lhs, const AliHFJetDefinition& rhs){ return !(lhs < rhs) ; }

    friend bool        operator==(const AliHFJetDefinition& lhs, const AliHFJetDefinition& rhs);
    friend inline bool operator!=(const AliHFJetDefinition& lhs, const AliHFJetDefinition& rhs){ return !(lhs == rhs); }

    EJetType_t                fJetType       ; 
    Double_t                  fRadius        ; 
    EJetAlgo_t                fJetAlgo       ; 
    ERecoScheme_t             fRecoScheme    ; 
    Double_t                  fMinJetPt      ; 
    Double_t                  fMaxJetPt      ; 
    Double_t                  fMinJetPhi     ; 
    Double_t                  fMaxJetPhi     ; 
    Double_t                  fMinJetEta     ; 
    Double_t                  fMaxJetEta     ; 
    Double_t                  fMinChargedPt  ; 
    Double_t                  fMaxChargedPt  ; 
    Double_t                  fMinNeutralPt  ; 
    Double_t                  fMaxNeutralPt  ; 
    TString                   fRhoName       ; 
    AliRhoParameter          *fRho           ; 
    std::vector<AliJetInfo>   fJets          ; 

  private:
    ClassDef(AliHFJetDefinition, 5);
  };

  class AnalysisEngine;

  class OutputHandler {
  public:
    OutputHandler();
    OutputHandler(AnalysisEngine* eng);

    virtual ~OutputHandler() {;}

    virtual Bool_t FillOutput(Bool_t applyKinCuts);

    virtual void BuildOutputObject(const char* taskName) {;}

    virtual TObject* GetOutputObject() const           { return 0 ; }
    virtual Int_t  GetDataSlotNumber() const           { return -1; }

  protected:
    ECandidateType_t                   fCandidateType         ; 
    EMCMode_t                          fMCMode                ; 
    Int_t                              fNMassBins             ; 
    Double_t                           fMinMass               ; 
    Double_t                           fMaxMass               ; 
    std::vector<AliHFJetDefinition>   *fJetDefinitions        ; 
    Float_t                            fPtBinWidth            ; 
    Float_t                            fMaxPt                 ; 
    Bool_t                             fD0Extended            ; 
    EventInfo                         *fEventInfo             ; 
    std::map<int, AliDmesonJetInfo>   *fDmesonJets            ; 
    THistManager                      *fHistManager           ; 
    TString                            fName                  ; 

  private:
    OutputHandler(const OutputHandler &source);  // not implemented
    OutputHandler& operator=(const OutputHandler& source);  // not implemented
  };

  class OutputHandlerTHnSparse : public OutputHandler {
  public:
    OutputHandlerTHnSparse();
    OutputHandlerTHnSparse(AnalysisEngine* eng);

    virtual ~OutputHandlerTHnSparse() {;}

    virtual Bool_t FillOutput(Bool_t applyKinCuts);

    virtual void BuildOutputObject(const char* taskName);

    Bool_t FillHnSparse(THnSparse* h, const AliDmesonJetInfo& DmesonJet, std::string n);

    Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo& dMesonJet) const;
    void SetEnabledAxis(UInt_t axis) { fEnabledAxis = axis; }

  protected:
    UInt_t               fEnabledAxis               ; 

  private:
    OutputHandlerTHnSparse(const OutputHandlerTHnSparse &source);  // not implemented
    OutputHandlerTHnSparse& operator=(const OutputHandlerTHnSparse& source);  // not implemented
  };

  class OutputHandlerTTree : public OutputHandler {
  public:
    OutputHandlerTTree();
    OutputHandlerTTree(AnalysisEngine* eng);

    virtual ~OutputHandlerTTree() {;}

    virtual Bool_t FillOutput(Bool_t applyKinCuts);

    virtual void BuildOutputObject(const char* taskName);
    virtual TObject* GetOutputObject() const           { return fTree; }
    virtual Int_t  GetDataSlotNumber() const           { return fDataSlotNumber    ; }

    TTree* GetTree() const { return fTree; }
    void   AssignDataSlot(Int_t n)             { fDataSlotNumber        = n; }

  protected:
    Int_t                              fDataSlotNumber        ; 
    TTree                             *fTree                  ; 
    AliDmesonInfoSummary              *fCurrentDmesonJetInfo  ; 
    AliJetInfoSummary                **fCurrentJetInfo        ; 

  private:
    OutputHandlerTTree(const OutputHandlerTTree &source);  // not implemented
    OutputHandlerTTree& operator=(const OutputHandlerTTree& source);  // not implemented
  };

  class OutputHandlerTTreeExtendedBase : public OutputHandler {
  public:
    OutputHandlerTTreeExtendedBase();
    OutputHandlerTTreeExtendedBase(AnalysisEngine* eng);

    virtual ~OutputHandlerTTreeExtendedBase() {;}

    virtual TObject* GetOutputObject() const           { return fTree; }
    virtual Int_t  GetDataSlotNumber() const           { return fDataSlotNumber    ; }

    TTree* GetTree() const { return fTree; }
    void   AssignDataSlot(Int_t n)             { fDataSlotNumber        = n; }

    static OutputHandlerTTreeExtendedBase* GenerateOutputHandler(AnalysisEngine* eng);

  protected:
    Int_t                                                       fDataSlotNumber        ; 
    TTree                                                      *fTree                  ; 
    TString                                                     fEventClassName        ;
    TString                                                     fDMesonClassName       ;
    TString                                                     fJetClassName          ;

  private:
    OutputHandlerTTreeExtendedBase(const OutputHandlerTTreeExtendedBase &source);  // not implemented
    OutputHandlerTTreeExtendedBase& operator=(const OutputHandlerTTreeExtendedBase& source);  // not implemented
  };

  template<class EVENTTYPE, class DMESONTYPE, class JETTYPE>
  class OutputHandlerTTreeExtended : public OutputHandlerTTreeExtendedBase {
  public:
    OutputHandlerTTreeExtended();
    OutputHandlerTTreeExtended(AnalysisEngine* eng);

    virtual ~OutputHandlerTTreeExtended() {;}

    virtual Bool_t FillOutput(Bool_t applyKinCuts);
    virtual void BuildOutputObject(const char* taskName);

  protected:
    EVENTTYPE                                                   fCurrentEventInfo      ; 
    std::vector<DMESONTYPE>                                     fCurrentDmesonInfo     ; 
    std::map<std::string, std::vector<JETTYPE> >                fCurrentJetInfo        ; 

  private:
    OutputHandlerTTreeExtended(const OutputHandlerTTreeExtended &source);  // not implemented
    OutputHandlerTTreeExtended& operator=(const OutputHandlerTTreeExtended& source);  // not implemented
  };

  class AnalysisEngine : public TObject {
  public:
    typedef std::pair<AliJetInfo*, Double_t> jet_distance_pair;

    enum EFindParticleOriginMode_t {
      kFindFirst,    
      kFindLast      
    };

    static AliAODMCParticle* FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode, const std::set<UInt_t>& pdgSet);
    static AliAODMCParticle* FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode);
    static std::pair<AliAnalysisTaskDmesonJetsSub::EMesonOrigin_t, AliAODMCParticle*> IsPromptCharm(const AliAODMCParticle* part, TClonesArray* mcArray);
    static EMesonDecayChannel_t CheckDecayChannel(const AliAODMCParticle* part, TClonesArray* mcArray);

    AnalysisEngine();
    AnalysisEngine(ECandidateType_t type, EMCMode_t MCmode, AliRDHFCuts* cuts = 0, Int_t nBins=80, Double_t range = 0.50);
    AnalysisEngine(const AnalysisEngine &source);
    AnalysisEngine& operator=(const AnalysisEngine& source);

    virtual ~AnalysisEngine();

    void SetCandidateType(ECandidateType_t t)       { fCandidateType  = t    ; }
    void SetMCMode(EMCMode_t m)                     { fMCMode         = m    ; }
    void SetNMassBins(Int_t n)                      { fNMassBins      = n    ; }
    void SetMassRange(Double_t min, Double_t max)   { fMinMass        = min  ; fMaxMass        = max  ; }
    void AdoptRDHFCuts(AliRDHFCuts* cuts);
    void SetRDHFCuts(AliRDHFCuts* cuts);
    void SetRejectedOriginMap(UInt_t m)             { fRejectedOrigin = m    ; }
    void SetAcceptedDecayMap(UInt_t m)              { fAcceptedDecay  = m    ; }
    void SetRejectISR(Bool_t b)                     { fRejectISR      = b    ; }

    void SetD0Extended(Bool_t b)                    { fD0Extended     = b    ; }

    const char* GetCandidateName() const { return fCandidateName.Data(); }
    const char* GetName() const;
    const char* GetName(const AliHFJetDefinition& jetDef) const;

    EMCMode_t GetMCMode()                     const { return fMCMode         ; }
    ECandidateType_t GetCandidateType()       const { return fCandidateType  ; }

    OutputHandler* GetOutputHandler() const { return fOutputHandler; }

    AliHFJetDefinition* AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco);
    AliHFJetDefinition* AddJetDefinition(const AliHFJetDefinition& def);
    std::vector<AliAnalysisTaskDmesonJetsSub::AliHFJetDefinition>::iterator FindJetDefinition(const AliHFJetDefinition& eng);
    std::vector<AliAnalysisTaskDmesonJetsSub::AliHFJetDefinition>& GetJetDefinitions() { return fJetDefinitions; }

    void SetJetPhiRange(Double_t min, Double_t max);
    void SetJetEtaRange(Double_t min, Double_t max);
    void SetJetPtRange(Double_t min, Double_t max);
    void SetChargedPtRange(Double_t min, Double_t max);
    void SetNeutralPtRange(Double_t min, Double_t max);

#if !(defined(__CINT__) || defined(__MAKECINT__))
    std::map<int, AliDmesonJetInfo>& GetDmesons() { return fDmesonJets; }
#endif

    void Init(const AliEMCALGeometry* const geom, Int_t runNumber);

    void   SetTrackEfficiency(Double_t t)      { fTrackEfficiency       = t; }

    Bool_t IsInhibit() const { return fInhibit; }

    Bool_t IsD0Extended() const { return fD0Extended; }

    friend bool        operator< (const AnalysisEngine& lhs, const AnalysisEngine& rhs);
    friend inline bool operator> (const AnalysisEngine& lhs, const AnalysisEngine& rhs){ return rhs < lhs    ; }
    friend inline bool operator<=(const AnalysisEngine& lhs, const AnalysisEngine& rhs){ return !(lhs > rhs) ; }
    friend inline bool operator>=(const AnalysisEngine& lhs, const AnalysisEngine& rhs){ return !(lhs < rhs) ; }

    friend bool        operator==(const AnalysisEngine& lhs, const AnalysisEngine& rhs);
    friend inline bool operator!=(const AnalysisEngine& lhs, const AnalysisEngine& rhs){ return !(lhs == rhs); }

  protected:
    void RunAnalysis();

    ECandidateType_t                   fCandidateType         ; 
    TString                            fCandidateName         ; 
    UInt_t                             fCandidatePDG          ; 
    UChar_t                            fNDaughters            ; 
    TArrayI                            fPDGdaughters          ; 
    TString                            fBranchName            ; 
    EMCMode_t                          fMCMode                ; 
    Int_t                              fNMassBins             ; 
    Double_t                           fMinMass               ; 
    Double_t                           fMaxMass               ; 
    AliRDHFCuts                       *fRDHFCuts              ; 
    UInt_t                             fRejectedOrigin        ; 
    UInt_t                             fAcceptedDecay         ; 
    Bool_t                             fInhibit               ; 
    std::vector<AliHFJetDefinition>    fJetDefinitions        ; 
    Float_t                            fPtBinWidth            ; 
    Float_t                            fMaxPt                 ; 
    Bool_t                             fD0Extended            ; 
    OutputHandler                     *fOutputHandler         ; 
    TRandom                           *fRandomGen             ; 
    Double_t                           fTrackEfficiency       ; 
    Bool_t                             fRejectISR             ; 
    std::map<int, AliDmesonJetInfo>    fDmesonJets            ; 
    TClonesArray                      *fCandidateArray        ; 
    AliHFAODMCParticleContainer*       fMCContainer           ; 
    std::vector<AliTrackContainer*>    fTrackContainers       ; 
    std::vector<AliClusterContainer*>  fClusterContainers     ; 
    AliAODEvent                       *fAodEvent              ; 
    AliFJWrapper                      *fFastJetWrapper        ; 
    THistManager                      *fHistManager           ; 
    EventInfo                          fEventInfo             ; 
    mutable TString                    fName                  ; 

    friend class AliAnalysisTaskDmesonJetsSub;
    friend class OutputHandler;

  private:

    void                AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist=0, Double_t eff=0.);
    void                SetCandidateProperties(Double_t range);
    AliAODMCParticle*   MatchToMC() const;
    void                RunDetectorLevelAnalysis();
    void                RunParticleLevelAnalysis();

    Bool_t              ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong* Dcand, AliDmesonJetInfo& DmesonJet, UInt_t i);
    Bool_t              ExtractD0Attributes(const AliAODRecoDecayHF2Prong* Dcand, AliDmesonJetInfo& DmesonJet, UInt_t i);
    Bool_t              ExtractDstarAttributes(const AliAODRecoCascadeHF* DstarCand, AliDmesonJetInfo& DmesonJet, UInt_t i);
    Bool_t              FindJet(AliAODRecoDecayHF2Prong* Dcand, AliDmesonJetInfo& DmesonJet, AliHFJetDefinition& jetDef, Int_t numcand);
    void                IterativeDeclustering(Int_t jetnum,Double_t type, AliHFJetDefinition& jetDef,Double_t invm); 
    ClassDef(AnalysisEngine, 3);
  };

  AliAnalysisTaskDmesonJetsSub();
  AliAnalysisTaskDmesonJetsSub(const char* name, Int_t nOutputTrees=2);
  virtual ~AliAnalysisTaskDmesonJetsSub();

  AnalysisEngine* AddAnalysisEngine(ECandidateType_t type, TString cutfname, TString cuttype, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString rhoName = "");
  AnalysisEngine* AddAnalysisEngine(ECandidateType_t type, TString cutfname, TString cuttype, EMCMode_t bkgMode, const AliHFJetDefinition& jetDef, TString rhoName = "");
  std::list<AliAnalysisTaskDmesonJetsSub::AnalysisEngine>::iterator FindAnalysisEngine(const AnalysisEngine& eng);

  void SetShowPositionD(Bool_t b = kTRUE)         { fEnabledAxis = b ?  fEnabledAxis | kPositionD         : fEnabledAxis & ~kPositionD         ; }
  void SetShowInvMass(Bool_t b = kTRUE)           { fEnabledAxis = b ?  fEnabledAxis | kInvMass           : fEnabledAxis & ~kInvMass           ; }
  void SetShow2ProngInvMass(Bool_t b = kTRUE)     { fEnabledAxis = b ?  fEnabledAxis | k2ProngInvMass     : fEnabledAxis & ~k2ProngInvMass     ; }
  void SetShowSoftPionPt(Bool_t b = kTRUE)        { fEnabledAxis = b ?  fEnabledAxis | kSoftPionPt        : fEnabledAxis & ~kSoftPionPt        ; }
  void SetShowDeltaR(Bool_t b = kTRUE)            { fEnabledAxis = b ?  fEnabledAxis | kDeltaR            : fEnabledAxis & ~kDeltaR            ; }
  void SetShowDeltaEta(Bool_t b = kTRUE)          { fEnabledAxis = b ?  fEnabledAxis | kDeltaEta          : fEnabledAxis & ~kDeltaEta          ; }
  void SetShowDeltaPhi(Bool_t b = kTRUE)          { fEnabledAxis = b ?  fEnabledAxis | kDeltaPhi          : fEnabledAxis & ~kDeltaPhi          ; }
  void SetShowPositionJet(Bool_t b = kTRUE)       { fEnabledAxis = b ?  fEnabledAxis | kPositionJet       : fEnabledAxis & ~kPositionJet       ; }
  void SetShowJetConstituents(Bool_t b = kTRUE)   { fEnabledAxis = b ?  fEnabledAxis | kJetConstituents   : fEnabledAxis & ~kJetConstituents   ; }

  UInt_t GetEnabledAxis() const { return fEnabledAxis; }

  void SetApplyKinematicCuts(Bool_t b)            { fApplyKinematicCuts = b ; }
  void SetOutputType(EOutputType_t b)             { SetOutputTypeInternal(b); }
  void SetTrackEfficiency(Double_t t)             { fTrackEfficiency    = t ; }
  void SetRejectISR(Bool_t b)                     { fRejectISR          = b ; }
  void SetJetArea(Int_t type,
      Double_t garea = 0.005)                     { fJetAreaType        = type; fJetGhostArea = garea; }

  virtual void         UserCreateOutputObjects();
  virtual void         ExecOnce();
  virtual Bool_t       Run();
  virtual Bool_t       FillHistograms();

  static AliAnalysisTaskDmesonJetsSub* AddTaskDmesonJetsSub(TString ntracks = "usedefault", TString nclusters = "usedefault", TString nMCpart = "", Int_t nMaxTrees = 2, TString suffix = "");

 protected:

  virtual void SetOutputTypeInternal(EOutputType_t b)             { fOutputType         = b; }

  AliRDHFCuts*         LoadDMesonCutsFromFile(TString cutfname, TString cutsname);

  static const char*   GetHFEventRejectionReasonLabel(UInt_t& bitmap);
  static void          CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t& minMass, Double_t& maxMass);

  Int_t                PostDataFromAnalysisEngine(OutputHandler const* handler);

  void                 FillPartonLevelHistograms();

  std::list<AnalysisEngine>
                       fAnalysisEngines           ; 
  UInt_t               fEnabledAxis               ; 
  EOutputType_t        fOutputType                ; 
  THistManager         fHistManager               ; 
  Bool_t               fApplyKinematicCuts        ; 
  Int_t                fNOutputTrees              ; 
  Double_t             fTrackEfficiency           ; 
  Bool_t               fRejectISR                 ; 
  Int_t                fJetAreaType               ; 
  Double_t             fJetGhostArea              ; 
  AliHFAODMCParticleContainer* fMCContainer       ; 
  AliAODEvent         *fAodEvent                  ; 
  AliFJWrapper        *fFastJetWrapper            ; 

 private:

  AliAnalysisTaskDmesonJetsSub(const AliAnalysisTaskDmesonJetsSub &source);
  AliAnalysisTaskDmesonJetsSub& operator=(const AliAnalysisTaskDmesonJetsSub& source);

  ClassDef(AliAnalysisTaskDmesonJetsSub, 10);
};

#endif