[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliFJWrapper.h

#ifndef AliFJWrapper_H
#define AliFJWrapper_H

// $Id$

#if !defined(__CINT__) 

#include <vector>
#include <TString.h>
#include "AliLog.h"
#include "FJ_includes.h"

class AliFJWrapper
{
 public:
  AliFJWrapper(const char *name, const char *title);
  virtual ~AliFJWrapper();

  virtual void  AddInputVector (Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index = -99999);
  virtual void  AddInputVector (const fastjet::PseudoJet& vec,                Int_t index = -99999);
  virtual void  AddInputVectors(const std::vector<fastjet::PseudoJet>& vecs,  Int_t offsetIndex = -99999);
  virtual const char *ClassName()                            const { return "AliFJWrapper";              }
  virtual void  Clear(const Option_t* /*opt*/ = "");
  virtual void  CopySettingsFrom (const AliFJWrapper& wrapper);
  virtual void  GetMedianAndSigma(Double_t& median, Double_t& sigma, Int_t remove = 0) const;
  fastjet::ClusterSequenceArea*         GetClusterSequence() const { return fClustSeq;                   }
  const std::vector<fastjet::PseudoJet> GetInputVectors()    const { return fInputVectors;               }
  const std::vector<fastjet::PseudoJet> GetInclusiveJets()   const { return fInclusiveJets;              }
  std::vector<fastjet::PseudoJet>       GetJetConstituents(UInt_t idx) const;
  Double_t                              GetMedianUsedForBgSubtraction() const { return fMedUsedForBgSub; }
  const char*                           GetName()            const { return fName;                       }
  const char*                           GetTitle()           const { return fTitle;                      }
  Double_t                              GetJetArea         (UInt_t idx) const;
  fastjet::PseudoJet                    GetJetAreaVector   (UInt_t idx) const;
  Double_t                              GetJetSubtractedPt (UInt_t idx) const;
  virtual std::vector<double>           GetSubtractedJetsPts(Double_t median_pt = -1, Bool_t sorted = kFALSE);
  Bool_t                                GetLegacyMode()            { return fLegacyMode; }

  virtual Int_t Run();

  void SetStrategy(const fastjet::Strategy &strat)                 { fStrategy = strat;  }
  void SetAlgorithm(const fastjet::JetAlgorithm &algor)            { fAlgor    = algor;  }
  void SetRecombScheme(const fastjet::RecombinationScheme &scheme) { fScheme   = scheme; }
  void SetAreaType(const fastjet::AreaType &atype)                 { fAreaType = atype;  }
  void SetNRepeats(Int_t nrepeat)       { fNGhostRepeats  = nrepeat; }
  void SetGhostArea(Double_t gharea)    { fGhostArea      = gharea;  }
  void SetMaxRap(Double_t maxrap)       { fMaxRap         = maxrap;  }
  void SetR(Double_t r)                 { fR              = r;       }
  void SetGridScatter(Double_t gridSc)  { fGridScatter    = gridSc;  }
  void SetKtScatter(Double_t ktSc)      { fKtScatter      = ktSc;    }
  void SetMeanGhostKt(Double_t meankt)  { fMeanGhostKt    = meankt;  }
  void SetPluginAlgor(Int_t plugin)     { fPluginAlgor    = plugin;  }
  void SetUseArea4Vector(Bool_t useA4v) { fUseArea4Vector = useA4v;  }
  void SetupAlgorithmfromOpt(const char *option);
  void SetupAreaTypefromOpt(const char *option);
  void SetupSchemefromOpt(const char *option);
  void SetupStrategyfromOpt(const char *option);
  void SetLegacyMode (Bool_t mode)      { fLegacyMode ^= mode; }
  void SetLegacyFJ();

 protected:
  TString                                fName;             //!
  TString                                fTitle;            //!
  std::vector<fastjet::PseudoJet>        fInputVectors;     //!
  std::vector<fastjet::PseudoJet>        fInclusiveJets;    //!
  std::vector<double>                    fSubtractedJetsPt; //!
  fastjet::AreaDefinition               *fAreaDef;          //!
  fastjet::VoronoiAreaSpec              *fVorAreaSpec;      //!
  fastjet::GhostedAreaSpec              *fGhostedAreaSpec;  //!
  fastjet::JetDefinition                *fJetDef;           //!
  fastjet::JetDefinition::Plugin        *fPlugin;           //!
  fastjet::RangeDefinition              *fRange;            //!
  fastjet::ClusterSequenceArea          *fClustSeq;         //!
  fastjet::Strategy                      fStrategy;         //!
  fastjet::JetAlgorithm                  fAlgor;            //!
  fastjet::RecombinationScheme           fScheme;           //!
  fastjet::AreaType                      fAreaType;         //!
  Int_t                                  fNGhostRepeats;    //!
  Double_t                               fGhostArea;	    //!
  Double_t                               fMaxRap;	    //!
  Double_t                               fR;                //!
  // no setters for the moment - used default values in the constructor
  Double_t                               fGridScatter;      //!
  Double_t                               fKtScatter;	    //!
  Double_t                               fMeanGhostKt;      //!
  Int_t                                  fPluginAlgor;      //!
  // extra parameters
  Double_t                               fMedUsedForBgSub;  //!
  Bool_t                                 fUseArea4Vector;   //!
#ifdef FASTJET_VERSION
  fastjet::JetMedianBackgroundEstimator *fBkrdEstimator;    //!
#endif
  Bool_t                                 fLegacyMode;       //!

  virtual void   SubtractBackground(const Double_t median_pt = -1);

 private:
  AliFJWrapper();
  AliFJWrapper(const AliFJWrapper& wrapper);
  AliFJWrapper& operator = (const AliFJWrapper& wrapper);
};
#endif
#endif

#ifdef AliFJWrapper_CXX
#undef AliFJWrapper_CXX 

#if defined __GNUC__
#pragma GCC system_header
#endif

namespace fj = fastjet;

//_________________________________________________________________________________________________
AliFJWrapper::AliFJWrapper(const char *name, const char *title)
  : 
    fName              (name)
  , fTitle             (title)
  , fInputVectors      ( )
  , fInclusiveJets     ( )
  , fSubtractedJetsPt  ( )
  , fAreaDef           (0)
  , fVorAreaSpec       (0)
  , fGhostedAreaSpec   (0)
  , fJetDef            (0)
  , fPlugin            (0)
  , fRange             (0)
  , fClustSeq          (0)
  , fStrategy          (fj::Best)
  , fAlgor             (fj::kt_algorithm)
  , fScheme            (fj::BIpt_scheme)
  , fAreaType          (fj::active_area)
  , fNGhostRepeats     (1)
  , fGhostArea         (0.01)
  , fMaxRap            (1.)
  , fR                 (0.4)
  , fGridScatter       (1.0)
  , fKtScatter         (0.1)
  , fMeanGhostKt       (1e-100)
  , fPluginAlgor       (0)
  , fMedUsedForBgSub   (0)
  , fUseArea4Vector    (kFALSE)
#ifdef FASTJET_VERSION
  , fBkrdEstimator     (0)
#endif
  , fLegacyMode        (false)
{
  // Constructor.
}

//_________________________________________________________________________________________________
AliFJWrapper::~AliFJWrapper()
{  
  // Destructor.

  delete fAreaDef;
  delete fVorAreaSpec;
  delete fGhostedAreaSpec;
  delete fJetDef;
  delete fPlugin;
  delete fRange;
  delete fClustSeq;
#ifdef FASTJET_VERSION
  if (fBkrdEstimator) delete fBkrdEstimator;
#endif
}

//_________________________________________________________________________________________________
void AliFJWrapper::CopySettingsFrom(const AliFJWrapper& wrapper)
{
  // Copy some settings.
  // You very often want to keep most of the settings 
  // but change only the algorithm or R - do it after call to this function

  fStrategy       = wrapper.fStrategy;
  fAlgor          = wrapper.fAlgor;
  fScheme         = wrapper.fScheme;
  fAreaType       = wrapper.fAreaType;
  fNGhostRepeats  = wrapper.fNGhostRepeats;
  fGhostArea      = wrapper.fGhostArea;
  fMaxRap         = wrapper.fMaxRap;
  fR              = wrapper.fR;
  fGridScatter    = wrapper.fGridScatter;
  fKtScatter      = wrapper.fKtScatter;
  fMeanGhostKt    = wrapper.fMeanGhostKt;
  fPluginAlgor    = wrapper.fPluginAlgor;
  fUseArea4Vector = wrapper.fUseArea4Vector;
  fLegacyMode     = wrapper.fLegacyMode;
}

//_________________________________________________________________________________________________
void AliFJWrapper::Clear(const Option_t */*opt*/)
{
  // Simply clear the input vectors.
  // Make sure done on every event if the instance is reused
  // Reset the median to zero.

  fInputVectors.clear();
  fMedUsedForBgSub = 0;

  // for the moment brute force delete everything
  delete fAreaDef;         fAreaDef         = 0;
  delete fVorAreaSpec;     fVorAreaSpec     = 0;
  delete fGhostedAreaSpec; fGhostedAreaSpec = 0; 
  delete fJetDef;          fJetDef          = 0;
  delete fPlugin;          fPlugin          = 0;
  delete fRange;           fRange           = 0;
  delete fClustSeq;        fClustSeq        = 0;
#ifdef FASTJET_VERSION
  if (fBkrdEstimator) delete fBkrdEstimator; fBkrdEstimator = 0;
#endif
}

//_________________________________________________________________________________________________
void AliFJWrapper::AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index)
{
  // Make the input pseudojet.

  fastjet::PseudoJet inVec(px, py, pz, E);
  
  if (index > -99999) {
    inVec.set_user_index(index);
  } else {
    inVec.set_user_index(fInputVectors.size());
  }

  // add to the fj container of input vectors
  fInputVectors.push_back(inVec);
}

