PWGJE/EMCALJetTasks/AliFJWrapper.h

   1 #ifndef AliFJWrapper_H
   2 #define AliFJWrapper_H
   3
   4 // $Id$
   5
   6 #if !defined(__CINT__)
   7
   8 #include <vector>
   9 #include <TString.h>
  10 #include "AliLog.h"
  11 #include "FJ_includes.h"
  12 #include "AliJetShape.h"
  13
  14 class AliFJWrapper
  15 {
  16  public:
  17   AliFJWrapper(const char *name, const char *title);
  18   virtual ~AliFJWrapper();
  19
  20   virtual void  AddInputVector (Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index = -99999);
  21   virtual void  AddInputVector (const fastjet::PseudoJet& vec,                Int_t index = -99999);
  22   virtual void  AddInputVectors(const std::vector<fastjet::PseudoJet>& vecs,  Int_t offsetIndex = -99999);
  23   virtual const char *ClassName()                            const { return "AliFJWrapper";              }
  24   virtual void  Clear(const Option_t* /*opt*/ = "");
  25   virtual void  CopySettingsFrom (const AliFJWrapper& wrapper);
  26   virtual void  GetMedianAndSigma(Double_t& median, Double_t& sigma, Int_t remove = 0) const;
  27   fastjet::ClusterSequenceArea*           GetClusterSequence() const { return fClustSeq;                   }
  28   const std::vector<fastjet::PseudoJet>   GetInputVectors()    const { return fInputVectors;               }
  29   const std::vector<fastjet::PseudoJet>   GetInclusiveJets()   const { return fInclusiveJets;              }
  30   std::vector<fastjet::PseudoJet>         GetJetConstituents(UInt_t idx) const;
  31   Double_t                                GetMedianUsedForBgSubtraction() const { return fMedUsedForBgSub; }
  32   const char*                             GetName()            const { return fName;                       }
  33   const char*                             GetTitle()           const { return fTitle;                      }
  34   Double_t                                GetJetArea         (UInt_t idx) const;
  35   fastjet::PseudoJet                      GetJetAreaVector   (UInt_t idx) const;
  36   Double_t                                GetJetSubtractedPt (UInt_t idx) const;
  37   virtual std::vector<double>             GetSubtractedJetsPts(Double_t median_pt = -1, Bool_t sorted = kFALSE);
  38   Bool_t                                  GetLegacyMode()            { return fLegacyMode; }
  39 #ifdef FASTJET_CONTRIB
  40   const std::vector<fastjet::contrib::GenericSubtractorInfo> GetGenSubtractorInfoJetMass() const {return fGenSubtractorInfoJetMass;}
  41   const std::vector<fastjet::PseudoJet>   GetConstituentSubtrJets() const { return fConstituentSubtrJets;  }
  42 #endif
  43
  44   virtual Int_t Run();
  45   virtual Int_t DoGenericSubtractionJetMass();
  46   virtual Int_t DoConstituentSubtraction();
  47
  48   void SetStrategy(const fastjet::Strategy &strat)                 { fStrategy = strat;  }
  49   void SetAlgorithm(const fastjet::JetAlgorithm &algor)            { fAlgor    = algor;  }
  50   void SetRecombScheme(const fastjet::RecombinationScheme &scheme) { fScheme   = scheme; }
  51   void SetAreaType(const fastjet::AreaType &atype)                 { fAreaType = atype;  }
  52   void SetNRepeats(Int_t nrepeat)       { fNGhostRepeats  = nrepeat; }
  53   void SetGhostArea(Double_t gharea)    { fGhostArea      = gharea;  }
  54   void SetMaxRap(Double_t maxrap)       { fMaxRap         = maxrap;  }
  55   void SetR(Double_t r)                 { fR              = r;       }
  56   void SetGridScatter(Double_t gridSc)  { fGridScatter    = gridSc;  }
  57   void SetKtScatter(Double_t ktSc)      { fKtScatter      = ktSc;    }
  58   void SetMeanGhostKt(Double_t meankt)  { fMeanGhostKt    = meankt;  }
  59   void SetPluginAlgor(Int_t plugin)     { fPluginAlgor    = plugin;  }
  60   void SetUseArea4Vector(Bool_t useA4v) { fUseArea4Vector = useA4v;  }
  61   void SetupAlgorithmfromOpt(const char *option);
  62   void SetupAreaTypefromOpt(const char *option);
  63   void SetupSchemefromOpt(const char *option);
  64   void SetupStrategyfromOpt(const char *option);
  65   void SetLegacyMode (Bool_t mode)      { fLegacyMode ^= mode; }
  66   void SetLegacyFJ();
  67   void SetUseExternalBkg(Bool_t b, Double_t rho, Double_t rhom) { fUseExternalBkg = b; fRho = rho; fRhom = rhom;}
  68
  69  protected:
  70   TString                                fName;               //!
  71   TString                                fTitle;              //!
  72   std::vector<fastjet::PseudoJet>        fInputVectors;       //!
  73   std::vector<fastjet::PseudoJet>        fInclusiveJets;      //!
  74   std::vector<double>                    fSubtractedJetsPt;   //!
  75   std::vector<fastjet::PseudoJet>        fConstituentSubtrJets; //!
  76   fastjet::AreaDefinition               *fAreaDef;            //!
  77   fastjet::VoronoiAreaSpec              *fVorAreaSpec;        //!
  78   fastjet::GhostedAreaSpec              *fGhostedAreaSpec;    //!
  79   fastjet::JetDefinition                *fJetDef;             //!
  80   fastjet::JetDefinition::Plugin        *fPlugin;             //!
  81   fastjet::RangeDefinition              *fRange;              //!
  82   fastjet::ClusterSequenceArea          *fClustSeq;           //!
  83   fastjet::Strategy                      fStrategy;           //!
  84   fastjet::JetAlgorithm                  fAlgor;              //!
  85   fastjet::RecombinationScheme           fScheme;             //!
  86   fastjet::AreaType                      fAreaType;           //!
  87   Int_t                                  fNGhostRepeats;      //!
  88   Double_t                               fGhostArea;          //!
  89   Double_t                               fMaxRap;             //!
  90   Double_t                               fR;                  //!
  91   // no setters for the moment - used default values in the constructor
  92   Double_t                               fGridScatter;        //!
  93   Double_t                               fKtScatter;          //!
  94   Double_t                               fMeanGhostKt;        //!
  95   Int_t                                  fPluginAlgor;        //!
  96   // extra parameters
  97   Double_t                               fMedUsedForBgSub;    //!
  98   Bool_t                                 fUseArea4Vector;     //!
  99 #ifdef FASTJET_VERSION
 100   fastjet::JetMedianBackgroundEstimator   *fBkrdEstimator;        //!
 101 #endif
 102 #ifdef FASTJET_CONTRIB
 103   fastjet::contrib::GenericSubtractor     *fGenSubtractor;        //!
 104   std::vector<fastjet::contrib::GenericSubtractorInfo> fGenSubtractorInfoJetMass;    //!
 105 #endif
 106   Bool_t                                   fLegacyMode;           //!
 107   Bool_t                                   fUseExternalBkg;       //!
 108   Double_t                                 fRho;                  //  pT background density
 109   Double_t                                 fRhom;                 //  mT background density
 110
 111   virtual void   SubtractBackground(const Double_t median_pt = -1);
 112
 113  private:
 114   AliFJWrapper();
 115   AliFJWrapper(const AliFJWrapper& wrapper);
 116   AliFJWrapper& operator = (const AliFJWrapper& wrapper);
 117 };
 118 #endif
 119 #endif
 120
 121 #ifdef AliFJWrapper_CXX
 122 #undef AliFJWrapper_CXX
 123
 124 #if defined __GNUC__
 125 #pragma GCC system_header
 126 #endif
 127
 128 namespace fj = fastjet;
 129
 130 //_________________________________________________________________________________________________
 131 AliFJWrapper::AliFJWrapper(const char *name, const char *title)
 132   :
 133     fName              (name)
 134   , fTitle             (title)
 135   , fInputVectors      ( )
 136   , fInclusiveJets     ( )
 137   , fSubtractedJetsPt  ( )
 138   , fConstituentSubtrJets ( )
 139   , fAreaDef           (0)
 140   , fVorAreaSpec       (0)
 141   , fGhostedAreaSpec   (0)
 142   , fJetDef            (0)
 143   , fPlugin            (0)
 144   , fRange             (0)
 145   , fClustSeq          (0)
 146   , fStrategy          (fj::Best)
 147   , fAlgor             (fj::kt_algorithm)
 148   , fScheme            (fj::BIpt_scheme)
 149   , fAreaType          (fj::active_area)
 150   , fNGhostRepeats     (1)
 151   , fGhostArea         (0.005)
 152   , fMaxRap            (1.)
 153   , fR                 (0.4)
 154   , fGridScatter       (1.0)
 155   , fKtScatter         (0.1)
 156   , fMeanGhostKt       (1e-100)
 157   , fPluginAlgor       (0)
 158   , fMedUsedForBgSub   (0)
 159   , fUseArea4Vector    (kFALSE)
 160 #ifdef FASTJET_VERSION
 161   , fBkrdEstimator     (0)
 162 #endif
 163 #ifdef FASTJET_CONTRIB
 164   , fGenSubtractor     (0)
 165   , fGenSubtractorInfoJetMass ( )
 166 #endif
 167   , fLegacyMode        (false)
 168   , fUseExternalBkg    (false)
 169   , fRho               (0)
 170   , fRhom              (0)
 171 {
 172   // Constructor.
 173 }
 174
 175 //_________________________________________________________________________________________________
 176 AliFJWrapper::~AliFJWrapper()
 177 {
 178   // Destructor.
 179
 180   delete fAreaDef;
 181   delete fVorAreaSpec;
 182   delete fGhostedAreaSpec;
 183   delete fJetDef;
 184   delete fPlugin;
 185   delete fRange;
 186   delete fClustSeq;
 187 #ifdef FASTJET_VERSION
 188   if (fBkrdEstimator)     delete fBkrdEstimator;
 189 #endif
 190 #ifdef FASTJET_CONTRIB
 191   if (fGenSubtractor)     delete fGenSubtractor;
 192 #endif
 193 }
 194
 195 //_________________________________________________________________________________________________
 196 void AliFJWrapper::CopySettingsFrom(const AliFJWrapper& wrapper)
 197 {
 198   // Copy some settings.
 199   // You very often want to keep most of the settings
 200   // but change only the algorithm or R - do it after call to this function
 201
 202   fStrategy         = wrapper.fStrategy;
 203   fAlgor            = wrapper.fAlgor;
 204   fScheme           = wrapper.fScheme;
 205   fAreaType         = wrapper.fAreaType;
 206   fNGhostRepeats    = wrapper.fNGhostRepeats;
 207   fGhostArea        = wrapper.fGhostArea;
 208   fMaxRap           = wrapper.fMaxRap;
 209   fR                = wrapper.fR;
 210   fGridScatter      = wrapper.fGridScatter;
 211   fKtScatter        = wrapper.fKtScatter;
 212   fMeanGhostKt      = wrapper.fMeanGhostKt;
 213   fPluginAlgor      = wrapper.fPluginAlgor;
 214   fUseArea4Vector   = wrapper.fUseArea4Vector;
 215   fLegacyMode       = wrapper.fLegacyMode;
 216   fUseExternalBkg   = wrapper.fUseExternalBkg;
 217   fRho              = wrapper.fRho;
 218   fRhom             = wrapper.fRhom;
 219 }
 220
 221 //_________________________________________________________________________________________________
 222 void AliFJWrapper::Clear(const Option_t */*opt*/)
 223 {
 224   // Simply clear the input vectors.
 225   // Make sure done on every event if the instance is reused
 226   // Reset the median to zero.
 227
 228   fInputVectors.clear();
 229   fMedUsedForBgSub = 0;
 230
 231   // for the moment brute force delete everything
 232   delete fAreaDef;         fAreaDef         = 0;
 233   delete fVorAreaSpec;     fVorAreaSpec     = 0;
 234   delete fGhostedAreaSpec; fGhostedAreaSpec = 0;
 235   delete fJetDef;          fJetDef          = 0;
 236   delete fPlugin;          fPlugin          = 0;
 237   delete fRange;           fRange           = 0;
 238   delete fClustSeq;        fClustSeq        = 0;
 239 #ifdef FASTJET_VERSION
 240   if (fBkrdEstimator)     delete fBkrdEstimator     ;  fBkrdEstimator     = 0;
 241 #endif
 242 #ifdef FASTJET_CONTRIB
 243   if (fGenSubtractor)     delete fGenSubtractor     ;  fGenSubtractor     = 0;
 244 #endif
 245 }
 246
 247 //_________________________________________________________________________________________________
 248 void AliFJWrapper::AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index)
 249 {
 250   // Make the input pseudojet.
 251
 252   fastjet::PseudoJet inVec(px, py, pz, E);
 253
 254   if (index > -99999) {
 255     inVec.set_user_index(index);
 256   } else {
 257     inVec.set_user_index(fInputVectors.size());
 258   }
 259
 260   // add to the fj container of input vectors
 261   fInputVectors.push_back(inVec);
 262 }
 263
 264 //_________________________________________________________________________________________________
 265 void AliFJWrapper::AddInputVector(const fj::PseudoJet& vec, Int_t index)
 266 {
 267   // Add an input pseudojet.
 268
 269   fj::PseudoJet inVec = vec;
 270
 271   if (index > -99999) {
 272     inVec.set_user_index(index);
 273   } else {
 274     inVec.set_user_index(fInputVectors.size());
 275   }
 276
 277   // add to the fj container of input vectors
 278   fInputVectors.push_back(inVec);
 279 }
 280
 281 //_________________________________________________________________________________________________
 282 void AliFJWrapper::AddInputVectors(const std::vector<fj::PseudoJet>& vecs, Int_t offsetIndex)
 283 {
 284   // Add the input from vector of pseudojets.
 285
 286   for (UInt_t i = 0; i < vecs.size(); ++i) {
 287     fj::PseudoJet inVec = vecs[i];
 288     if (offsetIndex > -99999)
 289       inVec.set_user_index(fInputVectors.size() + offsetIndex);
 290     // add to the fj container of input vectors
 291     fInputVectors.push_back(inVec);
 292   }
 293 }
 294
 295 //_________________________________________________________________________________________________
 296 Double_t AliFJWrapper::GetJetArea(UInt_t idx) const
 297 {
 298   // Get the jet area.
 299
 300   Double_t retval = -1; // really wrong area..
 301   if ( idx < fInclusiveJets.size() ) {
 302     retval = fClustSeq->area(fInclusiveJets[idx]);
 303   } else {
 304     AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
 305   }
 306   return retval;
 307 }
 308
 309 //_________________________________________________________________________________________________
 310 fastjet::PseudoJet AliFJWrapper::GetJetAreaVector(UInt_t idx) const
 311 {
 312   // Get the jet area as vector.
 313   fastjet::PseudoJet retval;
 314   if ( idx < fInclusiveJets.size() ) {
 315     retval = fClustSeq->area_4vector(fInclusiveJets[idx]);
 316   } else {
 317     AliError(Form("[e] ::GetJetArea wrong index: %d",idx));
 318   }
 319   return retval;
 320 }
 321
 322 //_________________________________________________________________________________________________
 323 std::vector<double> AliFJWrapper::GetSubtractedJetsPts(Double_t median_pt, Bool_t sorted)
 324 {
 325   // Get subtracted jets pTs, returns vector.
 326
 327   SubtractBackground(median_pt);
 328
 329   if (kTRUE == sorted) {
 330     std::sort(fSubtractedJetsPt.begin(), fSubtractedJetsPt.begin());
 331   }
 332   return fSubtractedJetsPt;
 333 }
 334
 335 //_________________________________________________________________________________________________
 336 Double_t AliFJWrapper::GetJetSubtractedPt(UInt_t idx) const
 337 {
 338   // Get subtracted jets pTs, returns Double_t.
 339
 340   Double_t retval = -99999.; // really wrong pt..
 341   if ( idx < fSubtractedJetsPt.size() ) {
 342     retval = fSubtractedJetsPt[idx];
 343   }
 344   return retval;
 345 }
 346
 347 //_________________________________________________________________________________________________
 348 std::vector<fastjet::PseudoJet>
 349 AliFJWrapper::GetJetConstituents(UInt_t idx) const
 350 {
 351   // Get jets constituents.
 352
 353   std::vector<fastjet::PseudoJet> retval;
 354
 355   if ( idx < fInclusiveJets.size() ) {
 356     retval = fClustSeq->constituents(fInclusiveJets[idx]);
 357   } else {
 358     AliError(Form("[e] ::GetJetConstituents wrong index: %d",idx));
 359   }
 360
 361   return retval;
 362 }
 363
 364 //_________________________________________________________________________________________________
 365 void AliFJWrapper::GetMedianAndSigma(Double_t &median, Double_t &sigma, Int_t remove) const
 366 {
 367   // Get the median and sigma from fastjet.
 368   // User can also do it on his own because the cluster sequence is exposed (via a getter)
 369
 370   if (!fClustSeq) {
 371     AliError("[e] Run the jfinder first.");
 372     return;
 373   }
 374
 375   Double_t mean_area = 0;
 376   try {
 377     if(0 == remove) {
 378       fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
 379     }  else {
 380       std::vector<fastjet::PseudoJet> input_jets = sorted_by_pt(fClustSeq->inclusive_jets());
 381       input_jets.erase(input_jets.begin(), input_jets.begin() + remove);
 382       fClustSeq->get_median_rho_and_sigma(input_jets, *fRange, fUseArea4Vector, median, sigma, mean_area);
 383       input_jets.clear();
 384     }
 385   } catch (fj::Error) {
 386     AliError(" [w] FJ Exception caught.");
 387     median = -1.;
 388     sigma = -1;
 389   }
 390 }
 391
 392 //_________________________________________________________________________________________________
 393 Int_t AliFJWrapper::Run()
 394 {
 395   // Run the actual jet finder.
 396
 397   if (fAreaType == fj::voronoi_area) {
 398     // Rfact - check dependence - default is 1.
 399     // NOTE: hardcoded variable!
 400     fVorAreaSpec = new fj::VoronoiAreaSpec(1.);
 401     fAreaDef     = new fj::AreaDefinition(*fVorAreaSpec);
 402   } else {
 403     fGhostedAreaSpec = new fj::GhostedAreaSpec(fMaxRap,
 404                                                fNGhostRepeats,
 405                                                fGhostArea,
 406                                                fGridScatter,
 407                                                fKtScatter,
 408                                                fMeanGhostKt);
 409
 410     fAreaDef = new fj::AreaDefinition(*fGhostedAreaSpec, fAreaType);
 411   }
 412
 413   // this is acceptable by fastjet:
 414   fRange = new fj::RangeDefinition(fMaxRap - 0.95 * fR);
 415
 416   if (fAlgor == fj::plugin_algorithm) {
 417     if (fPluginAlgor == 0) {
 418       // SIS CONE ALGOR
 419       // NOTE: hardcoded split parameter
 420       Double_t overlap_threshold = 0.75; // NOTE: this actually splits a lot: thr/min(pt1,pt2)
 421       fPlugin = new fj::SISConePlugin(fR,
 422                                       overlap_threshold,
 423                                       0,    //search of stable cones - zero = until no more
 424                                       1.0); // this should be seed effectively for proto jets
 425       fJetDef = new fastjet::JetDefinition(fPlugin);
 426     } else if (fPluginAlgor == 1) {
 427       // CDF cone
 428       // NOTE: hardcoded split parameter
 429       Double_t overlap_threshold = 0.75; // NOTE: this actually splits a lot: thr/min(pt1,pt2)
 430       fPlugin = new fj::CDFMidPointPlugin(fR,
 431                                       overlap_threshold,
 432                                       1.0,    //search of stable cones - zero = until no more
 433                                       1.0); // this should be seed effectively for proto jets
 434       fJetDef = new fastjet::JetDefinition(fPlugin);
 435     } else {
 436       AliError("[e] Unrecognized plugin number!");
 437     }
 438   } else {
 439     fJetDef = new fj::JetDefinition(fAlgor, fR, fScheme, fStrategy);
 440   }
 441
 442   try {
 443     fClustSeq = new fj::ClusterSequenceArea(fInputVectors, *fJetDef, *fAreaDef);
 444   } catch (fj::Error) {
 445     AliError(" [w] FJ Exception caught.");
 446     return -1;
 447   }
 448
 449   // FJ3 :: Define an JetMedianBackgroundEstimator just in case it will be used
 450 #ifdef FASTJET_VERSION
 451   fBkrdEstimator     = new fj::JetMedianBackgroundEstimator(fj::SelectorAbsRapMax(fMaxRap));
 452 #endif
 453
 454   if (fLegacyMode) { SetLegacyFJ(); } // for FJ 2.x even if fLegacyMode is set, SetLegacyFJ is dummy
 455
 456   // inclusive jets:
 457   fInclusiveJets.clear();
 458   fInclusiveJets = fClustSeq->inclusive_jets(0.0);
 459
 460   return 0;
 461 }
 462
 463 //_________________________________________________________________________________________________
 464 void AliFJWrapper::SetLegacyFJ()
 465 {
 466   // This methods enable legacy behaviour (FastJet 2.x) when AliROOT is compiled with FastJet 3.x
 467 #ifdef FASTJET_VERSION
 468     std::cout << "WARNING! Setting FastJet in legacy mode" << std::endl;
 469     if (fGhostedAreaSpec) { fGhostedAreaSpec->set_fj2_placement(kTRUE); }
 470     if (fBkrdEstimator) {
 471       fBkrdEstimator->set_provide_fj2_sigma(kTRUE);
 472       fBkrdEstimator->set_use_area_4vector(kFALSE);
 473     }
 474 #endif
 475 }
 476
 477 //_________________________________________________________________________________________________
 478 void AliFJWrapper::SubtractBackground(Double_t median_pt)
 479 {
 480   // Subtract the background (specify the value - see below the meaning).
 481   // Negative argument means the bg will be determined with the current algorithm
 482   // this is the default behavior. Zero is allowed
 483   // Note: user may set the switch for area4vector based subtraction.
 484
 485   Double_t median    = 0;
 486   Double_t sigma     = 0;
 487   Double_t mean_area = 0;
 488
 489   // clear the subtracted jet pt's vector<double>
 490   fSubtractedJetsPt.clear();
 491
 492   // check what was specified (default is -1)
 493   if (median_pt < 0) {
 494     try {
 495       fClustSeq->get_median_rho_and_sigma(*fRange, fUseArea4Vector, median, sigma, mean_area);
 496     }
 497
 498     catch (fj::Error) {
 499       AliError(" [w] FJ Exception caught.");
 500       median = -9999.;
 501       sigma = -1;
 502       fMedUsedForBgSub = median;
 503       return;
 504     }
 505     fMedUsedForBgSub = median;
 506   } else {
 507     // we do not know the sigma in this case
 508     sigma = -1;
 509     if (0.0 == median_pt) {
 510       // AliWarning(Form(" [w] Median specified for bg subtraction is ZERO: %f", median_pt )) << endl;
 511       fMedUsedForBgSub = 0.;
 512     } else {
 513       fMedUsedForBgSub = median_pt;
 514     }
 515   }
 516
 517   // subtract:
 518   for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
 519     if ( fUseArea4Vector ) {
 520       // subtract the background using the area4vector
 521       fj::PseudoJet area4v = fClustSeq->area_4vector(fInclusiveJets[i]);
 522       fj::PseudoJet jet_sub = fInclusiveJets[i] - area4v * fMedUsedForBgSub;
 523       fSubtractedJetsPt.push_back(jet_sub.perp()); // here we put only the pt of the jet - note: this can be negative
 524     } else {
 525       // subtract the background using scalars
 526       // fj::PseudoJet jet_sub = fInclusiveJets[i] - area * fMedUsedForBgSub_;
 527       Double_t area = fClustSeq->area(fInclusiveJets[i]);
 528       // standard subtraction
 529       Double_t pt_sub = fInclusiveJets[i].perp() - fMedUsedForBgSub * area;
 530       fSubtractedJetsPt.push_back(pt_sub); // here we put only the pt of the jet - note: this can be negative
 531     }
 532   }
 533 }
 534
 535 //_________________________________________________________________________________________________
 536 Int_t AliFJWrapper::DoGenericSubtractionJetMass() {
 537   //Do generic subtraction for jet mass
 538 #ifdef FASTJET_CONTRIB
 539   if(fUseExternalBkg)   fGenSubtractor     = new fj::contrib::GenericSubtractor(fRho,fRhom);
 540   else                  fGenSubtractor     = new fj::contrib::GenericSubtractor(fBkrdEstimator);
 541
 542   // Define jet shape
 543   AliJetShapeMass shapeMass;
 544
 545   // clear the generic subtractor info vector
 546   fGenSubtractorInfoJetMass.clear();
 547   for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
 548     fj::contrib::GenericSubtractorInfo info;
 549     if(fInclusiveJets[i].perp()>0.)
 550       double subtracted_shape = (*fGenSubtractor)(shapeMass, fInclusiveJets[i], info);
 551     fGenSubtractorInfoJetMass.push_back(info);
 552   }
 553
 554 #endif
 555   return 0;
 556 }
 557
 558 //_________________________________________________________________________________________________
 559 Int_t AliFJWrapper::DoConstituentSubtraction() {
 560   //Do constituent subtraction
 561 #ifdef FASTJET_CONTRIB
 562   fj::contrib::ConstituentSubtractor *subtractor;
 563   if(fUseExternalBkg)
 564     subtractor     = new fj::contrib::ConstituentSubtractor(fRho,fRhom,kFALSE,kTRUE);
 565   else                 subtractor     = new fj::contrib::ConstituentSubtractor(fBkrdEstimator);
 566
 567   //clear constituent subtracted jets
 568   fConstituentSubtrJets.clear();
 569   for (unsigned i = 0; i < fInclusiveJets.size(); i++) {
 570     fj::PseudoJet subtracted_jet(0.,0.,0.,0.);
 571     if(fInclusiveJets[i].perp()>0.)
 572       subtracted_jet = (*subtractor)(fInclusiveJets[i]);
 573     fConstituentSubtrJets.push_back(subtracted_jet);
 574   }
 575   if(subtractor) delete subtractor;
 576
 577 #endif
 578   return 0;
 579 }
 580
 581 //_________________________________________________________________________________________________
 582 void AliFJWrapper::SetupAlgorithmfromOpt(const char *option)
 583 {
 584   // Setup algorithm from char.
 585
 586   std::string opt(option);
 587
 588   if (!opt.compare("kt"))                fAlgor    = fj::kt_algorithm;
 589   if (!opt.compare("antikt"))            fAlgor    = fj::antikt_algorithm;
 590   if (!opt.compare("cambridge"))         fAlgor    = fj::cambridge_algorithm;
 591   if (!opt.compare("genkt"))             fAlgor    = fj::genkt_algorithm;
 592   if (!opt.compare("cambridge_passive")) fAlgor    = fj::cambridge_for_passive_algorithm;
 593   if (!opt.compare("genkt_passive"))     fAlgor    = fj::genkt_for_passive_algorithm;
 594   if (!opt.compare("ee_kt"))             fAlgor    = fj::ee_kt_algorithm;
 595   if (!opt.compare("ee_genkt"))          fAlgor    = fj::ee_genkt_algorithm;
 596   if (!opt.compare("plugin"))            fAlgor    = fj::plugin_algorithm;
 597 }
 598
 599 //_________________________________________________________________________________________________
 600 void AliFJWrapper::SetupAreaTypefromOpt(const char *option)
 601 {
 602   // Setup area type from char.
 603
 604   std::string opt(option);
 605
 606   if (!opt.compare("active"))                      fAreaType = fj::active_area;
 607   if (!opt.compare("invalid"))                     fAreaType = fj::invalid_area;
 608   if (!opt.compare("active_area_explicit_ghosts")) fAreaType = fj::active_area_explicit_ghosts;
 609   if (!opt.compare("one_ghost_passive"))           fAreaType = fj::one_ghost_passive_area;
 610   if (!opt.compare("passive"))                     fAreaType = fj::passive_area;
 611   if (!opt.compare("voronoi"))                     fAreaType = fj::voronoi_area;
 612 }
 613
 614 //_________________________________________________________________________________________________
 615 void AliFJWrapper::SetupSchemefromOpt(const char *option)
 616 {
 617   //
 618   // setup scheme from char
 619   //
 620
 621   std::string opt(option);
 622
 623   if (!opt.compare("BIpt"))   fScheme   = fj::BIpt_scheme;
 624   if (!opt.compare("BIpt2"))  fScheme   = fj::BIpt2_scheme;
 625   if (!opt.compare("E"))      fScheme   = fj::E_scheme;
 626   if (!opt.compare("pt"))     fScheme   = fj::pt_scheme;
 627   if (!opt.compare("pt2"))    fScheme   = fj::pt2_scheme;
 628   if (!opt.compare("Et"))     fScheme   = fj::Et_scheme;
 629   if (!opt.compare("Et2"))    fScheme   = fj::Et2_scheme;
 630 }
 631
 632 //_________________________________________________________________________________________________
 633 void AliFJWrapper::SetupStrategyfromOpt(const char *option)
 634 {
 635   // Setup strategy from char.
 636
 637   std::string opt(option);
 638
 639   if (!opt.compare("Best"))            fStrategy = fj::Best;
 640   if (!opt.compare("N2MinHeapTiled"))  fStrategy = fj::N2MinHeapTiled;
 641   if (!opt.compare("N2Tiled"))         fStrategy = fj::N2Tiled;
 642   if (!opt.compare("N2PoorTiled"))     fStrategy = fj::N2PoorTiled;
 643   if (!opt.compare("N2Plain"))         fStrategy = fj::N2Plain;
 644   if (!opt.compare("N3Dumb"))          fStrategy = fj::N3Dumb;
 645   if (!opt.compare("NlnN"))            fStrategy = fj::NlnN;
 646   if (!opt.compare("NlnN3pi"))         fStrategy = fj::NlnN3pi;
 647   if (!opt.compare("NlnN4pi"))         fStrategy = fj::NlnN4pi;
 648   if (!opt.compare("NlnNCam4pi"))      fStrategy = fj::NlnNCam4pi;
 649   if (!opt.compare("NlnNCam2pi2R"))    fStrategy = fj::NlnNCam2pi2R;
 650   if (!opt.compare("NlnNCam"))         fStrategy = fj::NlnNCam;
 651   if (!opt.compare("plugin"))          fStrategy = fj::plugin_strategy;
 652 }
 653 #endif