AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliFJWrapper.h
Go to the documentation of this file.
1 #ifndef AliFJWrapper_H
2 #define AliFJWrapper_H
3 
4 #if !defined(__CINT__)
5 
6 #include <vector>
7 #include <TString.h>
8 #include "AliLog.h"
9 #include "FJ_includes.h"
10 #include "AliJetShape.h"
11 
12 
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_casd() const {return fGenSubtractorInfoJet1subjettiness_casd ; }
74  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJet2subjettiness_casd() const {return fGenSubtractorInfoJet2subjettiness_casd ; }
75  const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetOpeningAngle_casd() const {return fGenSubtractorInfoJetOpeningAngle_casd ; }
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) {fUseMaxDelR = kTRUE; fMaxDelR = r;}
145 
146  protected:
149  std::vector<fastjet::PseudoJet> fInputVectors;
150  std::vector<fastjet::PseudoJet> fEventSubInputVectors;
151  std::vector<fastjet::PseudoJet> fEventSubCorrectedVectors;
152  std::vector<fastjet::PseudoJet> fInputGhosts;
153  std::vector<fastjet::PseudoJet> fInclusiveJets;
154  std::vector<fastjet::PseudoJet> fEventSubJets;
155  std::vector<fastjet::PseudoJet> fFilteredJets;
156  std::vector<double> fSubtractedJetsPt;
157  std::vector<fastjet::PseudoJet> fConstituentSubtrJets;
158  std::vector<fastjet::PseudoJet> fGroomedJets;
159  fastjet::AreaDefinition *fAreaDef;
160  fastjet::VoronoiAreaSpec *fVorAreaSpec;
161  fastjet::GhostedAreaSpec *fGhostedAreaSpec;
162  fastjet::JetDefinition *fJetDef;
163  fastjet::JetDefinition::Plugin *fPlugin;
164 #ifndef FASTJET_VERSION
165  fastjet::RangeDefinition *fRange;
166 #else
167  fastjet::Selector *fRange;
168 #endif
169  fastjet::ClusterSequenceArea *fClustSeq;
170  fastjet::ClusterSequenceArea *fClustSeqES;
171  fastjet::ClusterSequence *fClustSeqSA;
172  fastjet::ClusterSequenceActiveAreaExplicitGhosts *fClustSeqActGhosts;
173  fastjet::Strategy fStrategy;
174  fastjet::JetAlgorithm fAlgor;
175  fastjet::RecombinationScheme fScheme;
176  fastjet::AreaType fAreaType;
182  // no setters for the moment - used default values in the constructor
187  // extra parameters
190  // condition to stop the grooming (rejection of soft splitting) z > fZcut theta^fBeta
191  Double_t fZcut; // fZcut = 0.1
192  Double_t fBeta; // fBeta = 0
196 #ifdef FASTJET_VERSION
197  fastjet::JetMedianBackgroundEstimator *fBkrdEstimator;
198  //from contrib package
199  fastjet::contrib::GenericSubtractor *fGenSubtractor;
200  fastjet::contrib::ConstituentSubtractor *fConstituentSubtractor;
201  fastjet::contrib::ConstituentSubtractor *fEventConstituentSubtractor;
202  fastjet::contrib::SoftDrop *fSoftDrop;
203  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetMass;
204  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoGRNum;
205  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoGRDen;
206  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetAngularity;
207  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetpTD;
208  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetCircularity;
209  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetSigma2;
210  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetConstituent;
211  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetLeSub;
212  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet1subjettiness_kt;
213  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet2subjettiness_kt;
214  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet3subjettiness_kt;
215  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetOpeningAngle_kt;
216  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet1subjettiness_ca;
217  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet2subjettiness_ca;
218  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetOpeningAngle_ca;
219  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet1subjettiness_akt02;
220  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet2subjettiness_akt02;
221  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetOpeningAngle_akt02;
222  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet1subjettiness_casd;
223  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJet2subjettiness_casd;
224  std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetOpeningAngle_casd;
225 #endif
229  Double_t fRho; // pT background density
230  Double_t fRhom; // mT background density
233  std::vector<double> fGRNumerator;
234  std::vector<double> fGRDenominator;
235  std::vector<double> fGRNumeratorSub;
236  std::vector<double> fGRDenominatorSub;
237 
238  virtual void SubtractBackground(const Double_t median_pt = -1);
239 
240  private:
241  AliFJWrapper();
242  AliFJWrapper(const AliFJWrapper& wrapper);
243  AliFJWrapper& operator = (const AliFJWrapper& wrapper);
244 };
245 #endif
246 #endif
247 
248 #ifdef AliFJWrapper_CXX
249 #undef AliFJWrapper_CXX
250 
251 #if defined __GNUC__
252 #pragma GCC system_header
253 #endif
254 
255 namespace fj = fastjet;
256 
257 //_________________________________________________________________________________________________
258 AliFJWrapper::AliFJWrapper(const char *name, const char *title)
259  :
260  fName (name)
261  , fTitle (title)
262  , fInputVectors ( )
263  , fEventSubInputVectors ( )
264  , fEventSubCorrectedVectors ( )
265  , fInputGhosts ( )
266  , fInclusiveJets ( )
267  , fEventSubJets ( )
268  , fFilteredJets ( )
269  , fSubtractedJetsPt ( )
270  , fConstituentSubtrJets ( )
271  , fSoftDrop ( )
272  , fAreaDef (0)
273  , fVorAreaSpec (0)
274  , fGhostedAreaSpec (0)
275  , fJetDef (0)
276  , fPlugin (0)
277  , fRange (0)
278  , fClustSeq (0)
279  , fClustSeqES (0)
280  , fClustSeqSA (0)
281  , fClustSeqActGhosts (0)
282  , fStrategy (fj::Best)
283  , fAlgor (fj::kt_algorithm)
284  , fScheme (fj::BIpt_scheme)
285  , fAreaType (fj::active_area)
286  , fNGhostRepeats (1)
287  , fGhostArea (0.005)
288  , fMaxRap (1.)
289  , fR (0.4)
290  , fGridScatter (1.0)
291  , fKtScatter (0.1)
292  , fMeanGhostKt (1e-100)
293  , fPluginAlgor (0)
294  , fMedUsedForBgSub (0)
295  , fUseArea4Vector (kFALSE)
296  , fZcut(0.1)
297  , fBeta(0)
298  , fEventSub (kFALSE)
299  , fUseMaxDelR (kFALSE)
300  , fMaxDelR (0.4)
301 #ifdef FASTJET_VERSION
302  , fBkrdEstimator (0)
303  , fGenSubtractor (0)
304  , fConstituentSubtractor (0)
305  , fEventConstituentSubtractor (0)
306  , fGenSubtractorInfoJetMass ( )
307  , fGenSubtractorInfoGRNum ( )
308  , fGenSubtractorInfoGRDen ( )
309  , fGenSubtractorInfoJetAngularity ( )
310  , fGenSubtractorInfoJetpTD ( )
311  , fGenSubtractorInfoJetCircularity( )
312  , fGenSubtractorInfoJetSigma2()
313  , fGenSubtractorInfoJetConstituent ( )
314  , fGenSubtractorInfoJetLeSub ( )
315  , fGenSubtractorInfoJet1subjettiness_kt ( )
316  , fGenSubtractorInfoJet2subjettiness_kt ( )
317  , fGenSubtractorInfoJet3subjettiness_kt ( )
318  , fGenSubtractorInfoJetOpeningAngle_kt ( )
319  , fGenSubtractorInfoJet1subjettiness_ca ( )
320  , fGenSubtractorInfoJet2subjettiness_ca ( )
321  , fGenSubtractorInfoJetOpeningAngle_ca ( )
322  , fGenSubtractorInfoJet1subjettiness_akt02 ( )
323  , fGenSubtractorInfoJet2subjettiness_akt02 ( )
324  , fGenSubtractorInfoJetOpeningAngle_akt02 ( )
325  , fGenSubtractorInfoJet1subjettiness_casd ( )
326  , fGenSubtractorInfoJet2subjettiness_casd ( )
327  , fGenSubtractorInfoJetOpeningAngle_casd ( )
328 
329 #endif
330  , fDoFilterArea (false)
331  , fLegacyMode (false)
332  , fUseExternalBkg (false)
333  , fRho (0)
334  , fRhom (0)
335  , fRMax(2.)
336  , fDRStep(0.04)
337  , fGRNumerator()
338  , fGRDenominator()
339  , fGRNumeratorSub()
340  , fGRDenominatorSub()
341 {
342  // Constructor.
343 }
344 
345 //_________________________________________________________________________________________________
347 {
348  // Destructor.
349  ClearMemory();
350 }
351 
352 //_________________________________________________________________________________________________
354 {
355  // Destructor.
356  if (fAreaDef) { delete fAreaDef; fAreaDef = NULL; }
357  if (fVorAreaSpec) { delete fVorAreaSpec; fVorAreaSpec = NULL; }
358  if (fGhostedAreaSpec) { delete fGhostedAreaSpec; fGhostedAreaSpec = NULL; }
359  if (fJetDef) { delete fJetDef; fJetDef = NULL; }
360  if (fPlugin) { delete fPlugin; fPlugin = NULL; }
361  if (fRange) { delete fRange; fRange = NULL; }
362  if (fClustSeq) { delete fClustSeq; fClustSeq = NULL; }
363  if (fClustSeqES) { delete fClustSeqES; fClustSeqES = NULL; }
364  if (fClustSeqSA) { delete fClustSeqSA; fClustSeqSA = NULL; }
366  #ifdef FASTJET_VERSION
367  if (fBkrdEstimator) { delete fBkrdEstimator; fBkrdEstimator = NULL; }
368  if (fGenSubtractor) { delete fGenSubtractor; fGenSubtractor = NULL; }
369  if (fConstituentSubtractor) { delete fConstituentSubtractor; fConstituentSubtractor = NULL; }
370  if (fSoftDrop) { delete fSoftDrop; fSoftDrop = NULL;}
371  #endif
372 }
373 
374 //_________________________________________________________________________________________________
375 void AliFJWrapper::CopySettingsFrom(const AliFJWrapper& wrapper)
376 {
377  // Copy some settings.
378  // You very often want to keep most of the settings
379  // but change only the algorithm or R - do it after call to this function
380 
381  fStrategy = wrapper.fStrategy;
382  fAlgor = wrapper.fAlgor;
383  fScheme = wrapper.fScheme;
384  fAreaType = wrapper.fAreaType;
385  fNGhostRepeats = wrapper.fNGhostRepeats;
386  fGhostArea = wrapper.fGhostArea;
387  fMaxRap = wrapper.fMaxRap;
388  fR = wrapper.fR;
389  fGridScatter = wrapper.fGridScatter;
390  fKtScatter = wrapper.fKtScatter;
391  fMeanGhostKt = wrapper.fMeanGhostKt;
392  fPluginAlgor = wrapper.fPluginAlgor;
394  fZcut = wrapper.fZcut;
395  fBeta = wrapper.fBeta;
396  fLegacyMode = wrapper.fLegacyMode;
398  fRho = wrapper.fRho;
399  fRhom = wrapper.fRhom;
400 }
401 
402 //_________________________________________________________________________________________________
403 void AliFJWrapper::Clear(const Option_t */*opt*/)
404 {
405  // Simply clear the input vectors.
406  // Make sure done on every event if the instance is reused
407  // Reset the median to zero.
408 
409  fInputVectors.clear();
410  fEventSubInputVectors.clear();
411  fInputGhosts.clear();
412  fMedUsedForBgSub = 0;
413 
414  // for the moment brute force delete everything
415  ClearMemory();
416 }
417 
418 //_________________________________________________________________________________________________
420 {
421  // Remove last input vector
422 
423  fInputVectors.pop_back();
424  fEventSubInputVectors.pop_back();
425 }
426 
427 //_________________________________________________________________________________________________
429 {
430  // Make the input pseudojet.
431 
432  fastjet::PseudoJet inVec(px, py, pz, E);
433 
434  // Salvatore Aiola: not sure why this was done...
435  //if (index > -99999) {
436  inVec.set_user_index(index);
437  //} else {
438  //inVec.set_user_index(fInputVectors.size());
439  //}
440 
441  // add to the fj container of input vectors
442  fInputVectors.push_back(inVec);
443  if(fEventSub) fEventSubInputVectors.push_back(inVec);
444 
445 }
446 
447 //_________________________________________________________________________________________________
448 void AliFJWrapper::AddInputVector(const fj::PseudoJet& vec, Int_t index)
449 {
450  // Add an input pseudojet.
451 
452  fj::PseudoJet inVec = vec;
453 
454  // Salvatore Aiola: not sure why this was done...
456  inVec.set_user_index(index);
457  //} else {
458  //inVec.set_user_index(fInputVectors.size());
459  //}
460 
461  // add to the fj container of input vectors
462  fInputVectors.push_back(inVec);
463  //if(fEventSub) fEventSubInputVectors.push_back(inVec);
464 }
465 
466 //_________________________________________________________________________________________________
467 void AliFJWrapper::AddInputVectors(const std::vector<fj::PseudoJet>& vecs, Int_t offsetIndex)
468 {
469  // Add the input from vector of pseudojets.
470 
471  for (UInt_t i = 0; i < vecs.size(); ++i) {
472  fj::PseudoJet inVec = vecs[i];
473  if (offsetIndex > -99999)
474  inVec.set_user_index(fInputVectors.size() + offsetIndex);
475  // add to the fj container of input vectors
476  fInputVectors.push_back(inVec);
477  }
478  //Is this correct?
479  /* if(fEventSub){
480  for (UInt_t i = 0; i < vecs.size(); ++i) {
481  fj::PseudoJet inVec = vecs[i];
482  if (offsetIndex > -99999)
483  inVec.set_user_index(fEventSubInputVectors.size() + offsetIndex);
484  // add to the fj container of input vectors
485  fEventSubInputVectors.push_back(inVec);
486  }
487  }*/
488 }
489 
490 //_________________________________________________________________________________________________
492 {
493  // Make the input pseudojet.
494 
495  fastjet::PseudoJet inVec(px, py, pz, E);
496 
497  if (index > -99999) {
498  inVec.set_user_index(index);
499  } else {
500  inVec.set_user_index(fInputGhosts.size());
501  }
502 
503  // add to the fj container of input vectors
504  fInputGhosts.push_back(inVec);
505  if (!fDoFilterArea) fDoFilterArea = kTRUE;
506 }
507 
508 //_________________________________________________________________________________________________
510 {
511  // Get the jet area.
512 
513  Double_t retval = -1; // really wrong area..
514  if ( idx < fInclusiveJets.size() ) {
515  retval = fClustSeq->area(fInclusiveJets[idx]);
516  } else {
517  AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
518  }
519  return retval;
520 }
521 
522 //_________________________________________________________________________________________________
524 {
525  // Get the jet area.
526 
527  Double_t retval = -1; // really wrong area..
528  if ( idx < fEventSubJets.size() ) {
529  retval = fClustSeqES->area(fEventSubJets[idx]);
530  } else {
531  AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
532  }
533  return retval;
534 }
535 //_________________________________________________________________________________________________
537 {
538  // Get the filtered jet area.
539 
540  Double_t retval = -1; // really wrong area..
541  if (fDoFilterArea && fClustSeqActGhosts && (idx<fFilteredJets.size())) {
542  retval = fClustSeqActGhosts->area(fFilteredJets[idx]);
543  } else {
544  AliError(Form("[e] ::GetFilteredJetArea wrong index: %d",idx));
545  }
546  return retval;
547 }
548 
549 //_________________________________________________________________________________________________
550 fastjet::PseudoJet AliFJWrapper::GetJetAreaVector(UInt_t idx) const
551 {
552  // Get the jet area as vector.
553  fastjet::PseudoJet retval;
554  if ( idx < fInclusiveJets.size() ) {
555  retval = fClustSeq->area_4vector(fInclusiveJets[idx]);
556  } else {
557  AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
558  }
559  return retval;
560 }
561 
562 //_________________________________________________________________________________________________
563 fastjet::PseudoJet AliFJWrapper::GetEventSubJetAreaVector(UInt_t idx) const
564 {
565  // Get the jet area as vector.
566  fastjet::PseudoJet retval;
567  if ( idx < fEventSubJets.size() ) {
568  retval = fClustSeqES->area_4vector(fEventSubJets[idx]);
569  } else {
570  AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
571  }
572  return retval;
573 }
574 
575 //_________________________________________________________________________________________________
576 fastjet::PseudoJet AliFJWrapper::GetFilteredJetAreaVector(UInt_t idx) const
577 {
578  // Get the jet area as vector.
579  fastjet::PseudoJet retval;
580  if (fDoFilterArea && fClustSeqActGhosts && (idx<fFilteredJets.size())) {
581  retval = fClustSeqActGhosts->area_4vector(fFilteredJets[idx]);
582  } else {
583  AliError(Form("[e] ::GetFilteredJetArea wrong index: %d",idx));
584  }
585  return retval;
586 }
587 
588 //_________________________________________________________________________________________________
589 std::vector<double> AliFJWrapper::GetSubtractedJetsPts(Double_t median_pt, Bool_t sorted)
590 {
591  // Get subtracted jets pTs, returns vector.
592 
593  SubtractBackground(median_pt);
594 
595  if (kTRUE == sorted) {
596  std::sort(fSubtractedJetsPt.begin(), fSubtractedJetsPt.begin());
597  }
598  return fSubtractedJetsPt;
599 }
600 
601 //_________________________________________________________________________________________________
603 {
604  // Get subtracted jets pTs, returns Double_t.
605 
606  Double_t retval = -99999.; // really wrong pt..
607  if ( idx < fSubtractedJetsPt.size() ) {
608  retval = fSubtractedJetsPt[idx];
609  }
610  return retval;
611 }
612 
613 //_________________________________________________________________________________________________
614 std::vector<fastjet::PseudoJet>
616 {
617  // Get jets constituents.
618 
619  std::vector<fastjet::PseudoJet> retval;
620 
621  if ( idx < fInclusiveJets.size() ) {
622  retval = fClustSeq->constituents(fInclusiveJets[idx]);
623  } else {
624  AliError(Form("[e] ::GetJetConstituents wrong index: %d",idx));
625  }
626 
627  return retval;
628 }
629 
630 //_________________________________________________________________________________________________
631 std::vector<fastjet::PseudoJet>
633 {
634  // Get jets constituents.
635 
636  std::vector<fastjet::PseudoJet> retval;
637 
638  if ( idx < fEventSubJets.size() ) {
639  retval = fClustSeqES->constituents(fEventSubJets[idx]);
640  } else {
641  AliError(Form("[e] ::GetJetConstituents wrong index: %d",idx));
642  }
643 
644  return retval;
645 }
646 
647 //_________________________________________________________________________________________________
648 std::vector<fastjet::PseudoJet>
650 {
651  // Get jets constituents.
652 
653  std::vector<fastjet::PseudoJet> retval;
654 
655  if ( idx < fFilteredJets.size() ) {
656  if (fClustSeqSA) retval = fClustSeqSA->constituents(fFilteredJets[idx]);
657  if (fClustSeqActGhosts) retval = fClustSeqActGhosts->constituents(fFilteredJets[idx]);
658  } else {
659  AliError(Form("[e] ::GetFilteredJetConstituents wrong index: %d",idx));
660  }
661 
662  return retval;
663 }
664 
665 //_________________________________________________________________________________________________
666 void AliFJWrapper::GetMedianAndSigma(Double_t &median, Double_t &sigma, Int_t remove) const
667 {
668  // Get the median and sigma from fastjet.
669  // User can also do it on his own because the cluster sequence is exposed (via a getter)
670 
671  if (!fClustSeq) {
672  AliError("[e] Run the jfinder first.");
673  return;
674  }
675 
676  Double_t mean_area = 0;
677  try {
678  if(0 == remove) {
679  fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
680  } else {
681  std::vector<fastjet::PseudoJet> input_jets = sorted_by_pt(fClustSeq->inclusive_jets());
682  input_jets.erase(input_jets.begin(), input_jets.begin() + remove);
683  fClustSeq->get_median_rho_and_sigma(input_jets, *fRange, fUseArea4Vector, median, sigma, mean_area);
684  input_jets.clear();
685  }
686  } catch (fj::Error) {
687  AliError(" [w] FJ Exception caught.");
688  median = -1.;
689  sigma = -1;
690  }
691 }
692 
693 //_________________________________________________________________________________________________
695 {
696  // Run the actual jet finder.
697 
698  if (fAreaType == fj::voronoi_area) {
699  // Rfact - check dependence - default is 1.
700  // NOTE: hardcoded variable!
701  fVorAreaSpec = new fj::VoronoiAreaSpec(1.);
702  fAreaDef = new fj::AreaDefinition(*fVorAreaSpec);
703  } else {
704  fGhostedAreaSpec = new fj::GhostedAreaSpec(fMaxRap,
706  fGhostArea,
707  fGridScatter,
708  fKtScatter,
709  fMeanGhostKt);
710 
711  fAreaDef = new fj::AreaDefinition(*fGhostedAreaSpec, fAreaType);
712  }
713 
714  // this is acceptable by fastjet:
715 #ifndef FASTJET_VERSION
716  fRange = new fj::RangeDefinition(fMaxRap - 0.95 * fR);
717 #else
718  fRange = new fj::Selector(fj::SelectorAbsRapMax(fMaxRap - 0.95 * fR));
719 #endif
720 
721  if (fAlgor == fj::plugin_algorithm) {
722  if (fPluginAlgor == 0) {
723  // SIS CONE ALGOR
724  // NOTE: hardcoded split parameter
725  Double_t overlap_threshold = 0.75; // NOTE: this actually splits a lot: thr/min(pt1,pt2)
726  fPlugin = new fj::SISConePlugin(fR,
727  overlap_threshold,
728  0, //search of stable cones - zero = until no more
729  1.0); // this should be seed effectively for proto jets
730  fJetDef = new fastjet::JetDefinition(fPlugin);
731  } else if (fPluginAlgor == 1) {
732  // CDF cone
733  // NOTE: hardcoded split parameter
734  Double_t overlap_threshold = 0.75; // NOTE: this actually splits a lot: thr/min(pt1,pt2)
735  fPlugin = new fj::CDFMidPointPlugin(fR,
736  overlap_threshold,
737  1.0, //search of stable cones - zero = until no more
738  1.0); // this should be seed effectively for proto jets
739  fJetDef = new fastjet::JetDefinition(fPlugin);
740  } else {
741  AliError("[e] Unrecognized plugin number!");
742  }
743  } else {
744  fJetDef = new fj::JetDefinition(fAlgor, fR, fScheme, fStrategy);
745  }
746 
747  try {
748  fClustSeq = new fj::ClusterSequenceArea(fInputVectors, *fJetDef, *fAreaDef);
749  if(fEventSub){
751  fClustSeqES = new fj::ClusterSequenceArea(fEventSubCorrectedVectors, *fJetDef, *fAreaDef);
752  }
753  } catch (fj::Error) {
754  AliError(" [w] FJ Exception caught.");
755  return -1;
756  }
757 
758  // FJ3 :: Define an JetMedianBackgroundEstimator just in case it will be used
759 #ifdef FASTJET_VERSION
760  fBkrdEstimator = new fj::JetMedianBackgroundEstimator(fj::SelectorAbsRapMax(fMaxRap));
761 #endif
762 
763  if (fLegacyMode) { SetLegacyFJ(); } // for FJ 2.x even if fLegacyMode is set, SetLegacyFJ is dummy
764 
765  // inclusive jets:
766  fInclusiveJets.clear();
767  fEventSubJets.clear();
768  fInclusiveJets = fClustSeq->inclusive_jets(0.0);
769  if(fEventSub) fEventSubJets = fClustSeqES->inclusive_jets(0.0);
770 
771  return 0;
772 }
773 
774 //_________________________________________________________________________________________________
776 {
777 //
778 // AliFJWrapper::Filter
779 //
780 
781  fJetDef = new fj::JetDefinition(fAlgor, fR, fScheme, fStrategy);
782 
783  if (fDoFilterArea) {
784  if (fInputGhosts.size()>0) {
785  try {
786  fClustSeqActGhosts = new fj::ClusterSequenceActiveAreaExplicitGhosts(fInputVectors,
787  *fJetDef,
788  fInputGhosts,
789  fGhostArea);
790  } catch (fj::Error) {
791  AliError(" [w] FJ Exception caught.");
792  return -1;
793  }
794 
795  fFilteredJets.clear();
796  fFilteredJets = fClustSeqActGhosts->inclusive_jets(0.0);
797  } else {
798  return -1;
799  }
800  } else {
801  try {
802  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
803  } catch (fj::Error) {
804  AliError(" [w] FJ Exception caught.");
805  return -1;
806  }
807 
808  fFilteredJets.clear();
809  fFilteredJets = fClustSeqSA->inclusive_jets(0.0);
810  }
811 
812  return 0;
813 }
814 
815 //_________________________________________________________________________________________________
817 {
818  // This methods enable legacy behaviour (FastJet 2.x) when AliROOT is compiled with FastJet 3.x
819 #ifdef FASTJET_VERSION
820  std::cout << "WARNING! Setting FastJet in legacy mode" << std::endl;
821  if (fGhostedAreaSpec) { fGhostedAreaSpec->set_fj2_placement(kTRUE); }
822  if (fBkrdEstimator) {
823  fBkrdEstimator->set_provide_fj2_sigma(kTRUE);
824  fBkrdEstimator->set_use_area_4vector(kFALSE);
825  }
826 #endif
827 }
828 
829 //_________________________________________________________________________________________________
831 {
832  // Subtract the background (specify the value - see below the meaning).
833  // Negative argument means the bg will be determined with the current algorithm
834  // this is the default behavior. Zero is allowed
835  // Note: user may set the switch for area4vector based subtraction.
836 
837  Double_t median = 0;
838  Double_t sigma = 0;
839  Double_t mean_area = 0;
840 
841  // clear the subtracted jet pt's vector<double>
842  fSubtractedJetsPt.clear();
843 
844  // check what was specified (default is -1)
845  if (median_pt < 0) {
846  try {
847  fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
848  }
849 
850  catch (fj::Error) {
851  AliError(" [w] FJ Exception caught.");
852  median = -9999.;
853  sigma = -1;
854  fMedUsedForBgSub = median;
855  return;
856  }
857  fMedUsedForBgSub = median;
858  } else {
859  // we do not know the sigma in this case
860  sigma = -1;
861  if (0.0 == median_pt) {
862  // AliWarning(Form(" [w] Median specified for bg subtraction is ZERO: %f", median_pt )) << endl;
863  fMedUsedForBgSub = 0.;
864  } else {
865  fMedUsedForBgSub = median_pt;
866  }
867  }
868 
869  // subtract:
870  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
871  if ( fUseArea4Vector ) {
872  // subtract the background using the area4vector
873  fj::PseudoJet area4v = fClustSeq->area_4vector(fInclusiveJets[i]);
874  fj::PseudoJet jet_sub = fInclusiveJets[i] - area4v * fMedUsedForBgSub;
875  fSubtractedJetsPt.push_back(jet_sub.perp()); // here we put only the pt of the jet - note: this can be negative
876  } else {
877  // subtract the background using scalars
878  // fj::PseudoJet jet_sub = fInclusiveJets[i] - area * fMedUsedForBgSub_;
879  Double_t area = fClustSeq->area(fInclusiveJets[i]);
880  // standard subtraction
881  Double_t pt_sub = fInclusiveJets[i].perp() - fMedUsedForBgSub * area;
882  fSubtractedJetsPt.push_back(pt_sub); // here we put only the pt of the jet - note: this can be negative
883  }
884  }
885 }
886 
887 //_________________________________________________________________________________________________
888 void AliFJWrapper::DoGenericSubtraction(const fastjet::FunctionOfPseudoJet<Double32_t>& jetshape, std::vector<fastjet::contrib::GenericSubtractorInfo>& output) {
889  //Do generic subtraction for 1subjettiness
890 #ifdef FASTJET_VERSION
891  CreateGenSub();
892 
893  // clear the generic subtractor info vector
894  output.clear();
895  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
896  fj::contrib::GenericSubtractorInfo info_jetshape;
897  if(fInclusiveJets[i].perp()>1.e-4)
898  double subtracted_shape = (*fGenSubtractor)(jetshape, fInclusiveJets[i], info_jetshape);
899  output.push_back(info_jetshape);
900  }
901 #endif
902 }
903 
904 //_________________________________________________________________________________________________
906  //Do generic subtraction for jet mass
907 #ifdef FASTJET_VERSION
908  CreateGenSub();
909 
910  // Define jet shape
911  AliJetShapeMass shapeMass;
912 
913  // clear the generic subtractor info vector
914  fGenSubtractorInfoJetMass.clear();
915  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
916  fj::contrib::GenericSubtractorInfo info;
917  if(fInclusiveJets[i].perp()>1.e-4)
918  double subtracted_shape = (*fGenSubtractor)(shapeMass, fInclusiveJets[i], info);
919  fGenSubtractorInfoJetMass.push_back(info);
920  }
921 #endif
922  return 0;
923 }
924 
925 //_________________________________________________________________________________________________
927  //Do generic subtraction for jet mass
928 #ifdef FASTJET_VERSION
929  CreateGenSub();
930 
931  if(ijet>fInclusiveJets.size()) return 0;
932 
933  fGRNumerator.clear();
934  fGRDenominator.clear();
935  fGRNumeratorSub.clear();
936  fGRDenominatorSub.clear();
937 
938  // Define jet shape
939  for(Double_t r = 0.; r<fRMax; r+=fDRStep) {
940  AliJetShapeGRNum shapeGRNum(r,fDRStep);
941  AliJetShapeGRDen shapeGRDen(r,fDRStep);
942 
943  // clear the generic subtractor info vector
944  fGenSubtractorInfoGRNum.clear();
945  fGenSubtractorInfoGRDen.clear();
946  fj::contrib::GenericSubtractorInfo infoNum;
947  fj::contrib::GenericSubtractorInfo infoDen;
948  if(fInclusiveJets[ijet].perp()>1.e-4) {
949  double sub_num = (*fGenSubtractor)(shapeGRNum, fInclusiveJets[ijet], infoNum);
950  double sub_den = (*fGenSubtractor)(shapeGRDen, fInclusiveJets[ijet], infoDen);
951  }
952  fGenSubtractorInfoGRNum.push_back(infoNum);
953  fGenSubtractorInfoGRDen.push_back(infoDen);
954  fGRNumerator.push_back(infoNum.unsubtracted());
955  fGRDenominator.push_back(infoDen.unsubtracted());
956  fGRNumeratorSub.push_back(infoNum.second_order_subtracted());
957  fGRDenominatorSub.push_back(infoDen.second_order_subtracted());
958  }
959 #endif
960  return 0;
961 }
962 //_________________________________________________________________________________________________
964  //Do generic subtraction for jet mass
965 #ifdef FASTJET_VERSION
966  CreateGenSub();
967 
968  // Define jet shape
969  AliJetShapeAngularity shapeAngularity;
970 
971  // clear the generic subtractor info vector
972  fGenSubtractorInfoJetAngularity.clear();
973  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
974  fj::contrib::GenericSubtractorInfo infoAng;
975  if(fInclusiveJets[i].perp()>1.e-4)
976  double subtracted_shape = (*fGenSubtractor)(shapeAngularity, fInclusiveJets[i], infoAng);
977  fGenSubtractorInfoJetAngularity.push_back(infoAng);
978  }
979 #endif
980  return 0;
981 }
982 //_________________________________________________________________________________________________
984  //Do generic subtraction for jet mass
985 #ifdef FASTJET_VERSION
986  CreateGenSub();
987 
988  // Define jet shape
989  AliJetShapepTD shapepTD;
990 
991  // clear the generic subtractor info vector
992  fGenSubtractorInfoJetpTD.clear();
993  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
994  fj::contrib::GenericSubtractorInfo infopTD;
995  if(fInclusiveJets[i].perp()>1.e-4)
996  double subtracted_shape = (*fGenSubtractor)(shapepTD, fInclusiveJets[i], infopTD);
997  fGenSubtractorInfoJetpTD.push_back(infopTD);
998  }
999 #endif
1000  return 0;
1001 }
1002 //_________________________________________________________________________________________________
1004  //Do generic subtraction for jet mass
1005 #ifdef FASTJET_VERSION
1006  CreateGenSub();
1007 
1008  // Define jet shape
1009  AliJetShapeCircularity shapecircularity;
1010 
1011  // clear the generic subtractor info vector
1012  fGenSubtractorInfoJetCircularity.clear();
1013  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1014  fj::contrib::GenericSubtractorInfo infoCirc;
1015  if(fInclusiveJets[i].perp()>1.e-4)
1016  double subtracted_shape = (*fGenSubtractor)(shapecircularity, fInclusiveJets[i], infoCirc);
1017  fGenSubtractorInfoJetCircularity.push_back(infoCirc);
1018  }
1019 #endif
1020  return 0;
1021 }
1022 //_________________________________________________________________________________________________
1024  //Do generic subtraction for jet mass
1025 #ifdef FASTJET_VERSION
1026  CreateGenSub();
1027 
1028  // Define jet shape
1029  AliJetShapeSigma2 shapesigma2;
1030 
1031  // clear the generic subtractor info vector
1032  fGenSubtractorInfoJetSigma2.clear();
1033  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1034  fj::contrib::GenericSubtractorInfo infoSigma;
1035  if(fInclusiveJets[i].perp()>1.e-4)
1036  double subtracted_shape = (*fGenSubtractor)(shapesigma2, fInclusiveJets[i], infoSigma);
1037  fGenSubtractorInfoJetSigma2.push_back(infoSigma);
1038  }
1039 #endif
1040  return 0;
1041 }
1042 //_________________________________________________________________________________________________
1044  //Do generic subtraction for jet mass
1045 #ifdef FASTJET_VERSION
1046  CreateGenSub();
1047 
1048  // Define jet shape
1049  AliJetShapeConstituent shapeconst;
1050 
1051  // clear the generic subtractor info vector
1052  fGenSubtractorInfoJetConstituent.clear();
1053  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1054  fj::contrib::GenericSubtractorInfo infoConst;
1055  if(fInclusiveJets[i].perp()>1.e-4)
1056  double subtracted_shape = (*fGenSubtractor)(shapeconst, fInclusiveJets[i], infoConst);
1057  fGenSubtractorInfoJetConstituent.push_back(infoConst);
1058  }
1059 #endif
1060  return 0;
1061 }
1062 
1063 //_________________________________________________________________________________________________
1065  //Do generic subtraction for jet mass
1066 #ifdef FASTJET_VERSION
1067  CreateGenSub();
1068 
1069  // Define jet shape
1070  AliJetShapeLeSub shapeLeSub;
1071 
1072  // clear the generic subtractor info vector
1073  fGenSubtractorInfoJetLeSub.clear();
1074  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1075  fj::contrib::GenericSubtractorInfo infoLeSub;
1076  if(fInclusiveJets[i].perp()>1.e-4)
1077  double subtracted_shape = (*fGenSubtractor)(shapeLeSub, fInclusiveJets[i], infoLeSub);
1078  fGenSubtractorInfoJetLeSub.push_back(infoLeSub);
1079  }
1080 #endif
1081  return 0;
1082 }
1083 
1084 //_________________________________________________________________________________________________
1086  //Do generic subtraction for 1subjettiness
1087 #ifdef FASTJET_VERSION
1088  // Define jet shape
1089  AliJetShape1subjettiness_kt shape1subjettiness_kt;
1090  DoGenericSubtraction(shape1subjettiness_kt, fGenSubtractorInfoJet1subjettiness_kt);
1091 #endif
1092  return 0;
1093 }
1094 
1095 //_________________________________________________________________________________________________
1097  //Do generic subtraction for 2subjettiness
1098 #ifdef FASTJET_VERSION
1099  // Define jet shape
1100  AliJetShape2subjettiness_kt shape2subjettiness_kt;
1101  DoGenericSubtraction(shape2subjettiness_kt, fGenSubtractorInfoJet2subjettiness_kt);
1102 #endif
1103  return 0;
1104 }
1105 
1106 //_________________________________________________________________________________________________
1108  //Do generic subtraction for 3subjettiness
1109 #ifdef FASTJET_VERSION
1110  // Define jet shape
1111  AliJetShape3subjettiness_kt shape3subjettiness_kt;
1112  DoGenericSubtraction(shape3subjettiness_kt, fGenSubtractorInfoJet3subjettiness_kt);
1113 #endif
1114  return 0;
1115 }
1116 
1117 //_________________________________________________________________________________________________
1119  //Do generic subtraction for 2subjettiness axes opening angle
1120 #ifdef FASTJET_VERSION
1121  // Define jet shape
1122  AliJetShapeOpeningAngle_kt shapeOpeningAngle_kt;
1123  DoGenericSubtraction(shapeOpeningAngle_kt, fGenSubtractorInfoJetOpeningAngle_kt);
1124 #endif
1125  return 0;
1126 }
1127 
1128 //_________________________________________________________________________________________________
1130  //Do generic subtraction for 1subjettiness
1131 #ifdef FASTJET_VERSION
1132  // Define jet shape
1133  AliJetShape1subjettiness_ca shape1subjettiness_ca;
1134  DoGenericSubtraction(shape1subjettiness_ca, fGenSubtractorInfoJet1subjettiness_ca);
1135 #endif
1136  return 0;
1137 }
1138 
1139 //_________________________________________________________________________________________________
1141  //Do generic subtraction for 2subjettiness
1142 #ifdef FASTJET_VERSION
1143  // Define jet shape
1144  AliJetShape2subjettiness_ca shape2subjettiness_ca;
1145  DoGenericSubtraction(shape2subjettiness_ca, fGenSubtractorInfoJet2subjettiness_ca);
1146 #endif
1147  return 0;
1148 }
1149 
1150 //_________________________________________________________________________________________________
1152  //Do generic subtraction for 2subjettiness axes opening angle
1153 #ifdef FASTJET_VERSION
1154  // Define jet shape
1155  AliJetShapeOpeningAngle_ca shapeOpeningAngle_ca;
1156  DoGenericSubtraction(shapeOpeningAngle_ca, fGenSubtractorInfoJetOpeningAngle_ca);
1157 #endif
1158  return 0;
1159 }
1160 
1161 //_________________________________________________________________________________________________
1163  //Do generic subtraction for 1subjettiness
1164 #ifdef FASTJET_VERSION
1165  // Define jet shape
1166  AliJetShape1subjettiness_akt02 shape1subjettiness_akt02;
1167  DoGenericSubtraction(shape1subjettiness_akt02, fGenSubtractorInfoJet1subjettiness_akt02);
1168 #endif
1169  return 0;
1170 }
1171 
1172 //_________________________________________________________________________________________________
1174  //Do generic subtraction for 2subjettiness
1175 #ifdef FASTJET_VERSION
1176  // Define jet shape
1177  AliJetShape2subjettiness_akt02 shape2subjettiness_akt02;
1178  DoGenericSubtraction(shape2subjettiness_akt02, fGenSubtractorInfoJet2subjettiness_akt02);
1179 #endif
1180  return 0;
1181 }
1182 
1183 //_________________________________________________________________________________________________
1185  //Do generic subtraction for 2subjettiness axes opening angle
1186 #ifdef FASTJET_VERSION
1187  // Define jet shape
1188  AliJetShapeOpeningAngle_akt02 shapeOpeningAngle_akt02;
1189  DoGenericSubtraction(shapeOpeningAngle_akt02, fGenSubtractorInfoJetOpeningAngle_akt02);
1190 #endif
1191  return 0;
1192 }
1193 
1194 //_________________________________________________________________________________________________
1196  //Do generic subtraction for 1subjettiness
1197 #ifdef FASTJET_VERSION
1198  // Define jet shape
1199  AliJetShape1subjettiness_casd shape1subjettiness_casd;
1200  DoGenericSubtraction(shape1subjettiness_casd, fGenSubtractorInfoJet1subjettiness_casd);
1201 #endif
1202  return 0;
1203 }
1204 
1205 //_________________________________________________________________________________________________
1207  //Do generic subtraction for 2subjettiness
1208 #ifdef FASTJET_VERSION
1209  // Define jet shape
1210  AliJetShape2subjettiness_casd shape2subjettiness_casd;
1211  DoGenericSubtraction(shape2subjettiness_casd, fGenSubtractorInfoJet2subjettiness_casd);
1212 #endif
1213  return 0;
1214 }
1215 
1216 
1217 //_________________________________________________________________________________________________
1219  //Do generic subtraction for 2subjettiness axes opening angle
1220 #ifdef FASTJET_VERSION
1221  // Define jet shape
1222  AliJetShapeOpeningAngle_casd shapeOpeningAngle_casd;
1223  DoGenericSubtraction(shapeOpeningAngle_casd, fGenSubtractorInfoJetOpeningAngle_casd);
1224 #endif
1225  return 0;
1226 }
1227 
1228 //_________________________________________________________________________________________________
1230  //Do constituent subtraction
1231 #ifdef FASTJET_VERSION
1232  CreateConstituentSub();
1233  // fConstituentSubtractor->set_alpha(/* double alpha */);
1234  // fConstituentSubtractor->set_max_deltaR(/* double max_deltaR */);
1235 
1236  //clear constituent subtracted jets
1237  fConstituentSubtrJets.clear();
1238  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
1239  fj::PseudoJet subtracted_jet(0.,0.,0.,0.);
1240  if(fInclusiveJets[i].perp()>0.)
1241  subtracted_jet = (*fConstituentSubtractor)(fInclusiveJets[i]);
1242  fConstituentSubtrJets.push_back(subtracted_jet);
1243  }
1244  if(fConstituentSubtractor) { delete fConstituentSubtractor; fConstituentSubtractor = NULL; }
1245 
1246 #endif
1247  return 0;
1248 }
1249 
1250 //_________________________________________________________________________________________________
1252  //Do constituent subtraction
1253 #ifdef FASTJET_VERSION
1254  CreateEventConstituentSub();
1255  if(fUseMaxDelR) fEventConstituentSubtractor->set_max_standardDeltaR(fMaxDelR);
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); }
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);
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 
1472  if (Algorithm==0){
1473  Beta = 1.0;
1474  fR=0.4;
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  Result= nSub.result(fFilteredJets[0]);
1545  SubJet_Axes=nSub.currentAxes();
1546  SubJets=nSub.currentSubjets();
1547  }
1548 
1549 
1550 
1551  SubJet1_Axis=SubJet_Axes[0];
1552  Double_t SubJet1_Eta=SubJet1_Axis.pseudorapidity();
1553  Double_t SubJet2_Eta=0.0;
1554  Double_t SubJet1_Phi=SubJet1_Axis.phi();
1555  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
1556  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
1557  Double_t SubJet2_Phi=0.0;
1558  Double_t DeltaPhi=-5;
1559  if (SubJet_Axes.size()>1){
1560  SubJet2_Axis=SubJet_Axes[1];
1561  SubJet2_Eta=SubJet2_Axis.pseudorapidity();
1562  SubJet2_Phi=SubJet2_Axis.phi();
1563  if(SubJet2_Phi < -1*TMath::Pi()) SubJet2_Phi += (2*TMath::Pi());
1564  else if (SubJet2_Phi > TMath::Pi()) SubJet2_Phi -= (2*TMath::Pi());
1565  DeltaPhi=SubJet1_Phi-SubJet2_Phi;
1566  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1567  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1568  }
1569 
1570 
1571  SubJet1=SubJets[0];
1572  Double_t SubJet1Eta=SubJet1.pseudorapidity();
1573  Double_t SubJet2Eta=0.0;
1574  Double_t SubJet1Phi=SubJet1.phi();
1575  if(SubJet1Phi < -1*TMath::Pi()) SubJet1Phi += (2*TMath::Pi());
1576  else if (SubJet1Phi > TMath::Pi()) SubJet1Phi -= (2*TMath::Pi());
1577  Double_t SubJet2Phi=0.0;
1578  Double_t DeltaPhiSubJets=-5;
1579  Double_t SubJet1LeadingTrackPt=-3.0;
1580  Double_t SubJet2LeadingTrackPt=-3.0;
1581  std::vector<fj::PseudoJet> SubJet1Tracks = SubJet1.constituents();
1582  for (Int_t i=0; i<SubJet1Tracks.size(); i++){
1583  if (SubJet1Tracks[i].perp() > SubJet1LeadingTrackPt) SubJet1LeadingTrackPt=SubJet1Tracks[i].perp();
1584  }
1585  if (SubJet_Axes.size()>1){
1586  SubJet2=SubJets[1];
1587  SubJet2Eta=SubJet2.pseudorapidity();
1588  SubJet2Phi=SubJet2.phi();
1589  if(SubJet2Phi < -1*TMath::Pi()) SubJet2Phi += (2*TMath::Pi());
1590  else if (SubJet2Phi > TMath::Pi()) SubJet2Phi -= (2*TMath::Pi());
1591  DeltaPhiSubJets=SubJet1Phi-SubJet2Phi;
1592  if(DeltaPhiSubJets < -1*TMath::Pi()) DeltaPhiSubJets += (2*TMath::Pi());
1593  else if (DeltaPhiSubJets > TMath::Pi()) DeltaPhiSubJets -= (2*TMath::Pi());
1594  std::vector<fj::PseudoJet> SubJet2Tracks = SubJet2.constituents();
1595  for (Int_t i=0; i<SubJet2Tracks.size(); i++){
1596  if (SubJet2Tracks[i].perp() > SubJet2LeadingTrackPt) SubJet2LeadingTrackPt=SubJet2Tracks[i].perp();
1597  }
1598  }
1599 
1600 
1601 
1602 
1603  if (Option==0) return Result;
1604  else if (Option==1 && SubJet_Axes.size()>1 && N==2) return TMath::Sqrt(TMath::Power(SubJet1_Eta-SubJet2_Eta,2)+TMath::Power(DeltaPhi,2));
1605  else if (Option==2 && SubJet_Axes.size()>1 && N==2) return TMath::Sqrt(TMath::Power(SubJet1_Eta-SubJet2_Eta,2)+TMath::Power(DeltaPhi,2));
1606  else if ((Option==3 || Option==4) && N==2) return Result_SoftDrop;
1607  else if (Option==5 && SubJets.size()>1 && N==2) return SubJet1.perp();
1608  else if (Option==6 && SubJets.size()>1 && N==2) return SubJet2.perp();
1609  else if (Option==7 && SubJets.size()>1 && N==2) return TMath::Sqrt(TMath::Power(SubJet1Eta-SubJet2Eta,2)+TMath::Power(DeltaPhiSubJets,2));
1610  else if (Option==8 && SubJets.size()>1 && N==2) return SubJet1LeadingTrackPt;
1611  else if (Option==9 && SubJets.size()>1 && N==2) return SubJet2LeadingTrackPt;
1612  else return -2;
1613 }
1614 
1615 
1616 
1617 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
1618 
1619  Double_t Result=-1;
1620  std::vector<fastjet::PseudoJet> SubJet_Axes;
1621  fj::PseudoJet SubJet1_Axis;
1622  fj::PseudoJet SubJet2_Axis;
1623 
1624 
1625  if (SoftDropOn==1){
1626  // Beta_SD=0.0; //change these later so they are actually variable. currently not variable due to Kine trains
1627  // ZCut=0.1;
1628  fj::contrib::SoftDrop Soft_Drop(Beta_SD,ZCut);
1629  //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
1630  fj::PseudoJet Soft_Dropped_Jet=Soft_Drop(jet);
1631  if (Soft_Dropped_Jet==0) return -3;
1632  else{
1633  if (Soft_Dropped_Jet.constituents().size()>=N) jet=Soft_Dropped_Jet;
1634  else if (Soft_Dropped_Jet.constituents().size()<N) return -1;
1635  }
1636  }
1637 
1638 
1639  if (Algorithm==0){
1640  fj::contrib::Nsubjettiness nSub(N, fj::contrib::KT_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1641  Result= nSub.result(jet);
1642  SubJet_Axes=nSub.currentAxes();
1643  }
1644  else if (Algorithm==1) {
1645  fj::contrib::Nsubjettiness nSub(N, fj::contrib::CA_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1646  Result= nSub.result(jet);
1647  SubJet_Axes=nSub.currentAxes();
1648  }
1649  else if (Algorithm==2){
1650  fj::contrib::Nsubjettiness nSub(N, fj::contrib::AntiKT_Axes(Radius), fj::contrib::NormalizedMeasure(Beta,JetR));
1651  Result= nSub.result(jet);
1652  SubJet_Axes=nSub.currentAxes();
1653  }
1654  else if (Algorithm==3) {
1655  fj::contrib::Nsubjettiness nSub(N, fj::contrib::WTA_KT_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1656  Result= nSub.result(jet);
1657  SubJet_Axes=nSub.currentAxes();
1658  }
1659  else if (Algorithm==4) {
1660  fj::contrib::Nsubjettiness nSub(N, fj::contrib::WTA_CA_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1661  Result= nSub.result(jet);
1662  SubJet_Axes=nSub.currentAxes();
1663  }
1664  else if (Algorithm==5){
1665  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_KT_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1666  Result= nSub.result(jet);
1667  SubJet_Axes=nSub.currentAxes();
1668  }
1669  else if (Algorithm==6){
1670  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_CA_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1671  Result= nSub.result(jet);
1672  SubJet_Axes=nSub.currentAxes();
1673  }
1674  else if (Algorithm==7){
1675  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_AntiKT_Axes(Radius), fj::contrib::NormalizedMeasure(Beta,JetR));
1676  Result= nSub.result(jet);
1677  SubJet_Axes=nSub.currentAxes();
1678  }
1679  else if (Algorithm==8){
1680  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_WTA_KT_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1681  Result= nSub.result(jet);
1682  SubJet_Axes=nSub.currentAxes();
1683  }
1684  else if (Algorithm==9){
1685  fj::contrib::Nsubjettiness nSub(N, fj::contrib::OnePass_WTA_CA_Axes(), fj::contrib::NormalizedMeasure(Beta,JetR));
1686  Result= nSub.result(jet);
1687  SubJet_Axes=nSub.currentAxes();
1688  }
1689  else if (Algorithm==10){
1690  fj::contrib::Nsubjettiness nSub(N, fj::contrib::MultiPass_Axes(100), fj::contrib::NormalizedMeasure(Beta,JetR));
1691  Result= nSub.result(jet);
1692  SubJet_Axes=nSub.currentAxes();
1693  }
1694 
1695  SubJet1_Axis=SubJet_Axes[0];
1696  Double_t SubJet1_Eta=SubJet1_Axis.pseudorapidity();
1697  Double_t SubJet2_Eta;
1698  Double_t SubJet1_Phi=SubJet1_Axis.phi();
1699  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
1700  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
1701  Double_t SubJet2_Phi;
1702  Double_t DeltaPhi=-5;
1703  if (SubJet_Axes.size()>1){
1704  SubJet2_Axis=SubJet_Axes[1];
1705  SubJet2_Eta=SubJet2_Axis.pseudorapidity();
1706  SubJet2_Phi=SubJet2_Axis.phi();
1707  if(SubJet2_Phi < -1*TMath::Pi()) SubJet2_Phi += (2*TMath::Pi());
1708  else if (SubJet2_Phi > TMath::Pi()) SubJet2_Phi -= (2*TMath::Pi());
1709  DeltaPhi=SubJet1_Phi-SubJet2_Phi;
1710  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1711  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1712  }
1713 
1714 
1715  if (Option==0) return Result;
1716  else if (Option==1 && SubJet_Axes.size()>1 && N==2) return TMath::Sqrt(TMath::Power(SubJet1_Eta-SubJet2_Eta,2)+TMath::Power(DeltaPhi,2));
1717  else if (Option==2 && SubJet_Axes.size()>1 && N==2) return TMath::Sqrt(TMath::Power(SubJet1_Eta-SubJet2_Eta,2)+TMath::Power(DeltaPhi,2));
1718  else return -2;
1719 
1720 }
1721 
1722 
1723 
1724 
1725 
1726 
1727 
1728 
1729 
1730 #endif
std::vector< fastjet::PseudoJet > fEventSubInputVectors
Definition: AliFJWrapper.h:150
fastjet::PseudoJet GetJetAreaVector(UInt_t idx) const
fastjet::PseudoJet GetEventSubJetAreaVector(UInt_t idx) const
std::vector< double > fGRNumeratorSub
Definition: AliFJWrapper.h:235
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
Double_t GetFilteredJetArea(UInt_t idx) const
std::vector< fastjet::PseudoJet > fInputVectors
Definition: AliFJWrapper.h:149
Double_t fR
Definition: AliFJWrapper.h:180
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:171
virtual Int_t DoEventConstituentSubtraction()
virtual Int_t DoSoftDrop()
Int_t fPluginAlgor
Definition: AliFJWrapper.h:186
fastjet::ClusterSequenceActiveAreaExplicitGhosts * fClustSeqActGhosts
Definition: AliFJWrapper.h:172
virtual Int_t DoGenericSubtractionJetOpeningAngle_kt()
Double_t fKtScatter
Definition: AliFJWrapper.h:184
Int_t fNGhostRepeats
Definition: AliFJWrapper.h:177
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:153
fastjet::JetAlgorithm fAlgor
Definition: AliFJWrapper.h:174
const char * GetTitle() const
Definition: AliFJWrapper.h:42
fastjet::RecombinationScheme fScheme
Definition: AliFJWrapper.h:175
std::vector< fastjet::PseudoJet > GetEventSubJetConstituents(UInt_t idx) const
Double_t fMaxDelR
Definition: AliFJWrapper.h:195
Double_t fRMax
Definition: AliFJWrapper.h:231
virtual void CopySettingsFrom(const AliFJWrapper &wrapper)
Bool_t fUseMaxDelR
Definition: AliFJWrapper.h:194
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:181
fastjet::PseudoJet GetFilteredJetAreaVector(UInt_t idx) const
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
fastjet::RangeDefinition * fRange
Definition: AliFJWrapper.h:165
Double_t fZcut
Definition: AliFJWrapper.h:191
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
virtual Int_t DoGenericSubtractionJet2subjettiness_casd()
std::vector< fastjet::PseudoJet > fFilteredJets
Definition: AliFJWrapper.h:155
void SetupAreaTypefromOpt(const char *option)
virtual void ClearMemory()
fastjet::AreaType fAreaType
Definition: AliFJWrapper.h:176
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
Double_t fMaxRap
Definition: AliFJWrapper.h:179
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
virtual Int_t DoGenericSubtractionJetOpeningAngle_casd()
int Int_t
Definition: External.C:63
fastjet::ClusterSequenceArea * fClustSeqES
Definition: AliFJWrapper.h:170
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:156
fastjet::JetDefinition::Plugin * fPlugin
Definition: AliFJWrapper.h:163
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:226
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
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:158
void SetUseArea4Vector(Bool_t useA4v)
Definition: AliFJWrapper.h:132
fastjet::ClusterSequenceArea * fClustSeq
Definition: AliFJWrapper.h:169
void SetName(const char *name)
Definition: AliFJWrapper.h:118
Double_t fBeta
Definition: AliFJWrapper.h:192
Definition: Option.C:68
virtual void SubtractBackground(const Double_t median_pt=-1)
std::vector< fastjet::PseudoJet > fInputGhosts
Definition: AliFJWrapper.h:152
Double_t fGridScatter
Definition: AliFJWrapper.h:183
virtual void Clear(const Option_t *="")
Bool_t fLegacyMode
Definition: AliFJWrapper.h:227
fastjet::VoronoiAreaSpec * fVorAreaSpec
Definition: AliFJWrapper.h:160
Int_t mode
Definition: anaM.C:41
Double_t fDRStep
Definition: AliFJWrapper.h:232
std::vector< double > fGRDenominatorSub
Definition: AliFJWrapper.h:236
Bool_t GetLegacyMode()
Definition: AliFJWrapper.h:51
fastjet::GhostedAreaSpec * fGhostedAreaSpec
Definition: AliFJWrapper.h:161
Bool_t fUseArea4Vector
Definition: AliFJWrapper.h:189
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:234
virtual Int_t DoGenericSubtractionJet1subjettiness_kt()
virtual Int_t DoGenericSubtractionJet1subjettiness_casd()
Double_t fGhostArea
Definition: AliFJWrapper.h:178
virtual void GetMedianAndSigma(Double_t &median, Double_t &sigma, Int_t remove=0) const
TString fName
Definition: AliFJWrapper.h:147
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:188
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
fastjet::JetDefinition * fJetDef
Definition: AliFJWrapper.h:162
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:159
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:228
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:230
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:157
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:154
std::vector< fastjet::PseudoJet > fEventSubCorrectedVectors
Definition: AliFJWrapper.h:151
void SetMeanGhostKt(Double_t meankt)
Definition: AliFJWrapper.h:130
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:193
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:229
Bool_t GetDoFilterArea()
Definition: AliFJWrapper.h:52
virtual Int_t DoGenericSubtractionJetLeSub()
virtual Int_t DoGenericSubtractionJetConstituent()
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:233
Double_t fMeanGhostKt
Definition: AliFJWrapper.h:185
fastjet::Strategy fStrategy
Definition: AliFJWrapper.h:173
TString fTitle
Definition: AliFJWrapper.h:148