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