diff --git a/analyses/pluginALICE/ALICE_2016_I1507157.cc b/analyses/pluginALICE/ALICE_2016_I1507157.cc --- a/analyses/pluginALICE/ALICE_2016_I1507157.cc +++ b/analyses/pluginALICE/ALICE_2016_I1507157.cc @@ -1,201 +1,186 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" -#include "Rivet/Tools/AliceCommon.hh" +#include "Rivet/Projections/AliceCommon.hh" #include "Rivet/Projections/PrimaryParticles.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/EventMixingFinalState.hh" namespace Rivet { /// @brief Correlations of identified particles in pp. /// Also showcasing use of EventMixingFinalState. + class ALICE_2016_I1507157 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ALICE_2016_I1507157); /// @name Analysis methods //@{ - + /// @brief Calculate angular distance between particles. double phaseDif(double a1, double a2){ double dif = a1 - a2; while (dif < -M_PI/2) dif += 2*M_PI; while (dif > 3*M_PI/2) dif -= 2*M_PI; return dif; } + /// Book histograms and initialise projections before the run void init() { double etamax = 0.8; - double pTmin = 0.5; // GeV - + double pTmin = 0.5; // GeV + // Trigger declare(ALICE::V0AndTrigger(), "V0-AND"); // Charged tracks used to manage the mixing observable. ChargedFinalState cfsMult(Cuts::abseta < etamax); declare(cfsMult, "CFSMult"); - + // Primary particles. - PrimaryParticles pp({Rivet::PID::PIPLUS, Rivet::PID::KPLUS, - Rivet::PID::K0S, Rivet::PID::K0L, Rivet::PID::PROTON, + PrimaryParticles pp({Rivet::PID::PIPLUS, Rivet::PID::KPLUS, + Rivet::PID::K0S, Rivet::PID::K0L, Rivet::PID::PROTON, Rivet::PID::NEUTRON, Rivet::PID::LAMBDA, Rivet::PID::SIGMAMINUS, - Rivet::PID::SIGMAPLUS, Rivet::PID::XIMINUS, Rivet::PID::XI0, + Rivet::PID::SIGMAPLUS, Rivet::PID::XIMINUS, Rivet::PID::XI0, Rivet::PID::OMEGAMINUS},Cuts::abseta < etamax && Cuts::pT > pTmin*GeV); declare(pp,"APRIM"); // The event mixing projection - declare(EventMixingFinalState(cfsMult, pp, 5, 0, 100, 10),"EVM"); + declare(EventMixingFinalState(&cfsMult, pp, 5, 0, 100, 10),"EVM"); // The particle pairs. pid = {{211, -211}, {321, -321}, {2212, -2212}, {3122, -3122}, {211, 211}, {321, 321}, {2212, 2212}, {3122, 3122}, {2212, 3122}, {2212, -3122}}; // The associated histograms in the data file. vector refdata = {"d04-x01-y01","d04-x01-y02","d04-x01-y03", "d06-x01-y02","d05-x01-y01","d05-x01-y02","d05-x01-y03","d06-x01-y01", "d01-x01-y02","d02-x01-y02"}; ratio.resize(refdata.size()); signal.resize(refdata.size()); background.resize(refdata.size()); for (int i = 0, N = refdata.size(); i < N; ++i) { // The ratio plots. - book(ratio[i], refdata[i], true); - // Signal and mixed background. - book(signal[i], "/TMP/" + refdata[i] + "-s", *ratio[i], refdata[i] + "-s"); - book(background[i], "/TMP/" + refdata[i] + "-b", *ratio[i], refdata[i] + "-b"); + book(ratio[i],refdata[i], true); + // Signal and mixed background. + book(signal[i],"/TMP/" + refdata[i] + "-s", *ratio[i], refdata[i] + "-s"); + book(background[i],"/TMP/" + refdata[i] + "-b", *ratio[i], refdata[i] + "-b"); // Number of signal and mixed pairs. - nsp.push_back(0.); + nsp.push_back(0.); nmp.push_back(0.); - } + } } /// Perform the per-event analysis void analyze(const Event& event) { - // Triggering if (!apply(event, "V0-AND")()) return; // The projections - const PrimaryParticles& pp = applyProjection(event,"APRIM"); - const EventMixingFinalState& evm = + const PrimaryParticles& pp = + applyProjection(event,"APRIM"); + const EventMixingFinalState& evm = applyProjection(event, "EVM"); - // Get all mixing events - vector mixEvents = evm.getMixingEvents(); - - // If there are not enough events to mix, don't fill any histograms - if(mixEvents.size() == 0) - return; - - // Make a vector of mixed event particles - vector mixParticles; - size_t pSize = 0; - for(size_t i = 0; i < mixEvents.size(); ++i) - pSize+=mixEvents[i].size(); - mixParticles.reserve(pSize); - for(size_t i = 0; i < mixEvents.size(); ++i) - mixParticles.insert(mixParticles.end(), mixEvents[i].begin(), - mixEvents[i].end()); - - // Shuffle the particles in the mixing events - random_shuffle(mixParticles.begin(), mixParticles.end()); + // Test if we have enough mixing events available to continue. + if (!evm.hasMixingEvents()) return; for(const Particle& p1 : pp.particles()) { // Start by doing the signal distributions for(const Particle& p2 : pp.particles()) { if(p1 == p2) continue; double dEta = abs(p1.eta() - p2.eta()); double dPhi = phaseDif(p1.phi(), p2.phi()); if(dEta < 1.3) { for (int i = 0, N = pid.size(); i < N; ++i) { int pid1 = pid[i].first; int pid2 = pid[i].second; bool samesign = (pid1 * pid2 > 0); - if (samesign && ((pid1 == p1.pid() && pid2 == p2.pid()) || + if (samesign && ((pid1 == p1.pid() && pid2 == p2.pid()) || (pid1 == -p1.pid() && pid2 == -p2.pid()))) { signal[i]->fill(dPhi); nsp[i] += 1.0; } if (!samesign && abs(pid1) == abs(pid2) && pid1 == p1.pid() && pid2 == p2.pid()) { signal[i]->fill(dPhi); nsp[i] += 1.0; } - if (!samesign && abs(pid1) != abs(pid2) && + if (!samesign && abs(pid1) != abs(pid2) && ( (pid1 == p1.pid() && pid2 == p2.pid()) || (pid2 == p1.pid() && pid1 == p2.pid()) ) ) { signal[i]->fill(dPhi); nsp[i] += 1.0; } } } } // Then do the background distribution - for(const Particle& pMix : mixParticles){ + for(const Particle& pMix : evm.particles()){ double dEta = abs(p1.eta() - pMix.eta()); double dPhi = phaseDif(p1.phi(), pMix.phi()); if(dEta < 1.3) { for (int i = 0, N = pid.size(); i < N; ++i) { int pid1 = pid[i].first; int pid2 = pid[i].second; bool samesign = (pid1 * pid2 > 0); - if (samesign && ((pid1 == p1.pid() && pid2 == pMix.pid()) || + if (samesign && ((pid1 == p1.pid() && pid2 == pMix.pid()) || (pid1 == -p1.pid() && pid2 == -pMix.pid()))) { background[i]->fill(dPhi); nmp[i] += 1.0; } if (!samesign && abs(pid1) == abs(pid2) && pid1 == p1.pid() && pid2 == pMix.pid()) { background[i]->fill(dPhi); nmp[i] += 1.0; } - if (!samesign && abs(pid1) != abs(pid2) && + if (!samesign && abs(pid1) != abs(pid2) && ( (pid1 == p1.pid() && pid2 == pMix.pid()) || (pid2 == p1.pid() && pid1 == pMix.pid()) ) ) { background[i]->fill(dPhi); nmp[i] += 1.0; } } } } } } /// Normalise histograms etc., after the run void finalize() { for (int i = 0, N = pid.size(); i < N; ++i) { double sc = nmp[i] / nsp[i]; - signal[i]->scaleW(sc); - divide(signal[i],background[i],ratio[i]); + signal[i]->scaleW(sc); + divide(signal[i],background[i],ratio[i]); } } //@} /// @name Histograms //@{ vector > pid; vector signal; vector background; vector ratio; vector nsp; vector nmp; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2016_I1507157); } diff --git a/include/Rivet/Projections/EventMixingFinalState.hh b/include/Rivet/Projections/EventMixingFinalState.hh --- a/include/Rivet/Projections/EventMixingFinalState.hh +++ b/include/Rivet/Projections/EventMixingFinalState.hh @@ -1,104 +1,225 @@ // -*- C++ -*- #ifndef RIVET_EventMixingFinalState_HH #define RIVET_EventMixingFinalState_HH #include "Rivet/Projection.hh" #include "Rivet/Projections/ParticleFinder.hh" +#include "Rivet/Tools/Random.hh" #include #include + namespace Rivet { - using std::deque; + // @brief Projects out an event mixed of several events, given + // a mixing observable (eg. number of final state particles), + // defining what should qualify as a mixable event. + // Binning in the mixing observable is defined in the constructor, + // as is the number of events one wants to mix with. + // The method calculateMixingObs() must can be overloaded + // in derived classes, to provide the definition of the mixing observable, + // on the provided projection, eg. centrality or something more elaborate. + // + // The implementation can correcly handle mixing of weighted events, but + // not multi-weighted events. Nor does the implementation attemt to handle + // calculation of one (or more) event weights for the mixed events. For most + // common use cases, like calculating a background sample, this is sufficient. + // If a more elaborate use case ever turns up, this must be reevaluated. + // + // + // @author Christian Bierlich + // Weighted random shuffle, similar to std::random_shuffle, which + // allows the passing of a weight for each element to be shuffled. + template + void weighted_shuffle(RandomAccessIterator first, RandomAccessIterator last, + WeightIterator fw, WeightIterator lw, RandomNumberGenerator& g) { + while(first != last && fw != lw) { + std::discrete_distribution weightDist(fw, lw); + int i = weightDist(g); + if(i){ + std::iter_swap(first, next(first, i)); + std::iter_swap(fw, next(fw, i)); + } + ++first; + ++fw; + } + } + // A MixEvent is a vector of particles with and associated weight. + typedef pair MixEvent; + typedef map > MixMap; + // EventMixingBase is the base class for event mixing projections. + // Most methods are defined in this base class as they should. + // In order to use it, a derived class should be implemented where: + // - The constructor is reimplmented, giving the derived projection type + // from which the mixing observable is calculated. The constructor must also + // be declared public in the derived class. + // - The calculateMixingObs is implemented. + // To examples of such derived classes are given below, + // 1) EventMixingFinalState, where the mixing observable are calculated + // on a multiplicity of a charged final state, and: + // 2) EventMixingCentrality, where the mixing observable is centrality. - /// Projects out an event mixed of several events, given - /// a mixing observable (eg. number of final state particles), - /// defining what should qualify as a mixable event. - /// Binning in the mixing observable is defined in the constructor, - /// as is the number of events one wants to mix with. - /// The protected method calculateMixingObs() can be overloaded - /// in derived classes, to define other mixing observables, eg. - /// centrality or something even more elaborate. - /// - /// @warning !!!DISCLAIMER!!! The projection makes no attempt at correct handling - /// of event weights - ie. what one should use as event weight for several - /// mixed events. Proceed with caution if you do not use an MC with - /// unit weights. - /// - /// @author Christian Bierlich - /// - typedef map > MixMap; - class EventMixingFinalState : public Projection { - public: + class EventMixingBase : public Projection { + protected: // Constructor - EventMixingFinalState(const ParticleFinder& fsp, const ParticleFinder& mix, size_t nMixIn, - double oMin, double oMax, double deltao ) : nMix(nMixIn){ - setName("EventMixingFinalState"); - declare(fsp, "FS"); - declare(mix, "MIX"); - MSG_WARNING("EventMixingFinalState is a naive implementation, not currently " << - "validated. Use with caution."); + EventMixingBase(const Projection* mixObsProjPtr, const ParticleFinder& mix, + size_t nMixIn, double oMin, double oMax, double deltao) : nMix(nMixIn), + unitWeights(true) { + // The base class contructor should be called explicitly in derived classes + // to add projections below. + setName("EventMixingBase"); + declare(*mixObsProjPtr,"OBS"); + declare(mix,"MIX"); + MSG_WARNING("EventMixing is not fully validated. Use with caution."); - // Set up the map for mixing events - for (double o = oMin; o < oMax; o+=deltao) - mixEvents[o] = deque(); + // Set up the map for mixing events. + for(double o = oMin; o < oMax; o+=deltao ) + mixEvents[o] = std::deque(); } - // Clone on the heap - DEFAULT_RIVET_PROJ_CLONE(EventMixingFinalState); - - + public: + + // Test if we have enough mixing events available for projected, + // current mixing observable. + bool hasMixingEvents() const { + MixMap::const_iterator mixItr = mixEvents.lower_bound(mObs); + if(mixItr == mixEvents.end() || mixItr->second.size() < nMix + 1) + return false; + return true; + } + // Return a vector of mixing events. - vector getMixingEvents() const { - MixMap::const_iterator mixItr = mixEvents.lower_bound(mObs); - if(mixItr == mixEvents.end() || mixItr->second.size() < nMix + 1) - return vector(); - return vector(mixItr->second.begin(), mixItr->second.end() - 1); + vector getMixingEvents() const { + if (!hasMixingEvents()) + return vector(); + MixMap::const_iterator mixItr = mixEvents.lower_bound(mObs); + return vector(mixItr->second.begin(), mixItr->second.end() - 1); } + // Return a vector of particles from the mixing events. Can + // be overloaded in derived classes, though normally not neccesary. + virtual const Particles particles() const { + // Test if we have enough mixing events. + if(!hasMixingEvents()) + return vector(); + // Get mixing events for the current, projected mixing observable. + MixMap::const_iterator mixItr = mixEvents.lower_bound(mObs); + vector mixEvents = + vector(mixItr->second.begin(), mixItr->second.end() - 1); + // Make the vector of mixed particles. + Particles mixParticles; + vector weights; + size_t pSize = 0; + for(size_t i = 0; i < mixEvents.size(); ++i) + pSize+=mixEvents[i].first.size(); + mixParticles.reserve(pSize); + weights.reserve(pSize); + // Put the particles in the vector. + for(size_t i = 0; i < mixEvents.size(); ++i) { + mixParticles.insert(mixParticles.end(), mixEvents[i].first.begin(), + mixEvents[i].first.end()); + vector tmp = vector(mixEvents[i].first.size(), + mixEvents[i].second); + weights.insert(weights.end(), tmp.begin(), tmp.end()); + } + + // Shuffle the particles. + if (unitWeights) { + // Use the thread safe random number generator. + auto rnd = [&] (int i) {return rng()()%i;}; + random_shuffle(mixParticles.begin(), mixParticles.end(), rnd); + return mixParticles; + } + else { + weighted_shuffle(mixParticles.begin(), mixParticles.end(), + weights.begin(), weights.end(), rng()); + Particles tmp = vector(mixParticles.begin(), + mixParticles.begin() + size_t(ceil(mixParticles.size() / 2))); + return tmp; + } + } + protected: // Calulate mixing observable. - // Can be overloaded to define more elaborate observables. - virtual void calculateMixingObs(const Particles& parts){ - mObs = parts.size(); + // Must be overloaded in derived classes. + virtual void calculateMixingObs(const Projection* mProj) = 0; + + /// Perform the projection on the Event. + void project(const Event& e){ + const Projection* mixObsProjPtr = &applyProjection(e, "OBS"); + calculateMixingObs(mixObsProjPtr); + MixMap::iterator mixItr = mixEvents.lower_bound(mObs); + if(mixItr == mixEvents.end()){ + // We are out of bounds. + MSG_DEBUG("Mixing observable out of bounds."); + return; + } + const Particles mix = applyProjection(e, "MIX").particles(); + mixItr->second.push_back(make_pair(mix,e.weights()[0])); + // Assume unit weights until we see otherwise. + if (unitWeights && e.weights()[0] != 1.0 ) { + unitWeights = false; + nMix *= 2; + } + if(mixItr->second.size() > nMix + 1) + mixItr->second.pop_front(); } - /// Perform the projection on the Event - void project(const Event& e){ - const Particles parts = applyProjection(e, "FS").particles(); - - calculateMixingObs(parts); - - MixMap::iterator mixItr = mixEvents.lower_bound(mObs); - if(mixItr == mixEvents.end()){ - // We are out of bounds. - MSG_DEBUG("Mixing observable out of bounds."); - return; - } - const Particles mix = applyProjection(e, "MIX").particles(); - - mixItr->second.push_back(mix); - if(mixItr->second.size() > nMix + 1) - mixItr->second.pop_front(); + /// Compare with other projections. + CmpState compare(const Projection& p) const { + return mkNamedPCmp(p,"OBS"); } - - /// Compare with other projections - CmpState compare(const Projection& p) const { - return mkNamedPCmp(p, "FS"); - } - + + // The mixing observable of the current event. + double mObs; + private: - // The number of event to mix with + // The number of event to mix with. size_t nMix; - // The mixing observable of the current event - double mObs; - // The event map; + // The event map. MixMap mixEvents; + // Using unit weights or not. + bool unitWeights; }; + // EventMixingFinalState has multiplicity in the mixing projection. + // as the mixing observable. + class EventMixingFinalState : public EventMixingBase { + public: + EventMixingFinalState(const ParticleFinder* mixObsProjPtr, + const ParticleFinder& mix, size_t nMixIn, double oMin, double oMax, + double deltao) : + EventMixingBase(mixObsProjPtr, mix, nMixIn, oMin, oMax, deltao) { + setName("EventMixingFinalState"); + } + DEFAULT_RIVET_PROJ_CLONE(EventMixingFinalState); + + protected: + // Calulate mixing observable. + virtual void calculateMixingObs(const Projection* mProj) { + mObs = ((ParticleFinder*) mProj)->particles().size(); + } + }; + + // EventMixingCentrality has centrality as the mixing observable. + class EventMixingCentrality : public EventMixingBase { + public: + EventMixingCentrality(const CentralityProjection* mixObsProjPtr, + const ParticleFinder& mix, size_t nMixIn, double oMin, double oMax, + double deltao) : + EventMixingBase(mixObsProjPtr, mix, nMixIn, oMin, oMax, deltao) { + setName("EventMixingCentrality"); + } + + DEFAULT_RIVET_PROJ_CLONE(EventMixingCentrality); + protected: + virtual void calculateMixingObs(const Projection* mProj) { + mObs = ((CentralityProjection*) mProj)->operator()(); + } + }; } - #endif