Analysis.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
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     //setBeams(ANY, ANY);
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   // Metadata
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     /// @todo Need to also check internal consistency of the analysis'
00204     /// beam requirements with those of the projections it uses.
00205   }
00206 
00207   bool Analysis::isCompatible(const BeamPair& beams) const {
00208     return compatible(beams, requiredBeams());
00209     /// @todo Need to also check internal consistency of the analysis'
00210     /// beam requirements with those of the projections it uses.
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   // Histogramming
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     // Get the bin edges (only read the AIDA file once)
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     // Get the bin edges (only read the AIDA file once)
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     // Get the bin edges (only read the AIDA file once)
00449     _cacheXAxisData();
00450     // Build the axis code
00451     const string axisCode = makeAxisCode(datasetId, xAxisId, yAxisId);
00452     //const map<string, vector<DPSXPoint> > xpoints = getDPSXValsErrs(papername);
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         // vector<string> dirs;
00475         // split(dirs, histoDir(), "/");
00476         // string pathpart;
00477         // foreach (const string& d, dirs) {
00478         //tree().mkdir();
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       // Leaving out factor of binWidth because AIDA's "height" already includes a width factor.
00500       oldintg += histo->binHeight(iBin); // * histo->axis().binWidth(iBin);
00501     }
00502     if (oldintg == 0.0) {
00503       getLog() << Log::WARN << "Histo " << hpath << " has null integral during normalisation" << endl;
00504       return;
00505     }
00506 
00507     // Scale by the normalisation factor.
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       // "Bin height" is a misnomer in the AIDA spec: width is neglected.
00527       // We'd like to do this: y.push_back(histo->binHeight(i) * scale);
00528       y.push_back(histo->binHeight(i)*scale/histo->axis().binWidth(i));
00529 
00530       // "Bin error" is a misnomer in the AIDA spec: width is neglected.
00531       // We'd like to do this: ey.push_back(histo->binError(i) * scale);
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     // Set histo pointer to null - it can no longer be used.
00550     histo = 0;
00551   }
00552 
00553 
00554 }