[u/mrichter/AliRoot.git] / HLT / FJWrapper / AliFJWrapper.cxx

#include "AliFJWrapper.h"

#include "TString.h"

#include "AliLog.h"
#include "AliFastJetHeader.h"

#include <cmath>
#include <iostream>
/*
@
@
@
@
@
@
@
 */
using namespace std;
namespace fj = fastjet;

AliFJWrapper::AliFJWrapper(const char *name, const char *title)
  : TNamed(name, 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         (1.0)
  , fMeanGhostKt       (1e-100)
  , fPluginAlgor       (0)

  , fMedianUsedForBgSubtraction (0)
  , fUseArea4Vector    (false)
{
  ;
}

void AliFJWrapper::CopySettingsFrom (const AliFJWrapper& wrapper)
{
  //
  // copy some settings - avoid double typing...
  // 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;
  
}

AliFJWrapper::~AliFJWrapper()
{  
  //
  // destructor
  //
  delete fAreaDef;
  delete fVorAreaSpec;
  delete fGhostedAreaSpec;
  delete fJetDef;
  delete fPlugin;
  delete fRange;
  delete fClustSeq;  
}

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();
  fMedianUsedForBgSubtraction = 0;

  // for the moment brute force delete everything
  delete fAreaDef;
  delete fVorAreaSpec;
  delete fGhostedAreaSpec;
  delete fJetDef;
  delete fPlugin;
  delete fRange;
  delete fClustSeq;  

  TNamed::Clear(opt);
}

void AliFJWrapper::SetupFromHeader(const AliFastJetHeader *header)
{
  //
  // feed the wrapper from the parameters of the header
  //
  
  fR             = header->GetRparam();
  fStrategy      = header->GetStrategy();
  fScheme        = header->GetRecombScheme();
  fAlgor         = header->GetAlgorithm(); 
  fAreaType      = header->GetAreaType();
  fMaxRap        = header->GetGhostEtaMax(); 
  fGhostArea     = header->GetGhostArea(); 
  fNGhostRepeats = header->GetActiveAreaRepeats(); 

}

void AliFJWrapper::AddInputVector(double px, double py, double pz, double E, int 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 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 vector<fastjet::PseudoJet>& vecs, int offsetIndex)
{
  //
  // add the input from vector of pseudojets
  //

  for (unsigned int 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);
    }
}

void AliFJWrapper::GetMedianAndSigma(double &median, double &sigma, int remove)
{
  //
  // 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 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, false, median, sigma, mean_area);
	  input_jets.clear();
	}
    }
  
  catch (fj::Error)
    {
      AliError(" [w] FJ Exception caught.");
      median = -1.;
      sigma = -1;
    }
}

int 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 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
	  //__INFO("Define siscone");
	  fJetDef = new fastjet::JetDefinition(fPlugin);
	  //__INFO("SIS cone initialized");
	}
      else
	{
	  AliError("[e] Unrecognized plugin number!");
	}
    }
  else
    {
      //fJetDef = new fj::JetDefinition(fj::kt_algorithm, 0.4, fj::pt_scheme, fj::Best);
      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); //sorted_by_pt(fClustSeq->inclusive_jets(0.0));

  return 0;
}

void AliFJWrapper::SubtractBackground(const double 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 median = 0;
  double sigma = 0;
  double 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, false, median, sigma, mean_area);
	}
      
      catch (fj::Error)
	{
	  AliError(" [w] FJ Exception caught.");
	  median = -9999.;
	  sigma = -1;
	  fMedianUsedForBgSubtraction = median;
	  return;
	}
      fMedianUsedForBgSubtraction = 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;
	  fMedianUsedForBgSubtraction = 0.;
	}
      else
	{
	  fMedianUsedForBgSubtraction = median_pt;
	}
    }

  // subtract:
  for (unsigned i = 0; i < fInclusiveJets.size(); i++) 
    {
      //fj::PseudoJet jet_sub = fClustSeq->subtracted_jet(fInclusiveJets[i], fMedianUsedForBgSubtraction_);
      //cout << fInclusiveJets[i].perp() << " (" << fMedianUsedForBgSubtraction_ << ") -> " << jet_sub.perp() << endl;
      
      if (true == fUseArea4Vector)
	{
	  // sutract the background using the area4vector
	  fj::PseudoJet area4v = fClustSeq->area_4vector(fInclusiveJets[i]);
	  fj::PseudoJet jet_sub = fInclusiveJets[i] - area4v * fMedianUsedForBgSubtraction;
	  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 * fMedianUsedForBgSubtraction_;

	  double area = fClustSeq->area(fInclusiveJets[i]);

	  // standard subtraction
	  double pt_sub = fInclusiveJets[i].perp() - fMedianUsedForBgSubtraction * area;
	  fSubtractedJetsPt.push_back(pt_sub); // here we put only the pt of the jet - note: this can be negative
	  //cout << "phi incl : " << fInclusiveJets[i].phi() << " subtracted : " << jet_sub.phi() << endl;
	}
    }
}

double AliFJWrapper::GetJetArea(const unsigned int idx)
{
  //
  // get the jet area
  //
  double retval = -1; // really wrong area..
  if ( idx < fInclusiveJets.size() )
    {
      retval = fClustSeq->area(fInclusiveJets[idx]);
    }
  else
    {
      cerr << "[e] ::GetJetArea wrong index: " << idx << endl;
    }
  return retval;
}

std::vector<double> AliFJWrapper::GetSubtractedJetsPts(const double median_pt, const bool sorted)
{ 
  //
  // get subtracted jets' pTs returns vector
  //
  SubtractBackground(median_pt);
  
  if (true == sorted)
    {
      std::sort(fSubtractedJetsPt.begin(), fSubtractedJetsPt.begin());
    }
  return fSubtractedJetsPt;
}

double AliFJWrapper::GetJetSubtractedPt(const unsigned int idx)
{
  //
  // get subtracted jets' pTs returns double (pT)
  //
  double retval = -99999.; // really wrong pt..
  if ( idx < fSubtractedJetsPt.size() )
    {
      retval = fSubtractedJetsPt[idx];
    }
  return retval;
}

std::vector<fastjet::PseudoJet> 
AliFJWrapper::GetJetConstituents(const unsigned int idx)
{
  //
  // get jets' constituents
  //
  std::vector<fastjet::PseudoJet> retval;
  
  if ( idx < fInclusiveJets.size() )
    {
      retval = fClustSeq->constituents(fInclusiveJets[idx]);
    }  
  else
    {
      cerr << "[e] ::GetJetConstituents wrong index: " << idx << endl;
    }
  
  return retval;
}

void AliFJWrapper::SetupAlgorithmfromOpt(const char *option)
{
  //
  // setup algorithm from char
  //
  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::SetupStrategyfromOpt(const char *option)
{
  //
  // setup strategy from char
  //
  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;

}
void AliFJWrapper::SetupSchemefromOpt(const char *option)
{
  //
  // setup scheme from char
  //
  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::SetupAreaTypefromOpt(const char *option)
{
  //
  // setup area type from char
  //
  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;

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