PWGJE/EMCALJetTasks/AliFJWrapper.h

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