#include <Onia2MuMuPAT.h>

Inheritance diagram for Onia2MuMuPAT:

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

	~Onia2MuMuPAT () override

Public Member Functions inherited from edm::EDProducer
	EDProducer ()

SerialTaskQueue *	globalLuminosityBlocksQueue ()

SerialTaskQueue *	globalRunsQueue ()

ModuleDescription const &	moduleDescription () const

	~EDProducer () override

Public Member Functions inherited from edm::ProducerBase
void	callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)

std::vector< edm::ProductResolverIndex > const &	indiciesForPutProducts (BranchType iBranchType) const

	ProducerBase ()

std::vector< edm::ProductResolverIndex > const &	putTokenIndexToProductResolverIndex () const

std::vector< bool > const &	recordProvenanceList () const

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

std::function< void(BranchDescription const &)>	registrationCallback () const
	used by the fwk to register list of products More...

void	resolvePutIndicies (BranchType iBranchType, ModuleToResolverIndicies const &iIndicies, std::string const &moduleLabel)

TypeLabelList const &	typeLabelList () const
	used by the fwk to register the list of products of this module More...

	~ProducerBase () noexcept(false) override

Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfo >	consumesInfo () const

void	convertCurrentProcessAlias (std::string const &processName)
	Convert "@currentProcess" in InputTag process names to the actual current process name. More...

	EDConsumerBase ()

	EDConsumerBase (EDConsumerBase &&)=default

	EDConsumerBase (EDConsumerBase const &)=delete

ESProxyIndex const *	esGetTokenIndices (edm::Transition iTrans) const

ProductResolverIndexAndSkipBit	indexFrom (EDGetToken, BranchType, TypeID const &) const

void	itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

void	itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const

std::vector< ProductResolverIndexAndSkipBit > const &	itemsToGetFrom (BranchType iType) const

void	labelsForToken (EDGetToken iToken, Labels &oLabels) const

void	modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const

EDConsumerBase &	operator= (EDConsumerBase &&)=default

EDConsumerBase const &	operator= (EDConsumerBase const &)=delete

bool	registeredToConsume (ProductResolverIndex, bool, BranchType) const

bool	registeredToConsumeMany (TypeID const &, BranchType) const

ProductResolverIndexAndSkipBit	uncheckedIndexFrom (EDGetToken) const

void	updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)

void	updateLookup (eventsetup::ESRecordsToProxyIndices const &)

virtual	~EDConsumerBase () noexcept(false)

Private Member Functions
void	beginJob () override

void	endJob () override

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

bool	isAbHadron (int pdgID)

bool	isAMixedbHadron (int pdgID, int momPdgID)

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

Additional Inherited Members
Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

Public Types inherited from edm::ProducerBase
using	ModuleToResolverIndicies = std::unordered_multimap< std::string, std::tuple< edm::TypeID const , const char , edm::ProductResolverIndex > >

typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Public Types inherited from edm::EDConsumerBase
typedef ProductLabels	Labels

Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

static bool	wantsGlobalLuminosityBlocks ()

static bool	wantsGlobalRuns ()

static bool	wantsStreamLuminosityBlocks ()

static bool	wantsStreamRuns ()

Protected Member Functions inherited from edm::ProducerBase
template<class ProductType >
BranchAliasSetterT< ProductType >	produces ()
	declare what type of product will make and with which optional label More...

template<typename ProductType , BranchType B>
BranchAliasSetterT< ProductType >	produces ()

template<typename ProductType , Transition B>
BranchAliasSetterT< ProductType >	produces ()

BranchAliasSetter	produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)

template<BranchType B>
BranchAliasSetter	produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)

template<Transition B>
BranchAliasSetter	produces (const TypeID &id, std::string instanceName=std::string(), bool recordProvenance=true)

template<typename ProductType , Transition B>
BranchAliasSetterT< ProductType >	produces (std::string instanceName)

template<class ProductType >
BranchAliasSetterT< ProductType >	produces (std::string instanceName)

template<typename ProductType , BranchType B>
BranchAliasSetterT< ProductType >	produces (std::string instanceName)

ProducesCollector	producesCollector ()

