HLT/FJWrapper/AliFJWrapper.cxx

   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 }