//_________________________________________________________________________________________________
void AliFJWrapper::AddInputVector(const fj::PseudoJet& vec, Int_t index)
{
  // Add an input pseudojet.

  fj::PseudoJet inVec = vec;
  
  if (index > -99999) {
    inVec.set_user_index(index);
  } else {
    inVec.set_user_index(fInputVectors.size());
  }

  // add to the fj container of input vectors
  fInputVectors.push_back(inVec);
}

//_________________________________________________________________________________________________
void AliFJWrapper::AddInputVectors(const std::vector<fj::PseudoJet>& vecs, Int_t offsetIndex)
{
  // Add the input from vector of pseudojets.

  for (UInt_t i = 0; i < vecs.size(); ++i) {
    fj::PseudoJet inVec = vecs[i];
    if (offsetIndex > -99999)
      inVec.set_user_index(fInputVectors.size() + offsetIndex);
    // add to the fj container of input vectors
    fInputVectors.push_back(inVec);
  }
}

//_________________________________________________________________________________________________
Double_t AliFJWrapper::GetJetArea(UInt_t idx) const
{
  // Get the jet area.

  Double_t retval = -1; // really wrong area..
  if ( idx < fInclusiveJets.size() ) {
    retval = fClustSeq->area(fInclusiveJets[idx]);
  } else {
    AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
  }
  return retval;
}

//_________________________________________________________________________________________________
fastjet::PseudoJet AliFJWrapper::GetJetAreaVector(UInt_t idx) const
{
  // Get the jet area as vector.
  fastjet::PseudoJet retval;
  if ( idx < fInclusiveJets.size() ) {
    retval = fClustSeq->area_4vector(fInclusiveJets[idx]);
  } else {
    AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
  }
  return retval;
}

//_________________________________________________________________________________________________
std::vector<double> AliFJWrapper::GetSubtractedJetsPts(Double_t median_pt, Bool_t sorted)
{ 
  // Get subtracted jets pTs, returns vector.

  SubtractBackground(median_pt);
  
  if (kTRUE == sorted) {
    std::sort(fSubtractedJetsPt.begin(), fSubtractedJetsPt.begin());
  }
  return fSubtractedJetsPt;
}

//_________________________________________________________________________________________________
Double_t AliFJWrapper::GetJetSubtractedPt(UInt_t idx) const
{
  // Get subtracted jets pTs, returns Double_t.

  Double_t retval = -99999.; // really wrong pt..
  if ( idx < fSubtractedJetsPt.size() ) {
    retval = fSubtractedJetsPt[idx];
  }
  return retval;
}

