#include <Onia2MuMuPAT.h>

Inheritance diagram for Onia2MuMuPAT:

Public Member Functions
	Onia2MuMuPAT (const edm::ParameterSet &)

Public Member Functions inherited from edm::stream::EDProducer<>
	EDProducer ()=default

	EDProducer (const EDProducer &)=delete

bool	hasAbilityToProduceInBeginLumis () const final

bool	hasAbilityToProduceInBeginProcessBlocks () const final

bool	hasAbilityToProduceInBeginRuns () const final

bool	hasAbilityToProduceInEndLumis () const final

bool	hasAbilityToProduceInEndProcessBlocks () const final

bool	hasAbilityToProduceInEndRuns () const final

const EDProducer &	operator= (const EDProducer &)=delete

Static Public Member Functions
static void	fillDescriptions (edm::ConfigurationDescriptions &)

Private Member Functions
std::pair< int, float >	findJpsiMCInfo (reco::GenParticleRef genJpsi) const

bool	isAbHadron (int pdgID) const

bool	isAMixedbHadron (int pdgID, int momPdgID) const

void	produce (edm::Event &, const edm::EventSetup &) override

Private Attributes
bool	addCommonVertex_

bool	addMCTruth_

bool	addMuonlessPrimaryVertex_

StringCutObjectSelector< reco::Candidate, true >	dimuonSelection_

StringCutObjectSelector< pat::Muon >	higherPuritySelection_

StringCutObjectSelector< pat::Muon >	lowerPuritySelection_

edm::ESGetToken< MagneticField, IdealMagneticFieldRecord >	magneticFieldToken_

InvariantMassFromVertex	massCalculator

edm::EDGetTokenT< edm::View< pat::Muon > >	muons_

bool	resolveAmbiguity_

edm::EDGetTokenT< reco::BeamSpot >	revtxbs_

edm::EDGetTokenT< reco::TrackCollection >	revtxtrks_

edm::EDGetTokenT< reco::BeamSpot >	thebeamspot_

edm::EDGetTokenT< reco::VertexCollection >	thePVs_

edm::ESGetToken< TransientTrackBuilder, TransientTrackRecord >	theTTBuilderToken_

GreaterByVProb< pat::CompositeCandidate >	vPComparator_

Additional Inherited Members
Public Types inherited from edm::stream::EDProducer<>
using	CacheTypes = CacheContexts< T... >

using	GlobalCache = typename CacheTypes::GlobalCache

using	HasAbility = AbilityChecker< T... >

using	InputProcessBlockCache = typename CacheTypes::InputProcessBlockCache

using	LuminosityBlockCache = typename CacheTypes::LuminosityBlockCache

using	LuminosityBlockContext = LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >

using	LuminosityBlockSummaryCache = typename CacheTypes::LuminosityBlockSummaryCache

using	RunCache = typename CacheTypes::RunCache

using	RunContext = RunContextT< RunCache, GlobalCache >

using	RunSummaryCache = typename CacheTypes::RunSummaryCache

Detailed Description

Definition at line 39 of file Onia2MuMuPAT.h.

Constructor & Destructor Documentation

◆ Onia2MuMuPAT()

Onia2MuMuPAT::Onia2MuMuPAT ( const edm::ParameterSet & iConfig )

explicit

Definition at line 26 of file Onia2MuMuPAT.cc.

References ProducerED_cfi::InputTag, revtxbs_, and revtxtrks_.

     : muons_(consumes<edm::View<pat::Muon>>(iConfig.getParameter<edm::InputTag>("muons"))),
       thebeamspot_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpotTag"))),
       thePVs_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("primaryVertexTag"))),
       magneticFieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
       theTTBuilderToken_(
           esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))),
       higherPuritySelection_(iConfig.getParameter<std::string>("higherPuritySelection")),
       lowerPuritySelection_(iConfig.getParameter<std::string>("lowerPuritySelection")),
       dimuonSelection_(iConfig.getParameter<std::string>("dimuonSelection")),
       addCommonVertex_(iConfig.getParameter<bool>("addCommonVertex")),
       addMuonlessPrimaryVertex_(iConfig.getParameter<bool>("addMuonlessPrimaryVertex")),
       resolveAmbiguity_(iConfig.getParameter<bool>("resolvePileUpAmbiguity")),
       addMCTruth_(iConfig.getParameter<bool>("addMCTruth")) {
   revtxtrks_ = consumes<reco::TrackCollection>(
       (edm::InputTag) "generalTracks");  //if that is not true, we will raise an exception
   revtxbs_ = consumes<reco::BeamSpot>((edm::InputTag) "offlineBeamSpot");
   produces<pat::CompositeCandidateCollection>();
 }

