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     /// @todo This is all fragile and application-specific: remove!
00022 
00023     foreach (const GenParticle* p, Rivet::particles(e.genEvent())) {
00024       const GenVertex* pv = p->production_vertex();
00025       const GenVertex* dv = p->end_vertex();
00026       const PdgId pid = abs(p->pdg_id());
00027       bool passed = (inRange(pid, 1, 5));
00028       if (passed) {
00029         if (pv != 0) {
00030           foreach (const GenParticle* pp, particles_in(pv)) {
00031             // Only accept if parent is electron or Z0
00032             const PdgId pid = abs(pp->pdg_id());
00033             passed = (pid == ELECTRON || abs(pp->pdg_id()) == ZBOSON);
00034           }
00035         } else {
00036           passed = false;
00037         }
00038       }
00039 
00040       if (getLog().isActive(Log::TRACE)) {
00041         const int st = p->status();
00042         const double pT = p->momentum().perp();
00043         const double eta = p->momentum().eta();
00044         getLog() << Log::TRACE << std::boolalpha
00045                  << "ID = " << p->pdg_id() << ", status = " << st << ", pT = " << pT
00046                  << ", eta = " << eta << ": result = " << passed << endl;
00047         if (pv != 0) {
00048           foreach (const GenParticle* pp, particles_in(pv)) {
00049             getLog() << Log::TRACE << std::boolalpha
00050                      << "     parent ID = " << pp->pdg_id() << endl;
00051           }
00052         }
00053         if (dv != 0) {
00054           foreach (const GenParticle* pp, particles_out(dv)) {
00055             getLog() << Log::TRACE << std::boolalpha
00056                      << "     child ID  = " << pp->pdg_id() << endl;
00057           }
00058         }
00059       }
00060       if (passed) _theParticles.push_back(Particle(*p));
00061     }
00062     getLog() << Log::DEBUG << "Number of initial quarks = "
00063              << _theParticles.size() << endl;
00064     if (! _theParticles.empty())
00065       for (size_t i=0 ; i < _theParticles.size() ; i++)
00066         getLog() << Log::DEBUG << "Initial quark[" << i << "] = "
00067                  << _theParticles[i].pdgId() << std::endl;
00068   }
00069 
00070 
00071 }
InitialQuarks.cc
