|
|
Rivet 4.0.0
|
2#ifndef RIVET_CENTRALITYBINNER_HH
3#define RIVET_CENTRALITYBINNER_HH
5#include "Rivet/Config/RivetCommon.hh"
6#include "Rivet/Tools/RivetYODA.hh"
7#include "Rivet/Projections/HepMCHeavyIon.hh"
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();
80 struct CentralityBinTraits {
83 static T clone( const T & t) {
84 return T(t->newclone());
88 static void add(T & t, const T & o) {
93 static void scale(T & t, double f) {
99 static void normalize(T & t, double sumw) {
100 if ( t->sumW() > 0.0 ) t->normalize(t->sumW()/sumw);
104 static string path(T t) {
113 struct MergeDistance {
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));
138 template < typename T = Histo1DPtr, typename MDist = MergeDistance>
147 : _currentCEst(-1.0), _maxBins(maxbins), _warnlimit(wlim), _weightsum(0.0) {
148 _percentiles.insert(0.0);
149 _percentiles.insert(1.0);
160 virtual std::string name() const {
161 return "Rivet::CentralityBinner";
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);
183 _unfilled.push_back(Bin(t, 1.0 - cmax/100.0, 1.0 - cmin/100.0));
185 _ready[t] = Bin(t, 1.0 - cmax/100.0, 1.0 - cmin/100.0, cestmin, cestmax);
197 return select(applyProjection<CentralityEstimator>
198 (event, _estimator).estimate(), weight);
207 T select( double cest, double weight = 1.0);
220 for ( auto & b : _ready ) b.second.normalizePerEvent();
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;
245 vector<T> allObjects() {
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);
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) {}
264 FlexiBin( double cest)
265 : _cestLo(cest), _cestHi(cest), _weightsum(0.0), _n(0), _m(0) {}
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);
278 bool operator< ( const FlexiBin & fb) const {
279 return _cestLo < fb._cestLo;
284 bool inRange( double cest) const {
285 return cest == _cestLo || ( _cestLo < cest && cest < _cestHi );
293 double _cestLo, _cestHi;
297 mutable double _weightsum;
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) {}
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) {}
328 bool inRange( double cest) const {
329 return _cestLo >= 0 && _cestLo <= cest &&
330 ( _cestHi < 0.0 || cest <= _cestHi );
334 void normalizePerEvent() {
335 CentralityBinTraits<T>::normalize(_t, _weightsum);
342 double _centLo, _centHi;
345 double _cestLo, _cestHi;
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();
402 vector<Bin> _unfilled;
411 set<double> _percentiles;
431 struct CentralityBinTraits<Profile1DPtr> {
433 typedef Profile1DPtr T;
436 static T clone( const T & t) {
437 return Profile1DPtr(t->newclone());
441 static void add(T & t, const T & o) {
446 static void scale(T & t, double f) {
450 static void normalize(T & t, double sumw) {}
453 static string path(T t) {
462 struct CentralityBinTraits<Profile2DPtr> {
464 typedef Profile2DPtr T;
467 static T clone( const T & t) {
468 return Profile2DPtr(t->newclone());
472 static void add(T & t, const T & o) {
477 static void scale(T & t, double f) {
481 static void normalize(T & t, double sumw) {}
484 static string path(T t) {
490 template < typename T>
491 struct CentralityBinTraits< vector<T> > {
494 static vector<T> clone( const vector<T> & tv) {
496 for ( auto t : tv ) rtv.push_back(CentralityBinTraits<T>::clone(t));
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]);
507 static void scale(vector<T> & tv, double f) {
508 for ( auto t : tv ) CentralityBinTraits<T>::scale(t, f);
511 static void normalize(vector<T> & tv, double sumw) {
512 for ( auto t : tv ) CentralityBinTraits<T>::normalize(t, sumw);
516 static string path( const vector<T> & tv) {
517 string ret = "(vector:";
518 for ( auto t : tv ) {
520 ret += CentralityBinTraits<T>::path(t);
528 template < size_t I, typename... Types>
529 struct TupleCentralityBinTraitsHelper {
531 typedef tuple<Types...> Tuple;
532 typedef typename tuple_element<I-1,Tuple>::type T;
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);
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);
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);
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);
554 static string path( const Tuple & tup) {
555 return " " + CentralityBinTraits<T>::path(get<I-1>(tup))
556 + TupleCentralityBinTraitsHelper<I-1,Types...>::path(tup);
560 template < typename... Types>
561 struct TupleCentralityBinTraitsHelper<0,Types...> {
563 typedef tuple<Types...> Tuple;
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 "";}
573 template < typename... Types>
574 struct CentralityBinTraits< tuple<Types...> > {
576 typedef tuple<Types...> Tuple;
577 static const size_t N = tuple_size<Tuple>::value;
580 static Tuple clone( const Tuple & tup) {
582 TupleCentralityBinTraitsHelper<N,Types...>::clone(ret, tup);
587 static void add(Tuple & tup, const Tuple & otup) {
588 TupleCentralityBinTraitsHelper<N,Types...>::add(tup, otup);
592 static void scale(Tuple & tup, double f) {
593 TupleCentralityBinTraitsHelper<N,Types...>::scale(tup, f);
596 static void normalize(Tuple & tup, double sumw) {
597 TupleCentralityBinTraitsHelper<N,Types...>::normalize(tup, sumw);
601 static string path( const Tuple & tup) {
602 string ret = "(tuple:";
603 ret += TupleCentralityBinTraitsHelper<N,Types...>::path(tup);
610 template < typename T, typename MDist>
614 _weightsum += weight;
617 if ( _currentCEst < 0.0 ) return _currenT;
622 if ( _unfilled.empty() ) {
623 for ( auto & b : _ready ) if ( b.second.inRange(_currentCEst) ) {
624 b.second._weightsum += weight;
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;
635 it->_weightsum += weight;
640 if ( ( int)_flexiBins.size() <= _maxBins ) return _currenT;
643 set<double>::iterator citn = _percentiles.begin();
644 set<double>::iterator cit0 = citn++;
645 auto selectit = _flexiBins.end();
646 double mindist = -1.0;
648 auto next = _flexiBins.begin();
650 for ( ; next != _flexiBins.end(); prev = next++ ) {
651 acc += prev->_weightsum/_weightsum;
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 ) {
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);
682 template < typename T, typename MDist>
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;
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);
706 bin._underflow = clo;
707 bin._ambiguous += fb._n*frac;
708 bin._ambweight += fb._weightsum*frac*(1.0 - frac);
712 fb._cestLo + (fb._cestHi - fb._cestLo)*(ohi - clo)/(chi - clo);
714 bin._ambiguous += fb._n*frac;
715 bin._ambweight += fb._weightsum*frac*(1.0 - frac);
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.");
740 template < typename T, typename MDist>
744 set<double>::iterator citn = _percentiles.begin();
745 set<double>::iterator cit0 = citn++;
747 for ( auto it = _flexiBins.begin(); it != _flexiBins.end(); ) {
750 double w = curr->_weightsum/_weightsum;
752 if ( curr == _flexiBins.begin() || it == _flexiBins.end() || acc > *citn )
756 if ( acc > *citn ) cit0 = citn++;
758 << setw(12) << acc - w
760 << setw(8) << curr->_n
761 << setw(8) << curr->_m
762 << setw(12) << curr->_cestLo
763 << setw(12) << curr->_cestHi << endl;
765 cerr << "Number of sampler bins: " << _flexiBins.size() << endl;
768 template < typename T, typename MDist>
773 set<double>::iterator citn = _percentiles.begin();
774 set<double>::iterator cit0 = citn++;
775 for ( auto it = _flexiBins.begin(); it != _flexiBins.end(); ) {
778 double w = curr->_weightsum/_weightsum;
780 if ( curr == _flexiBins.begin() || it == _flexiBins.end() || acc > *citn ) {
781 if ( acc > *citn ) cit0 = citn++;
783 << setw(12) << acc - w
785 << setw(8) << curr->_n
786 << setw(8) << curr->_m
787 << setw(12) << curr->_cestLo
788 << setw(12) << curr->_cestHi << endl;
792 cerr << "Number of sampler bins: " << _flexiBins.size() << endl;
813 _estimate = apply<HepMCHeavyIon>(e, "HI").centrality();
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
|