Member Function Documentation

◆ fillDescriptions()

void Onia2MuMuPAT::fillDescriptions ( edm::ConfigurationDescriptions & iDescriptions )

static

Definition at line 486 of file Onia2MuMuPAT.cc.

References edm::ConfigurationDescriptions::addDefault(), submitPVResolutionJobs::desc, and AlCaHLTBitMon_QueryRunRegistry::string.

                                                                                {
   edm::ParameterSetDescription desc;
   desc.add<edm::InputTag>("muons");
   desc.add<edm::InputTag>("beamSpotTag");
   desc.add<edm::InputTag>("primaryVertexTag");
   desc.add<std::string>("higherPuritySelection");
   desc.add<std::string>("lowerPuritySelection");
   desc.add<std::string>("dimuonSelection", "");
   desc.add<bool>("addCommonVertex");
   desc.add<bool>("addMuonlessPrimaryVertex");
   desc.add<bool>("resolvePileUpAmbiguity");
   desc.add<bool>("addMCTruth");
 
   iDescriptions.addDefault(desc);
 }

◆ findJpsiMCInfo()

std::pair< int, float > Onia2MuMuPAT::findJpsiMCInfo ( reco::GenParticleRef genJpsi ) const

private

Definition at line 430 of file Onia2MuMuPAT.cc.

References isAbHadron(), isAMixedbHadron(), edm::Ref< C, T, F >::isNonnull(), edm::Ref< C, T, F >::isNull(), and mps_fire::result.

Referenced by produce().

                                                                                  {
   int momJpsiID = 0;
   float trueLife = -99.;
 
   if (genJpsi->numberOfMothers() > 0) {
     TVector3 trueVtx(0.0, 0.0, 0.0);
     TVector3 trueP(0.0, 0.0, 0.0);
     TVector3 trueVtxMom(0.0, 0.0, 0.0);
 
     trueVtx.SetXYZ(genJpsi->vertex().x(), genJpsi->vertex().y(), genJpsi->vertex().z());
     trueP.SetXYZ(genJpsi->momentum().x(), genJpsi->momentum().y(), genJpsi->momentum().z());
 
     bool aBhadron = false;
     reco::GenParticleRef Jpsimom = genJpsi->motherRef();  // find mothers
     if (Jpsimom.isNull()) {
       std::pair<int, float> result = std::make_pair(momJpsiID, trueLife);
       return result;
     } else {
       reco::GenParticleRef Jpsigrandmom = Jpsimom->motherRef();
       if (isAbHadron(Jpsimom->pdgId())) {
         if (Jpsigrandmom.isNonnull() && isAMixedbHadron(Jpsimom->pdgId(), Jpsigrandmom->pdgId())) {
           momJpsiID = Jpsigrandmom->pdgId();
           trueVtxMom.SetXYZ(Jpsigrandmom->vertex().x(), Jpsigrandmom->vertex().y(), Jpsigrandmom->vertex().z());
         } else {
           momJpsiID = Jpsimom->pdgId();
           trueVtxMom.SetXYZ(Jpsimom->vertex().x(), Jpsimom->vertex().y(), Jpsimom->vertex().z());
         }
         aBhadron = true;
       } else {
         if (Jpsigrandmom.isNonnull() && isAbHadron(Jpsigrandmom->pdgId())) {
           reco::GenParticleRef JpsiGrandgrandmom = Jpsigrandmom->motherRef();
           if (JpsiGrandgrandmom.isNonnull() && isAMixedbHadron(Jpsigrandmom->pdgId(), JpsiGrandgrandmom->pdgId())) {
             momJpsiID = JpsiGrandgrandmom->pdgId();
             trueVtxMom.SetXYZ(
                 JpsiGrandgrandmom->vertex().x(), JpsiGrandgrandmom->vertex().y(), JpsiGrandgrandmom->vertex().z());
           } else {
             momJpsiID = Jpsigrandmom->pdgId();
             trueVtxMom.SetXYZ(Jpsigrandmom->vertex().x(), Jpsigrandmom->vertex().y(), Jpsigrandmom->vertex().z());
           }
           aBhadron = true;
         }
       }
       if (!aBhadron) {
         momJpsiID = Jpsimom->pdgId();
         trueVtxMom.SetXYZ(Jpsimom->vertex().x(), Jpsimom->vertex().y(), Jpsimom->vertex().z());
       }
     }
 
     TVector3 vdiff = trueVtx - trueVtxMom;
     //trueLife = vdiff.Perp()*3.09688/trueP.Perp();
     trueLife = vdiff.Perp() * genJpsi->mass() / trueP.Perp();
   }
   std::pair<int, float> result = std::make_pair(momJpsiID, trueLife);
   return result;
 }

