00001 // -*- C++ -*-
00002 #include "Rivet/Projections/InitialQuarks.hh"
00003 #include "Rivet/Tools/Logging.hh"
00004 #include "Rivet/Cmp.hh"
00005 
00006 
00007 #define IS_PARTON_PDGID(id) ( abs(id) <= 100 && abs(id) != 22 && (abs(id) < 11 || abs(id) > 18) )
00008 
00009 
00010 namespace Rivet {
00011 
00012 
00013   int InitialQuarks::compare(const Projection& p) const {
00014     return EQUIVALENT;
00015   }
00016 
00017 
00018   void InitialQuarks::project(const Event& e) {
00019     _theParticles.clear();
00020 
00021     foreach (const GenParticle* p, Rivet::particles(e.genEvent())) {
00022       const GenVertex* pv = p->production_vertex();
00023       const GenVertex* dv = p->end_vertex();
00024       const PdgId pid = abs(p->pdg_id());
00025       bool passed = (inRange(pid, 1, 5));
00026       if (passed) {
00027         if (pv != 0) {
00028           foreach (const GenParticle* pp, particles_in(pv)) {
00029             // Only accept if parent is electron or Z0
00030             const PdgId pid = abs(pp->pdg_id());
00031             passed = (pid == ELECTRON || abs(pp->pdg_id()) == ZBOSON);
00032           }
00033         } else {
00034           passed = false;
00035         }
00036       }
00037 
00038       if (getLog().isActive(Log::TRACE)) {
00039         const int st = p->status();
00040         const double pT = p->momentum().perp();
00041         const double eta = p->momentum().eta();
00042         getLog() << Log::TRACE << std::boolalpha
00043                  << "ID = " << p->pdg_id() << ", status = " << st << ", pT = " << pT
00044                  << ", eta = " << eta << ": result = " << passed << endl;
00045         if (pv != 0) {
00046           foreach (const GenParticle* pp, particles_in(pv)) {
00047             getLog() << Log::TRACE << std::boolalpha
00048                      << "     parent ID = " << pp->pdg_id() << endl;
00049           }
00050         }
00051         if (dv != 0) {
00052           foreach (const GenParticle* pp, particles_out(dv)) {
00053             getLog() << Log::TRACE << std::boolalpha
00054                      << "     child ID  = " << pp->pdg_id() << endl;
00055           }
00056         }
00057       }
00058       if (passed) _theParticles.push_back(Particle(*p));
00059     }
00060     getLog() << Log::DEBUG << "Number of initial quarks = "
00061              << _theParticles.size() << endl;
00062     if (! _theParticles.empty())
00063       for (size_t i=0 ; i < _theParticles.size() ; i++)
00064         getLog() << Log::DEBUG << "Initial quark[" << i << "] = "
00065                  << _theParticles[i].pdgId() << std::endl;
00066   }
00067 
00068 
00069 }
InitialQuarks.cc
