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