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