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,202 +1,187 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Tools/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 // Trigger declare(ALICE::V0AndTrigger(), "V0-AND"); // Charged tracks used to manage the mixing observable. ChargedFinalState cfsMult(Cuts::abseta < etamax); addProjection(cfsMult, "CFSMult"); // Primary particles. 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::OMEGAMINUS},Cuts::abseta < etamax && Cuts::pT > pTmin*GeV); addProjection(pp,"APRIM"); // The event mixing projection 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"}; for (int i = 0, N = refdata.size(); i < N; ++i) { // The ratio plots. ratio.push_back(bookScatter2D(refdata[i], true)); // Signal and mixed background. signal.push_back(bookHisto1D("/TMP/" + refdata[i] + "-s", *ratio[i], refdata[i] + "-s")); background.push_back(bookHisto1D("/TMP/" + refdata[i] + "-b", *ratio[i], refdata[i] + "-b")); // Number of signal and mixed pairs. nsp.push_back(0.); nmp.push_back(0.); } } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); - + // Triggering if (!apply(event, "V0-AND")()) return; // The projections - const PrimaryParticles& pp = applyProjection(event,"APRIM"); + 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; + // Test if we have enough mixing events available to continue. + if (!evm.hasMixingEvents()) 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()); - 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()) || (pid1 == -p1.pid() && pid2 == -p2.pid()))) { signal[i]->fill(dPhi, weight); nsp[i] += 1.0; } if (!samesign && abs(pid1) == abs(pid2) && pid1 == p1.pid() && pid2 == p2.pid()) { signal[i]->fill(dPhi, weight); nsp[i] += 1.0; } if (!samesign && abs(pid1) != abs(pid2) && ( (pid1 == p1.pid() && pid2 == p2.pid()) || (pid2 == p1.pid() && pid1 == p2.pid()) ) ) { signal[i]->fill(dPhi, weight); 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()) || (pid1 == -p1.pid() && pid2 == -pMix.pid()))) { background[i]->fill(dPhi, weight); nmp[i] += 1.0; } if (!samesign && abs(pid1) == abs(pid2) && pid1 == p1.pid() && pid2 == pMix.pid()) { background[i]->fill(dPhi, weight); nmp[i] += 1.0; } if (!samesign && abs(pid1) != abs(pid2) && ( (pid1 == p1.pid() && pid2 == pMix.pid()) || (pid2 == p1.pid() && pid1 == pMix.pid()) ) ) { background[i]->fill(dPhi, weight); 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]); } } //@} /// @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,141 +1,218 @@ // -*- 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 { // @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 observables, // on the provided projection, eg. centrality or something more elaborate. // - // !!!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; + // 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) { + discrete_distribution weightDist(fw, lw); + int i = weightDist(g); + if(i){ + iter_swap(first, next(first, i)); + 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. class EventMixingBase : public Projection { protected: // Constructor - EventMixingBase(const Projection* mixObsProjPtr, const ParticleFinder& mix, size_t nMixIn, - double oMin, double oMax, double deltao ) : nMix(nMixIn){ - setName("EventMixingBase"); - addProjection(*mixObsProjPtr,"OBS"); - addProjection(mix,"MIX"); - MSG_WARNING("EventMixing 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"); + addProjection(*mixObsProjPtr,"OBS"); + addProjection(mix,"MIX"); + MSG_WARNING("EventMixing is not currently 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] = deque(); } 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 mixParticles; } protected: // Calulate mixing observable. // 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(mix); - if(mixItr->second.size() > nMix + 1) - mixItr->second.pop_front(); + 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.weight())); + // Assume unit weights until we see otherwise. + if (true && unitWeights && e.weight() != 1.0 ) { + unitWeights = false; + nMix *= 2.0; + } + if(mixItr->second.size() > nMix + 1) + mixItr->second.pop_front(); } /// Compare with other projections. int compare(const Projection& p) const { - return mkNamedPCmp(p,"OBS"); + return mkNamedPCmp(p,"OBS"); } // The mixing observable of the current event. double mObs; private: // The number of event to mix with. size_t nMix; // 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 ) : + 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 ) : + 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