ExhumeInterface/src/Weight.cpp

0001 //-*-c++-*-
0002 //-*-Weighting.cpp-*-
0003 //   Written by James Monk and Andrew Pilkington
0004 /////////////////////////////////////////////////////////////////////////////
0005 #include "GeneratorInterface/ExhumeInterface/interface/Weight.h"
0006
0007 /////////////////////////////////////////////////////////////////////////////
0008 void Exhume::Weight::WeightInit(const double &min, const double &max) {
0009   Max_ = max;
0010   FuncMap.insert(std::pair<double, double>(min, WeightFunc(min)));
0011
0012   std::map<double, double> UNLineShape;
0013
0014   UNLineShape.insert(std::pair<double, double>(0.0, min));
0015   UNLineShape.insert(std::pair<double, double>(1.0, max));
0016
0017   unsigned int jj10 = 1;
0018   double incr, y, x;
0019   double NewF, OldF, Oldx;
0020
0021   while (jj10 < NPoints) {
0022     jj10 = jj10 * 10;
0023     LineShape = UNLineShape;
0024     UNLineShape.clear();
0025     UNLineShape.insert(std::pair<double, double>(0.0, min));
0026     incr = (LineShape.rbegin()->first) / double(jj10 - 1.0);
0027     y = 0.0;
0028     x = GetValue(y);
0029     NewF = WeightFunc(x);
0030
0031     for (unsigned int kk = 1; kk < jj10; kk++) {
0032       y = y + incr;
0033       Oldx = x;
0034       x = GetValue(y);
0035       OldF = NewF;
0036       NewF = WeightFunc(x);
0037       FuncMap.insert(std::pair<double, double>(x, NewF));
0038       UNLineShape.insert(std::pair<double, double>(UNLineShape.rbegin()->first + 0.5 * (NewF + OldF) * (x - Oldx), x));
0039     }
0040   }
0041
0042   LineShape.clear();
0043   UNLineShape.clear();
0044   UNLineShape.insert(std::pair<double, double>(0.0, min));
0045
0046   std::map<double, double>::iterator jj = FuncMap.begin();
0047   jj++;
0048   for (std::map<double, double>::iterator ii = FuncMap.begin(); jj != FuncMap.end(); ii++) {
0049     UNLineShape.insert(std::pair<double, double>(
0050         UNLineShape.rbegin()->first + 0.5 * (ii->second + jj->second) * (jj->first - ii->first), jj->first));
0051     jj++;
0052   }
0053
0054   TotalIntegral = UNLineShape.rbegin()->first;
0055   double C = 1.0 / TotalIntegral;
0056
0057   for (jj = UNLineShape.begin(); jj != UNLineShape.end(); jj++) {
0058     LineShape.insert(std::pair<double, double>(C * (jj->first), jj->second));
0059   }
0060
0061   UNLineShape.clear();
0062
0063   std::cout << std::endl
0064             << "  " << LineShape.size() << " points calculated mapping {" << LineShape.begin()->first << ","
0065             << LineShape.rbegin()->first << "} to {" << (LineShape.begin())->second << ","
0066             << (LineShape.rbegin())->second << "}" << std::endl
0067             << std::endl;
0068
0069   std::cout << "  Event weighting initialised" << std::endl;
0070
0071   std::cout << std::endl << " ........................................................................." << std::endl;
0072
0073   return;
0074 }
0075 /////////////////////////////////////////////////////////////////////////////
0076 void Exhume::Weight::AddPoint(const double &xx_, const double &f_) {
0077   std::map<double, double>::iterator high_, low_;
0078   high_ = FuncMap.upper_bound(xx_);
0079   low_ = high_;
0080   low_--;
0081
0082   double OldSegment = 0.5 * (high_->second + low_->second) * (high_->first - low_->first);
0083
0084   double NewSegment = 0.5 * ((low_->second + f_) * (xx_ - low_->first) + (high_->second + f_) * (high_->first - xx_));
0085
0086   std::cout << "   Adding point to weight function map" << std::endl;
0087   std::cout << "   Updating TotalIntegral from " << TotalIntegral;
0088
0089   TotalIntegral = TotalIntegral - OldSegment + NewSegment;
0090
0091   std::cout << " to " << TotalIntegral << std::endl;
0092
0093   FuncMap.insert(std::pair<double, double>(xx_, f_));
0094
0095   return;
0096 }
0097 /////////////////////////////////////////////////////////////////////////////