//_________________________________________________________________________________________________
std::vector<fastjet::PseudoJet> 
AliFJWrapper::GetJetConstituents(UInt_t idx) const
{
  // Get jets constituents.

  std::vector<fastjet::PseudoJet> retval;
  
  if ( idx < fInclusiveJets.size() ) {
    retval = fClustSeq->constituents(fInclusiveJets[idx]);
  } else {
    AliError(Form("[e] ::GetJetConstituents wrong index: %d",idx));
  }
  
  return retval;
}

//_________________________________________________________________________________________________
void AliFJWrapper::GetMedianAndSigma(Double_t &median, Double_t &sigma, Int_t remove) const
{
  // Get the median and sigma from fastjet.
  // User can also do it on his own because the cluster sequence is exposed (via a getter)

  if (!fClustSeq) {
    AliError("[e] Run the jfinder first.");
    return;
  }

  Double_t mean_area = 0;
  try {
    if(0 == remove) {
      fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
    }  else {
      std::vector<fastjet::PseudoJet> input_jets = sorted_by_pt(fClustSeq->inclusive_jets());
      input_jets.erase(input_jets.begin(), input_jets.begin() + remove);
      fClustSeq->get_median_rho_and_sigma(input_jets, *fRange, fUseArea4Vector, median, sigma, mean_area);
      input_jets.clear();
    }
  } catch (fj::Error) {
    AliError(" [w] FJ Exception caught.");
    median = -1.;
    sigma = -1;
  }
}

//_________________________________________________________________________________________________
Int_t AliFJWrapper::Run()
{
  // Run the actual jet finder.

  if (fAreaType == fj::voronoi_area) {
    // Rfact - check dependence - default is 1.
    // NOTE: hardcoded variable!
    fVorAreaSpec = new fj::VoronoiAreaSpec(1.); 
    fAreaDef     = new fj::AreaDefinition(*fVorAreaSpec);      
  } else {
    fGhostedAreaSpec = new fj::GhostedAreaSpec(fMaxRap,
                                               fNGhostRepeats, 
                                               fGhostArea,
                                               fGridScatter,
                                               fKtScatter,
                                               fMeanGhostKt);

    fAreaDef = new fj::AreaDefinition(*fGhostedAreaSpec, fAreaType);
  }
  
  // this is acceptable by fastjet:
  fRange = new fj::RangeDefinition(fMaxRap - 0.95 * fR);

  if (fAlgor == fj::plugin_algorithm) {
    if (fPluginAlgor == 0) {
      // SIS CONE ALGOR
      // NOTE: hardcoded split parameter
      Double_t overlap_threshold = 0.75; // NOTE: this actually splits a lot: thr/min(pt1,pt2)
      fPlugin = new fj::SISConePlugin(fR, 
                                      overlap_threshold,
                                      0,    //search of stable cones - zero = until no more
                                      1.0); // this should be seed effectively for proto jets
      fJetDef = new fastjet::JetDefinition(fPlugin);
    } else {
      AliError("[e] Unrecognized plugin number!");
    }
  } else {
    fJetDef = new fj::JetDefinition(fAlgor, fR, fScheme, fStrategy);
  }
  
  try {
    fClustSeq = new fj::ClusterSequenceArea(fInputVectors, *fJetDef, *fAreaDef);
  } catch (fj::Error) {
    AliError(" [w] FJ Exception caught.");
    return -1;
  }

  // FJ3 :: Define an JetMedianBackgroundEstimator just in case it will be used 
#ifdef FASTJET_VERSION
  //fBkrdEstimator = new fj::JetMedianBackgroundEstimator(*fRange, *fJetDef, *fAreaDef) ;
#endif

  if (fLegacyMode) { SetLegacyFJ(); } // for FJ 2.x even if fLegacyMode is set, SetLegacyFJ is dummy

  // inclusive jets:
  fInclusiveJets.clear();
  fInclusiveJets = fClustSeq->inclusive_jets(0.0); 

  return 0;
}

//_________________________________________________________________________________________________
void AliFJWrapper::SetLegacyFJ()
{
  // This methods enable legacy behaviour (FastJet 2.x) when AliROOT is compiled with FastJet 3.x
#ifdef FASTJET_VERSION
    std::cout << "WARNING! Setting FastJet in legacy mode" << std::endl;
    if (fGhostedAreaSpec) { fGhostedAreaSpec->set_fj2_placement(kTRUE); }
    if (fBkrdEstimator) {
      fBkrdEstimator->set_provide_fj2_sigma(kTRUE);
      fBkrdEstimator->set_use_area_4vector(kFALSE);
    } 
#endif
}

//_________________________________________________________________________________________________
void AliFJWrapper::SubtractBackground(Double_t median_pt)
{
  // Subtract the background (specify the value - see below the meaning).
  // Negative argument means the bg will be determined with the current algorithm
  // this is the default behavior. Zero is allowed
  // Note: user may set the switch for area4vector based subtraction.

  Double_t median    = 0;
  Double_t sigma     = 0;
  Double_t mean_area = 0;

  // clear the subtracted jet pt's vector<double>
  fSubtractedJetsPt.clear();

  // check what was specified (default is -1)
  if (median_pt < 0) {
    try {
      fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
    }

    catch (fj::Error) {
      AliError(" [w] FJ Exception caught.");
      median = -9999.;
      sigma = -1;
      fMedUsedForBgSub = median;
      return;
    }
    fMedUsedForBgSub = median;
  } else {
    // we do not know the sigma in this case
    sigma = -1;
    if (0.0 == median_pt) {
      // AliWarning(Form(" [w] Median specified for bg subtraction is ZERO: %f", median_pt )) << endl;
      fMedUsedForBgSub = 0.;
    } else {
      fMedUsedForBgSub = median_pt;
    }
  }

  // subtract:
  for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
    if ( fUseArea4Vector ) {
      // subtract the background using the area4vector
      fj::PseudoJet area4v = fClustSeq->area_4vector(fInclusiveJets[i]);
      fj::PseudoJet jet_sub = fInclusiveJets[i] - area4v * fMedUsedForBgSub;
      fSubtractedJetsPt.push_back(jet_sub.perp()); // here we put only the pt of the jet - note: this can be negative
    } else {
      // subtract the background using scalars
      // fj::PseudoJet jet_sub = fInclusiveJets[i] - area * fMedUsedForBgSub_;
      Double_t area = fClustSeq->area(fInclusiveJets[i]);
      // standard subtraction
      Double_t pt_sub = fInclusiveJets[i].perp() - fMedUsedForBgSub * area;
      fSubtractedJetsPt.push_back(pt_sub); // here we put only the pt of the jet - note: this can be negative
    }
  }
}

