PWGJE/EMCALJetTasks/AliEmcalJetFinder.cxx

   1 // $Id$
   2 //
   3 // Standalone jet finder
   4 //   CINT-compatible wrapper for AliFJWrapper, can be used in macros and from the ROOT command line.
   5 //   Compiled code can use AliFJWrapper directly
   6 //
   7 // Authors: R.Haake
   8
   9 #include "AliFJWrapper.h"
  10 #include "AliEmcalJet.h"
  11 #include "AliEmcalJetFinder.h"
  12 #include "AliLog.h"
  13 #include <vector>
  14 #include "TH1.h"
  15
  16 //________________________________________________________________________
  17 AliEmcalJetFinder::AliEmcalJetFinder() :
  18   TNamed("EmcalJetFinder","EmcalJetFinder"), fFastjetWrapper(0), fInputVectorIndex(0), fJetCount(0), fJetArray(), fGhostArea(0.005), fRadius(0.4), fJetAlgorithm(0), fRecombScheme(-1), fTrackMaxEta(0.9), fJetMaxEta(0.5), fJetMinPt(0), fJetMinArea(0)
  19 {
  20   // Constructor
  21   fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
  22 }
  23
  24 //________________________________________________________________________
  25 AliEmcalJetFinder::AliEmcalJetFinder(const char* name) :
  26   TNamed(name, name), fFastjetWrapper(0), fInputVectorIndex(0), fJetCount(0), fJetArray(), fGhostArea(0.005), fRadius(0.4), fJetAlgorithm(0), fRecombScheme(-1), fTrackMaxEta(0.9), fJetMaxEta(0.5), fJetMinPt(0), fJetMinArea(0)
  27 {
  28   // Constructor
  29   fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
  30 }
  31
  32 //________________________________________________________________________
  33 AliEmcalJetFinder::~AliEmcalJetFinder()
  34 {
  35   if(fFastjetWrapper)
  36     delete fFastjetWrapper;
  37 }
  38
  39
  40 //________________________________________________________________________
  41 Bool_t AliEmcalJetFinder::FindJets()
  42 {
  43   // Tidy up and check input
  44   for (UInt_t i=0; i<fJetArray.size(); i++) {
  45     delete fJetArray[i];
  46     fJetArray[i] = 0;
  47   }
  48   fJetArray.clear();
  49   fJetCount = 0;
  50   if(!fInputVectorIndex)
  51   {
  52     AliError("No input vectors added to jet finder!");
  53     return kFALSE;
  54   }
  55
  56   // Pass settings to fastjet
  57   fFastjetWrapper->SetAreaType(fastjet::active_area_explicit_ghosts);
  58   fFastjetWrapper->SetGhostArea(fGhostArea);
  59   fFastjetWrapper->SetR(fRadius);
  60   if(fJetAlgorithm == 0)
  61     fFastjetWrapper->SetAlgorithm(fastjet::antikt_algorithm);
  62   if(fJetAlgorithm == 1)
  63     fFastjetWrapper->SetAlgorithm(fastjet::kt_algorithm);
  64   if(fRecombScheme>=0)
  65     fFastjetWrapper->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme));
  66
  67   fFastjetWrapper->SetMaxRap(fTrackMaxEta);
  68
  69   // Run jet finding
  70   fFastjetWrapper->Run();
  71
  72   // Save the found jets as light-weight objects
  73   std::vector<fastjet::PseudoJet> fastjets = fFastjetWrapper->GetInclusiveJets();
  74   fJetArray.resize(fastjets.size());
  75
  76   for (UInt_t i=0; i<fastjets.size(); i++)
  77   {
  78     // Apply jet cuts
  79     if (fastjets[i].perp()<fJetMinPt)
  80       continue;
  81     if (fFastjetWrapper->GetJetArea(i)<fJetMinArea)
  82       continue;
  83     if (TMath::Abs(fastjets[i].eta())>fJetMaxEta)
  84       continue;
  85
  86     AliEmcalJet* jet = new AliEmcalJet(fastjets[i].perp(), fastjets[i].eta(), fastjets[i].phi(), fastjets[i].m());
  87
  88     // Set the most important properties of the jet
  89     Int_t nConstituents(fFastjetWrapper->GetJetConstituents(i).size());
  90     jet->SetNumberOfTracks(nConstituents);
  91     jet->SetNumberOfClusters(nConstituents);
  92     fJetArray[fJetCount] = jet;
  93     fJetCount++;
  94   }
  95
  96   fJetArray.resize(fJetCount);
  97
  98   //fastjets.clear(); // will be done by the destructor at the end of the function
  99   fFastjetWrapper->Clear();
 100   fInputVectorIndex = 0;
 101   return kTRUE;
 102 }
 103
 104 //________________________________________________________________________
 105 void AliEmcalJetFinder::AddInputVector(Double_t px, Double_t py, Double_t pz)
 106 {
 107   fFastjetWrapper->AddInputVector(px, py, pz, TMath::Sqrt(px*px + py*py + pz*pz), fInputVectorIndex + 100);
 108   fInputVectorIndex++;
 109 }
 110
 111 //________________________________________________________________________
 112 void AliEmcalJetFinder::AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E)
 113 {
 114   fFastjetWrapper->AddInputVector(px, py, pz, E, fInputVectorIndex + 100);
 115   fInputVectorIndex++;
 116 }
 117
 118 //________________________________________________________________________
 119 void AliEmcalJetFinder::FillPtHistogram(TH1* histogram)
 120 {
 121   if(!histogram)
 122     return;
 123   for (std::size_t i=0; i<fJetArray.size(); i++)
 124   {
 125     histogram->Fill(fJetArray[i]->Pt());
 126   }
 127 }
 128
 129 //________________________________________________________________________
 130 void AliEmcalJetFinder::FillPhiHistogram(TH1* histogram)
 131 {
 132   if(!histogram)
 133     return;
 134   for (std::size_t i=0; i<fJetArray.size(); i++)
 135   {
 136     histogram->Fill(fJetArray[i]->Phi());
 137   }
 138 }
 139
 140 //________________________________________________________________________
 141 void AliEmcalJetFinder::FillEtaHistogram(TH1* histogram)
 142 {
 143   if(!histogram)
 144     return;
 145   for (std::size_t i=0; i<fJetArray.size(); i++)
 146   {
 147     histogram->Fill(fJetArray[i]->Eta());
 148   }
 149 }