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