rivet is hosted by Hepforge, IPPP Durham
InvMassFinalState.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 #include "Rivet/Projections/InvMassFinalState.hh"
00003 #include "Rivet/Tools/Logging.hh"
00004 #include "Rivet/Cmp.hh"
00005 
00006 namespace Rivet {
00007 
00008 
00009   InvMassFinalState::InvMassFinalState(const FinalState& fsp,
00010                                        const pair<PdgId, PdgId>& idpair, // pair of decay products
00011                                        double minmass, // min inv mass
00012                                        double maxmass, // max inv mass
00013                                        double masstarget)
00014     : _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
00015   {
00016     setName("InvMassFinalState");
00017     addProjection(fsp, "FS");
00018     _decayids.push_back(idpair);
00019   }
00020 
00021 
00022   InvMassFinalState::InvMassFinalState(const FinalState& fsp,
00023                                        const vector<pair<PdgId, PdgId> >& idpairs,  // vector of pairs of decay products
00024                                        double minmass, // min inv mass
00025                                        double maxmass, // max inv mass
00026                                        double masstarget)
00027     : _decayids(idpairs), _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
00028   {
00029     setName("InvMassFinalState");
00030     addProjection(fsp, "FS");
00031   }
00032 
00033 
00034   InvMassFinalState::InvMassFinalState(const pair<PdgId, PdgId>& idpair, // pair of decay products
00035                                        double minmass, // min inv mass
00036                                        double maxmass, // max inv mass
00037                                        double masstarget)
00038     : _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
00039   {
00040     setName("InvMassFinalState");
00041     _decayids.push_back(idpair);
00042   }
00043 
00044 
00045   InvMassFinalState::InvMassFinalState(const vector<pair<PdgId, PdgId> >& idpairs,  // vector of pairs of decay products
00046                                        double minmass, // min inv mass
00047                                        double maxmass, // max inv mass
00048                                        double masstarget)
00049     : _decayids(idpairs), _minmass(minmass), _maxmass(maxmass), _masstarget(masstarget), _useTransverseMass(false)
00050   {
00051     setName("InvMassFinalState");
00052   }
00053 
00054 
00055   int InvMassFinalState::compare(const Projection& p) const {
00056     // First compare the final states we are running on
00057     int fscmp = mkNamedPCmp(p, "FS");
00058     if (fscmp != EQUIVALENT) return fscmp;
00059 
00060     // Then compare the two as final states
00061     const InvMassFinalState& other = dynamic_cast<const InvMassFinalState&>(p);
00062     fscmp = FinalState::compare(other);
00063     if (fscmp != EQUIVALENT) return fscmp;
00064 
00065     // Compare the mass limits
00066     int masstypecmp = cmp(_useTransverseMass, other._useTransverseMass);
00067     if (masstypecmp != EQUIVALENT) return masstypecmp;
00068     int massllimcmp = cmp(_minmass, other._minmass);
00069     if (massllimcmp != EQUIVALENT) return massllimcmp;
00070     int masshlimcmp = cmp(_maxmass, other._maxmass);
00071     if (masshlimcmp != EQUIVALENT) return masshlimcmp;
00072 
00073     // Compare the decay species
00074     int decaycmp = cmp(_decayids, other._decayids);
00075     if (decaycmp != EQUIVALENT) return decaycmp;
00076 
00077     // Finally compare them as final states
00078     return FinalState::compare(other);
00079   }
00080 
00081 
00082   void InvMassFinalState::project(const Event& e) {
00083     const FinalState& fs = applyProjection<FinalState>(e, "FS");
00084     calc(fs.particles());
00085   }
00086 
00087 
00088   void InvMassFinalState::calc(const ParticleVector& inparticles) {
00089     _theParticles.clear();
00090     _particlePairs.clear();
00091 
00092     // Containers for the particles of type specified in the pair
00093     vector<const Particle*> type1;
00094     vector<const Particle*> type2;
00095     // Get all the particles of the type specified in the pair from the particle list
00096     foreach (const Particle& ipart, inparticles) {
00097       // Loop around possible particle pairs (typedef needed to keep BOOST_FOREACH happy)
00098       foreach (const PdgIdPair& ipair, _decayids) {
00099         if (ipart.pdgId() == ipair.first) {
00100           if (accept(ipart)) {
00101             type1 += &ipart;
00102           }
00103         } else if (ipart.pdgId() == ipair.second) {
00104           if (accept(ipart)) {
00105             type2 += &ipart;
00106           }
00107         }
00108       }
00109     }
00110     if (type1.empty() || type2.empty()) return;
00111 
00112     // Temporary container of selected particles iterators
00113     // Useful to compare iterators and avoid double occurrences of the same
00114     // particle in case it matches with more than another particle
00115     vector<const Particle*> tmp;
00116 
00117     // Now calculate the inv mass
00118     pair<double, pair<Particle, Particle> > closestPair;
00119     closestPair.first = 1e30;
00120     foreach (const Particle* i1, type1) {
00121       foreach (const Particle* i2, type2) {
00122         // Check this is actually a pair
00123         // (if more than one pair in vector particles can be unrelated)
00124         bool found = false;
00125         foreach (const PdgIdPair& ipair, _decayids) {
00126           if (i1->pdgId() == ipair.first && i2->pdgId() == ipair.second) {
00127             found = true;
00128             break;
00129           }
00130         }
00131         if (!found) continue;
00132 
00133         FourMomentum v4 = i1->momentum() + i2->momentum();
00134         if (v4.mass2() < 0) {
00135           MSG_DEBUG("Constructed negative inv mass2: skipping!");
00136           continue;
00137         }
00138         bool passedMassCut = false;
00139         if (_useTransverseMass) {
00140           passedMassCut = inRange(massT(i1->momentum(), i2->momentum()), _minmass, _maxmass);
00141         } else {
00142           passedMassCut = inRange(v4.mass(), _minmass, _maxmass);
00143         }
00144 
00145         if (passedMassCut) {
00146           MSG_DEBUG("Selecting particles with IDs " << i1->pdgId() << " & " << i2->pdgId()
00147                     << " and mass = " << v4.mass()/GeV << " GeV");
00148           // Store accepted particles, avoiding duplicates
00149           if (find(tmp.begin(), tmp.end(), i1) == tmp.end()) {
00150             tmp.push_back(i1);
00151             _theParticles += *i1;
00152           }
00153           if (find(tmp.begin(), tmp.end(), i2) == tmp.end()) {
00154             tmp.push_back(i2);
00155             _theParticles += *i2;
00156           }
00157           // Store accepted particle pairs
00158           _particlePairs += make_pair(*i1, *i2);
00159           if (_masstarget>0.0) {
00160             double diff=fabs(v4.mass()-_masstarget);
00161             if (diff<closestPair.first) {
00162               closestPair.first=diff;
00163               closestPair.second=make_pair(*i1, *i2);
00164             }
00165           }
00166         }
00167       }
00168     }
00169     if (_masstarget > 0.0 && closestPair.first < 1e30) {
00170       _theParticles.clear();
00171       _particlePairs.clear();
00172       _theParticles += closestPair.second.first;
00173       _theParticles += closestPair.second.second;
00174       _particlePairs += closestPair.second;
00175     }
00176 
00177     MSG_DEBUG("Selected " << _theParticles.size() << " particles " << "(" << _particlePairs.size() << " pairs)");
00178     if (getLog().isActive(Log::TRACE)) {
00179       foreach (const Particle& p, _theParticles) {
00180         MSG_TRACE("ID: " << p.pdgId() << ", barcode: " << p.genParticle().barcode());
00181       }
00182     }
00183   }
00184 
00185 
00186   /// Constituent pairs
00187   const std::vector<std::pair<Particle, Particle> >& InvMassFinalState::particlePairs() const {
00188     return _particlePairs;
00189   }
00190 
00191 
00192 }