AliAnalysisTaskJetExtractor.h

Go to the documentation of this file.

#ifndef ALIANALYSISTASKJETEXTRACTOR_H
#define ALIANALYSISTASKJETEXTRACTOR_H

/* Copyright(c) 1998-2018, ALICE Experiment at CERN, All rights reserved. *
 * See cxx source for full Copyright notice                               */


//###############################################################################################################################################3
class AliRDHFJetsCutsVertex;
class AliEmcalJetTree;
class AliHFJetsTaggingVertex;

// 

class AliAnalysisTaskJetExtractor : public AliAnalysisTaskEmcalJet {
 public:
  // ## Define class to store and compare secondary vertices (needed for matching of jet splittings & sec. vertices)
  struct SimpleSecondaryVertex
  {
    Int_t           fIndex;         // Index of sec. vertices in the array that we save to the tree
    Double_t        fLxy;           // decay length
    AliVParticle*   fDaughter1;     // daughter particle
    AliVParticle*   fDaughter2;     // daughter particle
    AliVParticle*   fDaughter3;     // daughter particle
  };

  AliAnalysisTaskJetExtractor();
  AliAnalysisTaskJetExtractor(const char *name);
  virtual ~AliAnalysisTaskJetExtractor();
  static AliAnalysisTaskJetExtractor* AddTaskJetExtractor(TString trackArray, TString clusterArray, TString jetArray, TString rhoObject, Double_t jetRadius, AliRDHFJetsCutsVertex* vertexerCuts, const char* taskNameSuffix);
  void                        UserCreateOutputObjects();
  void                        Terminate(Option_t *option);
  AliEmcalJetTree*            GetJetTree() {return fJetTree;}

  void                        SetSaveConstituents(Bool_t val) {fSaveConstituents = val; fInitialized = kFALSE;}
  void                        SetSaveConstituentsIP(Bool_t val) {fSaveConstituentsIP = val; fInitialized = kFALSE;}
  void                        SetSaveConstituentPID(Bool_t val) {fSaveConstituentPID = val; fInitialized = kFALSE;}
  void                        SetSaveJetShapes(Bool_t val) {fSaveJetShapes = val; fInitialized = kFALSE;}
  void                        SetSaveJetSplittings(Bool_t val) {fSaveJetSplittings = val; fInitialized = kFALSE;}
  void                        SetSaveMCInformation(Bool_t val) {fSaveMCInformation = val; fInitialized = kFALSE;}
  void                        SetSaveSecondaryVertices(Bool_t val) {fSaveSecondaryVertices = val; fInitialized = kFALSE;}
  void                        SetSaveTriggerTracks(Bool_t val) {fSaveTriggerTracks = val; fInitialized = kFALSE;}
  void                        SetSaveCaloClusters(Bool_t val) {fSaveCaloClusters = val; fInitialized = kFALSE;}

  void                        ActivateJetMatching(const char* arrayName, const char* rhoName = 0, const char* rhoMassName = 0)
                                {fMatchedJetsArrayName = arrayName; fMatchedJetsRhoName = rhoName ? rhoName : ""; fMatchedJetsRhoMassName = rhoMassName ? rhoMassName : "";}
  void                        SetMCParticleArrayName(const char* name)            { fMCParticleArrayName = name; }
  void                        SetHadronMatchingRadius(Double_t val)               { fHadronMatchingRadius  = val; }
  void                        SetSecondaryVertexMaxChi2(Double_t val   )          { fSecondaryVertexMaxChi2 = val; }
  void                        SetSecondaryVertexMaxDispersion(Double_t val)       { fSecondaryVertexMaxDispersion = val; }
  void                        SetCustomStartupScript(const char* path)            { fCustomStartupScript = path; }
  void                        SetVertexerCuts(AliRDHFJetsCutsVertex* val)         { fVertexerCuts = val; }
  void                        SetSetEmcalJetFlavour(Bool_t val)                   { fSetEmcalJetFlavour = val; }
  void                        SetEventPercentage(Double_t val)                    { fEventPercentage  = val; }
  void                        SetTruthLabelRange(Int_t min, Int_t max)            { fTruthMinLabel = min; fTruthMaxLabel = max; }
  void                        SetSaveTrackPDGCode(Bool_t val)                     { fSaveTrackPDGCode = val; }
  void                        SetRandomSeed(ULong_t val)                          { fRandomSeed  = val; }
  void                        SetRandomSeedCones(ULong_t val)                     { fRandomSeedCones  = val; }

