diff --git a/analyses/pluginALICE/ALICE_2012_I1127497.cc b/analyses/pluginALICE/ALICE_2012_I1127497.cc --- a/analyses/pluginALICE/ALICE_2012_I1127497.cc +++ b/analyses/pluginALICE/ALICE_2012_I1127497.cc @@ -1,201 +1,221 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Tools/Cuts.hh" #include "Rivet/Projections/SingleValueProjection.hh" #include "Rivet/Tools/AliceCommon.hh" #include "Rivet/Projections/AliceCommon.hh" #include "Rivet/Projections/HepMCHeavyIon.hh" #include #define _USE_MATH_DEFINES #include namespace Rivet { /// @brief ALICE PbPb at 2.76 TeV R_AA analysis. class ALICE_2012_I1127497 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ALICE_2012_I1127497); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Access the HepMC heavy ion info declare(HepMCHeavyIon(), "HepMC"); // Declare centrality projection declareCentrality(ALICE::V0MMultiplicity(), "ALICE_2015_PBPBCentrality", "V0M", "V0M"); // Charged, primary particles with |eta| < 0.5 and pT > 150 MeV declare(ALICE::PrimaryParticles(Cuts::abseta < 0.5 && Cuts::pT > 150*MeV && Cuts::abscharge > 0), "APRIM"); // Loop over all histograms for (size_t ihist = 0; ihist < NHISTOS; ++ihist) { // Initialize PbPb objects _histNch[PBPB][ihist] = bookHisto1D(ihist+1, 1, 1); std::string nameCounterPbPb = "counter.pbpb." + std::to_string(ihist); _counterSOW[PBPB][ihist] = bookCounter(nameCounterPbPb, "Sum of weights counter for PbPb"); std::string nameCounterNcoll = "counter.ncoll." + std::to_string(ihist); _counterNcoll[ihist] = bookCounter(nameCounterNcoll, "Ncoll counter for PbPb"); // Initialize pp objects. In principle, only one pp histogram would be // needed since centrality does not make any difference here. However, // in some cases in this analysis the binning differ from each other, // so this is easy-to-implement way to account for that. std::string namePP = _histNch[PBPB][ihist]->name() + "-pp"; // The binning is taken from the reference data _histNch[PP][ihist] = bookHisto1D(namePP, refData(ihist+1, 1, 1)); std::string nameCounterpp = "counter.pp." + std::to_string(ihist); _counterSOW[PP][ihist] = bookCounter(nameCounterpp, "Sum of weights counter for pp"); } // Centrality regions keeping boundaries for a certain region. // Note, that some regions overlap with other regions. _centrRegions.clear(); _centrRegions = {{0., 5.}, {5., 10.}, {10., 20.}, {20., 30.}, {30., 40.}, {40., 50.}, {50., 60.}, {60., 70.}, {70., 80.}, {0., 10.}, {0., 20.}, {20., 40.}, {40., 60.}, {40., 80.}, {60., 80.}}; + // Find out the beam type, also specified from option. + string beamOpt = getOption("beam","NONE"); + if (beamOpt != "NONE") { + MSG_WARNING("You are using a specified beam type, instead of using what" + "is provided by the generator. " + "Only do this if you are completely sure what you are doing."); + if (beamOpt=="PP") isHI = false; + else if (beamOpt=="HI") isHI = true; + else { + MSG_ERROR("Beam error (option)!"); + return; + } + } + else { + const ParticlePair& beam = beams(); + if (beam.first.pid() == 2212 && beam.second.pid() == 2212) isHI = false; + else if (beam.first.pid() == 1000822080 && beam.second.pid() == 1000822080) + isHI = true; + else { + MSG_ERROR("Beam error (found)!"); + return; + } + } } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); // Charged, primary particles with at least pT = 150 MeV // in eta range of |eta| < 0.5 Particles chargedParticles = applyProjection(event,"APRIM").particlesByPt(); - // 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 - if (event.genEvent()->heavy_ion()) { + // Check type of event. + if ( isHI ) { const HepMCHeavyIon & hi = apply(event, "HepMC"); // Prepare centrality projection and value const CentralityProjection& centrProj = apply(event, "V0M"); double centr = centrProj(); // Veto event for too large centralities since those are not used // in the analysis at all if ((centr < 0.) || (centr > 80.)) vetoEvent; // Fill PbPb histograms and add weights based on centrality value for (size_t ihist = 0; ihist < NHISTOS; ++ihist) { if (inRange(centr, _centrRegions[ihist].first, _centrRegions[ihist].second)) { _counterSOW[PBPB][ihist]->fill(weight); _counterNcoll[ihist]->fill(hi.Ncoll(), weight); foreach (const Particle& p, chargedParticles) { float pT = p.pT()/GeV; if (pT < 50.) { double pTAtBinCenter = _histNch[PBPB][ihist]->binAt(pT).xMid(); _histNch[PBPB][ihist]->fill(pT, weight/pTAtBinCenter); } } } } } else { // Fill all pp histograms and add weights for (size_t ihist = 0; ihist < NHISTOS; ++ihist) { _counterSOW[PP][ihist]->fill(weight); foreach (const Particle& p, chargedParticles) { float pT = p.pT()/GeV; if (pT < 50.) { double pTAtBinCenter = _histNch[PP][ihist]->binAt(pT).xMid(); _histNch[PP][ihist]->fill(pT, weight/pTAtBinCenter); } } } } } /// Normalise histograms etc., after the run void finalize() { // Right scaling of the histograms with their individual weights. for (size_t itype = 0; itype < EVENT_TYPES; ++itype ) { for (size_t ihist = 0; ihist < NHISTOS; ++ihist) { if (_counterSOW[itype][ihist]->sumW() > 0.) { scale(_histNch[itype][ihist], (1. / _counterSOW[itype][ihist]->sumW() / 2. / M_PI)); } } } // Postprocessing of the histograms for (size_t ihist = 0; ihist < NHISTOS; ++ihist) { // If there are entires in histograms for both beam types if (_histNch[PP][ihist]->numEntries() > 0 && _histNch[PBPB][ihist]->numEntries() > 0) { // Initialize and fill R_AA histograms _histRAA[ihist] = bookScatter2D(ihist+16, 1, 1); divide(_histNch[PBPB][ihist], _histNch[PP][ihist], _histRAA[ihist]); // Scale by Ncoll. Unfortunately some generators does not provide // Ncoll value (eg. JEWEL), so the following scaling will be done // only if there are entries in the counters double ncoll = _counterNcoll[ihist]->sumW(); double sow = _counterSOW[PBPB][ihist]->sumW(); if (ncoll > 1e-6 && sow > 1e-6) _histRAA[ihist]->scaleY(1. / (ncoll / sow)); } } } //@} private: + bool isHI; static const int NHISTOS = 15; static const int EVENT_TYPES = 2; static const int PP = 0; static const int PBPB = 1; /// @name Histograms //@{ Histo1DPtr _histNch[EVENT_TYPES][NHISTOS]; CounterPtr _counterSOW[EVENT_TYPES][NHISTOS]; CounterPtr _counterNcoll[NHISTOS]; Scatter2DPtr _histRAA[NHISTOS]; //@} std::vector> _centrRegions; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2012_I1127497); } diff --git a/analyses/pluginALICE/ALICE_2012_I1127497.info b/analyses/pluginALICE/ALICE_2012_I1127497.info --- a/analyses/pluginALICE/ALICE_2012_I1127497.info +++ b/analyses/pluginALICE/ALICE_2012_I1127497.info @@ -1,41 +1,42 @@ Name: ALICE_2012_I1127497 Year: 2012 Summary: Centrality dependence of charged particle production at large transverse momentum in Pb-Pb collisions at $\sqrt{s_{\rm{NN}}} = 2.76$ TeV Experiment: ALICE Collider: LHC SpireID: InspireID: 1127497 Status: VALIDATED REENTRANT Authors: - Przemyslaw Karczmarczyk - Jan Fiete Grosse-Oetringhaus - Jochen Klein References: - arXiv:1208.2711 [hep-ex] RunInfo: NumEvents: Options: - cent=REF,GEN,IMP,USR + - beam=HI,PP Beams: [[p, p], [1000822080, 1000822080]] # This is _total_ energy of beams, so this becomes 208*2760=574080 Energies: [2760, 574080] Description: 'The inclusive transverse momentum ($p_T$) distributions of primary charged particles are measured in the pseudo-rapidity range $|\eta| < 0.8$ as a function of event centrality in Pb--Pb collisions at $\sqrt{s_{nn}} = 2.76$ TeV with ALICE at the LHC. The data are presented in the $p_T$ range $0.15 30$ GeV/c. In peripheral collisions (70-80%), the suppression is weaker with $R_{AA} \approx 0.7$ almost independently of $p_T$. The measured nuclear modification factors are compared to other measurements and model calculations.' BibKey: Abelev:2012hxa BibTeX: '@article{Abelev:2012hxa, author = "Abelev, Betty and others", title = "{Centrality Dependence of Charged Particle Production at Large Transverse Momentum in Pb--Pb Collisions at $\sqrt{s_{\rm{NN}}} = 2.76$ TeV}", collaboration = "ALICE", journal = "Phys. Lett.", volume = "B720", year = "2013", pages = "52-62", doi = "10.1016/j.physletb.2013.01.051", eprint = "1208.2711", archivePrefix = "arXiv", primaryClass = "hep-ex", reportNumber = "CERN-PH-EP-2012-233", SLACcitation = "%%CITATION = ARXIV:1208.2711;%%" }' 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,319 +1,344 @@ // -*- 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); } } + // Find out the beam type, also specified from option. + string beamOpt = getOption("beam","NONE"); + if (beamOpt != "NONE") { + MSG_WARNING("You are using a specified beam type, instead of using what" + "is provided by the generator. " + "Only do this if you are completely sure what you are doing."); + if (beamOpt=="PP") isHI = false; + else if (beamOpt=="HI") isHI = true; + else { + MSG_ERROR("Beam error (option)!"); + return; + } + } + else { + const ParticlePair& beam = beams(); + if (beam.first.pid() == 2212 && beam.second.pid() == 2212) isHI = false; + else if (beam.first.pid() == 1000822080 && beam.second.pid() == 1000822080) + isHI = true; + else { + MSG_ERROR("Beam error (found)!"); + return; + } + } + } /// 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 - if (event.genEvent()->heavy_ion()) { + if ( isHI ) { // 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()*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: + bool isHI; 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.info b/analyses/pluginALICE/ALICE_2012_I930312.info --- a/analyses/pluginALICE/ALICE_2012_I930312.info +++ b/analyses/pluginALICE/ALICE_2012_I930312.info @@ -1,41 +1,42 @@ Name: ALICE_2012_I930312 Year: 2012 Summary: Particle-yield modification in jet-like azimuthal di-hadron correlations in Pb-Pb collisions at $\sqrt{s_\rm{NN}} = 2.76$ TeV Experiment: ALICE Collider: LHC SpireID: InspireID: 930312 Status: VALIDATED REENTRANT Authors: - Przemyslaw Karczmarczyk - Jan Fiete Grosse-Oetringhaus - Jochen Klein References: - arXiv:1110.0121 [nucl-ex] RunInfo: NumEvents: Options: - cent=REF,GEN,IMP,USR + - beam=HI,PP Beams: [[p, p], [1000822080, 1000822080]] # This is _total_ energy of beams, so this becomes 208*2760=574080 Energies: [2760, 574080] Description: 'The yield of charged particles associated with high-pT trigger particles ($8 < p_{\perp} < 15$ GeV/c) is measured with the ALICE detector in Pb-Pb collisions at $\sqrt{s_{NN}} = 2.76$ TeV relative to proton-proton collisions at the same energy. The conditional per-trigger yields are extracted from the narrow jet-like correlation peaks in azimuthal di-hadron correlations. In the 5\% most central collisions, we observe that the yield of associated charged particles with transverse momenta $p_\perp > 3$ GeV/c on the away-side drops to about 60\% of that observed in pp collisions, while on the near-side a moderate enhancement of 20-30\% is found.' BibKey: Aamodt:2011vg BibTeX: '@article{Aamodt:2011vg, author = "Aamodt, K. and others", title = "{Particle-yield modification in jet-like azimuthal di-hadron correlations in Pb-Pb collisions at $\sqrt{s_{NN}} = 2.76$ TeV}", collaboration = "ALICE", journal = "Phys. Rev. Lett.", volume = "108", year = "2012", pages = "092301", doi = "10.1103/PhysRevLett.108.092301", eprint = "1110.0121", archivePrefix = "arXiv", primaryClass = "nucl-ex", reportNumber = "CERN-PH-EP-2011-161", SLACcitation = "%%CITATION = ARXIV:1110.0121;%%" }' diff --git a/bin/rivet-cmphistos b/bin/rivet-cmphistos --- a/bin/rivet-cmphistos +++ b/bin/rivet-cmphistos @@ -1,495 +1,498 @@ #! /usr/bin/env python """\ Generate histogram comparison plots USAGE: %(prog)s [options] yodafile1[:'PlotOption1=Value':'PlotOption2=Value':...] [path/to/yodafile2 ...] [PLOT:Key1=Val1:...] where the plot options are described in the make-plots manual in the HISTOGRAM section. ENVIRONMENT: * RIVET_ANALYSIS_PATH: list of paths to be searched for plugin analysis libraries at runtime * RIVET_DATA_PATH: list of paths to be searched for data files """ from __future__ import print_function import rivet, yoda, sys, os rivet.util.check_python_version() rivet.util.set_process_name(os.path.basename(__file__)) class Plot(dict): "A tiny Plot object to help writing out the head in the .dat file" def __repr__(self): return "# BEGIN PLOT\n" + "\n".join("%s=%s" % (k,v) for k,v in self.items()) + "\n# END PLOT\n\n" def sanitiseString(s): #s = s.replace('_','\\_') #s = s.replace('^','\\^{}') #s = s.replace('$','\\$') s = s.replace('#','\\#') s = s.replace('%','\\%') return s def getHistos(filelist): """Loop over all input files. Only use the first occurrence of any REF-histogram and the first occurrence in each MC file for every MC-histogram.""" refhistos, mchistos = {}, {} for infile in filelist: mchistos.setdefault(infile, {}) analysisobjects = yoda.read(infile, patterns=args.PATHPATTERNS, unpatterns=args.PATHUNPATTERNS) #print(analysisobjects) for path, ao in analysisobjects.items(): ## We can't plot non-histograms yet # TODO: support counter plotting with a faked x (or y) position and forced plot width/height if ao.type not in ("Counter", "Histo1D", "Histo2D", "Profile1D", "Profile2D", "Scatter1D", "Scatter2D", "Scatter3D"): continue ## Make a path object and ensure the path is in standard form try: aop = rivet.AOPath(path) except Exception as e: #print(e) print("Found analysis object with non-standard path structure:", path, "... skipping") continue ## We don't plot data objects with path components hidden by an underscore prefix if aop.istmp() or aop.israw(): continue ## Add it to the ref or mc paths, if this path isn't already known basepath = aop.basepath(keepref=False) if aop.isref() and basepath not in refhistos: ao.path = aop.varpath(keepref=False, defaultvarid=0) refhistos[basepath] = ao else: #if basepath not in mchistos[infile]: mchistos[infile].setdefault(basepath, {})[aop.varid(0)] = ao return refhistos, mchistos def getRivetRefData(anas=None): "Find all Rivet reference data files" refhistos = {} rivet_data_dirs = rivet.getAnalysisRefPaths() dirlist = [] for d in rivet_data_dirs: if anas is None: import glob dirlist.append(glob.glob(os.path.join(d, '*.yoda'))) else: dirlist.append([os.path.join(d, a+'.yoda') for a in anas]) for filelist in dirlist: # TODO: delegate to getHistos? for infile in filelist: analysisobjects = yoda.read(infile, patterns=args.PATHPATTERNS, unpatterns=args.PATHUNPATTERNS) for path, ao in analysisobjects.items(): aop = rivet.AOPath(ao.path) if aop.isref(): ao.path = aop.basepath(keepref=False) refhistos[ao.path] = ao return refhistos def parseArgs(args): """Look at the argument list and split it at colons, in order to separate the file names from the plotting options. Store the file names and file specific plotting options.""" filelist = [] plotoptions = {} for a in args: asplit = a.split(':') path = asplit[0] filelist.append(path) plotoptions[path] = [] has_title = False for i in range(1, len(asplit)): ## Add 'Title' if there is no = sign before math mode if '=' not in asplit[i] or ('$' in asplit[i] and asplit[i].index('$') < asplit[i].index('=')): asplit[i] = 'Title=%s' % asplit[i] if asplit[i].startswith('Title='): has_title = True plotoptions[path].append(asplit[i]) if path != "PLOT" and not has_title: plotoptions[path].append('Title=%s' % sanitiseString(os.path.basename( os.path.splitext(path)[0] )) ) return filelist, plotoptions def setStyle(ao, istyle, variation=False): """Set default plot styles (color and line width) colors borrowed from Google Ngrams""" # LINECOLORS = ['{[HTML]{EE3311}}', # red (Google uses 'DC3912') # '{[HTML]{3366FF}}', # blue # '{[HTML]{109618}}', # green # '{[HTML]{FF9900}}', # orange # '{[HTML]{990099}}'] # lilac LINECOLORS = ['red', 'blue', 'green', 'orange', 'lilac'] LINESTYLES = ['solid', 'dashed', 'dashdotted', 'dotted'] if args.STYLE == 'talk': ao.setAnnotation('LineWidth', '1pt') if args.STYLE == 'bw': LINECOLORS = ['black!90', 'black!50', 'black!30'] jc = istyle % len(LINECOLORS) c = LINECOLORS[jc] js = (istyle // len(LINECOLORS)) % len(LINESTYLES) s = LINESTYLES[js] ## If plotting a variation (i.e. band), fade the colour if variation: c += "!30" ao.setAnnotation('LineStyle', '%s' % s) ao.setAnnotation('LineColor', '%s' % c) def setOptions(ao, options): "Set arbitrary annotations" for opt in options: key, val = opt.split('=', 1) ao.setAnnotation(key, val) # TODO: move to rivet.utils def mkoutdir(outdir): "Function to make output directories" if not os.path.exists(outdir): try: os.makedirs(outdir) except: msg = "Can't make output directory '%s'" % outdir raise Exception(msg) if not os.access(outdir, os.W_OK): msg = "Can't write to output directory '%s'" % outdir raise Exception(msg) def mkOutput(hpath, aos, plot=None, special=None): """ Make the .dat file string. We can't use "yoda.writeFLAT(anaobjects, 'foobar.dat')" because the PLOT and SPECIAL blocks don't have a corresponding analysis object. """ output = '' if plot is not None: output += str(plot) if special is not None: output += "\n" output += "# BEGIN SPECIAL %s\n" % hpath output += special output += "# END SPECIAL\n\n" from io import StringIO sio = StringIO() yoda.writeFLAT(aos, sio) output += sio.getvalue() return output def writeOutput(output, h): "Choose output file name and dir" if args.HIER_OUTPUT: hparts = h.strip("/").split("/", 1) ana = "_".join(hparts[:-1]) if len(hparts) > 1 else "ANALYSIS" outdir = os.path.join(args.OUTDIR, ana) outfile = '%s.dat' % hparts[-1].replace("/", "_") else: hparts = h.strip("/").split("/") outdir = args.OUTDIR outfile = '%s.dat' % "_".join(hparts) mkoutdir(outdir) outfilepath = os.path.join(outdir, outfile) f = open(outfilepath, 'w') f.write(output) f.close() #-------------------------------------------------------------------------------------------- if __name__ == '__main__': ## Command line parsing import argparse parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("ARGS", nargs="*") parser.add_argument('-o', '--outdir', dest='OUTDIR', default='.', help='write data files into this directory') parser.add_argument("--hier-out", action="store_true", dest="HIER_OUTPUT", default=False, help="write output dat files into a directory hierarchy which matches the analysis paths") parser.add_argument('--plotinfodir', dest='PLOTINFODIRS', action='append', default=['.'], help='directory which may contain plot header information (in addition ' 'to standard Rivet search paths)') parser.add_argument("--no-rivet-refs", dest="RIVETREFS", action="store_false", default=True, help="don't use Rivet reference data files") # parser.add_argument("--refid", dest="REF_ID", # default="REF", help="ID of reference data set (file path for non-REF data)") parser.add_argument("--reftitle", dest="REFTITLE", default='Data', help="Reference data legend entry") parser.add_argument("--pwd", dest="PATH_PWD", action="store_true", default=False, help="append the current directory (pwd) to the analysis/data search paths (cf. $RIVET_ANALYSIS/DATA_PATH)") parser.add_argument("-v", "--verbose", dest="VERBOSE", action="store_true", default=False, help="produce debug output to the terminal") stygroup = parser.add_argument_group("Plot style") stygroup.add_argument("--linear", action="store_true", dest="LINEAR", default=False, help="plot with linear scale") stygroup.add_argument("--errs", "--mcerrs", "--mc-errs", action="store_true", dest="MC_ERRS", default=False, help="show vertical error bars on the MC lines") stygroup.add_argument("--no-ratio", action="store_false", dest="RATIO", default=True, help="disable the ratio plot") stygroup.add_argument("--rel-ratio", action="store_true", dest="RATIO_DEVIATION", default=False, help="show the ratio plots scaled to the ref error") stygroup.add_argument("--no-plottitle", action="store_true", dest="NOPLOTTITLE", default=False, help="don't show the plot title on the plot " "(useful when the plot description should only be given in a caption)") stygroup.add_argument("--style", dest="STYLE", default="default", help="change plotting style: default|bw|talk") stygroup.add_argument("-c", "--config", dest="CONFIGFILES", action="append", default=["~/.make-plots"], help="additional plot config file(s). Settings will be included in the output configuration.") stygroup.add_argument("--remove-options", help="remove options label from legend", dest="REMOVE_OPTIONS", action="store_true", default=False) selgroup = parser.add_argument_group("Selective plotting") # selgroup.add_argument("--show-single", dest="SHOW_SINGLE", choices=("no", "ref", "mc", "all"), # default="mc", help="control if a plot file is made if there is only one dataset to be plotted " # "[default=%(default)s]. If the value is 'no', single plots are always skipped, for 'ref' and 'mc', " # "the plot will be written only if the single plot is a reference plot or an MC " # "plot respectively, and 'all' will always create single plot files.\n The 'ref' and 'all' values " # "should be used with great care, as they will also write out plot files for all reference " # "histograms without MC traces: combined with the -R/--rivet-refs flag, this is a great way to " # "write out several thousand irrelevant reference data histograms!") # selgroup.add_argument("--show-mc-only", "--all", action="store_true", dest="SHOW_IF_MC_ONLY", # default=False, help="make a plot file even if there is only one dataset to be plotted and " # "it is an MC one. Deprecated and will be removed: use --show-single instead, which overrides this.") # # selgroup.add_argument("-l", "--histogram-list", dest="HISTOGRAMLIST", # # default=None, help="specify a file containing a list of histograms to plot, in the format " # # "/ANALYSIS_ID/histoname, one per line, e.g. '/DELPHI_1996_S3430090/d01-x01-y01'.") selgroup.add_argument("-m", "--match", action="append", help="only write out histograms whose $path/$name string matches these regexes. The argument " "may also be a text file.", dest="PATHPATTERNS") selgroup.add_argument("-M", "--unmatch", action="append", help="exclude histograms whose $path/$name string matches these regexes", dest="PATHUNPATTERNS") args = parser.parse_args() ## Add pwd to search paths if args.PATH_PWD: rivet.addAnalysisLibPath(os.path.abspath(".")) rivet.addAnalysisDataPath(os.path.abspath(".")) ## Split the input file names and the associated plotting options ## given on the command line into two separate lists filelist, plotoptions = parseArgs(args.ARGS) ## Remove the PLOT dummy file from the file list if "PLOT" in filelist: filelist.remove("PLOT") ## Check that the files exist for f in filelist: if not os.access(f, os.R_OK): print("Error: cannot read from %s" % f) sys.exit(1) ## Read the .plot files plotdirs = args.PLOTINFODIRS + [os.path.abspath(os.path.dirname(f)) for f in filelist] plotparser = rivet.mkStdPlotParser(plotdirs, args.CONFIGFILES) ## Create a list of all histograms to be plotted, and identify if they are 2D histos (which need special plotting) try: refhistos, mchistos = getHistos(filelist) except IOError as e: print("File reading error: ", e.strerror) exit(1) hpaths, h2ds = [], [] for aos in mchistos.values(): for p in aos.keys(): ps = rivet.stripOptions(p) if ps and ps not in hpaths: hpaths.append(ps) firstaop = aos[p][sorted(aos[p].keys())[0]] # TODO: Would be nicer to test via isHisto and dim or similar, or yoda.Scatter/Histo/Profile base classes if type(firstaop) in (yoda.Histo2D, yoda.Profile2D) and ps not in h2ds: h2ds.append(ps) + # Unique list of analyses + anas = list(set([x.split("/")[1] for x in hpaths])) + ## Take reference data from the Rivet search paths, if there is not already if args.RIVETREFS: try: - refhistos2 = getRivetRefData() + refhistos2 = getRivetRefData(anas) except IOError as e: print("File reading error: ", e.strerror) exit(1) refhistos2.update(refhistos) refhistos = refhistos2 ## Purge unmatched ref data entries to save memory keylist = list(refhistos.keys()) # can't modify for-loop target for refhpath in keylist: if refhpath not in hpaths: del refhistos[refhpath] ## Now loop over all MC histograms and plot them # TODO: factorize much of this into a rivet.utils mkplotfile(mchists, refhist, kwargs, is2d=False) function for hpath in hpaths: #print('Currently looking at', hpath) ## The analysis objects to be plotted anaobjects = [] ## List of histos to be drawn, to sync the legend and plotted lines mainlines = [] varlines = [] ## Is this a 2D histo? is2d = (hpath in h2ds) ## Will we be drawing a ratio plot? showratio = args.RATIO and not is2d ## A Plot object to represent the PLOT section in the .dat file plot = Plot() if not is2d: plot['Legend'] = '1' plot['LogY'] = '1' headers = plotparser.getHeaders(hpath) if headers: plot.update(headers) # for key, val in headers.items(): # plot[key] = val if "PLOT" in plotoptions: for key_val in plotoptions["PLOT"]: key, val = [s.strip() for s in key_val.split("=", 1)] if 'ReplaceOption' in key_val: opt_group = key_val.split("=", 2) key, val = "=".join(opt_group[:2]), opt_group[2] plot[key] = val if args.REMOVE_OPTIONS: plot['RemoveOptions'] = '1' if args.LINEAR: plot['LogY'] = '0' if args.NOPLOTTITLE: plot['Title'] = '' if showratio and args.RATIO_DEVIATION: plot['RatioPlotMode'] = 'deviation' if args.STYLE == 'talk': plot['PlotSize'] = '8,6' elif args.STYLE == 'bw' and showratio: plot['RatioPlotErrorBandColor'] = 'black!10' ## Get a special object, if there is one for this path special = plotparser.getSpecial(hpath) ## Handle reference data histogram, if there is one ratioreference, hasdataref = None, False if hpath in refhistos: hasdataref = True refdata = refhistos[hpath] refdata.setAnnotation('Title', args.REFTITLE) if not is2d: refdata.setAnnotation('ErrorBars', '1') refdata.setAnnotation('PolyMarker', '*') refdata.setAnnotation('ConnectBins', '0') if showratio: ratioreference = hpath ## For 1D anaobjects.append(refdata) mainlines.append(hpath) ## For 2D if is2d: s = mkOutput(hpath, [refdata], plot, special) writeOutput(s, hpath) ## Loop over the MC files to plot all instances of the histogram styleidx = 0 for infile in filelist: if infile in mchistos: for xpath in sorted(mchistos[infile]): if rivet.stripOptions(xpath) != hpath: continue hmcs = mchistos[infile][xpath] ## For now, just plot all the different variation histograms (reversed, so [0] is on top) # TODO: calculate and plot an appropriate error band, somehow... for i in sorted(hmcs.keys(), reverse=True): iscanonical = (str(i) == "0") hmc = hmcs[i] ## Default linecolor, linestyle if not is2d: setStyle(hmc, styleidx, not iscanonical) if args.MC_ERRS: hmc.setAnnotation('ErrorBars', '1') ## Plot defaults from .plot files histopts = plotparser.getHistogramOptions(hpath) if histopts: for key, val in histargs.items(): hmc.setAnnotation(key, val) ## Command line plot options setOptions(hmc, plotoptions[infile]) ## Set path attribute fullpath = "/"+infile+xpath if not iscanonical: fullpath += "["+str(i)+"]" hmc.setAnnotation('Path', fullpath) ## Add object / path to appropriate lists #if hmc.hasAnnotation("Title"): # hmc.setAnnotation("Title", hmc.annotation("Title") + # rivet.extractOptionString(xpath)) anaobjects.append(hmc) if iscanonical: mainlines.append(fullpath) else: varlines.append(fullpath) if showratio and ratioreference is None and iscanonical: ratioreference = fullpath ## For 2D, plot each histo now (since overlay makes no sense) if is2d: s = mkOutput(hpath, [hmc], plot, special) writeOutput(s, fullpath) styleidx += 1 ## Finally render the combined plots; only show the first one if it's 2D # TODO: Only show the first *MC* one if 2D? if is2d: anaobjects = anaobjects[:1] ## Add final attrs to Plot - if not plot.has_key('DrawOnly'): + if 'DrawOnly' not in plot: plot['DrawOnly'] = ' '.join(varlines + mainlines).strip() - if not plot.has_key('LegendOnly'): + if 'LegendOnly' not in plot: plot['LegendOnly'] = ' '.join(mainlines).strip() - if showratio and not plot.has_key('RatioPlot') and len(varlines + mainlines) > 1: + if showratio and 'RatioPlot' not in plot and len(varlines + mainlines) > 1: plot['RatioPlot'] = '1' - if showratio and plot.has_key('RatioPlot') and plot['RatioPlot'] == '1' and len(varlines + mainlines) > 0: - if not plot.has_key('RatioPlotReference'): + if showratio and 'RatioPlot' in plot and plot['RatioPlot'] == '1' and len(varlines + mainlines) > 0: + if 'RatioPlotReference' not in plot: plot['RatioPlotReference'] = ratioreference if not hasdataref and "RatioPlotYLabel" not in plot: if plot.get('RatioPlotMode', '') == 'deviation': plot['RatioPlotYLabel'] = 'Deviation' #r'$\text{MC}-\text{MC}_\text{ref}$' else: plot['RatioPlotYLabel'] = 'Ratio' #r'$\text{MC}/\text{MC}_\text{ref}$' ## Make the output and write to file o = mkOutput(hpath, anaobjects, plot, special) writeOutput(o, hpath)