00001
00002 #include "Rivet/Rivet.hh"
00003 #include "Rivet/RivetAIDA.hh"
00004 #include "Rivet/AnalysisHandler.hh"
00005 #include "Rivet/AnalysisInfo.hh"
00006 #include "Rivet/Analysis.hh"
00007 #include "Rivet/Tools/Logging.hh"
00008 #include "LWH/AIManagedObject.h"
00009 using namespace AIDA;
00010
00011 namespace Rivet {
00012
00013
00014 namespace {
00015 string makeAxisCode(const size_t datasetId, const size_t xAxisId, const size_t yAxisId) {
00016 stringstream axisCode;
00017 axisCode << "d";
00018 if (datasetId < 10) axisCode << 0;
00019 axisCode << datasetId;
00020 axisCode << "-x";
00021 if (xAxisId < 10) axisCode << 0;
00022 axisCode << xAxisId;
00023 axisCode << "-y";
00024 if (yAxisId < 10) axisCode << 0;
00025 axisCode << yAxisId;
00026 return axisCode.str();
00027 }
00028 }
00029
00030
00031
00032
00033
00034 Analysis::Analysis(const string& name)
00035 : _crossSection(-1.0),
00036 _gotCrossSection(false),
00037 _needsCrossSection(false),
00038 _analysishandler(NULL),
00039 _madeHistoDir(false)
00040 {
00041 ProjectionApplier::_allowProjReg = false;
00042 _defaultname = name;
00043 AnalysisInfo* ai = AnalysisInfo::make(name);
00044 assert(ai != 0);
00045 _info.reset(ai);
00046 assert(_info.get() != 0);
00047
00048 }
00049
00050
00051 Analysis::~Analysis()
00052 { }
00053
00054
00055 IAnalysisFactory& Analysis::analysisFactory() {
00056 return handler().analysisFactory();
00057 }
00058
00059
00060 ITree& Analysis::tree() {
00061 return handler().tree();
00062 }
00063
00064
00065 IHistogramFactory& Analysis::histogramFactory() {
00066 return handler().histogramFactory();
00067 }
00068
00069
00070 IDataPointSetFactory& Analysis::datapointsetFactory() {
00071 return handler().datapointsetFactory();
00072 }
00073
00074
00075 double Analysis::sqrtS() const {
00076 return handler().sqrtS();
00077 }
00078
00079 const ParticlePair& Analysis::beams() const {
00080 return handler().beams();
00081 }
00082
00083 const BeamPair Analysis::beamIds() const {
00084 return handler().beamIds();
00085 }
00086
00087
00088 const string Analysis::histoDir() const {
00089 string path = "/" + name();
00090 if (handler().runName().length() > 0) {
00091 path = "/" + handler().runName() + path;
00092 }
00093 while (find_first(path, "//")) {
00094 replace_all(path, "//", "/");
00095 }
00096 return path;
00097 }
00098
00099
00100 const string Analysis::histoPath(const string& hname) const {
00101 const string path = histoDir() + "/" + hname;
00102 return path;
00103 }
00104
00105
00106 Log& Analysis::getLog() const {
00107 string logname = "Rivet.Analysis." + name();
00108 return Log::getLog(logname);
00109 }
00110
00111
00112 size_t Analysis::numEvents() const {
00113 return handler().numEvents();
00114 }
00115
00116
00117 double Analysis::sumOfWeights() const {
00118 return handler().sumOfWeights();
00119 }
00120
00121
00122
00123
00124
00125 const AnalysisInfo& Analysis::info() const {
00126 assert(_info.get() != 0);
00127 return *_info;
00128 }
00129
00130 string Analysis::name() const {
00131 if (_info && !_info->name().empty()) return _info->name();
00132 return _defaultname;
00133 }
00134
00135 string Analysis::spiresId() const {
00136 if (!_info) return "NONE";
00137 return _info->spiresId();
00138 }
00139
00140 vector<string> Analysis::authors() const {
00141 if (!_info) return std::vector<std::string>();
00142 return _info->authors();
00143 }
00144
00145 string Analysis::summary() const {
00146 if (!_info) return "NONE";
00147 return _info->summary();
00148 }
00149
00150 string Analysis::description() const {
00151 if (!_info) return "NONE";
00152 return _info->description();
00153 }
00154
00155 string Analysis::runInfo() const {
00156 if (!_info) return "NONE";
00157 return _info->runInfo();
00158 }
00159
00160 const std::vector<std::pair<double,double> >& Analysis::energies() const {
00161 return info().energies();
00162 }
00163
00164 string Analysis::experiment() const {
00165 if (!_info) return "NONE";
00166 return _info->experiment();
00167 }
00168
00169 string Analysis::collider() const {
00170 if (!_info) return "NONE";
00171 return _info->collider();
00172 }
00173
00174 string Analysis::year() const {
00175 if (!_info) return "NONE";
00176 return _info->year();
00177 }
00178
00179 vector<string> Analysis::references() const {
00180 if (!_info) return vector<string>();
00181 return _info->references();
00182 }
00183
00184 string Analysis::status() const {
00185 if (!_info) return "UNVALIDATED";
00186 return _info->status();
00187 }
00188
00189 const BeamPair Analysis::requiredBeams() const {
00190 return make_pdgid_pair(info().beams());
00191 }
00192
00193 Analysis& Analysis::setBeams(const ParticleName& beam1, const ParticleName& beam2) {
00194 assert(_info.get() != 0);
00195 _info->_beams = make_pair(beam1, beam2);
00196 return *this;
00197 }
00198
00199
00200 bool Analysis::isCompatible(const ParticleName& beam1, const ParticleName& beam2) const {
00201 BeamPair beams(beam1, beam2);
00202 return compatible(beams, requiredBeams());
00203
00204
00205 }
00206
00207 bool Analysis::isCompatible(const BeamPair& beams) const {
00208 return compatible(beams, requiredBeams());
00209
00210
00211 }
00212
00213
00214 Analysis& Analysis::setCrossSection(double xs) {
00215 _crossSection = xs;
00216 _gotCrossSection = true;
00217 return *this;
00218 }
00219
00220 bool Analysis::needsCrossSection() const {
00221 return _needsCrossSection;
00222 }
00223
00224 Analysis& Analysis::setNeedsCrossSection(bool needed) {
00225 _needsCrossSection = needed;
00226 return *this;
00227 }
00228
00229 double Analysis::crossSection() const {
00230 if (!_gotCrossSection || _crossSection < 0) {
00231 string errMsg = "You did not set the cross section for the analysis " + name();
00232 throw Error(errMsg);
00233 }
00234 return _crossSection;
00235 }
00236
00237 double Analysis::crossSectionPerEvent() const {
00238 const double sumW = sumOfWeights();
00239 assert(sumW > 0);
00240 return _crossSection / sumW;
00241 }
00242
00243
00244 AnalysisHandler& Analysis::handler() const {
00245 return *_analysishandler;
00246 }
00247
00248
00249
00250
00251
00252
00253
00254 void Analysis::_cacheBinEdges() const {
00255 _cacheXAxisData();
00256 if (_histBinEdges.empty()) {
00257 getLog() << Log::TRACE << "Getting histo bin edges from AIDA for paper " << name() << endl;
00258 _histBinEdges = getBinEdges(_dpsData);
00259 }
00260 }
00261
00262
00263 void Analysis::_cacheXAxisData() const {
00264 if (_dpsData.empty()) {
00265 getLog() << Log::TRACE << "Getting DPS x-axis data from AIDA for paper " << name() << endl;
00266 _dpsData = getDPSXValsErrs(name());
00267 }
00268 }
00269
00270
00271 const BinEdges& Analysis::binEdges(const string& hname) const {
00272 _cacheBinEdges();
00273 getLog() << Log::TRACE << "Using histo bin edges for " << name() << ":" << hname << endl;
00274 const BinEdges& edges = _histBinEdges.find(hname)->second;
00275 if (getLog().isActive(Log::TRACE)) {
00276 stringstream edges_ss;
00277 foreach (const double be, edges) {
00278 edges_ss << " " << be;
00279 }
00280 getLog() << Log::TRACE << "Edges:" << edges_ss.str() << endl;
00281 }
00282 return edges;
00283 }
00284
00285
00286 const BinEdges& Analysis::binEdges(size_t datasetId, size_t xAxisId, size_t yAxisId) const {
00287 const string hname = makeAxisCode(datasetId, xAxisId, yAxisId);
00288 return binEdges(hname);
00289 }
00290
00291
00292 BinEdges Analysis::logBinEdges(size_t nbins, double lower, double upper) {
00293 assert(lower>0.0);
00294 assert(upper>lower);
00295 double loglower=log10(lower);
00296 double logupper=log10(upper);
00297 vector<double> binedges;
00298 double stepwidth=(logupper-loglower)/double(nbins);
00299 for (size_t i=0; i<=nbins; ++i) {
00300 binedges.push_back(pow(10.0, loglower+double(i)*stepwidth));
00301 }
00302 return binedges;
00303 }
00304
00305 IHistogram1D* Analysis::bookHistogram1D(size_t datasetId, size_t xAxisId,
00306 size_t yAxisId, const string& title,
00307 const string& xtitle, const string& ytitle)
00308 {
00309 const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
00310 return bookHistogram1D(axisCode, title, xtitle, ytitle);
00311 }
00312
00313
00314 IHistogram1D* Analysis::bookHistogram1D(const string& hname, const string& title,
00315 const string& xtitle, const string& ytitle)
00316 {
00317
00318 const BinEdges edges = binEdges(hname);
00319 _makeHistoDir();
00320 const string path = histoPath(hname);
00321 IHistogram1D* hist = histogramFactory().createHistogram1D(path, title, edges);
00322 getLog() << Log::TRACE << "Made histogram " << hname << " for " << name() << endl;
00323 hist->setXTitle(xtitle);
00324 hist->setYTitle(ytitle);
00325 return hist;
00326 }
00327
00328
00329 IHistogram1D* Analysis::bookHistogram1D(const string& hname,
00330 size_t nbins, double lower, double upper,
00331 const string& title,
00332 const string& xtitle, const string& ytitle) {
00333 _makeHistoDir();
00334 const string path = histoPath(hname);
00335 IHistogram1D* hist = histogramFactory().createHistogram1D(path, title, nbins, lower, upper);
00336 getLog() << Log::TRACE << "Made histogram " << hname << " for " << name() << endl;
00337 hist->setXTitle(xtitle);
00338 hist->setYTitle(ytitle);
00339 return hist;
00340 }
00341
00342
00343 IHistogram1D* Analysis::bookHistogram1D(const string& hname,
00344 const vector<double>& binedges,
00345 const string& title,
00346 const string& xtitle, const string& ytitle) {
00347 _makeHistoDir();
00348 const string path = histoPath(hname);
00349 IHistogram1D* hist = histogramFactory().createHistogram1D(path, title, binedges);
00350 getLog() << Log::TRACE << "Made histogram " << hname << " for " << name() << endl;
00351 hist->setXTitle(xtitle);
00352 hist->setYTitle(ytitle);
00353 return hist;
00354 }
00355
00356
00357
00358
00359
00360 IProfile1D* Analysis::bookProfile1D(size_t datasetId, size_t xAxisId,
00361 size_t yAxisId, const string& title,
00362 const string& xtitle, const string& ytitle) {
00363 const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
00364 return bookProfile1D(axisCode, title, xtitle, ytitle);
00365 }
00366
00367
00368 IProfile1D* Analysis::bookProfile1D(const string& hname, const string& title,
00369 const string& xtitle, const string& ytitle)
00370 {
00371
00372 const BinEdges edges = binEdges(hname);
00373 _makeHistoDir();
00374 const string path = histoPath(hname);
00375 IProfile1D* prof = histogramFactory().createProfile1D(path, title, edges);
00376 getLog() << Log::TRACE << "Made profile histogram " << hname << " for " << name() << endl;
00377 prof->setXTitle(xtitle);
00378 prof->setYTitle(ytitle);
00379 return prof;
00380 }
00381
00382
00383 IProfile1D* Analysis::bookProfile1D(const string& hname,
00384 size_t nbins, double lower, double upper,
00385 const string& title,
00386 const string& xtitle, const string& ytitle) {
00387 _makeHistoDir();
00388 const string path = histoPath(hname);
00389 IProfile1D* prof = histogramFactory().createProfile1D(path, title, nbins, lower, upper);
00390 getLog() << Log::TRACE << "Made profile histogram " << hname << " for " << name() << endl;
00391 prof->setXTitle(xtitle);
00392 prof->setYTitle(ytitle);
00393 return prof;
00394 }
00395
00396
00397 IProfile1D* Analysis::bookProfile1D(const string& hname,
00398 const vector<double>& binedges,
00399 const string& title,
00400 const string& xtitle, const string& ytitle) {
00401 _makeHistoDir();
00402 const string path = histoPath(hname);
00403 IProfile1D* prof = histogramFactory().createProfile1D(path, title, binedges);
00404 getLog() << Log::TRACE << "Made profile histogram " << hname << " for " << name() << endl;
00405 prof->setXTitle(xtitle);
00406 prof->setYTitle(ytitle);
00407 return prof;
00408 }
00409
00410
00411
00412
00413
00414
00415 IDataPointSet* Analysis::bookDataPointSet(const string& hname, const string& title,
00416 const string& xtitle, const string& ytitle) {
00417 _makeHistoDir();
00418 const string path = histoPath(hname);
00419 IDataPointSet* dps = datapointsetFactory().create(path, title, 2);
00420 getLog() << Log::TRACE << "Made data point set " << hname << " for " << name() << endl;
00421 dps->setXTitle(xtitle);
00422 dps->setYTitle(ytitle);
00423 return dps;
00424 }
00425
00426
00427 IDataPointSet* Analysis::bookDataPointSet(const string& hname,
00428 size_t npts, double lower, double upper,
00429 const string& title,
00430 const string& xtitle, const string& ytitle) {
00431 IDataPointSet* dps = bookDataPointSet(hname, title, xtitle, ytitle);
00432 for (size_t pt = 0; pt < npts; ++pt) {
00433 const double binwidth = (upper-lower)/npts;
00434 const double bincentre = lower + (pt + 0.5) * binwidth;
00435 dps->addPoint();
00436 IMeasurement* meas = dps->point(pt)->coordinate(0);
00437 meas->setValue(bincentre);
00438 meas->setErrorPlus(binwidth/2.0);
00439 meas->setErrorMinus(binwidth/2.0);
00440 }
00441 return dps;
00442 }
00443
00444
00445 IDataPointSet* Analysis::bookDataPointSet(size_t datasetId, size_t xAxisId,
00446 size_t yAxisId, const string& title,
00447 const string& xtitle, const string& ytitle) {
00448
00449 _cacheXAxisData();
00450
00451 const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
00452
00453 getLog() << Log::TRACE << "Using DPS x-positions for " << name() << ":" << axisCode << endl;
00454 IDataPointSet* dps = bookDataPointSet(axisCode, title, xtitle, ytitle);
00455 const vector<DPSXPoint> xpts = _dpsData.find(axisCode)->second;
00456 for (size_t pt = 0; pt < xpts.size(); ++pt) {
00457 dps->addPoint();
00458 IMeasurement* meas = dps->point(pt)->coordinate(0);
00459 meas->setValue(xpts[pt].val);
00460 meas->setErrorPlus(xpts[pt].errplus);
00461 meas->setErrorMinus(xpts[pt].errminus);
00462 }
00463 getLog() << Log::TRACE << "Made DPS " << axisCode << " for " << name() << endl;
00464 return dps;
00465 }
00466
00467
00468
00469
00470
00471 void Analysis::_makeHistoDir() {
00472 if (!_madeHistoDir) {
00473 if (! name().empty()) {
00474
00475
00476
00477
00478
00479
00480 tree().mkdirs(histoDir());
00481 }
00482 _madeHistoDir = true;
00483 }
00484 }
00485
00486
00487 void Analysis::normalize(AIDA::IHistogram1D*& histo, double norm) {
00488 if (!histo) {
00489 getLog() << Log::ERROR << "Failed to normalise histo=NULL in analysis "
00490 << name() << "(norm=" << norm << ")" << endl;
00491 return;
00492 }
00493 const string hpath = tree().findPath(dynamic_cast<const AIDA::IManagedObject&>(*histo));
00494 getLog() << Log::TRACE << "Normalizing histo " << hpath << " to " << norm << endl;
00495
00496 double oldintg = 0.0;
00497 int nBins = histo->axis().bins();
00498 for (int iBin = 0; iBin != nBins; ++iBin) {
00499
00500 oldintg += histo->binHeight(iBin);
00501 }
00502 if (oldintg == 0.0) {
00503 getLog() << Log::WARN << "Histo " << hpath << " has null integral during normalisation" << endl;
00504 return;
00505 }
00506
00507
00508 scale(histo, norm/oldintg);
00509 }
00510
00511
00512 void Analysis::scale(AIDA::IHistogram1D*& histo, double scale) {
00513 if (!histo) {
00514 getLog() << Log::ERROR << "Failed to scale histo=NULL in analysis "
00515 << name() << "(scale=" << scale << ")" << endl;
00516 return;
00517 }
00518 const string hpath = tree().findPath(dynamic_cast<const AIDA::IManagedObject&>(*histo));
00519 getLog() << Log::TRACE << "Scaling histo " << hpath << endl;
00520
00521 vector<double> x, y, ex, ey;
00522 for (size_t i = 0, N = histo->axis().bins(); i < N; ++i) {
00523 x.push_back(0.5 * (histo->axis().binLowerEdge(i) + histo->axis().binUpperEdge(i)));
00524 ex.push_back(histo->axis().binWidth(i)*0.5);
00525
00526
00527
00528 y.push_back(histo->binHeight(i)*scale/histo->axis().binWidth(i));
00529
00530
00531
00532 ey.push_back(histo->binError(i)*scale/(0.5*histo->axis().binWidth(i)));
00533 }
00534
00535 string title = histo->title();
00536 string xtitle = histo->xtitle();
00537 string ytitle = histo->ytitle();
00538
00539 tree().mkdir("/tmpnormalize");
00540 tree().mv(hpath, "/tmpnormalize");
00541
00542 AIDA::IDataPointSet* dps = datapointsetFactory().createXY(hpath, title, x, y, ex, ey);
00543 dps->setXTitle(xtitle);
00544 dps->setYTitle(ytitle);
00545
00546 tree().rm(tree().findPath(dynamic_cast<AIDA::IManagedObject&>(*histo)));
00547 tree().rmdir("/tmpnormalize");
00548
00549
00550 histo = 0;
00551 }
00552
00553
00554 }