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