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), fRecombScheme(-1), 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), fRecombScheme(-1), 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   if(fRecombScheme>0)
  59     fFastjetWrapper->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme));
  60
  61   fFastjetWrapper->SetMaxRap(fTrackMaxEta);
  62
  63   // Run jet finding
  64   fFastjetWrapper->Run();
  65
  66   // Save the found jets as light-weight objects
  67   std::vector<fastjet::PseudoJet> fastjets = fFastjetWrapper->GetInclusiveJets();
  68   fJetArray.resize(fastjets.size());
  69
  70   for (UInt_t i=0; i<fastjets.size(); i++)
  71   {
  72     // Apply jet cuts
  73     if (fastjets[i].perp()<fJetMinPt)
  74       continue;
  75     if (fFastjetWrapper->GetJetArea(i)<fJetMinArea)
  76       continue;
  77     if (TMath::Abs(fastjets[i].eta())>fJetMaxEta)
  78       continue;
  79
  80     AliEmcalJet* jet = new AliEmcalJet(fastjets[i].perp(), fastjets[i].eta(), fastjets[i].phi(), fastjets[i].m());
  81
  82     // Set the most important properties of the jet
  83     Int_t nConstituents(fFastjetWrapper->GetJetConstituents(i).size());
  84     jet->SetNumberOfTracks(nConstituents);
  85     jet->SetNumberOfClusters(nConstituents);
  86     fJetArray[fJetCount] = jet;
  87     fJetCount++;
  88   }
  89
  90   fJetArray.resize(fJetCount);
  91
  92   fastjets.clear();
  93   fFastjetWrapper->Clear();
  94   fInputVectorIndex = 0;
  95   return kTRUE;
  96 }
  97
  98 //________________________________________________________________________
  99 void AliEmcalJetFinder::AddInputVector(Double_t px, Double_t py, Double_t pz)
 100 {
 101   fFastjetWrapper->AddInputVector(px, py, pz, TMath::Sqrt(px*px + py*py + pz*pz), fInputVectorIndex + 100);
 102   fInputVectorIndex++;
 103 }
 104
 105 //________________________________________________________________________
 106 void AliEmcalJetFinder::AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E)
 107 {
 108   fFastjetWrapper->AddInputVector(px, py, pz, E, fInputVectorIndex + 100);
 109   fInputVectorIndex++;
 110 }
 111
 112 //________________________________________________________________________
 113 void AliEmcalJetFinder::FillPtHistogram(TH1* histogram)
 114 {
 115   if(!histogram)
 116     return;
 117   for (Int_t i=0; i<fJetArray.size(); i++)
 118   {
 119     histogram->Fill(fJetArray[i]->Pt());
 120   }
 121 }
 122
 123 //________________________________________________________________________
 124 void AliEmcalJetFinder::FillPhiHistogram(TH1* histogram)
 125 {
 126   if(!histogram)
 127     return;
 128   for (Int_t i=0; i<fJetArray.size(); i++)
 129   {
 130     histogram->Fill(fJetArray[i]->Phi());
 131   }
 132 }
 133
 134 //________________________________________________________________________
 135 void AliEmcalJetFinder::FillEtaHistogram(TH1* histogram)
 136 {
 137   if(!histogram)
 138     return;
 139   for (Int_t i=0; i<fJetArray.size(); i++)
 140   {
 141     histogram->Fill(fJetArray[i]->Eta());
 142   }
 143 }