]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/EMCALJetTasks/AliFJWrapper.h
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Adapt add macro and particle/cluster containers for the analysis on jets
[u/mrichter/AliRoot.git] / PWGJE / 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"
8082a80b 12#include "AliJetShape.h"
48d874ff 13
a9c8e539 14class AliFJWrapper
48d874ff 15{
16 public:
17 AliFJWrapper(const char *name, const char *title);
18 virtual ~AliFJWrapper();
a9c8e539 19
48d874ff 20 virtual void AddInputVector (Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index = -99999);
21 virtual void AddInputVector (const fastjet::PseudoJet& vec, Int_t index = -99999);
22 virtual void AddInputVectors(const std::vector<fastjet::PseudoJet>& vecs, Int_t offsetIndex = -99999);
a9c8e539 23 virtual const char *ClassName() const { return "AliFJWrapper"; }
24 virtual void Clear(const Option_t* /*opt*/ = "");
48d874ff 25 virtual void CopySettingsFrom (const AliFJWrapper& wrapper);
26 virtual void GetMedianAndSigma(Double_t& median, Double_t& sigma, Int_t remove = 0) const;
8082a80b 27 fastjet::ClusterSequenceArea* GetClusterSequence() const { return fClustSeq; }
28 const std::vector<fastjet::PseudoJet> GetInputVectors() const { return fInputVectors; }
29 const std::vector<fastjet::PseudoJet> GetInclusiveJets() const { return fInclusiveJets; }
30 std::vector<fastjet::PseudoJet> GetJetConstituents(UInt_t idx) const;
31 Double_t GetMedianUsedForBgSubtraction() const { return fMedUsedForBgSub; }
32 const char* GetName() const { return fName; }
33 const char* GetTitle() const { return fTitle; }
34 Double_t GetJetArea (UInt_t idx) const;
35 fastjet::PseudoJet GetJetAreaVector (UInt_t idx) const;
36 Double_t GetJetSubtractedPt (UInt_t idx) const;
37 virtual std::vector<double> GetSubtractedJetsPts(Double_t median_pt = -1, Bool_t sorted = kFALSE);
38 Bool_t GetLegacyMode() { return fLegacyMode; }
93daf3f7 39#ifdef FASTJET_VERSION
3ec5da8e 40 const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetMass() const {return fGenSubtractorInfoJetMass ; }
41 const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetAngularity() const {return fGenSubtractorInfoJetAngularity ; }
42 const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetpTD() const {return fGenSubtractorInfoJetpTD ; }
43 const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetCircularity() const {return fGenSubtractorInfoJetCircularity ; }
df20822b 44 const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetSigma2() const {return fGenSubtractorInfoJetSigma2 ; }
3ec5da8e 45 const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetConstituent() const {return fGenSubtractorInfoJetConstituent ; }
46 const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetLeSub() const {return fGenSubtractorInfoJetLeSub ; }
47 const std::vector<fastjet::PseudoJet> GetConstituentSubtrJets() const {return fConstituentSubtrJets ; }
8082a80b 48#endif
3ec5da8e 49 virtual std::vector<double> GetGRNumerator() const { return fGRNumerator ; }
50 virtual std::vector<double> GetGRDenominator() const { return fGRDenominator ; }
51 virtual std::vector<double> GetGRNumeratorSub() const { return fGRNumeratorSub ; }
52 virtual std::vector<double> GetGRDenominatorSub() const { return fGRDenominatorSub ; }
48d874ff 53
54 virtual Int_t Run();
8082a80b 55 virtual Int_t DoGenericSubtractionJetMass();
726c9488 56 virtual Int_t DoGenericSubtractionGR(Int_t ijet);
0d13a63c 57 virtual Int_t DoGenericSubtractionJetAngularity();
58 virtual Int_t DoGenericSubtractionJetpTD();
59 virtual Int_t DoGenericSubtractionJetCircularity();
b9ccc456 60 virtual Int_t DoGenericSubtractionJetSigma2();
0d13a63c 61 virtual Int_t DoGenericSubtractionJetConstituent();
62 virtual Int_t DoGenericSubtractionJetLeSub();
8082a80b 63 virtual Int_t DoConstituentSubtraction();
48d874ff 64
65 void SetStrategy(const fastjet::Strategy &strat) { fStrategy = strat; }
66 void SetAlgorithm(const fastjet::JetAlgorithm &algor) { fAlgor = algor; }
67 void SetRecombScheme(const fastjet::RecombinationScheme &scheme) { fScheme = scheme; }
68 void SetAreaType(const fastjet::AreaType &atype) { fAreaType = atype; }
69 void SetNRepeats(Int_t nrepeat) { fNGhostRepeats = nrepeat; }
70 void SetGhostArea(Double_t gharea) { fGhostArea = gharea; }
71 void SetMaxRap(Double_t maxrap) { fMaxRap = maxrap; }
72 void SetR(Double_t r) { fR = r; }
73 void SetGridScatter(Double_t gridSc) { fGridScatter = gridSc; }
74 void SetKtScatter(Double_t ktSc) { fKtScatter = ktSc; }
75 void SetMeanGhostKt(Double_t meankt) { fMeanGhostKt = meankt; }
76 void SetPluginAlgor(Int_t plugin) { fPluginAlgor = plugin; }
77 void SetUseArea4Vector(Bool_t useA4v) { fUseArea4Vector = useA4v; }
78 void SetupAlgorithmfromOpt(const char *option);
79 void SetupAreaTypefromOpt(const char *option);
80 void SetupSchemefromOpt(const char *option);
81 void SetupStrategyfromOpt(const char *option);
5d0a5007 82 void SetLegacyMode (Bool_t mode) { fLegacyMode ^= mode; }
83 void SetLegacyFJ();
8082a80b 84 void SetUseExternalBkg(Bool_t b, Double_t rho, Double_t rhom) { fUseExternalBkg = b; fRho = rho; fRhom = rhom;}
726c9488 85 void SetRMaxAndStep(Double_t rmax, Double_t dr) {fRMax = rmax; fDRStep = dr; }
48d874ff 86
87 protected:
8082a80b 88 TString fName; //!
89 TString fTitle; //!
90 std::vector<fastjet::PseudoJet> fInputVectors; //!
91 std::vector<fastjet::PseudoJet> fInclusiveJets; //!
92 std::vector<double> fSubtractedJetsPt; //!
93 std::vector<fastjet::PseudoJet> fConstituentSubtrJets; //!
94 fastjet::AreaDefinition *fAreaDef; //!
95 fastjet::VoronoiAreaSpec *fVorAreaSpec; //!
96 fastjet::GhostedAreaSpec *fGhostedAreaSpec; //!
97 fastjet::JetDefinition *fJetDef; //!
98 fastjet::JetDefinition::Plugin *fPlugin; //!
99 fastjet::RangeDefinition *fRange; //!
100 fastjet::ClusterSequenceArea *fClustSeq; //!
101 fastjet::Strategy fStrategy; //!
102 fastjet::JetAlgorithm fAlgor; //!
103 fastjet::RecombinationScheme fScheme; //!
104 fastjet::AreaType fAreaType; //!
105 Int_t fNGhostRepeats; //!
106 Double_t fGhostArea; //!
107 Double_t fMaxRap; //!
108 Double_t fR; //!
48d874ff 109 // no setters for the moment - used default values in the constructor
8082a80b 110 Double_t fGridScatter; //!
111 Double_t fKtScatter; //!
112 Double_t fMeanGhostKt; //!
113 Int_t fPluginAlgor; //!
48d874ff 114 // extra parameters
8082a80b 115 Double_t fMedUsedForBgSub; //!
116 Bool_t fUseArea4Vector; //!
5d0a5007 117#ifdef FASTJET_VERSION
726c9488 118 fastjet::JetMedianBackgroundEstimator *fBkrdEstimator; //!
93daf3f7 119 //from contrib package
726c9488 120 fastjet::contrib::GenericSubtractor *fGenSubtractor; //!
3ec5da8e 121 std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetMass; //!
122 std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoGRNum; //!
123 std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoGRDen; //!
124 std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetAngularity; //!
125 std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetpTD; //!
126 std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetCircularity; //!
df20822b 127 std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetSigma2; //!
3ec5da8e 128 std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetConstituent; //!
129 std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetLeSub; //!
5d0a5007 130#endif
8082a80b 131 Bool_t fLegacyMode; //!
132 Bool_t fUseExternalBkg; //!
133 Double_t fRho; // pT background density
134 Double_t fRhom; // mT background density
726c9488 135 Double_t fRMax; //!
136 Double_t fDRStep; //!
137 std::vector<double> fGRNumerator; //!
138 std::vector<double> fGRDenominator; //!
139 std::vector<double> fGRNumeratorSub; //!
140 std::vector<double> fGRDenominatorSub; //!
48d874ff 141
142 virtual void SubtractBackground(const Double_t median_pt = -1);
143
144 private:
145 AliFJWrapper();
146 AliFJWrapper(const AliFJWrapper& wrapper);
147 AliFJWrapper& operator = (const AliFJWrapper& wrapper);
148};
149#endif
150#endif
151
152#ifdef AliFJWrapper_CXX
153#undef AliFJWrapper_CXX
154
155#if defined __GNUC__
156#pragma GCC system_header
157#endif
158
159namespace fj = fastjet;
160
161//_________________________________________________________________________________________________
162AliFJWrapper::AliFJWrapper(const char *name, const char *title)
a9c8e539 163 :
164 fName (name)
165 , fTitle (title)
48d874ff 166 , fInputVectors ( )
167 , fInclusiveJets ( )
168 , fSubtractedJetsPt ( )
8082a80b 169 , fConstituentSubtrJets ( )
48d874ff 170 , fAreaDef (0)
171 , fVorAreaSpec (0)
172 , fGhostedAreaSpec (0)
173 , fJetDef (0)
174 , fPlugin (0)
175 , fRange (0)
176 , fClustSeq (0)
177 , fStrategy (fj::Best)
178 , fAlgor (fj::kt_algorithm)
179 , fScheme (fj::BIpt_scheme)
180 , fAreaType (fj::active_area)
181 , fNGhostRepeats (1)
6d75c7ee 182 , fGhostArea (0.005)
48d874ff 183 , fMaxRap (1.)
184 , fR (0.4)
185 , fGridScatter (1.0)
7dc91d51 186 , fKtScatter (0.1)
48d874ff 187 , fMeanGhostKt (1e-100)
188 , fPluginAlgor (0)
189 , fMedUsedForBgSub (0)
190 , fUseArea4Vector (kFALSE)
5d0a5007 191#ifdef FASTJET_VERSION
192 , fBkrdEstimator (0)
8082a80b 193 , fGenSubtractor (0)
194 , fGenSubtractorInfoJetMass ( )
fde4b180 195 , fGenSubtractorInfoGRNum ( )
196 , fGenSubtractorInfoGRDen ( )
0d13a63c 197 , fGenSubtractorInfoJetAngularity ( )
198 , fGenSubtractorInfoJetpTD ( )
199 , fGenSubtractorInfoJetCircularity( )
b9ccc456 200 , fGenSubtractorInfoJetSigma2()
0d13a63c 201 , fGenSubtractorInfoJetConstituent ( )
202 , fGenSubtractorInfoJetLeSub ( )
5d0a5007 203#endif
204 , fLegacyMode (false)
8082a80b 205 , fUseExternalBkg (false)
206 , fRho (0)
207 , fRhom (0)
726c9488 208 , fRMax(2.)
209 , fDRStep(0.04)
726c9488 210 , fGRNumerator()
211 , fGRDenominator()
212 , fGRNumeratorSub()
213 , fGRDenominatorSub()
48d874ff 214{
215 // Constructor.
216}
217
218//_________________________________________________________________________________________________
219AliFJWrapper::~AliFJWrapper()
220{
221 // Destructor.
222
223 delete fAreaDef;
224 delete fVorAreaSpec;
225 delete fGhostedAreaSpec;
226 delete fJetDef;
227 delete fPlugin;
228 delete fRange;
5d0a5007 229 delete fClustSeq;
230#ifdef FASTJET_VERSION
8082a80b 231 if (fBkrdEstimator) delete fBkrdEstimator;
232 if (fGenSubtractor) delete fGenSubtractor;
5d0a5007 233#endif
48d874ff 234}
235
236//_________________________________________________________________________________________________
237void AliFJWrapper::CopySettingsFrom(const AliFJWrapper& wrapper)
238{
239 // Copy some settings.
240 // You very often want to keep most of the settings
241 // but change only the algorithm or R - do it after call to this function
242
8082a80b 243 fStrategy = wrapper.fStrategy;
244 fAlgor = wrapper.fAlgor;
245 fScheme = wrapper.fScheme;
246 fAreaType = wrapper.fAreaType;
247 fNGhostRepeats = wrapper.fNGhostRepeats;
248 fGhostArea = wrapper.fGhostArea;
249 fMaxRap = wrapper.fMaxRap;
250 fR = wrapper.fR;
251 fGridScatter = wrapper.fGridScatter;
252 fKtScatter = wrapper.fKtScatter;
253 fMeanGhostKt = wrapper.fMeanGhostKt;
254 fPluginAlgor = wrapper.fPluginAlgor;
255 fUseArea4Vector = wrapper.fUseArea4Vector;
256 fLegacyMode = wrapper.fLegacyMode;
257 fUseExternalBkg = wrapper.fUseExternalBkg;
258 fRho = wrapper.fRho;
259 fRhom = wrapper.fRhom;
48d874ff 260}
261
262//_________________________________________________________________________________________________
a9c8e539 263void AliFJWrapper::Clear(const Option_t */*opt*/)
48d874ff 264{
265 // Simply clear the input vectors.
266 // Make sure done on every event if the instance is reused
267 // Reset the median to zero.
268
269 fInputVectors.clear();
270 fMedUsedForBgSub = 0;
271
272 // for the moment brute force delete everything
273 delete fAreaDef; fAreaDef = 0;
274 delete fVorAreaSpec; fVorAreaSpec = 0;
275 delete fGhostedAreaSpec; fGhostedAreaSpec = 0;
276 delete fJetDef; fJetDef = 0;
277 delete fPlugin; fPlugin = 0;
278 delete fRange; fRange = 0;
279 delete fClustSeq; fClustSeq = 0;
5d0a5007 280#ifdef FASTJET_VERSION
8082a80b 281 if (fBkrdEstimator) delete fBkrdEstimator ; fBkrdEstimator = 0;
282 if (fGenSubtractor) delete fGenSubtractor ; fGenSubtractor = 0;
5d0a5007 283#endif
48d874ff 284}
285
286//_________________________________________________________________________________________________
287void AliFJWrapper::AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index)
288{
289 // Make the input pseudojet.
290
291 fastjet::PseudoJet inVec(px, py, pz, E);
292
293 if (index > -99999) {
294 inVec.set_user_index(index);
295 } else {
296 inVec.set_user_index(fInputVectors.size());
297 }
298
299 // add to the fj container of input vectors
300 fInputVectors.push_back(inVec);
301}
302
303//_________________________________________________________________________________________________
304void AliFJWrapper::AddInputVector(const fj::PseudoJet& vec, Int_t index)
305{
306 // Add an input pseudojet.
307
308 fj::PseudoJet inVec = vec;
309
310 if (index > -99999) {
311 inVec.set_user_index(index);
312 } else {
313 inVec.set_user_index(fInputVectors.size());
314 }
315
316 // add to the fj container of input vectors
317 fInputVectors.push_back(inVec);
318}
319
320//_________________________________________________________________________________________________
321void AliFJWrapper::AddInputVectors(const std::vector<fj::PseudoJet>& vecs, Int_t offsetIndex)
322{
323 // Add the input from vector of pseudojets.
324
325 for (UInt_t i = 0; i < vecs.size(); ++i) {
326 fj::PseudoJet inVec = vecs[i];
48d874ff 327 if (offsetIndex > -99999)
328 inVec.set_user_index(fInputVectors.size() + offsetIndex);
329 // add to the fj container of input vectors
48d874ff 330 fInputVectors.push_back(inVec);
331 }
332}
333
334//_________________________________________________________________________________________________
335Double_t AliFJWrapper::GetJetArea(UInt_t idx) const
336{
337 // Get the jet area.
338
339 Double_t retval = -1; // really wrong area..
340 if ( idx < fInclusiveJets.size() ) {
7dc91d51 341 retval = fClustSeq->area(fInclusiveJets[idx]);
48d874ff 342 } else {
343 AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
344 }
345 return retval;
346}
347
4e22c11c 348//_________________________________________________________________________________________________
349fastjet::PseudoJet AliFJWrapper::GetJetAreaVector(UInt_t idx) const
350{
351 // Get the jet area as vector.
352 fastjet::PseudoJet retval;
353 if ( idx < fInclusiveJets.size() ) {
354 retval = fClustSeq->area_4vector(fInclusiveJets[idx]);
355 } else {
356 AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
357 }
358 return retval;
359}
360
48d874ff 361//_________________________________________________________________________________________________
362std::vector<double> AliFJWrapper::GetSubtractedJetsPts(Double_t median_pt, Bool_t sorted)
363{
364 // Get subtracted jets pTs, returns vector.
365
366 SubtractBackground(median_pt);
367
368 if (kTRUE == sorted) {
369 std::sort(fSubtractedJetsPt.begin(), fSubtractedJetsPt.begin());
370 }
371 return fSubtractedJetsPt;
372}
373
374//_________________________________________________________________________________________________
375Double_t AliFJWrapper::GetJetSubtractedPt(UInt_t idx) const
376{
377 // Get subtracted jets pTs, returns Double_t.
378
379 Double_t retval = -99999.; // really wrong pt..
380 if ( idx < fSubtractedJetsPt.size() ) {
381 retval = fSubtractedJetsPt[idx];
382 }
383 return retval;
384}
385
386//_________________________________________________________________________________________________
387std::vector<fastjet::PseudoJet>
388AliFJWrapper::GetJetConstituents(UInt_t idx) const
389{
390 // Get jets constituents.
391
392 std::vector<fastjet::PseudoJet> retval;
393
394 if ( idx < fInclusiveJets.size() ) {
395 retval = fClustSeq->constituents(fInclusiveJets[idx]);
396 } else {
397 AliError(Form("[e] ::GetJetConstituents wrong index: %d",idx));
398 }
399
400 return retval;
401}
402
403//_________________________________________________________________________________________________
404void AliFJWrapper::GetMedianAndSigma(Double_t &median, Double_t &sigma, Int_t remove) const
405{
406 // Get the median and sigma from fastjet.
407 // User can also do it on his own because the cluster sequence is exposed (via a getter)
408
409 if (!fClustSeq) {
410 AliError("[e] Run the jfinder first.");
b0a20945 411 return;
48d874ff 412 }
413
414 Double_t mean_area = 0;
415 try {
416 if(0 == remove) {
417 fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
418 } else {
419 std::vector<fastjet::PseudoJet> input_jets = sorted_by_pt(fClustSeq->inclusive_jets());
420 input_jets.erase(input_jets.begin(), input_jets.begin() + remove);
5d0a5007 421 fClustSeq->get_median_rho_and_sigma(input_jets, *fRange, fUseArea4Vector, median, sigma, mean_area);
48d874ff 422 input_jets.clear();
423 }
424 } catch (fj::Error) {
425 AliError(" [w] FJ Exception caught.");
426 median = -1.;
427 sigma = -1;
428 }
429}
430
431//_________________________________________________________________________________________________
432Int_t AliFJWrapper::Run()
433{
434 // Run the actual jet finder.
435
436 if (fAreaType == fj::voronoi_area) {
437 // Rfact - check dependence - default is 1.
438 // NOTE: hardcoded variable!
439 fVorAreaSpec = new fj::VoronoiAreaSpec(1.);
440 fAreaDef = new fj::AreaDefinition(*fVorAreaSpec);
441 } else {
442 fGhostedAreaSpec = new fj::GhostedAreaSpec(fMaxRap,
443 fNGhostRepeats,
444 fGhostArea,
445 fGridScatter,
446 fKtScatter,
447 fMeanGhostKt);
448
449 fAreaDef = new fj::AreaDefinition(*fGhostedAreaSpec, fAreaType);
450 }
451
452 // this is acceptable by fastjet:
453 fRange = new fj::RangeDefinition(fMaxRap - 0.95 * fR);
454
455 if (fAlgor == fj::plugin_algorithm) {
456 if (fPluginAlgor == 0) {
457 // SIS CONE ALGOR
458 // NOTE: hardcoded split parameter
459 Double_t overlap_threshold = 0.75; // NOTE: this actually splits a lot: thr/min(pt1,pt2)
460 fPlugin = new fj::SISConePlugin(fR,
461 overlap_threshold,
462 0, //search of stable cones - zero = until no more
463 1.0); // this should be seed effectively for proto jets
464 fJetDef = new fastjet::JetDefinition(fPlugin);
c515ed6d 465 } else if (fPluginAlgor == 1) {
466 // CDF cone
467 // NOTE: hardcoded split parameter
468 Double_t overlap_threshold = 0.75; // NOTE: this actually splits a lot: thr/min(pt1,pt2)
469 fPlugin = new fj::CDFMidPointPlugin(fR,
470 overlap_threshold,
471 1.0, //search of stable cones - zero = until no more
472 1.0); // this should be seed effectively for proto jets
473 fJetDef = new fastjet::JetDefinition(fPlugin);
48d874ff 474 } else {
475 AliError("[e] Unrecognized plugin number!");
476 }
477 } else {
478 fJetDef = new fj::JetDefinition(fAlgor, fR, fScheme, fStrategy);
479 }
480
481 try {
482 fClustSeq = new fj::ClusterSequenceArea(fInputVectors, *fJetDef, *fAreaDef);
483 } catch (fj::Error) {
484 AliError(" [w] FJ Exception caught.");
485 return -1;
486 }
487
5d0a5007 488 // FJ3 :: Define an JetMedianBackgroundEstimator just in case it will be used
489#ifdef FASTJET_VERSION
8082a80b 490 fBkrdEstimator = new fj::JetMedianBackgroundEstimator(fj::SelectorAbsRapMax(fMaxRap));
5d0a5007 491#endif
492
493 if (fLegacyMode) { SetLegacyFJ(); } // for FJ 2.x even if fLegacyMode is set, SetLegacyFJ is dummy
494
48d874ff 495 // inclusive jets:
496 fInclusiveJets.clear();
7dc91d51 497 fInclusiveJets = fClustSeq->inclusive_jets(0.0);
48d874ff 498
499 return 0;
500}
501
5d0a5007 502//_________________________________________________________________________________________________
503void AliFJWrapper::SetLegacyFJ()
504{
505 // This methods enable legacy behaviour (FastJet 2.x) when AliROOT is compiled with FastJet 3.x
506#ifdef FASTJET_VERSION
507 std::cout << "WARNING! Setting FastJet in legacy mode" << std::endl;
508 if (fGhostedAreaSpec) { fGhostedAreaSpec->set_fj2_placement(kTRUE); }
509 if (fBkrdEstimator) {
510 fBkrdEstimator->set_provide_fj2_sigma(kTRUE);
511 fBkrdEstimator->set_use_area_4vector(kFALSE);
512 }
513#endif
514}
515
48d874ff 516//_________________________________________________________________________________________________
517void AliFJWrapper::SubtractBackground(Double_t median_pt)
518{
519 // Subtract the background (specify the value - see below the meaning).
520 // Negative argument means the bg will be determined with the current algorithm
521 // this is the default behavior. Zero is allowed
522 // Note: user may set the switch for area4vector based subtraction.
523
524 Double_t median = 0;
525 Double_t sigma = 0;
526 Double_t mean_area = 0;
527
528 // clear the subtracted jet pt's vector<double>
529 fSubtractedJetsPt.clear();
5d0a5007 530
48d874ff 531 // check what was specified (default is -1)
532 if (median_pt < 0) {
533 try {
5d0a5007 534 fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
48d874ff 535 }
5d0a5007 536
48d874ff 537 catch (fj::Error) {
538 AliError(" [w] FJ Exception caught.");
539 median = -9999.;
540 sigma = -1;
541 fMedUsedForBgSub = median;
542 return;
543 }
544 fMedUsedForBgSub = median;
545 } else {
546 // we do not know the sigma in this case
547 sigma = -1;
548 if (0.0 == median_pt) {
549 // AliWarning(Form(" [w] Median specified for bg subtraction is ZERO: %f", median_pt )) << endl;
550 fMedUsedForBgSub = 0.;
551 } else {
552 fMedUsedForBgSub = median_pt;
553 }
554 }
555
556 // subtract:
557 for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
5d0a5007 558 if ( fUseArea4Vector ) {
48d874ff 559 // subtract the background using the area4vector
560 fj::PseudoJet area4v = fClustSeq->area_4vector(fInclusiveJets[i]);
561 fj::PseudoJet jet_sub = fInclusiveJets[i] - area4v * fMedUsedForBgSub;
562 fSubtractedJetsPt.push_back(jet_sub.perp()); // here we put only the pt of the jet - note: this can be negative
563 } else {
564 // subtract the background using scalars
565 // fj::PseudoJet jet_sub = fInclusiveJets[i] - area * fMedUsedForBgSub_;
566 Double_t area = fClustSeq->area(fInclusiveJets[i]);
567 // standard subtraction
568 Double_t pt_sub = fInclusiveJets[i].perp() - fMedUsedForBgSub * area;
569 fSubtractedJetsPt.push_back(pt_sub); // here we put only the pt of the jet - note: this can be negative
570 }
571 }
572}
573
8082a80b 574//_________________________________________________________________________________________________
575Int_t AliFJWrapper::DoGenericSubtractionJetMass() {
576 //Do generic subtraction for jet mass
93daf3f7 577#ifdef FASTJET_VERSION
8082a80b 578 if(fUseExternalBkg) fGenSubtractor = new fj::contrib::GenericSubtractor(fRho,fRhom);
579 else fGenSubtractor = new fj::contrib::GenericSubtractor(fBkrdEstimator);
580
581 // Define jet shape
582 AliJetShapeMass shapeMass;
583
584 // clear the generic subtractor info vector
585 fGenSubtractorInfoJetMass.clear();
586 for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
587 fj::contrib::GenericSubtractorInfo info;
588 if(fInclusiveJets[i].perp()>0.)
589 double subtracted_shape = (*fGenSubtractor)(shapeMass, fInclusiveJets[i], info);
590 fGenSubtractorInfoJetMass.push_back(info);
591 }
8082a80b 592#endif
593 return 0;
594}
595
726c9488 596//_________________________________________________________________________________________________
597Int_t AliFJWrapper::DoGenericSubtractionGR(Int_t ijet) {
598 //Do generic subtraction for jet mass
599#ifdef FASTJET_VERSION
600 if(fUseExternalBkg) fGenSubtractor = new fj::contrib::GenericSubtractor(fRho,fRhom);
601 else fGenSubtractor = new fj::contrib::GenericSubtractor(fBkrdEstimator);
602
603 if(ijet>fInclusiveJets.size()) return 0;
604
605 fGRNumerator.clear();
606 fGRDenominator.clear();
607 fGRNumeratorSub.clear();
608 fGRDenominatorSub.clear();
609
610 // Define jet shape
611 for(Double_t r = 0.; r<fRMax; r+=fDRStep) {
612 AliJetShapeGRNum shapeGRNum(r,fDRStep);
613 AliJetShapeGRDen shapeGRDen(r,fDRStep);
614
615 // clear the generic subtractor info vector
616 fGenSubtractorInfoGRNum.clear();
617 fGenSubtractorInfoGRDen.clear();
618 fj::contrib::GenericSubtractorInfo infoNum;
619 fj::contrib::GenericSubtractorInfo infoDen;
620 if(fInclusiveJets[ijet].perp()>0.) {
621 double sub_num = (*fGenSubtractor)(shapeGRNum, fInclusiveJets[ijet], infoNum);
622 double sub_den = (*fGenSubtractor)(shapeGRDen, fInclusiveJets[ijet], infoDen);
623 }
624 fGenSubtractorInfoGRNum.push_back(infoNum);
625 fGenSubtractorInfoGRDen.push_back(infoDen);
626 fGRNumerator.push_back(infoNum.unsubtracted());
627 fGRDenominator.push_back(infoDen.unsubtracted());
628 fGRNumeratorSub.push_back(infoNum.second_order_subtracted());
629 fGRDenominatorSub.push_back(infoDen.second_order_subtracted());
630 }
631#endif
632 return 0;
633}
0d13a63c 634//_________________________________________________________________________________________________
635Int_t AliFJWrapper::DoGenericSubtractionJetAngularity() {
636 //Do generic subtraction for jet mass
637#ifdef FASTJET_VERSION
638 if(fUseExternalBkg) fGenSubtractor = new fj::contrib::GenericSubtractor(fRho,fRhom);
639 else fGenSubtractor = new fj::contrib::GenericSubtractor(fBkrdEstimator);
640
641 // Define jet shape
642 AliJetShapeAngularity shapeAngularity;
643
644 // clear the generic subtractor info vector
645 fGenSubtractorInfoJetAngularity.clear();
646 for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
647 fj::contrib::GenericSubtractorInfo infoAng;
648 if(fInclusiveJets[i].perp()>0.)
649 double subtracted_shape = (*fGenSubtractor)(shapeAngularity, fInclusiveJets[i], infoAng);
650 fGenSubtractorInfoJetAngularity.push_back(infoAng);
651 }
0d13a63c 652#endif
653 return 0;
654}
655//_________________________________________________________________________________________________
656Int_t AliFJWrapper::DoGenericSubtractionJetpTD() {
657 //Do generic subtraction for jet mass
658#ifdef FASTJET_VERSION
659 if(fUseExternalBkg) fGenSubtractor = new fj::contrib::GenericSubtractor(fRho,fRhom);
660 else fGenSubtractor = new fj::contrib::GenericSubtractor(fBkrdEstimator);
661
662 // Define jet shape
663 AliJetShapepTD shapepTD;
664
665 // clear the generic subtractor info vector
666 fGenSubtractorInfoJetpTD.clear();
667 for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
668 fj::contrib::GenericSubtractorInfo infopTD;
669 if(fInclusiveJets[i].perp()>0.)
670 double subtracted_shape = (*fGenSubtractor)(shapepTD, fInclusiveJets[i], infopTD);
671 fGenSubtractorInfoJetpTD.push_back(infopTD);
672 }
0d13a63c 673#endif
674 return 0;
675}
676//_________________________________________________________________________________________________
677Int_t AliFJWrapper::DoGenericSubtractionJetCircularity() {
678 //Do generic subtraction for jet mass
679#ifdef FASTJET_VERSION
680 if(fUseExternalBkg) fGenSubtractor = new fj::contrib::GenericSubtractor(fRho,fRhom);
681 else fGenSubtractor = new fj::contrib::GenericSubtractor(fBkrdEstimator);
682
683 // Define jet shape
684 AliJetShapeCircularity shapecircularity;
685
686 // clear the generic subtractor info vector
687 fGenSubtractorInfoJetCircularity.clear();
688 for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
689 fj::contrib::GenericSubtractorInfo infoCirc;
690 if(fInclusiveJets[i].perp()>0.)
691 double subtracted_shape = (*fGenSubtractor)(shapecircularity, fInclusiveJets[i], infoCirc);
692 fGenSubtractorInfoJetCircularity.push_back(infoCirc);
693 }
0d13a63c 694#endif
695 return 0;
696}
b9ccc456 697//_________________________________________________________________________________________________
698Int_t AliFJWrapper::DoGenericSubtractionJetSigma2() {
699 //Do generic subtraction for jet mass
700#ifdef FASTJET_VERSION
701 if(fUseExternalBkg) fGenSubtractor = new fj::contrib::GenericSubtractor(fRho,fRhom);
702 else fGenSubtractor = new fj::contrib::GenericSubtractor(fBkrdEstimator);
703
704 // Define jet shape
705 AliJetShapeSigma2 shapesigma2;
0d13a63c 706
b9ccc456 707 // clear the generic subtractor info vector
708 fGenSubtractorInfoJetSigma2.clear();
709 for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
710 fj::contrib::GenericSubtractorInfo infoSigma;
711 if(fInclusiveJets[i].perp()>0.)
712 double subtracted_shape = (*fGenSubtractor)(shapesigma2, fInclusiveJets[i], infoSigma);
713 fGenSubtractorInfoJetSigma2.push_back(infoSigma);
714 }
715#endif
df20822b 716 return 0;
b9ccc456 717}
0d13a63c 718//_________________________________________________________________________________________________
719Int_t AliFJWrapper::DoGenericSubtractionJetConstituent() {
720 //Do generic subtraction for jet mass
721#ifdef FASTJET_VERSION
722 if(fUseExternalBkg) fGenSubtractor = new fj::contrib::GenericSubtractor(fRho,fRhom);
723 else fGenSubtractor = new fj::contrib::GenericSubtractor(fBkrdEstimator);
724
725 // Define jet shape
726 AliJetShapeConstituent shapeconst;
727
728 // clear the generic subtractor info vector
729 fGenSubtractorInfoJetConstituent.clear();
730 for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
731 fj::contrib::GenericSubtractorInfo infoConst;
732 if(fInclusiveJets[i].perp()>0.)
733 double subtracted_shape = (*fGenSubtractor)(shapeconst, fInclusiveJets[i], infoConst);
734 fGenSubtractorInfoJetConstituent.push_back(infoConst);
735 }
3ec5da8e 736#endif
0d13a63c 737 return 0;
738}
739
740//_________________________________________________________________________________________________
741Int_t AliFJWrapper::DoGenericSubtractionJetLeSub() {
742 //Do generic subtraction for jet mass
743#ifdef FASTJET_VERSION
744 if(fUseExternalBkg) fGenSubtractor = new fj::contrib::GenericSubtractor(fRho,fRhom);
745 else fGenSubtractor = new fj::contrib::GenericSubtractor(fBkrdEstimator);
746
747 // Define jet shape
748 AliJetShapeLeSub shapeLeSub;
749
750 // clear the generic subtractor info vector
751 fGenSubtractorInfoJetLeSub.clear();
752 for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
753 fj::contrib::GenericSubtractorInfo infoLeSub;
754 if(fInclusiveJets[i].perp()>0.)
755 double subtracted_shape = (*fGenSubtractor)(shapeLeSub, fInclusiveJets[i], infoLeSub);
756 fGenSubtractorInfoJetLeSub.push_back(infoLeSub);
757 }
0d13a63c 758#endif
759 return 0;
760}
726c9488 761
8082a80b 762//_________________________________________________________________________________________________
763Int_t AliFJWrapper::DoConstituentSubtraction() {
764 //Do constituent subtraction
93daf3f7 765#ifdef FASTJET_VERSION
8082a80b 766 fj::contrib::ConstituentSubtractor *subtractor;
8ed8c5bb 767 if(fUseExternalBkg)
8082a80b 768 subtractor = new fj::contrib::ConstituentSubtractor(fRho,fRhom,kFALSE,kTRUE);
8082a80b 769 else subtractor = new fj::contrib::ConstituentSubtractor(fBkrdEstimator);
770
771 //clear constituent subtracted jets
772 fConstituentSubtrJets.clear();
773 for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
774 fj::PseudoJet subtracted_jet(0.,0.,0.,0.);
775 if(fInclusiveJets[i].perp()>0.)
776 subtracted_jet = (*subtractor)(fInclusiveJets[i]);
777 fConstituentSubtrJets.push_back(subtracted_jet);
778 }
779 if(subtractor) delete subtractor;
780
781#endif
782 return 0;
783}
784
48d874ff 785//_________________________________________________________________________________________________
786void AliFJWrapper::SetupAlgorithmfromOpt(const char *option)
787{
788 // Setup algorithm from char.
789
790 std::string opt(option);
791
792 if (!opt.compare("kt")) fAlgor = fj::kt_algorithm;
793 if (!opt.compare("antikt")) fAlgor = fj::antikt_algorithm;
794 if (!opt.compare("cambridge")) fAlgor = fj::cambridge_algorithm;
795 if (!opt.compare("genkt")) fAlgor = fj::genkt_algorithm;
796 if (!opt.compare("cambridge_passive")) fAlgor = fj::cambridge_for_passive_algorithm;
797 if (!opt.compare("genkt_passive")) fAlgor = fj::genkt_for_passive_algorithm;
798 if (!opt.compare("ee_kt")) fAlgor = fj::ee_kt_algorithm;
799 if (!opt.compare("ee_genkt")) fAlgor = fj::ee_genkt_algorithm;
800 if (!opt.compare("plugin")) fAlgor = fj::plugin_algorithm;
801}
802
803//_________________________________________________________________________________________________
804void AliFJWrapper::SetupAreaTypefromOpt(const char *option)
805{
806 // Setup area type from char.
807
808 std::string opt(option);
809
810 if (!opt.compare("active")) fAreaType = fj::active_area;
811 if (!opt.compare("invalid")) fAreaType = fj::invalid_area;
812 if (!opt.compare("active_area_explicit_ghosts")) fAreaType = fj::active_area_explicit_ghosts;
813 if (!opt.compare("one_ghost_passive")) fAreaType = fj::one_ghost_passive_area;
814 if (!opt.compare("passive")) fAreaType = fj::passive_area;
815 if (!opt.compare("voronoi")) fAreaType = fj::voronoi_area;
816}
817
818//_________________________________________________________________________________________________
819void AliFJWrapper::SetupSchemefromOpt(const char *option)
820{
821 //
822 // setup scheme from char
823 //
824
825 std::string opt(option);
826
827 if (!opt.compare("BIpt")) fScheme = fj::BIpt_scheme;
828 if (!opt.compare("BIpt2")) fScheme = fj::BIpt2_scheme;
829 if (!opt.compare("E")) fScheme = fj::E_scheme;
830 if (!opt.compare("pt")) fScheme = fj::pt_scheme;
831 if (!opt.compare("pt2")) fScheme = fj::pt2_scheme;
832 if (!opt.compare("Et")) fScheme = fj::Et_scheme;
833 if (!opt.compare("Et2")) fScheme = fj::Et2_scheme;
834}
835
836//_________________________________________________________________________________________________
837void AliFJWrapper::SetupStrategyfromOpt(const char *option)
838{
839 // Setup strategy from char.
840
841 std::string opt(option);
842
843 if (!opt.compare("Best")) fStrategy = fj::Best;
844 if (!opt.compare("N2MinHeapTiled")) fStrategy = fj::N2MinHeapTiled;
845 if (!opt.compare("N2Tiled")) fStrategy = fj::N2Tiled;
846 if (!opt.compare("N2PoorTiled")) fStrategy = fj::N2PoorTiled;
847 if (!opt.compare("N2Plain")) fStrategy = fj::N2Plain;
848 if (!opt.compare("N3Dumb")) fStrategy = fj::N3Dumb;
849 if (!opt.compare("NlnN")) fStrategy = fj::NlnN;
850 if (!opt.compare("NlnN3pi")) fStrategy = fj::NlnN3pi;
851 if (!opt.compare("NlnN4pi")) fStrategy = fj::NlnN4pi;
852 if (!opt.compare("NlnNCam4pi")) fStrategy = fj::NlnNCam4pi;
853 if (!opt.compare("NlnNCam2pi2R")) fStrategy = fj::NlnNCam2pi2R;
854 if (!opt.compare("NlnNCam")) fStrategy = fj::NlnNCam;
855 if (!opt.compare("plugin")) fStrategy = fj::plugin_strategy;
856}
857#endif