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,202 @@ // -*- 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"); + 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 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()); 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){ 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/analyses/pluginALICE/EvMixing_TestAnalysis.cc b/analyses/pluginALICE/EvMixing_TestAnalysis.cc new file mode 100644 --- /dev/null +++ b/analyses/pluginALICE/EvMixing_TestAnalysis.cc @@ -0,0 +1,149 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" +#include "Rivet/Projections/FinalState.hh" +#include "Rivet/Projections/FastJets.hh" +#include "Rivet/Projections/EventMixingFinalState.hh" +#include "Rivet/Math/Constants.hh" +#include "Rivet/Tools/AliceCommon.hh" +#include "Rivet/Projections/AliceCommon.hh" + +namespace Rivet { + + /// @brief Add a short analysis description here + class EvMixing_TestAnalysis : public Analysis { + public: + + /// Constructor + DEFAULT_RIVET_ANALYSIS_CTOR(EvMixing_TestAnalysis); + + + /// @name Analysis methods + //@{ + + /// Book histograms and initialise projections before the run + void init() { + + // Initialise and register projections + const FinalState fs(Cuts::abseta < 0.8 && Cuts::pT > 150*MeV); + declare(fs, "FS"); + const FinalState fsmix(Cuts::abseta < 0.8 && Cuts::pT > 150*MeV); + declare(fsmix, "MIX"); + + // Declare centrality projection + const CentralityProjection& cp = declareCentrality(ALICE::V0MMultiplicity(), + "ALICE_2015_PBPBCentrality", "V0M", "V0M"); + + // Declare event mixing projection + declare(EventMixingCentrality(&cp, fsmix, 10, 0, 100, 10), "EVMIXFS"); + + // Book histograms + for (int i = 0; i < NBINS; i++) { + _hDEtaDPhiSame[i] = bookHisto2D("detadphisame" + toString(i), 36, -0.5*pi, 1.5*pi, 32, -1.6, 1.6, + "Same", "$\\Delta\\phi$ (rad)", "$\\Delta\\eta$", + "$1 / N_{same} {\\rm d}N_{same} / {\\rm d}\\Delta\\eta{\\rm d}\\Delta\\phi$ (rad$^-1$)"); + _hDEtaDPhiMixed[i] = bookHisto2D("detadphimixed" + toString(i), 36, -0.5*pi, 1.5*pi, 32, -1.6, 1.6, + "Mixed", "$\\Delta\\phi$ (rad)", "$\\Delta\\eta$", + "$1 / N_{mixed} {\\rm d}N_{mixed} / {\\rm d}\\Delta\\eta{\\rm d}\\Delta\\phi$ (rad$^-1$)"); + _hDEtaDPhi[i] = bookHisto2D("detadphi" + toString(i), 36, -0.5*pi, 1.5*pi, 32, -1.6, 1.6, + "DEtaDPhi", "$\\Delta\\phi$ (rad)", "$\\Delta\\eta$", + "$1 / N {\\rm d}N / {\\rm d}\\Delta\\eta{\\rm d}\\Delta\\phi$ (rad$^-1$)"); + } + + } + /// @brief Calculate angular distance between particles. + double deltaPhi(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; + } + + + + /// Perform the per-event analysis + void analyze(const Event& event) { + + // Get centrality bin number for current event + int nr = (int)(apply(event, "V0M")() / 10); + if (nr > NBINS) vetoEvent; + + // Get particles from current event + const FinalState& fs = applyProjection(event, "FS"); + Particles particles = fs.particles(); + + // Fill 'same' histograms + for (const Particle& particle1 : particles) { + for (const Particle& particle2 : particles) { + if (particle1 != particle2) { + double dPhi = deltaPhi(particle1.phi(), particle2.phi()); + if (dPhi < -0.5 * pi) dPhi += 2 * pi; + double dEta = abs(particle1.eta()-particle2.eta()); + _hDEtaDPhiSame[nr]->fill(dPhi, dEta, event.weight()); + } + } + } + + // Get events for mixing + const EventMixingCentrality& evmixfs = applyProjection(event, "EVMIXFS"); + vector mixedEvents = evmixfs.getMixingEvents(); + + // Fill 'mixed' histograms + if (mixedEvents.size() > 0) { // not really needed + for (const auto& mixedParticles : mixedEvents) { + //std::cout << particles.size() << "\t" << mixedParticles.size() << std::endl; + for (const Particle& particle1 : particles) { + for (const Particle& particle2 : mixedParticles) { + if (particle1 != particle2) { + double dPhi = deltaPhi(particle1.phi(), particle2.phi()); + if (dPhi < -0.5 * pi) dPhi += 2 * pi; + double dEta = abs(particle1.eta()-particle2.eta()); + _hDEtaDPhiMixed[nr]->fill(dPhi, dEta, 1. / mixedEvents.size()); + } + } + } + } + //std::cout << std::endl; + } + + } + + + /// Normalise histograms etc., after the run + void finalize() { + + for (int ibin = 0; ibin < NBINS; ibin++) { + if (_hDEtaDPhiMixed[ibin]->numEntries() < 1) + continue; + std::cout << "Creating final histogram " << ibin << std::endl; + // We have to do that manually (bin by bin) because there is no way to use divide() method on Histo2DPtr + // and it is not possible to use Scatter3DPtr because there is no bookScatter3DPtr method + for (unsigned int ixy = 0; ixy < _hDEtaDPhiSame[ibin]->numBins(); ixy++) { + _hDEtaDPhi[ibin]->fillBin(ixy, _hDEtaDPhiSame[ibin]->bin(ixy).sumW() / _hDEtaDPhiMixed[ibin]->bin(ixy).sumW()); + } + } + + } + + //@} + + private: + + static const int NBINS = 10; + + /// @name Histograms + //@{ + Histo2DPtr _hDEtaDPhiSame[NBINS]; + Histo2DPtr _hDEtaDPhiMixed[NBINS]; + Histo2DPtr _hDEtaDPhi[NBINS]; + //@} + + }; + + + // The hook for the plugin system + DECLARE_RIVET_PLUGIN(EvMixing_TestAnalysis); + + +} 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,139 @@ // -*- C++ -*- #ifndef RIVET_EventMixingFinalState_HH #define RIVET_EventMixingFinalState_HH #include "Rivet/Projection.hh" #include "Rivet/Projections/ParticleFinder.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. // // @author Christian Bierlich typedef map > MixMap; - class EventMixingFinalState : public Projection { + // 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 { public: // Constructor - EventMixingFinalState(const ParticleFinder& fsp, const ParticleFinder& mix, size_t nMixIn, + EventMixingBase(const Projection* mixObsProjPtr, const ParticleFinder& mix, size_t nMixIn, double oMin, double oMax, double deltao ) : nMix(nMixIn){ - setName("EventMixingFinalState"); - addProjection(fsp,"FS"); + setName("EventMixingBase"); + addProjection(*mixObsProjPtr,"OBS"); addProjection(mix,"MIX"); - MSG_WARNING("EventMixingFinalState is a naive implementation, not currently " << + MSG_WARNING("EventMixing is a naive implementation, not currently " << "validated. Use with caution."); - // Set up the map for mixing events + // 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 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); } protected: // Calulate mixing observable. // Can be overloaded to define more elaborate observables. - virtual void calculateMixingObs(const Particles& parts){ - mObs = parts.size(); - } + virtual void calculateMixingObs(const Projection* mProj) = 0; /// Perform the projection on the Event void project(const Event& e){ - const Particles parts = applyProjection(e, "FS").particles(); - - calculateMixingObs(parts); - + 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(); } /// Compare with other projections int compare(const Projection& p) const { - return mkNamedPCmp(p,"FS"); + 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 mixing observable of the current event - double mObs; // The event map; MixMap mixEvents; }; + + // 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))(); + } + }; } #endif