//_________________________________________________________________________________________________
void AliFJWrapper::SetupAlgorithmfromOpt(const char *option)
{
  // Setup algorithm from char.

  std::string opt(option);
  
  if (!opt.compare("kt"))                fAlgor    = fj::kt_algorithm;
  if (!opt.compare("antikt"))            fAlgor    = fj::antikt_algorithm;
  if (!opt.compare("cambridge"))         fAlgor    = fj::cambridge_algorithm;
  if (!opt.compare("genkt"))             fAlgor    = fj::genkt_algorithm;
  if (!opt.compare("cambridge_passive")) fAlgor    = fj::cambridge_for_passive_algorithm;
  if (!opt.compare("genkt_passive"))     fAlgor    = fj::genkt_for_passive_algorithm;
  if (!opt.compare("ee_kt"))             fAlgor    = fj::ee_kt_algorithm;
  if (!opt.compare("ee_genkt"))          fAlgor    = fj::ee_genkt_algorithm;
  if (!opt.compare("plugin"))            fAlgor    = fj::plugin_algorithm;
}

//_________________________________________________________________________________________________
void AliFJWrapper::SetupAreaTypefromOpt(const char *option)
{
  // Setup area type from char.

  std::string opt(option);

  if (!opt.compare("active"))                      fAreaType = fj::active_area;
  if (!opt.compare("invalid"))                     fAreaType = fj::invalid_area;
  if (!opt.compare("active_area_explicit_ghosts")) fAreaType = fj::active_area_explicit_ghosts;
  if (!opt.compare("one_ghost_passive"))           fAreaType = fj::one_ghost_passive_area;
  if (!opt.compare("passive"))                     fAreaType = fj::passive_area;
  if (!opt.compare("voronoi"))                     fAreaType = fj::voronoi_area;
}

//_________________________________________________________________________________________________
void AliFJWrapper::SetupSchemefromOpt(const char *option)
{
  //
  // setup scheme from char
  //

  std::string opt(option);

  if (!opt.compare("BIpt"))   fScheme   = fj::BIpt_scheme;
  if (!opt.compare("BIpt2"))  fScheme   = fj::BIpt2_scheme;
  if (!opt.compare("E"))      fScheme   = fj::E_scheme;
  if (!opt.compare("pt"))     fScheme   = fj::pt_scheme;
  if (!opt.compare("pt2"))    fScheme   = fj::pt2_scheme;
  if (!opt.compare("Et"))     fScheme   = fj::Et_scheme;
  if (!opt.compare("Et2"))    fScheme   = fj::Et2_scheme;
}

//_________________________________________________________________________________________________
void AliFJWrapper::SetupStrategyfromOpt(const char *option)
{
  // Setup strategy from char.

  std::string opt(option);
  
  if (!opt.compare("Best"))            fStrategy = fj::Best;
  if (!opt.compare("N2MinHeapTiled"))  fStrategy = fj::N2MinHeapTiled;
  if (!opt.compare("N2Tiled"))         fStrategy = fj::N2Tiled;
  if (!opt.compare("N2PoorTiled"))     fStrategy = fj::N2PoorTiled;
  if (!opt.compare("N2Plain"))         fStrategy = fj::N2Plain;
  if (!opt.compare("N3Dumb"))          fStrategy = fj::N3Dumb;
  if (!opt.compare("NlnN"))            fStrategy = fj::NlnN;
  if (!opt.compare("NlnN3pi"))         fStrategy = fj::NlnN3pi;
  if (!opt.compare("NlnN4pi"))         fStrategy = fj::NlnN4pi;
  if (!opt.compare("NlnNCam4pi"))      fStrategy = fj::NlnNCam4pi;
  if (!opt.compare("NlnNCam2pi2R"))    fStrategy = fj::NlnNCam2pi2R;
  if (!opt.compare("NlnNCam"))         fStrategy = fj::NlnNCam;
  if (!opt.compare("plugin"))          fStrategy = fj::plugin_strategy;
}
#endif
Commit	Line	Data
48d874ff	1	#ifndef AliFJWrapper_H
	2	#define AliFJWrapper_H
	3
	4	// $Id$
	5
a9c8e539	6	#if !defined(__CINT__)
48d874ff	7
48d874ff	8	#include <vector>
a9c8e539	9	#include <TString.h>
48d874ff	10	#include "AliLog.h"
	11	#include "FJ_includes.h"
	12
a9c8e539	13	class AliFJWrapper
48d874ff	14	{
	15	public:
	16	AliFJWrapper(const char name, const char title);
	17	virtual ~AliFJWrapper();
a9c8e539	18
48d874ff	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);
a9c8e539	22	virtual const char *ClassName() const { return "AliFJWrapper"; }
a9c8e539	23	virtual void Clear(const Option_t* /opt/ = "");
48d874ff	24	virtual void CopySettingsFrom (const AliFJWrapper& wrapper);
	25	virtual void GetMedianAndSigma(Double_t& median, Double_t& sigma, Int_t remove = 0) const;
	26	fastjet::ClusterSequenceArea* GetClusterSequence() const { return fClustSeq; }
	27	const std::vector<fastjet::PseudoJet> GetInputVectors() const { return fInputVectors; }
	28	const std::vector<fastjet::PseudoJet> GetInclusiveJets() const { return fInclusiveJets; }
	29	std::vector<fastjet::PseudoJet> GetJetConstituents(UInt_t idx) const;
	30	Double_t GetMedianUsedForBgSubtraction() const { return fMedUsedForBgSub; }
