| 20 |
|
|
| 21 |
|
#include "DataFormats/Math/interface/deltaR.h" |
| 22 |
|
#include "DataFormats/Common/interface/View.h" |
| 23 |
+ |
#include "DataFormats/MuonReco/interface/MuonFwd.h" |
| 24 |
|
#include "DataFormats/MuonReco/interface/Muon.h" |
| 25 |
|
#include "DataFormats/TrackReco/interface/Track.h" |
| 26 |
+ |
#include "DataFormats/VertexReco/interface/Vertex.h" |
| 27 |
+ |
|
| 28 |
+ |
#include "DataFormats/PatCandidates/interface/GenericParticle.h" |
| 29 |
|
|
| 30 |
|
#include "CommonTools/Utils/interface/StringCutObjectSelector.h" |
| 31 |
|
|
| 39 |
|
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" |
| 40 |
|
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h" |
| 41 |
|
#include <TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h> |
| 42 |
+ |
#include <TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h> |
| 43 |
+ |
#include "TrackingTools/IPTools/interface/IPTools.h" |
| 44 |
+ |
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" |
| 45 |
+ |
#include "TrackingTools/TransientTrack/interface/TransientTrack.h" |
| 46 |
+ |
#include "TrackingTools/Records/interface/TransientTrackRecord.h" |
| 47 |
|
|
| 48 |
|
#include "Math/GenVector/Boost.h" |
| 49 |
|
|
| 63 |
|
edm::InputTag muons_; |
| 64 |
|
/// input tracker tracks |
| 65 |
|
edm::InputTag tracks_; |
| 66 |
+ |
/// input primary vertex |
| 67 |
+ |
edm::InputTag vertices_; |
| 68 |
|
/// preselection for muons |
| 69 |
|
StringCutObjectSelector<reco::Muon> muonCut_; |
| 70 |
|
/// preselection for candidate "inner" tracks (hadrons that decay) |
| 78 |
|
/// maximum distance between the two tracks |
| 79 |
|
double maxTwoTrackDist_; |
| 80 |
|
|
| 81 |
+ |
/// Use only variables available in AOD |
| 82 |
+ |
bool aodDataOnly_; |
| 83 |
+ |
|
| 84 |
|
/// Magnetic field for extrapolation |
| 85 |
|
edm::ESHandle<MagneticField> magfield_; |
| 86 |
|
|
| 90 |
|
/// Global tracking geometry |
| 91 |
|
edm::ESHandle<GlobalTrackingGeometry> geometry_; |
| 92 |
|
|
| 93 |
+ |
/// Transient track builder, for better d0 determination |
| 94 |
+ |
edm::ESHandle<TransientTrackBuilder> ttBuilder_; |
| 95 |
+ |
|
| 96 |
|
struct EqByTrack { |
| 97 |
|
EqByTrack() : tk() {} |
| 98 |
|
EqByTrack(const reco::TrackRef & ref) : tk(ref) {} |
| 121 |
|
innerTrack(inner), outerTrack(outer), distance(-1) {} |
| 122 |
|
reco::TrackRef innerTrack; // never null |
| 123 |
|
reco::TrackRef outerTrack; // never null |
| 124 |
< |
edm::RefToBase<reco::Muon> innerMu; // may be null |
| 125 |
< |
edm::RefToBase<reco::Muon> outerMu; // may be null |
| 124 |
> |
edm::Ptr<reco::Muon> innerMu; // may be null |
| 125 |
> |
edm::Ptr<reco::Muon> outerMu; // may be null |
| 126 |
|
double distance; |
| 127 |
|
GlobalPoint decayPoint; |
| 128 |
|
reco::Particle::LorentzVector p4mu, p4nu, p4had; |
| 139 |
|
InFlightDecayRecoFinder::InFlightDecayRecoFinder(const edm::ParameterSet & iConfig) : |
| 140 |
|
muons_(iConfig.getParameter<edm::InputTag>("muons")), |
| 141 |
|
tracks_(iConfig.getParameter<edm::InputTag>("tracks")), |
| 142 |
+ |
vertices_(iConfig.getParameter<edm::InputTag>("vertices")), |
| 143 |
|
muonCut_(iConfig.getParameter<std::string>("muonCut")), |
| 144 |
|
innerTrackCut_(iConfig.getParameter<std::string>("innerTrackCut")), |
| 145 |
|
outerTrackCut_(iConfig.getParameter<std::string>("outerTrackCut")), |
| 146 |
|
maxDeltaR2_(std::pow(iConfig.getParameter<double>("maxDeltaR"),2)), |
| 147 |
|
maxDeltaPtRel_(iConfig.getParameter<double>("maxDeltaPtRel")), |
| 148 |
< |
maxTwoTrackDist_(iConfig.getParameter<double>("maxTwoTrackDist")) |
| 148 |
> |
maxTwoTrackDist_(iConfig.getParameter<double>("maxTwoTrackDist")), |
| 149 |
> |
aodDataOnly_(iConfig.getParameter<bool>("aodDataOnly")) |
| 150 |
|
{ |
| 151 |
|
produces<reco::TrackCollection>("inner"); |
| 152 |
|
produces<reco::TrackCollection>("outer"); |
| 153 |
+ |
produces<std::vector<reco::Muon> >("inner"); |
| 154 |
+ |
produces<std::vector<reco::Muon> >("outer"); |
| 155 |
+ |
produces<pat::GenericParticleCollection>(); |
| 156 |
|
} |
| 157 |
|
|
| 158 |
|
void |
| 163 |
|
// Get EventSetup data |
| 164 |
|
iSetup.get<IdealMagneticFieldRecord>().get(magfield_); |
| 165 |
|
//iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator_); |
| 166 |
< |
iSetup.get<GlobalTrackingGeometryRecord>().get(geometry_); |
| 166 |
> |
if (!aodDataOnly_) iSetup.get<GlobalTrackingGeometryRecord>().get(geometry_); |
| 167 |
> |
iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", ttBuilder_); |
| 168 |
> |
|
| 169 |
|
|
| 170 |
|
// Get Event data |
| 171 |
|
Handle<View<reco::Muon> > muons; |
| 148 |
– |
Handle<reco::TrackCollection> tracks; |
| 172 |
|
iEvent.getByLabel(muons_, muons); |
| 173 |
+ |
Handle<reco::TrackCollection> tracks; |
| 174 |
|
iEvent.getByLabel(tracks_, tracks); |
| 175 |
+ |
Handle<std::vector<reco::Vertex> > vertices; |
| 176 |
+ |
iEvent.getByLabel(vertices_, vertices); |
| 177 |
|
|
| 178 |
< |
auto_ptr<reco::TrackCollection> innerSels(new reco::TrackCollection()); |
| 179 |
< |
auto_ptr<reco::TrackCollection> outerSels(new reco::TrackCollection()); |
| 178 |
> |
auto_ptr<reco::TrackCollection> innerSelTks(new reco::TrackCollection()); |
| 179 |
> |
auto_ptr<reco::TrackCollection> outerSelTks(new reco::TrackCollection()); |
| 180 |
> |
auto_ptr<reco::MuonCollection> innerSelMus(new reco::MuonCollection()); |
| 181 |
> |
auto_ptr<reco::MuonCollection> outerSelMus(new reco::MuonCollection()); |
| 182 |
> |
auto_ptr<pat::GenericParticleCollection> output(new pat::GenericParticleCollection()); |
| 183 |
|
|
| 184 |
|
typedef std::vector<reco::TrackRef> TkRefs; |
| 185 |
|
EqByTrack eq; |
| 199 |
|
DecayCandidate cand(*it1,*it2); |
| 200 |
|
if (std::find(decayCandidates.begin(), decayCandidates.end(), cand) != decayCandidates.end()) continue; // can be seeded by both muons |
| 201 |
|
|
| 202 |
< |
edm::View<reco::Muon>::const_iterator innerMuIt = std::find_if(muBegin, muEnd, EqByTrack(*it1)); |
| 202 |
> |
edm::View<reco::Muon>::const_iterator innerMuIt = std::find_if(muBegin, muEnd, EqByTrack(*it1)); |
| 203 |
|
edm::View<reco::Muon>::const_iterator outerMuIt = std::find_if(muBegin, muEnd, EqByTrack(*it2)); |
| 204 |
< |
if (innerMuIt != muEnd) { cand.innerMu = muons->refAt(innerMuIt - muBegin); } |
| 205 |
< |
if (outerMuIt != muEnd) { cand.outerMu = muons->refAt(outerMuIt - muBegin); } |
| 204 |
> |
if (innerMuIt != muEnd) { cand.innerMu = muons->ptrAt(innerMuIt - muBegin); } |
| 205 |
> |
if (outerMuIt != muEnd) { cand.outerMu = muons->ptrAt(outerMuIt - muBegin); } |
| 206 |
|
|
| 207 |
|
TwoTrackMinimumDistance ttmd; |
| 208 |
|
TrajectoryStateTransform transformer; |
| 180 |
– |
TrajectoryStateOnSurface innerInnerTSOS = transformer.innerStateOnSurface(**it1, *geometry_, &*magfield_); // start from the two far states |
| 181 |
– |
TrajectoryStateOnSurface outerOuterTSOS = transformer.outerStateOnSurface(**it2, *geometry_, &*magfield_); // so that distance is less biased by shared hits |
| 209 |
|
|
| 210 |
< |
if (!ttmd.calculate(innerInnerTSOS, outerOuterTSOS)) { std::cerr << " TTMD failed :-( " << std::endl; continue; } |
| 210 |
> |
if (!aodDataOnly_) { |
| 211 |
> |
TrajectoryStateOnSurface innerInnerTSOS = transformer.innerStateOnSurface(**it1, *geometry_, &*magfield_); // start from the two far states |
| 212 |
> |
TrajectoryStateOnSurface outerOuterTSOS = transformer.outerStateOnSurface(**it2, *geometry_, &*magfield_); // so that distance is less biased by shared hits |
| 213 |
> |
if (!ttmd.calculate(innerInnerTSOS, outerOuterTSOS)) { std::cerr << " TTMD failed :-( " << std::endl; continue; } |
| 214 |
> |
} else { |
| 215 |
> |
FreeTrajectoryState innerFTS = transformer.initialFreeState(**it1, &*magfield_); // start from the two far states |
| 216 |
> |
FreeTrajectoryState outerFTS = transformer.initialFreeState(**it2, &*magfield_); // so that distance is less biased by shared hits |
| 217 |
> |
if (!ttmd.calculate(innerFTS, outerFTS)) { std::cerr << " TTMD failed :-( " << std::endl; continue; } |
| 218 |
> |
} |
| 219 |
|
|
| 220 |
|
cand.decayPoint = ttmd.crossingPoint(); |
| 221 |
|
cand.distance = ttmd.distance(); |
| 222 |
+ |
|
| 223 |
+ |
|
| 224 |
|
if (cand.distance > maxTwoTrackDist_) continue; |
| 225 |
|
|
| 226 |
|
typedef reco::Candidate::Vector Vector; |
| 227 |
+ |
typedef reco::Candidate::Point Point; |
| 228 |
|
typedef reco::Candidate::LorentzVector LorentzVector; |
| 229 |
< |
Vector ph = cand.innerTrack->outerMomentum(); // start from nearest states |
| 230 |
< |
Vector pm = cand.outerTrack->innerMomentum(); // so the material is correct |
| 229 |
> |
TSCPBuilderNoMaterial tscpBuilder; |
| 230 |
> |
Vector ph, pm; |
| 231 |
> |
if (!aodDataOnly_) { |
| 232 |
> |
TrajectoryStateOnSurface innerOuterTSOS = transformer.outerStateOnSurface(**it1, *geometry_, &*magfield_); // start from the two near states |
| 233 |
> |
TrajectoryStateOnSurface outerInnerTSOS = transformer.innerStateOnSurface(**it2, *geometry_, &*magfield_); // so the material is more correct |
| 234 |
> |
TrajectoryStateClosestToPoint innerTSCP = tscpBuilder(innerOuterTSOS, ttmd.crossingPoint()); |
| 235 |
> |
TrajectoryStateClosestToPoint outerTSCP = tscpBuilder(outerInnerTSOS, ttmd.crossingPoint()); |
| 236 |
> |
ph = Vector(innerTSCP.momentum()); |
| 237 |
> |
pm = Vector(outerTSCP.momentum()); |
| 238 |
> |
} else { |
| 239 |
> |
FreeTrajectoryState innerFTS = transformer.initialFreeState(**it1, &*magfield_); // start from the two far states |
| 240 |
> |
FreeTrajectoryState outerFTS = transformer.initialFreeState(**it2, &*magfield_); // so that distance is less biased by shared hits |
| 241 |
> |
TrajectoryStateClosestToPoint innerTSCP = tscpBuilder(innerFTS, ttmd.crossingPoint()); |
| 242 |
> |
TrajectoryStateClosestToPoint outerTSCP = tscpBuilder(outerFTS, ttmd.crossingPoint()); |
| 243 |
> |
ph = Vector(innerTSCP.momentum()); |
| 244 |
> |
pm = Vector(outerTSCP.momentum()); |
| 245 |
> |
} |
| 246 |
|
Vector pn = ph-pm; |
| 247 |
|
cand.p4mu = LorentzVector(pm.x(), pm.y(), pm.z(), hypot(pm.r(), 0.10566)); |
| 248 |
|
cand.p4nu = LorentzVector(pn.x(), pn.y(), pn.z(), pn.r()); |
| 281 |
|
std::cerr << " theta* = " << cand.thetaStar << std::endl; |
| 282 |
|
|
| 283 |
|
decayCandidates.push_back(cand); |
| 284 |
+ |
innerSelTks->push_back(*cand.innerTrack); |
| 285 |
+ |
outerSelTks->push_back(*cand.outerTrack); |
| 286 |
+ |
if (cand.innerMu.isNonnull()) innerSelMus->push_back(*cand.innerMu); |
| 287 |
+ |
if (cand.outerMu.isNonnull()) outerSelMus->push_back(*cand.outerMu); |
| 288 |
+ |
|
| 289 |
+ |
output->push_back(pat::GenericParticle()); |
| 290 |
+ |
pat::GenericParticle &gp = output->back(); |
| 291 |
+ |
/// ----- Base kinematics ----- |
| 292 |
+ |
gp.setCharge(cand.innerTrack->charge()); |
| 293 |
+ |
gp.setP4(math::PtEtaPhiMLorentzVector(cand.innerTrack->pt(), |
| 294 |
+ |
cand.innerTrack->eta(), |
| 295 |
+ |
cand.innerTrack->phi(), |
| 296 |
+ |
(cand.p4had.mass2() > 0 ? cand.p4had.mass() : -sqrt(-cand.p4had.mass2())))); |
| 297 |
+ |
gp.setVertex(cand.innerTrack->vertex()); |
| 298 |
+ |
/// ----- References ------ |
| 299 |
+ |
gp.setTrack(cand.innerTrack); |
| 300 |
+ |
reco::TrackRefVector tks; |
| 301 |
+ |
tks.push_back(cand.innerTrack); |
| 302 |
+ |
tks.push_back(cand.outerTrack); |
| 303 |
+ |
gp.setTracks(tks); |
| 304 |
+ |
if (cand.outerMu.isNonnull() && cand.outerMu->isStandAloneMuon()) { |
| 305 |
+ |
gp.setStandAloneMuon(cand.outerMu->standAloneMuon()); |
| 306 |
+ |
} |
| 307 |
+ |
if (cand.outerMu.isNonnull() && cand.outerMu->isGlobalMuon()) { |
| 308 |
+ |
gp.setCombinedMuon(cand.outerMu->combinedMuon()); |
| 309 |
+ |
} |
| 310 |
+ |
gp.addUserCand("innerMu", cand.innerMu); |
| 311 |
+ |
gp.addUserCand("outerMu", cand.outerMu); |
| 312 |
+ |
/// ---- Decay kinematics ---- |
| 313 |
+ |
gp.addUserFloat("distance", cand.distance); |
| 314 |
+ |
gp.addUserFloat("decayRho", cand.decayPoint.perp()); |
| 315 |
+ |
gp.addUserFloat("decayZ", cand.decayPoint.z()); |
| 316 |
+ |
gp.addUserData( "decayPoint", Point(cand.decayPoint)); |
| 317 |
+ |
gp.addUserFloat("theta", cand.theta); |
| 318 |
+ |
gp.addUserFloat("thetaStar", cand.thetaStar); |
| 319 |
+ |
/// ---- Impact parameter stuff ---- |
| 320 |
+ |
// first the basics |
| 321 |
+ |
reco::VertexRef vtx(vertices, 0); // primary vertex |
| 322 |
+ |
gp.addUserFloat("innerDxyVtx",cand.innerTrack->dxy(vtx->position())); |
| 323 |
+ |
gp.addUserFloat("outerDxyVtx",cand.outerTrack->dxy(vtx->position())); |
| 324 |
+ |
gp.addUserFloat("innerDszVtx",cand.innerTrack->dsz(vtx->position())); |
| 325 |
+ |
gp.addUserFloat("outerDszVtx",cand.outerTrack->dsz(vtx->position())); |
| 326 |
+ |
// the we get fancy |
| 327 |
+ |
GlobalVector p4atVtx(cand.innerTrack->px(), cand.innerTrack->py(), cand.innerTrack->pz()); |
| 328 |
+ |
reco::TransientTrack innerTT = ttBuilder_->build(*cand.innerTrack); |
| 329 |
+ |
reco::TransientTrack outerTT = ttBuilder_->build(*cand.outerTrack); |
| 330 |
+ |
std::pair<bool,Measurement1D> innerTIP = IPTools::signedTransverseImpactParameter(innerTT, p4atVtx, *vtx); |
| 331 |
+ |
std::pair<bool,Measurement1D> inner3DIP = IPTools::signedImpactParameter3D(innerTT, p4atVtx, *vtx); |
| 332 |
+ |
std::pair<bool,Measurement1D> innerDL3D = IPTools::signedDecayLength3D(innerTT, p4atVtx, *vtx); |
| 333 |
+ |
std::pair<bool,Measurement1D> outerTIP = IPTools::signedTransverseImpactParameter(outerTT, p4atVtx, *vtx); |
| 334 |
+ |
std::pair<bool,Measurement1D> outer3DIP = IPTools::signedImpactParameter3D(outerTT, p4atVtx, *vtx); |
| 335 |
+ |
std::pair<bool,Measurement1D> outerDL3D = IPTools::signedDecayLength3D(outerTT, p4atVtx, *vtx); |
| 336 |
+ |
if (innerTIP.first ) { gp.addUserFloat("innerTIP", innerTIP.second.value()); gp.addUserFloat("innerTIPsig", innerTIP.second.significance()); } |
| 337 |
+ |
if (inner3DIP.first) { gp.addUserFloat("inner3DIP",inner3DIP.second.value()); gp.addUserFloat("inner3DIPsig",inner3DIP.second.significance()); } |
| 338 |
+ |
if (innerDL3D.first) { gp.addUserFloat("innerDL3D",innerDL3D.second.value()); gp.addUserFloat("innerDL3Dsig",innerDL3D.second.significance()); } |
| 339 |
+ |
if (outerTIP.first ) { gp.addUserFloat("outerTIP", outerTIP.second.value()); gp.addUserFloat("outerTIPsig", outerTIP.second.significance()); } |
| 340 |
+ |
if (outer3DIP.first) { gp.addUserFloat("outer3DIP",outer3DIP.second.value()); gp.addUserFloat("outer3DIPsig",outer3DIP.second.significance()); } |
| 341 |
+ |
if (outerDL3D.first) { gp.addUserFloat("outerDL3D",outerDL3D.second.value()); gp.addUserFloat("outerDL3Dsig",outerDL3D.second.significance()); } |
| 342 |
|
} |
| 343 |
|
} |
| 344 |
|
} |
| 345 |
+ |
|
| 346 |
+ |
iEvent.put(innerSelTks, "inner"); |
| 347 |
+ |
iEvent.put(outerSelTks, "outer"); |
| 348 |
+ |
iEvent.put(innerSelMus, "inner"); |
| 349 |
+ |
iEvent.put(outerSelMus, "outer"); |
| 350 |
+ |
iEvent.put(output); |
| 351 |
|
} |
| 352 |
|
|
| 353 |
|
#include "FWCore/Framework/interface/MakerMacros.h" |