  void                        SetEventCutTriggerTrack(Double_t minPt, Double_t maxPt, Int_t minLabel=-9999999, Int_t maxLabel=+9999999)
                                { fEventCut_TriggerTrackMinPt = minPt; fEventCut_TriggerTrackMaxPt = maxPt; fEventCut_TriggerTrackMinLabel = minLabel;
                                  fEventCut_TriggerTrackMaxLabel = maxLabel;}

 protected:
  Bool_t                      CreateControlHistograms();
  void                        ExecOnce();
  Bool_t                      Run();
  void                        FillTrackControlHistograms(AliVTrack* track);
  void                        FillEventControlHistograms();
  void                        FillJetControlHistograms(AliEmcalJet* jet);
  void                        CalculateJetShapes(AliEmcalJet* jet, Double_t& leSub_noCorr, Double_t& angularity, Double_t& momentumDispersion, Double_t& trackPtMean, Double_t& trackPtMedian);
  void                        GetTrueJetPtFraction(AliEmcalJet* jet, Double_t& truePtFraction, Double_t& truePtFraction_mcparticles);
  void                        GetMatchedJetObservables(AliEmcalJet* jet, Double_t& matchedJetPt, Double_t& matchedJetMass, Double_t& matchedJetDistance);
  void                        GetJetType(AliEmcalJet* jet, Int_t& typeHM, Int_t& typePM, Int_t& typeIC);
  Bool_t                      IsTriggerTrackInEvent();
  Bool_t                      IsTrackInCone(const AliVParticle* track, Double_t eta, Double_t phi, Double_t radius);
  Bool_t                      IsStrangeJet(AliEmcalJet* jet);
  void                        PrintConfig();
  void                        AddPIDInformation(const AliVParticle* particle, Float_t& sigITS, Float_t& sigTPC, Float_t& sigTOF, Float_t& sigTRD, Short_t& recoPID, Int_t& truePID);
  void                        GetTrackImpactParameters(const AliVVertex* vtx, const AliAODTrack* track, Float_t& d0, Float_t& d0cov, Float_t& z0, Float_t& z0cov);
  void                        ReconstructSecondaryVertices(const AliVVertex* primVtx, const AliEmcalJet* jet, std::vector<Float_t>& secVtx_X, std::vector<Float_t>& secVtx_Y, std::vector<Float_t>& secVtx_Z, std::vector<Float_t>& secVtx_Mass, std::vector<Float_t>& secVtx_Lxy, std::vector<Float_t>& secVtx_SigmaLxy, std::vector<Float_t>& secVtx_Chi2, std::vector<Float_t>& secVtx_Dispersion);
  void                        GetJetSplittings(AliEmcalJet* jet, std::vector<Float_t>& splittings_radiatorE, std::vector<Float_t>& splittings_kT, std::vector<Float_t>& splittings_theta, std::vector<Int_t>& splittings_secVtx_rank, std::vector<Int_t>& splittings_secVtx_index);

  // ################## CUTS AND SETTINGS: Should be set during task initialization
  Bool_t                      fSaveConstituents;                        
  Bool_t                      fSaveConstituentsIP;                      
  Bool_t                      fSaveConstituentPID;                      
  Bool_t                      fSaveJetShapes;                           
  Bool_t                      fSaveJetSplittings;                       
  Bool_t                      fSaveMCInformation;                       
  Bool_t                      fSaveSecondaryVertices;                   
  Bool_t                      fSaveTriggerTracks;                       
  Bool_t                      fSaveCaloClusters;                        

  Double_t                    fEventPercentage;                         
  Double_t                    fEventCut_TriggerTrackMinPt;              
  Double_t                    fEventCut_TriggerTrackMaxPt;              
  Int_t                       fEventCut_TriggerTrackMinLabel;           
  Int_t                       fEventCut_TriggerTrackMaxLabel;           
  Int_t                       fEventCut_TriggerTrackOrigin;             
  Int_t                       fTruthMinLabel;                           
  Int_t                       fTruthMaxLabel;                           
  Double_t                    fHadronMatchingRadius;                    

