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