◆ isAbHadron()

bool Onia2MuMuPAT::isAbHadron ( int pdgID ) const

private

Definition at line 417 of file Onia2MuMuPAT.cc.

References funct::abs().

Referenced by findJpsiMCInfo().

                                              {
   if (abs(pdgID) == 511 || abs(pdgID) == 521 || abs(pdgID) == 531 || abs(pdgID) == 5122)
     return true;
   return false;
 }

◆ isAMixedbHadron()

bool Onia2MuMuPAT::isAMixedbHadron	(	int	pdgID,
		int	momPdgID
	)		const

private

Definition at line 423 of file Onia2MuMuPAT.cc.

References funct::abs().

Referenced by findJpsiMCInfo().

                                                                 {
   if ((abs(pdgID) == 511 && abs(momPdgID) == 511 && pdgID * momPdgID < 0) ||
       (abs(pdgID) == 531 && abs(momPdgID) == 531 && pdgID * momPdgID < 0))
     return true;
   return false;
 }

◆ produce()

void Onia2MuMuPAT::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

overrideprivate

DCA

end DCA

Definition at line 51 of file Onia2MuMuPAT.cc.

References IPTools::absoluteImpactParameter3D(), addCommonVertex_, reco::CompositeCandidate::addDaughter(), addMCTruth_, addMuonlessPrimaryVertex_, pat::PATObject< ObjectType >::addUserData(), pat::PATObject< ObjectType >::addUserFloat(), pat::PATObject< ObjectType >::addUserInt(), align::BeamSpot, cms::cuda::bs, TransientTrackBuilder::build(), TwoTrackMinimumDistance::calculate(), ClosestApproachInRPhi::calculate(), gather_cfg::cout, TransientVertex::degreesOfFreedom(), l1tTrackerHTMiss_cfi::deltaZ, dimuonSelection_, VertexDistanceXY::distance(), ClosestApproachInRPhi::distance(), pat::PATObject< ObjectType >::embedGenParticle(), submitPVResolutionJobs::err, Measurement1D::error(), relativeConstraints::error, reco::Vertex::error(), cppFunctionSkipper::exception, findJpsiMCInfo(), edm::EventSetup::getData(), edm::EventSetup::getHandle(), reco::Vertex::hasRefittedTracks(), higherPuritySelection_, reco::TrackBase::highPurity, edm::HandleBase::id(), edm::Ref< C, T, F >::id(), iEvent, reco::Muon::innerTrack(), ProducerED_cfi::InputTag, InvariantMassFromVertex::invariantMass(), edm::Ref< C, T, F >::isNonnull(), edm::HandleBase::isValid(), TrajectoryStateClosestToPoint::isValid(), TransientVertex::isValid(), ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, edm::RefToBase< T >::key(), edm::Ref< C, T, F >::key(), lowerPuritySelection_, magneticFieldToken_, OniaVtxReProducer::makeVertices(), reco::LeafCandidate::mass(), massCalculator, GlobalErrorBase< T, ErrorWeightType >::matrix(), pfMETCorrectionType0_cfi::minDz, eostools::move(), DiMuonV_cfg::muons, muons_, MagneticField::nominalValue(), reco::Vertex::originalTrack(), reco::Candidate::pdgId(), TwoTrackMinimumDistance::points(), reco::Vertex::position(), TransientVertex::position(), FSQDQM_cfi::pvs, reco::LeafCandidate::px(), reco::LeafCandidate::py(), reco::LeafCandidate::pz(), reco::Vertex::refittedTracks(), reco::Vertex::reserve(), resolveAmbiguity_, revtxbs_, revtxtrks_, reco::LeafCandidate::setCharge(), pat::PATObject< ObjectType >::setGenParticleRef(), reco::LeafCandidate::setP4(), jetUpdater_cfi::sort, mathSSE::sqrt(), ClosestApproachInRPhi::status(), mps_update::status, thebeamspot_, thePVs_, TrajectoryStateClosestToPoint::theState(), theTTBuilderToken_, TransientVertex::totalChiSquared(), reco::Muon::track(), HLT_2023v12_cff::track, reco::Vertex::tracks_begin(), reco::Vertex::tracks_end(), reco::Vertex::tracksSize(), reco::Vertex::trackWeight(), Measurement1D::value(), HltBtagValidation_cff::Vertex, KalmanVertexFitter::vertex(), vPComparator_, L1BJetProducer_cff::vtx, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                                         {
   using namespace edm;
   using namespace std;
   using namespace reco;
   typedef Candidate::LorentzVector LorentzVector;
 
   vector<double> muMasses;
   muMasses.push_back(0.1056583715);
   muMasses.push_back(0.1056583715);
 
   std::unique_ptr<pat::CompositeCandidateCollection> oniaOutput(new pat::CompositeCandidateCollection);
 
   Vertex thePrimaryV;
   Vertex theBeamSpotV;
 
   const MagneticField &field = iSetup.getData(magneticFieldToken_);
 
   Handle<BeamSpot> theBeamSpot;
   iEvent.getByToken(thebeamspot_, theBeamSpot);
   BeamSpot bs = *theBeamSpot;
   theBeamSpotV = Vertex(bs.position(), bs.covariance3D());
 
   Handle<VertexCollection> priVtxs;
   iEvent.getByToken(thePVs_, priVtxs);
   if (priVtxs->begin() != priVtxs->end()) {
     thePrimaryV = Vertex(*(priVtxs->begin()));
   } else {
     thePrimaryV = Vertex(bs.position(), bs.covariance3D());
   }
 
   Handle<View<pat::Muon>> muons;
   iEvent.getByToken(muons_, muons);
 
   edm::ESHandle<TransientTrackBuilder> theTTBuilder = iSetup.getHandle(theTTBuilderToken_);
   KalmanVertexFitter vtxFitter(true);
   TrackCollection muonLess;
 
   // JPsi candidates only from muons
   for (View<pat::Muon>::const_iterator it = muons->begin(), itend = muons->end(); it != itend; ++it) {
     // both must pass low quality
     if (!lowerPuritySelection_(*it))
       continue;
     for (View<pat::Muon>::const_iterator it2 = it + 1; it2 != itend; ++it2) {
       // both must pass low quality
       if (!lowerPuritySelection_(*it2))
         continue;
       // one must pass tight quality
       if (!(higherPuritySelection_(*it) || higherPuritySelection_(*it2)))
         continue;
 
       pat::CompositeCandidate myCand;
       vector<TransientVertex> pvs;
 
       // ---- no explicit order defined ----
       myCand.addDaughter(*it, "muon1");
       myCand.addDaughter(*it2, "muon2");
 
       // ---- define and set candidate's 4momentum  ----
       LorentzVector jpsi = it->p4() + it2->p4();
       myCand.setP4(jpsi);
       myCand.setCharge(it->charge() + it2->charge());
 
       // ---- apply the dimuon cut ----
       if (!dimuonSelection_(myCand))
         continue;
 
       // ---- fit vertex using Tracker tracks (if they have tracks) ----
       if (it->track().isNonnull() && it2->track().isNonnull()) {
         //build the dimuon secondary vertex
         vector<TransientTrack> t_tks;
         t_tks.push_back(theTTBuilder->build(
             *it->track()));  // pass the reco::Track, not  the reco::TrackRef (which can be transient)
         t_tks.push_back(theTTBuilder->build(*it2->track()));  // otherwise the vertex will have transient refs inside.
         TransientVertex myVertex = vtxFitter.vertex(t_tks);
 
         CachingVertex<5> VtxForInvMass = vtxFitter.vertex(t_tks);
 
         Measurement1D MassWErr(jpsi.M(), -9999.);
         if (field.nominalValue() > 0) {
           MassWErr = massCalculator.invariantMass(VtxForInvMass, muMasses);
         } else {
           myVertex = TransientVertex();  // with no arguments it is invalid
         }
 
         myCand.addUserFloat("MassErr", MassWErr.error());
 
         if (myVertex.isValid()) {
           float vChi2 = myVertex.totalChiSquared();
           float vNDF = myVertex.degreesOfFreedom();
           float vProb(TMath::Prob(vChi2, (int)vNDF));
 
           myCand.addUserFloat("vNChi2", vChi2 / vNDF);
           myCand.addUserFloat("vProb", vProb);
 
           TVector3 vtx;
           TVector3 pvtx;
           VertexDistanceXY vdistXY;
 
           vtx.SetXYZ(myVertex.position().x(), myVertex.position().y(), 0);
           TVector3 pperp(jpsi.px(), jpsi.py(), 0);
           AlgebraicVector3 vpperp(pperp.x(), pperp.y(), 0);
 
           if (resolveAmbiguity_) {
             float minDz = 999999.;
             TwoTrackMinimumDistance ttmd;
             bool status = ttmd.calculate(
                 GlobalTrajectoryParameters(
                     GlobalPoint(myVertex.position().x(), myVertex.position().y(), myVertex.position().z()),
                     GlobalVector(myCand.px(), myCand.py(), myCand.pz()),
                     TrackCharge(0),
                     &(field)),
                 GlobalTrajectoryParameters(GlobalPoint(bs.position().x(), bs.position().y(), bs.position().z()),
                                            GlobalVector(bs.dxdz(), bs.dydz(), 1.),
                                            TrackCharge(0),
                                            &(field)));
             float extrapZ = -9E20;
             if (status)
               extrapZ = ttmd.points().first.z();
 
             for (VertexCollection::const_iterator itv = priVtxs->begin(), itvend = priVtxs->end(); itv != itvend;
                  ++itv) {
               float deltaZ = fabs(extrapZ - itv->position().z());
               if (deltaZ < minDz) {
                 minDz = deltaZ;
                 thePrimaryV = Vertex(*itv);
               }
             }
           }
 
           Vertex theOriginalPV = thePrimaryV;
 
           muonLess.clear();
           muonLess.reserve(thePrimaryV.tracksSize());
           if (addMuonlessPrimaryVertex_ && thePrimaryV.tracksSize() > 2) {
             // Primary vertex matched to the dimuon, now refit it removing the two muons
             OniaVtxReProducer revertex(priVtxs, iEvent);
             edm::Handle<reco::TrackCollection> pvtracks;
             iEvent.getByToken(revtxtrks_, pvtracks);
             if (!pvtracks.isValid()) {
               std::cout << "pvtracks NOT valid " << std::endl;
             } else {
               edm::Handle<reco::BeamSpot> pvbeamspot;
               iEvent.getByToken(revtxbs_, pvbeamspot);
               if (pvbeamspot.id() != theBeamSpot.id())
                 edm::LogWarning("Inconsistency")
                     << "The BeamSpot used for PV reco is not the same used in this analyzer.";
               // I need to go back to the reco::Muon object, as the TrackRef in the pat::Muon can be an embedded ref.
               const reco::Muon *rmu1 = dynamic_cast<const reco::Muon *>(it->originalObject());
               const reco::Muon *rmu2 = dynamic_cast<const reco::Muon *>(it2->originalObject());
               // check that muons are truly from reco::Muons (and not, e.g., from PF Muons)
               // also check that the tracks really come from the track collection used for the BS
               if (rmu1 != nullptr && rmu2 != nullptr && rmu1->track().id() == pvtracks.id() &&
                   rmu2->track().id() == pvtracks.id()) {
                 // Save the keys of the tracks in the primary vertex
                 // std::vector<size_t> vertexTracksKeys;
                 // vertexTracksKeys.reserve(thePrimaryV.tracksSize());
                 if (thePrimaryV.hasRefittedTracks()) {
                   // Need to go back to the original tracks before taking the key
                   std::vector<reco::Track>::const_iterator itRefittedTrack = thePrimaryV.refittedTracks().begin();
                   std::vector<reco::Track>::const_iterator refittedTracksEnd = thePrimaryV.refittedTracks().end();
                   for (; itRefittedTrack != refittedTracksEnd; ++itRefittedTrack) {
                     if (thePrimaryV.originalTrack(*itRefittedTrack).key() == rmu1->track().key())
                       continue;
                     if (thePrimaryV.originalTrack(*itRefittedTrack).key() == rmu2->track().key())
                       continue;
                     // vertexTracksKeys.push_back(thePrimaryV.originalTrack(*itRefittedTrack).key());
                     muonLess.push_back(*(thePrimaryV.originalTrack(*itRefittedTrack)));
                   }
                 } else {
                   std::vector<reco::TrackBaseRef>::const_iterator itPVtrack = thePrimaryV.tracks_begin();
                   for (; itPVtrack != thePrimaryV.tracks_end(); ++itPVtrack)
                     if (itPVtrack->isNonnull()) {
                       if (itPVtrack->key() == rmu1->track().key())
                         continue;
                       if (itPVtrack->key() == rmu2->track().key())
                         continue;
                       // vertexTracksKeys.push_back(itPVtrack->key());
                       muonLess.push_back(**itPVtrack);
                     }
                 }
                 if (muonLess.size() > 1 && muonLess.size() < thePrimaryV.tracksSize()) {
                   pvs = revertex.makeVertices(muonLess, *pvbeamspot, *theTTBuilder);
                   if (!pvs.empty()) {
                     Vertex muonLessPV = Vertex(pvs.front());
                     thePrimaryV = muonLessPV;
                   }
                 }
               }
             }
           }
 
           // count the number of high Purity tracks with pT > 900 MeV attached to the chosen vertex
           double vertexWeight = -1., sumPTPV = -1.;
           int countTksOfPV = -1;
           const reco::Muon *rmu1 = dynamic_cast<const reco::Muon *>(it->originalObject());
           const reco::Muon *rmu2 = dynamic_cast<const reco::Muon *>(it2->originalObject());
           try {
             for (reco::Vertex::trackRef_iterator itVtx = theOriginalPV.tracks_begin();
                  itVtx != theOriginalPV.tracks_end();
                  itVtx++)
               if (itVtx->isNonnull()) {
                 const reco::Track &track = **itVtx;
                 if (!track.quality(reco::TrackBase::highPurity))
                   continue;
                 if (track.pt() < 0.5)
                   continue;  //reject all rejects from counting if less than 900 MeV
                 TransientTrack tt = theTTBuilder->build(track);
                 pair<bool, Measurement1D> tkPVdist = IPTools::absoluteImpactParameter3D(tt, thePrimaryV);
                 if (!tkPVdist.first)
                   continue;
                 if (tkPVdist.second.significance() > 3)
                   continue;
                 if (track.ptError() / track.pt() > 0.1)
                   continue;
                 // do not count the two muons
                 if (rmu1 != nullptr && rmu1->innerTrack().key() == itVtx->key())
                   continue;
                 if (rmu2 != nullptr && rmu2->innerTrack().key() == itVtx->key())
                   continue;
 
                 vertexWeight += theOriginalPV.trackWeight(*itVtx);
                 if (theOriginalPV.trackWeight(*itVtx) > 0.5) {
                   countTksOfPV++;
                   sumPTPV += track.pt();
                 }
               }
           } catch (std::exception &err) {
             std::cout << " muon Selection%G�%@failed " << std::endl;
             return;
           }
 
           myCand.addUserInt("countTksOfPV", countTksOfPV);
           myCand.addUserFloat("vertexWeight", (float)vertexWeight);
           myCand.addUserFloat("sumPTPV", (float)sumPTPV);
 
           TrajectoryStateClosestToPoint mu1TS = t_tks[0].impactPointTSCP();
           TrajectoryStateClosestToPoint mu2TS = t_tks[1].impactPointTSCP();
           float dca = 1E20;
           if (mu1TS.isValid() && mu2TS.isValid()) {
             ClosestApproachInRPhi cApp;
             cApp.calculate(mu1TS.theState(), mu2TS.theState());
             if (cApp.status())
               dca = cApp.distance();
           }
           myCand.addUserFloat("DCA", dca);
 
           if (addMuonlessPrimaryVertex_)
             myCand.addUserData("muonlessPV", Vertex(thePrimaryV));
           else
             myCand.addUserData("PVwithmuons", thePrimaryV);
 
           // lifetime using PV
           pvtx.SetXYZ(thePrimaryV.position().x(), thePrimaryV.position().y(), 0);
           TVector3 vdiff = vtx - pvtx;
           double cosAlpha = vdiff.Dot(pperp) / (vdiff.Perp() * pperp.Perp());
           Measurement1D distXY = vdistXY.distance(Vertex(myVertex), thePrimaryV);
           //double ctauPV = distXY.value()*cosAlpha*3.09688/pperp.Perp();
           double ctauPV = distXY.value() * cosAlpha * myCand.mass() / pperp.Perp();
           GlobalError v1e = (Vertex(myVertex)).error();
           GlobalError v2e = thePrimaryV.error();
           AlgebraicSymMatrix33 vXYe = v1e.matrix() + v2e.matrix();
           //double ctauErrPV = sqrt(vXYe.similarity(vpperp))*3.09688/(pperp.Perp2());
           double ctauErrPV = sqrt(ROOT::Math::Similarity(vpperp, vXYe)) * myCand.mass() / (pperp.Perp2());
 
           myCand.addUserFloat("ppdlPV", ctauPV);
           myCand.addUserFloat("ppdlErrPV", ctauErrPV);
           myCand.addUserFloat("cosAlpha", cosAlpha);
 
           // lifetime using BS
           pvtx.SetXYZ(theBeamSpotV.position().x(), theBeamSpotV.position().y(), 0);
           vdiff = vtx - pvtx;
           cosAlpha = vdiff.Dot(pperp) / (vdiff.Perp() * pperp.Perp());
           distXY = vdistXY.distance(Vertex(myVertex), theBeamSpotV);
           //double ctauBS = distXY.value()*cosAlpha*3.09688/pperp.Perp();
           double ctauBS = distXY.value() * cosAlpha * myCand.mass() / pperp.Perp();
           GlobalError v1eB = (Vertex(myVertex)).error();
           GlobalError v2eB = theBeamSpotV.error();
           AlgebraicSymMatrix33 vXYeB = v1eB.matrix() + v2eB.matrix();
           //double ctauErrBS = sqrt(vXYeB.similarity(vpperp))*3.09688/(pperp.Perp2());
           double ctauErrBS = sqrt(ROOT::Math::Similarity(vpperp, vXYeB)) * myCand.mass() / (pperp.Perp2());
 
           myCand.addUserFloat("ppdlBS", ctauBS);
           myCand.addUserFloat("ppdlErrBS", ctauErrBS);
 
           if (addCommonVertex_) {
             myCand.addUserData("commonVertex", Vertex(myVertex));
           }
         } else {
           myCand.addUserFloat("vNChi2", -1);
           myCand.addUserFloat("vProb", -1);
           myCand.addUserFloat("ppdlPV", -100);
           myCand.addUserFloat("ppdlErrPV", -100);
           myCand.addUserFloat("cosAlpha", -100);
           myCand.addUserFloat("ppdlBS", -100);
           myCand.addUserFloat("ppdlErrBS", -100);
           myCand.addUserFloat("DCA", -1);
           if (addCommonVertex_) {
             myCand.addUserData("commonVertex", Vertex());
           }
           if (addMuonlessPrimaryVertex_) {
             myCand.addUserData("muonlessPV", Vertex());
           } else {
             myCand.addUserData("PVwithmuons", Vertex());
           }
         }
       }
 
       // ---- MC Truth, if enabled ----
       if (addMCTruth_) {
         reco::GenParticleRef genMu1 = it->genParticleRef();
         reco::GenParticleRef genMu2 = it2->genParticleRef();
         if (genMu1.isNonnull() && genMu2.isNonnull()) {
           if (genMu1->numberOfMothers() > 0 && genMu2->numberOfMothers() > 0) {
             reco::GenParticleRef mom1 = genMu1->motherRef();
             reco::GenParticleRef mom2 = genMu2->motherRef();
             if (mom1.isNonnull() && (mom1 == mom2)) {
               myCand.setGenParticleRef(mom1);  // set
               myCand.embedGenParticle();       // and embed
               std::pair<int, float> MCinfo = findJpsiMCInfo(mom1);
               myCand.addUserInt("momPDGId", MCinfo.first);
               myCand.addUserFloat("ppdlTrue", MCinfo.second);
             } else {
               myCand.addUserInt("momPDGId", 0);
               myCand.addUserFloat("ppdlTrue", -99.);
             }
           } else {
             edm::Handle<reco::GenParticleCollection> theGenParticles;
             edm::EDGetTokenT<reco::GenParticleCollection> genCands_ =
                 consumes<reco::GenParticleCollection>((edm::InputTag) "genParticles");
             iEvent.getByToken(genCands_, theGenParticles);
             if (theGenParticles.isValid()) {
               for (size_t iGenParticle = 0; iGenParticle < theGenParticles->size(); ++iGenParticle) {
                 const Candidate &genCand = (*theGenParticles)[iGenParticle];
                 if (genCand.pdgId() == 443 || genCand.pdgId() == 100443 || genCand.pdgId() == 553 ||
                     genCand.pdgId() == 100553 || genCand.pdgId() == 200553) {
                   reco::GenParticleRef mom1(theGenParticles, iGenParticle);
                   myCand.setGenParticleRef(mom1);
                   myCand.embedGenParticle();
                   std::pair<int, float> MCinfo = findJpsiMCInfo(mom1);
                   myCand.addUserInt("momPDGId", MCinfo.first);
                   myCand.addUserFloat("ppdlTrue", MCinfo.second);
                 }
               }
             } else {
               myCand.addUserInt("momPDGId", 0);
               myCand.addUserFloat("ppdlTrue", -99.);
             }
           }
         } else {
           myCand.addUserInt("momPDGId", 0);
           myCand.addUserFloat("ppdlTrue", -99.);
         }
       }
 
       // ---- Push back output ----
       oniaOutput->push_back(myCand);
     }
   }
 
   std::sort(oniaOutput->begin(), oniaOutput->end(), vPComparator_);
 
   iEvent.put(std::move(oniaOutput));
 }