  TString                     fMatchedJetsArrayName;                    
  TString                     fMatchedJetsRhoName;                      
  TString                     fMatchedJetsRhoMassName;                  
  TString                     fMCParticleArrayName;                     

  ULong_t                     fRandomSeed;                              
  ULong_t                     fRandomSeedCones;                         
  AliRDHFJetsCutsVertex*      fVertexerCuts;                            

  Double_t                    fSecondaryVertexMaxChi2;                  
  Double_t                    fSecondaryVertexMaxDispersion;            
  TString                     fCustomStartupScript;                     
  Bool_t                      fSetEmcalJetFlavour;                      
  Bool_t                      fSaveTrackPDGCode;                        

  // ################## BUFFERS: Variables set during execution time
  AliEmcalJetTree*            fJetTree;
  Double_t                    fEventWeight;                             
  Double_t                    fImpactParameter;                         
  Int_t                       fMultiplicity;                            
  std::vector<Float_t>        fTriggerTracks_Pt;                        
  std::vector<Float_t>        fTriggerTracks_Eta;                       
  std::vector<Float_t>        fTriggerTracks_Phi;                       
  TClonesArray*               fMCParticleArray;                         
  TRandom3*                   fRandomGenerator;                         
  TRandom3*                   fRandomGeneratorCones;                    
  AliHFJetsTaggingVertex*     fVtxTagger;                               
  Bool_t                      fIsEmbeddedEvent;                         
  std::vector<SimpleSecondaryVertex> fSimpleSecVertices;  

  // ################## HELPER FUNCTIONS
  Double_t                    GetDistance(Double_t eta1, Double_t eta2, Double_t phi1, Double_t phi2)
  {
    Double_t deltaPhi = TMath::Min(TMath::Abs(phi1-phi2),TMath::TwoPi() - TMath::Abs(phi1-phi2));
    return TMath::Sqrt((eta1-eta2)*(eta1-eta2) + deltaPhi*deltaPhi);
  }
  void                        FillHistogram(const char * key, Double_t x);
  void                        FillHistogram(const char * key, Double_t x, Double_t y);
  void                        FillHistogram(const char * key, Double_t x, Double_t y, Double_t add);
  void                        FillHistogram3D(const char * key, Double_t x, Double_t y, Double_t z, Double_t add = 0);
  template <class T> T*       AddHistogram1D(const char* name = "CustomHistogram", const char* title = "NO_TITLE", const char* options = "", Int_t xBins = 100, Double_t xMin = 0.0, Double_t xMax = 20.0, const char* xTitle = "x axis", const char* yTitle = "y axis");
  template <class T> T*       AddHistogram2D(const char* name = "CustomHistogram", const char* title = "NO_TITLE", const char* options = "", Int_t xBins = 100, Double_t xMin = 0.0, Double_t xMax = 20.0, Int_t yBins = 100, Double_t yMin = 0.0, Double_t yMax = 20.0,  const char* xTitle = "x axis", const char* yTitle = "y axis", const char* zTitle = "z axis");
  template <class T> T*       AddHistogram3D(const char* name = "CustomHistogram", const char* title = "NO_TITLE", const char* options = "", Int_t xBins = 100, Double_t xMin = 0.0, Double_t xMax = 20.0, Int_t yBins = 100, Double_t yMin = 0.0, Double_t yMax = 20.0, Int_t zBins = 100, Double_t zMin = 0.0, Double_t zMax = 20.0, const char* xTitle = "x axis", const char* yTitle = "y axis", const char* zTitle = "z axis");

 private:
  AliAnalysisTaskJetExtractor(const AliAnalysisTaskJetExtractor&);            // not implemented
  AliAnalysisTaskJetExtractor &operator=(const AliAnalysisTaskJetExtractor&); // not implemented

  ClassDef(AliAnalysisTaskJetExtractor, 2) // Jet extraction task
};


//###############################################################################################################################################3
const Int_t kMaxNumConstituents = 500;
// 
class AliEmcalJetTree : public TNamed
{
  public:
    AliEmcalJetTree();
    AliEmcalJetTree(const char* name);

