rivet is hosted by Hepforge, IPPP Durham
ALEPH_2001_S4656318.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 #include "Rivet/Analysis.hh"
00003 #include "Rivet/RivetYODA.hh"
00004 #include "Rivet/Tools/ParticleIdUtils.hh"
00005 #include "Rivet/Projections/Beam.hh"
00006 #include "Rivet/Projections/FinalState.hh"
00007 #include "Rivet/Projections/ChargedFinalState.hh"
00008 
00009 
00010 /// @todo Use inline PID functions instead
00011 #define IS_PARTON_PDGID(id) ( abs(id) <= 100 && abs(id) != 22 && (abs(id) < 11 || abs(id) > 18) )
00012 #define IS_BHADRON_PDGID(id) ( ((abs(id)/100)%10 == 5) || (abs(id) >= 5000 && abs(id) <= 5999) )
00013 
00014 namespace Rivet {
00015 
00016 
00017   /// @brief DELPHI b-fragmentation measurement
00018   /// @author Hendrik Hoeth
00019   class ALEPH_2001_S4656318 : public Analysis {
00020   public:
00021 
00022     /// Constructor
00023    ALEPH_2001_S4656318()
00024       : Analysis("ALEPH_2001_S4656318")
00025     {
00026     }
00027 
00028 
00029     /// @name Analysis methods
00030     //@{
00031 
00032     /// Book projections and histograms
00033     void init() {
00034       addProjection(Beam(), "Beams");
00035       addProjection(ChargedFinalState(), "FS");
00036       _histXbweak     = bookHisto1D(8, 1, 1);
00037       _histXbprim     = bookHisto1D(8, 1, 2);
00038       _histMeanXbweak = bookProfile1D(6, 1, 1);
00039       _histMeanXbprim = bookProfile1D(7, 1, 1);
00040     }
00041 
00042 
00043     void analyze(const Event& e) {
00044       const FinalState& fs = applyProjection<FinalState>(e, "FS");
00045       const size_t numParticles = fs.particles().size();
00046 
00047       // Even if we only generate hadronic events, we still need a cut on numCharged >= 2.
00048       if (numParticles < 2) {
00049         MSG_DEBUG("Failed ncharged cut");
00050         vetoEvent;
00051       }
00052       MSG_DEBUG("Passed ncharged cut");
00053 
00054       // Get event weight for histo filling
00055       const double weight = e.weight();
00056 
00057       // Get beams and average beam momentum
00058       const ParticlePair& beams = applyProjection<Beam>(e, "Beams").beams();
00059       const double meanBeamMom = ( beams.first.momentum().vector3().mod() +
00060                                    beams.second.momentum().vector3().mod() ) / 2.0;
00061       MSG_DEBUG("Avg beam momentum = " << meanBeamMom);
00062 
00063 
00064       foreach (const GenParticle* p, particles(e.genEvent())) {
00065         const GenVertex* pv = p->production_vertex();
00066         const GenVertex* dv = p->end_vertex();
00067         if (IS_BHADRON_PDGID(p->pdg_id())) {
00068           const double xp = p->momentum().e()/meanBeamMom;
00069 
00070           // If the B-hadron has a parton as parent, call it primary B-hadron:
00071           if (pv) {
00072             bool is_primary = false;
00073             for (GenVertex::particles_in_const_iterator pp = pv->particles_in_const_begin(); pp != pv->particles_in_const_end() ; ++pp) {
00074               if (IS_PARTON_PDGID((*pp)->pdg_id())) is_primary = true;
00075             }
00076             if (is_primary) {
00077               _histXbprim->fill(xp, weight);
00078               _histMeanXbprim->fill(_histMeanXbprim->bin(0).xMean(), xp, weight);
00079             }
00080           }
00081 
00082           // If the B-hadron has no B-hadron as a child, it decayed weakly:
00083           if (dv) {
00084             bool is_weak = true;
00085             for (GenVertex::particles_out_const_iterator pp = dv->particles_out_const_begin() ;
00086                  pp != dv->particles_out_const_end() ; ++pp) {
00087               if (IS_BHADRON_PDGID((*pp)->pdg_id())) {
00088                 is_weak = false;
00089               }
00090             }
00091             if (is_weak) {
00092               _histXbweak->fill(xp, weight);
00093               _histMeanXbweak->fill(_histMeanXbweak->bin(0).xMean(), xp, weight);
00094             }
00095           }
00096 
00097         }
00098       }
00099     }
00100 
00101 
00102     // Finalize
00103     void finalize() {
00104       normalize(_histXbprim);
00105       normalize(_histXbweak);
00106     }
00107 
00108 
00109   private:
00110 
00111     /// Store the weighted sums of numbers of charged / charged+neutral
00112     /// particles - used to calculate average number of particles for the
00113     /// inclusive single particle distributions' normalisations.
00114 
00115     Histo1DPtr _histXbprim;
00116     Histo1DPtr _histXbweak;
00117 
00118     Profile1DPtr _histMeanXbprim;
00119     Profile1DPtr _histMeanXbweak;
00120 
00121     //@}
00122 
00123   };
00124 
00125 
00126 
00127   // The hook for the plugin system
00128   DECLARE_RIVET_PLUGIN(ALEPH_2001_S4656318);
00129 
00130 }