]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGGA/EMCALJetTasks/AliFJWrapper.h
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fix from Salvatore
[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliFJWrapper.h
CommitLineData
48d874ff 1#ifndef AliFJWrapper_H
2#define AliFJWrapper_H
3
4// $Id$
5
06185cfa 6#if !defined(__CINT__) && !defined(__MAKECINT__)
48d874ff 7
8#include <vector>
9#include <TNamed.h>
10#include "AliLog.h"
11#include "FJ_includes.h"
12
13class AliFJWrapper : public TNamed
14{
15 public:
16 AliFJWrapper(const char *name, const char *title);
17 virtual ~AliFJWrapper();
18
19 virtual void AddInputVector (Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index = -99999);
20 virtual void AddInputVector (const fastjet::PseudoJet& vec, Int_t index = -99999);
21 virtual void AddInputVectors(const std::vector<fastjet::PseudoJet>& vecs, Int_t offsetIndex = -99999);
22 virtual void Clear(const Option_t* opt = "");
23 virtual void CopySettingsFrom (const AliFJWrapper& wrapper);
24 virtual void GetMedianAndSigma(Double_t& median, Double_t& sigma, Int_t remove = 0) const;
25 fastjet::ClusterSequenceArea* GetClusterSequence() const { return fClustSeq; }
26 const std::vector<fastjet::PseudoJet> GetInputVectors() const { return fInputVectors; }
27 const std::vector<fastjet::PseudoJet> GetInclusiveJets() const { return fInclusiveJets; }
28 std::vector<fastjet::PseudoJet> GetJetConstituents(UInt_t idx) const;
29 Double_t GetMedianUsedForBgSubtraction() const { return fMedUsedForBgSub; }
30 Double_t GetJetArea (UInt_t idx) const;
31 Double_t GetJetSubtractedPt (UInt_t idx) const;
32 virtual std::vector<double> GetSubtractedJetsPts(Double_t median_pt = -1, Bool_t sorted = kFALSE);
33
34 virtual Int_t Run();
35
36 void SetStrategy(const fastjet::Strategy &strat) { fStrategy = strat; }
37 void SetAlgorithm(const fastjet::JetAlgorithm &algor) { fAlgor = algor; }
38 void SetRecombScheme(const fastjet::RecombinationScheme &scheme) { fScheme = scheme; }
39 void SetAreaType(const fastjet::AreaType &atype) { fAreaType = atype; }
40 void SetNRepeats(Int_t nrepeat) { fNGhostRepeats = nrepeat; }
41 void SetGhostArea(Double_t gharea) { fGhostArea = gharea; }
42 void SetMaxRap(Double_t maxrap) { fMaxRap = maxrap; }
43 void SetR(Double_t r) { fR = r; }
44 void SetGridScatter(Double_t gridSc) { fGridScatter = gridSc; }
45 void SetKtScatter(Double_t ktSc) { fKtScatter = ktSc; }
46 void SetMeanGhostKt(Double_t meankt) { fMeanGhostKt = meankt; }
47 void SetPluginAlgor(Int_t plugin) { fPluginAlgor = plugin; }
48 void SetUseArea4Vector(Bool_t useA4v) { fUseArea4Vector = useA4v; }
49 void SetupAlgorithmfromOpt(const char *option);
50 void SetupAreaTypefromOpt(const char *option);
51 void SetupSchemefromOpt(const char *option);
52 void SetupStrategyfromOpt(const char *option);
53
54 protected:
55 std::vector<fastjet::PseudoJet> fInputVectors; //!
56 std::vector<fastjet::PseudoJet> fInclusiveJets; //!
57 std::vector<double> fSubtractedJetsPt; //!
58 fastjet::AreaDefinition *fAreaDef; //!
59 fastjet::VoronoiAreaSpec *fVorAreaSpec; //!
60 fastjet::GhostedAreaSpec *fGhostedAreaSpec; //!
61 fastjet::JetDefinition *fJetDef; //!
62 fastjet::JetDefinition::Plugin *fPlugin; //!
63 fastjet::RangeDefinition *fRange; //!
64 fastjet::ClusterSequenceArea *fClustSeq; //!
65 fastjet::Strategy fStrategy; //!
66 fastjet::JetAlgorithm fAlgor; //!
67 fastjet::RecombinationScheme fScheme; //!
68 fastjet::AreaType fAreaType; //!
69 Int_t fNGhostRepeats; //!
70 Double_t fGhostArea; //!
71 Double_t fMaxRap; //!
72 Double_t fR; //!
73 // no setters for the moment - used default values in the constructor
74 Double_t fGridScatter; //!
75 Double_t fKtScatter; //!
76 Double_t fMeanGhostKt; //!
77 Int_t fPluginAlgor; //!
78 // extra parameters
79 Double_t fMedUsedForBgSub; //!
80 Bool_t fUseArea4Vector; //!
81
82 virtual void SubtractBackground(const Double_t median_pt = -1);
83
84 private:
85 AliFJWrapper();
86 AliFJWrapper(const AliFJWrapper& wrapper);
87 AliFJWrapper& operator = (const AliFJWrapper& wrapper);
88};
89#endif
90#endif
91
92#ifdef AliFJWrapper_CXX
93#undef AliFJWrapper_CXX
94
95#if defined __GNUC__
96#pragma GCC system_header
97#endif
98
99namespace fj = fastjet;
100
101//_________________________________________________________________________________________________
102AliFJWrapper::AliFJWrapper(const char *name, const char *title)
103 : TNamed(name, title)
104 , fInputVectors ( )
105 , fInclusiveJets ( )
106 , fSubtractedJetsPt ( )
107 , fAreaDef (0)
108 , fVorAreaSpec (0)
109 , fGhostedAreaSpec (0)
110 , fJetDef (0)
111 , fPlugin (0)
112 , fRange (0)
113 , fClustSeq (0)
114 , fStrategy (fj::Best)
115 , fAlgor (fj::kt_algorithm)
116 , fScheme (fj::BIpt_scheme)
117 , fAreaType (fj::active_area)
118 , fNGhostRepeats (1)
119 , fGhostArea (0.01)
120 , fMaxRap (1.)
121 , fR (0.4)
122 , fGridScatter (1.0)
123 , fKtScatter (1.0)
124 , fMeanGhostKt (1e-100)
125 , fPluginAlgor (0)
126 , fMedUsedForBgSub (0)
127 , fUseArea4Vector (kFALSE)
128{
129 // Constructor.
130}
131
132//_________________________________________________________________________________________________
133AliFJWrapper::~AliFJWrapper()
134{
135 // Destructor.
136
137 delete fAreaDef;
138 delete fVorAreaSpec;
139 delete fGhostedAreaSpec;
140 delete fJetDef;
141 delete fPlugin;
142 delete fRange;
143 delete fClustSeq;
144}
145
146//_________________________________________________________________________________________________
147void AliFJWrapper::CopySettingsFrom(const AliFJWrapper& wrapper)
148{
149 // Copy some settings.
150 // You very often want to keep most of the settings
151 // but change only the algorithm or R - do it after call to this function
152
153 fStrategy = wrapper.fStrategy;
154 fAlgor = wrapper.fAlgor;
155 fScheme = wrapper.fScheme;
156 fAreaType = wrapper.fAreaType;
157 fNGhostRepeats = wrapper.fNGhostRepeats;
158 fGhostArea = wrapper.fGhostArea;
159 fMaxRap = wrapper.fMaxRap;
160 fR = wrapper.fR;
161 fGridScatter = wrapper.fGridScatter;
162 fKtScatter = wrapper.fKtScatter;
163 fMeanGhostKt = wrapper.fMeanGhostKt;
164 fPluginAlgor = wrapper.fPluginAlgor;
165 fUseArea4Vector = wrapper.fUseArea4Vector;
166}
167
168//_________________________________________________________________________________________________
169void AliFJWrapper::Clear(const Option_t *opt)
170{
171 // Simply clear the input vectors.
172 // Make sure done on every event if the instance is reused
173 // Reset the median to zero.
174
175 fInputVectors.clear();
176 fMedUsedForBgSub = 0;
177
178 // for the moment brute force delete everything
179 delete fAreaDef; fAreaDef = 0;
180 delete fVorAreaSpec; fVorAreaSpec = 0;
181 delete fGhostedAreaSpec; fGhostedAreaSpec = 0;
182 delete fJetDef; fJetDef = 0;
183 delete fPlugin; fPlugin = 0;
184 delete fRange; fRange = 0;
185 delete fClustSeq; fClustSeq = 0;
186
187 TNamed::Clear(opt);
188}
189
190//_________________________________________________________________________________________________
191void AliFJWrapper::AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index)
192{
193 // Make the input pseudojet.
194
195 fastjet::PseudoJet inVec(px, py, pz, E);
196
197 if (index > -99999) {
198 inVec.set_user_index(index);
199 } else {
200 inVec.set_user_index(fInputVectors.size());
201 }
202
203 // add to the fj container of input vectors
204 fInputVectors.push_back(inVec);
205}
206
207//_________________________________________________________________________________________________
208void AliFJWrapper::AddInputVector(const fj::PseudoJet& vec, Int_t index)
209{
210 // Add an input pseudojet.
211
212 fj::PseudoJet inVec = vec;
213
214 if (index > -99999) {
215 inVec.set_user_index(index);
216 } else {
217 inVec.set_user_index(fInputVectors.size());
218 }
219
220 // add to the fj container of input vectors
221 fInputVectors.push_back(inVec);
222}
223
224//_________________________________________________________________________________________________
225void AliFJWrapper::AddInputVectors(const std::vector<fj::PseudoJet>& vecs, Int_t offsetIndex)
226{
227 // Add the input from vector of pseudojets.
228
229 for (UInt_t i = 0; i < vecs.size(); ++i) {
230 fj::PseudoJet inVec = vecs[i];
231
232 if (offsetIndex > -99999)
233 inVec.set_user_index(fInputVectors.size() + offsetIndex);
234 // add to the fj container of input vectors
235
236 fInputVectors.push_back(inVec);
237 }
238}
239
240//_________________________________________________________________________________________________
241Double_t AliFJWrapper::GetJetArea(UInt_t idx) const
242{
243 // Get the jet area.
244
245 Double_t retval = -1; // really wrong area..
246 if ( idx < fInclusiveJets.size() ) {
247 retval = fClustSeq->area(fInclusiveJets[idx]);
248 } else {
249 AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
250 }
251 return retval;
252}
253
254//_________________________________________________________________________________________________
255std::vector<double> AliFJWrapper::GetSubtractedJetsPts(Double_t median_pt, Bool_t sorted)
256{
257 // Get subtracted jets pTs, returns vector.
258
259 SubtractBackground(median_pt);
260
261 if (kTRUE == sorted) {
262 std::sort(fSubtractedJetsPt.begin(), fSubtractedJetsPt.begin());
263 }
264 return fSubtractedJetsPt;
265}
266
267//_________________________________________________________________________________________________
268Double_t AliFJWrapper::GetJetSubtractedPt(UInt_t idx) const
269{
270 // Get subtracted jets pTs, returns Double_t.
271
272 Double_t retval = -99999.; // really wrong pt..
273 if ( idx < fSubtractedJetsPt.size() ) {
274 retval = fSubtractedJetsPt[idx];
275 }
276 return retval;
277}
278
279//_________________________________________________________________________________________________
280std::vector<fastjet::PseudoJet>
281AliFJWrapper::GetJetConstituents(UInt_t idx) const
282{
283 // Get jets constituents.
284
285 std::vector<fastjet::PseudoJet> retval;
286
287 if ( idx < fInclusiveJets.size() ) {
288 retval = fClustSeq->constituents(fInclusiveJets[idx]);
289 } else {
290 AliError(Form("[e] ::GetJetConstituents wrong index: %d",idx));
291 }
292
293 return retval;
294}
295
296//_________________________________________________________________________________________________
297void AliFJWrapper::GetMedianAndSigma(Double_t &median, Double_t &sigma, Int_t remove) const
298{
299 // Get the median and sigma from fastjet.
300 // User can also do it on his own because the cluster sequence is exposed (via a getter)
301
302 if (!fClustSeq) {
303 AliError("[e] Run the jfinder first.");
304 }
305
306 Double_t mean_area = 0;
307 try {
308 if(0 == remove) {
309 fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
310 } else {
311 std::vector<fastjet::PseudoJet> input_jets = sorted_by_pt(fClustSeq->inclusive_jets());
312 input_jets.erase(input_jets.begin(), input_jets.begin() + remove);
313 fClustSeq->get_median_rho_and_sigma(input_jets, *fRange, kFALSE, median, sigma, mean_area);
314 input_jets.clear();
315 }
316 } catch (fj::Error) {
317 AliError(" [w] FJ Exception caught.");
318 median = -1.;
319 sigma = -1;
320 }
321}
322
323//_________________________________________________________________________________________________
324Int_t AliFJWrapper::Run()
325{
326 // Run the actual jet finder.
327
328 if (fAreaType == fj::voronoi_area) {
329 // Rfact - check dependence - default is 1.
330 // NOTE: hardcoded variable!
331 fVorAreaSpec = new fj::VoronoiAreaSpec(1.);
332 fAreaDef = new fj::AreaDefinition(*fVorAreaSpec);
333 } else {
334 fGhostedAreaSpec = new fj::GhostedAreaSpec(fMaxRap,
335 fNGhostRepeats,
336 fGhostArea,
337 fGridScatter,
338 fKtScatter,
339 fMeanGhostKt);
340
341 fAreaDef = new fj::AreaDefinition(*fGhostedAreaSpec, fAreaType);
342 }
343
344 // this is acceptable by fastjet:
345 fRange = new fj::RangeDefinition(fMaxRap - 0.95 * fR);
346
347 if (fAlgor == fj::plugin_algorithm) {
348 if (fPluginAlgor == 0) {
349 // SIS CONE ALGOR
350 // NOTE: hardcoded split parameter
351 Double_t overlap_threshold = 0.75; // NOTE: this actually splits a lot: thr/min(pt1,pt2)
352 fPlugin = new fj::SISConePlugin(fR,
353 overlap_threshold,
354 0, //search of stable cones - zero = until no more
355 1.0); // this should be seed effectively for proto jets
356 fJetDef = new fastjet::JetDefinition(fPlugin);
357 } else {
358 AliError("[e] Unrecognized plugin number!");
359 }
360 } else {
361 fJetDef = new fj::JetDefinition(fAlgor, fR, fScheme, fStrategy);
362 }
363
364 try {
365 fClustSeq = new fj::ClusterSequenceArea(fInputVectors, *fJetDef, *fAreaDef);
366 } catch (fj::Error) {
367 AliError(" [w] FJ Exception caught.");
368 return -1;
369 }
370
371 // inclusive jets:
372 fInclusiveJets.clear();
373 fInclusiveJets = fClustSeq->inclusive_jets(0.0); //sorted_by_pt(fClustSeq->inclusive_jets(0.0));
374
375 return 0;
376}
377
378//_________________________________________________________________________________________________
379void AliFJWrapper::SubtractBackground(Double_t median_pt)
380{
381 // Subtract the background (specify the value - see below the meaning).
382 // Negative argument means the bg will be determined with the current algorithm
383 // this is the default behavior. Zero is allowed
384 // Note: user may set the switch for area4vector based subtraction.
385
386 Double_t median = 0;
387 Double_t sigma = 0;
388 Double_t mean_area = 0;
389
390 // clear the subtracted jet pt's vector<double>
391 fSubtractedJetsPt.clear();
392
393 // check what was specified (default is -1)
394 if (median_pt < 0) {
395 try {
396 fClustSeq->get_median_rho_and_sigma(*fRange, kFALSE, median, sigma, mean_area);
397 }
398
399 catch (fj::Error) {
400 AliError(" [w] FJ Exception caught.");
401 median = -9999.;
402 sigma = -1;
403 fMedUsedForBgSub = median;
404 return;
405 }
406 fMedUsedForBgSub = median;
407 } else {
408 // we do not know the sigma in this case
409 sigma = -1;
410 if (0.0 == median_pt) {
411 // AliWarning(Form(" [w] Median specified for bg subtraction is ZERO: %f", median_pt )) << endl;
412 fMedUsedForBgSub = 0.;
413 } else {
414 fMedUsedForBgSub = median_pt;
415 }
416 }
417
418 // subtract:
419 for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
420 if (kTRUE == fUseArea4Vector) {
421 // subtract the background using the area4vector
422 fj::PseudoJet area4v = fClustSeq->area_4vector(fInclusiveJets[i]);
423 fj::PseudoJet jet_sub = fInclusiveJets[i] - area4v * fMedUsedForBgSub;
424 fSubtractedJetsPt.push_back(jet_sub.perp()); // here we put only the pt of the jet - note: this can be negative
425 } else {
426 // subtract the background using scalars
427 // fj::PseudoJet jet_sub = fInclusiveJets[i] - area * fMedUsedForBgSub_;
428 Double_t area = fClustSeq->area(fInclusiveJets[i]);
429 // standard subtraction
430 Double_t pt_sub = fInclusiveJets[i].perp() - fMedUsedForBgSub * area;
431 fSubtractedJetsPt.push_back(pt_sub); // here we put only the pt of the jet - note: this can be negative
432 }
433 }
434}
435
436//_________________________________________________________________________________________________
437void AliFJWrapper::SetupAlgorithmfromOpt(const char *option)
438{
439 // Setup algorithm from char.
440
441 std::string opt(option);
442
443 if (!opt.compare("kt")) fAlgor = fj::kt_algorithm;
444 if (!opt.compare("antikt")) fAlgor = fj::antikt_algorithm;
445 if (!opt.compare("cambridge")) fAlgor = fj::cambridge_algorithm;
446 if (!opt.compare("genkt")) fAlgor = fj::genkt_algorithm;
447 if (!opt.compare("cambridge_passive")) fAlgor = fj::cambridge_for_passive_algorithm;
448 if (!opt.compare("genkt_passive")) fAlgor = fj::genkt_for_passive_algorithm;
449 if (!opt.compare("ee_kt")) fAlgor = fj::ee_kt_algorithm;
450 if (!opt.compare("ee_genkt")) fAlgor = fj::ee_genkt_algorithm;
451 if (!opt.compare("plugin")) fAlgor = fj::plugin_algorithm;
452}
453
454//_________________________________________________________________________________________________
455void AliFJWrapper::SetupAreaTypefromOpt(const char *option)
456{
457 // Setup area type from char.
458
459 std::string opt(option);
460
461 if (!opt.compare("active")) fAreaType = fj::active_area;
462 if (!opt.compare("invalid")) fAreaType = fj::invalid_area;
463 if (!opt.compare("active_area_explicit_ghosts")) fAreaType = fj::active_area_explicit_ghosts;
464 if (!opt.compare("one_ghost_passive")) fAreaType = fj::one_ghost_passive_area;
465 if (!opt.compare("passive")) fAreaType = fj::passive_area;
466 if (!opt.compare("voronoi")) fAreaType = fj::voronoi_area;
467}
468
469//_________________________________________________________________________________________________
470void AliFJWrapper::SetupSchemefromOpt(const char *option)
471{
472 //
473 // setup scheme from char
474 //
475
476 std::string opt(option);
477
478 if (!opt.compare("BIpt")) fScheme = fj::BIpt_scheme;
479 if (!opt.compare("BIpt2")) fScheme = fj::BIpt2_scheme;
480 if (!opt.compare("E")) fScheme = fj::E_scheme;
481 if (!opt.compare("pt")) fScheme = fj::pt_scheme;
482 if (!opt.compare("pt2")) fScheme = fj::pt2_scheme;
483 if (!opt.compare("Et")) fScheme = fj::Et_scheme;
484 if (!opt.compare("Et2")) fScheme = fj::Et2_scheme;
485}
486
487//_________________________________________________________________________________________________
488void AliFJWrapper::SetupStrategyfromOpt(const char *option)
489{
490 // Setup strategy from char.
491
492 std::string opt(option);
493
494 if (!opt.compare("Best")) fStrategy = fj::Best;
495 if (!opt.compare("N2MinHeapTiled")) fStrategy = fj::N2MinHeapTiled;
496 if (!opt.compare("N2Tiled")) fStrategy = fj::N2Tiled;
497 if (!opt.compare("N2PoorTiled")) fStrategy = fj::N2PoorTiled;
498 if (!opt.compare("N2Plain")) fStrategy = fj::N2Plain;
499 if (!opt.compare("N3Dumb")) fStrategy = fj::N3Dumb;
500 if (!opt.compare("NlnN")) fStrategy = fj::NlnN;
501 if (!opt.compare("NlnN3pi")) fStrategy = fj::NlnN3pi;
502 if (!opt.compare("NlnN4pi")) fStrategy = fj::NlnN4pi;
503 if (!opt.compare("NlnNCam4pi")) fStrategy = fj::NlnNCam4pi;
504 if (!opt.compare("NlnNCam2pi2R")) fStrategy = fj::NlnNCam2pi2R;
505 if (!opt.compare("NlnNCam")) fStrategy = fj::NlnNCam;
506 if (!opt.compare("plugin")) fStrategy = fj::plugin_strategy;
507}
508#endif