AliPhysics  deb3cd0 (deb3cd0)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliFJWrapper.h
Go to the documentation of this file.
1 #ifndef AliFJWrapper_H
2 #define AliFJWrapper_H
3 
4 #if !defined(__CINT__)
5 
6 #include <vector>
7 #include <TString.h>
8 #include "AliLog.h"
9 #include "FJ_includes.h"
10 #include "AliJetShape.h"
11 
12 
13 class AliFJWrapper
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 AddInputGhost (Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index = -99999);
23  virtual const char *ClassName() const { return "AliFJWrapper"; }
24  virtual void Clear(const Option_t* /*opt*/ = "");
25  virtual void ClearMemory();
26  virtual void CopySettingsFrom (const AliFJWrapper& wrapper);
27  virtual void GetMedianAndSigma(Double_t& median, Double_t& sigma, Int_t remove = 0) const;
28  fastjet::ClusterSequenceArea* GetClusterSequence() const { return fClustSeq; }
29  fastjet::ClusterSequence* GetClusterSequenceSA() const { return fClustSeqSA; }
30  fastjet::ClusterSequenceActiveAreaExplicitGhosts* GetClusterSequenceGhosts() const { return fClustSeqActGhosts; }
31  const std::vector<fastjet::PseudoJet>& GetInputVectors() const { return fInputVectors; }
32  const std::vector<fastjet::PseudoJet>& GetEventSubInputVectors() const { return fEventSubInputVectors; }
33  const std::vector<fastjet::PseudoJet>& GetInputGhosts() const { return fInputGhosts; }
34  const std::vector<fastjet::PseudoJet>& GetInclusiveJets() const { return fInclusiveJets; }
35  const std::vector<fastjet::PseudoJet>& GetEventSubJets() const { return fEventSubJets; }
36  const std::vector<fastjet::PseudoJet>& GetFilteredJets() const { return fFilteredJets; }
37  std::vector<fastjet::PseudoJet> GetJetConstituents(UInt_t idx) const;
38  std::vector<fastjet::PseudoJet> GetEventSubJetConstituents(UInt_t idx) const;
39  std::vector<fastjet::PseudoJet> GetFilteredJetConstituents(UInt_t idx) const;
40  Double_t GetMedianUsedForBgSubtraction() const { return fMedUsedForBgSub; }
41  const char* GetName() const { return fName; }
42  const char* GetTitle() const { return fTitle; }
43  Double_t GetJetArea (UInt_t idx) const;
44  Double_t GetEventSubJetArea (UInt_t idx) const;
45  fastjet::PseudoJet GetJetAreaVector (UInt_t idx) const;
46  fastjet::PseudoJet GetEventSubJetAreaVector (UInt_t idx) const;
47  Double_t GetFilteredJetArea (UInt_t idx) const;
48  fastjet::PseudoJet GetFilteredJetAreaVector(UInt_t idx) const;
49  Double_t GetJetSubtractedPt (UInt_t idx) const;
50  virtual std::vector<double> GetSubtractedJetsPts(Double_t median_pt = -1, Bool_t sorted = kFALSE);
51  Bool_t GetLegacyMode() { return fLegacyMode; }
52  Bool_t GetDoFilterArea() { return fDoFilterArea; }
53  Double_t NSubjettiness(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);
54  Double32_t NSubjettinessDerivativeSub(Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Double_t JetR, fastjet::PseudoJet jet, Int_t Option=0, Int_t Measure=0);
55 #ifdef FASTJET_VERSION
56  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetMass() const {return fGenSubtractorInfoJetMass ; }
57  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetAngularity() const {return fGenSubtractorInfoJetAngularity ; }
58  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetpTD() const {return fGenSubtractorInfoJetpTD ; }
59  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetCircularity() const {return fGenSubtractorInfoJetCircularity ; }
60  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetSigma2() const {return fGenSubtractorInfoJetSigma2 ; }
61  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetConstituent() const {return fGenSubtractorInfoJetConstituent ; }
62  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetLeSub() const {return fGenSubtractorInfoJetLeSub ; }
63  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJet1subjettiness_kt() const {return fGenSubtractorInfoJet1subjettiness_kt ; }
64  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJet2subjettiness_kt() const {return fGenSubtractorInfoJet2subjettiness_kt ; }
65  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJet3subjettiness_kt() const {return fGenSubtractorInfoJet3subjettiness_kt ; }
66  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetOpeningAngle_kt() const {return fGenSubtractorInfoJetOpeningAngle_kt ; }
67  const std::vector<fastjet::PseudoJet> GetConstituentSubtrJets() const {return fConstituentSubtrJets ; }
68  const std::vector<fastjet::PseudoJet> GetGroomedJets() const {return fGroomedJets ; }
69  Int_t CreateGenSub(); // fastjet::contrib::GenericSubtractor
70  Int_t CreateConstituentSub(); // fastjet::contrib::ConstituentSubtractor
71  Int_t CreateEventConstituentSub(); //fastjet::contrib::ConstituentSubtractor
72  Int_t CreateSoftDrop();
73 #endif
74  virtual std::vector<double> GetGRNumerator() const { return fGRNumerator ; }
75  virtual std::vector<double> GetGRDenominator() const { return fGRDenominator ; }
76  virtual std::vector<double> GetGRNumeratorSub() const { return fGRNumeratorSub ; }
77  virtual std::vector<double> GetGRDenominatorSub() const { return fGRDenominatorSub ; }
78 
79  virtual void RemoveLastInputVector();
80 
81  virtual Int_t Run();
82  virtual Int_t Filter();
83  virtual Int_t DoGenericSubtractionJetMass();
84  virtual Int_t DoGenericSubtractionGR(Int_t ijet);
85  virtual Int_t DoGenericSubtractionJetAngularity();
86  virtual Int_t DoGenericSubtractionJetpTD();
87  virtual Int_t DoGenericSubtractionJetCircularity();
88  virtual Int_t DoGenericSubtractionJetSigma2();
89  virtual Int_t DoGenericSubtractionJetConstituent();
90  virtual Int_t DoGenericSubtractionJetLeSub();
91  virtual Int_t DoGenericSubtractionJet1subjettiness_kt();
92  virtual Int_t DoGenericSubtractionJet2subjettiness_kt();
93  virtual Int_t DoGenericSubtractionJet3subjettiness_kt();
94  virtual Int_t DoGenericSubtractionJetOpeningAngle_kt();
95  virtual Int_t DoConstituentSubtraction();
96  virtual Int_t DoEventConstituentSubtraction();
97  virtual Int_t DoSoftDrop();
98 
99  void SetName(const char* name) { fName = name; }
100  void SetTitle(const char* title) { fTitle = title; }
101  void SetStrategy(const fastjet::Strategy &strat) { fStrategy = strat; }
102  void SetAlgorithm(const fastjet::JetAlgorithm &algor) { fAlgor = algor; }
103  void SetRecombScheme(const fastjet::RecombinationScheme &scheme) { fScheme = scheme; }
104  void SetAreaType(const fastjet::AreaType &atype) { fAreaType = atype; }
105  void SetNRepeats(Int_t nrepeat) { fNGhostRepeats = nrepeat; }
106  void SetGhostArea(Double_t gharea) { fGhostArea = gharea; }
107  void SetMaxRap(Double_t maxrap) { fMaxRap = maxrap; }
108  void SetR(Double_t r) { fR = r; }
109  void SetGridScatter(Double_t gridSc) { fGridScatter = gridSc; }
110  void SetKtScatter(Double_t ktSc) { fKtScatter = ktSc; }
111  void SetMeanGhostKt(Double_t meankt) { fMeanGhostKt = meankt; }
112  void SetPluginAlgor(Int_t plugin) { fPluginAlgor = plugin; }
113  void SetUseArea4Vector(Bool_t useA4v) { fUseArea4Vector = useA4v; }
114  void SetupAlgorithmfromOpt(const char *option);
115  void SetupAreaTypefromOpt(const char *option);
116  void SetupSchemefromOpt(const char *option);
117  void SetupStrategyfromOpt(const char *option);
118  void SetLegacyMode (Bool_t mode) { fLegacyMode ^= mode; }
119  void SetLegacyFJ();
120  void SetUseExternalBkg(Bool_t b, Double_t rho, Double_t rhom) { fUseExternalBkg = b; fRho = rho; fRhom = rhom;}
121  void SetRMaxAndStep(Double_t rmax, Double_t dr) {fRMax = rmax; fDRStep = dr; }
122  void SetRhoRhom (Double_t rho, Double_t rhom) { fUseExternalBkg = kTRUE; fRho = rho; fRhom = rhom;} // if using rho,rhom then fUseExternalBkg is true
123  void SetMinJetPt(Double_t MinPt) {fMinJetPt=MinPt;}
124  void SetEventSub(Bool_t b) {fEventSub = b;}
125  void SetMaxDelR(Double_t r) {fUseMaxDelR = kTRUE; fMaxDelR = r;}
126 
127  protected:
128  TString fName;
129  TString fTitle;
130  std::vector<fastjet::PseudoJet> fInputVectors;
131  std::vector<fastjet::PseudoJet> fEventSubInputVectors;
132  std::vector<fastjet::PseudoJet> fEventSubCorrectedVectors;
133  std::vector<fastjet::PseudoJet> fInputGhosts;
134  std::vector<fastjet::PseudoJet> fInclusiveJets;
135  std::vector<fastjet::PseudoJet> fEventSubJets;
136  std::vector<fastjet::PseudoJet> fFilteredJets;
137  std::vector<double> fSubtractedJetsPt;
138  std::vector<fastjet::PseudoJet> fConstituentSubtrJets;
139  std::vector<fastjet::PseudoJet> fGroomedJets;
140  fastjet::AreaDefinition *fAreaDef;
141  fastjet::VoronoiAreaSpec *fVorAreaSpec;
142  fastjet::GhostedAreaSpec *fGhostedAreaSpec;
143  fastjet::JetDefinition *fJetDef;
144  fastjet::JetDefinition::Plugin *fPlugin;
145 #ifndef FASTJET_VERSION
146  fastjet::RangeDefinition *fRange;
147 #else
148  fastjet::Selector *fRange;
149 #endif
150  fastjet::ClusterSequenceArea *fClustSeq;
151  fastjet::ClusterSequenceArea *fClustSeqES;
152  fastjet::ClusterSequence *fClustSeqSA;
153  fastjet::ClusterSequenceActiveAreaExplicitGhosts *fClustSeqActGhosts;
154  fastjet::Strategy fStrategy;
155  fastjet::JetAlgorithm fAlgor;
156  fastjet::RecombinationScheme fScheme;
157  fastjet::AreaType fAreaType;
158  Int_t fNGhostRepeats;
159  Double_t fGhostArea;
160  Double_t fMaxRap;
161  Double_t fR;
162  Double_t fMinJetPt;
163  // no setters for the moment - used default values in the constructor
164  Double_t fGridScatter;
165  Double_t fKtScatter;
166  Double_t fMeanGhostKt;
167  Int_t fPluginAlgor;
168  // extra parameters
169  Double_t fMedUsedForBgSub;
170  Bool_t fUseArea4Vector;
171  // condition to stop the grooming (rejection of soft splitting) z > fZcut theta^fBeta
172  Double_t fZcut; // fZcut = 0.1
173  Double_t fBeta; // fBeta = 0
174  Bool_t fEventSub;
175  Bool_t fUseMaxDelR;
176  Double_t fMaxDelR;
177 #ifdef FASTJET_VERSION
178  fastjet::JetMedianBackgroundEstimator *fBkrdEstimator;
179  //from contrib package
180  fastjet::contrib::GenericSubtractor *fGenSubtractor;
181  fastjet::contrib::ConstituentSubtractor *fConstituentSubtractor;
182  fastjet::contrib::ConstituentSubtractor *fEventConstituentSubtractor;
183  fastjet::contrib::SoftDrop *fSoftDrop;
184  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetMass;
185  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoGRNum;
186  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoGRDen;
187  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetAngularity;
188  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetpTD;
189  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetCircularity;
190  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetSigma2;
191  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetConstituent;
192  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetLeSub;
193  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet1subjettiness_kt;
194  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet2subjettiness_kt;
195  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet3subjettiness_kt;
196  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetOpeningAngle_kt;
197 #endif
198  Bool_t fDoFilterArea;
199  Bool_t fLegacyMode;
200  Bool_t fUseExternalBkg;
201  Double_t fRho; // pT background density
202  Double_t fRhom; // mT background density
203  Double_t fRMax;
204  Double_t fDRStep;
205  std::vector<double> fGRNumerator;
206  std::vector<double> fGRDenominator;
207  std::vector<double> fGRNumeratorSub;
208  std::vector<double> fGRDenominatorSub;
209 
210  virtual void SubtractBackground(const Double_t median_pt = -1);
211 
212  private:
213  AliFJWrapper();
214  AliFJWrapper(const AliFJWrapper& wrapper);
215  AliFJWrapper& operator = (const AliFJWrapper& wrapper);
216 };
217 #endif
218 #endif
219 
220 #ifdef AliFJWrapper_CXX
221 #undef AliFJWrapper_CXX
222 
223 #if defined __GNUC__
224 #pragma GCC system_header
225 #endif
226 
227 namespace fj = fastjet;
228 
229 //_________________________________________________________________________________________________
230 AliFJWrapper::AliFJWrapper(const char *name, const char *title)
231  :
232  fName (name)
233  , fTitle (title)
234  , fInputVectors ( )
235  , fEventSubInputVectors ( )
236  , fEventSubCorrectedVectors ( )
237  , fInputGhosts ( )
238  , fInclusiveJets ( )
239  , fEventSubJets ( )
240  , fFilteredJets ( )
241  , fSubtractedJetsPt ( )
242  , fConstituentSubtrJets ( )
243  , fSoftDrop ( )
244  , fAreaDef (0)
245  , fVorAreaSpec (0)
246  , fGhostedAreaSpec (0)
247  , fJetDef (0)
248  , fPlugin (0)
249  , fRange (0)
250  , fClustSeq (0)
251  , fClustSeqES (0)
252  , fClustSeqSA (0)
253  , fClustSeqActGhosts (0)
254  , fStrategy (fj::Best)
255  , fAlgor (fj::kt_algorithm)
256  , fScheme (fj::BIpt_scheme)
257  , fAreaType (fj::active_area)
258  , fNGhostRepeats (1)
259  , fGhostArea (0.005)
260  , fMaxRap (1.)
261  , fR (0.4)
262  , fGridScatter (1.0)
263  , fKtScatter (0.1)
264  , fMeanGhostKt (1e-100)
265  , fPluginAlgor (0)
266  , fMedUsedForBgSub (0)
267  , fUseArea4Vector (kFALSE)
268  , fZcut(0.1)
269  , fBeta(0)
270  , fEventSub (kFALSE)
271  , fUseMaxDelR (kFALSE)
272  , fMaxDelR (0.4)
273 #ifdef FASTJET_VERSION
274  , fBkrdEstimator (0)
275  , fGenSubtractor (0)
276  , fConstituentSubtractor (0)
277  , fEventConstituentSubtractor (0)
278  , fGenSubtractorInfoJetMass ( )
279  , fGenSubtractorInfoGRNum ( )
280  , fGenSubtractorInfoGRDen ( )
281  , fGenSubtractorInfoJetAngularity ( )
282  , fGenSubtractorInfoJetpTD ( )
283  , fGenSubtractorInfoJetCircularity( )
284  , fGenSubtractorInfoJetSigma2()
285  , fGenSubtractorInfoJetConstituent ( )
286  , fGenSubtractorInfoJetLeSub ( )
287  , fGenSubtractorInfoJet1subjettiness_kt ( )
288  , fGenSubtractorInfoJet2subjettiness_kt ( )
289  , fGenSubtractorInfoJet3subjettiness_kt ( )
290  , fGenSubtractorInfoJetOpeningAngle_kt ( )
291 #endif
292  , fDoFilterArea (false)
293  , fLegacyMode (false)
294  , fUseExternalBkg (false)
295  , fRho (0)
296  , fRhom (0)
297  , fRMax(2.)
298  , fDRStep(0.04)
299  , fGRNumerator()
300  , fGRDenominator()
301  , fGRNumeratorSub()
302  , fGRDenominatorSub()
303 {
304  // Constructor.
305 }
306 
307 //_________________________________________________________________________________________________
308 AliFJWrapper::~AliFJWrapper()
309 {
310  // Destructor.
311  ClearMemory();
312 }
313 
314 //_________________________________________________________________________________________________
315 void AliFJWrapper::ClearMemory()
316 {
317  // Destructor.
318  if (fAreaDef) { delete fAreaDef; fAreaDef = NULL; }
319  if (fVorAreaSpec) { delete fVorAreaSpec; fVorAreaSpec = NULL; }
320  if (fGhostedAreaSpec) { delete fGhostedAreaSpec; fGhostedAreaSpec = NULL; }
321  if (fJetDef) { delete fJetDef; fJetDef = NULL; }
322  if (fPlugin) { delete fPlugin; fPlugin = NULL; }
323  if (fRange) { delete fRange; fRange = NULL; }
324  if (fClustSeq) { delete fClustSeq; fClustSeq = NULL; }
325  if (fClustSeqES) { delete fClustSeqES; fClustSeqES = NULL; }
326  if (fClustSeqSA) { delete fClustSeqSA; fClustSeqSA = NULL; }
327  if (fClustSeqActGhosts) { delete fClustSeqActGhosts; fClustSeqActGhosts = NULL; }
328  #ifdef FASTJET_VERSION
329  if (fBkrdEstimator) { delete fBkrdEstimator; fBkrdEstimator = NULL; }
330  if (fGenSubtractor) { delete fGenSubtractor; fGenSubtractor = NULL; }
331  if (fConstituentSubtractor) { delete fConstituentSubtractor; fConstituentSubtractor = NULL; }
332  if (fSoftDrop) { delete fSoftDrop; fSoftDrop = NULL;}
333  #endif
334 }
335 
336 //_________________________________________________________________________________________________
337 void AliFJWrapper::CopySettingsFrom(const AliFJWrapper& wrapper)
338 {
339  // Copy some settings.
340  // You very often want to keep most of the settings
341  // but change only the algorithm or R - do it after call to this function
342 
343  fStrategy = wrapper.fStrategy;
344  fAlgor = wrapper.fAlgor;
345  fScheme = wrapper.fScheme;
346  fAreaType = wrapper.fAreaType;
347  fNGhostRepeats = wrapper.fNGhostRepeats;
348  fGhostArea = wrapper.fGhostArea;
349  fMaxRap = wrapper.fMaxRap;
350  fR = wrapper.fR;
351  fGridScatter = wrapper.fGridScatter;
352  fKtScatter = wrapper.fKtScatter;
353  fMeanGhostKt = wrapper.fMeanGhostKt;
354  fPluginAlgor = wrapper.fPluginAlgor;
355  fUseArea4Vector = wrapper.fUseArea4Vector;
356  fZcut = wrapper.fZcut;
357  fBeta = wrapper.fBeta;
358  fLegacyMode = wrapper.fLegacyMode;
359  fUseExternalBkg = wrapper.fUseExternalBkg;
360  fRho = wrapper.fRho;
361  fRhom = wrapper.fRhom;
362 }
363 
364 //_________________________________________________________________________________________________
365 void AliFJWrapper::Clear(const Option_t */*opt*/)
366 {
367  // Simply clear the input vectors.
368  // Make sure done on every event if the instance is reused
369  // Reset the median to zero.
370 
371  fInputVectors.clear();
372  fEventSubInputVectors.clear();
373  fInputGhosts.clear();
374  fMedUsedForBgSub = 0;
375 
376  // for the moment brute force delete everything
377  ClearMemory();
378 }
379 
380 //_________________________________________________________________________________________________
381 void AliFJWrapper::RemoveLastInputVector()
382 {
383  // Remove last input vector
384 
385  fInputVectors.pop_back();
386  fEventSubInputVectors.pop_back();
387 }
388 
389 //_________________________________________________________________________________________________
390 void AliFJWrapper::AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index)
391 {
392  // Make the input pseudojet.
393 
394  fastjet::PseudoJet inVec(px, py, pz, E);
395 
396  // Salvatore Aiola: not sure why this was done...
397  //if (index > -99999) {
398  inVec.set_user_index(index);
399  //} else {
400  //inVec.set_user_index(fInputVectors.size());
401  //}
402 
403  // add to the fj container of input vectors
404  fInputVectors.push_back(inVec);
405  if(fEventSub) fEventSubInputVectors.push_back(inVec);
406 
407 }
408 
409 //_________________________________________________________________________________________________
410 void AliFJWrapper::AddInputVector(const fj::PseudoJet& vec, Int_t index)
411 {
412  // Add an input pseudojet.
413 
414  fj::PseudoJet inVec = vec;
415 
416  // Salvatore Aiola: not sure why this was done...
418  inVec.set_user_index(index);
419  //} else {
420  //inVec.set_user_index(fInputVectors.size());
421  //}
422 
423  // add to the fj container of input vectors
424  fInputVectors.push_back(inVec);
425  //if(fEventSub) fEventSubInputVectors.push_back(inVec);
426 }
427 
428 //_________________________________________________________________________________________________
429 void AliFJWrapper::AddInputVectors(const std::vector<fj::PseudoJet>& vecs, Int_t offsetIndex)
430 {
431  // Add the input from vector of pseudojets.
432 
433  for (UInt_t i = 0; i < vecs.size(); ++i) {
434  fj::PseudoJet inVec = vecs[i];
435  if (offsetIndex > -99999)
436  inVec.set_user_index(fInputVectors.size() + offsetIndex);
437  // add to the fj container of input vectors
438  fInputVectors.push_back(inVec);
439  }
440  //Is this correct?
441  /* if(fEventSub){
442  for (UInt_t i = 0; i < vecs.size(); ++i) {
443  fj::PseudoJet inVec = vecs[i];
444  if (offsetIndex > -99999)
445  inVec.set_user_index(fEventSubInputVectors.size() + offsetIndex);
446  // add to the fj container of input vectors
447  fEventSubInputVectors.push_back(inVec);
448  }
449  }*/
450 }
451 
452 //_________________________________________________________________________________________________
453 void AliFJWrapper::AddInputGhost(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index)
454 {
455  // Make the input pseudojet.
456 
457  fastjet::PseudoJet inVec(px, py, pz, E);
458 
459  if (index > -99999) {
460  inVec.set_user_index(index);
461  } else {
462  inVec.set_user_index(fInputGhosts.size());
463  }
464 
465  // add to the fj container of input vectors
466  fInputGhosts.push_back(inVec);
467  if (!fDoFilterArea) fDoFilterArea = kTRUE;
468 }
469 
470 //_________________________________________________________________________________________________
471 Double_t AliFJWrapper::GetJetArea(UInt_t idx) const
472 {
473  // Get the jet area.
474 
475  Double_t retval = -1; // really wrong area..
476  if ( idx < fInclusiveJets.size() ) {
477  retval = fClustSeq->area(fInclusiveJets[idx]);
478  } else {
479  AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
480  }
481  return retval;
482 }
483 
484 //_________________________________________________________________________________________________
485 Double_t AliFJWrapper::GetEventSubJetArea(UInt_t idx) const
486 {
487  // Get the jet area.
488 
489  Double_t retval = -1; // really wrong area..
490  if ( idx < fEventSubJets.size() ) {
491  retval = fClustSeqES->area(fEventSubJets[idx]);
492  } else {
493  AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
494  }
495  return retval;
496 }
497 //_________________________________________________________________________________________________
498 Double_t AliFJWrapper::GetFilteredJetArea(UInt_t idx) const
499 {
500  // Get the filtered jet area.
501 
502  Double_t retval = -1; // really wrong area..
503  if (fDoFilterArea && fClustSeqActGhosts && (idx<fFilteredJets.size())) {
504  retval = fClustSeqActGhosts->area(fFilteredJets[idx]);
505  } else {
506  AliError(Form("[e] ::GetFilteredJetArea wrong index: %d",idx));
507  }
508  return retval;
509 }
510 
511 //_________________________________________________________________________________________________
512 fastjet::PseudoJet AliFJWrapper::GetJetAreaVector(UInt_t idx) const
513 {
514  // Get the jet area as vector.
515  fastjet::PseudoJet retval;
516  if ( idx < fInclusiveJets.size() ) {
517  retval = fClustSeq->area_4vector(fInclusiveJets[idx]);
518  } else {
519  AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
520  }
521  return retval;
522 }
523 
524 //_________________________________________________________________________________________________
525 fastjet::PseudoJet AliFJWrapper::GetEventSubJetAreaVector(UInt_t idx) const
526 {
527  // Get the jet area as vector.
528  fastjet::PseudoJet retval;
529  if ( idx < fEventSubJets.size() ) {
530  retval = fClustSeqES->area_4vector(fEventSubJets[idx]);
531  } else {
532  AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
533  }
534  return retval;
535 }
536 
537 //_________________________________________________________________________________________________
538 fastjet::PseudoJet AliFJWrapper::GetFilteredJetAreaVector(UInt_t idx) const
539 {
540  // Get the jet area as vector.
541  fastjet::PseudoJet retval;
542  if (fDoFilterArea && fClustSeqActGhosts && (idx<fFilteredJets.size())) {
543  retval = fClustSeqActGhosts->area_4vector(fFilteredJets[idx]);
544  } else {
545  AliError(Form("[e] ::GetFilteredJetArea wrong index: %d",idx));
546  }
547  return retval;
548 }
549 
550 //_________________________________________________________________________________________________
551 std::vector<double> AliFJWrapper::GetSubtractedJetsPts(Double_t median_pt, Bool_t sorted)
552 {
553  // Get subtracted jets pTs, returns vector.
554 
555  SubtractBackground(median_pt);
556 
557  if (kTRUE == sorted) {
558  std::sort(fSubtractedJetsPt.begin(), fSubtractedJetsPt.begin());
559  }
560  return fSubtractedJetsPt;
561 }
562 
563 //_________________________________________________________________________________________________
564 Double_t AliFJWrapper::GetJetSubtractedPt(UInt_t idx) const
565 {
566  // Get subtracted jets pTs, returns Double_t.
567 
568  Double_t retval = -99999.; // really wrong pt..
569  if ( idx < fSubtractedJetsPt.size() ) {
570  retval = fSubtractedJetsPt[idx];
571  }
572  return retval;
573 }
574 
575 //_________________________________________________________________________________________________
576 std::vector<fastjet::PseudoJet>
577 AliFJWrapper::GetJetConstituents(UInt_t idx) const
578 {
579  // Get jets constituents.
580 
581  std::vector<fastjet::PseudoJet> retval;
582 
583  if ( idx < fInclusiveJets.size() ) {
584  retval = fClustSeq->constituents(fInclusiveJets[idx]);
585  } else {
586  AliError(Form("[e] ::GetJetConstituents wrong index: %d",idx));
587  }
588 
589  return retval;
590 }
591 
592 //_________________________________________________________________________________________________
593 std::vector<fastjet::PseudoJet>
594 AliFJWrapper::GetEventSubJetConstituents(UInt_t idx) const
595 {
596  // Get jets constituents.
597 
598  std::vector<fastjet::PseudoJet> retval;
599 
600  if ( idx < fEventSubJets.size() ) {
601  retval = fClustSeqES->constituents(fEventSubJets[idx]);
602  } else {
603  AliError(Form("[e] ::GetJetConstituents wrong index: %d",idx));
604  }
605 
606  return retval;
607 }
608 
609 //_________________________________________________________________________________________________
610 std::vector<fastjet::PseudoJet>
611 AliFJWrapper::GetFilteredJetConstituents(UInt_t idx) const
612 {
613  // Get jets constituents.
614 
615  std::vector<fastjet::PseudoJet> retval;
616 
617  if ( idx < fFilteredJets.size() ) {
618  if (fClustSeqSA) retval = fClustSeqSA->constituents(fFilteredJets[idx]);
619  if (fClustSeqActGhosts) retval = fClustSeqActGhosts->constituents(fFilteredJets[idx]);
620  } else {
621  AliError(Form("[e] ::GetFilteredJetConstituents wrong index: %d",idx));
622  }
623 
624  return retval;
625 }
626 
627 //_________________________________________________________________________________________________
628 void AliFJWrapper::GetMedianAndSigma(Double_t &median, Double_t &sigma, Int_t remove) const
629 {
630  // Get the median and sigma from fastjet.
631  // User can also do it on his own because the cluster sequence is exposed (via a getter)
632 
633  if (!fClustSeq) {
634  AliError("[e] Run the jfinder first.");
635  return;
636  }
637 
638  Double_t mean_area = 0;
639  try {
640  if(0 == remove) {
641  fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
642  } else {
643  std::vector<fastjet::PseudoJet> input_jets = sorted_by_pt(fClustSeq->inclusive_jets());
644  input_jets.erase(input_jets.begin(), input_jets.begin() + remove);
645  fClustSeq->get_median_rho_and_sigma(input_jets, *fRange, fUseArea4Vector, median, sigma, mean_area);
646  input_jets.clear();
647  }
648  } catch (fj::Error) {
649  AliError(" [w] FJ Exception caught.");
650  median = -1.;
651  sigma = -1;
652  }
653 }
654 
655 //_________________________________________________________________________________________________
656 Int_t AliFJWrapper::Run()
657 {
658  // Run the actual jet finder.
659 
660  if (fAreaType == fj::voronoi_area) {
661  // Rfact - check dependence - default is 1.
662  // NOTE: hardcoded variable!
663  fVorAreaSpec = new fj::VoronoiAreaSpec(1.);
664  fAreaDef = new fj::AreaDefinition(*fVorAreaSpec);
665  } else {
666  fGhostedAreaSpec = new fj::GhostedAreaSpec(fMaxRap,
667  fNGhostRepeats,
668  fGhostArea,
669  fGridScatter,
670  fKtScatter,
671  fMeanGhostKt);
672 
673  fAreaDef = new fj::AreaDefinition(*fGhostedAreaSpec, fAreaType);
674  }
675 
676  // this is acceptable by fastjet:
677 #ifndef FASTJET_VERSION
678  fRange = new fj::RangeDefinition(fMaxRap - 0.95 * fR);
679 #else
680  fRange = new fj::Selector(fj::SelectorAbsRapMax(fMaxRap - 0.95 * fR));
681 #endif
682 
683  if (fAlgor == fj::plugin_algorithm) {
684  if (fPluginAlgor == 0) {
685  // SIS CONE ALGOR
686  // NOTE: hardcoded split parameter
687  Double_t overlap_threshold = 0.75; // NOTE: this actually splits a lot: thr/min(pt1,pt2)
688  fPlugin = new fj::SISConePlugin(fR,
689  overlap_threshold,
690  0, //search of stable cones - zero = until no more
691  1.0); // this should be seed effectively for proto jets
692  fJetDef = new fastjet::JetDefinition(fPlugin);
693  } else if (fPluginAlgor == 1) {
694  // CDF cone
695  // NOTE: hardcoded split parameter
696  Double_t overlap_threshold = 0.75; // NOTE: this actually splits a lot: thr/min(pt1,pt2)
697  fPlugin = new fj::CDFMidPointPlugin(fR,
698  overlap_threshold,
699  1.0, //search of stable cones - zero = until no more
700  1.0); // this should be seed effectively for proto jets
701  fJetDef = new fastjet::JetDefinition(fPlugin);
702  } else {
703  AliError("[e] Unrecognized plugin number!");
704  }
705  } else {
706  fJetDef = new fj::JetDefinition(fAlgor, fR, fScheme, fStrategy);
707  }
708 
709  try {
710  fClustSeq = new fj::ClusterSequenceArea(fInputVectors, *fJetDef, *fAreaDef);
711  if(fEventSub){
712  DoEventConstituentSubtraction();
713  fClustSeqES = new fj::ClusterSequenceArea(fEventSubCorrectedVectors, *fJetDef, *fAreaDef);
714  }
715  } catch (fj::Error) {
716  AliError(" [w] FJ Exception caught.");
717  return -1;
718  }
719 
720  // FJ3 :: Define an JetMedianBackgroundEstimator just in case it will be used
721 #ifdef FASTJET_VERSION
722  fBkrdEstimator = new fj::JetMedianBackgroundEstimator(fj::SelectorAbsRapMax(fMaxRap));
723 #endif
724 
725  if (fLegacyMode) { SetLegacyFJ(); } // for FJ 2.x even if fLegacyMode is set, SetLegacyFJ is dummy
726 
727  // inclusive jets:
728  fInclusiveJets.clear();
729  fEventSubJets.clear();
730  fInclusiveJets = fClustSeq->inclusive_jets(0.0);
731  if(fEventSub) fEventSubJets = fClustSeqES->inclusive_jets(0.0);
732 
733  return 0;
734 }
735 
736 //_________________________________________________________________________________________________
737 Int_t AliFJWrapper::Filter()
738 {
739 //
740 // AliFJWrapper::Filter
741 //
742 
743  fJetDef = new fj::JetDefinition(fAlgor, fR, fScheme, fStrategy);
744 
745  if (fDoFilterArea) {
746  if (fInputGhosts.size()>0) {
747  try {
748  fClustSeqActGhosts = new fj::ClusterSequenceActiveAreaExplicitGhosts(fInputVectors,
749  *fJetDef,
750  fInputGhosts,
751  fGhostArea);
752  } catch (fj::Error) {
753  AliError(" [w] FJ Exception caught.");
754  return -1;
755  }
756 
757  fFilteredJets.clear();
758  fFilteredJets = fClustSeqActGhosts->inclusive_jets(0.0);
759  } else {
760  return -1;
761  }
762  } else {
763  try {
764  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
765  } catch (fj::Error) {
766  AliError(" [w] FJ Exception caught.");
767  return -1;
768  }
769 
770  fFilteredJets.clear();
771  fFilteredJets = fClustSeqSA->inclusive_jets(0.0);
772  }
773 
774  return 0;
775 }
776 
777 //_________________________________________________________________________________________________
778 void AliFJWrapper::SetLegacyFJ()
779 {
780  // This methods enable legacy behaviour (FastJet 2.x) when AliROOT is compiled with FastJet 3.x
781 #ifdef FASTJET_VERSION
782  std::cout << "WARNING! Setting FastJet in legacy mode" << std::endl;
783  if (fGhostedAreaSpec) { fGhostedAreaSpec->set_fj2_placement(kTRUE); }
784  if (fBkrdEstimator) {
785  fBkrdEstimator->set_provide_fj2_sigma(kTRUE);
786  fBkrdEstimator->set_use_area_4vector(kFALSE);
787  }
788 #endif
789 }
790 
791 //_________________________________________________________________________________________________
792 void AliFJWrapper::SubtractBackground(Double_t median_pt)
793 {
794  // Subtract the background (specify the value - see below the meaning).
795  // Negative argument means the bg will be determined with the current algorithm
796  // this is the default behavior. Zero is allowed
797  // Note: user may set the switch for area4vector based subtraction.
798 
799  Double_t median = 0;
800  Double_t sigma = 0;
801  Double_t mean_area = 0;
802 
803  // clear the subtracted jet pt's vector<double>
804  fSubtractedJetsPt.clear();
805 
806  // check what was specified (default is -1)
807  if (median_pt < 0) {
808  try {
809  fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
810  }
811 
812  catch (fj::Error) {
813  AliError(" [w] FJ Exception caught.");
814  median = -9999.;
815  sigma = -1;
816  fMedUsedForBgSub = median;
817  return;
818  }
819  fMedUsedForBgSub = median;
820  } else {
821  // we do not know the sigma in this case
822  sigma = -1;
823  if (0.0 == median_pt) {
824  // AliWarning(Form(" [w] Median specified for bg subtraction is ZERO: %f", median_pt )) << endl;
825  fMedUsedForBgSub = 0.;
826  } else {
827  fMedUsedForBgSub = median_pt;
828  }
829  }
830 
831  // subtract:
832  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
833  if ( fUseArea4Vector ) {
834  // subtract the background using the area4vector
835  fj::PseudoJet area4v = fClustSeq->area_4vector(fInclusiveJets[i]);
836  fj::PseudoJet jet_sub = fInclusiveJets[i] - area4v * fMedUsedForBgSub;
837  fSubtractedJetsPt.push_back(jet_sub.perp()); // here we put only the pt of the jet - note: this can be negative
838  } else {
839  // subtract the background using scalars
840  // fj::PseudoJet jet_sub = fInclusiveJets[i] - area * fMedUsedForBgSub_;
841  Double_t area = fClustSeq->area(fInclusiveJets[i]);
842  // standard subtraction
843  Double_t pt_sub = fInclusiveJets[i].perp() - fMedUsedForBgSub * area;
844  fSubtractedJetsPt.push_back(pt_sub); // here we put only the pt of the jet - note: this can be negative
845  }
846  }
847 }
848 
849 //_________________________________________________________________________________________________
850 Int_t AliFJWrapper::DoGenericSubtractionJetMass() {
851  //Do generic subtraction for jet mass
852 #ifdef FASTJET_VERSION
853  CreateGenSub();
854 
855  // Define jet shape
856  AliJetShapeMass shapeMass;
857 
858  // clear the generic subtractor info vector
859  fGenSubtractorInfoJetMass.clear();
860  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
861  fj::contrib::GenericSubtractorInfo info;
862  if(fInclusiveJets[i].perp()>1.e-4)
863  double subtracted_shape = (*fGenSubtractor)(shapeMass, fInclusiveJets[i], info);
864  fGenSubtractorInfoJetMass.push_back(info);
865  }
866 #endif
867  return 0;
868 }
869 
870 //_________________________________________________________________________________________________
871 Int_t AliFJWrapper::DoGenericSubtractionGR(Int_t ijet) {
872  //Do generic subtraction for jet mass
873 #ifdef FASTJET_VERSION
874  CreateGenSub();
875 
876  if(ijet>fInclusiveJets.size()) return 0;
877 
878  fGRNumerator.clear();
879  fGRDenominator.clear();
880  fGRNumeratorSub.clear();
881  fGRDenominatorSub.clear();
882 
883  // Define jet shape
884  for(Double_t r = 0.; r<fRMax; r+=fDRStep) {
885  AliJetShapeGRNum shapeGRNum(r,fDRStep);
886  AliJetShapeGRDen shapeGRDen(r,fDRStep);
887 
888  // clear the generic subtractor info vector
889  fGenSubtractorInfoGRNum.clear();
890  fGenSubtractorInfoGRDen.clear();
891  fj::contrib::GenericSubtractorInfo infoNum;
892  fj::contrib::GenericSubtractorInfo infoDen;
893  if(fInclusiveJets[ijet].perp()>1.e-4) {
894  double sub_num = (*fGenSubtractor)(shapeGRNum, fInclusiveJets[ijet], infoNum);
895  double sub_den = (*fGenSubtractor)(shapeGRDen, fInclusiveJets[ijet], infoDen);
896  }
897  fGenSubtractorInfoGRNum.push_back(infoNum);
898  fGenSubtractorInfoGRDen.push_back(infoDen);
899  fGRNumerator.push_back(infoNum.unsubtracted());
900  fGRDenominator.push_back(infoDen.unsubtracted());
901  fGRNumeratorSub.push_back(infoNum.second_order_subtracted());
902  fGRDenominatorSub.push_back(infoDen.second_order_subtracted());
903  }
904 #endif
905  return 0;
906 }
907 //_________________________________________________________________________________________________
908 Int_t AliFJWrapper::DoGenericSubtractionJetAngularity() {
909  //Do generic subtraction for jet mass
910 #ifdef FASTJET_VERSION
911  CreateGenSub();
912 
913  // Define jet shape
914  AliJetShapeAngularity shapeAngularity;
915 
916  // clear the generic subtractor info vector
917  fGenSubtractorInfoJetAngularity.clear();
918  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
919  fj::contrib::GenericSubtractorInfo infoAng;
920  if(fInclusiveJets[i].perp()>1.e-4)
921  double subtracted_shape = (*fGenSubtractor)(shapeAngularity, fInclusiveJets[i], infoAng);
922  fGenSubtractorInfoJetAngularity.push_back(infoAng);
923  }
924 #endif
925  return 0;
926 }
927 //_________________________________________________________________________________________________
928 Int_t AliFJWrapper::DoGenericSubtractionJetpTD() {
929  //Do generic subtraction for jet mass
930 #ifdef FASTJET_VERSION
931  CreateGenSub();
932 
933  // Define jet shape
934  AliJetShapepTD shapepTD;
935 
936  // clear the generic subtractor info vector
937  fGenSubtractorInfoJetpTD.clear();
938  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
939  fj::contrib::GenericSubtractorInfo infopTD;
940  if(fInclusiveJets[i].perp()>1.e-4)
941  double subtracted_shape = (*fGenSubtractor)(shapepTD, fInclusiveJets[i], infopTD);
942  fGenSubtractorInfoJetpTD.push_back(infopTD);
943  }
944 #endif
945  return 0;
946 }
947 //_________________________________________________________________________________________________
948 Int_t AliFJWrapper::DoGenericSubtractionJetCircularity() {
949  //Do generic subtraction for jet mass
950 #ifdef FASTJET_VERSION
951  CreateGenSub();
952 
953  // Define jet shape
954  AliJetShapeCircularity shapecircularity;
955 
956  // clear the generic subtractor info vector
957  fGenSubtractorInfoJetCircularity.clear();
958  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
959  fj::contrib::GenericSubtractorInfo infoCirc;
960  if(fInclusiveJets[i].perp()>1.e-4)
961  double subtracted_shape = (*fGenSubtractor)(shapecircularity, fInclusiveJets[i], infoCirc);
962  fGenSubtractorInfoJetCircularity.push_back(infoCirc);
963  }
964 #endif
965  return 0;
966 }
967 //_________________________________________________________________________________________________
968 Int_t AliFJWrapper::DoGenericSubtractionJetSigma2() {
969  //Do generic subtraction for jet mass
970 #ifdef FASTJET_VERSION
971  CreateGenSub();
972 
973  // Define jet shape
974  AliJetShapeSigma2 shapesigma2;
975 
976  // clear the generic subtractor info vector
977  fGenSubtractorInfoJetSigma2.clear();
978  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
979  fj::contrib::GenericSubtractorInfo infoSigma;
980  if(fInclusiveJets[i].perp()>1.e-4)
981  double subtracted_shape = (*fGenSubtractor)(shapesigma2, fInclusiveJets[i], infoSigma);
982  fGenSubtractorInfoJetSigma2.push_back(infoSigma);
983  }
984 #endif
985  return 0;
986 }
987 //_________________________________________________________________________________________________
988 Int_t AliFJWrapper::DoGenericSubtractionJetConstituent() {
989  //Do generic subtraction for jet mass
990 #ifdef FASTJET_VERSION
991  CreateGenSub();
992 
993  // Define jet shape
994  AliJetShapeConstituent shapeconst;
995 
996  // clear the generic subtractor info vector
997  fGenSubtractorInfoJetConstituent.clear();
998  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
999  fj::contrib::GenericSubtractorInfo infoConst;
1000  if(fInclusiveJets[i].perp()>1.e-4)
1001  double subtracted_shape = (*fGenSubtractor)(shapeconst, fInclusiveJets[i], infoConst);
1002  fGenSubtractorInfoJetConstituent.push_back(infoConst);
1003  }
1004 #endif
1005  return 0;
1006 }
1007 
1008 //_________________________________________________________________________________________________
1009 Int_t AliFJWrapper::DoGenericSubtractionJetLeSub() {
1010  //Do generic subtraction for jet mass
1011 #ifdef FASTJET_VERSION
1012  CreateGenSub();
1013 
1014  // Define jet shape
1015  AliJetShapeLeSub shapeLeSub;
1016 
1017  // clear the generic subtractor info vector
1018  fGenSubtractorInfoJetLeSub.clear();
1019  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1020  fj::contrib::GenericSubtractorInfo infoLeSub;
1021  if(fInclusiveJets[i].perp()>1.e-4)
1022  double subtracted_shape = (*fGenSubtractor)(shapeLeSub, fInclusiveJets[i], infoLeSub);
1023  fGenSubtractorInfoJetLeSub.push_back(infoLeSub);
1024  }
1025 #endif
1026  return 0;
1027 }
1028 
1029 //_________________________________________________________________________________________________
1030 Int_t AliFJWrapper::DoGenericSubtractionJet1subjettiness_kt() {
1031  //Do generic subtraction for 1subjettiness
1032 #ifdef FASTJET_VERSION
1033  CreateGenSub();
1034 
1035  // Define jet shape
1036  AliJetShape1subjettiness_kt shape1subjettiness_kt;
1037 
1038  // clear the generic subtractor info vector
1039  fGenSubtractorInfoJet1subjettiness_kt.clear();
1040  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1041  fj::contrib::GenericSubtractorInfo info1subjettiness_kt;
1042  if(fInclusiveJets[i].perp()>1.e-4)
1043  double subtracted_shape = (*fGenSubtractor)(shape1subjettiness_kt, fInclusiveJets[i], info1subjettiness_kt);
1044  fGenSubtractorInfoJet1subjettiness_kt.push_back(info1subjettiness_kt);
1045  }
1046 #endif
1047  return 0;
1048 }
1049 
1050 //_________________________________________________________________________________________________
1051 Int_t AliFJWrapper::DoGenericSubtractionJet2subjettiness_kt() {
1052  //Do generic subtraction for 2subjettiness
1053 #ifdef FASTJET_VERSION
1054  CreateGenSub();
1055 
1056  // Define jet shape
1057  AliJetShape2subjettiness_kt shape2subjettiness_kt;
1058 
1059  // clear the generic subtractor info vector
1060  fGenSubtractorInfoJet2subjettiness_kt.clear();
1061  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1062  fj::contrib::GenericSubtractorInfo info2subjettiness_kt;
1063  if(fInclusiveJets[i].perp()>1.e-4)
1064  double subtracted_shape = (*fGenSubtractor)(shape2subjettiness_kt, fInclusiveJets[i], info2subjettiness_kt);
1065  fGenSubtractorInfoJet2subjettiness_kt.push_back(info2subjettiness_kt);
1066  }
1067 #endif
1068  return 0;
1069 }
1070 
1071 //_________________________________________________________________________________________________
1072 Int_t AliFJWrapper::DoGenericSubtractionJet3subjettiness_kt() {
1073  //Do generic subtraction for 3subjettiness
1074 #ifdef FASTJET_VERSION
1075  CreateGenSub();
1076 
1077  // Define jet shape
1078  AliJetShape3subjettiness_kt shape3subjettiness_kt;
1079 
1080  // clear the generic subtractor info vector
1081  fGenSubtractorInfoJet3subjettiness_kt.clear();
1082  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1083  fj::contrib::GenericSubtractorInfo info3subjettiness_kt;
1084  if(fInclusiveJets[i].perp()>1.e-4)
1085  double subtracted_shape = (*fGenSubtractor)(shape3subjettiness_kt, fInclusiveJets[i], info3subjettiness_kt);
1086  fGenSubtractorInfoJet3subjettiness_kt.push_back(info3subjettiness_kt);
1087  }
1088 #endif
1089  return 0;
1090 }
1091 
1092 //_________________________________________________________________________________________________
1093 Int_t AliFJWrapper::DoGenericSubtractionJetOpeningAngle_kt() {
1094  //Do generic subtraction for 2subjettiness axes opening angle
1095 #ifdef FASTJET_VERSION
1096  CreateGenSub();
1097 
1098  // Define jet shape
1099  AliJetShapeOpeningAngle_kt shapeOpeningAngle_kt;
1100 
1101  // clear the generic subtractor info vector
1102  fGenSubtractorInfoJetOpeningAngle_kt.clear();
1103  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1104  fj::contrib::GenericSubtractorInfo infoOpeningAngle_kt;
1105  if(fInclusiveJets[i].perp()>1.e-4)
1106  double subtracted_shape = (*fGenSubtractor)(shapeOpeningAngle_kt, fInclusiveJets[i], infoOpeningAngle_kt);
1107  fGenSubtractorInfoJetOpeningAngle_kt.push_back(infoOpeningAngle_kt);
1108  }
1109 #endif
1110  return 0;
1111 }
1112 
1113 //_________________________________________________________________________________________________
1114 Int_t AliFJWrapper::DoConstituentSubtraction() {
1115  //Do constituent subtraction
1116 #ifdef FASTJET_VERSION
1117  CreateConstituentSub();
1118  // fConstituentSubtractor->set_alpha(/* double alpha */);
1119  // fConstituentSubtractor->set_max_deltaR(/* double max_deltaR */);
1120 
1121  //clear constituent subtracted jets
1122  fConstituentSubtrJets.clear();
1123  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1124  fj::PseudoJet subtracted_jet(0.,0.,0.,0.);
1125  if(fInclusiveJets[i].perp()>0.)
1126  subtracted_jet = (*fConstituentSubtractor)(fInclusiveJets[i]);
1127  fConstituentSubtrJets.push_back(subtracted_jet);
1128  }
1129  if(fConstituentSubtractor) { delete fConstituentSubtractor; fConstituentSubtractor = NULL; }
1130 
1131 #endif
1132  return 0;
1133 }
1134 
1135 //_________________________________________________________________________________________________
1136 Int_t AliFJWrapper::DoEventConstituentSubtraction() {
1137  //Do constituent subtraction
1138 #ifdef FASTJET_VERSION
1139  CreateEventConstituentSub();
1140  if(fUseMaxDelR) fEventConstituentSubtractor->set_max_standardDeltaR(fMaxDelR);
1141  fEventSubCorrectedVectors = fEventConstituentSubtractor->subtract_event(fEventSubInputVectors,fMaxRap); //second argument max rap?
1142  //clear constituent subtracted jets
1143  if(fEventConstituentSubtractor) { delete fEventConstituentSubtractor; fEventConstituentSubtractor = NULL; }
1144 
1145 #endif
1146  return 0;
1147 }
1148 
1149 //_________________________________________________________________________________________________
1150 Int_t AliFJWrapper::DoSoftDrop() {
1151  //Do grooming
1152 #ifdef FASTJET_VERSION
1153  CreateSoftDrop();
1154 
1155  //clear groomed jets
1156  fGroomedJets.clear();
1157  //fastjet::Subtractor fjsub (fBkrdEstimator);
1158  //fSoftDrop->set_subtractor(&fjsub);
1159  //fSoftDrop->set_input_jet_is_subtracted(false); //??
1160 
1161  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1162  fj::PseudoJet groomed_jet(0.,0.,0.,0.);
1163  if(fInclusiveJets[i].perp()>0.){
1164  groomed_jet = (*fSoftDrop)(fInclusiveJets[i]);
1165  groomed_jet.set_user_index(i); //index of the corresponding inclusve jet
1166  if(groomed_jet!=0) fGroomedJets.push_back(groomed_jet);
1167  }
1168 
1169  }
1170  if(fSoftDrop) { delete fSoftDrop; fSoftDrop = NULL; }
1171 
1172 #endif
1173  return 0;
1174 }
1175 
1176 //_________________________________________________________________________________________________
1177 Int_t AliFJWrapper::CreateSoftDrop() {
1178  //Do grooming
1179  #ifdef FASTJET_VERSION
1180  if (fSoftDrop) { delete fSoftDrop; } // protect against memory leaks
1181 
1182  fSoftDrop = new fj::contrib::SoftDrop(fBeta,fZcut);
1183  fSoftDrop->set_verbose_structure(kTRUE);
1184 
1185  #endif
1186  return 0;
1187 }
1188 
1189 
1190 //_________________________________________________________________________________________________
1191 Int_t AliFJWrapper::CreateGenSub() {
1192  //Do generic subtraction for jet mass
1193  #ifdef FASTJET_VERSION
1194  if (fGenSubtractor) { delete fGenSubtractor; } // protect against memory leaks
1195 
1196  if (fUseExternalBkg)
1197  { fGenSubtractor = new fj::contrib::GenericSubtractor(fRho,fRhom); }
1198  else
1199  {
1200  fGenSubtractor = new fj::contrib::GenericSubtractor(fBkrdEstimator);
1201  #if FASTJET_VERSION_NUMBER >= 30100
1202  fGenSubtractor->set_common_bge_for_rho_and_rhom(); // see contrib 1.020 GenericSubtractor.hh line 62
1203  #endif
1204  }
1205 
1206  #endif
1207  return 0;
1208 }
1209 
1210 //_________________________________________________________________________________________________
1211 Int_t AliFJWrapper::CreateConstituentSub() {
1212  //Do generic subtraction for jet mass
1213  #ifdef FASTJET_VERSION
1214  if (fConstituentSubtractor) { delete fConstituentSubtractor; } // protect against memory leaks
1215 
1216  // see ConstituentSubtractor.hh signatures
1217  // ConstituentSubtractor(double rho, double rhom=0, double alpha=0, double maxDeltaR=-1)
1218  if (fUseExternalBkg) { fConstituentSubtractor = new fj::contrib::ConstituentSubtractor(fRho,fRhom); }
1219  else { fConstituentSubtractor = new fj::contrib::ConstituentSubtractor(fBkrdEstimator); }
1220 
1221  #endif
1222  return 0;
1223 }
1224 
1225 //_________________________________________________________________________________________________
1226 Int_t AliFJWrapper::CreateEventConstituentSub() {
1227  //Do generic subtraction for jet mass
1228  #ifdef FASTJET_VERSION
1229  if (fEventConstituentSubtractor) { delete fEventConstituentSubtractor; } // protect against memory leaks
1230  // see ConstituentSubtractor.hh signatures
1231  // ConstituentSubtractor(double rho, double rhom=0, double alpha=0, double maxDeltaR=-1)
1232  if (fUseExternalBkg) fEventConstituentSubtractor = new fj::contrib::ConstituentSubtractor(fRho,fRhom);
1233  else fEventConstituentSubtractor = new fj::contrib::ConstituentSubtractor(fBkrdEstimator);
1234 
1235  #endif
1236  return 0;
1237 }
1238 
1239 //_________________________________________________________________________________________________
1240 void AliFJWrapper::SetupAlgorithmfromOpt(const char *option)
1241 {
1242  // Setup algorithm from char.
1243 
1244  std::string opt(option);
1245 
1246  if (!opt.compare("kt")) fAlgor = fj::kt_algorithm;
1247  if (!opt.compare("antikt")) fAlgor = fj::antikt_algorithm;
1248  if (!opt.compare("cambridge")) fAlgor = fj::cambridge_algorithm;
1249  if (!opt.compare("genkt")) fAlgor = fj::genkt_algorithm;
1250  if (!opt.compare("cambridge_passive")) fAlgor = fj::cambridge_for_passive_algorithm;
1251  if (!opt.compare("genkt_passive")) fAlgor = fj::genkt_for_passive_algorithm;
1252  if (!opt.compare("ee_kt")) fAlgor = fj::ee_kt_algorithm;
1253  if (!opt.compare("ee_genkt")) fAlgor = fj::ee_genkt_algorithm;
1254  if (!opt.compare("plugin")) fAlgor = fj::plugin_algorithm;
1255 }
1256 
1257 //_________________________________________________________________________________________________
1258 void AliFJWrapper::SetupAreaTypefromOpt(const char *option)
1259 {
1260  // Setup area type from char.
1261 
1262  std::string opt(option);
1263 
1264  if (!opt.compare("active")) fAreaType = fj::active_area;
1265  if (!opt.compare("invalid")) fAreaType = fj::invalid_area;
1266  if (!opt.compare("active_area_explicit_ghosts")) fAreaType = fj::active_area_explicit_ghosts;
1267  if (!opt.compare("one_ghost_passive")) fAreaType = fj::one_ghost_passive_area;
1268  if (!opt.compare("passive")) fAreaType = fj::passive_area;
1269  if (!opt.compare("voronoi")) fAreaType = fj::voronoi_area;
1270 }
1271 
1272 //_________________________________________________________________________________________________
1273 void AliFJWrapper::SetupSchemefromOpt(const char *option)
1274 {
1275  //
1276  // setup scheme from char
1277  //
1278 
1279  std::string opt(option);
1280 
1281  if (!opt.compare("BIpt")) fScheme = fj::BIpt_scheme;
1282  if (!opt.compare("BIpt2")) fScheme = fj::BIpt2_scheme;
1283  if (!opt.compare("E")) fScheme = fj::E_scheme;
1284  if (!opt.compare("pt")) fScheme = fj::pt_scheme;
1285  if (!opt.compare("pt2")) fScheme = fj::pt2_scheme;
1286  if (!opt.compare("Et")) fScheme = fj::Et_scheme;
1287  if (!opt.compare("Et2")) fScheme = fj::Et2_scheme;
1288 }
1289 
1290 //_________________________________________________________________________________________________
1291 void AliFJWrapper::SetupStrategyfromOpt(const char *option)
1292 {
1293  // Setup strategy from char.
1294 
1295  std::string opt(option);
1296 
1297  if (!opt.compare("Best")) fStrategy = fj::Best;
1298  if (!opt.compare("N2MinHeapTiled")) fStrategy = fj::N2MinHeapTiled;
1299  if (!opt.compare("N2Tiled")) fStrategy = fj::N2Tiled;
1300  if (!opt.compare("N2PoorTiled")) fStrategy = fj::N2PoorTiled;
1301  if (!opt.compare("N2Plain")) fStrategy = fj::N2Plain;
1302  if (!opt.compare("N3Dumb")) fStrategy = fj::N3Dumb;
1303  if (!opt.compare("NlnN")) fStrategy = fj::NlnN;
1304  if (!opt.compare("NlnN3pi")) fStrategy = fj::NlnN3pi;
1305  if (!opt.compare("NlnN4pi")) fStrategy = fj::NlnN4pi;
1306  if (!opt.compare("NlnNCam4pi")) fStrategy = fj::NlnNCam4pi;
1307  if (!opt.compare("NlnNCam2pi2R")) fStrategy = fj::NlnNCam2pi2R;
1308  if (!opt.compare("NlnNCam")) fStrategy = fj::NlnNCam;
1309  if (!opt.compare("plugin")) fStrategy = fj::plugin_strategy;
1310 }
1311 
1312 Double_t AliFJWrapper::NSubjettiness(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){
1313 
1314 
1315  //Option 0=Nsubjettiness result, 1=opening angle between axes in Eta-Phi plane, 2=Distance between axes in Eta-Phi plane
1316 
1317  fJetDef = new fj::JetDefinition(fAlgor, fR*2, fScheme, fStrategy ); //the *2 is becasue of a handful of jets that end up missing a track for some reason.
1318 
1319  try {
1320  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
1321  // ClustSeqSA = new fastjet::ClusterSequenceArea(fInputVectors, *fJetDef, *fAreaDef);
1322  } catch (fj::Error) {
1323  AliError(" [w] FJ Exception caught.");
1324  return -1;
1325  }
1326  fFilteredJets.clear();
1327  fFilteredJets = fClustSeqSA->inclusive_jets(fMinJetPt-0.1); //becasue this is < not <=
1328  Double_t Result=-1;
1329  Double_t Result_SoftDrop=-1;
1330  std::vector<fastjet::PseudoJet> SubJet_Axes; //this is actually not subjet axis...but just axis
1331  fj::PseudoJet SubJet1_Axis;
1332  fj::PseudoJet SubJet2_Axis;
1333  std::vector<fastjet::PseudoJet> SubJets;
1334  fj::PseudoJet SubJet1;
1335  fj::PseudoJet SubJet2;
1336  if (Algorithm==0){
1337  Beta = 1.0;
1338  fR=0.4;
1339  if (Measure==0){
1340  fj::contrib::Nsubjettiness nSub(N, fj::contrib::KT_Axes(), fj::contrib::NormalizedMeasure(Beta,fR));
1341  Result= nSub.result(fFilteredJets[0]);
1342  SubJet_Axes=nSub.currentAxes();
1343  SubJets=nSub.currentSubjets();
1344  }
1345  else if (Measure==1){
1346  fj::contrib::Nsubjettiness nSub(N, fj::contrib::KT_Axes(), fj::contrib::UnnormalizedMeasure(Beta));
1347  Result= nSub.result(fFilteredJets[0]);
1348  SubJet_Axes=nSub.currentAxes();
1349  SubJets=nSub.currentSubjets();
1350  }
1351  }
1352  else if (Algorithm==1) {
1353  fj::contrib::Nsubjettiness nSub(N, fj::contrib::CA_Axes(), fj::contrib::NormalizedMeasure(Beta,fR));
1354  Result= nSub.result(fFilteredJets[0]);
1355  SubJet_Axes=nSub.currentAxes();
1356  SubJets=nSub.currentSubjets();
1357  }
1358  else if (Algorithm==2){
1359  fj::contrib::Nsubjettiness nSub(N, fj::contrib::AntiKT_Axes(Radius), fj::contrib::NormalizedMeasure(Beta,fR));
1360  Result= nSub.result(fFilteredJets[0]);
1361  SubJet_Axes=nSub.currentAxes();
1362  SubJets=nSub.currentSubjets();
1363  }
1364  else if (Algorithm==3) {
1365  fj::contrib::Nsubjettiness nSub(N, fj::contrib::WTA_KT_Axes(), fj::contrib::NormalizedMeasure(Beta,fR));
1366  Result= nSub.result(fFilteredJets[0]);
1367  SubJet_Axes=nSub.currentAxes();
1368  SubJets=nSub.currentSubjets();
1369  }
1370  else if (Algorithm==4) {
1371  fj::contrib::Nsubjettiness nSub(N, fj::contrib::WTA_CA_Axes(), fj::contrib::NormalizedMeasure(Beta,fR));
1372  Result= nSub.result(fFilteredJets[0]);
1373  SubJet_Axes=nSub.currentAxes();
1374  SubJets=nSub.currentSubjets();
1375  }
1376  else if (Algorithm==5){
1377  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_KT_Axes(), fj::contrib::NormalizedMeasure(Beta,fR));
1378  Result= nSub.result(fFilteredJets[0]);
1379  SubJet_Axes=nSub.currentAxes();
1380  SubJets=nSub.currentSubjets();
1381  }
1382  else if (Algorithm==6){
1383  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_CA_Axes(), fj::contrib::NormalizedMeasure(Beta,fR));
1384  Result= nSub.result(fFilteredJets[0]);
1385  SubJet_Axes=nSub.currentAxes();
1386  SubJets=nSub.currentSubjets();
1387  }
1388  else if (Algorithm==7){
1389  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_AntiKT_Axes(Radius), fj::contrib::NormalizedMeasure(Beta,fR));
1390  Result= nSub.result(fFilteredJets[0]);
1391  SubJet_Axes=nSub.currentAxes();
1392  SubJets=nSub.currentSubjets();
1393  }
1394  else if (Algorithm==8){
1395  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_WTA_KT_Axes(), fj::contrib::NormalizedMeasure(Beta,fR));
1396  Result= nSub.result(fFilteredJets[0]);
1397  SubJet_Axes=nSub.currentAxes();
1398  SubJets=nSub.currentSubjets();
1399  }
1400  else if (Algorithm==9){
1401  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_WTA_CA_Axes(), fj::contrib::NormalizedMeasure(Beta,fR));
1402  Result= nSub.result(fFilteredJets[0]);
1403  SubJet_Axes=nSub.currentAxes();
1404  SubJets=nSub.currentSubjets();
1405  }
1406  else if (Algorithm==10){
1407  fj::contrib::Nsubjettiness nSub(N, fj::contrib::MultiPass_Axes(100), fj::contrib::NormalizedMeasure(Beta,fR));
1408  Result= nSub.result(fFilteredJets[0]);
1409  SubJet_Axes=nSub.currentAxes();
1410  SubJets=nSub.currentSubjets();
1411  }
1412 
1413 
1414 
1415  SubJet1_Axis=SubJet_Axes[0];
1416  Double_t SubJet1_Eta=SubJet1_Axis.pseudorapidity();
1417  Double_t SubJet2_Eta;
1418  Double_t SubJet1_Phi=SubJet1_Axis.phi();
1419  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
1420  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
1421  Double_t SubJet2_Phi;
1422  Double_t DeltaPhi=-5;
1423  if (SubJet_Axes.size()>1){
1424  SubJet2_Axis=SubJet_Axes[1];
1425  SubJet2_Eta=SubJet2_Axis.pseudorapidity();
1426  SubJet2_Phi=SubJet2_Axis.phi();
1427  if(SubJet2_Phi < -1*TMath::Pi()) SubJet2_Phi += (2*TMath::Pi());
1428  else if (SubJet2_Phi > TMath::Pi()) SubJet2_Phi -= (2*TMath::Pi());
1429  DeltaPhi=SubJet1_Phi-SubJet2_Phi;
1430  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1431  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1432  }
1433 
1434 
1435  SubJet1=SubJets[0];
1436  Double_t SubJet1Eta=SubJet1.pseudorapidity();
1437  Double_t SubJet2Eta;
1438  Double_t SubJet1Phi=SubJet1.phi();
1439  if(SubJet1Phi < -1*TMath::Pi()) SubJet1Phi += (2*TMath::Pi());
1440  else if (SubJet1Phi > TMath::Pi()) SubJet1Phi -= (2*TMath::Pi());
1441  Double_t SubJet2Phi;
1442  Double_t DeltaPhiSubJets=-5;
1443  Double_t SubJet1LeadingTrackPt=-3.0;
1444  Double_t SubJet2LeadingTrackPt=-3.0;
1445  std::vector<fj::PseudoJet> SubJet1Tracks = SubJet1.constituents();
1446  for (Int_t i=0; i<SubJet1Tracks.size(); i++){
1447  if (SubJet1Tracks[i].perp() > SubJet1LeadingTrackPt) SubJet1LeadingTrackPt=SubJet1Tracks[i].perp();
1448  }
1449  if (SubJet_Axes.size()>1){
1450  SubJet2=SubJets[1];
1451  SubJet2Eta=SubJet2.pseudorapidity();
1452  SubJet2Phi=SubJet2.phi();
1453  if(SubJet2Phi < -1*TMath::Pi()) SubJet2Phi += (2*TMath::Pi());
1454  else if (SubJet2Phi > TMath::Pi()) SubJet2Phi -= (2*TMath::Pi());
1455  DeltaPhiSubJets=SubJet1Phi-SubJet2Phi;
1456  if(DeltaPhiSubJets < -1*TMath::Pi()) DeltaPhiSubJets += (2*TMath::Pi());
1457  else if (DeltaPhiSubJets > TMath::Pi()) DeltaPhiSubJets -= (2*TMath::Pi());
1458  std::vector<fj::PseudoJet> SubJet2Tracks = SubJet2.constituents();
1459  for (Int_t i=0; i<SubJet2Tracks.size(); i++){
1460  if (SubJet2Tracks[i].perp() > SubJet2LeadingTrackPt) SubJet2LeadingTrackPt=SubJet2Tracks[i].perp();
1461  }
1462  }
1463 
1464  //Added for quality control of the DeltaR-Nsubjettiness variable (comparing Nsubjettiness and soft drop results)
1465  Beta_SD=0.0;
1466  ZCut=0.1;
1467  fj::contrib::SoftDrop Soft_Drop(Beta_SD,ZCut);
1468  Soft_Drop.set_tagging_mode(); //if the first two subjets fail the soft drop criteria a jet = 0 is returned
1469  fj::PseudoJet Soft_Dropped_Jet=Soft_Drop(fFilteredJets[0]);
1470  if (Soft_Dropped_Jet==0) Result_SoftDrop=-1;
1471  else{
1472  if (Option==3) Result_SoftDrop=Soft_Dropped_Jet.structure_of<fj::contrib::SoftDrop>().delta_R();
1473  if (Option==4) Result_SoftDrop=Soft_Dropped_Jet.structure_of<fj::contrib::SoftDrop>().symmetry();
1474  }
1475 
1476 
1477  if (Option==0) return Result;
1478  else if (Option==1 && SubJet_Axes.size()>1 && N==2) return TMath::Sqrt(TMath::Power(SubJet1_Eta-SubJet2_Eta,2)+TMath::Power(DeltaPhi,2));
1479  else if (Option==2 && SubJet_Axes.size()>1 && N==2) return TMath::Sqrt(TMath::Power(SubJet1_Eta-SubJet2_Eta,2)+TMath::Power(DeltaPhi,2));
1480  else if ((Option==3 || Option==4) && N==2) return Result_SoftDrop;
1481  else if (Option==5 && SubJets.size()>1 && N==2) return SubJet1.perp();
1482  else if (Option==6 && SubJets.size()>1 && N==2) return SubJet2.perp();
1483  else if (Option==7 && SubJets.size()>1 && N==2) return TMath::Sqrt(TMath::Power(SubJet1Eta-SubJet2Eta,2)+TMath::Power(DeltaPhiSubJets,2));
1484  else if (Option==8 && SubJets.size()>1 && N==2) return SubJet1LeadingTrackPt;
1485  else if (Option==9 && SubJets.size()>1 && N==2) return SubJet2LeadingTrackPt;
1486  else return -2;
1487 }
1488 
1489 
1490 
1491 Double32_t AliFJWrapper::NSubjettinessDerivativeSub(Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Double_t JetR, fastjet::PseudoJet jet, Int_t Option, Int_t Measure){ //For derivative subtraction
1492 
1493  Double_t Result=-1;
1494  std::vector<fastjet::PseudoJet> SubJet_Axes;
1495  fj::PseudoJet SubJet1_Axis;
1496  fj::PseudoJet SubJet2_Axis;
1497  if (Algorithm==0){
1498  fj::contrib::Nsubjettiness nSub(N, fj::contrib::KT_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1499  Result= nSub.result(jet);
1500  SubJet_Axes=nSub.currentAxes();
1501  }
1502  else if (Algorithm==1) {
1503  fj::contrib::Nsubjettiness nSub(N, fj::contrib::CA_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1504  Result= nSub.result(jet);
1505  SubJet_Axes=nSub.currentAxes();
1506  }
1507  else if (Algorithm==2){
1508  fj::contrib::Nsubjettiness nSub(N, fj::contrib::AntiKT_Axes(Radius), fj::contrib::NormalizedMeasure(Beta,JetR));
1509  Result= nSub.result(jet);
1510  SubJet_Axes=nSub.currentAxes();
1511  }
1512  else if (Algorithm==3) {
1513  fj::contrib::Nsubjettiness nSub(N, fj::contrib::WTA_KT_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1514  Result= nSub.result(jet);
1515  SubJet_Axes=nSub.currentAxes();
1516  }
1517  else if (Algorithm==4) {
1518  fj::contrib::Nsubjettiness nSub(N, fj::contrib::WTA_CA_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1519  Result= nSub.result(jet);
1520  SubJet_Axes=nSub.currentAxes();
1521  }
1522  else if (Algorithm==5){
1523  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_KT_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1524  Result= nSub.result(jet);
1525  SubJet_Axes=nSub.currentAxes();
1526  }
1527  else if (Algorithm==6){
1528  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_CA_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1529  Result= nSub.result(jet);
1530  SubJet_Axes=nSub.currentAxes();
1531  }
1532  else if (Algorithm==7){
1533  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_AntiKT_Axes(Radius), fj::contrib::NormalizedMeasure(Beta,JetR));
1534  Result= nSub.result(jet);
1535  SubJet_Axes=nSub.currentAxes();
1536  }
1537  else if (Algorithm==8){
1538  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_WTA_KT_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1539  Result= nSub.result(jet);
1540  SubJet_Axes=nSub.currentAxes();
1541  }
1542  else if (Algorithm==9){
1543  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_WTA_CA_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1544  Result= nSub.result(jet);
1545  SubJet_Axes=nSub.currentAxes();
1546  }
1547  else if (Algorithm==10){
1548  fj::contrib::Nsubjettiness nSub(N, fj::contrib::MultiPass_Axes(100), fj::contrib::NormalizedMeasure(Beta,JetR));
1549  Result= nSub.result(jet);
1550  SubJet_Axes=nSub.currentAxes();
1551  }
1552 
1553  SubJet1_Axis=SubJet_Axes[0];
1554  Double_t SubJet1_Eta=SubJet1_Axis.pseudorapidity();
1555  Double_t SubJet2_Eta;
1556  Double_t SubJet1_Phi=SubJet1_Axis.phi();
1557  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
1558  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
1559  Double_t SubJet2_Phi;
1560  Double_t DeltaPhi=-5;
1561  if (SubJet_Axes.size()>1){
1562  SubJet2_Axis=SubJet_Axes[1];
1563  SubJet2_Eta=SubJet2_Axis.pseudorapidity();
1564  SubJet2_Phi=SubJet2_Axis.phi();
1565  if(SubJet2_Phi < -1*TMath::Pi()) SubJet2_Phi += (2*TMath::Pi());
1566  else if (SubJet2_Phi > TMath::Pi()) SubJet2_Phi -= (2*TMath::Pi());
1567  DeltaPhi=SubJet1_Phi-SubJet2_Phi;
1568  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1569  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1570  }
1571 
1572 
1573  if (Option==0) return Result;
1574  else if (Option==1 && SubJet_Axes.size()>1 && N==2) return TMath::Sqrt(TMath::Power(SubJet1_Eta-SubJet2_Eta,2)+TMath::Power(DeltaPhi,2));
1575  else if (Option==2 && SubJet_Axes.size()>1 && N==2) return TMath::Sqrt(TMath::Power(SubJet1_Eta-SubJet2_Eta,2)+TMath::Power(DeltaPhi,2));
1576  else return -2;
1577 
1578 }
1579 
1580 
1581 
1582 
1583 
1584 
1585 
1586 
1587 
1588 #endif
AliFJWrapper::NSubjettinessDerivativeSub
Double32_t NSubjettinessDerivativeSub(Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Double_t JetR, fastjet::PseudoJet jet, Int_t Option=0, Int_t Measure=0)
AliFJWrapper::fEventSubInputVectors
std::vector< fastjet::PseudoJet > fEventSubInputVectors
Definition: AliFJWrapper.h:131
AliFJWrapper::GetJetAreaVector
fastjet::PseudoJet GetJetAreaVector(UInt_t idx) const
AliFJWrapper::GetEventSubJetAreaVector
fastjet::PseudoJet GetEventSubJetAreaVector(UInt_t idx) const
AliFJWrapper::fGRNumeratorSub
std::vector< double > fGRNumeratorSub
Definition: AliFJWrapper.h:207
AliFJWrapper::SetRMaxAndStep
void SetRMaxAndStep(Double_t rmax, Double_t dr)
Definition: AliFJWrapper.h:121
AliFJWrapper::SetUseExternalBkg
void SetUseExternalBkg(Bool_t b, Double_t rho, Double_t rhom)
Definition: AliFJWrapper.h:120
AliFJWrapper::GetFilteredJetArea
Double_t GetFilteredJetArea(UInt_t idx) const
AliFJWrapper::fInputVectors
std::vector< fastjet::PseudoJet > fInputVectors
Definition: AliFJWrapper.h:130
AliFJWrapper::fR
Double_t fR
Definition: AliFJWrapper.h:161
AliFJWrapper::Run
virtual Int_t Run()
Double_t
double Double_t
Definition: External.C:58
AliFJWrapper::Filter
virtual Int_t Filter()
title
const char * title
Definition: MakeQAPdf.C:26
AliFJWrapper::GetSubtractedJetsPts
virtual std::vector< double > GetSubtractedJetsPts(Double_t median_pt=-1, Bool_t sorted=kFALSE)
AliFJWrapper::fClustSeqSA
fastjet::ClusterSequence * fClustSeqSA
Definition: AliFJWrapper.h:152
AliFJWrapper::DoEventConstituentSubtraction
virtual Int_t DoEventConstituentSubtraction()
AliFJWrapper::DoSoftDrop
virtual Int_t DoSoftDrop()
AliFJWrapper::fPluginAlgor
Int_t fPluginAlgor
Definition: AliFJWrapper.h:167
AliFJWrapper::fClustSeqActGhosts
fastjet::ClusterSequenceActiveAreaExplicitGhosts * fClustSeqActGhosts
Definition: AliFJWrapper.h:153
AliFJWrapper::DoGenericSubtractionJetOpeningAngle_kt
virtual Int_t DoGenericSubtractionJetOpeningAngle_kt()
AliFJWrapper::fKtScatter
Double_t fKtScatter
Definition: AliFJWrapper.h:165
AliFJWrapper::fNGhostRepeats
Int_t fNGhostRepeats
Definition: AliFJWrapper.h:158
AliFJWrapper::DoGenericSubtractionJetpTD
virtual Int_t DoGenericSubtractionJetpTD()
AliFJWrapper::GetMedianUsedForBgSubtraction
Double_t GetMedianUsedForBgSubtraction() const
Definition: AliFJWrapper.h:40
AliFJWrapper::AddInputVectors
virtual void AddInputVectors(const std::vector< fastjet::PseudoJet > &vecs, Int_t offsetIndex=-99999)
AliFJWrapper::fInclusiveJets
std::vector< fastjet::PseudoJet > fInclusiveJets
Definition: AliFJWrapper.h:134
AliFJWrapper::fAlgor
fastjet::JetAlgorithm fAlgor
Definition: AliFJWrapper.h:155
AliFJWrapper::GetTitle
const char * GetTitle() const
Definition: AliFJWrapper.h:42
AliFJWrapper::fScheme
fastjet::RecombinationScheme fScheme
Definition: AliFJWrapper.h:156
AliFJWrapper::GetEventSubJetConstituents
std::vector< fastjet::PseudoJet > GetEventSubJetConstituents(UInt_t idx) const
AliFJWrapper::fMaxDelR
Double_t fMaxDelR
Definition: AliFJWrapper.h:176
AliFJWrapper::fRMax
Double_t fRMax
Definition: AliFJWrapper.h:203
AliFJWrapper::CopySettingsFrom
virtual void CopySettingsFrom(const AliFJWrapper &wrapper)
AliFJWrapper::fUseMaxDelR
Bool_t fUseMaxDelR
Definition: AliFJWrapper.h:175
AliFJWrapper::SetupAlgorithmfromOpt
void SetupAlgorithmfromOpt(const char *option)
AliFJWrapper::DoGenericSubtractionJetAngularity
virtual Int_t DoGenericSubtractionJetAngularity()
AliJetShape.h
AliFJWrapper::DoGenericSubtractionGR
virtual Int_t DoGenericSubtractionGR(Int_t ijet)
AliFJWrapper::DoGenericSubtractionJetMass
virtual Int_t DoGenericSubtractionJetMass()
AliFJWrapper::SetMinJetPt
void SetMinJetPt(Double_t MinPt)
Definition: AliFJWrapper.h:123
AliFJWrapper::SetMaxRap
void SetMaxRap(Double_t maxrap)
Definition: AliFJWrapper.h:107
AliFJWrapper::GetJetSubtractedPt
Double_t GetJetSubtractedPt(UInt_t idx) const
AliFJWrapper::fMinJetPt
Double_t fMinJetPt
Definition: AliFJWrapper.h:162
AliFJWrapper::GetFilteredJetAreaVector
fastjet::PseudoJet GetFilteredJetAreaVector(UInt_t idx) const
AliFJWrapper::SetRecombScheme
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:103
AliFJWrapper::fRange
fastjet::RangeDefinition * fRange
Definition: AliFJWrapper.h:146
AliFJWrapper::fZcut
Double_t fZcut
Definition: AliFJWrapper.h:172
AliFJWrapper::GetClusterSequenceGhosts
fastjet::ClusterSequenceActiveAreaExplicitGhosts * GetClusterSequenceGhosts() const
Definition: AliFJWrapper.h:30
AliFJWrapper::GetClusterSequenceSA
fastjet::ClusterSequence * GetClusterSequenceSA() const
Definition: AliFJWrapper.h:29
AliFJWrapper::SetStrategy
void SetStrategy(const fastjet::Strategy &strat)
Definition: AliFJWrapper.h:101
AliFJWrapper::GetEventSubJetArea
Double_t GetEventSubJetArea(UInt_t idx) const
AliFJWrapper::fFilteredJets
std::vector< fastjet::PseudoJet > fFilteredJets
Definition: AliFJWrapper.h:136
AliFJWrapper::SetupAreaTypefromOpt
void SetupAreaTypefromOpt(const char *option)
AliFJWrapper::ClearMemory
virtual void ClearMemory()
AliFJWrapper::fAreaType
fastjet::AreaType fAreaType
Definition: AliFJWrapper.h:157
FJ_includes.h
AliFJWrapper::SetR
void SetR(Double_t r)
Definition: AliFJWrapper.h:108
AliFJWrapper::fMaxRap
Double_t fMaxRap
Definition: AliFJWrapper.h:160
AliFJWrapper::GetClusterSequence
fastjet::ClusterSequenceArea * GetClusterSequence() const
Definition: AliFJWrapper.h:28
sigma
Double_t * sigma
Definition: CalculateAveragePt.C:58
AliFJWrapper::GetInclusiveJets
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
Int_t
int Int_t
Definition: External.C:63
AliFJWrapper::fClustSeqES
fastjet::ClusterSequenceArea * fClustSeqES
Definition: AliFJWrapper.h:151
AliFJWrapper::operator=
AliFJWrapper & operator=(const AliFJWrapper &wrapper)
UInt_t
unsigned int UInt_t
Definition: External.C:33
AliFJWrapper::GetJetArea
Double_t GetJetArea(UInt_t idx) const
AliFJWrapper::fSubtractedJetsPt
std::vector< double > fSubtractedJetsPt
Definition: AliFJWrapper.h:137
AliFJWrapper::fPlugin
fastjet::JetDefinition::Plugin * fPlugin
Definition: AliFJWrapper.h:144
AliFJWrapper::GetEventSubInputVectors
const std::vector< fastjet::PseudoJet > & GetEventSubInputVectors() const
Definition: AliFJWrapper.h:32
AliFJWrapper::SetLegacyMode
void SetLegacyMode(Bool_t mode)
Definition: AliFJWrapper.h:118
AliFJWrapper::fDoFilterArea
Bool_t fDoFilterArea
Definition: AliFJWrapper.h:198
AliFJWrapper::SetAlgorithm
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:102
AliFJWrapper::DoGenericSubtractionJetCircularity
virtual Int_t DoGenericSubtractionJetCircularity()
AliFJWrapper::SetupSchemefromOpt
void SetupSchemefromOpt(const char *option)
AliFJWrapper::GetInputGhosts
const std::vector< fastjet::PseudoJet > & GetInputGhosts() const
Definition: AliFJWrapper.h:33
AliFJWrapper::fGroomedJets
std::vector< fastjet::PseudoJet > fGroomedJets
Definition: AliFJWrapper.h:139
AliFJWrapper::SetUseArea4Vector
void SetUseArea4Vector(Bool_t useA4v)
Definition: AliFJWrapper.h:113
AliFJWrapper::fClustSeq
fastjet::ClusterSequenceArea * fClustSeq
Definition: AliFJWrapper.h:150
AliFJWrapper::SetName
void SetName(const char *name)
Definition: AliFJWrapper.h:99
AliFJWrapper::fBeta
Double_t fBeta
Definition: AliFJWrapper.h:173
Option
Definition: Option.C:68
AliFJWrapper::SubtractBackground
virtual void SubtractBackground(const Double_t median_pt=-1)
AliFJWrapper::fInputGhosts
std::vector< fastjet::PseudoJet > fInputGhosts
Definition: AliFJWrapper.h:133
AliFJWrapper::fGridScatter
Double_t fGridScatter
Definition: AliFJWrapper.h:164
AliFJWrapper::Clear
virtual void Clear(const Option_t *="")
AliFJWrapper::fLegacyMode
Bool_t fLegacyMode
Definition: AliFJWrapper.h:199
AliFJWrapper::fVorAreaSpec
fastjet::VoronoiAreaSpec * fVorAreaSpec
Definition: AliFJWrapper.h:141
mode
Int_t mode
Definition: anaM.C:41
AliFJWrapper::fDRStep
Double_t fDRStep
Definition: AliFJWrapper.h:204
AliFJWrapper::fGRDenominatorSub
std::vector< double > fGRDenominatorSub
Definition: AliFJWrapper.h:208
AliFJWrapper::GetLegacyMode
Bool_t GetLegacyMode()
Definition: AliFJWrapper.h:51
AliFJWrapper::fGhostedAreaSpec
fastjet::GhostedAreaSpec * fGhostedAreaSpec
Definition: AliFJWrapper.h:142
AliFJWrapper::fUseArea4Vector
Bool_t fUseArea4Vector
Definition: AliFJWrapper.h:170
AliFJWrapper::SetNRepeats
void SetNRepeats(Int_t nrepeat)
Definition: AliFJWrapper.h:105
AliFJWrapper::GetEventSubJets
const std::vector< fastjet::PseudoJet > & GetEventSubJets() const
Definition: AliFJWrapper.h:35
AliFJWrapper::GetGRNumeratorSub
virtual std::vector< double > GetGRNumeratorSub() const
Definition: AliFJWrapper.h:76
AliFJWrapper::~AliFJWrapper
virtual ~AliFJWrapper()
AliFJWrapper::fGRDenominator
std::vector< double > fGRDenominator
Definition: AliFJWrapper.h:206
AliFJWrapper::DoGenericSubtractionJet1subjettiness_kt
virtual Int_t DoGenericSubtractionJet1subjettiness_kt()
AliFJWrapper::fGhostArea
Double_t fGhostArea
Definition: AliFJWrapper.h:159
AliFJWrapper::GetMedianAndSigma
virtual void GetMedianAndSigma(Double_t &median, Double_t &sigma, Int_t remove=0) const
AliFJWrapper::fName
TString fName
Definition: AliFJWrapper.h:128
AliFJWrapper::SetGhostArea
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:106
AliFJWrapper::fMedUsedForBgSub
Double_t fMedUsedForBgSub
Definition: AliFJWrapper.h:169
AliFJWrapper::GetJetConstituents
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
AliFJWrapper::fJetDef
fastjet::JetDefinition * fJetDef
Definition: AliFJWrapper.h:143
AliFJWrapper::DoGenericSubtractionJetSigma2
virtual Int_t DoGenericSubtractionJetSigma2()
AliFJWrapper::SetupStrategyfromOpt
void SetupStrategyfromOpt(const char *option)
AliFJWrapper::GetInputVectors
const std::vector< fastjet::PseudoJet > & GetInputVectors() const
Definition: AliFJWrapper.h:31
AliFJWrapper::fAreaDef
fastjet::AreaDefinition * fAreaDef
Definition: AliFJWrapper.h:140
AliFJWrapper::GetGRDenominatorSub
virtual std::vector< double > GetGRDenominatorSub() const
Definition: AliFJWrapper.h:77
AliFJWrapper::SetMaxDelR
void SetMaxDelR(Double_t r)
Definition: AliFJWrapper.h:125
AliFJWrapper::fUseExternalBkg
Bool_t fUseExternalBkg
Definition: AliFJWrapper.h:200
AliFJWrapper::SetTitle
void SetTitle(const char *title)
Definition: AliFJWrapper.h:100
AliFJWrapper::NSubjettiness
Double_t NSubjettiness(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)
AliFJWrapper::SetPluginAlgor
void SetPluginAlgor(Int_t plugin)
Definition: AliFJWrapper.h:112
AliFJWrapper::DoConstituentSubtraction
virtual Int_t DoConstituentSubtraction()
AliFJWrapper::fRhom
Double_t fRhom
Definition: AliFJWrapper.h:202
AliFJWrapper::GetGRDenominator
virtual std::vector< double > GetGRDenominator() const
Definition: AliFJWrapper.h:75
AliFJWrapper::SetRhoRhom
void SetRhoRhom(Double_t rho, Double_t rhom)
Definition: AliFJWrapper.h:122
AliFJWrapper::SetGridScatter
void SetGridScatter(Double_t gridSc)
Definition: AliFJWrapper.h:109
AliFJWrapper::fConstituentSubtrJets
std::vector< fastjet::PseudoJet > fConstituentSubtrJets
Definition: AliFJWrapper.h:138
AliFJWrapper::fEventSubJets
std::vector< fastjet::PseudoJet > fEventSubJets
Definition: AliFJWrapper.h:135
AliFJWrapper::fEventSubCorrectedVectors
std::vector< fastjet::PseudoJet > fEventSubCorrectedVectors
Definition: AliFJWrapper.h:132
AliFJWrapper::SetMeanGhostKt
void SetMeanGhostKt(Double_t meankt)
Definition: AliFJWrapper.h:111
Option_t
const char Option_t
Definition: External.C:48
AliFJWrapper::AddInputVector
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
AliFJWrapper::AddInputGhost
virtual void AddInputGhost(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
AliFJWrapper::DoGenericSubtractionJet3subjettiness_kt
virtual Int_t DoGenericSubtractionJet3subjettiness_kt()
Bool_t
bool Bool_t
Definition: External.C:53
AliFJWrapper::fEventSub
Bool_t fEventSub
Definition: AliFJWrapper.h:174
AliFJWrapper::SetLegacyFJ
void SetLegacyFJ()
AliFJWrapper::SetEventSub
void SetEventSub(Bool_t b)
Definition: AliFJWrapper.h:124
AliFJWrapper::GetFilteredJetConstituents
std::vector< fastjet::PseudoJet > GetFilteredJetConstituents(UInt_t idx) const
AliFJWrapper::SetKtScatter
void SetKtScatter(Double_t ktSc)
Definition: AliFJWrapper.h:110
AliFJWrapper::GetGRNumerator
virtual std::vector< double > GetGRNumerator() const
Definition: AliFJWrapper.h:74
AliFJWrapper::ClassName
virtual const char * ClassName() const
Definition: AliFJWrapper.h:23
AliFJWrapper::RemoveLastInputVector
virtual void RemoveLastInputVector()
AliFJWrapper::SetAreaType
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:104
AliFJWrapper::fRho
Double_t fRho
Definition: AliFJWrapper.h:201
AliFJWrapper::GetDoFilterArea
Bool_t GetDoFilterArea()
Definition: AliFJWrapper.h:52
AliFJWrapper::DoGenericSubtractionJetLeSub
virtual Int_t DoGenericSubtractionJetLeSub()
AliFJWrapper::DoGenericSubtractionJetConstituent
virtual Int_t DoGenericSubtractionJetConstituent()
AliFJWrapper::GetName
const char * GetName() const
Definition: AliFJWrapper.h:41
AliFJWrapper::GetFilteredJets
const std::vector< fastjet::PseudoJet > & GetFilteredJets() const
Definition: AliFJWrapper.h:36
AliFJWrapper::DoGenericSubtractionJet2subjettiness_kt
virtual Int_t DoGenericSubtractionJet2subjettiness_kt()
AliFJWrapper::fGRNumerator
std::vector< double > fGRNumerator
Definition: AliFJWrapper.h:205
AliFJWrapper::fMeanGhostKt
Double_t fMeanGhostKt
Definition: AliFJWrapper.h:166
AliFJWrapper::fStrategy
fastjet::Strategy fStrategy
Definition: AliFJWrapper.h:154
AliFJWrapper::fTitle
TString fTitle
Definition: AliFJWrapper.h:129