rivet is hosted by Hepforge, IPPP Durham
Rivet 4.0.0
CentralityBinner.hh
1// -*- C++ -*-
2#ifndef RIVET_CENTRALITYBINNER_HH
3#define RIVET_CENTRALITYBINNER_HH
4#include <tuple>
5#include "Rivet/Config/RivetCommon.hh"
6#include "Rivet/Tools/RivetYODA.hh"
7#include "Rivet/Projections/HepMCHeavyIon.hh"
8
9namespace Rivet {
10
11
35 public:
36
38 CentralityEstimator(): _estimate(-1.0) {
39 setName("CentralityEstimator");
40 declare(HepMCHeavyIon(), "HepMC");
41 }
42
45
46 protected:
47
49 void project(const Event& e) {
50 _estimate = -1.0;
51 double imp = apply<HepMCHeavyIon>(e, "HepMC").impact_parameter();
52 if ( imp < 0.0 ) return;
53 _estimate = imp > 0.0? 1.0/imp: numeric_limits<double>::max();
54 }
55
57 CmpState compare(const Projection& p) const {
58 return mkNamedPCmp(p, "HepMC");
59 }
60
61
62 public:
63
65 double estimate() const { return _estimate; }
66
67
68 protected:
69
71 double _estimate;
72
73 };
74
75
79 template <typename T>
80 struct CentralityBinTraits {
81
83 static T clone(const T & t) {
84 return T(t->newclone());
85 }
86
88 static void add(T & t, const T & o) {
89 *t += *o;
90 }
91
93 static void scale(T & t, double f) {
94 t->scaleW(f);
95 }
96
99 static void normalize(T & t, double sumw) {
100 if ( t->sumW() > 0.0 ) t->normalize(t->sumW()/sumw);
101 }
102
104 static string path(T t) {
105 return t->path();
106 }
107
108 };
109
113 struct MergeDistance {
114
124 static double dist(double cestLo, double cestHi, double weight,
125 double clo, double chi, double, double) {
126 return (cestHi - cestLo)*weight/(cestHi*(chi - clo));
127 }
128 };
129
130
138 template <typename T = Histo1DPtr, typename MDist = MergeDistance>
140 public:
141
146 CentralityBinner(int maxbins = 200, double wlim = 0.02)
147 : _currentCEst(-1.0), _maxBins(maxbins), _warnlimit(wlim), _weightsum(0.0) {
148 _percentiles.insert(0.0);
149 _percentiles.insert(1.0);
150 }
151
154 void setProjection(const CentralityEstimator & p, string pname) {
155 declare(p, pname);
156 _estimator = pname;
157 }
158
160 virtual std::string name() const {
161 return "Rivet::CentralityBinner";
162 }
163
168
172
175
176 void add(T t, double cmin, double cmax,
177 double cestmin = -1.0, double cestmax = -1.0 ) {
178 _percentiles.insert(max(1.0 - cmax/100.0, 0.0));
179 _percentiles.insert(min(1.0 - cmin/100.0, 1.0));
180 if ( _unfilled.empty() && _ready.empty() )
181 _devnull = CentralityBinTraits<T>::clone(t);
182 if ( cestmin < 0.0 )
183 _unfilled.push_back(Bin(t, 1.0 - cmax/100.0, 1.0 - cmin/100.0));
184 else
185 _ready[t] = Bin(t, 1.0 - cmax/100.0, 1.0 - cmin/100.0, cestmin, cestmax);
186 }
187
196 T select(const Event & event, double weight = 1.0) {
197 return select(applyProjection<CentralityEstimator>
198 (event, _estimator).estimate(), weight);
199 }
200
207 T select(double cest, double weight = 1.0);
208
215 void finalize();
216
220 for ( auto & b : _ready ) b.second.normalizePerEvent();
221 }
222
225 map<double,double> edges() const {
226 map<double,double> ret;
227 for ( auto & b : _ready ) {
228 ret[1.0 - b.second._centLo] = b.second._cestLo;
229 ret[1.0 - b.second._centHi] = b.second._cestHi;
230 }
231 return ret;
232 }
233
235 const T & current() const {
236 return _currenT;
237 }
238
241 double estimator() const {
242 return _currentCEst;
243 }
244
245 vector<T> allObjects() {
246 vector<T> ret;
247 for ( auto & fb : _flexiBins ) ret.push_back(fb._t);
248 if ( !ret.empty() ) return ret;
249 for ( auto b : _ready ) ret.push_back(b.second._t);
250 return ret;
251 }
252
253 private:
254
256 struct FlexiBin {
257
260 FlexiBin(T & t, double cest = 0.0, double weight = 0.0)
261 : _t(t), _cestLo(cest), _cestHi(cest), _weightsum(weight), _n(1), _m(0) {}
262
264 FlexiBin(double cest)
265 : _cestLo(cest), _cestHi(cest), _weightsum(0.0), _n(0), _m(0) {}
266
268 void merge(const FlexiBin & fb) {
269 _cestLo = min(_cestLo, fb._cestLo);
270 _cestHi = max(_cestHi, fb._cestHi);
271 _weightsum += fb._weightsum;
272 CentralityBinTraits<T>::add(_t, fb._t);
273 _n += fb._n;
274 _m += fb._m + 1;
275 }
276
278 bool operator< (const FlexiBin & fb) const {
279 return _cestLo < fb._cestLo;
280 }
281
284 bool inRange(double cest) const {
285 return cest == _cestLo || ( _cestLo < cest && cest < _cestHi );
286 }
287
289 T _t;
290
293 double _cestLo, _cestHi;
294
297 mutable double _weightsum;
298
300 mutable int _n;
301
303 mutable int _m;
304
305 };
306
307 struct Bin {
308
310 Bin()
311 : _centLo(-1.0), _centHi(-1.0), _cestLo(-1.0), _cestHi(-1.0),
312 _weightsum(0.0), _underflow(0.0), _overflow(0.0),
313 _ambiguous(0), _ambweight(0.0) {}
314
319 Bin(T t, double centLo, double centHi,
320 double cestLo = -1.0, double cestHi = -1.0)
321 : _t(t), _centLo(centLo), _centHi(centHi),
322 _cestLo(cestLo), _cestHi(cestHi),
323 _weightsum(0.0), _underflow(0.0), _overflow(0.0),
324 _ambiguous(0.0), _ambweight(0.0) {}
325
328 bool inRange(double cest) const {
329 return _cestLo >= 0 && _cestLo <= cest &&
330 ( _cestHi < 0.0 || cest <= _cestHi );
331 }
332
334 void normalizePerEvent() {
335 CentralityBinTraits<T>::normalize(_t, _weightsum);
336 }
337
339 T _t;
340
342 double _centLo, _centHi;
343
345 double _cestLo, _cestHi;
346
348 double _weightsum;
349
352 double _underflow;
353
356 double _overflow;
357
359 double _ambiguous;
360
362 double _ambweight;
363
364 };
365
366 protected:
367
369 typedef set<FlexiBin> FlexiBinSet;
370
373 typename FlexiBinSet::iterator _findBin(double cest) {
374 if ( _flexiBins.empty() ) return _flexiBins.end();
375 auto it = _flexiBins.lower_bound(FlexiBin(cest));
376 if ( it->_cestLo == cest ) return it;
377 if ( it != _flexiBins.begin() ) --it;
378 if ( it->_cestLo < cest && cest < it->_cestHi ) return it;
379 return _flexiBins.end();
380 }
381
383 string _estimator;
384
387 T _currenT;
388
390 double _currentCEst;
391
394 int _maxBins;
395
398 double _warnlimit;
399
402 vector<Bin> _unfilled;
403
405 FlexiBinSet _flexiBins;
406
408 double _weightsum;
409
411 set<double> _percentiles;
412
414 map<T, Bin> _ready;
415
418 T _devnull;
419
420 public:
421
423 void debug();
424 void fulldebug();
425
426 };
427
428
430 template <>
431 struct CentralityBinTraits<Profile1DPtr> {
432
433 typedef Profile1DPtr T;
434
436 static T clone(const T & t) {
437 return Profile1DPtr(t->newclone());
438 }
439
441 static void add(T & t, const T & o) {
442 *t += *o;
443 }
444
446 static void scale(T & t, double f) {
447 t->scaleW(f);
448 }
449
450 static void normalize(T & t, double sumw) {}
451
453 static string path(T t) {
454 return t->path();
455 }
456
457 };
458
459
461 template <>
462 struct CentralityBinTraits<Profile2DPtr> {
463
464 typedef Profile2DPtr T;
465
467 static T clone(const T & t) {
468 return Profile2DPtr(t->newclone());
469 }
470
472 static void add(T & t, const T & o) {
473 *t += *o;
474 }
475
477 static void scale(T & t, double f) {
478 t->scaleW(f);
479 }
480
481 static void normalize(T & t, double sumw) {}
482
484 static string path(T t) {
485 return t->path();
486 }
487
488 };
489
490 template <typename T>
491 struct CentralityBinTraits< vector<T> > {
492
494 static vector<T> clone(const vector<T> & tv) {
495 vector<T> rtv;
496 for ( auto t : tv ) rtv.push_back(CentralityBinTraits<T>::clone(t));
497 return rtv;
498 }
499
501 static void add(vector<T> & tv, const vector<T> & ov) {
502 for ( int i = 0, N = tv.size(); i < N; ++i )
503 CentralityBinTraits::add(tv[i], ov[i]);
504 }
505
507 static void scale(vector<T> & tv, double f) {
508 for ( auto t : tv ) CentralityBinTraits<T>::scale(t, f);
509 }
510
511 static void normalize(vector<T> & tv, double sumw) {
512 for ( auto t : tv ) CentralityBinTraits<T>::normalize(t, sumw);
513 }
514
516 static string path(const vector<T> & tv) {
517 string ret = "(vector:";
518 for ( auto t : tv ) {
519 ret += " ";
520 ret += CentralityBinTraits<T>::path(t);
521 }
522 ret += ")";
523 return ret;
524 }
525
526 };
527
528 template <size_t I, typename... Types>
529 struct TupleCentralityBinTraitsHelper {
530
531 typedef tuple<Types...> Tuple;
532 typedef typename tuple_element<I-1,Tuple>::type T;
533
534 static void clone(Tuple & ret, const Tuple & tup) {
535 get<I-1>(ret) = CentralityBinTraits<T>::clone(get<I-1>(tup));
536 TupleCentralityBinTraitsHelper<I-1,Types...>::clone(ret, tup);
537 }
538
539 static void add(Tuple & tup, const Tuple & otup) {
540 CentralityBinTraits<T>::add(get<I-1>(tup),get<I-1>(otup));
541 TupleCentralityBinTraitsHelper<I-1,Types...>::add(tup, otup);
542 }
543
544 static void scale(Tuple & tup, double f) {
545 CentralityBinTraits<T>::scale(get<I-1>(tup), f);
546 TupleCentralityBinTraitsHelper<I-1,Types...>::scale(tup, f);
547 }
548
549 static void normalize(Tuple & tup, double sumw) {
550 CentralityBinTraits<T>::normalize(get<I-1>(tup), sumw);
551 TupleCentralityBinTraitsHelper<I-1,Types...>::normalize(tup, sumw);
552 }
553
554 static string path(const Tuple & tup) {
555 return " " + CentralityBinTraits<T>::path(get<I-1>(tup))
556 + TupleCentralityBinTraitsHelper<I-1,Types...>::path(tup);
557 }
558 };
559
560 template <typename... Types>
561 struct TupleCentralityBinTraitsHelper<0,Types...> {
562
563 typedef tuple<Types...> Tuple;
564
565 static void clone(Tuple &, const Tuple &) {}
566 static void add(Tuple & tup, const Tuple & otup) {}
567 static void scale(Tuple & tup, double f) {}
568 static void normalize(Tuple & tup, double sumw) {}
569 static string path(const Tuple & tup) {return "";}
570
571 };
572
573 template <typename... Types>
574 struct CentralityBinTraits< tuple<Types...> > {
575
576 typedef tuple<Types...> Tuple;
577 static const size_t N = tuple_size<Tuple>::value;
578
580 static Tuple clone(const Tuple & tup) {
581 Tuple ret;
582 TupleCentralityBinTraitsHelper<N,Types...>::clone(ret, tup);
583 return ret;
584 }
585
587 static void add(Tuple & tup, const Tuple & otup) {
588 TupleCentralityBinTraitsHelper<N,Types...>::add(tup, otup);
589 }
590
592 static void scale(Tuple & tup, double f) {
593 TupleCentralityBinTraitsHelper<N,Types...>::scale(tup, f);
594 }
595
596 static void normalize(Tuple & tup, double sumw) {
597 TupleCentralityBinTraitsHelper<N,Types...>::normalize(tup, sumw);
598 }
599
601 static string path(const Tuple & tup) {
602 string ret = "(tuple:";
603 ret += TupleCentralityBinTraitsHelper<N,Types...>::path(tup);
604 ret += ")";
605 return ret;
606 }
607
608 };
609
610 template <typename T, typename MDist>
611 T CentralityBinner<T,MDist>::select(double cest, double weight) {
612 _currenT = _devnull;
613 _currentCEst = cest;
614 _weightsum += weight;
615
616 // If estimator is negative, something has gone wrong.
617 if ( _currentCEst < 0.0 ) return _currenT;
618
619 // If we already have finalized the limits on the centrality
620 // estimator, we just add the weights to their bins and return the
621 // corresponding AnalysisObject.
622 if ( _unfilled.empty() ) {
623 for ( auto & b : _ready ) if ( b.second.inRange(_currentCEst) ) {
624 b.second._weightsum += weight;
625 return b.second._t;
626 }
627 return _currenT;
628 }
629
630 auto it = _findBin(cest);
631 if ( it == _flexiBins.end() ) {
632 _currenT = CentralityBinTraits<T>::clone(_unfilled.begin()->_t);
633 it = _flexiBins.insert(FlexiBin(_currenT, _currentCEst, weight)).first;
634 } else {
635 it->_weightsum += weight;
636 ++(it->_n);
637 _currenT = it->_t;
638 }
639
640 if ( (int)_flexiBins.size() <= _maxBins ) return _currenT;
641
642
643 set<double>::iterator citn = _percentiles.begin();
644 set<double>::iterator cit0 = citn++;
645 auto selectit = _flexiBins.end();
646 double mindist = -1.0;
647 double acc = 0.0;
648 auto next = _flexiBins.begin();
649 auto prev = next++;
650 for ( ; next != _flexiBins.end(); prev = next++ ) {
651 acc += prev->_weightsum/_weightsum;
652 if ( acc > *citn ) {
653 cit0 = citn++;
654 continue;
655 }
656 if ( acc + next->_weightsum/_weightsum > *citn ) continue;
657 double dist = MDist::dist(prev->_cestLo, next->_cestHi,
658 next->_weightsum + prev->_weightsum,
659 *cit0, *citn, next->_n + prev->_n,
660 next->_m + prev->_m);
661 if ( mindist < 0.0 || dist < mindist ) {
662 selectit = prev;
663 mindist = dist;
664 }
665 }
666
667 if ( selectit == _flexiBins.end() ) return _currenT;
668 auto mergeit = selectit++;
669 FlexiBin merged = *mergeit;
670 merged.merge(*selectit);
671 if ( merged.inRange(cest) || selectit->inRange(cest) )
672 _currenT = merged._t;
673 _flexiBins.erase(mergeit);
674 _flexiBins.erase(selectit);
675 _flexiBins.insert(merged);
676
677 return _currenT;
678
679 }
680
681
682 template <typename T, typename MDist>
684
685 // Take the contents of the dynamical binning and fill the original
686 // AnalysisObjects.
687
688 double clo = 0.0;
689 for ( const FlexiBin & fb : _flexiBins ) {
690 double chi = min(clo + fb._weightsum/_weightsum, 1.0);
691 for ( Bin & bin : _unfilled ) {
692 double olo = bin._centLo;
693 double ohi = bin._centHi;
694 if ( clo > ohi || chi <= olo ) continue;
695 // If we only have partial overlap we need to scale
696 double lo = max(olo, clo);
697 double hi = min(ohi, chi);
698 T t = CentralityBinTraits<T>::clone(fb._t);
699 double frac = (hi - lo)/(chi - clo);
700 CentralityBinTraits<T>::scale(t, frac);
701 CentralityBinTraits<T>::add(bin._t, t);
702 bin._weightsum += fb._weightsum*frac;
703 if ( clo <= olo ) bin._cestLo = fb._cestLo +
704 (fb._cestHi - fb._cestLo)*(olo - clo)/(chi - clo);
705 if ( clo < olo ) {
706 bin._underflow = clo;
707 bin._ambiguous += fb._n*frac;
708 bin._ambweight += fb._weightsum*frac*(1.0 - frac);
709 }
710 if ( chi > ohi ) {
711 bin._cestHi =
712 fb._cestLo + (fb._cestHi - fb._cestLo)*(ohi - clo)/(chi - clo);
713 bin._overflow = chi;
714 bin._ambiguous += fb._n*frac;
715 bin._ambweight += fb._weightsum*frac*(1.0 - frac);
716 }
717 }
718 clo = chi;
719 }
720 _flexiBins.clear();
721 for ( Bin & bin : _unfilled ) {
722 if ( bin._overflow == 0.0 ) bin._overflow = 1.0;
723 _ready[bin._t] = bin;
724 if ( bin._ambweight/bin._weightsum >_warnlimit )
725 MSG_WARNING("Analysis object \"" << CentralityBinTraits<T>::path(bin._t)
726 << "\", contains events with centralities between "
727 << bin._underflow*100.0
728 << " and " << bin._overflow*100.0 << "% ("
729 << int(bin._ambiguous + 0.5)
730 << " ambiguous events with effectively "
731 << 100.0*bin._ambweight/bin._weightsum
732 << "% of the weights)."
733 << "Consider increasing the number of bins.");
734
735 }
736 _unfilled.clear();
737
738 }
739
740 template <typename T, typename MDist>
742 cerr << endl;
743 double acc = 0.0;
744 set<double>::iterator citn = _percentiles.begin();
745 set<double>::iterator cit0 = citn++;
746 int i = 0;
747 for ( auto it = _flexiBins.begin(); it != _flexiBins.end(); ) {
748 ++i;
749 auto curr = it++;
750 double w = curr->_weightsum/_weightsum;
751 acc += w;
752 if ( curr == _flexiBins.begin() || it == _flexiBins.end() || acc > *citn )
753 cerr << "*";
754 else
755 cerr << " ";
756 if ( acc > *citn ) cit0 = citn++;
757 cerr << setw(6) << i
758 << setw(12) << acc - w
759 << setw(12) << acc
760 << setw(8) << curr->_n
761 << setw(8) << curr->_m
762 << setw(12) << curr->_cestLo
763 << setw(12) << curr->_cestHi << endl;
764 }
765 cerr << "Number of sampler bins: " << _flexiBins.size() << endl;
766 }
767
768 template <typename T, typename MDist>
770 cerr << endl;
771 double acc = 0.0;
772 int i = 0;
773 set<double>::iterator citn = _percentiles.begin();
774 set<double>::iterator cit0 = citn++;
775 for ( auto it = _flexiBins.begin(); it != _flexiBins.end(); ) {
776 auto curr = it++;
777 ++i;
778 double w = curr->_weightsum/_weightsum;
779 acc += w;
780 if ( curr == _flexiBins.begin() || it == _flexiBins.end() || acc > *citn ) {
781 if ( acc > *citn ) cit0 = citn++;
782 cerr << setw(6) << i
783 << setw(12) << acc - w
784 << setw(12) << acc
785 << setw(8) << curr->_n
786 << setw(8) << curr->_m
787 << setw(12) << curr->_cestLo
788 << setw(12) << curr->_cestHi << endl;
789
790 }
791 }
792 cerr << "Number of sampler bins: " << _flexiBins.size() << endl;
793 }
794
798
799 public:
800
805
808
809 protected:
810
812 void project(const Event& e) {
813 _estimate = apply<HepMCHeavyIon>(e, "HI").centrality();
814 }
815
817 CmpState compare(const Projection& p) const {
818 return mkNamedPCmp(p, "GeneratedCentrality");
819 }
820
821 };
822
823
824
825}
826
827#endif
Definition CentralityBinner.hh:139
double estimator() const
Definition CentralityBinner.hh:241
void normalizePerEvent()
Definition CentralityBinner.hh:219
map< double, double > edges() const
Definition CentralityBinner.hh:225
T select(const Event &event, double weight=1.0)
Definition CentralityBinner.hh:196
const T & current() const
Return the current AnalysisObject from the latest call to select().
Definition CentralityBinner.hh:235
void add(T t, double cmin, double cmax, double cestmin=-1.0, double cestmax=-1.0)
Definition CentralityBinner.hh:176
void finalize()
Definition CentralityBinner.hh:683
set< FlexiBin > FlexiBinSet
Convenient typedefs.
Definition CentralityBinner.hh:369
void setProjection(const CentralityEstimator &p, string pname)
Definition CentralityBinner.hh:154
CentralityBinner(int maxbins=200, double wlim=0.02)
Definition CentralityBinner.hh:146
void debug()
Print out the _flexiBins to cerr.
Definition CentralityBinner.hh:769
virtual std::string name() const
Return the class name.
Definition CentralityBinner.hh:160
Base class for projections profile observable value vs the collision centrality.
Definition CentralityBinner.hh:34
void project(const Event &e)
Perform the projection on the Event.
Definition CentralityBinner.hh:49
CentralityEstimator()
Constructor.
Definition CentralityBinner.hh:38
CmpState compare(const Projection &p) const
Compare projections.
Definition CentralityBinner.hh:57
RIVET_DEFAULT_PROJ_CLONE(CentralityEstimator)
Clone on the heap.
double estimate() const
The value of the centrality estimate.
Definition CentralityBinner.hh:65
Representation of a HepMC event, and enabler of Projection caching.
Definition Event.hh:22
Definition CentralityBinner.hh:797
GeneratedCentrality()
Constructor.
Definition CentralityBinner.hh:802
CmpState compare(const Projection &p) const
Compare projections.
Definition CentralityBinner.hh:817
void project(const Event &e)
Perform the projection on the Event.
Definition CentralityBinner.hh:812
RIVET_DEFAULT_PROJ_CLONE(GeneratedCentrality)
Clone on the heap.
Definition HepMCHeavyIon.hh:12
Common base class for Projection and Analysis, used for internal polymorphism.
Definition ProjectionApplier.hh:22
const PROJ & declare(const PROJ &proj, const std::string &name) const
Register a contained projection (user-facing version)
Definition ProjectionApplier.hh:175
Base class for all Rivet projections.
Definition Projection.hh:29
void setName(const std::string &name)
Used by derived classes to set their name.
Definition Projection.hh:148
Cmp< Projection > mkNamedPCmp(const Projection &otherparent, const std::string &pname) const
#define MSG_WARNING(x)
Warning messages for non-fatal bad things, using MSG_LVL.
Definition Logging.hh:187
Definition MC_CENT_PPB_Projections.hh:10
std::enable_if< std::is_arithmetic< N1 >::value &&std::is_arithmetic< N2 >::value, typenamestd::common_type< N1, N2 >::type >::type max(N1 a, N2 b)
Get the maximum of two numbers.
Definition MathUtils.hh:111
std::enable_if< std::is_arithmetic< N1 >::value &&std::is_arithmetic< N2 >::value, typenamestd::common_type< N1, N2 >::type >::type min(N1 a, N2 b)
Get the minimum of two numbers.
Definition MathUtils.hh:102
std::enable_if< std::is_arithmetic< N1 >::value &&std::is_arithmetic< N2 >::value &&std::is_arithmetic< N3 >::value, bool >::type inRange(N1 value, N2 low, N3 high, RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN)
Determine if value is in the range low to high, for floating point numbers.
Definition MathUtils.hh:133