AliPhysics  aaf9c62 (aaf9c62)
AliEmcalJetFinder.cxx
Go to the documentation of this file.
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 <vector>
10 #include <TLorentzVector.h>
11 
12 #include "AliLog.h"
13 #include "AliVTrack.h"
14 #include "AliVCluster.h"
15 #include "AliEmcalJet.h"
16 #include "AliParticleContainer.h"
17 #include "AliClusterContainer.h"
18 #include "AliJetContainer.h"
19 #include "AliFJWrapper.h"
20 #include "AliEmcalJetFinder.h"
21 #include "TH1.h"
22 
23 //________________________________________________________________________
25  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)
26 {
27  // Constructor
28  fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
29 }
30 
31 //________________________________________________________________________
34 {
35  // Constructor
36  fFastjetWrapper = new AliFJWrapper("FJWrapper", "FJWrapper");
37 }
38 
39 //________________________________________________________________________
41 {
42  if(fFastjetWrapper)
43  delete fFastjetWrapper;
44 }
45 
46 
47 //________________________________________________________________________
49 {
50  // Tidy up and check input
51  for (UInt_t i=0; i<fJetArray.size(); i++) {
52  delete fJetArray[i];
53  fJetArray[i] = 0;
54  }
55  fJetArray.clear();
56  fJetCount = 0;
58  {
59  AliError("No input vectors added to jet finder!");
60  return kFALSE;
61  }
62 
63  // Pass settings to fastjet
64  fFastjetWrapper->SetAreaType(fastjet::active_area_explicit_ghosts);
67  if(fJetAlgorithm == 0)
68  fFastjetWrapper->SetAlgorithm(fastjet::antikt_algorithm);
69  if(fJetAlgorithm == 1)
70  fFastjetWrapper->SetAlgorithm(fastjet::kt_algorithm);
71  if(fRecombScheme>=0)
72  fFastjetWrapper->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme));
73 
75 
76  // Run jet finding
78 
79  // Save the found jets as light-weight objects
80  std::vector<fastjet::PseudoJet> fastjets = fFastjetWrapper->GetInclusiveJets();
81  fJetArray.resize(fastjets.size());
82 
83  for (UInt_t i=0; i<fastjets.size(); i++)
84  {
85  // Apply jet cuts
86  if (fastjets[i].perp()<fJetMinPt)
87  continue;
89  continue;
90  if (TMath::Abs(fastjets[i].eta())>fJetMaxEta)
91  continue;
92 
93  AliEmcalJet* jet = new AliEmcalJet(fastjets[i].perp(), fastjets[i].eta(), fastjets[i].phi(), fastjets[i].m());
94 
95  // Set the most important properties of the jet
96  Int_t nConstituents(fFastjetWrapper->GetJetConstituents(i).size());
98  jet->SetNumberOfTracks(nConstituents);
99  jet->SetNumberOfClusters(nConstituents);
100  fJetArray[fJetCount] = jet;
101  fJetCount++;
102  }
103 
104  fJetArray.resize(fJetCount);
105 
106  //fastjets.clear(); // will be done by the destructor at the end of the function
108  fInputVectorIndex = 0;
109  return kTRUE;
110 }
111 
112 //________________________________________________________________________
114 {
115 //
116 // AliEmcalJetFinder::Filter
117 //
118 
119  fJetCount = 0;
120  for (UInt_t j=0; j<fJetArray.size(); j++) {
121  if (fJetArray[j]) { delete fJetArray[j]; fJetArray[j] = 0; }
122  } fJetArray.clear();
123 //=============================================================================
124 
125  if ((!pJet) || (!pContJets)) return kFALSE;
126 //=============================================================================
127 
128  AliParticleContainer *pContTrks = pContJets->GetParticleContainer();
129  if (pContTrks) for (Int_t i=0; i<pJet->GetNumberOfTracks(); i++) {
130  AliVParticle *pTrk = pJet->TrackAt(i, pContTrks->GetArray()); if (!pTrk) continue;
131  AddInputVector(pTrk->Px(), pTrk->Py(), pTrk->Pz(), pTrk->E(), pJet->TrackAt(i)+100);
132  }
133 
134  AliClusterContainer *pContClus = pContJets->GetClusterContainer();
135  if (pContClus) for (Int_t i=0; i<pJet->GetNumberOfClusters(); i++) {
136  AliVCluster *pClu = pJet->ClusterAt(i, pContClus->GetArray()); if (!pClu) continue;
137 
138  TLorentzVector vClu; pClu->GetMomentum(vClu, dVtx);
139  AddInputVector(vClu.Px(), vClu.Py(), vClu.Pz(), vClu.P(), -1*pJet->ClusterAt(i)-100);
140  }
141 //=============================================================================
142 
143  if(!fInputVectorIndex) {
144  AliError("No input vectors added to jet finder!");
145  return kFALSE;
146  }
147 //=============================================================================
148 
149  if (pJet->HasGhost()) {
150  const std::vector<TLorentzVector> aGhosts = pJet->GetGhosts();
151  for (UInt_t i=0; i<aGhosts.size(); i++) AddInputGhost(aGhosts[i].Px(),
152  aGhosts[i].Py(),
153  aGhosts[i].Pz(),
154  aGhosts[i].E());
155  }
156 //=============================================================================
157 
161  if(fJetAlgorithm==0) fFastjetWrapper->SetAlgorithm(fastjet::antikt_algorithm);
162  if(fJetAlgorithm==1) fFastjetWrapper->SetAlgorithm(fastjet::kt_algorithm);
163  if(fRecombScheme>=0) fFastjetWrapper->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme));
164 //=============================================================================
165 
167  std::vector<fastjet::PseudoJet> aFilteredJets = fFastjetWrapper->GetFilteredJets();
168 
169  fJetArray.resize(aFilteredJets.size());
170  for (UInt_t j=0; j<aFilteredJets.size(); j++) {
171  if (aFilteredJets[j].perp()<fJetMinPt) continue;
172  if (TMath::Abs(aFilteredJets[j].eta())>fJetMaxEta) continue;
174 
175  AliEmcalJet *piece = new AliEmcalJet(aFilteredJets[j].perp(),
176  aFilteredJets[j].eta(),
177  aFilteredJets[j].phi(),
178  aFilteredJets[j].m());
179 
180  piece->SetLabel(j);
182  fastjet::PseudoJet area(fFastjetWrapper->GetFilteredJetAreaVector(j));
183  piece->SetArea(area.perp());
184  piece->SetAreaEta(area.eta());
185  piece->SetAreaPhi(area.phi());
186  piece->SetAreaE(area.E());
187  }
188 
189  std::vector<fastjet::PseudoJet> aConstis(fFastjetWrapper->GetFilteredJetConstituents(j));
190 
191  UInt_t nConstis = aConstis.size();
192  piece->SetNumberOfTracks(nConstis);
193  piece->SetNumberOfClusters(nConstis);
194 
195  Int_t nt = 0;
196  Int_t nc = 0;
197  for (UInt_t i=0; i<nConstis; i++) {
198  Int_t uid = aConstis[i].user_index();
199  if (uid>= 100) { piece->AddTrackAt( 1*uid-100, nt); ++nt; }
200  if (uid<=-100) { piece->AddClusterAt(-1*uid-100, nc); ++nc; }
201  }
202 
203  piece->SetNumberOfTracks(nt);
204  piece->SetNumberOfClusters(nc);
205 
206  fJetArray[fJetCount] = piece;
207  fJetCount++;
208  }
209 
210  fJetArray.resize(fJetCount);
212  fInputVectorIndex = 0;
213  return kTRUE;
214 }
215 
216 //________________________________________________________________________
218 {
219  fFastjetWrapper->AddInputVector(px, py, pz, TMath::Sqrt(px*px + py*py + pz*pz), fInputVectorIndex + 100);
221 }
222 
223 //________________________________________________________________________
225 {
226  fFastjetWrapper->AddInputVector(px, py, pz, E, fInputVectorIndex + 100);
228 }
229 
230 //________________________________________________________________________
232 {
233  fFastjetWrapper->AddInputVector(px, py, pz, TMath::Sqrt(px*px + py*py + pz*pz), index);
235 }
236 
237 //________________________________________________________________________
239 {
240  fFastjetWrapper->AddInputVector(px, py, pz, E, index);
242 }
243 
244 //________________________________________________________________________
246 {
247  fFastjetWrapper->AddInputGhost(px, py, pz, E, -1);
248 }
249 
250 //________________________________________________________________________
252 {
253  if(!histogram)
254  return;
255  for (std::size_t i=0; i<fJetArray.size(); i++)
256  {
257  histogram->Fill(fJetArray[i]->Pt());
258  }
259 }
260 
261 //________________________________________________________________________
263 {
264  if(!histogram)
265  return;
266  for (std::size_t i=0; i<fJetArray.size(); i++)
267  {
268  histogram->Fill(fJetArray[i]->Phi());
269  }
270 }
271 
272 //________________________________________________________________________
274 {
275  if(!histogram)
276  return;
277  for (std::size_t i=0; i<fJetArray.size(); i++)
278  {
279  histogram->Fill(fJetArray[i]->Eta());
280  }
281 }
282 
283 
284 
285 //________________________________________________________________________
286 Double_t AliEmcalJetFinder::Nsubjettiness(AliEmcalJet *pJet, AliJetContainer *pContJets, Double_t dVtx[3], Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Int_t Option, Int_t Measure, Double_t Beta_SD, Double_t ZCut, Int_t SoftDropOn){
287 
288  fJetCount = 0;
289  for (UInt_t j=0; j<fJetArray.size(); j++) {
290  if (fJetArray[j]) { delete fJetArray[j]; fJetArray[j] = 0; }
291  } fJetArray.clear();
292 //=============================================================================
293 
294  if ((!pJet) || (!pContJets)) return kFALSE;
295 //=============================================================================
296 
297  AliParticleContainer *pContTrks = pContJets->GetParticleContainer();
298  if (pContTrks) for (Int_t i=0; i<pJet->GetNumberOfTracks(); i++) {
299  AliVParticle *pTrk = pJet->TrackAt(i, pContTrks->GetArray()); if (!pTrk) continue;
300  AddInputVector(pTrk->Px(), pTrk->Py(), pTrk->Pz(), pTrk->E(), pJet->TrackAt(i)+100);
301  }
302  AliClusterContainer *pContClus = pContJets->GetClusterContainer();
303  if (pContClus) for (Int_t i=0; i<pJet->GetNumberOfClusters(); i++) {
304  AliVCluster *pClu = pJet->ClusterAt(i, pContClus->GetArray()); if (!pClu) continue;
305 
306  TLorentzVector vClu; pClu->GetMomentum(vClu, dVtx);
307  AddInputVector(vClu.Px(), vClu.Py(), vClu.Pz(), vClu.P(), -1*pJet->ClusterAt(i)-100);
308  }
309 //=============================================================================
310 
311  if(!fInputVectorIndex) {
312  AliError("No input vectors added to jet finder!");
313  }
314 
315 //=============================================================================
316 
317  if (pJet->HasGhost()) {
318  const std::vector<TLorentzVector> aGhosts = pJet->GetGhosts();
319  for (UInt_t i=0; i<aGhosts.size(); i++) AddInputGhost(aGhosts[i].Px(),
320  aGhosts[i].Py(),
321  aGhosts[i].Pz(),
322  aGhosts[i].E());
323  }
324  //these are all for the jet not subjets
329  if(fJetAlgorithm==0) fFastjetWrapper->SetAlgorithm(fastjet::antikt_algorithm); //this is for the jet clustering not the subjet reclustering.
330  // if(fJetAlgorithm==1) fFastjetWrapper->SetAlgorithm(fastjet::kt_algorithm);
331  if(fRecombScheme>=0) fFastjetWrapper->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme));
332  if (Measure==1) return fFastjetWrapper->AliFJWrapper::NSubjettiness(N,Algorithm,Radius, Beta, Option,1);
333  else return fFastjetWrapper->AliFJWrapper::NSubjettiness(N,Algorithm,Radius, Beta, Option,0,Beta_SD,ZCut,SoftDropOn);
334 
335 }
336 
337 
338 
339 
340 
Double_t GetFilteredJetArea(UInt_t idx) const
void AddTrackAt(Int_t track, Int_t idx)
Definition: AliEmcalJet.h:275
void AddInputVector(Double_t px, Double_t py, Double_t pz)
virtual Int_t Run()
double Double_t
Definition: External.C:58
virtual Int_t Filter()
Bool_t Filter(AliEmcalJet *pJet, AliJetContainer *pContJets, Double_t dVtx[3])
void AddInputGhost(Double_t px, Double_t py, Double_t pz, Double_t E)
void AddClusterAt(Int_t clus, Int_t idx)
Definition: AliEmcalJet.h:274
Int_t ClusterAt(Int_t idx) const
Definition: AliEmcalJet.h:137
AliClusterContainer * GetClusterContainer() const
void SetArea(Double_t a)
Definition: AliEmcalJet.h:256
void SetMinJetPt(Double_t MinPt)
Definition: AliFJWrapper.h:142
void SetMaxRap(Double_t maxrap)
Definition: AliFJWrapper.h:126
fastjet::PseudoJet GetFilteredJetAreaVector(UInt_t idx) const
void SetAreaE(Double_t a)
Definition: AliEmcalJet.h:259
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
void FillPhiHistogram(TH1 *histogram)
void SetLabel(Int_t l)
Definition: AliEmcalJet.h:255
void FillPtHistogram(TH1 *histogram)
Container for particles within the EMCAL framework.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
Double_t Nsubjettiness(AliEmcalJet *pJet, AliJetContainer *pContJets, Double_t dVtx[3], Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Int_t Option=0, Int_t Measure=0, Double_t Beta_SD=0, Double_t ZCut=0.1, Int_t SoftDropOn=0)
AliParticleContainer * GetParticleContainer() const
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
int Int_t
Definition: External.C:63
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:138
unsigned int UInt_t
Definition: External.C:33
Double_t GetJetArea(UInt_t idx) const
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
Definition: Option.C:68
virtual void Clear(const Option_t *="")
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:125
Bool_t HasGhost() const
Definition: AliEmcalJet.h:342
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
void SetAreaEta(Double_t a)
Definition: AliEmcalJet.h:257
AliFJWrapper * fFastjetWrapper
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
void SetNumberOfTracks(Int_t n)
Definition: AliEmcalJet.h:266
const std::vector< TLorentzVector > GetGhosts() const
Definition: AliEmcalJet.h:343
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
virtual void AddInputGhost(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
std::vector< AliEmcalJet * > fJetArray
bool Bool_t
Definition: External.C:53
std::vector< fastjet::PseudoJet > GetFilteredJetConstituents(UInt_t idx) const
void SetAreaPhi(Double_t a)
Definition: AliEmcalJet.h:258
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:123
Bool_t GetDoFilterArea()
Definition: AliFJWrapper.h:52
Container structure for EMCAL clusters.
Container for jet within the EMCAL jet framework.
Definition: External.C:196
void SetNumberOfClusters(Int_t n)
Definition: AliEmcalJet.h:265
const std::vector< fastjet::PseudoJet > & GetFilteredJets() const
Definition: AliFJWrapper.h:36
void FillEtaHistogram(TH1 *histogram)