AliAnalysisTaskDmesonJets.h

Go to the documentation of this file.

#ifndef ALIANALYSISTASKDMESONJETS_H
#define ALIANALYSISTASKDMESONJETS_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;

#include <list>
#include <vector>
#include <map>

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

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

class AliAnalysisTaskDmesonJets : public AliAnalysisTaskEmcalLight
{
 public:

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

  enum EOutputType_t { kNoOutput, kTreeOutput, kTHnOutput };
  enum ECandidateType_t  { kD0toKpi, kDstartoKpipi, kD0toKpiLikeSign };
  enum EMCMode_t { kNoMC, kSignalOnly, kBackgroundOnly, kMCTruth };
  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)
  };

  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 AliJetInfo {
  public:
    AliJetInfo() : fMomentum(), fNConstituents(0), fNEF(0), fMaxChargedPt(0), fMaxNeutralPt(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) {}

    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 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         ; 

    ClassDef(AliJetInfo, 1);
  };

  class AliDmesonJetInfo {
  public:
    AliDmesonJetInfo() : fDmesonParticle(0), fD(), fSoftPionPt(0), fInvMass2Prong(0), fJets(), fMCLabel(-1), fReconstructed(kFALSE) {}

    virtual ~AliDmesonJetInfo() {;}

    AliVParticle      *fDmesonParticle          ; 
    AliTLorentzVector  fD                       ; 
    Double_t           fSoftPionPt              ; 
    Double_t           fInvMass2Prong           ; 
    std::map<std::string, AliJetInfo>
                       fJets                    ; 
    Int_t              fMCLabel                 ; 
    Bool_t             fReconstructed           ; 

    const AliJetInfo* GetJet(std::string n) const;
    AliJetInfo* GetJet(std::string n);
    void Reset();
    Double_t GetZ(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, 1);
  };

  class AliJetInfoSummary {
  public:
    AliJetInfoSummary() : fPt(0), fEta(0), fPhi(0), fR(0), fZ(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]

    ClassDef(AliJetInfoSummary, 3);
  };

  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 AliD0InfoSummary : public AliDmesonInfoSummary {
  public:
    AliD0InfoSummary() : AliDmesonInfoSummary(), fInvMass(0) {}
    AliD0InfoSummary(const AliDmesonJetInfo& source);
    virtual ~AliD0InfoSummary() {}

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

    Double32_t   fInvMass   ; //[0,6.5536,16]

    ClassDef(AliD0InfoSummary, 2);
  };

  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(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; }
    Double_t GetRadius() const { return fRadius; }

    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); }

  protected:
    friend class AliAnalysisTaskDmesonJets;
    friend class AnalysisEngine;

    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  ; 
    std::vector<AliJetInfo>   fJets          ; 

  private:
    ClassDef(AliHFJetDefinition, 4);
  };

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

    static EMesonOrigin_t CheckOrigin(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    ; }

    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  ; }

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

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

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

    TTree* BuildTree(const char* taskName);
    TTree* GetTree() const { return fTree; }
    Bool_t FillTree(Bool_t applyKinCuts);

    void AssignDataSlot(Int_t n) { fDataSlotNumber = n; }
    Int_t GetDataSlotNumber() const { return fDataSlotNumber; }

    void BuildHnSparse(UInt_t enabledAxis);
    Bool_t FillHnSparse(Bool_t applyKinCuts);
    Bool_t FillHnSparse(THnSparse* h, const AliDmesonJetInfo& DmesonJet, std::string n);

    Bool_t FillQA(Bool_t applyKinCuts);

    Bool_t IsInhibit() const { return fInhibit; }

    jet_distance_pair FindJetMacthedToGeneratedDMeson(const AliDmesonJetInfo& dmeson, AliHFJetDefinition& jetDef, Double_t dMax, Bool_t applyKinCuts);

    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();
    void FindJets(AliHFJetDefinition& jetDef);

    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                 ; 
    Int_t                              fDataSlotNumber        ; 
    TTree                             *fTree                  ; 
    AliDmesonInfoSummary              *fCurrentDmesonJetInfo  ; 
    AliJetInfoSummary                **fCurrentJetInfo        ; 
    std::map<int, AliDmesonJetInfo>    fDmesonJets            ; 
    TClonesArray                      *fCandidateArray        ; 
    AliHFAODMCParticleContainer       *fMCContainer           ; 
    AliHFTrackContainer               *fTrackContainer        ; 
    AliClusterContainer               *fClusterContainer      ; 
    AliAODEvent                       *fAodEvent              ; 
    AliFJWrapper                      *fFastJetWrapper        ; 
    THistManager                      *fHistManager           ; 

    friend class AliAnalysisTaskDmesonJets;

  private:

    void                AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist=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);

    ClassDef(AnalysisEngine, 2);
  };

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

  AnalysisEngine* AddAnalysisEngine(ECandidateType_t type, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString cutfname = "");
  AnalysisEngine* AddAnalysisEngine(ECandidateType_t type, EMCMode_t bkgMode, const AliHFJetDefinition& jetDef, TString cutfname = "");
  std::list<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   ; }

  void SetApplyKinematicCuts(Bool_t b)            { fApplyKinematicCuts = b ; }
  void SetOutputType(EOutputType_t b)             { SetOutputTypeInternal(b); }

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

 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(const AnalysisEngine& eng);

  std::list<AnalysisEngine>
                       fAnalysisEngines           ; 
  UInt_t               fEnabledAxis               ; 
  EOutputType_t        fOutputType                ; 
  THistManager         fHistManager               ; 
  Bool_t               fApplyKinematicCuts        ; 
  Int_t                fNOutputTrees              ; 
  AliAODEvent         *fAodEvent                  ; 
  AliFJWrapper        *fFastJetWrapper            ; 

 private:
   
  AliAnalysisTaskDmesonJets(const AliAnalysisTaskDmesonJets &source);
  AliAnalysisTaskDmesonJets& operator=(const AliAnalysisTaskDmesonJets& source);

  ClassDef(AliAnalysisTaskDmesonJets, 6);
};

#endif