a9c8e539	31	const char* GetName() const { return fName; }
a9c8e539	32	const char* GetTitle() const { return fTitle; }
48d874ff	33	Double_t GetJetArea (UInt_t idx) const;
4e22c11c	34	fastjet::PseudoJet GetJetAreaVector (UInt_t idx) const;
48d874ff	35	Double_t GetJetSubtractedPt (UInt_t idx) const;
48d874ff	36	virtual std::vector<double> GetSubtractedJetsPts(Double_t median_pt = -1, Bool_t sorted = kFALSE);
5d0a5007	37	Bool_t GetLegacyMode() { return fLegacyMode; }
48d874ff	38
	39	virtual Int_t Run();
	40
	41	void SetStrategy(const fastjet::Strategy &strat) { fStrategy = strat; }
	42	void SetAlgorithm(const fastjet::JetAlgorithm &algor) { fAlgor = algor; }
	43	void SetRecombScheme(const fastjet::RecombinationScheme &scheme) { fScheme = scheme; }
	44	void SetAreaType(const fastjet::AreaType &atype) { fAreaType = atype; }
	45	void SetNRepeats(Int_t nrepeat) { fNGhostRepeats = nrepeat; }
	46	void SetGhostArea(Double_t gharea) { fGhostArea = gharea; }
	47	void SetMaxRap(Double_t maxrap) { fMaxRap = maxrap; }
	48	void SetR(Double_t r) { fR = r; }
	49	void SetGridScatter(Double_t gridSc) { fGridScatter = gridSc; }
	50	void SetKtScatter(Double_t ktSc) { fKtScatter = ktSc; }
	51	void SetMeanGhostKt(Double_t meankt) { fMeanGhostKt = meankt; }
	52	void SetPluginAlgor(Int_t plugin) { fPluginAlgor = plugin; }
	53	void SetUseArea4Vector(Bool_t useA4v) { fUseArea4Vector = useA4v; }
	54	void SetupAlgorithmfromOpt(const char *option);
	55	void SetupAreaTypefromOpt(const char *option);
	56	void SetupSchemefromOpt(const char *option);
	57	void SetupStrategyfromOpt(const char *option);
5d0a5007	58	void SetLegacyMode (Bool_t mode) { fLegacyMode ^= mode; }
5d0a5007	59	void SetLegacyFJ();
48d874ff	60
48d874ff	61	protected:
5d0a5007	62	TString fName; //!
	63	TString fTitle; //!
	64	std::vector<fastjet::PseudoJet> fInputVectors; //!
	65	std::vector<fastjet::PseudoJet> fInclusiveJets; //!
	66	std::vector<double> fSubtractedJetsPt; //!
	67	fastjet::AreaDefinition *fAreaDef; //!
	68	fastjet::VoronoiAreaSpec *fVorAreaSpec; //!
	69	fastjet::GhostedAreaSpec *fGhostedAreaSpec; //!
	70	fastjet::JetDefinition *fJetDef; //!
	71	fastjet::JetDefinition::Plugin *fPlugin; //!
	72	fastjet::RangeDefinition *fRange; //!
	73	fastjet::ClusterSequenceArea *fClustSeq; //!
	74	fastjet::Strategy fStrategy; //!
	75	fastjet::JetAlgorithm fAlgor; //!
	76	fastjet::RecombinationScheme fScheme; //!
	77	fastjet::AreaType fAreaType; //!
	78	Int_t fNGhostRepeats; //!
	79	Double_t fGhostArea; //!
	80	Double_t fMaxRap; //!
	81	Double_t fR; //!
48d874ff	82	// no setters for the moment - used default values in the constructor
5d0a5007	83	Double_t fGridScatter; //!
	84	Double_t fKtScatter; //!
	85	Double_t fMeanGhostKt; //!
	86	Int_t fPluginAlgor; //!
48d874ff	87	// extra parameters
5d0a5007	88	Double_t fMedUsedForBgSub; //!
	89	Bool_t fUseArea4Vector; //!
	90	#ifdef FASTJET_VERSION
	91	fastjet::JetMedianBackgroundEstimator *fBkrdEstimator; //!
	92	#endif
	93	Bool_t fLegacyMode; //!
48d874ff	94
	95	virtual void SubtractBackground(const Double_t median_pt = -1);
	96
	97	private:
	98	AliFJWrapper();
	99	AliFJWrapper(const AliFJWrapper& wrapper);
	100	AliFJWrapper& operator = (const AliFJWrapper& wrapper);
	101	};
	102	#endif
	103	#endif
	104
	105	#ifdef AliFJWrapper_CXX
	106	#undef AliFJWrapper_CXX
	107
	108	#if defined __GNUC__
	109	#pragma GCC system_header
	110	#endif
	111
	112	namespace fj = fastjet;
	113
	114	//_________________________________________________________________________________________________
	115	AliFJWrapper::AliFJWrapper(const char name, const char title)
a9c8e539	116	:
	117	fName (name)
	118	, fTitle (title)
48d874ff	119	, fInputVectors ( )
	120	, fInclusiveJets ( )
	121	, fSubtractedJetsPt ( )
	122	, fAreaDef (0)
	123	, fVorAreaSpec (0)
	124	, fGhostedAreaSpec (0)
	125	, fJetDef (0)
	126	, fPlugin (0)
	127	, fRange (0)
	128	, fClustSeq (0)
	129	, fStrategy (fj::Best)
	130	, fAlgor (fj::kt_algorithm)
	131	, fScheme (fj::BIpt_scheme)
	132	, fAreaType (fj::active_area)
	133	, fNGhostRepeats (1)
	134	, fGhostArea (0.01)
	135	, fMaxRap (1.)
	136	, fR (0.4)
	137	, fGridScatter (1.0)
7dc91d51	138	, fKtScatter (0.1)
48d874ff	139	, fMeanGhostKt (1e-100)
	140	, fPluginAlgor (0)
	141	, fMedUsedForBgSub (0)
	142	, fUseArea4Vector (kFALSE)
5d0a5007	143	#ifdef FASTJET_VERSION
	144	, fBkrdEstimator (0)
	145	#endif
	146	, fLegacyMode (false)
