diff --git a/analyses/pluginALICE/ALICE_2016_I1394676.cc b/analyses/pluginALICE/ALICE_2016_I1394676.cc --- a/analyses/pluginALICE/ALICE_2016_I1394676.cc +++ b/analyses/pluginALICE/ALICE_2016_I1394676.cc @@ -1,108 +1,107 @@ // -*- 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 eta distributions, peripheral events. class ALICE_2016_I1394676 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ALICE_2016_I1394676); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Initialise and register projections // Centrality projection. declareCentrality(ALICE::V0MMultiplicity(), "ALICE_2015_PBPBCentrality", "V0M","V0M"); // Projections for the 2-out-of-3 trigger. declare(ChargedFinalState( (Cuts::eta > 2.8 && Cuts::eta < 5.1) && Cuts::pT > 0.1*GeV), "VZERO1"); declare(ChargedFinalState( (Cuts::eta > -3.7 && Cuts::eta < -1.7) && Cuts::pT > 0.1*GeV), "VZERO2"); declare(ChargedFinalState(Cuts::abseta < 1. && Cuts::pT > 0.15*GeV), "SPD"); // Primary particles. declare(ALICE::PrimaryParticles(Cuts::abseta < 5.6),"APRIM"); // The centrality bins and the corresponding histograms and sow counters. centralityBins = { 40, 50, 60, 70, 80, 90 }; - vector< std::tuple > refData; for (int i = 0, N = centralityBins.size(); i < N; ++i) { histEta[centralityBins[i]] = bookHisto1D(1, 1, i + 1); sow[centralityBins[i]] = bookCounter("sow_" + toString(i)); } } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); // Trigger projections. const ChargedFinalState& vz1 = applyProjection(event,"VZERO1"); const ChargedFinalState& vz2 = applyProjection(event,"VZERO2"); const ChargedFinalState& spd = applyProjection(event,"SPD"); int fwdTrig = (vz1.particles().size() > 0 ? 1 : 0); int bwdTrig = (vz2.particles().size() > 0 ? 1 : 0); int cTrig = (spd.particles().size() > 0 ? 1 : 0); if (fwdTrig + bwdTrig + cTrig < 2) vetoEvent; const CentralityProjection& cent = apply(event,"V0M"); double c = cent(); // No centralities below 30 % if (c < 30.) return; // Find the correct centrality histogram auto hItr = histEta.upper_bound(c); if (hItr == histEta.end()) return; // Find the correct sow. auto sItr = sow.upper_bound(c); if (sItr == sow.end()) return; sItr->second->fill(weight); // Fill the histograms. for ( const auto& p : applyProjection(event,"APRIM").particles() ) if(p.abscharge() > 0) hItr->second->fill(p.eta(), weight); } /// Normalise histograms etc., after the run void finalize() { for (int i = 0, N = centralityBins.size(); i < N; ++i) histEta[centralityBins[i]]->scaleW(1./sow[centralityBins[i]]->sumW()); } //@} /// @name Histograms //@{ vector centralityBins; map histEta; map sow; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2016_I1394676); } 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,143 +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/MixedFinalState.hh" +#include "Rivet/Projections/EventMixingFinalState.hh" namespace Rivet { - /// @brief Add a short analysis description here + /// @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.1; // 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); addProjection(pp,"APRIM"); // The event mixing projection - declare(MixedFinalState(cfsMult, pp, 5, 0, 100, 10),"EVM"); - - // Book histograms - _h_protonRatio = bookScatter2D(2, 1, 2, true); - _h_protonSignal = bookHisto1D("protonSignal",*_h_protonRatio,"protonSignal"); - _h_protonBackground = bookHisto1D("protonBackground",*_h_protonRatio,"protonBackground"); - nmp = 0; - nsp = 0; - + 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 MixedFinalState& evm = applyProjection(event, "EVM"); + 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; + 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(); + 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()); + 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 distribution - for(const Particle& p2 : pp.particles() ){ - if(p1 == p2) - continue; - nsp+=1.0; - double dEta = abs(p1.eta() - p2.eta()); - double dPhi = phaseDif(p1.phi(), p2.phi()); - if(dEta < 1.3 && p1.pid() == 2212 && p2.pid() == -2212){ - _h_protonSignal->fill(dPhi,weight); - } + 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; } - // Then do the background distribution - for(const Particle& pMix : mixParticles){ - nmp+=1.0; - double dEta = abs(p1.eta() - pMix.eta()); - double dPhi = phaseDif(p1.phi(), pMix.phi()); - if(dEta < 1.3 && p1.pid() == 2212 && pMix.pid() == -2212){ - _h_protonBackground->fill(dPhi,weight); - } + 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() { - double sc = nmp / nsp; - scale(_h_protonSignal,sc); - divide(_h_protonSignal,_h_protonBackground,_h_protonRatio); - scale(_h_protonSignal,1.0/sc); + 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 //@{ - Histo1DPtr _h_protonSignal; - Histo1DPtr _h_protonBackground; - Scatter2DPtr _h_protonRatio; - double nmp, nsp; + 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/ALICE_2016_I1507157.info b/analyses/pluginALICE/ALICE_2016_I1507157.info --- a/analyses/pluginALICE/ALICE_2016_I1507157.info +++ b/analyses/pluginALICE/ALICE_2016_I1507157.info @@ -1,48 +1,51 @@ Name: ALICE_2016_I1507157 Year: 2016 -Summary: +Summary: Angular correlations of identified particles at 7 TeV. Experiment: ALICE Collider: LHC InspireID: 1507157 Status: UNVALIDATED Authors: - - Your Name -#References: -#- '' -#- '' -#- '' -RunInfo: + - Christian Bierlich +References: +- 'Eur.Phys.J.C77(2017)no.8,569' +- 'DOI:10.1140/epjc/s10052-017-5129-6' +- 'arXiv:1612.08975' +RunInfo: Proton-proton minium bias events at 7 TeV. NeedCrossSection: no -#Beams: -#Energies: +Beams: [p+,p+] +Energies: [7000] #Luminosity_fb: Description: - ' 50\;\GeV$.>' + 'Angular correlations between like-sign and opposite-sign identified particles, + integrated over $\Delta \eta < 1.3$. The analysis makes use of event mixing + to remove background. The current implementation of the event micing is not + validated by experiment, and should be used with caution. Note in particular + that for the event mixing to behave sensibly, event weights are assumed to be + unity. Do not run this analysis with weighted events' Keywords: [] BibKey: Adam:2016iwf BibTeX: '%%% contains utf-8, see: http://inspirehep.net/info/faq/general#utf8 %%% add \usepackage[utf8]{inputenc} to your latex preamble @article{Adam:2016iwf, author = "Adam, Jaroslav and others", title = "{Insight into particle production mechanisms via angular correlations of identified particles in pp collisions at $\sqrt{\mathrm{s}}=7$  TeV}", collaboration = "ALICE", journal = "Eur. Phys. J.", volume = "C77", year = "2017", number = "8", pages = "569", doi = "10.1140/epjc/s10052-017-5129-6", eprint = "1612.08975", archivePrefix = "arXiv", primaryClass = "nucl-ex", reportNumber = "CERN-EP-2016-322", SLACcitation = "%%CITATION = ARXIV:1612.08975;%%" }' ToDo: - Implement the analysis, test it, remove this ToDo, and mark as VALIDATED :-) diff --git a/analyses/pluginALICE/ALICE_2016_I1507157.plot b/analyses/pluginALICE/ALICE_2016_I1507157.plot --- a/analyses/pluginALICE/ALICE_2016_I1507157.plot +++ b/analyses/pluginALICE/ALICE_2016_I1507157.plot @@ -1,8 +1,61 @@ -BEGIN PLOT /ALICE_2016_I1507157/d01-x01-y01 -Title=[Insert title for histogram d01-x01-y01 here] -XLabel=[Insert $x$-axis label for histogram d01-x01-y01 here] -YLabel=[Insert $y$-axis label for histogram d01-x01-y01 here] -# + any additional plot settings you might like, see make-plots documentation +BEGIN PLOT /ALICE_2016_I1507157/d04-x01-y01 +Title=Integrated correlation function for $\pi^+ \pi^-$ +XLabel=$\Delta \phi$ (rad) +YLabel=$C(\Delta \phi)$ +LogY=0 +END PLOT +BEGIN PLOT /ALICE_2016_I1507157/d04-x01-y02 +Title=Integrated correlation function for $K^+ K^-$ +XLabel=$\Delta \phi$ (rad) +YLabel=$C(\Delta \phi)$ +LogY=0 +END PLOT +BEGIN PLOT /ALICE_2016_I1507157/d04-x01-y03 +Title=Integrated correlation function for $p\bar{p}$ +XLabel=$\Delta \phi$ (rad) +YLabel=$C(\Delta \phi)$ +LogY=0 +END PLOT +BEGIN PLOT /ALICE_2016_I1507157/d06-x01-y02 +Title=Integrated correlation function for $\Lambda \bar{\Lambda}$ +XLabel=$\Delta \phi$ (rad) +YLabel=$C(\Delta \phi)$ +LogY=0 END PLOT -# ... add more histograms as you need them ... +BEGIN PLOT /ALICE_2016_I1507157/d05-x01-y01 +Title=Integrated correlation function for $\pi^+ \pi^+ + \pi^- \pi^-$ +XLabel=$\Delta \phi$ (rad) +YLabel=$C(\Delta \phi)$ +LogY=0 +END PLOT +BEGIN PLOT /ALICE_2016_I1507157/d05-x01-y02 +Title=Integrated correlation function for $K^+ K^+ + K^- \K-$ +XLabel=$\Delta \phi$ (rad) +YLabel=$C(\Delta \phi)$ +LogY=0 +END PLOT +BEGIN PLOT /ALICE_2016_I1507157/d05-x01-y03 +Title=Integrated correlation function for $pp + \bar{p}\bar{p}$ +XLabel=$\Delta \phi$ (rad) +YLabel=$C(\Delta \phi)$ +LogY=0 +END PLOT +BEGIN PLOT /ALICE_2016_I1507157/d06-x01-y01 +Title=Integrated correlation function for $\Lambda \Lambda + \bar{\Lambda}\bar{\Lambda}$ +XLabel=$\Delta \phi$ (rad) +YLabel=$C(\Delta \phi)$ +LogY=0 +END PLOT +BEGIN PLOT /ALICE_2016_I1507157/d01-x01-y02 +Title=Integrated correlation function for $p \Lambda + \bar{p}\bar{\Lambda}$ +XLabel=$\Delta \phi$ (rad) +YLabel=$C(\Delta \phi)$ +LogY=0 +END PLOT +BEGIN PLOT /ALICE_2016_I1507157/d02-x01-y02 +Title=Integrated correlation function for $p \bar{\Lambda} + \bar{p}\Lambda$ +XLabel=$\Delta \phi$ (rad) +YLabel=$C(\Delta \phi)$ +LogY=0 +END PLOT diff --git a/include/Rivet/Makefile.am b/include/Rivet/Makefile.am --- a/include/Rivet/Makefile.am +++ b/include/Rivet/Makefile.am @@ -1,187 +1,187 @@ ## Internal headers - not to be installed nobase_dist_noinst_HEADERS = ## Public headers - to be installed nobase_pkginclude_HEADERS = ## Rivet interface nobase_pkginclude_HEADERS += \ Rivet.hh \ Run.hh \ Event.hh \ ParticleBase.hh \ Particle.fhh Particle.hh \ Jet.fhh Jet.hh \ Projection.fhh Projection.hh \ ProjectionApplier.hh \ ProjectionHandler.hh \ Analysis.hh \ AnalysisHandler.hh \ AnalysisInfo.hh \ AnalysisBuilder.hh \ AnalysisLoader.hh ## Build config stuff nobase_pkginclude_HEADERS += \ Config/RivetCommon.hh \ Config/RivetConfig.hh \ Config/BuildOptions.hh ## Projections nobase_pkginclude_HEADERS += \ Projections/AliceCommon.hh \ Projections/AxesDefinition.hh \ Projections/Beam.hh \ Projections/BeamThrust.hh \ Projections/CentralEtHCM.hh \ Projections/CentralityProjection.hh \ Projections/ChargedFinalState.hh \ Projections/ChargedLeptons.hh \ Projections/ConstLossyFinalState.hh \ Projections/DirectFinalState.hh \ Projections/DISFinalState.hh \ Projections/DISKinematics.hh \ Projections/DISLepton.hh \ Projections/DressedLeptons.hh \ + Projections/EventMixingFinalState.hh \ Projections/FastJets.hh \ Projections/PxConePlugin.hh \ Projections/FinalPartons.hh \ Projections/FinalState.hh \ Projections/FoxWolframMoments.hh \ Projections/FParameter.hh \ Projections/GeneratedPercentileProjection.hh \ Projections/HadronicFinalState.hh \ Projections/HeavyHadrons.hh \ Projections/Hemispheres.hh \ Projections/IdentifiedFinalState.hh \ Projections/ImpactParameterProjection.hh \ Projections/IndirectFinalState.hh \ Projections/InitialQuarks.hh \ Projections/InvMassFinalState.hh \ Projections/JetAlg.hh \ Projections/JetShape.hh \ Projections/LeadingParticlesFinalState.hh \ Projections/LossyFinalState.hh \ Projections/MergedFinalState.hh \ Projections/MissingMomentum.hh \ - Projections/MixedFinalState.hh \ Projections/NeutralFinalState.hh \ Projections/NonHadronicFinalState.hh \ Projections/NonPromptFinalState.hh \ Projections/ParisiTensor.hh \ Projections/ParticleFinder.hh \ Projections/PartonicTops.hh \ Projections/PercentileProjection.hh \ Projections/PrimaryHadrons.hh \ Projections/PrimaryParticles.hh \ Projections/PromptFinalState.hh \ Projections/SingleValueProjection.hh \ Projections/SmearedParticles.hh \ Projections/SmearedJets.hh \ Projections/SmearedMET.hh \ Projections/Sphericity.hh \ Projections/Spherocity.hh \ Projections/TauFinder.hh \ Projections/Thrust.hh \ Projections/TriggerCDFRun0Run1.hh \ Projections/TriggerCDFRun2.hh \ Projections/TriggerProjection.hh \ Projections/TriggerUA5.hh \ Projections/UnstableFinalState.hh \ Projections/UserCentEstimate.hh \ Projections/VetoedFinalState.hh \ Projections/VisibleFinalState.hh \ Projections/WFinder.hh \ Projections/ZFinder.hh ## Meta-projection convenience headers nobase_pkginclude_HEADERS += \ Projections/FinalStates.hh \ Projections/Smearing.hh ## Analysis base class headers # TODO: Move to Rivet/AnalysisTools header dir? nobase_pkginclude_HEADERS += \ Analyses/MC_Cent_pPb.hh \ Analyses/MC_ParticleAnalysis.hh \ Analyses/MC_JetAnalysis.hh \ Analyses/MC_JetSplittings.hh ## Tools nobase_pkginclude_HEADERS += \ Tools/AliceCommon.hh \ Tools/AtlasCommon.hh \ Tools/BeamConstraint.hh \ Tools/BinnedHistogram.hh \ Tools/CentralityBinner.hh \ Tools/Cmp.fhh \ Tools/Cmp.hh \ Tools/Correlators.hh \ Tools/Cutflow.hh \ Tools/Cuts.fhh \ Tools/Cuts.hh \ Tools/Exceptions.hh \ Tools/JetUtils.hh \ Tools/Logging.hh \ Tools/Random.hh \ Tools/ParticleBaseUtils.hh \ Tools/ParticleIdUtils.hh \ Tools/ParticleUtils.hh \ Tools/ParticleName.hh \ Tools/Percentile.hh \ Tools/PrettyPrint.hh \ Tools/RivetPaths.hh \ Tools/RivetSTL.hh \ Tools/RivetFastJet.hh \ Tools/RivetHepMC.hh \ Tools/RivetYODA.hh \ Tools/RivetMT2.hh \ Tools/SmearingFunctions.hh \ Tools/MomentumSmearingFunctions.hh \ Tools/ParticleSmearingFunctions.hh \ Tools/JetSmearingFunctions.hh \ Tools/TypeTraits.hh \ Tools/Utils.hh \ Tools/Nsubjettiness/AxesDefinition.hh \ Tools/Nsubjettiness/ExtraRecombiners.hh \ Tools/Nsubjettiness/MeasureDefinition.hh \ Tools/Nsubjettiness/Njettiness.hh \ Tools/Nsubjettiness/NjettinessPlugin.hh \ Tools/Nsubjettiness/Nsubjettiness.hh \ Tools/Nsubjettiness/TauComponents.hh \ Tools/Nsubjettiness/XConePlugin.hh nobase_dist_noinst_HEADERS += \ Tools/osdir.hh ## Maths nobase_pkginclude_HEADERS += \ Math/Matrices.hh \ Math/Vector3.hh \ Math/VectorN.hh \ Math/MatrixN.hh \ Math/MatrixDiag.hh \ Math/MathHeader.hh \ Math/Vectors.hh \ Math/LorentzTrans.hh \ Math/Matrix3.hh \ Math/MathUtils.hh \ Math/Vector4.hh \ Math/Math.hh \ Math/Units.hh \ Math/Constants.hh \ Math/eigen/util.h \ Math/eigen/regressioninternal.h \ Math/eigen/regression.h \ Math/eigen/vector.h \ Math/eigen/ludecompositionbase.h \ Math/eigen/ludecomposition.h \ Math/eigen/linearsolver.h \ Math/eigen/linearsolverbase.h \ Math/eigen/matrix.h \ Math/eigen/vectorbase.h \ Math/eigen/projective.h \ Math/eigen/matrixbase.h diff --git a/include/Rivet/Projections/EventMixingFinalState.hh b/include/Rivet/Projections/EventMixingFinalState.hh new file mode 100644 --- /dev/null +++ b/include/Rivet/Projections/EventMixingFinalState.hh @@ -0,0 +1,99 @@ +// -*- 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 { + 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."); + + // 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(); + } + + /// 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 + int compare(const Projection& p) const { + return mkNamedPCmp(p,"FS"); + } + + 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; + }; +} +#endif