    // ############ SETTERS
    void            AddExtractionPercentage(Float_t minPt, Float_t maxPt, Float_t percentage)
    {
      fExtractionPercentagePtBins.push_back(minPt);
      fExtractionPercentagePtBins.push_back(maxPt);
      fExtractionPercentages.push_back(percentage);
    }
    void            ResetExtractionPercentages() {fExtractionPercentages.clear(); fExtractionPercentagePtBins.clear();}
    void            AddExtractionJetTypeHM(Int_t type) {fExtractionJetTypes_HM.push_back(type);}
    void            AddExtractionJetTypePM(Int_t type) {fExtractionJetTypes_PM.push_back(type);}

    void            InitializeTree(Bool_t saveCaloClusters, Bool_t saveMCInformation, Bool_t saveConstituents, Bool_t saveConstituentsIP, Bool_t saveConstituentPID, Bool_t saveJetShapes, Bool_t saveSplittings, Bool_t saveSecondaryVertices, Bool_t saveTriggerTracks);

    // ######################################
    Bool_t          AddJetToTree(AliEmcalJet* jet, Bool_t saveConstituents, Bool_t saveConstituentsIP, Bool_t saveCaloClusters, Double_t* vertex, Float_t rho, Float_t rhoMass, Float_t centrality, Int_t multiplicity, Long64_t eventID, Float_t magField);
    void            FillBuffer_SecVertices(std::vector<Float_t>& secVtx_X, std::vector<Float_t>& secVtx_Y, std::vector<Float_t>& secVtx_Z, std::vector<Float_t>& secVtx_Mass, std::vector<Float_t>& secVtx_Lxy, std::vector<Float_t>& secVtx_SigmaLxy, std::vector<Float_t>& secVtx_Chi2, std::vector<Float_t>& secVtx_Dispersion);
    void            FillBuffer_JetShapes(AliEmcalJet* jet, Double_t leSub_noCorr, Double_t angularity, Double_t momentumDispersion, Double_t trackPtMean, Double_t trackPtMedian);
    void            FillBuffer_Splittings(std::vector<Float_t>& splittings_radiatorE, std::vector<Float_t>& splittings_kT, std::vector<Float_t>& splittings_theta, Bool_t saveSecondaryVertices, std::vector<Int_t>& splittings_secVtx_rank, std::vector<Int_t>& splittings_secVtx_index);
    void            FillBuffer_PID(std::vector<Float_t>& trackPID_ITS, std::vector<Float_t>& trackPID_TPC, std::vector<Float_t>& trackPID_TOF, std::vector<Float_t>& trackPID_TRD, std::vector<Short_t>& trackPID_Reco, std::vector<Int_t>& trackPID_Truth);
    void            FillBuffer_MonteCarlo(Int_t motherParton, Int_t motherHadron, Int_t partonInitialCollision, Float_t matchJetDistance, Float_t matchedJetPt, Float_t matchedJetMass, Float_t truePtFraction, Float_t truePtFraction_mcparticles, Float_t ptHard, Float_t eventWeight, Float_t impactParameter);
    void            FillBuffer_ImpactParameters(std::vector<Float_t>& trackIP_d0, std::vector<Float_t>& trackIP_z0, std::vector<Float_t>& trackIP_d0cov, std::vector<Float_t>& trackIP_z0cov);
    void            FillBuffer_TriggerTracks(std::vector<Float_t>& triggerTrackPt, std::vector<Float_t>& triggerTrackDeltaEta, std::vector<Float_t>& triggerTrackDeltaPhi);
    // ######################################

    void            SetRandomGenerator(TRandom3* gen) {fRandomGenerator = gen;}
    std::vector<Float_t> GetExtractionPercentagePtBins() {return fExtractionPercentagePtBins;}
    std::vector<Float_t> GetExtractionPercentages() {return fExtractionPercentages;}
    std::vector<Int_t> GetExtractionJetTypes_HM() {return fExtractionJetTypes_HM;}
    std::vector<Int_t> GetExtractionJetTypes_PM() {return fExtractionJetTypes_PM;}
    TTree*          GetTreePointer() {return fJetTree;}

  private:
    TTree*          fJetTree;                             
    Bool_t          fInitialized;                         
    TRandom3*       fRandomGenerator;                     

    // Option flags
    std::vector<Float_t> fExtractionPercentages;          
    std::vector<Float_t> fExtractionPercentagePtBins;     
    std::vector<Int_t>   fExtractionJetTypes_HM;          
    std::vector<Int_t>   fExtractionJetTypes_PM;          