48d874ff	147	{
	148	// Constructor.
	149	}
	150
	151	//_________________________________________________________________________________________________
	152	AliFJWrapper::~AliFJWrapper()
	153	{
	154	// Destructor.
	155
	156	delete fAreaDef;
	157	delete fVorAreaSpec;
	158	delete fGhostedAreaSpec;
	159	delete fJetDef;
	160	delete fPlugin;
	161	delete fRange;
5d0a5007	162	delete fClustSeq;
	163	#ifdef FASTJET_VERSION
	164	if (fBkrdEstimator) delete fBkrdEstimator;
	165	#endif
48d874ff	166	}
	167
	168	//_________________________________________________________________________________________________
	169	void AliFJWrapper::CopySettingsFrom(const AliFJWrapper& wrapper)
	170	{
	171	// Copy some settings.
	172	// You very often want to keep most of the settings
	173	// but change only the algorithm or R - do it after call to this function
	174
	175	fStrategy = wrapper.fStrategy;
	176	fAlgor = wrapper.fAlgor;
	177	fScheme = wrapper.fScheme;
	178	fAreaType = wrapper.fAreaType;
	179	fNGhostRepeats = wrapper.fNGhostRepeats;
	180	fGhostArea = wrapper.fGhostArea;
	181	fMaxRap = wrapper.fMaxRap;
	182	fR = wrapper.fR;
	183	fGridScatter = wrapper.fGridScatter;
	184	fKtScatter = wrapper.fKtScatter;
	185	fMeanGhostKt = wrapper.fMeanGhostKt;
	186	fPluginAlgor = wrapper.fPluginAlgor;
	187	fUseArea4Vector = wrapper.fUseArea4Vector;
5d0a5007	188	fLegacyMode = wrapper.fLegacyMode;
48d874ff	189	}
	190
	191	//_________________________________________________________________________________________________
a9c8e539	192	void AliFJWrapper::Clear(const Option_t /opt*/)
48d874ff	193	{
	194	// Simply clear the input vectors.
	195	// Make sure done on every event if the instance is reused
	196	// Reset the median to zero.
	197
	198	fInputVectors.clear();
	199	fMedUsedForBgSub = 0;
	200
	201	// for the moment brute force delete everything
	202	delete fAreaDef; fAreaDef = 0;
	203	delete fVorAreaSpec; fVorAreaSpec = 0;
	204	delete fGhostedAreaSpec; fGhostedAreaSpec = 0;
	205	delete fJetDef; fJetDef = 0;
	206	delete fPlugin; fPlugin = 0;
	207	delete fRange; fRange = 0;
	208	delete fClustSeq; fClustSeq = 0;
5d0a5007	209	#ifdef FASTJET_VERSION
	210	if (fBkrdEstimator) delete fBkrdEstimator; fBkrdEstimator = 0;
	211	#endif
48d874ff	212	}
	213
	214	//_________________________________________________________________________________________________
	215	void AliFJWrapper::AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index)
	216	{
	217	// Make the input pseudojet.
	218
	219	fastjet::PseudoJet inVec(px, py, pz, E);
	220
	221	if (index > -99999) {
	222	inVec.set_user_index(index);
	223	} else {
	224	inVec.set_user_index(fInputVectors.size());
	225	}
	226
	227	// add to the fj container of input vectors
	228	fInputVectors.push_back(inVec);
	229	}
	230
	231	//_________________________________________________________________________________________________
	232	void AliFJWrapper::AddInputVector(const fj::PseudoJet& vec, Int_t index)
	233	{
	234	// Add an input pseudojet.
	235
	236	fj::PseudoJet inVec = vec;
	237
	238	if (index > -99999) {
	239	inVec.set_user_index(index);
	240	} else {
	241	inVec.set_user_index(fInputVectors.size());
	242	}
	243
	244	// add to the fj container of input vectors
	245	fInputVectors.push_back(inVec);
	246	}
	247
	248	//_________________________________________________________________________________________________
	249	void AliFJWrapper::AddInputVectors(const std::vector<fj::PseudoJet>& vecs, Int_t offsetIndex)
	250	{
	251	// Add the input from vector of pseudojets.
	252
	253	for (UInt_t i = 0; i < vecs.size(); ++i) {
	254	fj::PseudoJet inVec = vecs[i];
48d874ff	255	if (offsetIndex > -99999)
	256	inVec.set_user_index(fInputVectors.size() + offsetIndex);
	257	// add to the fj container of input vectors
48d874ff	258	fInputVectors.push_back(inVec);
	259	}
	260	}
	261
	262	//_________________________________________________________________________________________________
	263	Double_t AliFJWrapper::GetJetArea(UInt_t idx) const
	264	{
	265	// Get the jet area.
	266
	267	Double_t retval = -1; // really wrong area..
	268	if ( idx < fInclusiveJets.size() ) {
7dc91d51	269	retval = fClustSeq->area(fInclusiveJets[idx]);
48d874ff	270	} else {
	271	AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
	272	}
	273	return retval;
	274	}
	275
4e22c11c	276	//_________________________________________________________________________________________________
	277	fastjet::PseudoJet AliFJWrapper::GetJetAreaVector(UInt_t idx) const
	278	{
	279	// Get the jet area as vector.
	280	fastjet::PseudoJet retval;
	281	if ( idx < fInclusiveJets.size() ) {
	282	retval = fClustSeq->area_4vector(fInclusiveJets[idx]);
	283	} else {
	284	AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
	285	}
	286	return retval;
	287	}
	288
