Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:15:33

0001 #include "HeavyFlavorAnalysis/Onia2MuMu/interface/Onia2MuMuPAT.h"
0002 
0003 //Headers for the data items
0004 #include <DataFormats/TrackReco/interface/TrackFwd.h>
0005 #include <DataFormats/TrackReco/interface/Track.h>
0006 #include <DataFormats/MuonReco/interface/MuonFwd.h>
0007 #include <DataFormats/MuonReco/interface/Muon.h>
0008 #include <DataFormats/Common/interface/View.h>
0009 #include <DataFormats/HepMCCandidate/interface/GenParticle.h>
0010 #include <DataFormats/PatCandidates/interface/Muon.h>
0011 #include <DataFormats/VertexReco/interface/VertexFwd.h>
0012 
0013 //Headers for services and tools
0014 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0015 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
0016 #include "TMath.h"
0017 #include "Math/VectorUtil.h"
0018 #include "TVector3.h"
0019 #include "HeavyFlavorAnalysis/Onia2MuMu/interface/OniaVtxReProducer.h"
0020 
0021 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
0022 #include "TrackingTools/IPTools/interface/IPTools.h"
0023 #include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
0024 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0025 
0026 Onia2MuMuPAT::Onia2MuMuPAT(const edm::ParameterSet &iConfig)
0027     : muons_(consumes<edm::View<pat::Muon>>(iConfig.getParameter<edm::InputTag>("muons"))),
0028       thebeamspot_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpotTag"))),
0029       thePVs_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("primaryVertexTag"))),
0030       magneticFieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0031       theTTBuilderToken_(
0032           esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))),
0033       higherPuritySelection_(iConfig.getParameter<std::string>("higherPuritySelection")),
0034       lowerPuritySelection_(iConfig.getParameter<std::string>("lowerPuritySelection")),
0035       dimuonSelection_(iConfig.getParameter<std::string>("dimuonSelection")),
0036       addCommonVertex_(iConfig.getParameter<bool>("addCommonVertex")),
0037       addMuonlessPrimaryVertex_(iConfig.getParameter<bool>("addMuonlessPrimaryVertex")),
0038       resolveAmbiguity_(iConfig.getParameter<bool>("resolvePileUpAmbiguity")),
0039       addMCTruth_(iConfig.getParameter<bool>("addMCTruth")) {
0040   revtxtrks_ = consumes<reco::TrackCollection>(
0041       (edm::InputTag) "generalTracks");  //if that is not true, we will raise an exception
0042   revtxbs_ = consumes<reco::BeamSpot>((edm::InputTag) "offlineBeamSpot");
0043   produces<pat::CompositeCandidateCollection>();
0044 }
0045 
0046 //
0047 // member functions
0048 //
0049 
0050 // ------------ method called to produce the data  ------------
0051 void Onia2MuMuPAT::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0052   using namespace edm;
0053   using namespace std;
0054   using namespace reco;
0055   typedef Candidate::LorentzVector LorentzVector;
0056 
0057   vector<double> muMasses;
0058   muMasses.push_back(0.1056583715);
0059   muMasses.push_back(0.1056583715);
0060 
0061   std::unique_ptr<pat::CompositeCandidateCollection> oniaOutput(new pat::CompositeCandidateCollection);
0062 
0063   Vertex thePrimaryV;
0064   Vertex theBeamSpotV;
0065 
0066   const MagneticField &field = iSetup.getData(magneticFieldToken_);
0067 
0068   Handle<BeamSpot> theBeamSpot;
0069   iEvent.getByToken(thebeamspot_, theBeamSpot);
0070   BeamSpot bs = *theBeamSpot;
0071   theBeamSpotV = Vertex(bs.position(), bs.covariance3D());
0072 
0073   Handle<VertexCollection> priVtxs;
0074   iEvent.getByToken(thePVs_, priVtxs);
0075   if (priVtxs->begin() != priVtxs->end()) {
0076     thePrimaryV = Vertex(*(priVtxs->begin()));
0077   } else {
0078     thePrimaryV = Vertex(bs.position(), bs.covariance3D());
0079   }
0080 
0081   Handle<View<pat::Muon>> muons;
0082   iEvent.getByToken(muons_, muons);
0083 
0084   edm::ESHandle<TransientTrackBuilder> theTTBuilder = iSetup.getHandle(theTTBuilderToken_);
0085   KalmanVertexFitter vtxFitter(true);
0086   TrackCollection muonLess;
0087 
0088   // JPsi candidates only from muons
0089   for (View<pat::Muon>::const_iterator it = muons->begin(), itend = muons->end(); it != itend; ++it) {
0090     // both must pass low quality
0091     if (!lowerPuritySelection_(*it))
0092       continue;
0093     for (View<pat::Muon>::const_iterator it2 = it + 1; it2 != itend; ++it2) {
0094       // both must pass low quality
0095       if (!lowerPuritySelection_(*it2))
0096         continue;
0097       // one must pass tight quality
0098       if (!(higherPuritySelection_(*it) || higherPuritySelection_(*it2)))
0099         continue;
0100 
0101       pat::CompositeCandidate myCand;
0102       vector<TransientVertex> pvs;
0103 
0104       // ---- no explicit order defined ----
0105       myCand.addDaughter(*it, "muon1");
0106       myCand.addDaughter(*it2, "muon2");
0107 
0108       // ---- define and set candidate's 4momentum  ----
0109       LorentzVector jpsi = it->p4() + it2->p4();
0110       myCand.setP4(jpsi);
0111       myCand.setCharge(it->charge() + it2->charge());
0112 
0113       // ---- apply the dimuon cut ----
0114       if (!dimuonSelection_(myCand))
0115         continue;
0116 
0117       // ---- fit vertex using Tracker tracks (if they have tracks) ----
0118       if (it->track().isNonnull() && it2->track().isNonnull()) {
0119         //build the dimuon secondary vertex
0120         vector<TransientTrack> t_tks;
0121         t_tks.push_back(theTTBuilder->build(
0122             *it->track()));  // pass the reco::Track, not  the reco::TrackRef (which can be transient)
0123         t_tks.push_back(theTTBuilder->build(*it2->track()));  // otherwise the vertex will have transient refs inside.
0124         TransientVertex myVertex = vtxFitter.vertex(t_tks);
0125 
0126         CachingVertex<5> VtxForInvMass = vtxFitter.vertex(t_tks);
0127 
0128         Measurement1D MassWErr(jpsi.M(), -9999.);
0129         if (field.nominalValue() > 0) {
0130           MassWErr = massCalculator.invariantMass(VtxForInvMass, muMasses);
0131         } else {
0132           myVertex = TransientVertex();  // with no arguments it is invalid
0133         }
0134 
0135         myCand.addUserFloat("MassErr", MassWErr.error());
0136 
0137         if (myVertex.isValid()) {
0138           float vChi2 = myVertex.totalChiSquared();
0139           float vNDF = myVertex.degreesOfFreedom();
0140           float vProb(TMath::Prob(vChi2, (int)vNDF));
0141 
0142           myCand.addUserFloat("vNChi2", vChi2 / vNDF);
0143           myCand.addUserFloat("vProb", vProb);
0144 
0145           TVector3 vtx;
0146           TVector3 pvtx;
0147           VertexDistanceXY vdistXY;
0148 
0149           vtx.SetXYZ(myVertex.position().x(), myVertex.position().y(), 0);
0150           TVector3 pperp(jpsi.px(), jpsi.py(), 0);
0151           AlgebraicVector3 vpperp(pperp.x(), pperp.y(), 0);
0152 
0153           if (resolveAmbiguity_) {
0154             float minDz = 999999.;
0155             TwoTrackMinimumDistance ttmd;
0156             bool status = ttmd.calculate(
0157                 GlobalTrajectoryParameters(
0158                     GlobalPoint(myVertex.position().x(), myVertex.position().y(), myVertex.position().z()),
0159                     GlobalVector(myCand.px(), myCand.py(), myCand.pz()),
0160                     TrackCharge(0),
0161                     &(field)),
0162                 GlobalTrajectoryParameters(GlobalPoint(bs.position().x(), bs.position().y(), bs.position().z()),
0163                                            GlobalVector(bs.dxdz(), bs.dydz(), 1.),
0164                                            TrackCharge(0),
0165                                            &(field)));
0166             float extrapZ = -9E20;
0167             if (status)
0168               extrapZ = ttmd.points().first.z();
0169 
0170             for (VertexCollection::const_iterator itv = priVtxs->begin(), itvend = priVtxs->end(); itv != itvend;
0171                  ++itv) {
0172               float deltaZ = fabs(extrapZ - itv->position().z());
0173               if (deltaZ < minDz) {
0174                 minDz = deltaZ;
0175                 thePrimaryV = Vertex(*itv);
0176               }
0177             }
0178           }
0179 
0180           Vertex theOriginalPV = thePrimaryV;
0181 
0182           muonLess.clear();
0183           muonLess.reserve(thePrimaryV.tracksSize());
0184           if (addMuonlessPrimaryVertex_ && thePrimaryV.tracksSize() > 2) {
0185             // Primary vertex matched to the dimuon, now refit it removing the two muons
0186             OniaVtxReProducer revertex(priVtxs, iEvent);
0187             edm::Handle<reco::TrackCollection> pvtracks;
0188             iEvent.getByToken(revtxtrks_, pvtracks);
0189             if (!pvtracks.isValid()) {
0190               std::cout << "pvtracks NOT valid " << std::endl;
0191             } else {
0192               edm::Handle<reco::BeamSpot> pvbeamspot;
0193               iEvent.getByToken(revtxbs_, pvbeamspot);
0194               if (pvbeamspot.id() != theBeamSpot.id())
0195                 edm::LogWarning("Inconsistency")
0196                     << "The BeamSpot used for PV reco is not the same used in this analyzer.";
0197               // I need to go back to the reco::Muon object, as the TrackRef in the pat::Muon can be an embedded ref.
0198               const reco::Muon *rmu1 = dynamic_cast<const reco::Muon *>(it->originalObject());
0199               const reco::Muon *rmu2 = dynamic_cast<const reco::Muon *>(it2->originalObject());
0200               // check that muons are truly from reco::Muons (and not, e.g., from PF Muons)
0201               // also check that the tracks really come from the track collection used for the BS
0202               if (rmu1 != nullptr && rmu2 != nullptr && rmu1->track().id() == pvtracks.id() &&
0203                   rmu2->track().id() == pvtracks.id()) {
0204                 // Save the keys of the tracks in the primary vertex
0205                 // std::vector<size_t> vertexTracksKeys;
0206                 // vertexTracksKeys.reserve(thePrimaryV.tracksSize());
0207                 if (thePrimaryV.hasRefittedTracks()) {
0208                   // Need to go back to the original tracks before taking the key
0209                   std::vector<reco::Track>::const_iterator itRefittedTrack = thePrimaryV.refittedTracks().begin();
0210                   std::vector<reco::Track>::const_iterator refittedTracksEnd = thePrimaryV.refittedTracks().end();
0211                   for (; itRefittedTrack != refittedTracksEnd; ++itRefittedTrack) {
0212                     if (thePrimaryV.originalTrack(*itRefittedTrack).key() == rmu1->track().key())
0213                       continue;
0214                     if (thePrimaryV.originalTrack(*itRefittedTrack).key() == rmu2->track().key())
0215                       continue;
0216                     // vertexTracksKeys.push_back(thePrimaryV.originalTrack(*itRefittedTrack).key());
0217                     muonLess.push_back(*(thePrimaryV.originalTrack(*itRefittedTrack)));
0218                   }
0219                 } else {
0220                   std::vector<reco::TrackBaseRef>::const_iterator itPVtrack = thePrimaryV.tracks_begin();
0221                   for (; itPVtrack != thePrimaryV.tracks_end(); ++itPVtrack)
0222                     if (itPVtrack->isNonnull()) {
0223                       if (itPVtrack->key() == rmu1->track().key())
0224                         continue;
0225                       if (itPVtrack->key() == rmu2->track().key())
0226                         continue;
0227                       // vertexTracksKeys.push_back(itPVtrack->key());
0228                       muonLess.push_back(**itPVtrack);
0229                     }
0230                 }
0231                 if (muonLess.size() > 1 && muonLess.size() < thePrimaryV.tracksSize()) {
0232                   pvs = revertex.makeVertices(muonLess, *pvbeamspot, *theTTBuilder);
0233                   if (!pvs.empty()) {
0234                     Vertex muonLessPV = Vertex(pvs.front());
0235                     thePrimaryV = muonLessPV;
0236                   }
0237                 }
0238               }
0239             }
0240           }
0241 
0242           // count the number of high Purity tracks with pT > 900 MeV attached to the chosen vertex
0243           double vertexWeight = -1., sumPTPV = -1.;
0244           int countTksOfPV = -1;
0245           const reco::Muon *rmu1 = dynamic_cast<const reco::Muon *>(it->originalObject());
0246           const reco::Muon *rmu2 = dynamic_cast<const reco::Muon *>(it2->originalObject());
0247           try {
0248             for (reco::Vertex::trackRef_iterator itVtx = theOriginalPV.tracks_begin();
0249                  itVtx != theOriginalPV.tracks_end();
0250                  itVtx++)
0251               if (itVtx->isNonnull()) {
0252                 const reco::Track &track = **itVtx;
0253                 if (!track.quality(reco::TrackBase::highPurity))
0254                   continue;
0255                 if (track.pt() < 0.5)
0256                   continue;  //reject all rejects from counting if less than 900 MeV
0257                 TransientTrack tt = theTTBuilder->build(track);
0258                 pair<bool, Measurement1D> tkPVdist = IPTools::absoluteImpactParameter3D(tt, thePrimaryV);
0259                 if (!tkPVdist.first)
0260                   continue;
0261                 if (tkPVdist.second.significance() > 3)
0262                   continue;
0263                 if (track.ptError() / track.pt() > 0.1)
0264                   continue;
0265                 // do not count the two muons
0266                 if (rmu1 != nullptr && rmu1->innerTrack().key() == itVtx->key())
0267                   continue;
0268                 if (rmu2 != nullptr && rmu2->innerTrack().key() == itVtx->key())
0269                   continue;
0270 
0271                 vertexWeight += theOriginalPV.trackWeight(*itVtx);
0272                 if (theOriginalPV.trackWeight(*itVtx) > 0.5) {
0273                   countTksOfPV++;
0274                   sumPTPV += track.pt();
0275                 }
0276               }
0277           } catch (std::exception &err) {
0278             std::cout << " muon Selection%G�%@failed " << std::endl;
0279             return;
0280           }
0281 
0282           myCand.addUserInt("countTksOfPV", countTksOfPV);
0283           myCand.addUserFloat("vertexWeight", (float)vertexWeight);
0284           myCand.addUserFloat("sumPTPV", (float)sumPTPV);
0285 
0286           ///DCA
0287           TrajectoryStateClosestToPoint mu1TS = t_tks[0].impactPointTSCP();
0288           TrajectoryStateClosestToPoint mu2TS = t_tks[1].impactPointTSCP();
0289           float dca = 1E20;
0290           if (mu1TS.isValid() && mu2TS.isValid()) {
0291             ClosestApproachInRPhi cApp;
0292             cApp.calculate(mu1TS.theState(), mu2TS.theState());
0293             if (cApp.status())
0294               dca = cApp.distance();
0295           }
0296           myCand.addUserFloat("DCA", dca);
0297           ///end DCA
0298 
0299           if (addMuonlessPrimaryVertex_)
0300             myCand.addUserData("muonlessPV", Vertex(thePrimaryV));
0301           else
0302             myCand.addUserData("PVwithmuons", thePrimaryV);
0303 
0304           // lifetime using PV
0305           pvtx.SetXYZ(thePrimaryV.position().x(), thePrimaryV.position().y(), 0);
0306           TVector3 vdiff = vtx - pvtx;
0307           double cosAlpha = vdiff.Dot(pperp) / (vdiff.Perp() * pperp.Perp());
0308           Measurement1D distXY = vdistXY.distance(Vertex(myVertex), thePrimaryV);
0309           //double ctauPV = distXY.value()*cosAlpha*3.09688/pperp.Perp();
0310           double ctauPV = distXY.value() * cosAlpha * myCand.mass() / pperp.Perp();
0311           GlobalError v1e = (Vertex(myVertex)).error();
0312           GlobalError v2e = thePrimaryV.error();
0313           AlgebraicSymMatrix33 vXYe = v1e.matrix() + v2e.matrix();
0314           //double ctauErrPV = sqrt(vXYe.similarity(vpperp))*3.09688/(pperp.Perp2());
0315           double ctauErrPV = sqrt(ROOT::Math::Similarity(vpperp, vXYe)) * myCand.mass() / (pperp.Perp2());
0316 
0317           myCand.addUserFloat("ppdlPV", ctauPV);
0318           myCand.addUserFloat("ppdlErrPV", ctauErrPV);
0319           myCand.addUserFloat("cosAlpha", cosAlpha);
0320 
0321           // lifetime using BS
0322           pvtx.SetXYZ(theBeamSpotV.position().x(), theBeamSpotV.position().y(), 0);
0323           vdiff = vtx - pvtx;
0324           cosAlpha = vdiff.Dot(pperp) / (vdiff.Perp() * pperp.Perp());
0325           distXY = vdistXY.distance(Vertex(myVertex), theBeamSpotV);
0326           //double ctauBS = distXY.value()*cosAlpha*3.09688/pperp.Perp();
0327           double ctauBS = distXY.value() * cosAlpha * myCand.mass() / pperp.Perp();
0328           GlobalError v1eB = (Vertex(myVertex)).error();
0329           GlobalError v2eB = theBeamSpotV.error();
0330           AlgebraicSymMatrix33 vXYeB = v1eB.matrix() + v2eB.matrix();
0331           //double ctauErrBS = sqrt(vXYeB.similarity(vpperp))*3.09688/(pperp.Perp2());
0332           double ctauErrBS = sqrt(ROOT::Math::Similarity(vpperp, vXYeB)) * myCand.mass() / (pperp.Perp2());
0333 
0334           myCand.addUserFloat("ppdlBS", ctauBS);
0335           myCand.addUserFloat("ppdlErrBS", ctauErrBS);
0336 
0337           if (addCommonVertex_) {
0338             myCand.addUserData("commonVertex", Vertex(myVertex));
0339           }
0340         } else {
0341           myCand.addUserFloat("vNChi2", -1);
0342           myCand.addUserFloat("vProb", -1);
0343           myCand.addUserFloat("ppdlPV", -100);
0344           myCand.addUserFloat("ppdlErrPV", -100);
0345           myCand.addUserFloat("cosAlpha", -100);
0346           myCand.addUserFloat("ppdlBS", -100);
0347           myCand.addUserFloat("ppdlErrBS", -100);
0348           myCand.addUserFloat("DCA", -1);
0349           if (addCommonVertex_) {
0350             myCand.addUserData("commonVertex", Vertex());
0351           }
0352           if (addMuonlessPrimaryVertex_) {
0353             myCand.addUserData("muonlessPV", Vertex());
0354           } else {
0355             myCand.addUserData("PVwithmuons", Vertex());
0356           }
0357         }
0358       }
0359 
0360       // ---- MC Truth, if enabled ----
0361       if (addMCTruth_) {
0362         reco::GenParticleRef genMu1 = it->genParticleRef();
0363         reco::GenParticleRef genMu2 = it2->genParticleRef();
0364         if (genMu1.isNonnull() && genMu2.isNonnull()) {
0365           if (genMu1->numberOfMothers() > 0 && genMu2->numberOfMothers() > 0) {
0366             reco::GenParticleRef mom1 = genMu1->motherRef();
0367             reco::GenParticleRef mom2 = genMu2->motherRef();
0368             if (mom1.isNonnull() && (mom1 == mom2)) {
0369               myCand.setGenParticleRef(mom1);  // set
0370               myCand.embedGenParticle();       // and embed
0371               std::pair<int, float> MCinfo = findJpsiMCInfo(mom1);
0372               myCand.addUserInt("momPDGId", MCinfo.first);
0373               myCand.addUserFloat("ppdlTrue", MCinfo.second);
0374             } else {
0375               myCand.addUserInt("momPDGId", 0);
0376               myCand.addUserFloat("ppdlTrue", -99.);
0377             }
0378           } else {
0379             edm::Handle<reco::GenParticleCollection> theGenParticles;
0380             edm::EDGetTokenT<reco::GenParticleCollection> genCands_ =
0381                 consumes<reco::GenParticleCollection>((edm::InputTag) "genParticles");
0382             iEvent.getByToken(genCands_, theGenParticles);
0383             if (theGenParticles.isValid()) {
0384               for (size_t iGenParticle = 0; iGenParticle < theGenParticles->size(); ++iGenParticle) {
0385                 const Candidate &genCand = (*theGenParticles)[iGenParticle];
0386                 if (genCand.pdgId() == 443 || genCand.pdgId() == 100443 || genCand.pdgId() == 553 ||
0387                     genCand.pdgId() == 100553 || genCand.pdgId() == 200553) {
0388                   reco::GenParticleRef mom1(theGenParticles, iGenParticle);
0389                   myCand.setGenParticleRef(mom1);
0390                   myCand.embedGenParticle();
0391                   std::pair<int, float> MCinfo = findJpsiMCInfo(mom1);
0392                   myCand.addUserInt("momPDGId", MCinfo.first);
0393                   myCand.addUserFloat("ppdlTrue", MCinfo.second);
0394                 }
0395               }
0396             } else {
0397               myCand.addUserInt("momPDGId", 0);
0398               myCand.addUserFloat("ppdlTrue", -99.);
0399             }
0400           }
0401         } else {
0402           myCand.addUserInt("momPDGId", 0);
0403           myCand.addUserFloat("ppdlTrue", -99.);
0404         }
0405       }
0406 
0407       // ---- Push back output ----
0408       oniaOutput->push_back(myCand);
0409     }
0410   }
0411 
0412   std::sort(oniaOutput->begin(), oniaOutput->end(), vPComparator_);
0413 
0414   iEvent.put(std::move(oniaOutput));
0415 }
0416 
0417 bool Onia2MuMuPAT::isAbHadron(int pdgID) const {
0418   if (abs(pdgID) == 511 || abs(pdgID) == 521 || abs(pdgID) == 531 || abs(pdgID) == 5122)
0419     return true;
0420   return false;
0421 }
0422 
0423 bool Onia2MuMuPAT::isAMixedbHadron(int pdgID, int momPdgID) const {
0424   if ((abs(pdgID) == 511 && abs(momPdgID) == 511 && pdgID * momPdgID < 0) ||
0425       (abs(pdgID) == 531 && abs(momPdgID) == 531 && pdgID * momPdgID < 0))
0426     return true;
0427   return false;
0428 }
0429 
0430 std::pair<int, float> Onia2MuMuPAT::findJpsiMCInfo(reco::GenParticleRef genJpsi) const {
0431   int momJpsiID = 0;
0432   float trueLife = -99.;
0433 
0434   if (genJpsi->numberOfMothers() > 0) {
0435     TVector3 trueVtx(0.0, 0.0, 0.0);
0436     TVector3 trueP(0.0, 0.0, 0.0);
0437     TVector3 trueVtxMom(0.0, 0.0, 0.0);
0438 
0439     trueVtx.SetXYZ(genJpsi->vertex().x(), genJpsi->vertex().y(), genJpsi->vertex().z());
0440     trueP.SetXYZ(genJpsi->momentum().x(), genJpsi->momentum().y(), genJpsi->momentum().z());
0441 
0442     bool aBhadron = false;
0443     reco::GenParticleRef Jpsimom = genJpsi->motherRef();  // find mothers
0444     if (Jpsimom.isNull()) {
0445       std::pair<int, float> result = std::make_pair(momJpsiID, trueLife);
0446       return result;
0447     } else {
0448       reco::GenParticleRef Jpsigrandmom = Jpsimom->motherRef();
0449       if (isAbHadron(Jpsimom->pdgId())) {
0450         if (Jpsigrandmom.isNonnull() && isAMixedbHadron(Jpsimom->pdgId(), Jpsigrandmom->pdgId())) {
0451           momJpsiID = Jpsigrandmom->pdgId();
0452           trueVtxMom.SetXYZ(Jpsigrandmom->vertex().x(), Jpsigrandmom->vertex().y(), Jpsigrandmom->vertex().z());
0453         } else {
0454           momJpsiID = Jpsimom->pdgId();
0455           trueVtxMom.SetXYZ(Jpsimom->vertex().x(), Jpsimom->vertex().y(), Jpsimom->vertex().z());
0456         }
0457         aBhadron = true;
0458       } else {
0459         if (Jpsigrandmom.isNonnull() && isAbHadron(Jpsigrandmom->pdgId())) {
0460           reco::GenParticleRef JpsiGrandgrandmom = Jpsigrandmom->motherRef();
0461           if (JpsiGrandgrandmom.isNonnull() && isAMixedbHadron(Jpsigrandmom->pdgId(), JpsiGrandgrandmom->pdgId())) {
0462             momJpsiID = JpsiGrandgrandmom->pdgId();
0463             trueVtxMom.SetXYZ(
0464                 JpsiGrandgrandmom->vertex().x(), JpsiGrandgrandmom->vertex().y(), JpsiGrandgrandmom->vertex().z());
0465           } else {
0466             momJpsiID = Jpsigrandmom->pdgId();
0467             trueVtxMom.SetXYZ(Jpsigrandmom->vertex().x(), Jpsigrandmom->vertex().y(), Jpsigrandmom->vertex().z());
0468           }
0469           aBhadron = true;
0470         }
0471       }
0472       if (!aBhadron) {
0473         momJpsiID = Jpsimom->pdgId();
0474         trueVtxMom.SetXYZ(Jpsimom->vertex().x(), Jpsimom->vertex().y(), Jpsimom->vertex().z());
0475       }
0476     }
0477 
0478     TVector3 vdiff = trueVtx - trueVtxMom;
0479     //trueLife = vdiff.Perp()*3.09688/trueP.Perp();
0480     trueLife = vdiff.Perp() * genJpsi->mass() / trueP.Perp();
0481   }
0482   std::pair<int, float> result = std::make_pair(momJpsiID, trueLife);
0483   return result;
0484 }
0485 
0486 void Onia2MuMuPAT::fillDescriptions(edm::ConfigurationDescriptions &iDescriptions) {
0487   edm::ParameterSetDescription desc;
0488   desc.add<edm::InputTag>("muons");
0489   desc.add<edm::InputTag>("beamSpotTag");
0490   desc.add<edm::InputTag>("primaryVertexTag");
0491   desc.add<std::string>("higherPuritySelection");
0492   desc.add<std::string>("lowerPuritySelection");
0493   desc.add<std::string>("dimuonSelection", "");
0494   desc.add<bool>("addCommonVertex");
0495   desc.add<bool>("addMuonlessPrimaryVertex");
0496   desc.add<bool>("resolvePileUpAmbiguity");
0497   desc.add<bool>("addMCTruth");
0498 
0499   iDescriptions.addDefault(desc);
0500 }
0501 
0502 //define this as a plug-in
0503 DEFINE_FWK_MODULE(Onia2MuMuPAT);