1 #include "AliFJWrapper.h"
6 #include "AliFastJetHeader.h"
20 namespace fj = fastjet;
22 AliFJWrapper::AliFJWrapper(const char *name, const char *title)
28 , fSubtractedJetsPt ( )
32 , fGhostedAreaSpec (0)
38 , fStrategy (fj::Best)
39 , fAlgor (fj::kt_algorithm)
40 , fScheme (fj::BIpt_scheme)
41 , fAreaType (fj::active_area)
50 , fMeanGhostKt (1e-100)
53 , fMedianUsedForBgSubtraction (0)
54 , fUseArea4Vector (false)
59 void AliFJWrapper::CopySettingsFrom (const AliFJWrapper& wrapper)
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
67 fStrategy = wrapper.fStrategy;
68 fAlgor = wrapper.fAlgor;
69 fScheme = wrapper.fScheme;
70 fAreaType = wrapper.fAreaType;
72 fNGhostRepeats = wrapper.fNGhostRepeats;
73 fGhostArea = wrapper.fGhostArea;
74 fMaxRap = wrapper.fMaxRap;
77 fGridScatter = wrapper.fGridScatter;
78 fKtScatter = wrapper.fKtScatter;
79 fMeanGhostKt = wrapper.fMeanGhostKt;
80 fPluginAlgor = wrapper.fPluginAlgor;
82 fUseArea4Vector = wrapper.fUseArea4Vector;
86 AliFJWrapper::~AliFJWrapper()
93 delete fGhostedAreaSpec;
100 void AliFJWrapper::Clear(const Option_t *opt)
103 // simply clear the input vectors
104 // make sure done on every event if the instance is reused
105 // reset the median to zero
108 fInputVectors.clear();
109 fMedianUsedForBgSubtraction = 0;
111 // for the moment brute force delete everything
114 delete fGhostedAreaSpec;
123 void AliFJWrapper::SetupFromHeader(const AliFastJetHeader *header)
126 // feed the wrapper from the parameters of the header
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();
140 void AliFJWrapper::AddInputVector(double px, double py, double pz, double E, int index)
142 // make the input pseudojet
143 fastjet::PseudoJet inVec(px, py, pz, E);
147 inVec.set_user_index(index);
151 inVec.set_user_index(fInputVectors.size());
154 // add to the fj container of input vectors
155 fInputVectors.push_back(inVec);
158 void AliFJWrapper::AddInputVector(const fj::PseudoJet& vec, int index)
160 // add an input pseudojet
161 fj::PseudoJet inVec = vec;
165 inVec.set_user_index(index);
169 inVec.set_user_index(fInputVectors.size());
172 // add to the fj container of input vectors
173 fInputVectors.push_back(inVec);
176 void AliFJWrapper::AddInputVectors(const vector<fastjet::PseudoJet>& vecs, int offsetIndex)
179 // add the input from vector of pseudojets
182 for (unsigned int i = 0; i < vecs.size(); i++)
184 fj::PseudoJet inVec = vecs[i];
186 if (offsetIndex > -99999)
187 inVec.set_user_index(fInputVectors.size() + offsetIndex);
188 // add to the fj container of input vectors
190 fInputVectors.push_back(inVec);
194 void AliFJWrapper::GetMedianAndSigma(double &median, double &sigma, int remove)
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)
203 AliError("[e] Run the jfinder first.");
206 double mean_area = 0;
211 fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
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);
224 AliError(" [w] FJ Exception caught.");
230 int AliFJWrapper::Run()
233 // run the actual jet finder
236 if (fAreaType == fj::voronoi_area)
238 // Rfact - check dependence - default is 1.
239 // NOTE: hardcoded variable!
240 fVorAreaSpec = new fj::VoronoiAreaSpec(1.);
241 fAreaDef = new fj::AreaDefinition(*fVorAreaSpec);
245 fGhostedAreaSpec = new fj::GhostedAreaSpec(fMaxRap,
252 fAreaDef = new fj::AreaDefinition(*fGhostedAreaSpec, fAreaType);
255 // this is acceptable by fastjet:
256 fRange = new fj::RangeDefinition(fMaxRap - 0.95 * fR);
258 if (fAlgor == fj::plugin_algorithm)
260 if (fPluginAlgor == 0)
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,
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");
275 AliError("[e] Unrecognized plugin number!");
280 //fJetDef = new fj::JetDefinition(fj::kt_algorithm, 0.4, fj::pt_scheme, fj::Best);
281 fJetDef = new fj::JetDefinition(fAlgor, fR, fScheme, fStrategy);
286 fClustSeq = new fj::ClusterSequenceArea(fInputVectors, *fJetDef, *fAreaDef);
290 AliError(" [w] FJ Exception caught.");
295 fInclusiveJets.clear();
296 fInclusiveJets = fClustSeq->inclusive_jets(0.0); //sorted_by_pt(fClustSeq->inclusive_jets(0.0));
301 void AliFJWrapper::SubtractBackground(const double median_pt)
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
308 // Note: user may set the switch for area4vector based subtraction
313 double mean_area = 0;
315 // clear the subtracted jet pt's vector<double>
316 fSubtractedJetsPt.clear();
318 // check what was specified (default is -1)
323 fClustSeq->get_median_rho_and_sigma(*fRange, false, median, sigma, mean_area);
328 AliError(" [w] FJ Exception caught.");
331 fMedianUsedForBgSubtraction = median;
334 fMedianUsedForBgSubtraction = median;
338 // we do not know the sigma in this case
340 if (0.0 == median_pt)
342 // AliWarning(Form(" [w] Median specified for bg subtraction is ZERO: %f", median_pt )) << endl;
343 fMedianUsedForBgSubtraction = 0.;
347 fMedianUsedForBgSubtraction = median_pt;
352 for (unsigned i = 0; i < fInclusiveJets.size(); i++)
354 //fj::PseudoJet jet_sub = fClustSeq->subtracted_jet(fInclusiveJets[i], fMedianUsedForBgSubtraction_);
355 //cout << fInclusiveJets[i].perp() << " (" << fMedianUsedForBgSubtraction_ << ") -> " << jet_sub.perp() << endl;
357 if (true == fUseArea4Vector)
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
366 // subtract the background using scalars
367 // fj::PseudoJet jet_sub = fInclusiveJets[i] - area * fMedianUsedForBgSubtraction_;
369 double area = fClustSeq->area(fInclusiveJets[i]);
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;
379 double AliFJWrapper::GetJetArea(const unsigned int idx)
384 double retval = -1; // really wrong area..
385 if ( idx < fInclusiveJets.size() )
387 retval = fClustSeq->area(fInclusiveJets[idx]);
391 cerr << "[e] ::GetJetArea wrong index: " << idx << endl;
396 std::vector<double> AliFJWrapper::GetSubtractedJetsPts(const double median_pt, const bool sorted)
399 // get subtracted jets' pTs returns vector
401 SubtractBackground(median_pt);
405 std::sort(fSubtractedJetsPt.begin(), fSubtractedJetsPt.begin());
407 return fSubtractedJetsPt;
410 double AliFJWrapper::GetJetSubtractedPt(const unsigned int idx)
413 // get subtracted jets' pTs returns double (pT)
415 double retval = -99999.; // really wrong pt..
416 if ( idx < fSubtractedJetsPt.size() )
418 retval = fSubtractedJetsPt[idx];
423 std::vector<fastjet::PseudoJet>
424 AliFJWrapper::GetJetConstituents(const unsigned int idx)
427 // get jets' constituents
429 std::vector<fastjet::PseudoJet> retval;
431 if ( idx < fInclusiveJets.size() )
433 retval = fClustSeq->constituents(fInclusiveJets[idx]);
437 cerr << "[e] ::GetJetConstituents wrong index: " << idx << endl;
443 void AliFJWrapper::SetupAlgorithmfromOpt(const char *option)
446 // setup algorithm from char
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;
462 void AliFJWrapper::SetupStrategyfromOpt(const char *option)
465 // setup strategy from char
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;
484 void AliFJWrapper::SetupSchemefromOpt(const char *option)
487 // setup scheme from char
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;
500 void AliFJWrapper::SetupAreaTypefromOpt(const char *option)
503 // setup area type from char
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;