|
[Rivet-svn] r2113 - in trunk: . data/plotinfo include/Rivet include/Rivet/Tools src src/Analyses src/Projections src/Toolsblackhole at projects.hepforge.org blackhole at projects.hepforge.orgMon Nov 30 16:56:40 GMT 2009
Author: buckley Date: Mon Nov 30 16:56:40 2009 New Revision: 2113 Log: Many improvements to W analysis, addition of PID function wrappers which take Particle refs as arguments, new crossSectionPerEvent() function on Analysis Added: trunk/src/Particle.cc Modified: trunk/ChangeLog trunk/data/plotinfo/MC_LHC_WANALYSIS.plot trunk/include/Rivet/Analysis.hh trunk/include/Rivet/Particle.hh trunk/include/Rivet/Tools/ParticleIdUtils.hh trunk/src/Analyses/MC_LHC_WANALYSIS.cc trunk/src/Analysis.cc trunk/src/Makefile.am trunk/src/Projections/WFinder.cc trunk/src/Tools/ParticleIdUtils.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Mon Nov 30 12:57:23 2009 (r2112) +++ trunk/ChangeLog Mon Nov 30 16:56:40 2009 (r2113) @@ -1,3 +1,14 @@ +2009-11-30 Andy Buckley <andy at insectnation.org> + + * Adding crossSectionPerEvent() method [== + crossSection()/sumOfWeights()] to Analysis. Useful for histogram + scaling since numerator of sumW_passed/sumW_total (to calculate + pass-cuts xsec) is cancelled by dividing histo by sumW_passed. + + * Clean-up of Particle class and provision of inline PID:: + functions which take a Particle as an argument to avoid having to + explicitly call the Particle::pdgId() method. + 2009-11-30 Hendrik Hoeth <hendrik.hoeth at cern.ch> * Fixing division by zero in Profile1D bin errors for @@ -12,8 +23,8 @@ * Adding missing CDF_2001_S4751469 plots to uemerge * New "ShowZero" option in make-plots * Improving lots of plot defaults - * Fixing typos / non-existing bins in CDF_2001_S4751469 - and CDF_2004_S5839831 reference data + * Fixing typos / non-existing bins in CDF_2001_S4751469 and + CDF_2004_S5839831 reference data 2009-11-19 Hendrik Hoeth <hendrik.hoeth at cern.ch> @@ -21,8 +32,8 @@ 2009-11-17 Hendrik Hoeth <hendrik.hoeth at cern.ch> - * Zeroth version of STAR_2006_S6860818 analysis (identified strange - particles). Not working yet for unstable particles. + * Zeroth version of STAR_2006_S6860818 analysis (identified + strange particles). Not working yet for unstable particles. 2009-11-11 Andy Buckley <andy at insectnation.org> @@ -72,9 +83,9 @@ 2009-11-04 Hendrik Hoeth <hendrik.hoeth at cern.ch> - * New analysis STAR_2006_S6500200 (pion and proton pT spectra - in pp collisions at 200 GeV). It is still unclear if they used - a cut in rapidity or pseudorapidity, thus the analysis is declared + * New analysis STAR_2006_S6500200 (pion and proton pT spectra in + pp collisions at 200 GeV). It is still unclear if they used a cut + in rapidity or pseudorapidity, thus the analysis is declared "UNDER DEVELOPMENT" and "DO NOT USE". * Fixing compare() in NeutralFinalState and MergedFinalState @@ -89,12 +100,13 @@ 2009-11-02 Hendrik Hoeth <hendrik.hoeth at cern.ch> - * Fixing normalisation issue with stacked histograms in make-plots. + * Fixing normalisation issue with stacked histograms in + make-plots. 2009-10-30 Hendrik Hoeth <hendrik.hoeth at cern.ch> - * CDF_2009_S8233977: Updating data and axes labels to match - final paper. Normalise to cross-section instead of data. + * CDF_2009_S8233977: Updating data and axes labels to match final + paper. Normalise to cross-section instead of data. 2009-10-23 Andy Buckley <andy at insectnation.org> @@ -119,16 +131,16 @@ expected (E_T comparison still missing). * Adding new MergedFinalState projection. This merges two final - states, removing duplicate particles. Duplicates are identified - by looking at the genParticle(), so users need to take care of - any manually added particles themselves. + states, removing duplicate particles. Duplicates are identified by + looking at the genParticle(), so users need to take care of any + manually added particles themselves. * Fixing most open issues with the STAR_2009_UE_HELEN analysis. There is only one question left, regarding the away region. * Set the default split-merge value for SISCone in our FastJets - projection to the recommended (but not Fastjet-default!) value - of 0.75. + projection to the recommended (but not Fastjet-default!) value of + 0.75. 2009-10-17 Andy Buckley <andy at insectnation.org> @@ -424,8 +436,8 @@ * New analyses STAR_2006_S6870392 and STAR_2008_S7993412. In STAR_2008_S7993412 only the first distribution is filled at the moment. STAR_2006_S6870392 is normalised to data instead of the - Monte Carlo cross-section, since we don't have that available - in the HepMC stream yet. + Monte Carlo cross-section, since we don't have that available in + the HepMC stream yet. 2009-05-05 Andy Buckley <andy at insectnation.org> @@ -525,7 +537,8 @@ * Moving CDF_2001 constructor into implementation file. - * Cleaning up CDF_2008_LEADINGJETS a bit, e.g. using foreach loops. + * Cleaning up CDF_2008_LEADINGJETS a bit, e.g. using foreach + loops. * Adding function interface for specifying axis labels in histo bookings. Currently has no effect, since AIDA doesn't seem to have @@ -565,7 +578,8 @@ 2009-03-04 Andy Buckley <andy at insectnation.org> - * Fixing use of fastjet-config to not use the user's PATH variable. + * Fixing use of fastjet-config to not use the user's PATH + variable. * Fixing SWIG type table for HepMC object interchange. @@ -604,7 +618,8 @@ 2009-01-19 Andy Buckley <andy at insectnation.org> - * Added psyco Python optimiser to rivet, make-plots and compare-histos. + * Added psyco Python optimiser to rivet, make-plots and + compare-histos. * bin/aida2root: Added "-" -> "_" mangling, following requests. Modified: trunk/data/plotinfo/MC_LHC_WANALYSIS.plot ============================================================================== --- trunk/data/plotinfo/MC_LHC_WANALYSIS.plot Mon Nov 30 12:57:23 2009 (r2112) +++ trunk/data/plotinfo/MC_LHC_WANALYSIS.plot Mon Nov 30 16:56:40 2009 (r2113) @@ -2,82 +2,145 @@ Title=Charged multiplicity XLabel=$n_\text{ch}$ YLabel=$\mathrm{d}N_\text{evt}/\mathrm{d}n_\text{ch}$ +LogY=0 # END PLOT # BEGIN PLOT /MC_LHC_WANALYSIS/n-w Title=W candidate multiplicity -XLabel=$n_\text{W}$ +XLabel=$n_\text{W}^\pm$ YLabel=$\mathrm{d}N_\text{evt}/\mathrm{d}n_\text{W}$ +LogY=0 +# END PLOT + +# BEGIN PLOT /MC_LHC_WANALYSIS/n-wplus +Title=\PWplus candidate multiplicity +XLabel=$n_\text{W}^+$ +YLabel=$\mathrm{d}N_\text{evt}/\mathrm{d}n_\text{W}$ +LogY=0 +# END PLOT + +# BEGIN PLOT /MC_LHC_WANALYSIS/n-wminus +Title=\PWminus candidate multiplicity +XLabel=$n_\text{W}^-$ +YLabel=$\mathrm{d}N_\text{evt}/\mathrm{d}n_\text{W}$ +LogY=0 # END PLOT # BEGIN PLOT /MC_LHC_WANALYSIS/n-jet Title=Jet multiplicity XLabel=$n_\text{jet}$ YLabel=$\mathrm{d}N_\text{evt}/\mathrm{d}n_\text{jet}$ +LogY=0 # END PLOT # BEGIN PLOT /MC_LHC_WANALYSIS/pt-ch Title=Charged track $p_\perp$ XLabel=$p_\perp$ / GeV YLabel=$\mathrm{d}N_\text{ch}/\mathrm{d}p_\perp^\text{ch}$ / GeV$^{-1}$ +LogY=1 # END PLOT # BEGIN PLOT /MC_LHC_WANALYSIS/ptavg-ch Title=Average charged track $p_\perp$ XLabel=$\langle p_\perp^\text{ch} \rangle$ / GeV YLabel=$\mathrm{d}N_\text{ch}/\mathrm{d}\langle p_\perp^\text{ch} \langle$ / GeV$^{-1}$ +LogY=1 # END PLOT # BEGIN PLOT /MC_LHC_WANALYSIS/ptrms-ch Title=Width of charged track $p_\perp$ XLabel=$\hat{\sigma}_{p_\perp^\text{ch}}$ / GeV YLabel=$\mathrm{d}N_\text{ch}/\mathrm{d}\hat{\sigma}_{p_\perp^\text{ch}}$ / GeV$^{-1}$ +LogY=1 # END PLOT # BEGIN PLOT /MC_LHC_WANALYSIS/pt-w Title=W candidate $p_\perp$ XLabel=$p_\perp(\text{W})$ / GeV YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}p_\perp(\text{W})$ / GeV$^{-1}$ +LogY=1 +# END PLOT + +# BEGIN PLOT /MC_LHC_WANALYSIS/pt-wplus +Title=\PWplus candidate $p_\perp$ +XLabel=$p_\perp(\text{\PWplus})$ / GeV +YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}p_\perp(\text{W})$ / GeV$^{-1}$ +LogY=1 # END PLOT -# BEGIN PLOT /MC_LHC_WANALYSIS/logpt-w -Title=W candidate $\log(p_\perp/\text{GeV})$ -XLabel=$\log(p_\perp(\text{W})/\text{GeV})$ -YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}\log(p_\perp(\text{W})/\text{GeV})$ +# BEGIN PLOT /MC_LHC_WANALYSIS/pt-wminus +Title=\PWminus candidate $p_\perp$ +XLabel=$p_\perp(\text{\PWminus})$ / GeV +YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}p_\perp(\text{W})$ / GeV$^{-1}$ +LogY=1 # END PLOT # BEGIN PLOT /MC_LHC_WANALYSIS/pt-jet Title=Jet $p_\perp$ XLabel=$p_\perp(\text{jet})$ / GeV YLabel=$\mathrm{d}N_\text{jet}/\mathrm{d}p_\perp(\text{jet})$ / GeV$^{-1}$ -# END PLOT - -# BEGIN PLOT /MC_LHC_WANALYSIS/logpt-jet -Title=Jet $\log(p_\perp/\text{GeV})$ -XLabel=$\log(p_\perp(\text{jet})/\text{GeV})$ -YLabel=$\mathrm{d}N_\text{jet}/\mathrm{d}\log(p_\perp(\text{jet})/\text{GeV})$ +LogY=1 # END PLOT # BEGIN PLOT /MC_LHC_WANALYSIS/eta-w Title=W candidate $\eta$ XLabel=$\eta(\text{W})$ YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}\eta(\text{W})$ +LogY=0 +# END PLOT + +# BEGIN PLOT /MC_LHC_WANALYSIS/eta-wplus +Title=\PWplus candidate $\eta$ +XLabel=$\eta(\text{W})$ +YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}\eta(\text{W})$ +LogY=0 +# END PLOT + +# BEGIN PLOT /MC_LHC_WANALYSIS/eta-wminus +Title=\PWminus candidate $\eta$ +XLabel=$\eta(\text{W})$ +YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}\eta(\text{W})$ +LogY=0 # END PLOT # BEGIN PLOT /MC_LHC_WANALYSIS/phi-w Title=W candidate $\phi$ XLabel=$\phi(\text{W})$ YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}\phi(\text{W})$ +LogY=0 +# END PLOT + +# BEGIN PLOT /MC_LHC_WANALYSIS/phi-wplus +Title=\PWplus candidate $\phi$ +XLabel=$\phi(\text{W})$ +YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}\phi(\text{W})$ +LogY=0 +# END PLOT + +# BEGIN PLOT /MC_LHC_WANALYSIS/phi-wminus +Title=\PWminus candidate $\phi$ +XLabel=$\phi(\text{W})$ +YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}\phi(\text{W})$ +LogY=0 # END PLOT # BEGIN PLOT /MC_LHC_WANALYSIS/m-w Title=W candidate mass XLabel=$m(\text{W})$ / GeV YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}m(\text{W})$ / GeV$^{-1}$ +LogY=1 # END PLOT -# BEGIN PLOT /MC_LHC_WANALYSIS/m-w -Title=Log of W candidate mass -XLabel=$\log(m(\text{W})/\text{GeV})$ -YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}\log(m(\text{W})/\text{GeV})$ +# BEGIN PLOT /MC_LHC_WANALYSIS/m-wplus +Title=\PWplus candidate mass +XLabel=$m(\text{W})$ / GeV +YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}m(\text{W})$ / GeV$^{-1}$ +LogY=1 +# END PLOT + +# BEGIN PLOT /MC_LHC_WANALYSIS/m-wminus +Title=\PWminus candidate mass +XLabel=$m(\text{W})$ / GeV +YLabel=$\mathrm{d}N_\text{W}/\mathrm{d}m(\text{W})$ / GeV$^{-1}$ +LogY=1 # END PLOT Modified: trunk/include/Rivet/Analysis.hh ============================================================================== --- trunk/include/Rivet/Analysis.hh Mon Nov 30 12:57:23 2009 (r2112) +++ trunk/include/Rivet/Analysis.hh Mon Nov 30 16:56:40 2009 (r2113) @@ -193,8 +193,11 @@ protected: /// Get the process cross-section in pb. Throws if this hasn't been set. - const double& crossSection() const; + double crossSection() const; + /// Get the process cross-section per generated event in pb. Throws if this + /// hasn't been set. + double crossSectionPerEvent() const; /// Get a Log object based on the name() property of the calling analysis object. Log& getLog() const; Modified: trunk/include/Rivet/Particle.hh ============================================================================== --- trunk/include/Rivet/Particle.hh Mon Nov 30 12:57:23 2009 (r2112) +++ trunk/include/Rivet/Particle.hh Mon Nov 30 16:56:40 2009 (r2113) @@ -5,7 +5,9 @@ #include "Rivet/Rivet.hh" #include "Rivet/Particle.fhh" #include "Rivet/ParticleBase.hh" +#include "Rivet/ParticleName.hh" #include "Rivet/Math/Vectors.hh" +#include "Rivet/Tools/Logging.hh" namespace Rivet { @@ -15,17 +17,23 @@ public: /// Default constructor. + /// @deprecated A particle without info is useless. This only exists to keep STL containers happy. Particle() : ParticleBase(), - _original(0), _id(0), _momentum(), _mass(0.0) + _original(0), _id(0), _momentum() { } + /// Constructor without GenParticle. + Particle(PdgId pid, const FourMomentum& mom) : ParticleBase(), + _original(0), _id(pid), _momentum(mom) + { } /// Constructor from a HepMC GenParticle. Particle(const GenParticle& gp) : ParticleBase(), _original(&gp), _id(gp.pdg_id()), - _momentum(gp.momentum()), _mass(gp.momentum().m()) + _momentum(gp.momentum()) { } + public: /// Get a const reference to the original GenParticle. @@ -36,41 +44,48 @@ /// Check if the particle corresponds to a GenParticle. - bool hasGenParticle() const { return bool(_original); } + bool hasGenParticle() const { + return bool(_original); + } /// The PDG ID code for this Particle. - long pdgId() const { return _id; } - - - /// The momentum of this Particle. - FourMomentum& momentum() { return _momentum; } + long pdgId() const { + return _id; + } /// The momentum of this Particle. - const FourMomentum& momentum() const { return _momentum; } + const FourMomentum& momentum() const { + return _momentum; + } /// Set the momentum of this Particle. - Particle& setMomentum(const FourMomentum& momentum) { _momentum = momentum; return *this; } + Particle& setMomentum(const FourMomentum& momentum) { + _momentum = momentum; + return *this; + } /// The mass of this Particle. - double mass() const { return _mass; } + double mass() const { + return momentum().mass(); + } + // /// The charge of this Particle. + // double charge() const { + // return PID::charge(*this); + // } + + // /// Three times the charge of this Particle (i.e. integer multiple of smallest quark charge). + // int threeCharge() const { + // return PID::threeCharge(*this); + // } - bool hasAncestor(long pdg_id) const { - GenVertex* prodVtx = genParticle().production_vertex(); - if (prodVtx == 0) return false; - GenVertex::particle_iterator ibegin = prodVtx->particles_begin(HepMC::ancestors); - GenVertex::particle_iterator iend = prodVtx->particles_end(HepMC::ancestors); - for (GenVertex::particle_iterator ancestor = ibegin; ancestor != iend; ++ancestor) { - if ((*ancestor)->pdg_id() == pdg_id) { - return true; - } - } - return false; - } + + /// Check whether a given PID is found in the GenParticle's ancestor list + bool hasAncestor(PdgId pdg_id) const; private: @@ -83,9 +98,6 @@ /// The momentum of this projection of the Particle. FourMomentum _momentum; - - /// The mass of this Particle, stored for numerical hygiene. - double _mass; }; Modified: trunk/include/Rivet/Tools/ParticleIdUtils.hh ============================================================================== --- trunk/include/Rivet/Tools/ParticleIdUtils.hh Mon Nov 30 12:57:23 2009 (r2112) +++ trunk/include/Rivet/Tools/ParticleIdUtils.hh Mon Nov 30 16:56:40 2009 (r2113) @@ -1,7 +1,7 @@ // ---------------------------------------------------------------------- // // ParticleIDMethods.hh -// Author: Lynn Garren +// Author: Lynn Garren, Andy Buckley // // various utilities to extract information from the particle ID // @@ -17,88 +17,152 @@ #ifndef RIVET_PARTICLE_ID_METHODS_HH #define RIVET_PARTICLE_ID_METHODS_HH +#include "Rivet/Particle.hh" + + namespace Rivet { + namespace PID { -/// PID digits (base 10) are: n nr nl nq1 nq2 nq3 nj -/// The location enum provides a convenient index into the PID. -enum location { nj=1, nq3, nq2, nq1, nl, nr, n, n8, n9, n10 }; - -/// return the digit at a named location in the PID -unsigned short digit( location loc, const int & pid ); - -/// if this is a nucleus (ion), get A -/// Ion numbers are +/- 10LZZZAAAI. -int A(const int & pid ); - -/// if this is a nucleus (ion), get Z -/// Ion numbers are +/- 10LZZZAAAI. -int Z(const int & pid ); - -/// if this is a nucleus (ion), get nLambda -/// Ion numbers are +/- 10LZZZAAAI. -int lambda( const int & pid ); - -/// absolute value of particle ID -int abspid( const int & pid ); - -/// extract fundamental ID (1-100) if this is a "fundamental" particle -int fundamentalID( const int & pid ); -/// if this is a fundamental particle, does it have a valid antiparticle? -//bool hasFundamentalAnti( const int & pid ); - -/// returns everything beyond the 7th digit -/// (e.g. outside the standard numbering scheme) -int extraBits( const int & pid ); -// --- boolean methods: -// -/// is this a valid ID? -bool isValid( const int & pid ); -/// is this a valid meson ID? -bool isMeson( const int & pid ); -/// is this a valid baryon ID? -bool isBaryon( const int & pid ); -/// is this a valid diquark ID? -bool isDiQuark( const int & pid ); -/// is this a valid hadron ID? -bool isHadron( const int & pid ); -/// is this a valid lepton ID? -bool isLepton( const int & pid ); -/// is this a valid ion ID? -bool isNucleus( const int & pid ); -/// is this a valid pentaquark ID? -bool isPentaquark( const int & pid ); -/// is this a valid SUSY ID? -bool isSUSY( const int & pid ); -/// is this a valid R-hadron ID? -bool isRhadron( const int & pid ); - -/// does this particle contain an up quark? -bool hasUp( const int & pid ); -/// does this particle contain a down quark? -bool hasDown( const int & pid ); -/// does this particle contain a strange quark? -bool hasStrange( const int & pid ); -/// does this particle contain a charm quark? -bool hasCharm( const int & pid ); -/// does this particle contain a bottom quark? -bool hasBottom( const int & pid ); -/// does this particle contain a top quark? -bool hasTop( const int & pid ); + /// @name PID operations on Rivet::Particle wrapper + //@{ -// --- other information -// -/// jSpin returns 2J+1, where J is the total spin -int jSpin( const int & pid ); -/// sSpin returns 2S+1, where S is the spin -int sSpin( const int & pid ); -/// lSpin returns 2L+1, where L is the orbital angular momentum -int lSpin( const int & pid ); -/// return 3 times the charge (3 x quark charge is an int) -int threeCharge( const int & pid ); + // /// if this is a nucleus (ion), get A + // /// Ion numbers are +/- 10LZZZAAAI. + // int A(const int & pid ); + + // /// if this is a nucleus (ion), get Z + // /// Ion numbers are +/- 10LZZZAAAI. + // int Z(const int & pid ); + + // /// if this is a nucleus (ion), get nLambda + // /// Ion numbers are +/- 10LZZZAAAI. + // int lambda( const int & pid ); + + /// absolute value of particle ID + int abspid( const int & pid ); + + + /// is this a valid ID? + bool isValid( const int & pid ); + /// is this a valid meson ID? + bool isMeson( const int & pid ); + /// is this a valid baryon ID? + bool isBaryon( const int & pid ); + /// is this a valid diquark ID? + bool isDiQuark( const int & pid ); + /// is this a valid hadron ID? + bool isHadron( const int & pid ); + /// is this a valid lepton ID? + bool isLepton( const int & pid ); + /// is this a valid ion ID? + bool isNucleus( const int & pid ); + /// is this a valid pentaquark ID? + bool isPentaquark( const int & pid ); + /// is this a valid SUSY ID? + bool isSUSY( const int & pid ); + /// is this a valid R-hadron ID? + bool isRhadron( const int & pid ); + + /// does this particle contain an up quark? + bool hasUp( const int & pid ); + /// does this particle contain a down quark? + bool hasDown( const int & pid ); + /// does this particle contain a strange quark? + bool hasStrange( const int & pid ); + /// does this particle contain a charm quark? + bool hasCharm( const int & pid ); + /// does this particle contain a bottom quark? + bool hasBottom( const int & pid ); + /// does this particle contain a top quark? + bool hasTop( const int & pid ); + + /// jSpin returns 2J+1, where J is the total spin + int jSpin( const int & pid ); + /// sSpin returns 2S+1, where S is the spin + int sSpin( const int & pid ); + /// lSpin returns 2L+1, where L is the orbital angular momentum + int lSpin( const int & pid ); + + /// return 3 times the charge (3 x quark charge is an int) + int threeCharge( const int & pid ); + /// return 3 times the charge (3 x quark charge is an int) + inline double charge( const int & pid ) { return 3.0 * threeCharge(pid); } + + //@} + + + ///////////////////////// + + + + /// @name PID operations on Rivet::Particle wrapper + //@{ + + /// if this is a nucleus (ion), get A + /// Ion numbers are +/- 10LZZZAAAI. + // int A(const Particle& p) { return A(p.pdgId()); } + + /// if this is a nucleus (ion), get Z + /// Ion numbers are +/- 10LZZZAAAI. + // int Z(const Particle& p) { return Z(p.pdgId()); } + + /// if this is a nucleus (ion), get nLambda + /// Ion numbers are +/- 10LZZZAAAI. + // int lambda( const Particle& p) { return lambda(p.pdgId()); } + + /// absolute value of particle ID + inline int abspid( const Particle& p) { return abspid(p.pdgId()); } + + /// is this a valid meson ID? + inline bool isMeson( const Particle& p ) { return isMeson(p.pdgId()); } + /// is this a valid baryon ID? + inline bool isBaryon( const Particle& p ) { return isBaryon(p.pdgId()); } + /// is this a valid diquark ID? + inline bool isDiQuark( const Particle& p ) { return isDiQuark(p.pdgId()); } + /// is this a valid hadron ID? + inline bool isHadron( const Particle& p ) { return isHadron(p.pdgId()); } + /// is this a valid lepton ID? + inline bool isLepton( const Particle& p ) { return isLepton(p.pdgId()); } + /// is this a valid ion ID? + inline bool isNucleus( const Particle& p ) { return isNucleus(p.pdgId()); } + /// is this a valid pentaquark ID? + inline bool isPentaquark( const Particle& p ) { return isPentaquark(p.pdgId()); } + /// is this a valid SUSY ID? + inline bool isSUSY( const Particle& p ) { return isSUSY(p.pdgId()); } + /// is this a valid R-hadron ID? + inline bool isRhadron( const Particle& p ) { return isRhadron(p.pdgId()); } + + /// does this particle contain an up quark? + inline bool hasUp( const Particle& p ) { return hasUp(p.pdgId()); } + /// does this particle contain a down quark? + inline bool hasDown( const Particle& p ) { return hasDown(p.pdgId()); } + /// does this particle contain a strange quark? + inline bool hasStrange( const Particle& p ) { return hasStrange(p.pdgId()); } + /// does this particle contain a charm quark? + inline bool hasCharm( const Particle& p ) { return hasCharm(p.pdgId()); } + /// does this particle contain a bottom quark? + inline bool hasBottom( const Particle& p ) { return hasBottom(p.pdgId()); } + /// does this particle contain a top quark? + inline bool hasTop( const Particle& p ) { return hasTop(p.pdgId()); } + + /// jSpin returns 2J+1, where J is the total spin + inline int jSpin( const Particle& p ) { return jSpin(p.pdgId()); } + /// sSpin returns 2S+1, where S is the spin + inline int sSpin( const Particle& p ) { return sSpin(p.pdgId()); } + /// lSpin returns 2L+1, where L is the orbital angular momentum + inline int lSpin( const Particle& p ) { return lSpin(p.pdgId()); } + + /// return 3 times the charge (3 x quark charge is an int) + inline int threeCharge( const Particle& p ) { return threeCharge(p.pdgId()); } + /// return 3 times the charge (3 x quark charge is an int) + inline double charge( const Particle& p ) { return 3.0 * threeCharge(p); } + + //@} + } -}} +} #endif Modified: trunk/src/Analyses/MC_LHC_WANALYSIS.cc ============================================================================== --- trunk/src/Analyses/MC_LHC_WANALYSIS.cc Mon Nov 30 12:57:23 2009 (r2112) +++ trunk/src/Analyses/MC_LHC_WANALYSIS.cc Mon Nov 30 16:56:40 2009 (r2113) @@ -15,7 +15,7 @@ /// Default constructor MC_LHC_WANALYSIS() : Analysis("MC_LHC_WANALYSIS") { - // + _sumWPass = 0.0; } @@ -27,36 +27,53 @@ const ChargedFinalState cfs; addProjection(cfs, "CFS"); /// @todo Handle muon-decay Ws as well - const WFinder wf(-MAXRAPIDITY, MAXRAPIDITY, 0.0*GeV, ELECTRON, 30.0*GeV, 110.0*GeV, 0.2); + const WFinder wf(-MAXRAPIDITY, MAXRAPIDITY, 0.0*GeV, ELECTRON, 70.0*GeV, 90.0*GeV, 0.2); addProjection(wf, "WF"); FastJets fastjets(wf.remainingFinalState(), FastJets::KT, 0.7); addProjection(fastjets, "Jets"); - _hist_chargemulti = bookHistogram1D("n-ch", 50, -0.5, 199.5); + _hist_chargemulti = bookHistogram1D("n-ch", 50, -0.5, 399.5); _hist_chargept = bookHistogram1D("pt-ch", 25, 0, 25); /// @todo Use profile plots instead: - _hist_chargemeanpt = bookHistogram1D("ptavg-ch", 25, 0, 10); + _hist_chargemeanpt = bookHistogram1D("ptavg-ch", 20, 0, 10); /// @todo Use profile plots instead (and this isn't really an "RMS") - _hist_chargermspt = bookHistogram1D("ptrms-ch", 25, 0, 10); - _hist_wcount = bookHistogram1D("n-w", 16, -0.5, 15.5); - _hist_wpt = bookHistogram1D("pt-w", 25, 0, 25); - _hist_wlogpt = bookHistogram1D("logpt-w", 30, 0, 6); + _hist_chargermspt = bookHistogram1D("ptrms-ch", 20, 0, 10); + // + _hist_wcount = bookHistogram1D("n-w", 6, -0.5, 5.5); + _hist_wpluscount = bookHistogram1D("n-wplus", 6, -0.5, 5.5); + _hist_wminuscount = bookHistogram1D("n-wminus", 6, -0.5, 5.5); + _hist_wpt = bookHistogram1D("pt-w", 25, 0, 100); + _hist_wpluspt = bookHistogram1D("pt-wplus", 25, 0, 100); + _hist_wminuspt = bookHistogram1D("pt-wminus", 25, 0, 100); _hist_weta = bookHistogram1D("eta-w", 36, -6, 6); + _hist_wpluseta = bookHistogram1D("eta-wplus", 36, -6, 6); + _hist_wminuseta = bookHistogram1D("eta-wminus", 36, -6, 6); _hist_wphi = bookHistogram1D("phi-w", 25, 0, TWOPI); + _hist_wplusphi = bookHistogram1D("phi-wplus", 25, 0, TWOPI); + _hist_wminusphi = bookHistogram1D("phi-wminus", 25, 0, TWOPI); _hist_wmass = bookHistogram1D("m-w", 40, 60, 100); - _hist_wlogmass = bookHistogram1D("logm-w", 20, 0, 10); - _hist_jetcount = bookHistogram1D("n-jet", 16, -0.5, 15.5); + _hist_wplusmass = bookHistogram1D("m-wplus", 40, 60, 100); + _hist_wminusmass = bookHistogram1D("m-wminus", 40, 60, 100); + //_hist_weta_asymm = bookProfile1D("asymm-eta-w", 20, -5.0, 5.0); + // + _hist_jetcount = bookHistogram1D("n-jet", 11, -0.5, 10.5); _hist_jetpt = bookHistogram1D("pt-jet", 50, 20, 100); - _hist_jetlogpt = bookHistogram1D("logpt-jet", 20, 0, 20); } void analyze(const Event& event) { - const double weight = event.weight(); const FinalState& cfs = applyProjection<FinalState>(event, "CFS"); const WFinder& wf = applyProjection<WFinder>(event, "WF"); const FastJets& fastjets = applyProjection<FastJets>(event, "Jets"); const Jets jets = fastjets.jetsByPt(); + + // Veto if no Ws found + if (wf.size() == 0) { + getLog() << Log::DEBUG << "No W candidates found: vetoing" << endl; + vetoEvent; + } + const double weight = event.weight(); + _sumWPass += weight; // Charged particles part _hist_chargemulti->fill(cfs.particles().size(), weight); @@ -73,54 +90,102 @@ _hist_chargermspt->fill(rmspt/GeV, weight); // W part - _hist_wcount->fill(wf.particles().size(), weight); + uint n_wplus(0), n_wminus(0); foreach (const Particle& wp, wf.particles()) { const double pT = wp.momentum().pT(); - _hist_wpt->fill(pT/GeV, weight); - _hist_wlogpt->fill(log(pT/GeV), weight); - _hist_weta->fill(wp.momentum().pseudorapidity(), weight); - _hist_wphi->fill(wp.momentum().azimuthalAngle(), weight); + const double eta = wp.momentum().eta(); + const double phi = wp.momentum().phi(); const double m = wp.momentum().mass(); + /// @todo When histo handling is easier, build total histos by summing W+/- histos + _hist_wpt->fill(pT/GeV, weight); + _hist_weta->fill(eta, weight); + _hist_wphi->fill(phi, weight); _hist_wmass->fill(m/GeV, weight); - _hist_wlogmass->fill(log(m/GeV), weight); + if (wp.pdgId() == WPLUSBOSON) { + n_wplus += 1; + _hist_wpluspt->fill(pT/GeV, weight); + _hist_wpluseta->fill(eta, weight); + _hist_wplusphi->fill(phi, weight); + _hist_wplusmass->fill(m/GeV, weight); + } else if (wp.pdgId() == WMINUSBOSON) { + n_wminus += 1; + _hist_wminuspt->fill(pT/GeV, weight); + _hist_wminuseta->fill(eta, weight); + _hist_wminusphi->fill(phi, weight); + _hist_wminusmass->fill(m/GeV, weight); + } else { + // Just checking! + throw Error("There shouldn't be any W candidates without a W PID!"); + } } - + _hist_wcount->fill(n_wplus+n_wminus, weight); + _hist_wpluscount->fill(n_wplus, weight); + _hist_wminuscount->fill(n_wminus, weight); + // Jet part - _hist_jetcount->fill(fastjets.size(), weight); + _hist_jetcount->fill(fastjets.jets().size(), weight); foreach(const Jet& j, fastjets.jetsByPt()) { const double pT = j.momentum().pT(); _hist_jetpt->fill(pT/GeV, weight); - _hist_jetlogpt->fill(log(pT/GeV), weight); } + } void finalize() { - ///@todo Obtain cross-sections from generator and normalise histos to them. + const double xsec_sumw = crossSectionPerEvent()/picobarn; + /// @todo Actually "measure" separate W+ and W- xsecs + const double xsec_sumw_plus = xsec_sumw/2.0; + const double xsec_sumw_minus = xsec_sumw/2.0; + + scale(_hist_chargemulti, xsec_sumw); + scale(_hist_chargept, xsec_sumw); + scale(_hist_chargemeanpt, xsec_sumw); + scale(_hist_chargermspt, xsec_sumw); + scale(_hist_wcount, xsec_sumw); + scale(_hist_wpluscount, xsec_sumw_plus); + scale(_hist_wminuscount, xsec_sumw_minus); + scale(_hist_wpt, xsec_sumw); + scale(_hist_wpluspt, xsec_sumw_plus); + scale(_hist_wminuspt, xsec_sumw_minus); + scale(_hist_weta, xsec_sumw); + scale(_hist_wpluseta, xsec_sumw_plus); + scale(_hist_wminuseta, xsec_sumw_minus); + scale(_hist_wphi, xsec_sumw); + scale(_hist_wplusphi, xsec_sumw_plus); + scale(_hist_wminusphi, xsec_sumw_minus); + scale(_hist_wmass, xsec_sumw); + scale(_hist_wplusmass, xsec_sumw_plus); + scale(_hist_wminusmass, xsec_sumw_minus); + scale(_hist_jetcount, xsec_sumw); + scale(_hist_jetpt, xsec_sumw); } //@} + private: + + /// @name Counters + //@{ + double _sumWPass; + //@} + + /// @name Histograms //@{ AIDA::IHistogram1D* _hist_chargemulti; AIDA::IHistogram1D* _hist_chargept; AIDA::IHistogram1D* _hist_chargemeanpt; AIDA::IHistogram1D* _hist_chargermspt; - AIDA::IHistogram1D* _hist_wcount; - AIDA::IHistogram1D* _hist_wpt; - AIDA::IHistogram1D* _hist_wlogpt; - //AIDA::IHistogram1D* _hist_wpthigh; - //AIDA::IHistogram1D* _hist_wlogpthigh; - AIDA::IHistogram1D* _hist_weta; - AIDA::IHistogram1D* _hist_wphi; - AIDA::IHistogram1D* _hist_wmass; - AIDA::IHistogram1D* _hist_wlogmass; + AIDA::IHistogram1D *_hist_wcount, *_hist_wpluscount, *_hist_wminuscount; + AIDA::IHistogram1D *_hist_wpt, *_hist_wpluspt, *_hist_wminuspt; + AIDA::IHistogram1D *_hist_weta, *_hist_wpluseta, *_hist_wminuseta; + AIDA::IHistogram1D *_hist_wphi, *_hist_wplusphi, *_hist_wminusphi; + AIDA::IHistogram1D *_hist_wmass, *_hist_wplusmass, *_hist_wminusmass; AIDA::IHistogram1D* _hist_jetcount; AIDA::IHistogram1D* _hist_jetpt; - AIDA::IHistogram1D* _hist_jetlogpt; //@} }; Modified: trunk/src/Analysis.cc ============================================================================== --- trunk/src/Analysis.cc Mon Nov 30 12:57:23 2009 (r2112) +++ trunk/src/Analysis.cc Mon Nov 30 16:56:40 2009 (r2113) @@ -182,7 +182,7 @@ return *this; } - const double& Analysis::crossSection() const { + double Analysis::crossSection() const { if (!_gotCrossSection) { string errMsg = "You did not set the cross section for the analysis " + name(); throw Error(errMsg); @@ -190,6 +190,12 @@ return _crossSection; } + double Analysis::crossSectionPerEvent() const { + const double sumW = sumOfWeights(); + assert(sumW > 0); + return _crossSection / sumW; + } + AnalysisHandler& Analysis::handler() const { return *_analysishandler; } Modified: trunk/src/Makefile.am ============================================================================== --- trunk/src/Makefile.am Mon Nov 30 12:57:23 2009 (r2112) +++ trunk/src/Makefile.am Mon Nov 30 16:56:40 2009 (r2113) @@ -7,7 +7,7 @@ lib_LTLIBRARIES = libRivet.la libRivet_la_SOURCES = \ - Event.cc Jet.cc \ + Event.cc Jet.cc Particle.cc \ ProjectionApplier.cc Projection.cc \ Analysis.cc AnalysisLoader.cc AnalysisInfo.cc \ AnalysisHandler.cc Run.cc ProjectionHandler.cc HistoHandler.cc Added: trunk/src/Particle.cc ============================================================================== --- /dev/null 00:00:00 1970 (empty, because file is newly added) +++ trunk/src/Particle.cc Mon Nov 30 16:56:40 2009 (r2113) @@ -0,0 +1,16 @@ +#include "Rivet/Particle.hh" + +namespace Rivet { + + + bool Particle::hasAncestor(PdgId pdg_id) const { + GenVertex* prodVtx = genParticle().production_vertex(); + if (prodVtx == 0) return false; + foreach (const GenParticle* ancestor, particles(prodVtx, HepMC::ancestors)) { + if (ancestor->pdg_id() == pdg_id) return true; + } + return false; + } + + +} Modified: trunk/src/Projections/WFinder.cc ============================================================================== --- trunk/src/Projections/WFinder.cc Mon Nov 30 12:57:23 2009 (r2112) +++ trunk/src/Projections/WFinder.cc Mon Nov 30 16:56:40 2009 (r2113) @@ -128,10 +128,9 @@ msg << " = " << pW; getLog() << Log::DEBUG << msg.str() << endl; - Particle W; - W.setMomentum(pW); - - _theParticles.push_back(W); + // Make W Particle and insert into particles list + const PdgId wpid = (wcharge == 1) ? WPLUSBOSON : WPLUSBOSON; + _theParticles.push_back(Particle(wpid, pW)); getLog() << Log::DEBUG << name() << " found " << _theParticles.size() << " W candidates." << endl; } Modified: trunk/src/Tools/ParticleIdUtils.cc ============================================================================== --- trunk/src/Tools/ParticleIdUtils.cc Mon Nov 30 12:57:23 2009 (r2112) +++ trunk/src/Tools/ParticleIdUtils.cc Mon Nov 30 16:56:40 2009 (r2113) @@ -8,453 +8,473 @@ #include <cmath> // for pow() #include "Rivet/Tools/ParticleIdUtils.hh" -//#include "Rivet/Tools/ParticleName.hh" namespace Rivet { + namespace PID { -// absolute value -int abspid( const int & pid ) -{ - return (pid < 0) ? -pid : pid; -} -// returns everything beyond the 7th digit (e.g. outside the numbering scheme) -int extraBits( const int & pid ) -{ - return abspid(pid)/10000000; -} + /// PID digits (base 10) are: n nr nl nq1 nq2 nq3 nj + /// The location enum provides a convenient index into the PID. + enum location { nj=1, nq3, nq2, nq1, nl, nr, n, n8, n9, n10 }; + + /// return the digit at a named location in the PID + unsigned short digit( location loc, const int & pid ); + + /// extract fundamental ID (1-100) if this is a "fundamental" particle + int fundamentalID( const int & pid ); + /// if this is a fundamental particle, does it have a valid antiparticle? + //bool hasFundamentalAnti( const int & pid ); + + /// returns everything beyond the 7th digit + /// (e.g. outside the standard numbering scheme) + int extraBits( const int & pid ); + + + // absolute value + int abspid( const int & pid ) + { + return (pid < 0) ? -pid : pid; + } -// split the PID into constituent integers -unsigned short digit( location loc, const int & pid ) -{ - // PID digits (base 10) are: n nr nl nq1 nq2 nq3 nj - // the location enum provides a convenient index into the PID - int numerator = (int) std::pow(10.0,(loc-1)); - return (abspid(pid)/numerator)%10; -} + // returns everything beyond the 7th digit (e.g. outside the numbering scheme) + int extraBits( const int & pid ) + { + return abspid(pid)/10000000; + } -// return the first two digits if this is a "fundamental" particle -// ID = 100 is a special case (internal generator ID's are 81-100) -int fundamentalID( const int & pid ) -{ - if( extraBits(pid) > 0 ) return 0; - if( digit(nq2,pid) == 0 && digit(nq1,pid) == 0) { - return abspid(pid)%10000; - } else if( abspid(pid) <= 100 ) { - return abspid(pid); - } else { - return 0; + // split the PID into constituent integers + unsigned short digit( location loc, const int & pid ) + { + // PID digits (base 10) are: n nr nl nq1 nq2 nq3 nj + // the location enum provides a convenient index into the PID + int numerator = (int) std::pow(10.0,(loc-1)); + return (abspid(pid)/numerator)%10; } -} -// Ion numbers are +/- 10LZZZAAAI. -int Z( const int & pid ) -{ - // a proton can also be a Hydrogen nucleus - if( abspid(pid) == 2212 ) { return 1; } - if( isNucleus(pid) ) return (abspid(pid)/10000)%1000; - return 0; -} + // return the first two digits if this is a "fundamental" particle + // ID = 100 is a special case (internal generator ID's are 81-100) + int fundamentalID( const int & pid ) + { + if( extraBits(pid) > 0 ) return 0; + if( digit(nq2,pid) == 0 && digit(nq1,pid) == 0) { + return abspid(pid)%10000; + } else if( abspid(pid) <= 100 ) { + return abspid(pid); + } else { + return 0; + } + } -// Ion numbers are +/- 10LZZZAAAI. -int A( const int & pid ) -{ - // a proton can also be a Hydrogen nucleus - if( abspid(pid) == 2212 ) { return 1; } - if( isNucleus(pid) ) return (abspid(pid)/10)%1000; - return 0; -} + // Ion numbers are +/- 10LZZZAAAI. + int Z( const int & pid ) + { + // a proton can also be a Hydrogen nucleus + if( abspid(pid) == 2212 ) { return 1; } + if( isNucleus(pid) ) return (abspid(pid)/10000)%1000; + return 0; + } -// if this is a nucleus (ion), get nLambda -// Ion numbers are +/- 10LZZZAAAI. -int lambda( const int & pid ) -{ - // a proton can also be a Hydrogen nucleus - if( abspid(pid) == 2212 ) { return 0; } - if( isNucleus(pid) ) return digit(n8,pid); - return 0; -} + // Ion numbers are +/- 10LZZZAAAI. + int A( const int & pid ) + { + // a proton can also be a Hydrogen nucleus + if( abspid(pid) == 2212 ) { return 1; } + if( isNucleus(pid) ) return (abspid(pid)/10)%1000; + return 0; + } + // if this is a nucleus (ion), get nLambda + // Ion numbers are +/- 10LZZZAAAI. + int lambda( const int & pid ) + { + // a proton can also be a Hydrogen nucleus + if( abspid(pid) == 2212 ) { return 0; } + if( isNucleus(pid) ) return digit(n8,pid); + return 0; + } -// --- boolean methods: -// -// check to see if this is a valid PID -bool isValid( const int & pid ) -{ - if( extraBits(pid) > 0 ) { - if( isNucleus(pid) ) { return true; } - return false; - } - if( isSUSY(pid) ) { return true; } - if( isRhadron(pid) ) { return true; } - // Meson signature - if( isMeson(pid) ) { return true; } - // Baryon signature - if( isBaryon(pid) ) { return true; } - // DiQuark signature - if( isDiQuark(pid) ) { return true; } - // fundamental particle - if( fundamentalID(pid) > 0 ) { - if(pid > 0 ) { - return true; - } else { - // AB - disabled this to remove need for PID -> name lookup. - //if( hasFundamentalAnti(pid) ) { return true; } - return false; - } - } - // pentaquark - if( isPentaquark(pid) ) { return true; } - // don't recognize this number - return false; -} + // --- boolean methods: + // -// // if this is a fundamental particle, does it have a valid antiparticle? -// bool hasFundamentalAnti( const int & pid ) -// { -// // these are defined by the generator and therefore are always valid -// if( fundamentalID(pid) <= 100 && fundamentalID(pid) >= 80 ) { return true; } -// // check id's from 1 to 79 -// if( fundamentalID(pid) > 0 && fundamentalID(pid) < 80 ) { -// if( validParticleName(-pid) ) { return true; } -// } -// return false; -// } - -// check to see if this is a valid meson -bool isMeson( const int & pid ) -{ - if( extraBits(pid) > 0 ) { return false; } - if( abspid(pid) <= 100 ) { return false; } - if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) { return false; } - int aid = abspid(pid); - if( aid == 130 || aid == 310 || aid == 210 ) { return true; } - // EvtGen uses some odd numbers - if( aid == 150 || aid == 350 || aid == 510 || aid == 530 ) { return true; } - // pomeron, etc. - if( pid == 110 || pid == 990 || pid == 9990 ) { return true; } - if( digit(nj,pid) > 0 && digit(nq3,pid) > 0 - && digit(nq2,pid) > 0 && digit(nq1,pid) == 0 ) { - // check for illegal antiparticles - if( digit(nq3,pid) == digit(nq2,pid) && pid < 0 ) { + // check to see if this is a valid PID + bool isValid( const int & pid ) + { + if( extraBits(pid) > 0 ) { + if( isNucleus(pid) ) { return true; } return false; - } else { + } + if( isSUSY(pid) ) { return true; } + if( isRhadron(pid) ) { return true; } + // Meson signature + if( isMeson(pid) ) { return true; } + // Baryon signature + if( isBaryon(pid) ) { return true; } + // DiQuark signature + if( isDiQuark(pid) ) { return true; } + // fundamental particle + if( fundamentalID(pid) > 0 ) { + if(pid > 0 ) { return true; + } else { + // AB - disabled this to remove need for PID -> name lookup. + //if( hasFundamentalAnti(pid) ) { return true; } + return false; + } } + // pentaquark + if( isPentaquark(pid) ) { return true; } + // don't recognize this number + return false; } - return false; -} -// check to see if this is a valid baryon -bool isBaryon( const int & pid ) -{ - if( extraBits(pid) > 0 ) { return false; } - if( abspid(pid) <= 100 ) { return false; } - if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) { return false; } - if( abspid(pid) == 2110 || abspid(pid) == 2210 ) { return true; } - if( digit(nj,pid) > 0 && digit(nq3,pid) > 0 - && digit(nq2,pid) > 0 && digit(nq1,pid) > 0 ) { return true; } - return false; -} + // // if this is a fundamental particle, does it have a valid antiparticle? + // bool hasFundamentalAnti( const int & pid ) + // { + // // these are defined by the generator and therefore are always valid + // if( fundamentalID(pid) <= 100 && fundamentalID(pid) >= 80 ) { return true; } + // // check id's from 1 to 79 + // if( fundamentalID(pid) > 0 && fundamentalID(pid) < 80 ) { + // if( validParticleName(-pid) ) { return true; } + // } + // return false; + // } + + // check to see if this is a valid meson + bool isMeson( const int & pid ) + { + if( extraBits(pid) > 0 ) { return false; } + if( abspid(pid) <= 100 ) { return false; } + if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) { return false; } + int aid = abspid(pid); + if( aid == 130 || aid == 310 || aid == 210 ) { return true; } + // EvtGen uses some odd numbers + if( aid == 150 || aid == 350 || aid == 510 || aid == 530 ) { return true; } + // pomeron, etc. + if( pid == 110 || pid == 990 || pid == 9990 ) { return true; } + if( digit(nj,pid) > 0 && digit(nq3,pid) > 0 + && digit(nq2,pid) > 0 && digit(nq1,pid) == 0 ) { + // check for illegal antiparticles + if( digit(nq3,pid) == digit(nq2,pid) && pid < 0 ) { + return false; + } else { + return true; + } + } + return false; + } -// check to see if this is a valid diquark -bool isDiQuark( const int & pid ) -{ - if( extraBits(pid) > 0 ) { return false; } - if( abspid(pid) <= 100 ) { return false; } - if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) { return false; } - if( digit(nj,pid) > 0 && digit(nq3,pid) == 0 - && digit(nq2,pid) > 0 && digit(nq1,pid) > 0 ) { // diquark signature - // EvtGen uses the diquarks for quark pairs, so, for instance, - // 5501 is a valid "diquark" for EvtGen - //if( digit(nj) == 1 && digit(nq2) == digit(nq1) ) { // illegal - // return false; - //} else { - return true; - //} + // check to see if this is a valid baryon + bool isBaryon( const int & pid ) + { + if( extraBits(pid) > 0 ) { return false; } + if( abspid(pid) <= 100 ) { return false; } + if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) { return false; } + if( abspid(pid) == 2110 || abspid(pid) == 2210 ) { return true; } + if( digit(nj,pid) > 0 && digit(nq3,pid) > 0 + && digit(nq2,pid) > 0 && digit(nq1,pid) > 0 ) { return true; } + return false; } - return false; -} -// is this a valid hadron ID? -bool isHadron( const int & pid ) -{ - if( extraBits(pid) > 0 ) { return false; } - if( isMeson(pid) ) { return true; } - if( isBaryon(pid) ) { return true; } - if( isPentaquark(pid) ) { return true; } - return false; -} -// is this a valid lepton ID? -bool isLepton( const int & pid ) -{ - if( extraBits(pid) > 0 ) { return false; } - if( fundamentalID(pid) >= 11 && fundamentalID(pid) <= 18 ) { return true; } - return false; -} + // check to see if this is a valid diquark + bool isDiQuark( const int & pid ) + { + if( extraBits(pid) > 0 ) { return false; } + if( abspid(pid) <= 100 ) { return false; } + if( fundamentalID(pid) <= 100 && fundamentalID(pid) > 0 ) { return false; } + if( digit(nj,pid) > 0 && digit(nq3,pid) == 0 + && digit(nq2,pid) > 0 && digit(nq1,pid) > 0 ) { // diquark signature + // EvtGen uses the diquarks for quark pairs, so, for instance, + // 5501 is a valid "diquark" for EvtGen + //if( digit(nj) == 1 && digit(nq2) == digit(nq1) ) { // illegal + // return false; + //} else { + return true; + //} + } + return false; + } -// -// This implements the 2006 Monte Carlo nuclear code scheme. -// Ion numbers are +/- 10LZZZAAAI. -// AAA is A - total baryon number -// ZZZ is Z - total charge -// L is the total number of strange quarks. -// I is the isomer number, with I=0 corresponding to the ground state. -bool isNucleus( const int & pid ) -{ - // a proton can also be a Hydrogen nucleus - if( abspid(pid) == 2212 ) { return true; } - // new standard: +/- 10LZZZAAAI - if( ( digit(n10,pid) == 1 ) && ( digit(n9,pid) == 0 ) ) { - // charge should always be less than or equal to baryon number - // the following line is A >= Z - if( (abspid(pid)/10)%1000 >= (abspid(pid)/10000)%1000 ) { return true; } - } - return false; -} + // is this a valid hadron ID? + bool isHadron( const int & pid ) + { + if( extraBits(pid) > 0 ) { return false; } + if( isMeson(pid) ) { return true; } + if( isBaryon(pid) ) { return true; } + if( isPentaquark(pid) ) { return true; } + return false; + } + // is this a valid lepton ID? + bool isLepton( const int & pid ) + { + if( extraBits(pid) > 0 ) { return false; } + if( fundamentalID(pid) >= 11 && fundamentalID(pid) <= 18 ) { return true; } + return false; + } -// check to see if this is a valid pentaquark -bool isPentaquark( const int & pid ) -{ - // a pentaquark is of the form 9abcdej, - // where j is the spin and a, b, c, d, and e are quarks - if( extraBits(pid) > 0 ) { return false; } - if( digit(n,pid) != 9 ) { return false; } - if( digit(nr,pid) == 9 || digit(nr,pid) == 0 ) { return false; } - if( digit(nj,pid) == 9 || digit(nl,pid) == 0 ) { return false; } - if( digit(nq1,pid) == 0 ) { return false; } - if( digit(nq2,pid) == 0 ) { return false; } - if( digit(nq3,pid) == 0 ) { return false; } - if( digit(nj,pid) == 0 ) { return false; } - // check ordering - if( digit(nq2,pid) > digit(nq1,pid) ) { return false; } - if( digit(nq1,pid) > digit(nl,pid) ) { return false; } - if( digit(nl,pid) > digit(nr,pid) ) { return false; } - return true; -} + // + // This implements the 2006 Monte Carlo nuclear code scheme. + // Ion numbers are +/- 10LZZZAAAI. + // AAA is A - total baryon number + // ZZZ is Z - total charge + // L is the total number of strange quarks. + // I is the isomer number, with I=0 corresponding to the ground state. + bool isNucleus( const int & pid ) + { + // a proton can also be a Hydrogen nucleus + if( abspid(pid) == 2212 ) { return true; } + // new standard: +/- 10LZZZAAAI + if( ( digit(n10,pid) == 1 ) && ( digit(n9,pid) == 0 ) ) { + // charge should always be less than or equal to baryon number + // the following line is A >= Z + if( (abspid(pid)/10)%1000 >= (abspid(pid)/10000)%1000 ) { return true; } + } + return false; + } -// is this a SUSY? -bool isSUSY( const int & pid ) -{ - // fundamental SUSY particles have n = 1 or 2 - if( extraBits(pid) > 0 ) { return false; } - if( digit(n,pid) != 1 && digit(n,pid) != 2 ) { return false; } - if( digit(nr,pid) != 0 ) { return false; } - // check fundamental part - if( fundamentalID(pid) == 0 ) { return false; } - return true; -} + // check to see if this is a valid pentaquark + bool isPentaquark( const int & pid ) + { + // a pentaquark is of the form 9abcdej, + // where j is the spin and a, b, c, d, and e are quarks + if( extraBits(pid) > 0 ) { return false; } + if( digit(n,pid) != 9 ) { return false; } + if( digit(nr,pid) == 9 || digit(nr,pid) == 0 ) { return false; } + if( digit(nj,pid) == 9 || digit(nl,pid) == 0 ) { return false; } + if( digit(nq1,pid) == 0 ) { return false; } + if( digit(nq2,pid) == 0 ) { return false; } + if( digit(nq3,pid) == 0 ) { return false; } + if( digit(nj,pid) == 0 ) { return false; } + // check ordering + if( digit(nq2,pid) > digit(nq1,pid) ) { return false; } + if( digit(nq1,pid) > digit(nl,pid) ) { return false; } + if( digit(nl,pid) > digit(nr,pid) ) { return false; } + return true; + } -// is this an R-hadron? -bool isRhadron( const int & pid ) -{ - // an R-hadron is of the form 10abcdj, - // where j is the spin and a, b, c, and d are quarks or gluons - if( extraBits(pid) > 0 ) { return false; } - if( digit(n,pid) != 1 ) { return false; } - if( digit(nr,pid) != 0 ) { return false; } - // make sure this isn't a SUSY particle - if( isSUSY(pid) ) { return false; } - // All R-hadrons have at least 3 core digits - if( digit(nq2,pid) == 0 ) { return false; } - if( digit(nq3,pid) == 0 ) { return false; } - if( digit(nj,pid) == 0 ) { return false; } - return true; -} + // is this a SUSY? + bool isSUSY( const int & pid ) + { + // fundamental SUSY particles have n = 1 or 2 + if( extraBits(pid) > 0 ) { return false; } + if( digit(n,pid) != 1 && digit(n,pid) != 2 ) { return false; } + if( digit(nr,pid) != 0 ) { return false; } + // check fundamental part + if( fundamentalID(pid) == 0 ) { return false; } + return true; + } -// does this particle contain an up quark? -bool hasUp( const int & pid) -{ - if( extraBits(pid) > 0 ) { return false; } - if( fundamentalID(pid) > 0 ) { return false; } - if( digit(nq3,pid) == 2 || digit(nq2,pid) == 2 || digit(nq1,pid) == 2 ) { return true; } - return false; -} -// does this particle contain a down quark? -bool hasDown( const int & pid) -{ - if( extraBits(pid) > 0 ) { return false; } - if( fundamentalID(pid) > 0 ) { return false; } - if( digit(nq3,pid) == 1 || digit(nq2,pid) == 1 || digit(nq1,pid) == 1 ) { return true; } - return false; -} -// does this particle contain a strange quark? -bool hasStrange( const int & pid ) -{ - if( extraBits(pid) > 0 ) { return false; } - if( fundamentalID(pid) > 0 ) { return false; } - if( digit(nq3,pid) == 3 || digit(nq2,pid) == 3 || digit(nq1,pid) == 3 ) { return true; } - return false; -} -// does this particle contain a charm quark? -bool hasCharm( const int & pid ) -{ - if( extraBits(pid) > 0 ) { return false; } - if( fundamentalID(pid) > 0 ) { return false; } - if( digit(nq3,pid) == 4 || digit(nq2,pid) == 4 || digit(nq1,pid) == 4 ) { return true; } - return false; -} -// does this particle contain a bottom quark? -bool hasBottom( const int & pid ) -{ - if( extraBits(pid) > 0 ) { return false; } - if( fundamentalID(pid) > 0 ) { return false; } - if( digit(nq3,pid) == 5 || digit(nq2,pid) == 5 || digit(nq1,pid) == 5 ) { return true; } - return false; -} -// does this particle contain a top quark? -bool hasTop( const int & pid ) -{ - if( extraBits(pid) > 0 ) { return false; } - if( fundamentalID(pid) > 0 ) { return false; } - if( digit(nq3,pid) == 6 || digit(nq2,pid) == 6 || digit(nq1,pid) == 6 ) { return true; } - return false; -} + // is this an R-hadron? + bool isRhadron( const int & pid ) + { + // an R-hadron is of the form 10abcdj, + // where j is the spin and a, b, c, and d are quarks or gluons + if( extraBits(pid) > 0 ) { return false; } + if( digit(n,pid) != 1 ) { return false; } + if( digit(nr,pid) != 0 ) { return false; } + // make sure this isn't a SUSY particle + if( isSUSY(pid) ) { return false; } + // All R-hadrons have at least 3 core digits + if( digit(nq2,pid) == 0 ) { return false; } + if( digit(nq3,pid) == 0 ) { return false; } + if( digit(nj,pid) == 0 ) { return false; } + return true; + } -// --- other information -// -// jSpin returns 2J+1, where J is the total spin -int jSpin( const int & pid ) -{ - if( fundamentalID(pid) > 0 ) { - // some of these are known - int fund = fundamentalID(pid); - if( fund > 0 && fund < 7 ) return 2; - if( fund == 9 ) return 3; - if( fund > 10 && fund < 17 ) return 2; - if( fund > 20 && fund < 25 ) return 3; - return 0; - } else if( extraBits(pid) > 0 ) { - return 0; + // does this particle contain an up quark? + bool hasUp( const int & pid) + { + if( extraBits(pid) > 0 ) { return false; } + if( fundamentalID(pid) > 0 ) { return false; } + if( digit(nq3,pid) == 2 || digit(nq2,pid) == 2 || digit(nq1,pid) == 2 ) { return true; } + return false; } - return abspid(pid)%10; -} -// sSpin returns 2S+1, where S is the spin -int sSpin( const int & pid ) -{ - if( !isMeson(pid) ) { return 0; } - int inl = digit(nl,pid); - //int tent = digit(n,pid); - int js = digit(nj,pid); - if( digit(n,pid) == 9 ) { return 0; } // tentative ID - //if( tent == 9 ) { return 0; } // tentative assignment - if( inl == 0 && js >= 3 ) { - return 1; - } else if( inl == 0 && js == 1 ) { - return 0; - } else if( inl == 1 && js >= 3 ) { - return 0; - } else if( inl == 2 && js >= 3 ) { - return 1; - } else if( inl == 1 && js == 1 ) { - return 1; - } else if( inl == 3 && js >= 3 ) { - return 1; + // does this particle contain a down quark? + bool hasDown( const int & pid) + { + if( extraBits(pid) > 0 ) { return false; } + if( fundamentalID(pid) > 0 ) { return false; } + if( digit(nq3,pid) == 1 || digit(nq2,pid) == 1 || digit(nq1,pid) == 1 ) { return true; } + return false; } - // default to zero - return 0; -} -// lSpin returns 2L+1, where L is the orbital angular momentum -int lSpin( const int & pid ) -{ - if( !isMeson(pid) ) { return 0; } - int inl = digit(nl,pid); - //int tent = digit(n,pid); - int js = digit(nj,pid); - if( digit(n,pid) == 9 ) { return 0; } // tentative ID - if( inl == 0 && js == 3 ) { - return 0; - } else if( inl == 0 && js == 5 ) { - return 1; - } else if( inl == 0 && js == 7 ) { - return 2; - } else if( inl == 0 && js == 9 ) { - return 3; - } else if( inl == 0 && js == 1 ) { - return 0; - } else if( inl == 1 && js == 3 ) { - return 1; - } else if( inl == 1 && js == 5 ) { - return 2; - } else if( inl == 1 && js == 7 ) { - return 3; - } else if( inl == 1 && js == 9 ) { - return 4; - } else if( inl == 2 && js == 3 ) { - return 1; - } else if( inl == 2 && js == 5 ) { - return 2; - } else if( inl == 2 && js == 7 ) { - return 3; - } else if( inl == 2 && js == 9 ) { - return 4; - } else if( inl == 1 && js == 1 ) { - return 1; - } else if( inl == 3 && js == 3 ) { - return 2; - } else if( inl == 3 && js == 5 ) { - return 3; - } else if( inl == 3 && js == 7 ) { - return 4; - } else if( inl == 3 && js == 9 ) { - return 5; + // does this particle contain a strange quark? + bool hasStrange( const int & pid ) + { + if( extraBits(pid) > 0 ) { return false; } + if( fundamentalID(pid) > 0 ) { return false; } + if( digit(nq3,pid) == 3 || digit(nq2,pid) == 3 || digit(nq1,pid) == 3 ) { return true; } + return false; + } + // does this particle contain a charm quark? + bool hasCharm( const int & pid ) + { + if( extraBits(pid) > 0 ) { return false; } + if( fundamentalID(pid) > 0 ) { return false; } + if( digit(nq3,pid) == 4 || digit(nq2,pid) == 4 || digit(nq1,pid) == 4 ) { return true; } + return false; + } + // does this particle contain a bottom quark? + bool hasBottom( const int & pid ) + { + if( extraBits(pid) > 0 ) { return false; } + if( fundamentalID(pid) > 0 ) { return false; } + if( digit(nq3,pid) == 5 || digit(nq2,pid) == 5 || digit(nq1,pid) == 5 ) { return true; } + return false; + } + // does this particle contain a top quark? + bool hasTop( const int & pid ) + { + if( extraBits(pid) > 0 ) { return false; } + if( fundamentalID(pid) > 0 ) { return false; } + if( digit(nq3,pid) == 6 || digit(nq2,pid) == 6 || digit(nq1,pid) == 6 ) { return true; } + return false; } - // default to zero - return 0; -} -// 3 times the charge -int threeCharge( const int & pid ) -{ - int charge=0; - int ida, sid; - unsigned short q1, q2, q3; - static int ch100[100] = { -1, 2,-1, 2,-1, 2,-1, 2, 0, 0, - -3, 0,-3, 0,-3, 0,-3, 0, 0, 0, - 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, - 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, - 0, 6, 3, 6, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; - q1 = digit(nq1,pid); - q2 = digit(nq2,pid); - q3 = digit(nq3,pid); - ida = abspid(pid); - sid = fundamentalID(pid); - if( ida == 0 || extraBits(pid) > 0 ) { // ion or illegal - return 0; - } else if( sid > 0 && sid <= 100 ) { // use table - charge = ch100[sid-1]; - if(ida==1000017 || ida==1000018) { charge = 0; } - if(ida==1000034 || ida==1000052) { charge = 0; } - if(ida==1000053 || ida==1000054) { charge = 0; } - if(ida==5100061 || ida==5100062) { charge = 6; } - } else if( digit(nj,pid) == 0 ) { // KL, Ks, or undefined - return 0; - } else if( isMeson(pid) ) { // mesons - if( q2 == 3 || q2 == 5 ) { - charge = ch100[q3-1] - ch100[q2-1]; - } else { - charge = ch100[q2-1] - ch100[q3-1]; - } - } else if( isDiQuark(pid) ) { // diquarks - charge = ch100[q2-1] + ch100[q1-1]; - } else if( isBaryon(pid) ) { // baryons - charge = ch100[q3-1] + ch100[q2-1] + ch100[q1-1]; - } else { // unknown + // --- other information + // + // jSpin returns 2J+1, where J is the total spin + int jSpin( const int & pid ) + { + if( fundamentalID(pid) > 0 ) { + // some of these are known + int fund = fundamentalID(pid); + if( fund > 0 && fund < 7 ) return 2; + if( fund == 9 ) return 3; + if( fund > 10 && fund < 17 ) return 2; + if( fund > 20 && fund < 25 ) return 3; + return 0; + } else if( extraBits(pid) > 0 ) { + return 0; + } + return abspid(pid)%10; + } + // sSpin returns 2S+1, where S is the spin + int sSpin( const int & pid ) + { + if( !isMeson(pid) ) { return 0; } + int inl = digit(nl,pid); + //int tent = digit(n,pid); + int js = digit(nj,pid); + if( digit(n,pid) == 9 ) { return 0; } // tentative ID + //if( tent == 9 ) { return 0; } // tentative assignment + if( inl == 0 && js >= 3 ) { + return 1; + } else if( inl == 0 && js == 1 ) { + return 0; + } else if( inl == 1 && js >= 3 ) { + return 0; + } else if( inl == 2 && js >= 3 ) { + return 1; + } else if( inl == 1 && js == 1 ) { + return 1; + } else if( inl == 3 && js >= 3 ) { + return 1; + } + // default to zero return 0; } - if( charge == 0 ) { + // lSpin returns 2L+1, where L is the orbital angular momentum + int lSpin( const int & pid ) + { + if( !isMeson(pid) ) { return 0; } + int inl = digit(nl,pid); + //int tent = digit(n,pid); + int js = digit(nj,pid); + if( digit(n,pid) == 9 ) { return 0; } // tentative ID + if( inl == 0 && js == 3 ) { + return 0; + } else if( inl == 0 && js == 5 ) { + return 1; + } else if( inl == 0 && js == 7 ) { + return 2; + } else if( inl == 0 && js == 9 ) { + return 3; + } else if( inl == 0 && js == 1 ) { + return 0; + } else if( inl == 1 && js == 3 ) { + return 1; + } else if( inl == 1 && js == 5 ) { + return 2; + } else if( inl == 1 && js == 7 ) { + return 3; + } else if( inl == 1 && js == 9 ) { + return 4; + } else if( inl == 2 && js == 3 ) { + return 1; + } else if( inl == 2 && js == 5 ) { + return 2; + } else if( inl == 2 && js == 7 ) { + return 3; + } else if( inl == 2 && js == 9 ) { + return 4; + } else if( inl == 1 && js == 1 ) { + return 1; + } else if( inl == 3 && js == 3 ) { + return 2; + } else if( inl == 3 && js == 5 ) { + return 3; + } else if( inl == 3 && js == 7 ) { + return 4; + } else if( inl == 3 && js == 9 ) { + return 5; + } + // default to zero return 0; - } else if( pid < 0 ) { - charge = -charge; } - return charge; -} -}} + // 3 times the charge + int threeCharge( const int & pid ) + { + int charge=0; + int ida, sid; + unsigned short q1, q2, q3; + static int ch100[100] = { -1, 2,-1, 2,-1, 2,-1, 2, 0, 0, + -3, 0,-3, 0,-3, 0,-3, 0, 0, 0, + 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 3, 0, 0, 3, 0, 0, 0, + 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 6, 3, 6, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; + q1 = digit(nq1,pid); + q2 = digit(nq2,pid); + q3 = digit(nq3,pid); + ida = abspid(pid); + sid = fundamentalID(pid); + if( ida == 0 || extraBits(pid) > 0 ) { // ion or illegal + return 0; + } else if( sid > 0 && sid <= 100 ) { // use table + charge = ch100[sid-1]; + if(ida==1000017 || ida==1000018) { charge = 0; } + if(ida==1000034 || ida==1000052) { charge = 0; } + if(ida==1000053 || ida==1000054) { charge = 0; } + if(ida==5100061 || ida==5100062) { charge = 6; } + } else if( digit(nj,pid) == 0 ) { // KL, Ks, or undefined + return 0; + } else if( isMeson(pid) ) { // mesons + if( q2 == 3 || q2 == 5 ) { + charge = ch100[q3-1] - ch100[q2-1]; + } else { + charge = ch100[q2-1] - ch100[q3-1]; + } + } else if( isDiQuark(pid) ) { // diquarks + charge = ch100[q2-1] + ch100[q1-1]; + } else if( isBaryon(pid) ) { // baryons + charge = ch100[q3-1] + ch100[q2-1] + ch100[q1-1]; + } else { // unknown + return 0; + } + if( charge == 0 ) { + return 0; + } else if( pid < 0 ) { + charge = -charge; + } + return charge; + } + + + } +}
More information about the Rivet-svn mailing list |