    // Buffers that will be added to the tree
    Float_t         fBuffer_JetPt;                        
    Float_t         fBuffer_JetEta;                       
    Float_t         fBuffer_JetPhi;                       
    Float_t         fBuffer_JetArea;                      
    Int_t           fBuffer_NumTracks;                    
    Int_t           fBuffer_NumClusters;                  

    Float_t         fBuffer_Event_BackgroundDensityMass;  
    Float_t         fBuffer_Event_BackgroundDensity;      
    Float_t         fBuffer_Event_Vertex_X;               
    Float_t         fBuffer_Event_Vertex_Y;               
    Float_t         fBuffer_Event_Vertex_Z;               
    Float_t         fBuffer_Event_Centrality;             
    Int_t           fBuffer_Event_Multiplicity;           
    Long64_t        fBuffer_Event_ID;                     
    Float_t         fBuffer_Event_MagneticField;          
    Float_t         fBuffer_Event_PtHard;                 
    Float_t         fBuffer_Event_Weight;                 
    Float_t         fBuffer_Event_ImpactParameter;        

    Float_t*        fBuffer_Track_Pt;                     
    Float_t*        fBuffer_Track_Eta;                    
    Float_t*        fBuffer_Track_Phi;                    
    Float_t*        fBuffer_Track_Charge;                 
    Int_t*          fBuffer_Track_Label;                  
    Float_t*        fBuffer_Track_ProdVtx_X;              
    Float_t*        fBuffer_Track_ProdVtx_Y;              
    Float_t*        fBuffer_Track_ProdVtx_Z;              

    Float_t*        fBuffer_Cluster_Pt;                   
    Float_t*        fBuffer_Cluster_E;                    
    Float_t*        fBuffer_Cluster_Eta;                  
    Float_t*        fBuffer_Cluster_Phi;                  
    Float_t*        fBuffer_Cluster_M02;                  
    Float_t*        fBuffer_Cluster_Time;                 
    Int_t*          fBuffer_Cluster_Label;                

    Float_t         fBuffer_Shape_Mass_NoCorr;            
    Float_t         fBuffer_Shape_Mass_DerivCorr_1;       
    Float_t         fBuffer_Shape_Mass_DerivCorr_2;       
    Float_t         fBuffer_Shape_pTD_DerivCorr_1;        
    Float_t         fBuffer_Shape_pTD_DerivCorr_2;        
    Float_t         fBuffer_Shape_LeSub_NoCorr;           
    Float_t         fBuffer_Shape_LeSub_DerivCorr;        
    Float_t         fBuffer_Shape_Angularity_NoCorr;      
    Float_t         fBuffer_Shape_Angularity_DerivCorr_1; 
    Float_t         fBuffer_Shape_Angularity_DerivCorr_2; 
    Float_t         fBuffer_Shape_Circularity_DerivCorr_1;
    Float_t         fBuffer_Shape_Circularity_DerivCorr_2;
    Float_t         fBuffer_Shape_Sigma2_DerivCorr_1;     
    Float_t         fBuffer_Shape_Sigma2_DerivCorr_2;     
    Float_t         fBuffer_Shape_NumTracks_DerivCorr;    
    Float_t         fBuffer_Shape_MomentumDispersion;     
    Float_t         fBuffer_Shape_TrackPtMean;            
    Float_t         fBuffer_Shape_TrackPtMedian;          

    Int_t           fBuffer_Jet_MC_MotherParton;          
    Int_t           fBuffer_Jet_MC_MotherHadron;          
    Int_t           fBuffer_Jet_MC_MotherIC;              
    Float_t         fBuffer_Jet_MC_MatchedJet_Distance;   
    Float_t         fBuffer_Jet_MC_MatchedJet_Pt;         
    Float_t         fBuffer_Jet_MC_MatchedJet_Mass;       
    Float_t         fBuffer_Jet_MC_TruePtFraction;        
    Float_t         fBuffer_Jet_MC_TruePtFraction_mcparticles;

    Int_t           fBuffer_NumTriggerTracks;
    Int_t           fBuffer_NumSecVertices;
    Int_t           fBuffer_NumSplittings;

    ClassDef(AliEmcalJetTree, 8) // Jet tree class
};

#endif