PWGJE/EMCALJetTasks/AliEmcalJetFinder.cxx

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