diff --git a/analyses/pluginALICE/ALICE_2010_I880049.info.orig b/analyses/pluginALICE/ALICE_2010_I880049.info.orig new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/ALICE_2010_I880049.info.orig @@ -0,0 +1,41 @@ +Name: ALICE_2010_I880049 +Year: 2010 +Summary: Centrality dependence of the charged-particle multiplicity density at mid-rapidity in Pb--Pb collisions at $\sqrt{s_{\rm{NN}}} = 2.76$ TeV +Experiment: ALICE +Collider: LHC +SpireID: 032301 +InspireID: 880049 +Status: NEW +Reentrant: True +Authors: +References: + - arXiv:1012.1657 [nucl-ex] +RunInfo: +NumEvents: +Options: + - cent=REF,GEN,IMP,USR +Beams: [1000822080, 1000822080] +# This is _total_ energy of beams, so this becomes 208*2760=574080 +Energies: [574080] +Description: + 'The centrality dependence of the charged-particle multiplicity density at mid-rapidity in Pb-Pb collisions at $\sqrt{s_{NN}} = 2.76$ TeV is presented. + The charged-particle density normalized per participating nucleon pair increases by about a factor 2 from peripheral (70-80%) to central (0-5%) collisions. + The centrality dependence is found to be similar to that observed at lower collision energies. The data are compared with models based on different mechanisms for particle production in nuclear collisions.' +BibKey: Aamodt:2010cz +BibTeX: '@article{Aamodt:2010cz, + author = "Aamodt, Kenneth and others", + title = "{Centrality dependence of the charged-particle + multiplicity density at mid-rapidity in Pb-Pb collisions + at $\sqrt{s_{NN}}=2.76$ TeV}", + collaboration = "ALICE", + journal = "Phys. Rev. Lett.", + volume = "106", + year = "2011", + pages = "032301", + doi = "10.1103/PhysRevLett.106.032301", + eprint = "1012.1657", + archivePrefix = "arXiv", + primaryClass = "nucl-ex", + reportNumber = "CERN-PH-EP-2010-071", + SLACcitation = "%%CITATION = ARXIV:1012.1657;%%" +}' diff --git a/analyses/pluginALICE/ALICE_2012_I930312.cc b/analyses/pluginALICE/ALICE_2012_I930312.cc --- a/analyses/pluginALICE/ALICE_2012_I930312.cc +++ b/analyses/pluginALICE/ALICE_2012_I930312.cc @@ -1,320 +1,320 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Tools/AliceCommon.hh" #include "Rivet/Projections/AliceCommon.hh" namespace Rivet { /// @brief ALICE PbPb at 2.76 TeV azimuthal di-hadron correlations class ALICE_2012_I930312 : public Analysis { public: // Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ALICE_2012_I930312); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Declare centrality projection declareCentrality(ALICE::V0MMultiplicity(), "ALICE_2015_PBPBCentrality", "V0M", "V0M"); // Projection for trigger particles: charged, primary particles // with |eta| < 1.0 and 8 < pT < 15 GeV/c declare(ALICE::PrimaryParticles(Cuts::abseta < 1.0 && Cuts::abscharge > 0 && Cuts::ptIn(8.*GeV, 15.*GeV)), "APRIMTrig"); // pT bins edges vector ptBins = { 3., 4., 6., 8., 10. }; // Projections for associated particles: charged, primary particles // with |eta| < 1.0 and different pT bins for (int ipt = 0; ipt < PT_BINS; ++ipt) { Cut cut = Cuts::abseta < 1.0 && Cuts::abscharge > 0 && Cuts::ptIn(ptBins[ipt]*GeV, ptBins[ipt+1]*GeV); declare(ALICE::PrimaryParticles(cut), "APRIMAssoc" + toString(ipt)); } // Create event strings vector evString = { "pp", "central", "peripheral" }; // Initialize trigger counters and yield histograms string title = "Per trigger particle yield"; string xtitle = "$\\Delta\\eta$ (rad)"; string ytitle = "$1 / N_{trig} {\\rm d}N_{assoc} / {\\rm d}\\Delta\\eta$ (rad$^-1$)"; for (int itype = 0; itype < EVENT_TYPES; ++itype) { _counterTrigger[itype] = bookCounter("counter." + toString(itype)); for (int ipt = 0; ipt < PT_BINS; ++ipt) { string name = "yield." + evString[itype] + ".pt" + toString(ipt); _histYield[itype][ipt] = bookHisto1D(name, 36, -0.5*M_PI, 1.5*M_PI, title, xtitle, ytitle); } } } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); // Trigger particles Particles trigParticles = applyProjection(event,"APRIMTrig").particles(); // Associated particles Particles assocParticles[PT_BINS]; for (int ipt = 0; ipt < PT_BINS; ++ipt) { string pname = "APRIMAssoc" + toString(ipt); assocParticles[ipt] = applyProjection(event,pname).particles(); } // Check type of event. This may not be a perfect way to check for the // type of event as there might be some weird conditions hidden inside. // For example some HepMC versions check if number of hard collisions // is equal to 0 and assign 'false' in that case, which is usually wrong. // This might be changed in the future int ev_type = 0; // pp const HepMC::HeavyIon* hi = event.genEvent()->heavy_ion(); if (hi && hi->is_valid()) { // Prepare centrality projection and value const CentralityProjection& centrProj = apply(event, "V0M"); double centr = centrProj(); // Set the flag for the type of the event if (centr > 0.0 && centr < 5.0) ev_type = 1; // PbPb, central else if (centr > 60.0 && centr < 90.0) ev_type = 2; // PbPb, peripherial else vetoEvent; // PbPb, other, this is not used in the analysis at all } // Fill trigger histogram for a proper event type - _counterTrigger[ev_type]->fill(trigParticles.size()); + _counterTrigger[ev_type]->fill(trigParticles.size()*weight); // Loop over trigger particles for (const Particle& trigParticle : trigParticles) { // For each pt bin for (int ipt = 0; ipt < PT_BINS; ++ipt) { // Loop over associated particles for (const Particle& assocParticle : assocParticles[ipt]) { // If associated and trigger particle are not the same particles. if (!isSame(trigParticle, assocParticle)) { // Test trigger particle. if (trigParticle.pt() > assocParticle.pt()) { // Calculate delta phi in range (-0.5*PI, 1.5*PI). double dPhi = deltaPhi(trigParticle, assocParticle, true); if (dPhi < -0.5 * M_PI) dPhi += 2 * M_PI; // Fill yield histogram for calculated delta phi _histYield[ev_type][ipt]->fill(dPhi, weight); } } } } } } /// Normalise histograms etc., after the run void finalize() { // Check for the reentrant finalize bool pp_available = false, PbPb_available = false; for (int itype = 0; itype < EVENT_TYPES; ++itype) { for (int ipt = 0; ipt < PT_BINS; ++ipt) { if (_histYield[itype][ipt]->numEntries() > 0) itype == 0 ? pp_available = true : PbPb_available = true; } } // Skip postprocessing if pp or PbPb histograms are not available if (!(pp_available && PbPb_available)) return; // Initialize IAA and ICP histograms _histIAA[0] = bookScatter2D(1, 1, 1); _histIAA[1] = bookScatter2D(2, 1, 1); _histIAA[2] = bookScatter2D(5, 1, 1); _histIAA[3] = bookScatter2D(3, 1, 1); _histIAA[4] = bookScatter2D(4, 1, 1); _histIAA[5] = bookScatter2D(6, 1, 1); // Initialize background-subtracted yield histograms for (int itype = 0; itype < EVENT_TYPES; ++itype) { for (int ipt = 0; ipt < PT_BINS; ++ipt) { string newname = _histYield[itype][ipt]->name() + ".nobkg"; string newtitle = _histYield[itype][ipt]->title() + ", background subtracted"; _histYieldNoBkg[itype][ipt] = bookHisto1D(newname, 36, -0.5*M_PI, 1.5*M_PI, newtitle, _histYield[itype][ipt]->annotation("XLabel"), _histYield[itype][ipt]->annotation("YLabel")); } } // Variable for near and away side peak integral calculation double integral[EVENT_TYPES][PT_BINS][2] = { { {0.0} } }; // Variables for background calculation double bkg = 0.0; double bkgErr[EVENT_TYPES][PT_BINS] = { {0.0} }; // Variables for integration error calculation double norm[EVENT_TYPES] = {0.0}; double numEntries[EVENT_TYPES][PT_BINS][2] = { { {0.0} } }; int numBins[EVENT_TYPES][PT_BINS][2] = { { {0} } }; // For each event type for (int itype = 0; itype < EVENT_TYPES; ++itype) { // Get counter CounterPtr counter = _counterTrigger[itype]; // For each pT range for (int ipt = 0; ipt < PT_BINS; ++ipt) { // Get yield histograms Histo1DPtr hYield = _histYield[itype][ipt]; Histo1DPtr hYieldNoBkg = _histYieldNoBkg[itype][ipt]; // Check if histograms are fine if (counter->sumW() == 0 || hYield->numEntries() == 0) { MSG_WARNING("There are no entries in one of the histograms"); continue; } // Scale yield histogram norm[itype] = 1. / counter->sumW(); scale(hYield, norm[itype]); // Calculate background double sum = 0.0; int nbins = 0; for (size_t ibin = 0; ibin < hYield->numBins(); ++ibin) { double xmid = hYield->bin(ibin).xMid(); if (inRange(xmid, -0.5 * M_PI, -0.5 * M_PI + 0.4) || inRange(xmid, 0.5 * M_PI - 0.4, 0.5 * M_PI + 0.4) || inRange(xmid, 1.5 * M_PI - 0.4, 1.5 * M_PI)) { sum += hYield->bin(ibin).sumW(); nbins += 1; } } if (nbins == 0) { MSG_WARNING("Failed to estimate background!"); continue; } bkg = sum / nbins; // Calculate background error sum = 0.0; nbins = 0; for (size_t ibin = 0; ibin < hYield->numBins(); ++ibin) { double xmid = hYield->bin(ibin).xMid(); if (inRange(xmid, 0.5 * M_PI - 0.4, 0.5 * M_PI + 0.4)) { sum += (hYield->bin(ibin).sumW() - bkg) * (hYield->bin(ibin).sumW() - bkg); nbins++; } } if (nbins < 2) { MSG_WARNING("Failed to estimate background error!"); continue; } bkgErr[itype][ipt] = sqrt(sum / (nbins - 1)); // Fill histograms with removed background for (size_t ibin = 0; ibin < hYield->numBins(); ++ibin) { hYieldNoBkg->fillBin(ibin, hYield->bin(ibin).sumW() - bkg); } // Integrate near-side yield size_t lowerBin = hYield->binIndexAt(-0.7 + 0.02); size_t upperBin = hYield->binIndexAt( 0.7 - 0.02) + 1; nbins = upperBin - lowerBin; numBins[itype][ipt][NEAR] = nbins; integral[itype][ipt][NEAR] = hYield->integralRange(lowerBin, upperBin) - nbins * bkg; numEntries[itype][ipt][NEAR] = hYield->integralRange(lowerBin, upperBin) * counter->sumW(); // Integrate away-side yield lowerBin = hYield->binIndexAt(M_PI - 0.7 + 0.02); upperBin = hYield->binIndexAt(M_PI + 0.7 - 0.02) + 1; nbins = upperBin - lowerBin; numBins[itype][ipt][AWAY] = nbins; integral[itype][ipt][AWAY] = hYield->integralRange(lowerBin, upperBin) - nbins * bkg; numEntries[itype][ipt][AWAY] = hYield->integralRange(lowerBin, upperBin) * counter->sumW(); } } // Variables for IAA/ICP plots double yval[2] = { 0.0, 0.0 }; double yerr[2] = { 0.0, 0.0 }; double xval[PT_BINS] = { 3.5, 5.0, 7.0, 9.0 }; double xerr[PT_BINS] = { 0.5, 1.0, 1.0, 1.0 }; int types1[3] = {1, 2, 1}; int types2[3] = {0, 0, 2}; // Fill IAA/ICP plots for near and away side peak for (int ihist = 0; ihist < 3; ++ihist) { int type1 = types1[ihist]; int type2 = types2[ihist]; double norm1 = norm[type1]; double norm2 = norm[type2]; for (int ipt = 0; ipt < PT_BINS; ++ipt) { double bkgErr1 = bkgErr[type1][ipt]; double bkgErr2 = bkgErr[type2][ipt]; for (int ina = 0; ina < 2; ++ina) { double integ1 = integral[type1][ipt][ina]; double integ2 = integral[type2][ipt][ina]; double numEntries1 = numEntries[type1][ipt][ina]; double numEntries2 = numEntries[type2][ipt][ina]; double numBins1 = numBins[type1][ipt][ina]; double numBins2 = numBins[type2][ipt][ina]; yval[ina] = integ1 / integ2; yerr[ina] = norm1 * norm1 * numEntries1 + norm2 * norm2 * numEntries2 * integ1 * integ1 / (integ2 * integ2) + numBins1 * numBins1 * bkgErr1 * bkgErr1 + numBins2 * numBins2 * bkgErr2 * bkgErr2 * integ1 * integ1 / (integ2 * integ2); yerr[ina] = sqrt(yerr[ina])/integ2; } _histIAA[ihist]->addPoint(xval[ipt], yval[NEAR], xerr[ipt], yerr[NEAR]); _histIAA[ihist + 3]->addPoint(xval[ipt], yval[AWAY], xerr[ipt], yerr[AWAY]); } } } //@} private: static const int PT_BINS = 4; static const int EVENT_TYPES = 3; static const int NEAR = 0; static const int AWAY = 1; /// @name Histograms //@{ Histo1DPtr _histYield[EVENT_TYPES][PT_BINS]; Histo1DPtr _histYieldNoBkg[EVENT_TYPES][PT_BINS]; CounterPtr _counterTrigger[EVENT_TYPES]; Scatter2DPtr _histIAA[6]; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2012_I930312); } diff --git a/analyses/pluginALICE/ALICE_2012_I930312.plot b/analyses/pluginALICE/ALICE_2012_I930312.plot --- a/analyses/pluginALICE/ALICE_2012_I930312.plot +++ b/analyses/pluginALICE/ALICE_2012_I930312.plot @@ -1,47 +1,47 @@ # BEGIN PLOT /ALICE_2012_I930312/d01-x01-y01 -Title=$I_{AA}$ near-side +Title=$I_{AA}$ near-side (0-5\%) XLabel=$p_{t,assoc}$ [GeV/c] YLabel=$ I_{AA} $ LogX=0 LogY=0 # END PLOT # BEGIN PLOT /ALICE_2012_I930312/d02-x01-y01 -Title=$I_{AA}$ near-side +Title=$I_{AA}$ near-side (60-90\%) XLabel=$p_{t,assoc}$ [GeV/c] YLabel=$I_{AA}$ LogX=0 LogY=0 # END PLOT # BEGIN PLOT /ALICE_2012_I930312/d03-x01-y01 -Title=$I_{AA}$ away-side +Title=$I_{AA}$ away-side (0-5\%) XLabel=$p_{t,assoc}$ [GeV/c] YLabel=$I_{AA}$ LogX=0 LogY=0 # END PLOT # BEGIN PLOT /ALICE_2012_I930312/d04-x01-y01 -Title=$I_{AA}$ away-side +Title=$I_{AA}$ away-side (60-90\%) XLabel=$p_{t,assoc}$ [GeV/c] YLabel=$I_{AA}$ LogX=0 LogY=0 # END PLOT # BEGIN PLOT /ALICE_2012_I930312/d05-x01-y01 Title=$I_{CP}$ (0-5\% / 60-90\%) near-side XLabel=$p_{t,assoc}$ [GeV/c] YLabel=$I_{CP}$ (0-5\% / 60-90\%) LogX=0 LogY=0 # END PLOT # BEGIN PLOT /ALICE_2012_I930312/d06-x01-y01 Title=$I_{CP}$ (0-5\% / 60-90\%), away-side XLabel=$p_{t,assoc}$ [GeV/c] YLabel=$I_{CP}$ (0-5\% / 60-90\%) LogX=0 LogY=0 # END PLOT diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc --- a/src/Core/AnalysisHandler.cc +++ b/src/Core/AnalysisHandler.cc @@ -1,586 +1,585 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/AnalysisHandler.hh" #include "Rivet/Analysis.hh" #include "Rivet/Tools/ParticleName.hh" #include "Rivet/Tools/BeamConstraint.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/Beam.hh" #include "YODA/IO.h" namespace Rivet { AnalysisHandler::AnalysisHandler(const string& runname) : _runname(runname), _eventcounter("/_EVTCOUNT"), _xs(NAN), _xserr(NAN), _initialised(false), _ignoreBeams(false), _dumpPeriod(0), _dumping(0) { } AnalysisHandler::~AnalysisHandler() { } Log& AnalysisHandler::getLog() const { return Log::getLog("Rivet.Analysis.Handler"); } void AnalysisHandler::init(const GenEvent& ge) { if (_initialised) throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!"); setRunBeams(Rivet::beams(ge)); MSG_DEBUG("Initialising the analysis handler"); _eventcounter.reset(); // Check that analyses are beam-compatible, and remove those that aren't const size_t num_anas_requested = analysisNames().size(); vector anamestodelete; for (const AnaHandle a : analyses()) { if (!_ignoreBeams && !a->isCompatible(beams())) { //MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV"); anamestodelete.push_back(a->name()); } } for (const string& aname : anamestodelete) { MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing"); removeAnalysis(aname); } if (num_anas_requested > 0 && analysisNames().empty()) { cerr << "All analyses were incompatible with the first event's beams\n" << "Exiting, since this probably wasn't intentional!" << endl; exit(1); } // Warn if any analysis' status is not unblemished for (const AnaHandle a : analyses()) { if (toUpper(a->status()) == "PRELIMINARY") { MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!"); } else if (toUpper(a->status()) == "OBSOLETE") { MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!"); } else if (toUpper(a->status()).find("UNVALIDATED") != string::npos) { MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!"); } } // Initialize the remaining analyses for (AnaHandle a : analyses()) { MSG_DEBUG("Initialising analysis: " << a->name()); try { // Allow projection registration in the init phase onwards a->_allowProjReg = true; a->init(); //MSG_DEBUG("Checking consistency of analysis: " << a->name()); //a->checkConsistency(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::init method: " << err.what() << endl; exit(1); } MSG_DEBUG("Done initialising analysis: " << a->name()); } _initialised = true; MSG_DEBUG("Analysis handler initialised"); } void AnalysisHandler::analyze(const GenEvent& ge) { // Call init with event as template if not already initialised if (!_initialised) init(ge); assert(_initialised); // Ensure that beam details match those from the first event (if we're checking beams) if ( !_ignoreBeams ) { const PdgIdPair beams = Rivet::beamIds(ge); const double sqrts = Rivet::sqrtS(ge); if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) { cerr << "Event beams mismatch: " << PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams " << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << endl; exit(1); } } // Create the Rivet event wrapper /// @todo Filter/normalize the event here Event event(ge); // Weights /// @todo Drop this / just report first weight when we support multiweight events _eventcounter.fill(event.weight()); MSG_DEBUG("Event #" << _eventcounter.numEntries() << " weight = " << event.weight()); // Cross-section #ifdef HEPMC_HAS_CROSS_SECTION if (ge.cross_section()) { _xs = ge.cross_section()->cross_section(); _xserr = ge.cross_section()->cross_section_error(); } #endif // Run the analyses for (AnaHandle a : analyses()) { MSG_TRACE("About to run analysis " << a->name()); try { a->analyze(event); } catch (const Error& err) { cerr << "Error in " << a->name() << "::analyze method: " << err.what() << endl; exit(1); } MSG_TRACE("Finished running analysis " << a->name()); } if ( _dumpPeriod > 0 && numEvents()%_dumpPeriod == 0 ) { MSG_INFO("Dumping intermediate results to " << _dumpFile << "."); _dumping = numEvents()/_dumpPeriod; finalize(); writeData(_dumpFile); _dumping = 0; } } void AnalysisHandler::analyze(const GenEvent* ge) { if (ge == nullptr) { MSG_ERROR("AnalysisHandler received null pointer to GenEvent"); //throw Error("AnalysisHandler received null pointer to GenEvent"); } analyze(*ge); } void AnalysisHandler::finalize() { if (!_initialised) return; // First we make copies of all analysis objects. map backupAOs; for (auto ao : getData(false, true, false) ) backupAOs[ao->path()] = AnalysisObjectPtr(ao->newclone()); // Now we run the (re-entrant) finalize() functions for all analyses. MSG_INFO("Finalising analyses"); for (AnaHandle a : analyses()) { a->setCrossSection(_xs); try { if ( !_dumping || a->info().reentrant() ) a->finalize(); else if ( _dumping == 1 ) MSG_INFO("Skipping finalize in periodic dump of " << a->name() << " as it is not declared reentrant."); } catch (const Error& err) { cerr << "Error in " << a->name() << "::finalize method: " << err.what() << endl; exit(1); } } // Now we copy all analysis objects to the list of finalized // ones, and restore the value to their original ones. _finalizedAOs.clear(); for ( auto ao : getData(false, false, false) ) _finalizedAOs.push_back(AnalysisObjectPtr(ao->newclone())); for ( auto ao : getData(false, true, false) ) { // TODO: This should be possible to do in a nicer way, with a flag etc. if (ao->path().find("/FINAL") != std::string::npos) continue; auto aoit = backupAOs.find(ao->path()); if ( aoit == backupAOs.end() ) { AnaHandle ana = analysis(split(ao->path(), "/")[0]); if ( ana ) ana->removeAnalysisObject(ao->path()); } else copyao(aoit->second, ao); } // Print out number of events processed const int nevts = _eventcounter.numEntries(); MSG_INFO("Processed " << nevts << " event" << (nevts != 1 ? "s" : "")); // // Delete analyses // MSG_DEBUG("Deleting analyses"); // _analyses.clear(); // Print out MCnet boilerplate cout << endl; cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << endl; cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << endl; } AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname, std::map pars) { // Make an option handle. std::string parHandle = ""; for (map::iterator par = pars.begin(); par != pars.end(); ++par) { parHandle +=":"; parHandle += par->first + "=" + par->second; } return addAnalysis(analysisname + parHandle); } AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) { // Check for a duplicate analysis /// @todo Might we want to be able to run an analysis twice, with different params? /// Requires avoiding histo tree clashes, i.e. storing the histos on the analysis objects. string ananame = analysisname; vector anaopt = split(analysisname, ":"); if ( anaopt.size() > 1 ) ananame = anaopt[0]; AnaHandle analysis( AnalysisLoader::getAnalysis(ananame) ); if (analysis.get() != 0) { // < Check for null analysis. MSG_DEBUG("Adding analysis '" << analysisname << "'"); map opts; for ( int i = 1, N = anaopt.size(); i < N; ++i ) { vector opt = split(anaopt[i], "="); if ( opt.size() != 2 ) { MSG_WARNING("Error in option specification. Skipping analysis " << analysisname); return *this; } if ( !analysis->info().validOption(opt[0], opt[1]) ) { MSG_WARNING("Cannot set option '" << opt[0] << "' to '" << opt[1] << "'. Skipping analysis " << analysisname); return *this; } opts[opt[0]] = opt[1]; } for ( auto opt: opts) { analysis->_options[opt.first] = opt.second; analysis->_optstring += ":" + opt.first + "=" + opt.second; } for (const AnaHandle& a : analyses()) { if (a->name() == analysis->name() ) { MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate"); return *this; } } analysis->_analysishandler = this; _analyses[analysisname] = analysis; } else { MSG_WARNING("Analysis '" << analysisname << "' not found."); } // MSG_WARNING(_analyses.size()); // for (const AnaHandle& a : _analyses) MSG_WARNING(a->name()); return *this; } AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) { // std::shared_ptr toremove; // for (const AnaHandle a : _analyses) { // if (a->name() == analysisname) { // toremove = a; // break; // } // } // if (toremove.get() != 0) { MSG_DEBUG("Removing analysis '" << analysisname << "'"); if (_analyses.find(analysisname) != _analyses.end()) _analyses.erase(analysisname); // } return *this; } ///////////////////////////// void AnalysisHandler::addData(const std::vector& aos) { for (const AnalysisObjectPtr ao : aos) { string path = ao->path(); if ( path.substr(0, 5) != "/RAW/" ) { _orphanedPreloads.push_back(ao); continue; } path = path.substr(4); ao->setPath(path); if (path.size() > 1) { // path > "/" try { const string ananame = split(path, "/")[0]; AnaHandle a = analysis(ananame); a->addAnalysisObject(ao); /// @todo Need to statistically merge... } catch (const Error& e) { MSG_TRACE("Adding analysis object " << path << " to the list of orphans."); _orphanedPreloads.push_back(ao); } } } } void AnalysisHandler::stripOptions(AnalysisObjectPtr ao, const vector & delopts) const { string path = ao->path(); string ananame = split(path, "/")[0]; vector anaopts = split(ananame, ":"); for ( int i = 1, N = anaopts.size(); i < N; ++i ) for ( auto opt : delopts ) if ( opt == "*" || anaopts[i].find(opt + "=") == 0 ) path.replace(path.find(":" + anaopts[i]), (":" + anaopts[i]).length(), ""); ao->setPath(path); } void AnalysisHandler:: mergeYodas(const vector & aofiles, const vector & delopts, bool equiv) { vector< vector > aosv; vector xsecs; vector xsecerrs; vector sows; set ananames; _eventcounter.reset(); // First scan all files and extract analysis objects and add the // corresponding anayses.. for ( auto file : aofiles ) { Scatter1DPtr xsec; CounterPtr sow; // For each file make sure that cross section and sum-of-weights // objects are present and stor all RAW ones in a vector; vector aos; try { /// @todo Use new YODA SFINAE to fill the smart ptr vector directly vector aos_raw; YODA::read(file, aos_raw); for (AnalysisObject* aor : aos_raw) { AnalysisObjectPtr ao = AnalysisObjectPtr(aor); if ( ao->path().substr(0, 5) != "/RAW/" ) continue; ao->setPath(ao->path().substr(4)); if ( ao->path() == "/_XSEC" ) xsec = dynamic_pointer_cast(ao); else if ( ao->path() == "/_EVTCOUNT" ) sow = dynamic_pointer_cast(ao); else { stripOptions(ao, delopts); string ananame = split(ao->path(), "/")[0]; if ( ananames.insert(ananame).second ) addAnalysis(ananame); aos.push_back(ao); } } if ( !xsec || !sow ) { MSG_ERROR( "Error in AnalysisHandler::mergeYodas: The file " << file << " did not contain weights and cross section info."); exit(1); } xsecs.push_back(xsec->point(0).x()); - sows.push_back(sow); xsecerrs.push_back(sqr(xsec->point(0).xErrAvg())); _eventcounter += *sow; sows.push_back(sow); aosv.push_back(aos); } catch (...) { //< YODA::ReadError& throw UserError("Unexpected error in reading file: " + file); } } // Now calculate the scale to be applied for all bins in a file // and get the common cross section and sum of weights. _xs = _xserr = 0.0; for ( int i = 0, N = sows.size(); i < N; ++i ) { double effnent = sows[i]->effNumEntries(); _xs += (equiv? effnent: 1.0)*xsecs[i]; _xserr += (equiv? sqr(effnent): 1.0)*xsecerrs[i]; } vector scales(sows.size(), 1.0); if ( equiv ) { _xs /= _eventcounter.effNumEntries(); _xserr = sqrt(_xserr)/_eventcounter.effNumEntries(); } else { _xserr = sqrt(_xserr); for ( int i = 0, N = sows.size(); i < N; ++i ) scales[i] = (_eventcounter.sumW()/sows[i]->sumW())*(xsecs[i]/_xs); } // Initialize the analyses allowing them to book analysis objects. for (AnaHandle a : analyses()) { MSG_DEBUG("Initialising analysis: " << a->name()); if ( !a->info().reentrant() ) MSG_WARNING("Analysis " << a->name() << " has not been validated to have " << "a reentrant finalize method. The result is unpredictable."); try { // Allow projection registration in the init phase onwards a->_allowProjReg = true; cerr << "sqrtS " << sqrtS() << endl; a->init(); //MSG_DEBUG("Checking consistency of analysis: " << a->name()); //a->checkConsistency(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::init method: " << err.what() << endl; exit(1); } MSG_DEBUG("Done initialising analysis: " << a->name()); } _initialised = true; // Get a list of all anaysis objects to handle. map current; for ( auto ao : getData(false, true, false) ) current[ao->path()] = ao; // Go through all objects to be merged and add them to current // after appropriate scaling. for ( int i = 0, N = aosv.size(); i < N; ++i) for ( auto ao : aosv[i] ) { if ( ao->path() == "/_XSEC" || ao->path() == "_EVTCOUNT" ) continue; auto aoit = current.find(ao->path()); if ( aoit == current.end() ) { MSG_WARNING("" << ao->path() << " was not properly booked."); continue; } if ( !addaos(aoit->second, ao, scales[i]) ) MSG_WARNING("Cannot merge objects with path " << ao->path() <<" of type " << ao->annotation("Type") ); } // Now we can simply finalize() the analysis, leaving the // controlling program to write it out some yoda-file. finalize(); } void AnalysisHandler::readData(const string& filename) { vector aos; try { /// @todo Use new YODA SFINAE to fill the smart ptr vector directly vector aos_raw; YODA::read(filename, aos_raw); for (AnalysisObject* aor : aos_raw) aos.push_back(AnalysisObjectPtr(aor)); } catch (...) { //< YODA::ReadError& throw UserError("Unexpected error in reading file: " + filename); } if (!aos.empty()) addData(aos); } vector AnalysisHandler:: getData(bool includeorphans, bool includetmps, bool usefinalized) const { vector rtn; // Event counter rtn.push_back( make_shared(_eventcounter) ); // Cross-section + err as scatter YODA::Scatter1D::Points pts; pts.insert(YODA::Point1D(_xs, _xserr)); rtn.push_back( make_shared(pts, "/_XSEC") ); // Analysis histograms vector aos; if (usefinalized) aos = _finalizedAOs; else { for (const AnaHandle a : analyses()) { // MSG_WARNING(a->name() << " " << aos.size()); for (const AnalysisObjectPtr ao : a->analysisObjects()) { aos.push_back(ao); } } } for (const AnalysisObjectPtr ao : aos) { // Exclude paths from final write-out if they contain a "TMP" layer (i.e. matching "/TMP/") /// @todo This needs to be much more nuanced for re-entrant histogramming if ( !includetmps && ao->path().find("/TMP/" ) != string::npos) continue; rtn.push_back(ao); } // Sort histograms alphanumerically by path before write-out sort(rtn.begin(), rtn.end(), [](AnalysisObjectPtr a, AnalysisObjectPtr b) {return a->path() < b->path();}); if ( includeorphans ) rtn.insert(rtn.end(), _orphanedPreloads.begin(), _orphanedPreloads.end()); return rtn; } void AnalysisHandler::writeData(const string& filename) const { vector out = _finalizedAOs; set finalana; for ( auto ao : out) finalana.insert(ao->path()); out.reserve(2*out.size()); vector aos = getData(false, true, false); if ( _dumping ) { for ( auto ao : aos ) { if ( finalana.find(ao->path()) == finalana.end() ) out.push_back(AnalysisObjectPtr(ao->newclone())); } } for ( auto ao : aos ) { ao = AnalysisObjectPtr(ao->newclone()); ao->setPath("/RAW" + ao->path()); out.push_back(ao); } try { YODA::write(filename, out.begin(), out.end()); } catch (...) { //< YODA::WriteError& throw UserError("Unexpected error in writing file: " + filename); } } std::vector AnalysisHandler::analysisNames() const { std::vector rtn; for (AnaHandle a : analyses()) { rtn.push_back(a->name()); } return rtn; } AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { //MSG_DEBUG("Adding analysis '" << aname << "'"); addAnalysis(aname); } return *this; } AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { removeAnalysis(aname); } return *this; } bool AnalysisHandler::needCrossSection() const { bool rtn = false; for (const AnaHandle a : analyses()) { if (!rtn) rtn = a->needsCrossSection(); if (rtn) break; } return rtn; } AnalysisHandler& AnalysisHandler::setCrossSection(double xs, double xserr) { _xs = xs; _xserr = xserr; return *this; } bool AnalysisHandler::hasCrossSection() const { return (!std::isnan(crossSection())); } AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) { analysis->_analysishandler = this; // _analyses.insert(AnaHandle(analysis)); _analyses[analysis->name()] = AnaHandle(analysis); return *this; } PdgIdPair AnalysisHandler::beamIds() const { return Rivet::beamIds(beams()); } double AnalysisHandler::sqrtS() const { return Rivet::sqrtS(beams()); } void AnalysisHandler::setIgnoreBeams(bool ignore) { _ignoreBeams=ignore; } }