19 #include "Math/VectorUtil.h" 33 higherPuritySelection_(iConfig.getParameter<
std::
string>(
"higherPuritySelection")),
34 lowerPuritySelection_(iConfig.getParameter<
std::
string>(
"lowerPuritySelection")),
36 iConfig.existsAs<
std::
string>(
"dimuonSelection") ? iConfig.getParameter<
std::
string>(
"dimuonSelection") :
""),
37 addCommonVertex_(iConfig.getParameter<
bool>(
"addCommonVertex")),
38 addMuonlessPrimaryVertex_(iConfig.getParameter<
bool>(
"addMuonlessPrimaryVertex")),
39 resolveAmbiguity_(iConfig.getParameter<
bool>(
"resolvePileUpAmbiguity")),
40 addMCTruth_(iConfig.getParameter<
bool>(
"addMCTruth")) {
44 produces<pat::CompositeCandidateCollection>();
63 vector<double> muMasses;
64 muMasses.push_back(0.1056583715);
65 muMasses.push_back(0.1056583715);
79 theBeamSpotV =
Vertex(bs.position(), bs.covariance3D());
83 if (priVtxs->begin() != priVtxs->end()) {
84 thePrimaryV =
Vertex(*(priVtxs->begin()));
86 thePrimaryV =
Vertex(bs.position(), bs.covariance3D());
111 vector<TransientVertex>
pvs;
118 LorentzVector
jpsi = it->p4() + it2->p4();
120 myCand.
setCharge(it->charge() + it2->charge());
127 if (it->track().isNonnull() && it2->track().isNonnull()) {
129 vector<TransientTrack> t_tks;
130 t_tks.push_back(theTTBuilder->
build(
132 t_tks.push_back(theTTBuilder->
build(*it2->track()));
149 float vProb(TMath::Prob(vChi2, (
int)vNDF));
159 TVector3 pperp(jpsi.px(), jpsi.py(), 0);
163 float minDz = 999999.;
175 float extrapZ = -9E20;
177 extrapZ = ttmd.
points().first.z();
179 for (VertexCollection::const_iterator itv = priVtxs->begin(), itvend = priVtxs->end(); itv != itvend;
181 float deltaZ = fabs(extrapZ - itv->position().z());
182 if (deltaZ < minDz) {
184 thePrimaryV =
Vertex(*itv);
189 Vertex theOriginalPV = thePrimaryV;
199 std::cout <<
"pvtracks NOT valid " << std::endl;
203 if (pvbeamspot.
id() != theBeamSpot.
id())
205 <<
"The BeamSpot used for PV reco is not the same used in this analyzer.";
211 if (rmu1 !=
nullptr && rmu2 !=
nullptr && rmu1->
track().
id() == pvtracks.
id() &&
218 std::vector<reco::Track>::const_iterator itRefittedTrack = thePrimaryV.
refittedTracks().begin();
219 std::vector<reco::Track>::const_iterator refittedTracksEnd = thePrimaryV.
refittedTracks().end();
220 for (; itRefittedTrack != refittedTracksEnd; ++itRefittedTrack) {
226 muonLess.push_back(*(thePrimaryV.
originalTrack(*itRefittedTrack)));
229 std::vector<reco::TrackBaseRef>::const_iterator itPVtrack = thePrimaryV.
tracks_begin();
230 for (; itPVtrack != thePrimaryV.
tracks_end(); ++itPVtrack)
231 if (itPVtrack->isNonnull()) {
232 if (itPVtrack->key() == rmu1->
track().
key())
234 if (itPVtrack->key() == rmu2->
track().
key())
237 muonLess.push_back(**itPVtrack);
240 if (muonLess.size() > 1 && muonLess.size() < thePrimaryV.
tracksSize()) {
241 pvs = revertex.
makeVertices(muonLess, *pvbeamspot, iSetup);
244 thePrimaryV = muonLessPV;
252 double vertexWeight = -1., sumPTPV = -1.;
253 int countTksOfPV = -1;
260 if (itVtx->isNonnull()) {
264 if (track.
pt() < 0.5)
270 if (tkPVdist.second.significance() > 3)
275 if (rmu1 !=
nullptr && rmu1->
innerTrack().
key() == itVtx->key())
277 if (rmu2 !=
nullptr && rmu2->
innerTrack().
key() == itVtx->key())
283 sumPTPV += track.
pt();
287 std::cout <<
" muon Selection%G�%@failed " << std::endl;
291 myCand.
addUserInt(
"countTksOfPV", countTksOfPV);
292 myCand.
addUserFloat(
"vertexWeight", (
float)vertexWeight);
315 TVector3 vdiff = vtx - pvtx;
316 double cosAlpha = vdiff.Dot(pperp) / (vdiff.Perp() * pperp.Perp());
319 double ctauPV = distXY.
value() * cosAlpha * myCand.
mass() / pperp.Perp();
324 double ctauErrPV =
sqrt(ROOT::Math::Similarity(vpperp, vXYe)) * myCand.
mass() / (pperp.Perp2());
333 cosAlpha = vdiff.Dot(pperp) / (vdiff.Perp() * pperp.Perp());
336 double ctauBS = distXY.
value() * cosAlpha * myCand.
mass() / pperp.Perp();
341 double ctauErrBS =
sqrt(ROOT::Math::Similarity(vpperp, vXYeB)) * myCand.
mass() / (pperp.Perp2());
374 if (genMu1->numberOfMothers() > 0 && genMu2->numberOfMothers() > 0) {
377 if (mom1.
isNonnull() && (mom1 == mom2)) {
390 consumes<reco::GenParticleCollection>((
edm::InputTag)
"genParticles");
391 iEvent.
getByToken(genCands_, theGenParticles);
392 if (theGenParticles.
isValid()) {
393 for (
size_t iGenParticle = 0; iGenParticle < theGenParticles->size(); ++iGenParticle) {
394 const Candidate &genCand = (*theGenParticles)[iGenParticle];
395 if (genCand.
pdgId() == 443 || genCand.
pdgId() == 100443 || genCand.
pdgId() == 553 ||
396 genCand.
pdgId() == 100553 || genCand.
pdgId() == 200553) {
417 oniaOutput->push_back(myCand);
421 std::sort(oniaOutput->begin(), oniaOutput->end(),
vPComparator_);
427 if (
abs(pdgID) == 511 ||
abs(pdgID) == 521 ||
abs(pdgID) == 531 ||
abs(pdgID) == 5122)
433 if ((
abs(pdgID) == 511 &&
abs(momPdgID) == 511 && pdgID * momPdgID < 0) ||
434 (
abs(pdgID) == 531 &&
abs(momPdgID) == 531 && pdgID * momPdgID < 0))
441 float trueLife = -99.;
443 if (genJpsi->numberOfMothers() > 0) {
444 TVector3 trueVtx(0.0, 0.0, 0.0);
445 TVector3 trueP(0.0, 0.0, 0.0);
446 TVector3 trueVtxMom(0.0, 0.0, 0.0);
448 trueVtx.SetXYZ(genJpsi->vertex().x(), genJpsi->vertex().y(), genJpsi->vertex().z());
449 trueP.SetXYZ(genJpsi->momentum().x(), genJpsi->momentum().y(), genJpsi->momentum().z());
451 bool aBhadron =
false;
454 std::pair<int, float>
result = std::make_pair(momJpsiID, trueLife);
460 momJpsiID = Jpsigrandmom->pdgId();
461 trueVtxMom.SetXYZ(Jpsigrandmom->vertex().x(), Jpsigrandmom->vertex().y(), Jpsigrandmom->vertex().z());
463 momJpsiID = Jpsimom->pdgId();
464 trueVtxMom.SetXYZ(Jpsimom->vertex().x(), Jpsimom->vertex().y(), Jpsimom->vertex().z());
471 momJpsiID = JpsiGrandgrandmom->pdgId();
473 JpsiGrandgrandmom->vertex().x(), JpsiGrandgrandmom->vertex().y(), JpsiGrandgrandmom->vertex().z());
475 momJpsiID = Jpsigrandmom->pdgId();
476 trueVtxMom.SetXYZ(Jpsigrandmom->vertex().x(), Jpsigrandmom->vertex().y(), Jpsigrandmom->vertex().z());
482 momJpsiID = Jpsimom->pdgId();
483 trueVtxMom.SetXYZ(Jpsimom->vertex().x(), Jpsimom->vertex().y(), Jpsimom->vertex().z());
487 TVector3 vdiff = trueVtx - trueVtxMom;
489 trueLife = vdiff.Perp() * genJpsi->mass() / trueP.Perp();
491 std::pair<int, float>
result = std::make_pair(momJpsiID, trueLife);
Analysis-level particle class.
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
float distance() const override
TrackRef track() const override
reference to a Track
bool isNonnull() const
Checks for non-null.
std::pair< GlobalPoint, GlobalPoint > points() const override
trackRef_iterator tracks_end() const
last iterator over tracks
StringCutObjectSelector< pat::Muon > higherPuritySelection_
bool isAbHadron(int pdgID)
virtual TrackRef innerTrack() const
int nominalValue() const
The nominal field value for this map in kGauss.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
TrackBaseRef originalTrack(const Track &refTrack) const
const FreeTrajectoryState & theState() const
const AlgebraicSymMatrix33 matrix() const
Onia2MuMuPAT(const edm::ParameterSet &)
double px() const final
x coordinate of momentum vector
bool hasRefittedTracks() const
Checks whether refitted tracks are stored.
float totalChiSquared() const
Global3DPoint GlobalPoint
std::vector< Track > TrackCollection
collection of Tracks
reco::TransientTrack build(const reco::Track *p) const
std::pair< int, float > findJpsiMCInfo(reco::GenParticleRef genJpsi)
std::vector< Vertex > VertexCollection
collection of Vertex objects
InvariantMassFromVertex massCalculator
void addUserFloat(const std::string &label, float data, const bool overwrite=false)
Set user-defined float.
key_type key() const
Accessor for product key.
void produce(edm::Event &, const edm::EventSetup &) override
const std::vector< Track > & refittedTracks() const
Returns the container of refitted tracks.
const Point & position() const
position
void setGenParticleRef(const reco::GenParticleRef &ref, bool embed=false)
Set the generator level particle reference.
ProductID id() const
Accessor for product ID.
Measurement1D distance(const GlobalPoint &vtx1Position, const GlobalError &vtx1PositionError, const GlobalPoint &vtx2Position, const GlobalError &vtx2PositionError) const override
bool calculate(const TrajectoryStateOnSurface &sta, const TrajectoryStateOnSurface &stb) override
void setCharge(Charge q) final
set electric charge
#define DEFINE_FWK_MODULE(type)
bool status() const override
edm::EDGetTokenT< reco::BeamSpot > revtxbs_
float degreesOfFreedom() const
double pz() const final
z coordinate of momentum vector
bool calculate(const TrajectoryStateOnSurface &sta, const TrajectoryStateOnSurface &stb) override
virtual int pdgId() const =0
PDG identifier.
GlobalPoint position() const
double pt() const
track transverse momentum
double ptError() const
error on Pt (set to 1000 TeV if charge==0 for safety)
StringCutObjectSelector< reco::Candidate, true > dimuonSelection_
Abs< T >::type abs(const T &t)
std::vector< TransientVertex > makeVertices(const reco::TrackCollection &tracks, const reco::BeamSpot &bs, const edm::EventSetup &iSetup) const
Make the vertices.
CachingVertex< 5 > vertex(const std::vector< reco::TransientTrack > &tracks) const override
edm::EDGetTokenT< reco::TrackCollection > revtxtrks_
float trackWeight(const TREF &r) const
returns the weight with which a Track has contributed to the vertex-fit.
bool addMuonlessPrimaryVertex_
bool isNull() const
Checks for null.
void addUserInt(const std::string &label, int32_t data, const bool overwrite=false)
Set user-defined int.
void addDaughter(const Candidate &, const std::string &s="")
add a clone of the passed candidate as daughter
bool isAMixedbHadron(int pdgID, int momPdgID)
StringCutObjectSelector< pat::Muon > lowerPuritySelection_
std::vector< CompositeCandidate > CompositeCandidateCollection
double py() const final
y coordinate of momentum vector
Error error() const
return SMatrix
edm::EDGetTokenT< reco::BeamSpot > thebeamspot_
bool quality(const TrackQuality) const
Track quality.
math::XYZTLorentzVector LorentzVector
Lorentz vector.
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
ROOT::Math::SVector< double, 3 > AlgebraicVector3
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepSym< double, 3 > > AlgebraicSymMatrix33
edm::EDGetTokenT< reco::VertexCollection > thePVs_
GreaterByVProb< pat::CompositeCandidate > vPComparator_
trackRef_iterator tracks_begin() const
first iterator over tracks
T const * product() const
void setP4(const LorentzVector &p4) final
set 4-momentum
void addUserData(const std::string &label, const T &data, bool transientOnly=false, bool overwrite=false)
std::vector< TrackBaseRef >::const_iterator trackRef_iterator
The iteratator for the vector<TrackRef>
size_t tracksSize() const
number of tracks
Measurement1D invariantMass(const CachingVertex< 5 > &vertex, const std::vector< double > &masses) const
edm::EDGetTokenT< edm::View< pat::Muon > > muons_
math::PtEtaPhiELorentzVectorF LorentzVector
Global3DVector GlobalVector
double mass() const final
mass