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 } Generated on Fri Dec 21 2012 15:03:41 for The Rivet MC analysis system by ![]() |