Analysis.cc
Go to the documentation of this file.
00001 // -*- C++ -*- 00002 #include "Rivet/Config/RivetCommon.hh" 00003 #include "Rivet/Analysis.hh" 00004 #include "Rivet/AnalysisHandler.hh" 00005 #include "Rivet/AnalysisInfo.hh" 00006 #include "Rivet/BeamConstraint.hh" 00007 #include "Rivet/Tools/RivetYODA.hh" 00008 #include "Rivet/Tools/Logging.hh" 00009 00010 namespace Rivet { 00011 00012 00013 Analysis::Analysis(const string& name) 00014 : _crossSection(-1.0), 00015 _gotCrossSection(false), 00016 _analysishandler(NULL) 00017 { 00018 ProjectionApplier::_allowProjReg = false; 00019 _defaultname = name; 00020 00021 AnalysisInfo* ai = AnalysisInfo::make(name); 00022 assert(ai != 0); 00023 _info.reset(ai); 00024 assert(_info.get() != 0); 00025 } 00026 00027 double Analysis::sqrtS() const { 00028 return handler().sqrtS(); 00029 } 00030 00031 const ParticlePair& Analysis::beams() const { 00032 return handler().beams(); 00033 } 00034 00035 const PdgIdPair Analysis::beamIds() const { 00036 return handler().beamIds(); 00037 } 00038 00039 00040 const string Analysis::histoDir() const { 00041 /// @todo Cache in a member variable 00042 string _histoDir; 00043 if (_histoDir.empty()) { 00044 _histoDir = "/" + name(); 00045 if (handler().runName().length() > 0) { 00046 _histoDir = "/" + handler().runName() + _histoDir; 00047 } 00048 while (find_first(_histoDir, "//")) { 00049 replace_all(_histoDir, "//", "/"); 00050 } 00051 } 00052 return _histoDir; 00053 } 00054 00055 00056 const string Analysis::histoPath(const string& hname) const { 00057 const string path = histoDir() + "/" + hname; 00058 return path; 00059 } 00060 00061 00062 const string Analysis::histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { 00063 return histoDir() + "/" + makeAxisCode(datasetId, xAxisId, yAxisId); 00064 } 00065 00066 00067 const string Analysis::makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { 00068 stringstream axisCode; 00069 axisCode << "d"; 00070 if (datasetId < 10) axisCode << 0; 00071 axisCode << datasetId; 00072 axisCode << "-x"; 00073 if (xAxisId < 10) axisCode << 0; 00074 axisCode << xAxisId; 00075 axisCode << "-y"; 00076 if (yAxisId < 10) axisCode << 0; 00077 axisCode << yAxisId; 00078 return axisCode.str(); 00079 } 00080 00081 00082 Log& Analysis::getLog() const { 00083 string logname = "Rivet.Analysis." + name(); 00084 return Log::getLog(logname); 00085 } 00086 00087 00088 size_t Analysis::numEvents() const { 00089 return handler().numEvents(); 00090 } 00091 00092 00093 double Analysis::sumOfWeights() const { 00094 return handler().sumOfWeights(); 00095 } 00096 00097 00098 /////////////////////////////////////////// 00099 00100 00101 bool Analysis::isCompatible(const ParticlePair& beams) const { 00102 return isCompatible(beams.first.pid(), beams.second.pid(), 00103 beams.first.energy(), beams.second.energy()); 00104 } 00105 00106 00107 bool Analysis::isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const { 00108 PdgIdPair beams(beam1, beam2); 00109 pair<double,double> energies(e1, e2); 00110 return isCompatible(beams, energies); 00111 } 00112 00113 00114 bool Analysis::isCompatible(const PdgIdPair& beams, const pair<double,double>& energies) const { 00115 // First check the beam IDs 00116 bool beamIdsOk = false; 00117 foreach (const PdgIdPair& bp, requiredBeams()) { 00118 if (compatible(beams, bp)) { 00119 beamIdsOk = true; 00120 break; 00121 } 00122 } 00123 if (!beamIdsOk) return false; 00124 00125 // Next check that the energies are compatible (within 1% or 1 GeV, whichever is larger, for a bit of UI forgiveness) 00126 00127 /// @todo Use some sort of standard ordering to improve comparisons, esp. when the two beams are different particles 00128 bool beamEnergiesOk = requiredEnergies().size() > 0 ? false : true; 00129 typedef pair<double,double> DoublePair; 00130 foreach (const DoublePair& ep, requiredEnergies()) { 00131 if ((fuzzyEquals(ep.first, energies.first, 0.01) && fuzzyEquals(ep.second, energies.second, 0.01)) || 00132 (fuzzyEquals(ep.first, energies.second, 0.01) && fuzzyEquals(ep.second, energies.first, 0.01)) || 00133 (abs(ep.first - energies.first) < 1*GeV && abs(ep.second - energies.second) < 1*GeV) || 00134 (abs(ep.first - energies.second) < 1*GeV && abs(ep.second - energies.first) < 1*GeV)) { 00135 beamEnergiesOk = true; 00136 break; 00137 } 00138 } 00139 return beamEnergiesOk; 00140 00141 /// @todo Need to also check internal consistency of the analysis' 00142 /// beam requirements with those of the projections it uses. 00143 } 00144 00145 00146 /////////////////////////////////////////// 00147 00148 00149 Analysis& Analysis::setCrossSection(double xs) { 00150 _crossSection = xs; 00151 _gotCrossSection = true; 00152 return *this; 00153 } 00154 00155 double Analysis::crossSection() const { 00156 if (!_gotCrossSection || std::isnan(_crossSection)) { 00157 string errMsg = "You did not set the cross section for the analysis " + name(); 00158 throw Error(errMsg); 00159 } 00160 return _crossSection; 00161 } 00162 00163 double Analysis::crossSectionPerEvent() const { 00164 const double sumW = sumOfWeights(); 00165 assert(sumW != 0.0); 00166 return _crossSection / sumW; 00167 } 00168 00169 00170 00171 //////////////////////////////////////////////////////////// 00172 // Histogramming 00173 00174 00175 void Analysis::_cacheRefData() const { 00176 if (_refdata.empty()) { 00177 MSG_TRACE("Getting refdata cache for paper " << name()); 00178 _refdata = getRefData(name()); 00179 } 00180 } 00181 00182 00183 const Scatter2D& Analysis::refData(const string& hname) const { 00184 _cacheRefData(); 00185 MSG_TRACE("Using histo bin edges for " << name() << ":" << hname); 00186 if (!_refdata[hname]) { 00187 MSG_ERROR("Can't find reference histogram " << hname); 00188 throw Exception("Reference data " + hname + " not found."); 00189 } 00190 return dynamic_cast<Scatter2D&>(*_refdata[hname]); 00191 } 00192 00193 00194 const Scatter2D& Analysis::refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { 00195 const string hname = makeAxisCode(datasetId, xAxisId, yAxisId); 00196 return refData(hname); 00197 } 00198 00199 00200 Histo1DPtr Analysis::bookHisto1D(const string& hname, 00201 size_t nbins, double lower, double upper, 00202 const string& title, 00203 const string& xtitle, 00204 const string& ytitle) { 00205 const string path = histoPath(hname); 00206 Histo1DPtr hist( new Histo1D(nbins, lower, upper, path, title) ); 00207 addAnalysisObject(hist); 00208 MSG_TRACE("Made histogram " << hname << " for " << name()); 00209 hist->setAnnotation("XLabel", xtitle); 00210 hist->setAnnotation("YLabel", ytitle); 00211 return hist; 00212 } 00213 00214 00215 Histo1DPtr Analysis::bookHisto1D(const string& hname, 00216 const vector<double>& binedges, 00217 const string& title, 00218 const string& xtitle, 00219 const string& ytitle) { 00220 const string path = histoPath(hname); 00221 Histo1DPtr hist( new Histo1D(binedges, path, title) ); 00222 addAnalysisObject(hist); 00223 MSG_TRACE("Made histogram " << hname << " for " << name()); 00224 hist->setAnnotation("XLabel", xtitle); 00225 hist->setAnnotation("YLabel", ytitle); 00226 return hist; 00227 } 00228 00229 00230 Histo1DPtr Analysis::bookHisto1D(const string& hname, 00231 const Scatter2D& refscatter, 00232 const string& title, 00233 const string& xtitle, 00234 const string& ytitle) { 00235 const string path = histoPath(hname); 00236 Histo1DPtr hist( new Histo1D(refscatter, path) ); 00237 addAnalysisObject(hist); 00238 MSG_TRACE("Made histogram " << hname << " for " << name()); 00239 hist->setTitle(title); 00240 hist->setAnnotation("XLabel", xtitle); 00241 hist->setAnnotation("YLabel", ytitle); 00242 return hist; 00243 } 00244 00245 00246 Histo1DPtr Analysis::bookHisto1D(const string& hname, 00247 const string& title, 00248 const string& xtitle, 00249 const string& ytitle) { 00250 const Scatter2D& refdata = refData(hname); 00251 return bookHisto1D(hname, refdata, title, xtitle, ytitle); 00252 } 00253 00254 00255 Histo1DPtr Analysis::bookHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, 00256 const string& title, 00257 const string& xtitle, 00258 const string& ytitle) { 00259 const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId); 00260 return bookHisto1D(axisCode, title, xtitle, ytitle); 00261 } 00262 00263 00264 /// @todo Add booking methods which take a path, titles and *a reference Scatter from which to book* 00265 00266 00267 ///////////////// 00268 00269 00270 Histo2DPtr Analysis::bookHisto2D(const string& hname, 00271 size_t nxbins, double xlower, double xupper, 00272 size_t nybins, double ylower, double yupper, 00273 const string& title, 00274 const string& xtitle, 00275 const string& ytitle, 00276 const string& ztitle) 00277 { 00278 const string path = histoPath(hname); 00279 Histo2DPtr hist( new Histo2D(nxbins, xlower, xupper, nybins, ylower, yupper, path, title) ); 00280 addAnalysisObject(hist); 00281 MSG_TRACE("Made 2D histogram " << hname << " for " << name()); 00282 hist->setAnnotation("XLabel", xtitle); 00283 hist->setAnnotation("YLabel", ytitle); 00284 hist->setAnnotation("ZLabel", ztitle); 00285 return hist; 00286 } 00287 00288 00289 Histo2DPtr Analysis::bookHisto2D(const string& hname, 00290 const vector<double>& xbinedges, 00291 const vector<double>& ybinedges, 00292 const string& title, 00293 const string& xtitle, 00294 const string& ytitle, 00295 const string& ztitle) 00296 { 00297 const string path = histoPath(hname); 00298 Histo2DPtr hist( new Histo2D(xbinedges, ybinedges, path, title) ); 00299 addAnalysisObject(hist); 00300 MSG_TRACE("Made 2D histogram " << hname << " for " << name()); 00301 hist->setAnnotation("XLabel", xtitle); 00302 hist->setAnnotation("YLabel", ytitle); 00303 hist->setAnnotation("ZLabel", ztitle); 00304 return hist; 00305 } 00306 00307 00308 // Histo2DPtr Analysis::bookHisto2D(const string& hname, 00309 // const Scatter3D& refscatter, 00310 // const string& title="", 00311 // const string& xtitle="", 00312 // const string& ytitle="", 00313 // const string& ztitle="") { 00314 // const string path = histoPath(hname); 00315 // Histo2DPtr hist( new Histo2D(refscatter, path) ); 00316 // addAnalysisObject(hist); 00317 // MSG_TRACE("Made 2D histogram " << hname << " for " << name()); 00318 // hist->setTitle(title); 00319 // hist->setAnnotation("XLabel", xtitle); 00320 // hist->setAnnotation("YLabel", ytitle); 00321 // hist->setAnnotation("ZLabel", ztitle); 00322 // return hist; 00323 // } 00324 00325 00326 // Histo2DPtr Analysis::bookHisto2D(const string& hname, 00327 // const string& title, 00328 // const string& xtitle, 00329 // const string& ytitle, 00330 // const string& ztitle) { 00331 // const Scatter3D& refdata = refData(hname); 00332 // return bookHisto2D(hname, refdata, title, xtitle, ytitle, ztitle); 00333 // } 00334 00335 00336 // Histo2DPtr Analysis::bookHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, 00337 // const string& title, 00338 // const string& xtitle, 00339 // const string& ytitle, 00340 // const string& ztitle) { 00341 // const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId); 00342 // return bookHisto2D(axisCode, title, xtitle, ytitle, ztitle); 00343 // } 00344 00345 00346 ///////////////// 00347 00348 00349 Profile1DPtr Analysis::bookProfile1D(const string& hname, 00350 size_t nbins, double lower, double upper, 00351 const string& title, 00352 const string& xtitle, 00353 const string& ytitle) { 00354 const string path = histoPath(hname); 00355 Profile1DPtr prof( new Profile1D(nbins, lower, upper, path, title) ); 00356 addAnalysisObject(prof); 00357 MSG_TRACE("Made profile histogram " << hname << " for " << name()); 00358 prof->setAnnotation("XLabel", xtitle); 00359 prof->setAnnotation("YLabel", ytitle); 00360 return prof; 00361 } 00362 00363 00364 Profile1DPtr Analysis::bookProfile1D(const string& hname, 00365 const vector<double>& binedges, 00366 const string& title, 00367 const string& xtitle, 00368 const string& ytitle) { 00369 const string path = histoPath(hname); 00370 Profile1DPtr prof( new Profile1D(binedges, path, title) ); 00371 addAnalysisObject(prof); 00372 MSG_TRACE("Made profile histogram " << hname << " for " << name()); 00373 prof->setAnnotation("XLabel", xtitle); 00374 prof->setAnnotation("YLabel", ytitle); 00375 return prof; 00376 } 00377 00378 00379 Profile1DPtr Analysis::bookProfile1D(const string& hname, 00380 const Scatter2D& refscatter, 00381 const string& title, 00382 const string& xtitle, 00383 const string& ytitle) { 00384 const string path = histoPath(hname); 00385 Profile1DPtr prof( new Profile1D(refscatter, path) ); 00386 addAnalysisObject(prof); 00387 MSG_TRACE("Made profile histogram " << hname << " for " << name()); 00388 prof->setTitle(title); 00389 prof->setAnnotation("XLabel", xtitle); 00390 prof->setAnnotation("YLabel", ytitle); 00391 return prof; 00392 } 00393 00394 00395 Profile1DPtr Analysis::bookProfile1D(const string& hname, 00396 const string& title, 00397 const string& xtitle, 00398 const string& ytitle) { 00399 const Scatter2D& refdata = refData(hname); 00400 return bookProfile1D(hname, refdata, title, xtitle, ytitle); 00401 } 00402 00403 00404 Profile1DPtr Analysis::bookProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, 00405 const string& title, 00406 const string& xtitle, 00407 const string& ytitle) { 00408 const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId); 00409 return bookProfile1D(axisCode, title, xtitle, ytitle); 00410 } 00411 00412 00413 /////////////////// 00414 00415 00416 00417 Profile2DPtr Analysis::bookProfile2D(const string& hname, 00418 size_t nxbins, double xlower, double xupper, 00419 size_t nybins, double ylower, double yupper, 00420 const string& title, 00421 const string& xtitle, 00422 const string& ytitle, 00423 const string& ztitle) 00424 { 00425 const string path = histoPath(hname); 00426 Profile2DPtr prof( new Profile2D(nxbins, xlower, xupper, nybins, ylower, yupper, path, title) ); 00427 addAnalysisObject(prof); 00428 MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); 00429 prof->setAnnotation("XLabel", xtitle); 00430 prof->setAnnotation("YLabel", ytitle); 00431 prof->setAnnotation("ZLabel", ztitle); 00432 return prof; 00433 } 00434 00435 00436 Profile2DPtr Analysis::bookProfile2D(const string& hname, 00437 const vector<double>& xbinedges, 00438 const vector<double>& ybinedges, 00439 const string& title, 00440 const string& xtitle, 00441 const string& ytitle, 00442 const string& ztitle) 00443 { 00444 const string path = histoPath(hname); 00445 Profile2DPtr prof( new Profile2D(xbinedges, ybinedges, path, title) ); 00446 addAnalysisObject(prof); 00447 MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); 00448 prof->setAnnotation("XLabel", xtitle); 00449 prof->setAnnotation("YLabel", ytitle); 00450 prof->setAnnotation("ZLabel", ztitle); 00451 return prof; 00452 } 00453 00454 00455 // Profile2DPtr Analysis::bookProfile2D(const string& hname, 00456 // const Scatter3D& refscatter, 00457 // const string& title="", 00458 // const string& xtitle="", 00459 // const string& ytitle="", 00460 // const string& ztitle="") { 00461 // const string path = histoPath(hname); 00462 // Profile2DPtr prof( new Profile2D(refscatter, path) ); 00463 // addAnalysisObject(prof); 00464 // MSG_TRACE("Made 2D profile histogram " << hname << " for " << name()); 00465 // prof->setTitle(title); 00466 // prof->setAnnotation("XLabel", xtitle); 00467 // prof->setAnnotation("YLabel", ytitle); 00468 // prof->setAnnotation("ZLabel", ztitle); 00469 // return prof; 00470 // } 00471 00472 00473 // Profile2DPtr Analysis::bookProfile2D(const string& hname, 00474 // const string& title, 00475 // const string& xtitle, 00476 // const string& ytitle, 00477 // const string& ztitle) { 00478 // const Scatter3D& refdata = refData(hname); 00479 // return bookProfile2D(hname, refdata, title, xtitle, ytitle, ztitle); 00480 // } 00481 00482 00483 // Profile2DPtr Analysis::bookProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, 00484 // const string& title, 00485 // const string& xtitle, 00486 // const string& ytitle, 00487 // const string& ztitle) { 00488 // const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId); 00489 // return bookProfile2D(axisCode, title, xtitle, ytitle, ztitle); 00490 // } 00491 00492 00493 ///////////////// 00494 00495 00496 Scatter2DPtr Analysis::bookScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, 00497 bool copy_pts, 00498 const string& title, 00499 const string& xtitle, 00500 const string& ytitle) { 00501 const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId); 00502 return bookScatter2D(axisCode, copy_pts, title, xtitle, ytitle); 00503 } 00504 00505 00506 Scatter2DPtr Analysis::bookScatter2D(const string& hname, 00507 bool copy_pts, 00508 const string& title, 00509 const string& xtitle, 00510 const string& ytitle) { 00511 Scatter2DPtr s; 00512 const string path = histoPath(hname); 00513 if (copy_pts) { 00514 const Scatter2D& refdata = refData(hname); 00515 s.reset( new Scatter2D(refdata, path) ); 00516 foreach (Point2D& p, s->points()) p.setY(0, 0); 00517 } else { 00518 s.reset( new Scatter2D(path) ); 00519 } 00520 addAnalysisObject(s); 00521 MSG_TRACE("Made scatter " << hname << " for " << name()); 00522 s->setTitle(title); 00523 s->setAnnotation("XLabel", xtitle); 00524 s->setAnnotation("YLabel", ytitle); 00525 return s; 00526 } 00527 00528 00529 Scatter2DPtr Analysis::bookScatter2D(const string& hname, 00530 size_t npts, double lower, double upper, 00531 const string& title, 00532 const string& xtitle, 00533 const string& ytitle) { 00534 const string path = histoPath(hname); 00535 Scatter2DPtr s( new Scatter2D(path) ); 00536 const double binwidth = (upper-lower)/npts; 00537 for (size_t pt = 0; pt < npts; ++pt) { 00538 const double bincentre = lower + (pt + 0.5) * binwidth; 00539 s->addPoint(bincentre, 0, binwidth/2.0, 0); 00540 } 00541 addAnalysisObject(s); 00542 MSG_TRACE("Made scatter " << hname << " for " << name()); 00543 s->setTitle(title); 00544 s->setAnnotation("XLabel", xtitle); 00545 s->setAnnotation("YLabel", ytitle); 00546 return s; 00547 } 00548 00549 00550 Scatter2DPtr Analysis::bookScatter2D(const string& hname, 00551 const vector<double>& binedges, 00552 const string& title, 00553 const string& xtitle, 00554 const string& ytitle) { 00555 const string path = histoPath(hname); 00556 Scatter2DPtr s( new Scatter2D(path) ); 00557 for (size_t pt = 0; pt < binedges.size()-1; ++pt) { 00558 const double bincentre = (binedges[pt] + binedges[pt+1]) / 2.0; 00559 const double binwidth = binedges[pt+1] - binedges[pt]; 00560 s->addPoint(bincentre, 0, binwidth/2.0, 0); 00561 } 00562 addAnalysisObject(s); 00563 MSG_TRACE("Made scatter " << hname << " for " << name()); 00564 s->setTitle(title); 00565 s->setAnnotation("XLabel", xtitle); 00566 s->setAnnotation("YLabel", ytitle); 00567 return s; 00568 } 00569 00570 00571 ///////////////////// 00572 00573 00574 void Analysis::divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { 00575 const string path = s->path(); 00576 *s = *h1 / *h2; 00577 s->setPath(path); 00578 } 00579 00580 void Analysis::divide(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { 00581 const string path = s->path(); 00582 *s = h1 / h2; 00583 s->setPath(path); 00584 } 00585 00586 00587 void Analysis::divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const { 00588 const string path = s->path(); 00589 *s = *p1 / *p2; 00590 s->setPath(path); 00591 } 00592 00593 void Analysis::divide(const Profile1D& p1, const Profile1D& p2, Scatter2DPtr s) const { 00594 const string path = s->path(); 00595 *s = p1 / p2; 00596 s->setPath(path); 00597 } 00598 00599 00600 void Analysis::divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const { 00601 const string path = s->path(); 00602 *s = *h1 / *h2; 00603 s->setPath(path); 00604 } 00605 00606 void Analysis::divide(const Histo2D& h1, const Histo2D& h2, Scatter3DPtr s) const { 00607 const string path = s->path(); 00608 *s = h1 / h2; 00609 s->setPath(path); 00610 } 00611 00612 00613 void Analysis::divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const { 00614 const string path = s->path(); 00615 *s = *p1 / *p2; 00616 s->setPath(path); 00617 } 00618 00619 void Analysis::divide(const Profile2D& p1, const Profile2D& p2, Scatter3DPtr s) const { 00620 const string path = s->path(); 00621 *s = p1 / p2; 00622 s->setPath(path); 00623 } 00624 00625 00626 void Analysis::efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { 00627 const string path = s->path(); 00628 *s = YODA::efficiency(*h1, *h2); 00629 s->setPath(path); 00630 } 00631 00632 void Analysis::efficiency(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { 00633 const string path = s->path(); 00634 *s = YODA::efficiency(h1, h2); 00635 s->setPath(path); 00636 } 00637 00638 00639 void Analysis::asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const { 00640 const string path = s->path(); 00641 *s = YODA::asymm(*h1, *h2); 00642 s->setPath(path); 00643 } 00644 00645 void Analysis::asymm(const Histo1D& h1, const Histo1D& h2, Scatter2DPtr s) const { 00646 const string path = s->path(); 00647 *s = YODA::asymm(h1, h2); 00648 s->setPath(path); 00649 } 00650 00651 00652 void Analysis::normalize(Histo1DPtr histo, double norm, bool includeoverflows) { 00653 if (!histo) { 00654 MSG_ERROR("Failed to normalize histo=NULL in analysis " << name() << " (norm=" << norm << ")"); 00655 return; 00656 } 00657 MSG_TRACE("Normalizing histo " << histo->path() << " to " << norm); 00658 try { 00659 histo->normalize(norm, includeoverflows); 00660 } catch (YODA::Exception& we) { 00661 MSG_WARNING("Could not normalize histo " << histo->path()); 00662 return; 00663 } 00664 } 00665 00666 00667 void Analysis::scale(Histo1DPtr histo, double scale) { 00668 if (!histo) { 00669 MSG_ERROR("Failed to scale histo=NULL in analysis " << name() << " (scale=" << scale << ")"); 00670 return; 00671 } 00672 if (std::isnan(scale) || std::isinf(scale)) { 00673 MSG_ERROR("Failed to scale histo=" << histo->path() << " in analysis: " << name() << " (invalid scale factor = " << scale << ")"); 00674 scale = 0; 00675 } 00676 MSG_TRACE("Scaling histo " << histo->path() << " by factor " << scale); 00677 try { 00678 histo->scaleW(scale); 00679 } catch (YODA::Exception& we) { 00680 MSG_WARNING("Could not scale histo " << histo->path()); 00681 return; 00682 } 00683 // // Transforming the histo into a scatter after scaling 00684 // vector<double> x, y, ex, ey; 00685 // for (size_t i = 0, N = histo->numBins(); i < N; ++i) { 00686 // x.push_back( histo->bin(i).midpoint() ); 00687 // ex.push_back(histo->bin(i).width()*0.5); 00688 // y.push_back(histo->bin(i).height()*scale); 00689 // ey.push_back(histo->bin(i).heightErr()*scale); 00690 // } 00691 // string title = histo->title(); 00692 // Scatter2DPtr dps( new Scatter2D(x, y, ex, ey, hpath, title) ); 00693 // addAnalysisObject(dps); 00694 } 00695 00696 00697 /// @todo 2D versions of scale and normalize... 00698 00699 00700 void Analysis::integrate(Histo1DPtr h, Scatter2DPtr s) const { 00701 // preserve the path info 00702 const string path = s->path(); 00703 *s = toIntegralHisto(*h); 00704 s->setPath(path); 00705 } 00706 00707 void Analysis::integrate(const Histo1D& h, Scatter2DPtr s) const { 00708 // preserve the path info 00709 const string path = s->path(); 00710 *s = toIntegralHisto(h); 00711 s->setPath(path); 00712 } 00713 00714 00715 /// @todo 2D versions of integrate... defined how, exactly?!? 00716 00717 00718 ////////////////////////////////// 00719 00720 00721 void Analysis::addAnalysisObject(AnalysisObjectPtr ao) { 00722 _analysisobjects.push_back(ao); 00723 } 00724 00725 void Analysis::removeAnalysisObject(const string& path) { 00726 for (vector<AnalysisObjectPtr>::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) { 00727 if ((*it)->path() == path) { 00728 _analysisobjects.erase(it); 00729 break; 00730 } 00731 } 00732 } 00733 00734 void Analysis::removeAnalysisObject(AnalysisObjectPtr ao) { 00735 for (vector<AnalysisObjectPtr>::iterator it = _analysisobjects.begin(); it != _analysisobjects.end(); ++it) { 00736 if (*it == ao) { 00737 _analysisobjects.erase(it); 00738 break; 00739 } 00740 } 00741 } 00742 00743 00744 } Generated on Wed Oct 7 2015 12:09:09 for The Rivet MC analysis system by ![]() |