Protected Member Functions inherited from edm::EDConsumerBase
EDGetToken	consumes (const TypeToGet &id, edm::InputTag const &tag)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	consumes (edm::InputTag const &tag)

template<BranchType B>
EDGetToken	consumes (TypeToGet const &id, edm::InputTag const &tag)

ConsumesCollector	consumesCollector ()
	Use a ConsumesCollector to gather consumes information from helper functions. More...

template<typename ProductType , BranchType B = InEvent>
void	consumesMany ()

void	consumesMany (const TypeToGet &id)

template<BranchType B>
void	consumesMany (const TypeToGet &id)

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes ()

template<typename ESProduct , typename ESRecord , Transition Tr = Transition::Event>
auto	esConsumes (ESInputTag const &tag)

EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<BranchType B>
EDGetToken	mayConsume (const TypeToGet &id, edm::InputTag const &tag)

template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType >	mayConsume (edm::InputTag const &tag)

Detailed Description

Definition at line 33 of file Onia2MuMuPAT.h.

Constructor & Destructor Documentation

◆ Onia2MuMuPAT()

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

explicit

Definition at line 29 of file Onia2MuMuPAT.cc.

     : 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"))),
       higherPuritySelection_(iConfig.getParameter<std::string>("higherPuritySelection")),
       lowerPuritySelection_(iConfig.getParameter<std::string>("lowerPuritySelection")),
       dimuonSelection_(
           iConfig.existsAs<std::string>("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>();
 }

References revtxbs_, and revtxtrks_.

◆ ~Onia2MuMuPAT()

Onia2MuMuPAT::~Onia2MuMuPAT ( )

override

Definition at line 47 of file Onia2MuMuPAT.cc.

                             {
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)
 }

Member Function Documentation

◆ beginJob()

void Onia2MuMuPAT::beginJob ( void )

overrideprivatevirtual

Reimplemented from edm::EDProducer.

Definition at line 496 of file Onia2MuMuPAT.cc.

496 {}

◆ endJob()

void Onia2MuMuPAT::endJob ( void )

overrideprivatevirtual

Reimplemented from edm::EDProducer.

Definition at line 499 of file Onia2MuMuPAT.cc.

499 {}

◆ findJpsiMCInfo()

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

private

Definition at line 439 of file Onia2MuMuPAT.cc.

                                                                            {
   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;
 }

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

Referenced by produce().

◆ isAbHadron()

bool Onia2MuMuPAT::isAbHadron ( int pdgID )

private

Definition at line 426 of file Onia2MuMuPAT.cc.

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

References funct::abs().

Referenced by findJpsiMCInfo().

◆ isAMixedbHadron()

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

private

Definition at line 432 of file Onia2MuMuPAT.cc.

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

References funct::abs().

Referenced by findJpsiMCInfo().

◆ produce()

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

overrideprivatevirtual

DCA

end DCA

Implements edm::EDProducer.

Definition at line 57 of file Onia2MuMuPAT.cc.

                                                                         {
   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;
  
   ESHandle<MagneticField> magneticField;
   iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
   const MagneticField *field = magneticField.product();
  
   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.get<TransientTrackRecord>().get("TransientTrackBuilder", theTTBuilder);
   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),
                     &(*magneticField)),
                 GlobalTrajectoryParameters(GlobalPoint(bs.position().x(), bs.position().y(), bs.position().z()),
                                            GlobalVector(bs.dxdz(), bs.dydz(), 1.),
                                            TrackCharge(0),
                                            &(*magneticField)));
             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, iSetup);
                   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));
 }

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, dimuonSelection_, VertexDistanceXY::distance(), ClosestApproachInRPhi::distance(), pat::PATObject< ObjectType >::embedGenParticle(), runTheMatrix::err, Measurement1D::error(), relativeConstraints::error, reco::Vertex::error(), cppFunctionSkipper::exception, findJpsiMCInfo(), edm::EventSetup::get(), get, reco::Vertex::hasRefittedTracks(), higherPuritySelection_, reco::TrackBase::highPurity, edm::ProductID::id(), edm::HandleBase::id(), edm::Ref< C, T, F >::id(), iEvent, reco::Muon::innerTrack(), InvariantMassFromVertex::invariantMass(), edm::Ref< C, T, F >::isNonnull(), edm::HandleBase::isValid(), TrajectoryStateClosestToPoint::isValid(), TransientVertex::isValid(), jpsi, edm::RefToBase< T >::key(), edm::Ref< C, T, F >::key(), lowerPuritySelection_, HLT_2018_cff::magneticField, OniaVtxReProducer::makeVertices(), reco::LeafCandidate::mass(), massCalculator, GlobalErrorBase< T, ErrorWeightType >::matrix(), pfMETCorrectionType0_cfi::minDz, eostools::move(), PDWG_BPHSkim_cff::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(), resolveAmbiguity_, revtxbs_, revtxtrks_, reco::LeafCandidate::setCharge(), pat::PATObject< ObjectType >::setGenParticleRef(), reco::LeafCandidate::setP4(), mathSSE::sqrt(), ClosestApproachInRPhi::status(), mps_update::status, thebeamspot_, thePVs_, TrajectoryStateClosestToPoint::theState(), reco::Muon::track(), HLT_2018_cff::track, reco::Vertex::tracks_begin(), reco::Vertex::tracks_end(), reco::Vertex::tracksSize(), reco::Vertex::trackWeight(), groupFilesInBlocks::tt, Measurement1D::value(), HltBtagValidation_cff::Vertex, KalmanVertexFitter::vertex(), vPComparator_, badGlobalMuonTaggersAOD_cff::vtx, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Member Data Documentation

