c9baf2eea00910bb55227711ac5aed1c6032fe3b
[u/mrichter/AliRoot.git] / 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
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