48d874ff	289	//_________________________________________________________________________________________________
	290	std::vector<double> AliFJWrapper::GetSubtractedJetsPts(Double_t median_pt, Bool_t sorted)
	291	{
	292	// Get subtracted jets pTs, returns vector.
	293
	294	SubtractBackground(median_pt);
	295
	296	if (kTRUE == sorted) {
	297	std::sort(fSubtractedJetsPt.begin(), fSubtractedJetsPt.begin());
	298	}
	299	return fSubtractedJetsPt;
	300	}
	301
	302	//_________________________________________________________________________________________________
	303	Double_t AliFJWrapper::GetJetSubtractedPt(UInt_t idx) const
	304	{
	305	// Get subtracted jets pTs, returns Double_t.
	306
	307	Double_t retval = -99999.; // really wrong pt..
	308	if ( idx < fSubtractedJetsPt.size() ) {
	309	retval = fSubtractedJetsPt[idx];
	310	}
	311	return retval;
	312	}
	313
	314	//_________________________________________________________________________________________________
	315	std::vector<fastjet::PseudoJet>
	316	AliFJWrapper::GetJetConstituents(UInt_t idx) const
	317	{
	318	// Get jets constituents.
	319
	320	std::vector<fastjet::PseudoJet> retval;
	321
	322	if ( idx < fInclusiveJets.size() ) {
	323	retval = fClustSeq->constituents(fInclusiveJets[idx]);
	324	} else {
	325	AliError(Form("[e] ::GetJetConstituents wrong index: %d",idx));
	326	}
	327
	328	return retval;
	329	}
	330
	331	//_________________________________________________________________________________________________
	332	void AliFJWrapper::GetMedianAndSigma(Double_t &median, Double_t &sigma, Int_t remove) const
	333	{
	334	// Get the median and sigma from fastjet.
	335	// User can also do it on his own because the cluster sequence is exposed (via a getter)
	336
	337	if (!fClustSeq) {
	338	AliError("[e] Run the jfinder first.");
b0a20945	339	return;
48d874ff	340	}
	341
	342	Double_t mean_area = 0;
	343	try {
	344	if(0 == remove) {
	345	fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
	346	} else {
	347	std::vector<fastjet::PseudoJet> input_jets = sorted_by_pt(fClustSeq->inclusive_jets());
	348	input_jets.erase(input_jets.begin(), input_jets.begin() + remove);
5d0a5007	349	fClustSeq->get_median_rho_and_sigma(input_jets, *fRange, fUseArea4Vector, median, sigma, mean_area);
48d874ff	350	input_jets.clear();
	351	}
	352	} catch (fj::Error) {
	353	AliError(" [w] FJ Exception caught.");
	354	median = -1.;
	355	sigma = -1;
	356	}
	357	}
	358
	359	//_________________________________________________________________________________________________
	360	Int_t AliFJWrapper::Run()
	361	{
	362	// Run the actual jet finder.
	363
	364	if (fAreaType == fj::voronoi_area) {
	365	// Rfact - check dependence - default is 1.
	366	// NOTE: hardcoded variable!
	367	fVorAreaSpec = new fj::VoronoiAreaSpec(1.);
	368	fAreaDef = new fj::AreaDefinition(*fVorAreaSpec);
	369	} else {
	370	fGhostedAreaSpec = new fj::GhostedAreaSpec(fMaxRap,
	371	fNGhostRepeats,
	372	fGhostArea,
	373	fGridScatter,
	374	fKtScatter,
	375	fMeanGhostKt);
	376
	377	fAreaDef = new fj::AreaDefinition(*fGhostedAreaSpec, fAreaType);
	378	}
	379
	380	// this is acceptable by fastjet:
	381	fRange = new fj::RangeDefinition(fMaxRap - 0.95 * fR);
	382
	383	if (fAlgor == fj::plugin_algorithm) {
	384	if (fPluginAlgor == 0) {
	385	// SIS CONE ALGOR
	386	// NOTE: hardcoded split parameter
	387	Double_t overlap_threshold = 0.75; // NOTE: this actually splits a lot: thr/min(pt1,pt2)
	388	fPlugin = new fj::SISConePlugin(fR,
	389	overlap_threshold,
	390	0, //search of stable cones - zero = until no more
	391	1.0); // this should be seed effectively for proto jets
	392	fJetDef = new fastjet::JetDefinition(fPlugin);
	393	} else {
	394	AliError("[e] Unrecognized plugin number!");
	395	}
	396	} else {
	397	fJetDef = new fj::JetDefinition(fAlgor, fR, fScheme, fStrategy);
	398	}
	399
	400	try {
	401	fClustSeq = new fj::ClusterSequenceArea(fInputVectors, fJetDef, fAreaDef);
	402	} catch (fj::Error) {
	403	AliError(" [w] FJ Exception caught.");
	404	return -1;
	405	}
	406
5d0a5007	407	// FJ3 :: Define an JetMedianBackgroundEstimator just in case it will be used
	408	#ifdef FASTJET_VERSION
	409	//fBkrdEstimator = new fj::JetMedianBackgroundEstimator(fRange, fJetDef, *fAreaDef) ;
	410	#endif
	411
	412	if (fLegacyMode) { SetLegacyFJ(); } // for FJ 2.x even if fLegacyMode is set, SetLegacyFJ is dummy
	413
48d874ff	414	// inclusive jets:
48d874ff	415	fInclusiveJets.clear();
7dc91d51	416	fInclusiveJets = fClustSeq->inclusive_jets(0.0);
48d874ff	417
	418	return 0;
	419	}
	420
5d0a5007	421	//_________________________________________________________________________________________________
	422	void AliFJWrapper::SetLegacyFJ()
	423	{
	424	// This methods enable legacy behaviour (FastJet 2.x) when AliROOT is compiled with FastJet 3.x
	425	#ifdef FASTJET_VERSION
	426	std::cout << "WARNING! Setting FastJet in legacy mode" << std::endl;
	427	if (fGhostedAreaSpec) { fGhostedAreaSpec->set_fj2_placement(kTRUE); }
	428	if (fBkrdEstimator) {
	429	fBkrdEstimator->set_provide_fj2_sigma(kTRUE);
	430	fBkrdEstimator->set_use_area_4vector(kFALSE);
	431	}
	432	#endif
	433	}
	434
