[u/mrichter/AliRoot.git] / PWGGA / 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;
  Double_t                              GetJetSubtractedPt (UInt_t idx) const;
  virtual std::vector<double>           GetSubtractedJetsPts(Double_t median_pt = -1, Bool_t sorted = kFALSE);

  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);

 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;   //!

  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)
{
  // Constructor.
}

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

  delete fAreaDef;
  delete fVorAreaSpec;
  delete fGhostedAreaSpec;
  delete fJetDef;
  delete fPlugin;
  delete fRange;
  delete fClustSeq;  
}

//_________________________________________________________________________________________________
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;
}

//_________________________________________________________________________________________________
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;
}

//_________________________________________________________________________________________________
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;
}

//_________________________________________________________________________________________________
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.");
  }

  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, kFALSE, 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;
  }

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

  return 0;
}

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