Member Data Documentation

◆ addCommonVertex_

bool Onia2MuMuPAT::addCommonVertex_

private

Definition at line 63 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ addMCTruth_

bool Onia2MuMuPAT::addMCTruth_

private

Definition at line 65 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ addMuonlessPrimaryVertex_

bool Onia2MuMuPAT::addMuonlessPrimaryVertex_

private

Definition at line 63 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ dimuonSelection_

StringCutObjectSelector<reco::Candidate, true> Onia2MuMuPAT::dimuonSelection_

private

Definition at line 62 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ higherPuritySelection_

StringCutObjectSelector<pat::Muon> Onia2MuMuPAT::higherPuritySelection_

private

Definition at line 60 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ lowerPuritySelection_

StringCutObjectSelector<pat::Muon> Onia2MuMuPAT::lowerPuritySelection_

private

Definition at line 61 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ magneticFieldToken_

edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> Onia2MuMuPAT::magneticFieldToken_

private

Definition at line 58 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ massCalculator

InvariantMassFromVertex Onia2MuMuPAT::massCalculator

private

Definition at line 68 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ muons_

edm::EDGetTokenT<edm::View<pat::Muon> > Onia2MuMuPAT::muons_

private

Definition at line 53 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ resolveAmbiguity_

bool Onia2MuMuPAT::resolveAmbiguity_

