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