◆ addCommonVertex_

bool Onia2MuMuPAT::addCommonVertex_

private

Definition at line 56 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ addMCTruth_

bool Onia2MuMuPAT::addMCTruth_

private

Definition at line 58 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ addMuonlessPrimaryVertex_

bool Onia2MuMuPAT::addMuonlessPrimaryVertex_

private

Definition at line 56 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ dimuonSelection_

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

private

Definition at line 55 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ higherPuritySelection_

StringCutObjectSelector<pat::Muon> Onia2MuMuPAT::higherPuritySelection_

private

Definition at line 53 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ lowerPuritySelection_

StringCutObjectSelector<pat::Muon> Onia2MuMuPAT::lowerPuritySelection_

private

Definition at line 54 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ massCalculator

InvariantMassFromVertex Onia2MuMuPAT::massCalculator

private

Definition at line 61 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ muons_

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

private

Definition at line 48 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ resolveAmbiguity_

bool Onia2MuMuPAT::resolveAmbiguity_

private

Definition at line 57 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ revtxbs_

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

private

Definition at line 52 of file Onia2MuMuPAT.h.

Referenced by Onia2MuMuPAT(), and produce().

◆ revtxtrks_

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

private

Definition at line 51 of file Onia2MuMuPAT.h.

Referenced by Onia2MuMuPAT(), and produce().

◆ thebeamspot_

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

private

Definition at line 49 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ thePVs_

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

private

Definition at line 50 of file Onia2MuMuPAT.h.

Referenced by produce().

◆ vPComparator_

GreaterByVProb<pat::CompositeCandidate> Onia2MuMuPAT::vPComparator_

private

Definition at line 59 of file Onia2MuMuPAT.h.

Referenced by produce().

Private Attributes
bool	addCommonVertex_

bool	addMCTruth_

bool	addMuonlessPrimaryVertex_

StringCutObjectSelector< reco::Candidate, true >	dimuonSelection_

StringCutObjectSelector< pat::Muon >	higherPuritySelection_

StringCutObjectSelector< pat::Muon >	lowerPuritySelection_

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_

GreaterByVProb< pat::CompositeCandidate >	vPComparator_

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ Onia2MuMuPAT()

◆ ~Onia2MuMuPAT()

Member Function Documentation

◆ beginJob()

◆ endJob()

◆ findJpsiMCInfo()

◆ isAbHadron()

◆ isAMixedbHadron()

◆ produce()

Member Data Documentation

◆ addCommonVertex_

◆ addMCTruth_

◆ addMuonlessPrimaryVertex_

◆ dimuonSelection_

◆ higherPuritySelection_

◆ lowerPuritySelection_

◆ massCalculator

◆ muons_

◆ resolveAmbiguity_

◆ revtxbs_

◆ revtxtrks_

◆ thebeamspot_

◆ thePVs_

◆ vPComparator_