diff --git a/analyses/pluginCMS/CMS_2016_I1486238.cc b/analyses/pluginCMS/CMS_2016_I1486238.cc --- a/analyses/pluginCMS/CMS_2016_I1486238.cc +++ b/analyses/pluginCMS/CMS_2016_I1486238.cc @@ -1,124 +1,126 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" -#include "Rivet/Projections/InitialQuarks.hh" #include "Rivet/Projections/UnstableFinalState.hh" #include "Rivet/Projections/FastJets.hh" +#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT +#include "Rivet/Projections/InitialQuarks.hh" + namespace Rivet { /// Studies of 2 b-jet + 2 jet production in proton-proton collisions at 7 TeV class CMS_2016_I1486238 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2016_I1486238); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { FastJets akt(FinalState(), FastJets::ANTIKT, 0.5); declare(akt, "antikT"); book(_h_Deltaphi_newway, 1,1,1); book(_h_deltaphiafterlight, 9,1,1); book(_h_SumPLight, 5,1,1); book(_h_LeadingBJetpt, 11,1,1); book(_h_SubleadingBJetpt, 15,1,1); book(_h_LeadingLightJetpt, 13,1,1); book(_h_SubleadingLightJetpt, 17,1,1); book(_h_LeadingBJeteta, 10,1,1); book(_h_SubleadingBJeteta, 14,1,1); book(_h_LeadingLightJeteta, 12,1,1); book(_h_SubleadingLightJeteta, 16,1,1); } /// Perform the per-event analysis void analyze(const Event& event) { const Jets& jets = apply(event, "antikT").jetsByPt(Cuts::absrap < 4.7 && Cuts::pT > 20*GeV); if (jets.size() < 4) vetoEvent; // Initial quarks /// @note Quark-level tagging... Particles bquarks; for (const GenParticle* p : particles(event.genEvent())) { if (abs(p->pdg_id()) == PID::BQUARK) bquarks += Particle(p); } Jets bjets, ljets; for (const Jet& j : jets) { const bool btag = any(bquarks, deltaRLess(j, 0.3)); // for (const Particle& b : bquarks) if (deltaR(j, b) < 0.3) btag = true; (btag && j.abseta() < 2.4 ? bjets : ljets).push_back(j); } // Fill histograms const double weight = 1.0; if (bjets.size() >= 2 && ljets.size() >= 2) { _h_LeadingBJetpt->fill(bjets[0].pT()/GeV, weight); _h_SubleadingBJetpt->fill(bjets[1].pT()/GeV, weight); _h_LeadingLightJetpt->fill(ljets[0].pT()/GeV, weight); _h_SubleadingLightJetpt->fill(ljets[1].pT()/GeV, weight); // _h_LeadingBJeteta->fill(bjets[0].eta(), weight); _h_SubleadingBJeteta->fill(bjets[1].eta(), weight); _h_LeadingLightJeteta->fill(ljets[0].eta(), weight); _h_SubleadingLightJeteta->fill(ljets[1].eta(), weight); const double lightdphi = deltaPhi(ljets[0], ljets[1]); _h_deltaphiafterlight->fill(lightdphi, weight); const double vecsumlightjets = sqrt(sqr(ljets[0].px()+ljets[1].px()) + sqr(ljets[0].py()+ljets[1].py())); //< @todo Just (lj0+lj1).pT()? Or use add_quad const double term2 = vecsumlightjets/(sqrt(sqr(ljets[0].px()) + sqr(ljets[0].py())) + sqrt(sqr(ljets[1].px()) + sqr(ljets[1].py()))); //< @todo lj0.pT() + lj1.pT()? Or add_quad _h_SumPLight->fill(term2, weight); const double pxBsyst2 = bjets[0].px()+bjets[1].px(); // @todo (bj0+bj1).px() const double pyBsyst2 = bjets[0].py()+bjets[1].py(); // @todo (bj0+bj1).py() const double pxJetssyst2 = ljets[0].px()+ljets[1].px(); // @todo (lj0+lj1).px() const double pyJetssyst2 = ljets[0].py()+ljets[1].py(); // @todo (lj0+lj1).py() const double modulusB2 = sqrt(sqr(pxBsyst2)+sqr(pyBsyst2)); //< @todo add_quad const double modulusJets2 = sqrt(sqr(pxJetssyst2)+sqr(pyJetssyst2)); //< @todo add_quad const double cosphiBsyst2 = pxBsyst2/modulusB2; const double cosphiJetssyst2 = pxJetssyst2/modulusJets2; const double phiBsyst2 = ((pyBsyst2 > 0) ? 1 : -1) * acos(cosphiBsyst2); //< @todo sign(pyBsyst2) const double phiJetssyst2 = sign(pyJetssyst2) * acos(cosphiJetssyst2); const double Dphi2 = deltaPhi(phiBsyst2, phiJetssyst2); _h_Deltaphi_newway->fill(Dphi2,weight); } } /// Normalise histograms etc., after the run void finalize() { const double invlumi = crossSection()/picobarn/sumOfWeights(); normalize({_h_SumPLight, _h_deltaphiafterlight, _h_Deltaphi_newway}); scale({_h_LeadingLightJetpt, _h_SubleadingLightJetpt, _h_LeadingBJetpt, _h_SubleadingBJetpt}, invlumi); scale({_h_LeadingLightJeteta, _h_SubleadingLightJeteta, _h_LeadingBJeteta, _h_SubleadingBJeteta}, invlumi); } //@} private: /// @name Histograms //@{ Histo1DPtr _h_deltaphiafterlight, _h_Deltaphi_newway, _h_SumPLight; Histo1DPtr _h_LeadingBJetpt, _h_SubleadingBJetpt, _h_LeadingLightJetpt, _h_SubleadingLightJetpt; Histo1DPtr _h_LeadingBJeteta, _h_SubleadingBJeteta, _h_LeadingLightJeteta, _h_SubleadingLightJeteta; }; // Hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2016_I1486238); } diff --git a/analyses/pluginLEP/DELPHI_2000_S4328825.cc b/analyses/pluginLEP/DELPHI_2000_S4328825.cc --- a/analyses/pluginLEP/DELPHI_2000_S4328825.cc +++ b/analyses/pluginLEP/DELPHI_2000_S4328825.cc @@ -1,143 +1,145 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/Beam.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/Sphericity.hh" #include "Rivet/Projections/Thrust.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/ParisiTensor.hh" #include "Rivet/Projections/Hemispheres.hh" +#include + +#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT #include "Rivet/Projections/InitialQuarks.hh" -#include namespace Rivet { /// @brief OPAL multiplicities at various energies /// @author Peter Richardson class DELPHI_2000_S4328825 : public Analysis { public: /// Constructor DELPHI_2000_S4328825() : Analysis("DELPHI_2000_S4328825") {} /// @name Analysis methods //@{ void init() { // Projections declare(Beam(), "Beams"); declare(ChargedFinalState(), "CFS"); declare(InitialQuarks(), "IQF"); book(_weightedTotalChargedPartNumLight,"weight_totalch_light"); book(_weightedTotalChargedPartNumCharm,"weight_totalch_charm"); book(_weightedTotalChargedPartNumBottom,"weight_totalch_bottom"); book(_weightLight,"weight_light"); book(_weightCharm,"weight_charm"); book(_weightBottom,"weight_bottom"); book(h_bottom, 1, 1, 1); book(h_charm, 1, 1, 2); book(h_light, 1, 1, 3); book(h_diff, 1, 1, 4); // bottom minus light } void analyze(const Event& event) { // Even if we only generate hadronic events, we still need a cut on numCharged >= 2. const FinalState& cfs = apply(event, "CFS"); if (cfs.size() < 2) vetoEvent; int flavour = 0; const InitialQuarks& iqf = apply(event, "IQF"); // If we only have two quarks (qqbar), just take the flavour. // If we have more than two quarks, look for the highest energetic q-qbar pair. if (iqf.particles().size() == 2) { flavour = iqf.particles().front().abspid(); } else { map quarkmap; for (const Particle& p : iqf.particles()) { if (quarkmap[p.pid()] < p.E()) { quarkmap[p.pid()] = p.E(); } } double maxenergy = 0.; for (int i = 1; i <= 5; ++i) { if (quarkmap[i]+quarkmap[-i] > maxenergy) { flavour = i; } } } const size_t numParticles = cfs.particles().size(); switch (flavour) { case 1: case 2: case 3: _weightLight->fill(); _weightedTotalChargedPartNumLight ->fill(numParticles); break; case 4: _weightCharm->fill(); _weightedTotalChargedPartNumCharm ->fill(numParticles); break; case 5: _weightBottom->fill(); _weightedTotalChargedPartNumBottom->fill(numParticles); break; } } void finalize() { Histo1D temphisto(refData(1, 1, 1)); const double avgNumPartsBottom = dbl(*_weightedTotalChargedPartNumBottom / *_weightBottom); const double avgNumPartsCharm = dbl(*_weightedTotalChargedPartNumCharm / *_weightCharm); const double avgNumPartsLight = dbl(*_weightedTotalChargedPartNumLight / *_weightLight); for (size_t b = 0; b < temphisto.numBins(); b++) { const double x = temphisto.bin(b).xMid(); const double ex = temphisto.bin(b).xWidth()/2.; if (inRange(sqrtS()/GeV, x-ex, x+ex)) { // @TODO: Fix y-error: h_bottom->addPoint(x, avgNumPartsBottom, ex, 0.); h_charm->addPoint(x, avgNumPartsCharm, ex, 0.); h_light->addPoint(x, avgNumPartsLight, ex, 0.); h_diff->addPoint(x, avgNumPartsBottom-avgNumPartsLight, ex, 0.); } } } //@} private: Scatter2DPtr h_bottom, h_charm, h_light, h_diff; /// @name Multiplicities //@{ CounterPtr _weightedTotalChargedPartNumLight; CounterPtr _weightedTotalChargedPartNumCharm; CounterPtr _weightedTotalChargedPartNumBottom; //@} /// @name Weights //@{ CounterPtr _weightLight; CounterPtr _weightCharm; CounterPtr _weightBottom; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(DELPHI_2000_S4328825); } diff --git a/analyses/pluginLEP/L3_2004_I652683.cc b/analyses/pluginLEP/L3_2004_I652683.cc --- a/analyses/pluginLEP/L3_2004_I652683.cc +++ b/analyses/pluginLEP/L3_2004_I652683.cc @@ -1,210 +1,212 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/Beam.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" -#include "Rivet/Projections/InitialQuarks.hh" #include "Rivet/Projections/Thrust.hh" #include "Rivet/Projections/ParisiTensor.hh" #include "Rivet/Projections/Hemispheres.hh" +#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT +#include "Rivet/Projections/InitialQuarks.hh" + namespace Rivet { /// Jet rates and event shapes at LEP I+II class L3_2004_I652683 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(L3_2004_I652683); // L3_2004_I652683() : Analysis("L3_2004_I652683") // { } /// Book histograms and initialise projections before the run void init() { // Projections to use const FinalState FS; declare(FS, "FS"); declare(Beam(), "beams"); const ChargedFinalState CFS; declare(CFS, "CFS"); const Thrust thrust(FS); declare(thrust, "thrust"); declare(ParisiTensor(FS), "Parisi"); declare(Hemispheres(thrust), "Hemispheres"); declare(InitialQuarks(), "initialquarks"); // Book the histograms book(_h_Thrust_udsc , 47, 1, 1); book(_h_Thrust_bottom , 47, 1, 2); book(_h_heavyJetmass_udsc , 48, 1, 1); book(_h_heavyJetmass_bottom , 48, 1, 2); book(_h_totalJetbroad_udsc , 49, 1, 1); book(_h_totalJetbroad_bottom , 49, 1, 2); book(_h_wideJetbroad_udsc , 50, 1, 1); book(_h_wideJetbroad_bottom , 50, 1, 2); book(_h_Cparameter_udsc , 51, 1, 1); book(_h_Cparameter_bottom , 51, 1, 2); book(_h_Dparameter_udsc , 52, 1, 1); book(_h_Dparameter_bottom , 52, 1, 2); book(_h_Ncharged , 59, 1, 1); book(_h_Ncharged_udsc , 59, 1, 2); book(_h_Ncharged_bottom , 59, 1, 3); book(_h_scaledMomentum , 65, 1, 1); book(_h_scaledMomentum_udsc , 65, 1, 2); book(_h_scaledMomentum_bottom, 65, 1, 3); book(_sumW_udsc, "sumW_udsc"); book(_sumW_b, "sumW_b"); book(_sumW_ch, "sumW_ch"); book(_sumW_ch_udsc, "sumW_ch_udsc"); book(_sumW_ch_b, "sumW_ch_b"); } /// Perform the per-event analysis void analyze(const Event& event) { // Get beam average momentum const ParticlePair& beams = apply(event, "beams").beams(); const double beamMomentum = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0; // InitialQuarks projection to have udsc events separated from b events /// @todo Yuck!!! Eliminate when possible... int flavour = 0; const InitialQuarks& iqf = apply(event, "initialquarks"); Particles quarks; if ( iqf.particles().size() == 2 ) { flavour = iqf.particles().front().abspid(); quarks = iqf.particles(); } else { map quarkmap; for (const Particle& p : iqf.particles()) { if (quarkmap.find(p.pid()) == quarkmap.end()) quarkmap[p.pid()] = p; else if (quarkmap[p.pid()].E() < p.E()) quarkmap[p.pid()] = p; } double max_energy = 0.; for (int i = 1; i <= 5; ++i) { double energy = 0.; if (quarkmap.find(i) != quarkmap.end()) energy += quarkmap[ i].E(); if (quarkmap.find(-i) != quarkmap.end()) energy += quarkmap[-i].E(); if (energy > max_energy) flavour = i; } if (quarkmap.find(flavour) != quarkmap.end()) quarks.push_back(quarkmap[flavour]); if (quarkmap.find(-flavour) != quarkmap.end()) quarks.push_back(quarkmap[-flavour]); } // Flavour label /// @todo Change to a bool? const int iflav = (flavour == PID::DQUARK || flavour == PID::UQUARK || flavour == PID::SQUARK || flavour == PID::CQUARK) ? 1 : (flavour == PID::BQUARK) ? 5 : 0; // Update weight sums if (iflav == 1) { _sumW_udsc->fill(); } else if (iflav == 5) { _sumW_b->fill(); } _sumW_ch->fill(); // Charged multiplicity const FinalState& cfs = applyProjection(event, "CFS"); _h_Ncharged->fill(cfs.size()); if (iflav == 1) { _sumW_ch_udsc->fill(); _h_Ncharged_udsc->fill(cfs.size()); } else if (iflav == 5) { _sumW_ch_b->fill(); _h_Ncharged_bottom->fill(cfs.size()); } // Scaled momentum const Particles& chparticles = cfs.particlesByPt(); for (const Particle& p : chparticles) { const Vector3 momentum3 = p.p3(); const double mom = momentum3.mod(); const double scaledMom = mom/beamMomentum; const double logScaledMom = std::log(scaledMom); _h_scaledMomentum->fill(-logScaledMom); if (iflav == 1) { _h_scaledMomentum_udsc->fill(-logScaledMom); } else if (iflav == 5) { _h_scaledMomentum_bottom->fill(-logScaledMom); } } // Thrust const Thrust& thrust = applyProjection(event, "thrust"); if (iflav == 1) { _h_Thrust_udsc->fill(thrust.thrust()); } else if (iflav == 5) { _h_Thrust_bottom->fill(thrust.thrust()); } // C and D Parisi parameters const ParisiTensor& parisi = applyProjection(event, "Parisi"); if (iflav == 1) { _h_Cparameter_udsc->fill(parisi.C()); _h_Dparameter_udsc->fill(parisi.D()); } else if (iflav == 5) { _h_Cparameter_bottom->fill(parisi.C()); _h_Dparameter_bottom->fill(parisi.D()); } // The hemisphere variables const Hemispheres& hemisphere = applyProjection(event, "Hemispheres"); if (iflav == 1) { _h_heavyJetmass_udsc->fill(hemisphere.scaledM2high()); _h_totalJetbroad_udsc->fill(hemisphere.Bsum()); _h_wideJetbroad_udsc->fill(hemisphere.Bmax()); } else if (iflav == 5) { _h_heavyJetmass_bottom->fill(hemisphere.scaledM2high()); _h_totalJetbroad_bottom->fill(hemisphere.Bsum()); _h_wideJetbroad_bottom->fill(hemisphere.Bmax()); } } /// Normalise histograms etc., after the run void finalize() { scale({_h_Thrust_udsc, _h_heavyJetmass_udsc, _h_totalJetbroad_udsc, _h_wideJetbroad_udsc, _h_Cparameter_udsc, _h_Dparameter_udsc}, 1/ *_sumW_udsc); scale({_h_Thrust_bottom, _h_heavyJetmass_bottom, _h_totalJetbroad_bottom, _h_wideJetbroad_bottom, _h_Cparameter_bottom, _h_Dparameter_bottom}, 1./ *_sumW_b); scale(_h_Ncharged, 2/ *_sumW_ch); scale(_h_Ncharged_udsc, 2/ *_sumW_ch_udsc); scale(_h_Ncharged_bottom, 2/ *_sumW_ch_b); scale(_h_scaledMomentum, 1/ *_sumW_ch); scale(_h_scaledMomentum_udsc, 1/ *_sumW_ch_udsc); scale(_h_scaledMomentum_bottom, 1/ *_sumW_ch_b); } /// Weight counters CounterPtr _sumW_udsc, _sumW_b, _sumW_ch, _sumW_ch_udsc, _sumW_ch_b; /// @name Histograms //@{ Histo1DPtr _h_Thrust_udsc, _h_Thrust_bottom; Histo1DPtr _h_heavyJetmass_udsc, _h_heavyJetmass_bottom; Histo1DPtr _h_totalJetbroad_udsc, _h_totalJetbroad_bottom; Histo1DPtr _h_wideJetbroad_udsc, _h_wideJetbroad_bottom; Histo1DPtr _h_Cparameter_udsc, _h_Cparameter_bottom; Histo1DPtr _h_Dparameter_udsc, _h_Dparameter_bottom; Histo1DPtr _h_Ncharged, _h_Ncharged_udsc, _h_Ncharged_bottom; Histo1DPtr _h_scaledMomentum, _h_scaledMomentum_udsc, _h_scaledMomentum_bottom; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(L3_2004_I652683); } diff --git a/analyses/pluginLEP/OPAL_1998_S3780481.cc b/analyses/pluginLEP/OPAL_1998_S3780481.cc --- a/analyses/pluginLEP/OPAL_1998_S3780481.cc +++ b/analyses/pluginLEP/OPAL_1998_S3780481.cc @@ -1,193 +1,195 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/Beam.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" + +#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT #include "Rivet/Projections/InitialQuarks.hh" namespace Rivet { /// @brief OPAL flavour-dependent fragmentation paper /// @author Hendrik Hoeth class OPAL_1998_S3780481 : public Analysis { public: /// Constructor OPAL_1998_S3780481() : Analysis("OPAL_1998_S3780481") { } /// @name Analysis methods //@{ void analyze(const Event& e) { // First, veto on leptonic events by requiring at least 4 charged FS particles const FinalState& fs = apply(e, "FS"); const size_t numParticles = fs.particles().size(); // Even if we only generate hadronic events, we still need a cut on numCharged >= 2. if (numParticles < 2) { MSG_DEBUG("Failed ncharged cut"); vetoEvent; } MSG_DEBUG("Passed ncharged cut"); _weightedTotalPartNum->fill(numParticles); // Get beams and average beam momentum const ParticlePair& beams = apply(e, "Beams").beams(); const double meanBeamMom = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0; MSG_DEBUG("Avg beam momentum = " << meanBeamMom); int flavour = 0; const InitialQuarks& iqf = apply(e, "IQF"); // If we only have two quarks (qqbar), just take the flavour. // If we have more than two quarks, look for the highest energetic q-qbar pair. /// @todo Yuck... does this *really* have to be quark-based?!? if (iqf.particles().size() == 2) { flavour = iqf.particles().front().abspid(); } else { map quarkmap; for (const Particle& p : iqf.particles()) { if (quarkmap[p.pid()] < p.E()) { quarkmap[p.pid()] = p.E(); } } double maxenergy = 0.; for (int i = 1; i <= 5; ++i) { if (quarkmap[i]+quarkmap[-i] > maxenergy) { flavour = i; } } } switch (flavour) { case 1: case 2: case 3: _SumOfudsWeights->fill(); break; case 4: _SumOfcWeights->fill(); break; case 5: _SumOfbWeights->fill(); break; } for (const Particle& p : fs.particles()) { const double xp = p.p3().mod()/meanBeamMom; const double logxp = -std::log(xp); _histXpall->fill(xp); _histLogXpall->fill(logxp); _histMultiChargedall->fill(_histMultiChargedall->bin(0).xMid()); switch (flavour) { /// @todo Use PDG code enums case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _histXpuds->fill(xp); _histLogXpuds->fill(logxp); _histMultiChargeduds->fill(_histMultiChargeduds->bin(0).xMid()); break; case PID::CQUARK: _histXpc->fill(xp); _histLogXpc->fill(logxp); _histMultiChargedc->fill(_histMultiChargedc->bin(0).xMid()); break; case PID::BQUARK: _histXpb->fill(xp); _histLogXpb->fill(logxp); _histMultiChargedb->fill(_histMultiChargedb->bin(0).xMid()); break; } } } void init() { // Projections declare(Beam(), "Beams"); declare(ChargedFinalState(), "FS"); declare(InitialQuarks(), "IQF"); // Book histos book(_histXpuds ,1, 1, 1); book(_histXpc ,2, 1, 1); book(_histXpb ,3, 1, 1); book(_histXpall ,4, 1, 1); book(_histLogXpuds ,5, 1, 1); book(_histLogXpc ,6, 1, 1); book(_histLogXpb ,7, 1, 1); book(_histLogXpall ,8, 1, 1); book(_histMultiChargeduds ,9, 1, 1); book(_histMultiChargedc ,9, 1, 2); book(_histMultiChargedb ,9, 1, 3); book(_histMultiChargedall ,9, 1, 4); // Counters book(_weightedTotalPartNum, "TotalPartNum"); book(_SumOfudsWeights, "udsWeights"); book(_SumOfcWeights, "cWeights"); book(_SumOfbWeights, "bWeights"); } /// Finalize void finalize() { const double avgNumParts = dbl(*_weightedTotalPartNum) / sumOfWeights(); normalize(_histXpuds , avgNumParts); normalize(_histXpc , avgNumParts); normalize(_histXpb , avgNumParts); normalize(_histXpall , avgNumParts); normalize(_histLogXpuds , avgNumParts); normalize(_histLogXpc , avgNumParts); normalize(_histLogXpb , avgNumParts); normalize(_histLogXpall , avgNumParts); scale(_histMultiChargeduds, 1.0/ *_SumOfudsWeights); scale(_histMultiChargedc , 1.0/ *_SumOfcWeights); scale(_histMultiChargedb , 1.0/ *_SumOfbWeights); scale(_histMultiChargedall, 1.0/sumOfWeights()); } //@} private: /// Store the weighted sums of numbers of charged / charged+neutral /// particles - used to calculate average number of particles for the /// inclusive single particle distributions' normalisations. CounterPtr _weightedTotalPartNum; CounterPtr _SumOfudsWeights; CounterPtr _SumOfcWeights; CounterPtr _SumOfbWeights; Histo1DPtr _histXpuds; Histo1DPtr _histXpc; Histo1DPtr _histXpb; Histo1DPtr _histXpall; Histo1DPtr _histLogXpuds; Histo1DPtr _histLogXpc; Histo1DPtr _histLogXpb; Histo1DPtr _histLogXpall; Histo1DPtr _histMultiChargeduds; Histo1DPtr _histMultiChargedc; Histo1DPtr _histMultiChargedb; Histo1DPtr _histMultiChargedall; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(OPAL_1998_S3780481); } diff --git a/analyses/pluginLEP/OPAL_2002_S5361494.cc b/analyses/pluginLEP/OPAL_2002_S5361494.cc --- a/analyses/pluginLEP/OPAL_2002_S5361494.cc +++ b/analyses/pluginLEP/OPAL_2002_S5361494.cc @@ -1,145 +1,147 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/Beam.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/Sphericity.hh" #include "Rivet/Projections/Thrust.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/ParisiTensor.hh" #include "Rivet/Projections/Hemispheres.hh" +#include + +#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT #include "Rivet/Projections/InitialQuarks.hh" -#include namespace Rivet { /// @brief OPAL multiplicities at various energies /// @author Peter Richardson class OPAL_2002_S5361494 : public Analysis { public: /// Constructor OPAL_2002_S5361494() : Analysis("OPAL_2002_S5361494") {} /// @name Analysis methods //@{ void init() { // Projections declare(Beam(), "Beams"); declare(ChargedFinalState(), "CFS"); declare(InitialQuarks(), "IQF"); book(h_bottom, 1, 1, 1); book(h_charm, 1, 1, 2); book(h_light, 1, 1, 3); book(h_diff, 1, 1, 4); // bottom minus light book(_weightedTotalChargedPartNumLight, "TotalChargedPartNumLight"); book(_weightedTotalChargedPartNumCharm, "TotalChargedPartNumCharm"); book(_weightedTotalChargedPartNumBottom, "TotalChargedPartNumBottom"); book(_weightLight, "Light"); book(_weightCharm, "Charm"); book(_weightBottom, "Bottom"); } void analyze(const Event& event) { // Even if we only generate hadronic events, we still need a cut on numCharged >= 2. const FinalState& cfs = apply(event, "CFS"); if (cfs.size() < 2) vetoEvent; int flavour = 0; const InitialQuarks& iqf = apply(event, "IQF"); // If we only have two quarks (qqbar), just take the flavour. // If we have more than two quarks, look for the highest energetic q-qbar pair. if (iqf.particles().size() == 2) { flavour = iqf.particles().front().abspid(); } else { map quarkmap; for (const Particle& p : iqf.particles()) { if (quarkmap[p.pid()] < p.E()) { quarkmap[p.pid()] = p.E(); } } double maxenergy = 0.; for (int i = 1; i <= 5; ++i) { if (quarkmap[i]+quarkmap[-i] > maxenergy) { flavour = i; } } } const size_t numParticles = cfs.particles().size(); switch (flavour) { case 1: case 2: case 3: _weightLight ->fill(); _weightedTotalChargedPartNumLight ->fill(numParticles); break; case 4: _weightCharm ->fill(); _weightedTotalChargedPartNumCharm ->fill(numParticles); break; case 5: _weightBottom->fill(); _weightedTotalChargedPartNumBottom->fill(numParticles); break; } } void finalize() { Histo1D temphisto(refData(1, 1, 1)); const double avgNumPartsBottom = _weightBottom->val() != 0. ? dbl(*_weightedTotalChargedPartNumBottom / *_weightBottom) : 0.; const double avgNumPartsCharm = _weightCharm->val() != 0. ? dbl(*_weightedTotalChargedPartNumCharm / *_weightCharm ) : 0.; const double avgNumPartsLight = _weightLight->val() != 0. ? dbl(*_weightedTotalChargedPartNumLight / *_weightLight ) : 0.; for (size_t b = 0; b < temphisto.numBins(); b++) { const double x = temphisto.bin(b).xMid(); const double ex = temphisto.bin(b).xWidth()/2.; if (inRange(sqrtS()/GeV, x-ex, x+ex)) { // @TODO: Fix y-error: h_bottom->addPoint(x, avgNumPartsBottom, ex, 0.); h_charm->addPoint(x, avgNumPartsCharm, ex, 0.); h_light->addPoint(x, avgNumPartsLight, ex, 0.); h_diff->addPoint(x, avgNumPartsBottom-avgNumPartsLight, ex, 0.); } } } //@} private: Scatter2DPtr h_bottom; Scatter2DPtr h_charm ; Scatter2DPtr h_light ; Scatter2DPtr h_diff ; /// @name Multiplicities //@{ CounterPtr _weightedTotalChargedPartNumLight; CounterPtr _weightedTotalChargedPartNumCharm; CounterPtr _weightedTotalChargedPartNumBottom; //@} /// @name Weights //@{ CounterPtr _weightLight; CounterPtr _weightCharm; CounterPtr _weightBottom; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(OPAL_2002_S5361494); } diff --git a/analyses/pluginLEP/SLD_1996_S3398250.cc b/analyses/pluginLEP/SLD_1996_S3398250.cc --- a/analyses/pluginLEP/SLD_1996_S3398250.cc +++ b/analyses/pluginLEP/SLD_1996_S3398250.cc @@ -1,139 +1,141 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/Beam.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/Sphericity.hh" #include "Rivet/Projections/Thrust.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/ParisiTensor.hh" #include "Rivet/Projections/Hemispheres.hh" +#include + +#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT #include "Rivet/Projections/InitialQuarks.hh" -#include namespace Rivet { /// @brief SLD multiplicities at mZ /// @author Peter Richardson class SLD_1996_S3398250 : public Analysis { public: /// Constructor SLD_1996_S3398250() : Analysis("SLD_1996_S3398250") {} /// @name Analysis methods //@{ void init() { // Projections declare(Beam(), "Beams"); declare(ChargedFinalState(), "CFS"); declare(InitialQuarks(), "IQF"); book(_h_bottom ,1, 1, 1); book(_h_charm ,2, 1, 1); book(_h_light ,3, 1, 1); book(_weightLight, "weightLight"); book(_weightCharm, "weightCharm"); book(_weightBottom, "weightBottom"); book(scatter_c, 4,1,1); book(scatter_b, 5,1,1); } void analyze(const Event& event) { // Even if we only generate hadronic events, we still need a cut on numCharged >= 2. const FinalState& cfs = apply(event, "CFS"); if (cfs.size() < 2) vetoEvent; int flavour = 0; const InitialQuarks& iqf = apply(event, "IQF"); // If we only have two quarks (qqbar), just take the flavour. // If we have more than two quarks, look for the highest energetic q-qbar pair. if (iqf.particles().size() == 2) { flavour = iqf.particles().front().abspid(); } else { map quarkmap; for (const Particle& p : iqf.particles()) { if (quarkmap[p.pid()] < p.E()) { quarkmap[p.pid()] = p.E(); } } double maxenergy = 0.; for (int i = 1; i <= 5; ++i) { if (quarkmap[i]+quarkmap[-i] > maxenergy) { flavour = i; } } } const size_t numParticles = cfs.particles().size(); switch (flavour) { case 1: case 2: case 3: _weightLight ->fill(); _h_light->fillBin(0, numParticles); break; case 4: _weightCharm ->fill(); _h_charm->fillBin(0, numParticles); break; case 5: _weightBottom->fill(); _h_bottom->fillBin(0, numParticles); break; } } void multiplicity_subtract(const Histo1DPtr first, const Histo1DPtr second, Scatter2DPtr & scatter) { const double x = first->bin(0).xMid(); const double ex = first->bin(0).xWidth()/2.; const double y = first->bin(0).area() - second->bin(0).area(); const double ey = sqrt(sqr(first->bin(0).areaErr()) + sqr(second->bin(0).areaErr())); scatter->addPoint(x, y, ex, ey); } void finalize() { if (_weightBottom->val() != 0) scale(_h_bottom, 1./ *_weightBottom); if (_weightCharm->val() != 0) scale(_h_charm, 1./ *_weightCharm ); if (_weightLight->val() != 0) scale(_h_light, 1./ *_weightLight ); multiplicity_subtract(_h_charm, _h_light, scatter_c); multiplicity_subtract(_h_bottom, _h_light, scatter_b); } //@} private: Scatter2DPtr scatter_c, scatter_b; /// @name Weights //@{ CounterPtr _weightLight; CounterPtr _weightCharm; CounterPtr _weightBottom; //@} Histo1DPtr _h_bottom; Histo1DPtr _h_charm; Histo1DPtr _h_light; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(SLD_1996_S3398250); } diff --git a/analyses/pluginLEP/SLD_1999_S3743934.cc b/analyses/pluginLEP/SLD_1999_S3743934.cc --- a/analyses/pluginLEP/SLD_1999_S3743934.cc +++ b/analyses/pluginLEP/SLD_1999_S3743934.cc @@ -1,736 +1,738 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/Beam.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/UnstableFinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" +#include "Rivet/Projections/Thrust.hh" + +#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT #include "Rivet/Projections/InitialQuarks.hh" -#include "Rivet/Projections/Thrust.hh" namespace Rivet { /// @brief SLD flavour-dependent fragmentation paper /// @author Peter Richardson class SLD_1999_S3743934 : public Analysis { public: /// Constructor SLD_1999_S3743934() : Analysis("SLD_1999_S3743934"), _multPiPlus(4),_multKPlus(4),_multK0(4), _multKStar0(4),_multPhi(4), _multProton(4),_multLambda(4) { } /// @name Analysis methods //@{ void analyze(const Event& e) { // First, veto on leptonic events by requiring at least 4 charged FS particles const FinalState& fs = apply(e, "FS"); const size_t numParticles = fs.particles().size(); // Even if we only generate hadronic events, we still need a cut on numCharged >= 2. if (numParticles < 2) { MSG_DEBUG("Failed ncharged cut"); vetoEvent; } MSG_DEBUG("Passed ncharged cut"); // Get beams and average beam momentum const ParticlePair& beams = apply(e, "Beams").beams(); const double meanBeamMom = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0; MSG_DEBUG("Avg beam momentum = " << meanBeamMom); int flavour = 0; const InitialQuarks& iqf = apply(e, "IQF"); // If we only have two quarks (qqbar), just take the flavour. // If we have more than two quarks, look for the highest energetic q-qbar pair. /// @todo Can we make this based on hadron flavour instead? Particles quarks; if (iqf.particles().size() == 2) { flavour = iqf.particles().front().abspid(); quarks = iqf.particles(); } else { map quarkmap; for (const Particle& p : iqf.particles()) { if (quarkmap.find(p.pid()) == quarkmap.end()) quarkmap[p.pid()] = p; else if (quarkmap[p.pid()].E() < p.E()) quarkmap[p.pid()] = p; } double maxenergy = 0.; for (int i = 1; i <= 5; ++i) { double energy(0.); if (quarkmap.find( i) != quarkmap.end()) energy += quarkmap[ i].E(); if (quarkmap.find(-i) != quarkmap.end()) energy += quarkmap[-i].E(); if (energy > maxenergy) flavour = i; } if (quarkmap.find(flavour) != quarkmap.end()) quarks.push_back(quarkmap[flavour]); if (quarkmap.find(-flavour) != quarkmap.end()) quarks.push_back(quarkmap[-flavour]); } switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _SumOfudsWeights->fill(); break; case PID::CQUARK: _SumOfcWeights->fill(); break; case PID::BQUARK: _SumOfbWeights->fill(); break; } // thrust axis for projections Vector3 axis = apply(e, "Thrust").thrustAxis(); double dot(0.); if (!quarks.empty()) { dot = quarks[0].p3().dot(axis); if (quarks[0].pid() < 0) dot *= -1; } for (const Particle& p : fs.particles()) { const double xp = p.p3().mod()/meanBeamMom; // if in quark or antiquark hemisphere bool quark = p.p3().dot(axis)*dot > 0.; _h_XpChargedN->fill(xp); _temp_XpChargedN1->fill(xp); _temp_XpChargedN2->fill(xp); _temp_XpChargedN3->fill(xp); int id = p.abspid(); // charged pions if (id == PID::PIPLUS) { _h_XpPiPlusN->fill(xp); _multPiPlus[0]->fill(); switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _multPiPlus[1]->fill(); _h_XpPiPlusLight->fill(xp); if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 )) _h_RPiPlus->fill(xp); else _h_RPiMinus->fill(xp); break; case PID::CQUARK: _multPiPlus[2]->fill(); _h_XpPiPlusCharm->fill(xp); break; case PID::BQUARK: _multPiPlus[3]->fill(); _h_XpPiPlusBottom->fill(xp); break; } } else if (id == PID::KPLUS) { _h_XpKPlusN->fill(xp); _multKPlus[0]->fill(); switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _multKPlus[1]->fill(); _temp_XpKPlusLight->fill(xp); _h_XpKPlusLight->fill(xp); if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 )) _h_RKPlus->fill(xp); else _h_RKMinus->fill(xp); break; break; case PID::CQUARK: _multKPlus[2]->fill(); _h_XpKPlusCharm->fill(xp); _temp_XpKPlusCharm->fill(xp); break; case PID::BQUARK: _multKPlus[3]->fill(); _h_XpKPlusBottom->fill(xp); break; } } else if (id == PID::PROTON) { _h_XpProtonN->fill(xp); _multProton[0]->fill(); switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _multProton[1]->fill(); _temp_XpProtonLight->fill(xp); _h_XpProtonLight->fill(xp); if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 )) _h_RProton->fill(xp); else _h_RPBar ->fill(xp); break; break; case PID::CQUARK: _multProton[2]->fill(); _temp_XpProtonCharm->fill(xp); _h_XpProtonCharm->fill(xp); break; case PID::BQUARK: _multProton[3]->fill(); _h_XpProtonBottom->fill(xp); break; } } } const UnstableFinalState& ufs = apply(e, "UFS"); for (const Particle& p : ufs.particles()) { const double xp = p.p3().mod()/meanBeamMom; // if in quark or antiquark hemisphere bool quark = p.p3().dot(axis)*dot>0.; int id = p.abspid(); if (id == PID::LAMBDA) { _multLambda[0]->fill(); _h_XpLambdaN->fill(xp); switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _multLambda[1]->fill(); _h_XpLambdaLight->fill(xp); if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 )) _h_RLambda->fill(xp); else _h_RLBar ->fill(xp); break; case PID::CQUARK: _multLambda[2]->fill(); _h_XpLambdaCharm->fill(xp); break; case PID::BQUARK: _multLambda[3]->fill(); _h_XpLambdaBottom->fill(xp); break; } } else if (id == 313) { _multKStar0[0]->fill(); _h_XpKStar0N->fill(xp); switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _multKStar0[1]->fill(); _temp_XpKStar0Light->fill(xp); _h_XpKStar0Light->fill(xp); if ( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 )) _h_RKS0 ->fill(xp); else _h_RKSBar0->fill(xp); break; break; case PID::CQUARK: _multKStar0[2]->fill(); _temp_XpKStar0Charm->fill(xp); _h_XpKStar0Charm->fill(xp); break; case PID::BQUARK: _multKStar0[3]->fill(); _h_XpKStar0Bottom->fill(xp); break; } } else if (id == 333) { _multPhi[0]->fill(); _h_XpPhiN->fill(xp); switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _multPhi[1]->fill(); _h_XpPhiLight->fill(xp); break; case PID::CQUARK: _multPhi[2]->fill(); _h_XpPhiCharm->fill(xp); break; case PID::BQUARK: _multPhi[3]->fill(); _h_XpPhiBottom->fill(xp); break; } } else if (id == PID::K0S || id == PID::K0L) { _multK0[0]->fill(); _h_XpK0N->fill(xp); switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _multK0[1]->fill(); _h_XpK0Light->fill(xp); break; case PID::CQUARK: _multK0[2]->fill(); _h_XpK0Charm->fill(xp); break; case PID::BQUARK: _multK0[3]->fill(); _h_XpK0Bottom->fill(xp); break; } } } } void init() { // Projections declare(Beam(), "Beams"); declare(ChargedFinalState(), "FS"); declare(UnstableFinalState(), "UFS"); declare(InitialQuarks(), "IQF"); declare(Thrust(FinalState()), "Thrust"); book(_temp_XpChargedN1 ,"TMP/XpChargedN1", refData( 1, 1, 1)); book(_temp_XpChargedN2 ,"TMP/XpChargedN2", refData( 2, 1, 1)); book(_temp_XpChargedN3 ,"TMP/XpChargedN3", refData( 3, 1, 1)); book(_h_XpPiPlusN , 1, 1, 2); book(_h_XpKPlusN , 2, 1, 2); book(_h_XpProtonN , 3, 1, 2); book(_h_XpChargedN , 4, 1, 1); book(_h_XpK0N , 5, 1, 1); book(_h_XpLambdaN , 7, 1, 1); book(_h_XpKStar0N , 8, 1, 1); book(_h_XpPhiN , 9, 1, 1); book(_h_XpPiPlusLight ,10, 1, 1); book(_h_XpPiPlusCharm ,10, 1, 2); book(_h_XpPiPlusBottom ,10, 1, 3); book(_h_XpKPlusLight ,12, 1, 1); book(_h_XpKPlusCharm ,12, 1, 2); book(_h_XpKPlusBottom ,12, 1, 3); book(_h_XpKStar0Light ,14, 1, 1); book(_h_XpKStar0Charm ,14, 1, 2); book(_h_XpKStar0Bottom ,14, 1, 3); book(_h_XpProtonLight ,16, 1, 1); book(_h_XpProtonCharm ,16, 1, 2); book(_h_XpProtonBottom ,16, 1, 3); book(_h_XpLambdaLight ,18, 1, 1); book(_h_XpLambdaCharm ,18, 1, 2); book(_h_XpLambdaBottom ,18, 1, 3); book(_h_XpK0Light ,20, 1, 1); book(_h_XpK0Charm ,20, 1, 2); book(_h_XpK0Bottom ,20, 1, 3); book(_h_XpPhiLight ,22, 1, 1); book(_h_XpPhiCharm ,22, 1, 2); book(_h_XpPhiBottom ,22, 1, 3); book(_temp_XpKPlusCharm ,"TMP/XpKPlusCharm", refData(13, 1, 1)); book(_temp_XpKPlusLight ,"TMP/XpKPlusLight", refData(13, 1, 1)); book(_temp_XpKStar0Charm ,"TMP/XpKStar0Charm", refData(15, 1, 1)); book(_temp_XpKStar0Light ,"TMP/XpKStar0Light", refData(15, 1, 1)); book(_temp_XpProtonCharm ,"TMP/XpProtonCharm", refData(17, 1, 1)); book(_temp_XpProtonLight ,"TMP/XpProtonLight", refData(17, 1, 1)); book(_h_RPiPlus , 26, 1, 1); book(_h_RPiMinus , 26, 1, 2); book(_h_RKS0 , 28, 1, 1); book(_h_RKSBar0 , 28, 1, 2); book(_h_RKPlus , 30, 1, 1); book(_h_RKMinus , 30, 1, 2); book(_h_RProton , 32, 1, 1); book(_h_RPBar , 32, 1, 2); book(_h_RLambda , 34, 1, 1); book(_h_RLBar , 34, 1, 2); book(_s_Xp_PiPl_Ch , 1, 1, 1); book(_s_Xp_KPl_Ch , 2, 1, 1); book(_s_Xp_Pr_Ch , 3, 1, 1); book(_s_Xp_PiPlCh_PiPlLi, 11, 1, 1); book(_s_Xp_PiPlBo_PiPlLi, 11, 1, 2); book(_s_Xp_KPlCh_KPlLi , 13, 1, 1); book(_s_Xp_KPlBo_KPlLi , 13, 1, 2); book(_s_Xp_KS0Ch_KS0Li , 15, 1, 1); book(_s_Xp_KS0Bo_KS0Li , 15, 1, 2); book(_s_Xp_PrCh_PrLi , 17, 1, 1); book(_s_Xp_PrBo_PrLi , 17, 1, 2); book(_s_Xp_LaCh_LaLi , 19, 1, 1); book(_s_Xp_LaBo_LaLi , 19, 1, 2); book(_s_Xp_K0Ch_K0Li , 21, 1, 1); book(_s_Xp_K0Bo_K0Li , 21, 1, 2); book(_s_Xp_PhiCh_PhiLi , 23, 1, 1); book(_s_Xp_PhiBo_PhiLi , 23, 1, 2); book(_s_PiM_PiP , 27, 1, 1); book(_s_KSBar0_KS0, 29, 1, 1); book(_s_KM_KP , 31, 1, 1); book(_s_Pr_PBar , 33, 1, 1); book(_s_Lam_LBar , 35, 1, 1); book(_SumOfudsWeights, "SumOfudsWeights"); book(_SumOfcWeights, "SumOfcWeights"); book(_SumOfbWeights, "SumOfbWeights"); for ( size_t i=0; i<4; ++i) { book(_multPiPlus[i], "multPiPlus_"+to_str(i)); book(_multKPlus[i], "multKPlus_"+to_str(i)); book(_multK0[i], "multK0_"+to_str(i)); book(_multKStar0[i], "multKStar0_"+to_str(i)); book(_multPhi[i], "multPhi_"+to_str(i)); book(_multProton[i], "multProton_"+to_str(i)); book(_multLambda[i], "multLambda_"+to_str(i)); } book(tmp1, 24, 1, 1, true); book(tmp2, 24, 1, 2, true); book(tmp3, 24, 1, 3, true); book(tmp4, 24, 1, 4, true); book(tmp5, 25, 1, 1, true); book(tmp6, 25, 1, 2, true); book(tmp7, 24, 2, 1, true); book(tmp8, 24, 2, 2, true); book(tmp9, 24, 2, 3, true); book(tmp10, 24, 2, 4, true); book(tmp11, 25, 2, 1, true); book(tmp12, 25, 2, 2, true); book(tmp13, 24, 3, 1, true); book(tmp14, 24, 3, 2, true); book(tmp15, 24, 3, 3, true); book(tmp16, 24, 3, 4, true); book(tmp17, 25, 3, 1, true); book(tmp18, 25, 3, 2, true); book(tmp19, 24, 4, 1, true); book(tmp20, 24, 4, 2, true); book(tmp21, 24, 4, 3, true); book(tmp22, 24, 4, 4, true); book(tmp23, 25, 4, 1, true); book(tmp24, 25, 4, 2, true); book(tmp25, 24, 5, 1, true); book(tmp26, 24, 5, 2, true); book(tmp27, 24, 5, 3, true); book(tmp28, 24, 5, 4, true); book(tmp29, 25, 5, 1, true); book(tmp30, 25, 5, 2, true); book(tmp31, 24, 6, 1, true); book(tmp32, 24, 6, 2, true); book(tmp33, 24, 6, 3, true); book(tmp34, 24, 6, 4, true); book(tmp35, 25, 6, 1, true); book(tmp36, 25, 6, 2, true); book(tmp37, 24, 7, 1, true); book(tmp38, 24, 7, 2, true); book(tmp39, 24, 7, 3, true); book(tmp40, 24, 7, 4, true); book(tmp41, 25, 7, 1, true); book(tmp42, 25, 7, 2, true); } /// Finalize void finalize() { // Get the ratio plots sorted out first divide(_h_XpPiPlusN, _temp_XpChargedN1, _s_Xp_PiPl_Ch); divide(_h_XpKPlusN, _temp_XpChargedN2, _s_Xp_KPl_Ch); divide(_h_XpProtonN, _temp_XpChargedN3, _s_Xp_Pr_Ch); divide(_h_XpPiPlusCharm, _h_XpPiPlusLight, _s_Xp_PiPlCh_PiPlLi); _s_Xp_PiPlCh_PiPlLi->scale(1., dbl(*_SumOfudsWeights / *_SumOfcWeights)); divide(_h_XpPiPlusBottom, _h_XpPiPlusLight, _s_Xp_PiPlBo_PiPlLi); _s_Xp_PiPlBo_PiPlLi->scale(1., dbl(*_SumOfudsWeights / *_SumOfbWeights)); divide(_temp_XpKPlusCharm , _temp_XpKPlusLight, _s_Xp_KPlCh_KPlLi); _s_Xp_KPlCh_KPlLi->scale(1., dbl(*_SumOfudsWeights / *_SumOfcWeights)); divide(_h_XpKPlusBottom, _h_XpKPlusLight, _s_Xp_KPlBo_KPlLi); _s_Xp_KPlBo_KPlLi->scale(1., dbl(*_SumOfudsWeights / *_SumOfbWeights)); divide(_temp_XpKStar0Charm, _temp_XpKStar0Light, _s_Xp_KS0Ch_KS0Li); _s_Xp_KS0Ch_KS0Li->scale(1., dbl(*_SumOfudsWeights / *_SumOfcWeights)); divide(_h_XpKStar0Bottom, _h_XpKStar0Light, _s_Xp_KS0Bo_KS0Li); _s_Xp_KS0Bo_KS0Li->scale(1., dbl(*_SumOfudsWeights / *_SumOfbWeights)); divide(_temp_XpProtonCharm, _temp_XpProtonLight, _s_Xp_PrCh_PrLi); _s_Xp_PrCh_PrLi->scale(1., dbl(*_SumOfudsWeights / *_SumOfcWeights)); divide(_h_XpProtonBottom, _h_XpProtonLight, _s_Xp_PrBo_PrLi); _s_Xp_PrBo_PrLi->scale(1., dbl(*_SumOfudsWeights / *_SumOfbWeights)); divide(_h_XpLambdaCharm, _h_XpLambdaLight, _s_Xp_LaCh_LaLi); _s_Xp_LaCh_LaLi->scale(1., dbl(*_SumOfudsWeights / *_SumOfcWeights)); divide(_h_XpLambdaBottom, _h_XpLambdaLight, _s_Xp_LaBo_LaLi); _s_Xp_LaBo_LaLi->scale(1., dbl(*_SumOfudsWeights / *_SumOfbWeights)); divide(_h_XpK0Charm, _h_XpK0Light, _s_Xp_K0Ch_K0Li); _s_Xp_K0Ch_K0Li->scale(1., dbl(*_SumOfudsWeights / *_SumOfcWeights)); divide(_h_XpK0Bottom, _h_XpK0Light, _s_Xp_K0Bo_K0Li); _s_Xp_K0Bo_K0Li->scale(1., dbl(*_SumOfudsWeights / *_SumOfbWeights)); divide(_h_XpPhiCharm, _h_XpPhiLight, _s_Xp_PhiCh_PhiLi); _s_Xp_PhiCh_PhiLi->scale(1., dbl(*_SumOfudsWeights / *_SumOfcWeights)); divide(_h_XpPhiBottom, _h_XpPhiLight, _s_Xp_PhiBo_PhiLi); _s_Xp_PhiBo_PhiLi->scale(1., dbl(*_SumOfudsWeights / *_SumOfbWeights)); // Then the leading particles divide(*_h_RPiMinus - *_h_RPiPlus, *_h_RPiMinus + *_h_RPiPlus, _s_PiM_PiP); divide(*_h_RKSBar0 - *_h_RKS0, *_h_RKSBar0 + *_h_RKS0, _s_KSBar0_KS0); divide(*_h_RKMinus - *_h_RKPlus, *_h_RKMinus + *_h_RKPlus, _s_KM_KP); divide(*_h_RProton - *_h_RPBar, *_h_RProton + *_h_RPBar, _s_Pr_PBar); divide(*_h_RLambda - *_h_RLBar, *_h_RLambda + *_h_RLBar, _s_Lam_LBar); // Then the rest scale(_h_XpPiPlusN, 1/sumOfWeights()); scale(_h_XpKPlusN, 1/sumOfWeights()); scale(_h_XpProtonN, 1/sumOfWeights()); scale(_h_XpChargedN, 1/sumOfWeights()); scale(_h_XpK0N, 1/sumOfWeights()); scale(_h_XpLambdaN, 1/sumOfWeights()); scale(_h_XpKStar0N, 1/sumOfWeights()); scale(_h_XpPhiN, 1/sumOfWeights()); scale(_h_XpPiPlusLight, 1 / *_SumOfudsWeights); scale(_h_XpPiPlusCharm, 1 / *_SumOfcWeights); scale(_h_XpPiPlusBottom, 1 / *_SumOfbWeights); scale(_h_XpKPlusLight, 1 / *_SumOfudsWeights); scale(_h_XpKPlusCharm, 1 / *_SumOfcWeights); scale(_h_XpKPlusBottom, 1 / *_SumOfbWeights); scale(_h_XpKStar0Light, 1 / *_SumOfudsWeights); scale(_h_XpKStar0Charm, 1 / *_SumOfcWeights); scale(_h_XpKStar0Bottom, 1 / *_SumOfbWeights); scale(_h_XpProtonLight, 1 / *_SumOfudsWeights); scale(_h_XpProtonCharm, 1 / *_SumOfcWeights); scale(_h_XpProtonBottom, 1 / *_SumOfbWeights); scale(_h_XpLambdaLight, 1 / *_SumOfudsWeights); scale(_h_XpLambdaCharm, 1 / *_SumOfcWeights); scale(_h_XpLambdaBottom, 1 / *_SumOfbWeights); scale(_h_XpK0Light, 1 / *_SumOfudsWeights); scale(_h_XpK0Charm, 1 / *_SumOfcWeights); scale(_h_XpK0Bottom, 1 / *_SumOfbWeights); scale(_h_XpPhiLight, 1 / *_SumOfudsWeights); scale(_h_XpPhiCharm , 1 / *_SumOfcWeights); scale(_h_XpPhiBottom, 1 / *_SumOfbWeights); scale(_h_RPiPlus, 1 / *_SumOfudsWeights); scale(_h_RPiMinus, 1 / *_SumOfudsWeights); scale(_h_RKS0, 1 / *_SumOfudsWeights); scale(_h_RKSBar0, 1 / *_SumOfudsWeights); scale(_h_RKPlus, 1 / *_SumOfudsWeights); scale(_h_RKMinus, 1 / *_SumOfudsWeights); scale(_h_RProton, 1 / *_SumOfudsWeights); scale(_h_RPBar, 1 / *_SumOfudsWeights); scale(_h_RLambda, 1 / *_SumOfudsWeights); scale(_h_RLBar, 1 / *_SumOfudsWeights); // Multiplicities double avgNumPartsAll, avgNumPartsLight,avgNumPartsCharm, avgNumPartsBottom; // pi+/- // all avgNumPartsAll = dbl(*_multPiPlus[0])/sumOfWeights(); tmp1->point(0).setY(avgNumPartsAll); // light avgNumPartsLight = dbl(*_multPiPlus[1] / *_SumOfudsWeights); tmp2->point(0).setY(avgNumPartsLight); // charm avgNumPartsCharm = dbl(*_multPiPlus[2] / *_SumOfcWeights); tmp3->point(0).setY(avgNumPartsCharm); // bottom avgNumPartsBottom = dbl(*_multPiPlus[3] / *_SumOfbWeights); tmp4->point(0).setY(avgNumPartsBottom); // charm-light tmp5->point(0).setY(avgNumPartsCharm - avgNumPartsLight); // bottom-light tmp6->point(0).setY(avgNumPartsBottom - avgNumPartsLight); // K+/- // all avgNumPartsAll = dbl(*_multKPlus[0])/sumOfWeights(); tmp7->point(0).setY(avgNumPartsAll); // light avgNumPartsLight = dbl(*_multKPlus[1] / *_SumOfudsWeights); tmp8->point(0).setY(avgNumPartsLight); // charm avgNumPartsCharm = dbl(*_multKPlus[2] / *_SumOfcWeights); tmp9->point(0).setY(avgNumPartsCharm); // bottom avgNumPartsBottom = dbl(*_multKPlus[3] / *_SumOfbWeights); tmp10->point(0).setY(avgNumPartsBottom); // charm-light tmp11->point(0).setY(avgNumPartsCharm - avgNumPartsLight); // bottom-light tmp12->point(0).setY(avgNumPartsBottom - avgNumPartsLight); // K0 // all avgNumPartsAll = dbl(*_multK0[0])/sumOfWeights(); tmp13->point(0).setY(avgNumPartsAll); // light avgNumPartsLight = dbl(*_multK0[1] / *_SumOfudsWeights); tmp14->point(0).setY(avgNumPartsLight); // charm avgNumPartsCharm = dbl(*_multK0[2] / *_SumOfcWeights); tmp15->point(0).setY(avgNumPartsCharm); // bottom avgNumPartsBottom = dbl(*_multK0[3] / *_SumOfbWeights); tmp16->point(0).setY(avgNumPartsBottom); // charm-light tmp17->point(0).setY(avgNumPartsCharm - avgNumPartsLight); // bottom-light tmp18->point(0).setY(avgNumPartsBottom - avgNumPartsLight); // K*0 // all avgNumPartsAll = dbl(*_multKStar0[0])/sumOfWeights(); tmp19->point(0).setY(avgNumPartsAll); // light avgNumPartsLight = dbl(*_multKStar0[1] / *_SumOfudsWeights); tmp20->point(0).setY(avgNumPartsLight); // charm avgNumPartsCharm = dbl(*_multKStar0[2] / *_SumOfcWeights); tmp21->point(0).setY(avgNumPartsCharm); // bottom avgNumPartsBottom = dbl(*_multKStar0[3] / *_SumOfbWeights); tmp22->point(0).setY(avgNumPartsBottom); // charm-light tmp23->point(0).setY(avgNumPartsCharm - avgNumPartsLight); // bottom-light tmp24->point(0).setY(avgNumPartsBottom - avgNumPartsLight); // phi // all avgNumPartsAll = dbl(*_multPhi[0])/sumOfWeights(); tmp25->point(0).setY(avgNumPartsAll); // light avgNumPartsLight = dbl(*_multPhi[1] / *_SumOfudsWeights); tmp26->point(0).setY(avgNumPartsLight); // charm avgNumPartsCharm = dbl(*_multPhi[2] / *_SumOfcWeights); tmp27->point(0).setY(avgNumPartsCharm); // bottom avgNumPartsBottom = dbl(*_multPhi[3] / *_SumOfbWeights); tmp28->point(0).setY(avgNumPartsBottom); // charm-light tmp29->point(0).setY(avgNumPartsCharm - avgNumPartsLight); // bottom-light tmp30->point(0).setY(avgNumPartsBottom - avgNumPartsLight); // p // all avgNumPartsAll = dbl(*_multProton[0])/sumOfWeights(); tmp31->point(0).setY(avgNumPartsAll); // light avgNumPartsLight = dbl(*_multProton[1] / *_SumOfudsWeights); tmp32->point(0).setY(avgNumPartsLight); // charm avgNumPartsCharm = dbl(*_multProton[2] / *_SumOfcWeights); tmp33->point(0).setY(avgNumPartsCharm); // bottom avgNumPartsBottom = dbl(*_multProton[3] / *_SumOfbWeights); tmp34->point(0).setY(avgNumPartsBottom); // charm-light tmp35->point(0).setY(avgNumPartsCharm - avgNumPartsLight); // bottom-light tmp36->point(0).setY(avgNumPartsBottom - avgNumPartsLight); // Lambda // all avgNumPartsAll = dbl(*_multLambda[0])/sumOfWeights(); tmp37->point(0).setY(avgNumPartsAll); // light avgNumPartsLight = dbl(*_multLambda[1] / *_SumOfudsWeights); tmp38->point(0).setY(avgNumPartsLight); // charm avgNumPartsCharm = dbl(*_multLambda[2] / *_SumOfcWeights); tmp39->point(0).setY(avgNumPartsCharm); // bottom avgNumPartsBottom = dbl(*_multLambda[3] / *_SumOfbWeights); tmp40->point(0).setY(avgNumPartsBottom); // charm-light tmp41->point(0).setY(avgNumPartsCharm - avgNumPartsLight); // bottom-light tmp42->point(0).setY(avgNumPartsBottom - avgNumPartsLight); } //@} private: /// Store the weighted sums of numbers of charged / charged+neutral /// particles. Used to calculate average number of particles for the /// inclusive single particle distributions' normalisations. CounterPtr _SumOfudsWeights, _SumOfcWeights, _SumOfbWeights; vector _multPiPlus, _multKPlus, _multK0, _multKStar0, _multPhi, _multProton, _multLambda; Histo1DPtr _h_XpPiPlusSig, _h_XpPiPlusN; Histo1DPtr _h_XpKPlusSig, _h_XpKPlusN; Histo1DPtr _h_XpProtonSig, _h_XpProtonN; Histo1DPtr _h_XpChargedN; Histo1DPtr _h_XpK0N, _h_XpLambdaN; Histo1DPtr _h_XpKStar0N, _h_XpPhiN; Histo1DPtr _h_XpPiPlusLight, _h_XpPiPlusCharm, _h_XpPiPlusBottom; Histo1DPtr _h_XpKPlusLight, _h_XpKPlusCharm, _h_XpKPlusBottom; Histo1DPtr _h_XpKStar0Light, _h_XpKStar0Charm, _h_XpKStar0Bottom; Histo1DPtr _h_XpProtonLight, _h_XpProtonCharm, _h_XpProtonBottom; Histo1DPtr _h_XpLambdaLight, _h_XpLambdaCharm, _h_XpLambdaBottom; Histo1DPtr _h_XpK0Light, _h_XpK0Charm, _h_XpK0Bottom; Histo1DPtr _h_XpPhiLight, _h_XpPhiCharm, _h_XpPhiBottom; Histo1DPtr _temp_XpChargedN1, _temp_XpChargedN2, _temp_XpChargedN3; Histo1DPtr _temp_XpKPlusCharm , _temp_XpKPlusLight; Histo1DPtr _temp_XpKStar0Charm, _temp_XpKStar0Light; Histo1DPtr _temp_XpProtonCharm, _temp_XpProtonLight; Histo1DPtr _h_RPiPlus, _h_RPiMinus; Histo1DPtr _h_RKS0, _h_RKSBar0; Histo1DPtr _h_RKPlus, _h_RKMinus; Histo1DPtr _h_RProton, _h_RPBar; Histo1DPtr _h_RLambda, _h_RLBar; Scatter2DPtr _s_Xp_PiPl_Ch, _s_Xp_KPl_Ch, _s_Xp_Pr_Ch; Scatter2DPtr _s_Xp_PiPlCh_PiPlLi, _s_Xp_PiPlBo_PiPlLi; Scatter2DPtr _s_Xp_KPlCh_KPlLi, _s_Xp_KPlBo_KPlLi; Scatter2DPtr _s_Xp_KS0Ch_KS0Li, _s_Xp_KS0Bo_KS0Li; Scatter2DPtr _s_Xp_PrCh_PrLi, _s_Xp_PrBo_PrLi; Scatter2DPtr _s_Xp_LaCh_LaLi, _s_Xp_LaBo_LaLi; Scatter2DPtr _s_Xp_K0Ch_K0Li, _s_Xp_K0Bo_K0Li; Scatter2DPtr _s_Xp_PhiCh_PhiLi, _s_Xp_PhiBo_PhiLi; Scatter2DPtr _s_PiM_PiP, _s_KSBar0_KS0, _s_KM_KP, _s_Pr_PBar, _s_Lam_LBar; //@} Scatter2DPtr tmp1; Scatter2DPtr tmp2; Scatter2DPtr tmp3; Scatter2DPtr tmp4; Scatter2DPtr tmp5; Scatter2DPtr tmp6; Scatter2DPtr tmp7; Scatter2DPtr tmp8; Scatter2DPtr tmp9; Scatter2DPtr tmp10; Scatter2DPtr tmp11; Scatter2DPtr tmp12; Scatter2DPtr tmp13; Scatter2DPtr tmp14; Scatter2DPtr tmp15; Scatter2DPtr tmp16; Scatter2DPtr tmp17; Scatter2DPtr tmp18; Scatter2DPtr tmp19; Scatter2DPtr tmp20; Scatter2DPtr tmp21; Scatter2DPtr tmp22; Scatter2DPtr tmp23; Scatter2DPtr tmp24; Scatter2DPtr tmp25; Scatter2DPtr tmp26; Scatter2DPtr tmp27; Scatter2DPtr tmp28; Scatter2DPtr tmp29; Scatter2DPtr tmp30; Scatter2DPtr tmp31; Scatter2DPtr tmp32; Scatter2DPtr tmp33; Scatter2DPtr tmp34; Scatter2DPtr tmp35; Scatter2DPtr tmp36; Scatter2DPtr tmp37; Scatter2DPtr tmp38; Scatter2DPtr tmp39; Scatter2DPtr tmp40; Scatter2DPtr tmp41; Scatter2DPtr tmp42; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(SLD_1999_S3743934); } diff --git a/analyses/pluginLEP/SLD_2004_S5693039.cc b/analyses/pluginLEP/SLD_2004_S5693039.cc --- a/analyses/pluginLEP/SLD_2004_S5693039.cc +++ b/analyses/pluginLEP/SLD_2004_S5693039.cc @@ -1,377 +1,379 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/Beam.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" +#include "Rivet/Projections/Thrust.hh" + +#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT #include "Rivet/Projections/InitialQuarks.hh" -#include "Rivet/Projections/Thrust.hh" namespace Rivet { /// @brief SLD flavour-dependent fragmentation paper /// @author Peter Richardson class SLD_2004_S5693039 : public Analysis { public: /// Constructor SLD_2004_S5693039() : Analysis("SLD_2004_S5693039") {} /// @name Analysis methods //@{ void analyze(const Event& e) { // First, veto on leptonic events by requiring at least 2 charged FS particles const FinalState& fs = apply(e, "FS"); const size_t numParticles = fs.particles().size(); // Even if we only generate hadronic events, we still need a cut on numCharged >= 2. if (numParticles < 2) { MSG_DEBUG("Failed ncharged cut"); vetoEvent; } MSG_DEBUG("Passed ncharged cut"); // Get beams and average beam momentum const ParticlePair& beams = apply(e, "Beams").beams(); const double meanBeamMom = ( beams.first.p3().mod() + beams.second.p3().mod() ) / 2.0; MSG_DEBUG("Avg beam momentum = " << meanBeamMom); int flavour = 0; const InitialQuarks& iqf = apply(e, "IQF"); // If we only have two quarks (qqbar), just take the flavour. // If we have more than two quarks, look for the highest energetic q-qbar pair. Particles quarks; if (iqf.particles().size() == 2) { flavour = iqf.particles().front().abspid(); quarks = iqf.particles(); } else { map quarkmap; for (const Particle& p : iqf.particles()) { if (quarkmap.find(p.pid())==quarkmap.end()) quarkmap[p.pid()] = p; else if (quarkmap[p.pid()].E() < p.E()) quarkmap[p.pid()] = p; } double maxenergy = 0.; for (int i = 1; i <= 5; ++i) { double energy(0.); if(quarkmap.find( i)!=quarkmap.end()) energy += quarkmap[ i].E(); if(quarkmap.find(-i)!=quarkmap.end()) energy += quarkmap[-i].E(); if (energy > maxenergy) flavour = i; } if(quarkmap.find( flavour)!=quarkmap.end()) quarks.push_back(quarkmap[ flavour]); if(quarkmap.find(-flavour)!=quarkmap.end()) quarks.push_back(quarkmap[-flavour]); } // total multiplicities switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _weightLight ->fill(); _weightedTotalChargedPartNumLight ->fill(numParticles); break; case PID::CQUARK: _weightCharm ->fill(); _weightedTotalChargedPartNumCharm ->fill(numParticles); break; case PID::BQUARK: _weightBottom ->fill(); _weightedTotalChargedPartNumBottom ->fill(numParticles); break; } // thrust axis for projections Vector3 axis = apply(e, "Thrust").thrustAxis(); double dot(0.); if(!quarks.empty()) { dot = quarks[0].p3().dot(axis); if(quarks[0].pid()<0) dot *= -1.; } // spectra and individual multiplicities for (const Particle& p : fs.particles()) { double pcm = p.p3().mod(); const double xp = pcm/meanBeamMom; // if in quark or antiquark hemisphere bool quark = p.p3().dot(axis)*dot>0.; _h_PCharged ->fill(pcm ); // all charged switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _h_XpChargedL->fill(xp); break; case PID::CQUARK: _h_XpChargedC->fill(xp); break; case PID::BQUARK: _h_XpChargedB->fill(xp); break; } int id = p.abspid(); // charged pions if (id == PID::PIPLUS) { _h_XpPiPlus->fill(xp); _h_XpPiPlusTotal->fill(xp); switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _h_XpPiPlusL->fill(xp); _h_NPiPlusL->fill(sqrtS()); if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 )) _h_RPiPlus->fill(xp); else _h_RPiMinus->fill(xp); break; case PID::CQUARK: _h_XpPiPlusC->fill(xp); _h_NPiPlusC->fill(sqrtS()); break; case PID::BQUARK: _h_XpPiPlusB->fill(xp); _h_NPiPlusB->fill(sqrtS()); break; } } else if (id == PID::KPLUS) { _h_XpKPlus->fill(xp); _h_XpKPlusTotal->fill(xp); switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _h_XpKPlusL->fill(xp); _h_NKPlusL->fill(sqrtS()); if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 )) _h_RKPlus->fill(xp); else _h_RKMinus->fill(xp); break; case PID::CQUARK: _h_XpKPlusC->fill(xp); _h_NKPlusC->fill(sqrtS()); break; case PID::BQUARK: _h_XpKPlusB->fill(xp); _h_NKPlusB->fill(sqrtS()); break; } } else if (id == PID::PROTON) { _h_XpProton->fill(xp); _h_XpProtonTotal->fill(xp); switch (flavour) { case PID::DQUARK: case PID::UQUARK: case PID::SQUARK: _h_XpProtonL->fill(xp); _h_NProtonL->fill(sqrtS()); if( ( quark && p.pid()>0 ) || ( !quark && p.pid()<0 )) _h_RProton->fill(xp); else _h_RPBar ->fill(xp); break; case PID::CQUARK: _h_XpProtonC->fill(xp); _h_NProtonC->fill(sqrtS()); break; case PID::BQUARK: _h_XpProtonB->fill(xp); _h_NProtonB->fill(sqrtS()); break; } } } } void init() { // Projections declare(Beam(), "Beams"); declare(ChargedFinalState(), "FS"); declare(InitialQuarks(), "IQF"); declare(Thrust(FinalState()), "Thrust"); // Book histograms book(_h_PCharged , 1, 1, 1); book(_h_XpPiPlus , 2, 1, 2); book(_h_XpKPlus , 3, 1, 2); book(_h_XpProton , 4, 1, 2); book(_h_XpPiPlusTotal , 2, 2, 2); book(_h_XpKPlusTotal , 3, 2, 2); book(_h_XpProtonTotal , 4, 2, 2); book(_h_XpPiPlusL , 5, 1, 1); book(_h_XpPiPlusC , 5, 1, 2); book(_h_XpPiPlusB , 5, 1, 3); book(_h_XpKPlusL , 6, 1, 1); book(_h_XpKPlusC , 6, 1, 2); book(_h_XpKPlusB , 6, 1, 3); book(_h_XpProtonL , 7, 1, 1); book(_h_XpProtonC , 7, 1, 2); book(_h_XpProtonB , 7, 1, 3); book(_h_XpChargedL , 8, 1, 1); book(_h_XpChargedC , 8, 1, 2); book(_h_XpChargedB , 8, 1, 3); book(_h_NPiPlusL , 5, 2, 1); book(_h_NPiPlusC , 5, 2, 2); book(_h_NPiPlusB , 5, 2, 3); book(_h_NKPlusL , 6, 2, 1); book(_h_NKPlusC , 6, 2, 2); book(_h_NKPlusB , 6, 2, 3); book(_h_NProtonL , 7, 2, 1); book(_h_NProtonC , 7, 2, 2); book(_h_NProtonB , 7, 2, 3); book(_h_RPiPlus , 9, 1, 1); book(_h_RPiMinus , 9, 1, 2); book(_h_RKPlus ,10, 1, 1); book(_h_RKMinus ,10, 1, 2); book(_h_RProton ,11, 1, 1); book(_h_RPBar ,11, 1, 2); // Ratios: used as target of divide() later book(_s_PiM_PiP, 9, 1, 3); book(_s_KM_KP , 10, 1, 3); book(_s_Pr_PBar, 11, 1, 3); book(_weightedTotalChargedPartNumLight, "weightedTotalChargedPartNumLight"); book(_weightedTotalChargedPartNumCharm, "weightedTotalChargedPartNumCharm"); book(_weightedTotalChargedPartNumBottom, "weightedTotalChargedPartNumBottom"); book(_weightLight, "weightLight"); book(_weightCharm, "weightCharm"); book(_weightBottom, "weightBottom"); book(tmp1, 8, 2, 1, true); book(tmp2, 8, 2, 2, true); book(tmp3, 8, 2, 3, true); book(tmp4, 8, 3, 2, true); book(tmp5, 8, 3, 3, true); } /// Finalize void finalize() { // Multiplicities /// @todo Include errors const double avgNumPartsLight = _weightedTotalChargedPartNumLight->val() / _weightLight->val(); const double avgNumPartsCharm = _weightedTotalChargedPartNumCharm->val() / _weightCharm->val(); const double avgNumPartsBottom = _weightedTotalChargedPartNumBottom->val() / _weightBottom->val(); tmp1->point(0).setY(avgNumPartsLight); tmp2->point(0).setY(avgNumPartsCharm); tmp3->point(0).setY(avgNumPartsBottom); tmp4->point(0).setY(avgNumPartsCharm - avgNumPartsLight); tmp5->point(0).setY(avgNumPartsBottom - avgNumPartsLight); // Do divisions divide(*_h_RPiMinus - *_h_RPiPlus, *_h_RPiMinus + *_h_RPiPlus, _s_PiM_PiP); divide(*_h_RKMinus - *_h_RKPlus, *_h_RKMinus + *_h_RKPlus, _s_KM_KP); divide(*_h_RProton - *_h_RPBar, *_h_RProton + *_h_RPBar, _s_Pr_PBar); // Scale histograms scale(_h_PCharged, 1./sumOfWeights()); scale(_h_XpPiPlus, 1./sumOfWeights()); scale(_h_XpKPlus, 1./sumOfWeights()); scale(_h_XpProton, 1./sumOfWeights()); scale(_h_XpPiPlusTotal, 1./sumOfWeights()); scale(_h_XpKPlusTotal, 1./sumOfWeights()); scale(_h_XpProtonTotal, 1./sumOfWeights()); scale(_h_XpPiPlusL, 1. / *_weightLight); scale(_h_XpPiPlusC, 1. / *_weightCharm); scale(_h_XpPiPlusB, 1. / *_weightBottom); scale(_h_XpKPlusL, 1. / *_weightLight); scale(_h_XpKPlusC, 1. / *_weightCharm); scale(_h_XpKPlusB, 1. / *_weightBottom); scale(_h_XpProtonL, 1. / *_weightLight); scale(_h_XpProtonC, 1. / *_weightCharm); scale(_h_XpProtonB, 1. / *_weightBottom); scale(_h_XpChargedL, 1. / *_weightLight); scale(_h_XpChargedC, 1. / *_weightCharm); scale(_h_XpChargedB, 1. / *_weightBottom); scale(_h_NPiPlusL, 1. / *_weightLight); scale(_h_NPiPlusC, 1. / *_weightCharm); scale(_h_NPiPlusB, 1. / *_weightBottom); scale(_h_NKPlusL, 1. / *_weightLight); scale(_h_NKPlusC, 1. / *_weightCharm); scale(_h_NKPlusB, 1. / *_weightBottom); scale(_h_NProtonL, 1. / *_weightLight); scale(_h_NProtonC, 1. / *_weightCharm); scale(_h_NProtonB, 1. / *_weightBottom); // Paper suggests this should be 0.5/weight but it has to be 1.0 to get normalisations right... scale(_h_RPiPlus, 1. / *_weightLight); scale(_h_RPiMinus, 1. / *_weightLight); scale(_h_RKPlus, 1. / *_weightLight); scale(_h_RKMinus, 1. / *_weightLight); scale(_h_RProton, 1. / *_weightLight); scale(_h_RPBar, 1. / *_weightLight); // convert ratio to % _s_PiM_PiP->scale(1.,100.); _s_KM_KP ->scale(1.,100.); _s_Pr_PBar->scale(1.,100.); } //@} private: Scatter2DPtr tmp1; Scatter2DPtr tmp2; Scatter2DPtr tmp3; Scatter2DPtr tmp4; Scatter2DPtr tmp5; /// @name Multiplicities //@{ CounterPtr _weightedTotalChargedPartNumLight; CounterPtr _weightedTotalChargedPartNumCharm; CounterPtr _weightedTotalChargedPartNumBottom; //@} /// @name Weights //@{ CounterPtr _weightLight, _weightCharm, _weightBottom; //@} // Histograms //@{ Histo1DPtr _h_PCharged; Histo1DPtr _h_XpPiPlus, _h_XpKPlus, _h_XpProton; Histo1DPtr _h_XpPiPlusTotal, _h_XpKPlusTotal, _h_XpProtonTotal; Histo1DPtr _h_XpPiPlusL, _h_XpPiPlusC, _h_XpPiPlusB; Histo1DPtr _h_XpKPlusL, _h_XpKPlusC, _h_XpKPlusB; Histo1DPtr _h_XpProtonL, _h_XpProtonC, _h_XpProtonB; Histo1DPtr _h_XpChargedL, _h_XpChargedC, _h_XpChargedB; Histo1DPtr _h_NPiPlusL, _h_NPiPlusC, _h_NPiPlusB; Histo1DPtr _h_NKPlusL, _h_NKPlusC, _h_NKPlusB; Histo1DPtr _h_NProtonL, _h_NProtonC, _h_NProtonB; Histo1DPtr _h_RPiPlus, _h_RPiMinus, _h_RKPlus; Histo1DPtr _h_RKMinus, _h_RProton, _h_RPBar; Scatter2DPtr _s_PiM_PiP, _s_KM_KP, _s_Pr_PBar; //@} }; // Hook for the plugin system DECLARE_RIVET_PLUGIN(SLD_2004_S5693039); } diff --git a/analyses/pluginMisc/TPC_1987_I235694.cc b/analyses/pluginMisc/TPC_1987_I235694.cc --- a/analyses/pluginMisc/TPC_1987_I235694.cc +++ b/analyses/pluginMisc/TPC_1987_I235694.cc @@ -1,102 +1,104 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/Beam.hh" #include "Rivet/Projections/ChargedFinalState.hh" + +#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT #include "Rivet/Projections/InitialQuarks.hh" namespace Rivet { /// @brief TPC flavour separated N charged at 29 GeV class TPC_1987_I235694 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(TPC_1987_I235694); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Histograms book(_h_all , 5, 1, 4); book(_h_light , 4, 1, 4); book(_h_charm , 3, 1, 4); book(_h_bottom, 2, 1, 4); // Projections declare(Beam(), "Beams"); declare(ChargedFinalState(), "CFS"); declare(InitialQuarks(), "IQF"); } /// Perform the per-event analysis void analyze(const Event& event) { // Even if we only generate hadronic events, we still need a cut on numCharged >= 2. const FinalState& cfs = apply(event, "CFS"); if (cfs.size() < 2) vetoEvent; int flavour = 0; const InitialQuarks& iqf = apply(event, "IQF"); // If we only have two quarks (qqbar), just take the flavour. // If we have more than two quarks, look for the highest energetic q-qbar pair. if (iqf.particles().size() == 2) { flavour = iqf.particles().front().abspid(); } else { map quarkmap; for (const Particle& p : iqf.particles()) { if (quarkmap[p.pid()] < p.E()) { quarkmap[p.pid()] = p.E(); } } double maxenergy = 0.; for (int i = 1; i <= 5; ++i) { if (quarkmap[i]+quarkmap[-i] > maxenergy) { flavour = i; } } } const size_t numParticles = cfs.particles().size(); switch (flavour) { case 1: case 2: case 3: _h_light->fill(sqrtS()/GeV,numParticles); break; case 4: _h_charm->fill(sqrtS()/GeV,numParticles); break; case 5: _h_bottom->fill(sqrtS()/GeV,numParticles); break; } _h_all->fill(sqrtS()/GeV,numParticles); } /// Normalise histograms etc., after the run void finalize() { } //@} private: /// @name Multiplicities //@{ Profile1DPtr _h_all; Profile1DPtr _h_light; Profile1DPtr _h_charm; Profile1DPtr _h_bottom; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(TPC_1987_I235694); } diff --git a/include/Rivet/Projections/InitialQuarks.hh b/include/Rivet/Projections/InitialQuarks.hh --- a/include/Rivet/Projections/InitialQuarks.hh +++ b/include/Rivet/Projections/InitialQuarks.hh @@ -1,61 +1,63 @@ // -*- C++ -*- #ifndef RIVET_InitialQuarks_HH #define RIVET_InitialQuarks_HH +#ifndef I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT #warning "This is a dangerous projection for a few specific old analyses. Not for general use!" +#endif #include "Rivet/Projection.hh" #include "Rivet/Particle.hh" #include "Rivet/Event.hh" namespace Rivet { /// @brief Project out quarks from the hard process in \f$ e^+ e^- \to Z^0 \f$ events /// /// @deprecated We're not sure exactly when we'lll get rid of this, but it's going to happen... /// /// @warning This is a very dangerous and specific projection! class InitialQuarks : public Projection { public: /// @name Standard constructors and destructors. //@{ /// The default constructor. May specify the minimum and maximum /// pseudorapidity \f$ \eta \f$ and the min \f$ p_T \f$ (in GeV). InitialQuarks() { setName("InitialQuarks"); } /// Clone on the heap. DEFAULT_RIVET_PROJ_CLONE(InitialQuarks); //@} /// Access the projected final-state particles. virtual const Particles& particles() const { return _theParticles; } /// Is this final state empty? virtual bool empty() const { return _theParticles.empty(); } protected: /// Apply the projection to the event. virtual void project(const Event& e); /// Compare projections. virtual CmpState compare(const Projection& p) const; protected: /// The final-state particles. Particles _theParticles; }; } #endif diff --git a/include/Rivet/Projections/JetAlg.hh b/include/Rivet/Projections/JetAlg.hh --- a/include/Rivet/Projections/JetAlg.hh +++ b/include/Rivet/Projections/JetAlg.hh @@ -1,219 +1,214 @@ // -*- C++ -*- #ifndef RIVET_JetAlg_HH #define RIVET_JetAlg_HH #include "Rivet/Projection.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/VisibleFinalState.hh" #include "Rivet/Particle.hh" #include "Rivet/Jet.hh" namespace Rivet { /// Abstract base class for projections which can return a set of {@link Jet}s. class JetAlg : public Projection { public: /// Enum for the treatment of muons: whether to include all, some, or none in jet-finding enum class Muons { NONE, DECAY, ALL }; /// Enum for the treatment of invisible particles: whether to include all, some, or none in jet-finding enum class Invisibles { NONE, DECAY, ALL }; /// Constructor JetAlg(const FinalState& fs, Muons usemuons = Muons::ALL, Invisibles useinvis = Invisibles::NONE); /// Default constructor JetAlg() = default; /// Clone on the heap. virtual unique_ptr clone() const = 0; /// Destructor virtual ~JetAlg() = default; /// @name Control the treatment of muons and invisible particles /// /// Since MC-based jet calibration (and/or particle flow) can add back in /// particles that weren't seen in calorimeters/trackers. //@{ /// @brief Include (some) muons in jet construction. /// /// The default behaviour is that jets are only constructed from visible /// particles. Some jet studies, including those from ATLAS, use a definition /// in which neutrinos from hadron decays are included via MC-based calibrations. /// Setting this flag to true avoids the automatic restriction to a VisibleFinalState. void useMuons(Muons usemuons = Muons::ALL) { _useMuons = usemuons; } /// @brief Include (some) invisible particles in jet construction. /// /// The default behaviour is that jets are only constructed from visible /// particles. Some jet studies, including those from ATLAS, use a definition /// in which neutrinos from hadron decays are included via MC-based calibrations. /// Setting this flag to true avoids the automatic restriction to a VisibleFinalState. void useInvisibles(Invisibles useinvis = Invisibles::DECAY) { _useInvisibles = useinvis; } /// @brief obsolete chooser DEPRECATED("make an explicit choice from Invisibles::{NONE,DECAY,ALL}. This boolean call does not allow for ALL") void useInvisibles(bool useinvis) { _useInvisibles = useinvis ? Invisibles::DECAY : Invisibles::NONE; } //@} /// @name Access to jet objects //@{ /// Get jets in no guaranteed order, with an optional Cut /// @note Returns a copy rather than a reference, due to cuts virtual Jets jets(const Cut& c=Cuts::open()) const { return filter_select(_jets(), c); } /// Get jets in no guaranteed order, with a selection functor /// @note Returns a copy rather than a reference, due to cuts virtual Jets jets(const JetSelector& selector) const { return filter_select(_jets(), selector); } /// Get the jets with a Cut applied, and ordered by supplied sorting functor /// @note Returns a copy rather than a reference, due to cuts and sorting Jets jets(const Cut& c, const JetSorter& sorter) const { /// @todo Will the vector be efficiently std::move'd by value through this function chain? return sortBy(jets(c), sorter); } /// Get the jets, ordered by supplied sorting functor, with an optional Cut /// @note Returns a copy rather than a reference, due to cuts and sorting Jets jets(const JetSorter& sorter, const Cut& c=Cuts::open()) const { /// @todo Will the vector be efficiently std::move'd by value through this function chain? return jets(c, sorter); } /// Get the jets, ordered by supplied sorting function object, with optional cuts on \f$ p_\perp \f$ and rapidity. /// @note Returns a copy rather than a reference, due to cuts and sorting Jets jets(const JetSelector& selector, const JetSorter& sorter) const { /// @todo Will the vector be efficiently std::move'd by value through this function chain? return sortBy(jets(selector), sorter); } /// Get the jets, ordered by supplied sorting functor and with a selection functor applied /// @note Returns a copy rather than a reference, due to cuts and sorting Jets jets(const JetSorter& sorter, const JetSelector selector) const { - /// @todo Will the vector be efficiently std::move'd by value through this function chain? return jets(selector, sorter); } /// Get the jets, ordered by \f$ p_T \f$, with optional cuts. /// /// @note Returns a copy rather than a reference, due to cuts and sorting /// /// This is a very common use-case, so is available as syntatic sugar for jets(c, cmpMomByPt). - /// @todo The other sorted accessors should be removed in a cleanup. Jets jetsByPt(const Cut& c=Cuts::open()) const { return jets(c, cmpMomByPt); } /// Get the jets, ordered by \f$ p_T \f$, with cuts via a selection functor. /// /// @note Returns a copy rather than a reference, due to cuts and sorting /// /// This is a very common use-case, so is available as syntatic sugar for jets(c, cmpMomByPt). - /// @todo The other sorted accessors should be removed in a cleanup. Jets jetsByPt(const JetSelector& selector) const { return jets(selector, cmpMomByPt); } /// Get the jets, ordered by \f$ p_T \f$, with a cut on \f$ p_\perp \f$. /// /// @deprecated Use the version with a Cut argument /// @note Returns a copy rather than a reference, due to cuts and sorting /// /// This is a very common use-case, so is available as syntatic sugar for jets(Cuts::pT >= ptmin, cmpMomByPt). - /// @todo The other sorted accessors should be removed in a cleanup. - DEPRECATED("Use the version with a Cut argument") Jets jetsByPt(double ptmin) const { return jets(Cuts::pT >= ptmin, cmpMomByPt); } //@} protected: /// @brief Internal pure virtual method for getting jets in no guaranteed order. virtual Jets _jets() const = 0; public: /// Count the jets size_t size() const { return jets().size(); } /// Count the jets after a Cut is applied. size_t size(const Cut& c) const { return jets(c).size(); } /// Count the jets after a selection functor is applied. size_t size(const JetSelector& s) const { return jets(s).size(); } /// Is this jet finder empty? bool empty() const { return size() == 0; } /// Is this jet finder empty after a Cut is applied? bool empty(const Cut& c) const { return size(c) == 0; } /// Is this jet finder empty after a selection functor is applied? bool empty(const JetSelector& s) const { return size(s) == 0; } /// Clear the projection. virtual void reset() = 0; typedef Jet entity_type; typedef Jets collection_type; /// Template-usable interface common to FinalState. collection_type entities() const { return jets(); } // /// Do the calculation locally (no caching). // virtual void calc(const Particles& constituents, const Particles& tagparticles=Particles()) = 0; protected: /// Perform the projection on the Event. virtual void project(const Event& e) = 0; /// Compare projections. virtual CmpState compare(const Projection& p) const = 0; protected: /// Flag to determine whether or not to exclude (some) muons from the would-be constituents. Muons _useMuons; /// Flag to determine whether or not to exclude (some) invisible particles from the would-be constituents. Invisibles _useInvisibles; }; /// Compatibility typedef, for equivalence with ParticleFinder /// @todo Should we make this the canonical name? Would "require" a header filename change -> breakage or ugly. using JetFinder = JetAlg; } #endif diff --git a/src/Projections/InitialQuarks.cc b/src/Projections/InitialQuarks.cc --- a/src/Projections/InitialQuarks.cc +++ b/src/Projections/InitialQuarks.cc @@ -1,61 +1,62 @@ // -*- C++ -*- +#define I_KNOW_THE_INITIAL_QUARKS_PROJECTION_IS_DODGY_BUT_NEED_TO_USE_IT #include "Rivet/Projections/InitialQuarks.hh" namespace Rivet { CmpState InitialQuarks::compare(const Projection& p) const { return CmpState::EQ; } void InitialQuarks::project(const Event& e) { _theParticles.clear(); for (const GenParticle* p : Rivet::particles(e.genEvent())) { const GenVertex* pv = p->production_vertex(); const GenVertex* dv = p->end_vertex(); const PdgId pid = abs(p->pdg_id()); bool passed = inRange((long)pid, 1, 6); if (passed) { if (pv != 0) { for (const GenParticle* pp : particles_in(pv)) { // Only accept if parent is electron or Z0 const PdgId pid = abs(pp->pdg_id()); passed = (pid == PID::ELECTRON || abs(pp->pdg_id()) == PID::ZBOSON || abs(pp->pdg_id()) == PID::GAMMA); } } else { passed = false; } } if (getLog().isActive(Log::TRACE)) { const int st = p->status(); const double pT = p->momentum().perp(); const double eta = p->momentum().eta(); MSG_TRACE(std::boolalpha << "ID = " << p->pdg_id() << ", status = " << st << ", pT = " << pT << ", eta = " << eta << ": result = " << passed); if (pv != 0) { for (const GenParticle* pp : particles_in(pv)) { MSG_TRACE(std::boolalpha << " parent ID = " << pp->pdg_id()); } } if (dv != 0) { for (const GenParticle* pp : particles_out(dv)) { MSG_TRACE(std::boolalpha << " child ID = " << pp->pdg_id()); } } } if (passed) _theParticles.push_back(Particle(*p)); } MSG_DEBUG("Number of initial quarks = " << _theParticles.size()); if (!_theParticles.empty()) { for (size_t i = 0; i < _theParticles.size(); i++) { MSG_DEBUG("Initial quark[" << i << "] = " << _theParticles[i].pid()); } } } }