rivet is hosted by Hepforge, IPPP Durham
FinalState.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 #include "Rivet/Projections/FinalState.hh"
00003 
00004 namespace Rivet {
00005 
00006   // @deprecated, keep for backwards compatibility for now.
00007   FinalState::FinalState(double mineta, double maxeta, double minpt)
00008   {
00009     setName("FinalState");
00010     const bool openpt = isZero(minpt);
00011     const bool openeta = (mineta <= -MAXRAPIDITY && maxeta >= MAXRAPIDITY);
00012     MSG_TRACE("Check for open FS conditions:" << std::boolalpha
00013               << " eta=" << openeta
00014               << ", pt=" << openpt);
00015     if ( openpt && openeta ) {
00016       _cuts = Cuts::open();
00017     }
00018     else {
00019       addProjection(FinalState(), "OpenFS");
00020       if ( openeta ) 
00021     _cuts = Cuts::pT >= minpt;
00022       else if ( openpt )
00023     _cuts = EtaIn(mineta,maxeta);
00024       else
00025     _cuts = EtaIn(mineta,maxeta) & (Cuts::pT >= minpt);
00026     }
00027   }
00028 
00029   FinalState::FinalState(Cut c)
00030     : _cuts(c)
00031   {
00032     setName("FinalState");
00033     const bool open = ( c == Cuts::open() );
00034     MSG_TRACE("Check for open FS conditions: " << std::boolalpha << open);
00035     if ( !open ) {
00036       addProjection(FinalState(), "OpenFS");
00037     }
00038   }
00039 
00040 
00041   // FinalState::FinalState(const vector<pair<double, double> >& etaRanges, double minpt)
00042   //   : _etaRanges(etaRanges), _ptmin(minpt)
00043   // {
00044   //   setName("FinalState");
00045   //   const bool openpt = isZero(minpt);
00046   //   /// @todo Properly check whether any of these eta ranges (or their combination) are actually open
00047   //   const bool openeta = etaRanges.empty();
00048   //   MSG_TRACE("Check for open FS conditions:" << std::boolalpha
00049   //             << " eta=" << openeta
00050   //             << ", pt=" << openpt);
00051   //   if (!openeta || !openpt) {
00052   //     addProjection(FinalState(), "OpenFS");
00053   //   }
00054   // }
00055 
00056 
00057   /// @todo HOW DO WE COMPARE CUTS OBJECTS?
00058   int FinalState::compare(const Projection& p) const {
00059     const FinalState& other = dynamic_cast<const FinalState&>(p);
00060 
00061     //MSG_TRACE("FS::compare: " << 1 << " " << this << " " << &p);
00062     //    std::vector<std::pair<double, double> > eta1(_etaRanges);
00063     //std::vector<std::pair<double, double> > eta2(other._etaRanges);
00064     //std::sort(eta1.begin(), eta1.end());
00065     //std::sort(eta2.begin(), eta2.end());
00066 
00067     //MSG_TRACE("FS::compare: " << 2 << " " << this << " " << &p);
00068     //if (eta1 < eta2) return ORDERED;
00069     //else if (eta2 < eta1) return UNORDERED;
00070 
00071     //MSG_TRACE("FS::compare: " << 3 << " " << this << " " << &p);
00072     //return cmp(_ptmin, other._ptmin);
00073     return _cuts == other._cuts ? EQUIVALENT : UNDEFINED;
00074   }
00075 
00076 
00077 
00078   void FinalState::project(const Event& e) {
00079     _theParticles.clear();
00080 
00081     // Handle "open FS" special case
00082     if ( _cuts == Cuts::open() ) {
00083       //MSG_TRACE("Open FS processing: should only see this once per event ("
00084       //           << e.genEvent().event_number() << ")");
00085       foreach (const GenParticle* p, Rivet::particles(e.genEvent())) {
00086         if (p->status() == 1) {
00087           //MSG_TRACE("FS GV = " << p->production_vertex());
00088           _theParticles.push_back(Particle(*p));
00089         }
00090       }
00091       return;
00092     }
00093 
00094     // If this is not itself the "open" FS, base the calculations on the open FS' results
00095     /// @todo In general, we'd like to calculate a restrictive FS based on the most restricted superset FS.
00096     const Particles allstable = applyProjection<FinalState>(e, "OpenFS").particles();
00097     foreach (const Particle& p, allstable) {
00098       const bool passed = accept(p);
00099       MSG_TRACE("Choosing: ID = " << p.pdgId()
00100                 << ", pT = " << p.pT()
00101                 << ", eta = " << p.eta()
00102                 << ": result = " << std::boolalpha << passed);
00103       if (passed) _theParticles.push_back(p);
00104     }
00105     //MSG_DEBUG("Number of final-state particles = " << _theParticles.size());
00106   }
00107 
00108 
00109   /// Decide if a particle is to be accepted or not.
00110   bool FinalState::accept(const Particle& p) const {
00111     // Not having s.c. == 1 should never happen!
00112     assert(p.genParticle() == NULL || p.genParticle()->status() == 1);
00113 
00114     return _cuts->accept(p);
00115   }
00116 
00117 
00118 }