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

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

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