rivet is hosted by Hepforge, IPPP Durham
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 }