private

Definition at line 64 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ revtxbs_

edm::EDGetTokenT<reco::BeamSpot> Onia2MuMuPAT::revtxbs_

private

Definition at line 57 of file Onia2MuMuPAT.h.

Referenced by Onia2MuMuPAT(), and produce().

◆ revtxtrks_

edm::EDGetTokenT<reco::TrackCollection> Onia2MuMuPAT::revtxtrks_

private

Definition at line 56 of file Onia2MuMuPAT.h.

Referenced by Onia2MuMuPAT(), and produce().

◆ thebeamspot_

edm::EDGetTokenT<reco::BeamSpot> Onia2MuMuPAT::thebeamspot_

private

Definition at line 54 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ thePVs_

edm::EDGetTokenT<reco::VertexCollection> Onia2MuMuPAT::thePVs_

private

Definition at line 55 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ theTTBuilderToken_

edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> Onia2MuMuPAT::theTTBuilderToken_

private

Definition at line 59 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ vPComparator_

GreaterByVProb<pat::CompositeCandidate> Onia2MuMuPAT::vPComparator_

private

Definition at line 66 of file Onia2MuMuPAT.h.

Referenced by produce().

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ Onia2MuMuPAT()

Member Function Documentation

◆ fillDescriptions()

◆ findJpsiMCInfo()

◆ isAbHadron()

◆ isAMixedbHadron()

◆ produce()

Member Data Documentation

◆ addCommonVertex_

◆ addMCTruth_

◆ addMuonlessPrimaryVertex_

◆ dimuonSelection_

◆ higherPuritySelection_

◆ lowerPuritySelection_

◆ magneticFieldToken_

◆ massCalculator

◆ muons_

◆ resolveAmbiguity_

◆ revtxbs_

◆ revtxtrks_

◆ thebeamspot_

◆ thePVs_

◆ theTTBuilderToken_

◆ vPComparator_