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,167 +1,187 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" -#include "Rivet/Projections/AliceCommon.hh" +#include "Rivet/Tools/AliceCommon.hh" #include "Rivet/Projections/PrimaryParticles.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/EventMixingFinalState.hh" namespace Rivet { - /// @brief ALICE correlations of identified particles in pp + /// @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); 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); + 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"); + 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"}; + "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. + 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")); + "-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.); + 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 EventMixingFinalState& evm = + applyProjection(event, "EVM"); - // The primary particles - const PrimaryParticles& pp = apply(event, "APRIM"); - const Particles pparticles = pp.particles(); + // Test if we have enough mixing events available to continue. + if (!evm.hasMixingEvents()) return; - // The mixed events - const EventMixingFinalState& evm = apply(event, "EVM"); - const vector mixEvents = evm.getMixingEvents(); - if (mixEvents.empty()) vetoEvent; - - // 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()); - random_shuffle(mixParticles.begin(), mixParticles.end()); - - for (size_t ip1 = 0; ip1 < pparticles.size()-1; ++ip1) { - const Particle& p1 = pparticles[ip1]; - + for(const Particle& p1 : pp.particles()) { // Start by doing the signal distributions - for (size_t ip2 = 0; ip2 < pparticles.size(); ++ip1) { - if (ip1 == ip2) continue; - const Particle& p2 = pparticles[ip2]; - const double dEta = deltaEta(p1, p2); - const double dPhi = deltaPhi(p1, p2, true); - if (dEta > 1.3) continue; - for (int i = 0, N = pid.size(); i < N; ++i) { - const int pid1 = pid[i].first; - const int pid2 = pid[i].second; - const bool samesign = (pid1 * pid2 > 0); - const bool pidmatch1 = (pid1 == p1.pid() && pid2 == p2.pid()) || (pid1 == -p1.pid() && pid2 == -p2.pid()); - const bool pidmatch2 = abs(pid1) == abs(pid2) && pid1 == p1.pid() && pid2 == p2.pid(); - const bool pidmatch3 = abs(pid1) != abs(pid2) && ( (pid1 == p1.pid() && pid2 == p2.pid()) || (pid2 == p1.pid() && pid1 == p2.pid()) ); - if ((samesign && pidmatch1) || (!samesign && (pidmatch2 || pidmatch3))) { - signal[i]->fill(dPhi, weight); - nsp[i] += 1.0; - } - } - } - - // Then do the background distribution - for (const Particle& pMix : mixParticles){ - const double dEta = deltaEta(p1, pMix); - const double dPhi = deltaPhi(p1, pMix, true); - if (dEta > 1.3) continue; - for (int i = 0, N = pid.size(); i < N; ++i) { - const int pid1 = pid[i].first; - const int pid2 = pid[i].second; - const bool samesign = (pid1 * pid2 > 0); - const bool pidmatch1 = (pid1 == p1.pid() && pid2 == pMix.pid()) || (pid1 == -p1.pid() && pid2 == -pMix.pid()); - const bool pidmatch2 = abs(pid1) == abs(pid2) && pid1 == p1.pid() && pid2 == pMix.pid(); - const bool pidmatch3 = abs(pid1) != abs(pid2) && ( (pid1 == p1.pid() && pid2 == pMix.pid()) || (pid2 == p1.pid() && pid1 == pMix.pid()) ); - if ((samesign && pidmatch1) || (!samesign && (pidmatch2 || pidmatch3))) { - background[i]->fill(dPhi, weight); - nmp[i] += 1.0; - } - } - } + for(const Particle& p2 : pp.particles()) { + if(isSame(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 : 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) { - const double sc = nmp[i] / nsp[i]; - signal[i]->scaleW(sc); - divide(signal[i],background[i],ratio[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,99 +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 protected method calculateMixingObs() can be overloaded - // in derived classes, to define other mixing observables, eg. - // centrality or something even 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. + // 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. // // @author Christian Bierlich - typedef map > MixMap; - class EventMixingFinalState : public Projection { + // 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), + 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(); + } + public: - // Constructor - EventMixingFinalState(const ParticleFinder& fsp, const ParticleFinder& mix, size_t nMixIn, - double oMin, double oMax, double deltao ) : nMix(nMixIn){ - setName("EventMixingFinalState"); - addProjection(fsp,"FS"); - addProjection(mix,"MIX"); - MSG_WARNING("EventMixingFinalState is a naive implementation, not currently " << - "validated. Use with caution."); + + // 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 { + if (!hasMixingEvents()) + return vector(); + MixMap::const_iterator mixItr = mixEvents.lower_bound(mObs); + return vector(mixItr->second.begin(), mixItr->second.end() - 1); + } - // Set up the map for mixing events - for(double o = oMin; o < oMax; o+=deltao ) - mixEvents[o] = deque(); - - } - // Clone on the heap - DEFAULT_RIVET_PROJ_CLONE(EventMixingFinalState); + // 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()); + } - - // 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); + // 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.weight())); + // Assume unit weights until we see otherwise. + if (unitWeights && e.weight() != 1.0 ) { + unitWeights = false; + nMix *= 2; + } + if(mixItr->second.size() > nMix + 1) + mixItr->second.pop_front(); + } + + /// Compare with other projections. + int compare(const Projection& p) const { + return mkNamedPCmp(p,"OBS"); } - /// Perform the projection on the Event - virtual void project(const Event& e){ - const Particles parts = applyProjection(e, "FS").particles(); + // 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; + }; - calculateMixingObs(parts); + // 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"); + } - 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(); + DEFAULT_RIVET_PROJ_CLONE(EventMixingFinalState); + + protected: + // Calulate mixing observable. + virtual void calculateMixingObs(const Projection* mProj) { + mObs = ((ParticleFinder*) mProj)->particles().size(); } + }; - /// Compare with other projections - virtual int compare(const Projection& p) const { - return mkNamedPCmp(p,"FS"); - } + // 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"); + } - protected: - // The number of event to mix with - size_t nMix; - // The mixing observable of the current event - double mObs; - // The event map; - MixMap mixEvents; - }; + DEFAULT_RIVET_PROJ_CLONE(EventMixingCentrality); + protected: + virtual void calculateMixingObs(const Projection* mProj) { + mObs = ((CentralityProjection*) mProj)->operator()(); + } + }; } #endif