48d874ff	435	//_________________________________________________________________________________________________
	436	void AliFJWrapper::SubtractBackground(Double_t median_pt)
	437	{
	438	// Subtract the background (specify the value - see below the meaning).
	439	// Negative argument means the bg will be determined with the current algorithm
	440	// this is the default behavior. Zero is allowed
	441	// Note: user may set the switch for area4vector based subtraction.
	442
	443	Double_t median = 0;
	444	Double_t sigma = 0;
	445	Double_t mean_area = 0;
	446
	447	// clear the subtracted jet pt's vector<double>
	448	fSubtractedJetsPt.clear();
5d0a5007	449
48d874ff	450	// check what was specified (default is -1)
	451	if (median_pt < 0) {
	452	try {
5d0a5007	453	fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
48d874ff	454	}
5d0a5007	455
48d874ff	456	catch (fj::Error) {
	457	AliError(" [w] FJ Exception caught.");
	458	median = -9999.;
	459	sigma = -1;
	460	fMedUsedForBgSub = median;
	461	return;
	462	}
	463	fMedUsedForBgSub = median;
	464	} else {
	465	// we do not know the sigma in this case
	466	sigma = -1;
	467	if (0.0 == median_pt) {
	468	// AliWarning(Form(" [w] Median specified for bg subtraction is ZERO: %f", median_pt )) << endl;
	469	fMedUsedForBgSub = 0.;
	470	} else {
	471	fMedUsedForBgSub = median_pt;
	472	}
	473	}
	474
	475	// subtract:
	476	for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
5d0a5007	477	if ( fUseArea4Vector ) {
48d874ff	478	// subtract the background using the area4vector
	479	fj::PseudoJet area4v = fClustSeq->area_4vector(fInclusiveJets[i]);
	480	fj::PseudoJet jet_sub = fInclusiveJets[i] - area4v * fMedUsedForBgSub;
	481	fSubtractedJetsPt.push_back(jet_sub.perp()); // here we put only the pt of the jet - note: this can be negative
	482	} else {
	483	// subtract the background using scalars
	484	// fj::PseudoJet jet_sub = fInclusiveJets[i] - area * fMedUsedForBgSub_;
	485	Double_t area = fClustSeq->area(fInclusiveJets[i]);
	486	// standard subtraction
	487	Double_t pt_sub = fInclusiveJets[i].perp() - fMedUsedForBgSub * area;
	488	fSubtractedJetsPt.push_back(pt_sub); // here we put only the pt of the jet - note: this can be negative
	489	}
	490	}
	491	}
	492
	493	//_________________________________________________________________________________________________
	494	void AliFJWrapper::SetupAlgorithmfromOpt(const char *option)
	495	{
	496	// Setup algorithm from char.
	497
	498	std::string opt(option);
	499
	500	if (!opt.compare("kt")) fAlgor = fj::kt_algorithm;
	501	if (!opt.compare("antikt")) fAlgor = fj::antikt_algorithm;
	502	if (!opt.compare("cambridge")) fAlgor = fj::cambridge_algorithm;
	503	if (!opt.compare("genkt")) fAlgor = fj::genkt_algorithm;
	504	if (!opt.compare("cambridge_passive")) fAlgor = fj::cambridge_for_passive_algorithm;
	505	if (!opt.compare("genkt_passive")) fAlgor = fj::genkt_for_passive_algorithm;
	506	if (!opt.compare("ee_kt")) fAlgor = fj::ee_kt_algorithm;
	507	if (!opt.compare("ee_genkt")) fAlgor = fj::ee_genkt_algorithm;
	508	if (!opt.compare("plugin")) fAlgor = fj::plugin_algorithm;
	509	}
	510
	511	//_________________________________________________________________________________________________
	512	void AliFJWrapper::SetupAreaTypefromOpt(const char *option)
	513	{
	514	// Setup area type from char.
	515
	516	std::string opt(option);
	517
	518	if (!opt.compare("active")) fAreaType = fj::active_area;
	519	if (!opt.compare("invalid")) fAreaType = fj::invalid_area;
	520	if (!opt.compare("active_area_explicit_ghosts")) fAreaType = fj::active_area_explicit_ghosts;
	521	if (!opt.compare("one_ghost_passive")) fAreaType = fj::one_ghost_passive_area;
	522	if (!opt.compare("passive")) fAreaType = fj::passive_area;
	523	if (!opt.compare("voronoi")) fAreaType = fj::voronoi_area;
	524	}
	525
	526	//_________________________________________________________________________________________________
	527	void AliFJWrapper::SetupSchemefromOpt(const char *option)
	528	{
	529	//
	530	// setup scheme from char
	531	//
	532
	533	std::string opt(option);
	534
	535	if (!opt.compare("BIpt")) fScheme = fj::BIpt_scheme;
	536	if (!opt.compare("BIpt2")) fScheme = fj::BIpt2_scheme;
	537	if (!opt.compare("E")) fScheme = fj::E_scheme;
	538	if (!opt.compare("pt")) fScheme = fj::pt_scheme;
	539	if (!opt.compare("pt2")) fScheme = fj::pt2_scheme;
	540	if (!opt.compare("Et")) fScheme = fj::Et_scheme;
	541	if (!opt.compare("Et2")) fScheme = fj::Et2_scheme;
542	}
543
544	//_________________________________________________________________________________________________
545	void AliFJWrapper::SetupStrategyfromOpt(const char *option)
546	{
547	// Setup strategy from char.
548
549	std::string opt(option);
550
551	if (!opt.compare("Best")) fStrategy = fj::Best;
552	if (!opt.compare("N2MinHeapTiled")) fStrategy = fj::N2MinHeapTiled;
553	if (!opt.compare("N2Tiled")) fStrategy = fj::N2Tiled;
554	if (!opt.compare("N2PoorTiled")) fStrategy = fj::N2PoorTiled;
555	if (!opt.compare("N2Plain")) fStrategy = fj::N2Plain;
556	if (!opt.compare("N3Dumb")) fStrategy = fj::N3Dumb;
557	if (!opt.compare("NlnN")) fStrategy = fj::NlnN;
558	if (!opt.compare("NlnN3pi")) fStrategy = fj::NlnN3pi;
559	if (!opt.compare("NlnN4pi")) fStrategy = fj::NlnN4pi;
560	if (!opt.compare("NlnNCam4pi")) fStrategy = fj::NlnNCam4pi;
561	if (!opt.compare("NlnNCam2pi2R")) fStrategy = fj::NlnNCam2pi2R;
562	if (!opt.compare("NlnNCam")) fStrategy = fj::NlnNCam;
563	if (!opt.compare("plugin")) fStrategy = fj::plugin_strategy;
564	}
565	#endif