Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:13 

0001 #include "MuScleFitMuonSelector.h"
0002 #include "DataFormats/MuonReco/interface/CaloMuon.h"
0003 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0004 #include "MuonAnalysis/MomentumScaleCalibration/interface/Muon.h"
0005 
0006 const double MuScleFitMuonSelector::mMu2 = 0.011163612;
0007 const unsigned int MuScleFitMuonSelector::motherPdgIdArray[] = {23, 100553, 100553, 553, 100443, 443};
0008 
0009 const reco::Candidate* MuScleFitMuonSelector::getStatus1Muon(const reco::Candidate* status3Muon) {
0010   const reco::Candidate* tempMuon = status3Muon;
0011   //  bool lastCopy = (static_cast<const reco::GenParticle*>(tempMuon))->isLastCopy();                      //  isLastCopy() likely not enough robust
0012   bool isPromptFinalState = static_cast<const reco::GenParticle*>(tempMuon)
0013                                 ->isPromptFinalState();  //  pre-CMSSW_74X code: int status = tempStatus1Muon->status();
0014   while (tempMuon == nullptr || tempMuon->numberOfDaughters() != 0) {
0015     if (isPromptFinalState)
0016       break;  //  pre-CMSSW_74X code: if (status == 1) break;
0017     //std::vector<const reco::Candidate*> daughters;
0018     for (unsigned int i = 0; i < tempMuon->numberOfDaughters(); ++i) {
0019       if (tempMuon->daughter(i)->pdgId() == tempMuon->pdgId()) {
0020         tempMuon = tempMuon->daughter(i);
0021         isPromptFinalState = static_cast<const reco::GenParticle*>(tempMuon)
0022                                  ->isPromptFinalState();  //  pre-CMSSW_74X code: status = tempStatus1Muon->status();
0023         break;
0024       } else
0025         continue;
0026     }  //for loop
0027   }  //while loop
0028 
0029   return tempMuon;
0030 }
0031 
0032 const reco::Candidate* MuScleFitMuonSelector::getStatus3Muon(const reco::Candidate* status3Muon) {
0033   const reco::Candidate* tempMuon = status3Muon;
0034   bool lastCopy = static_cast<const reco::GenParticle*>(tempMuon)
0035                       ->isLastCopyBeforeFSR();  //  pre-CMSSW_74X code: int status = tempStatus1Muon->status();
0036   while (tempMuon == nullptr || tempMuon->numberOfDaughters() != 0) {
0037     if (lastCopy)
0038       break;  //  pre-CMSSW_74X code: if (status == 3) break;
0039     //std::vector<const reco::Candidate*> daughters;
0040     for (unsigned int i = 0; i < tempMuon->numberOfDaughters(); ++i) {
0041       if (tempMuon->daughter(i)->pdgId() == tempMuon->pdgId()) {
0042         tempMuon = tempMuon->daughter(i);
0043         lastCopy = static_cast<const reco::GenParticle*>(tempMuon)
0044                        ->isLastCopyBeforeFSR();  //  pre-CMSSW_74X code: status = tempStatus1Muon->status();
0045         break;
0046       } else
0047         continue;
0048     }  //for loop
0049   }  //while loop
0050 
0051   return tempMuon;
0052 }
0053 
0054 bool MuScleFitMuonSelector::selGlobalMuon(const pat::Muon* aMuon) {
0055   reco::TrackRef iTrack = aMuon->innerTrack();
0056   const reco::HitPattern& p = iTrack->hitPattern();
0057 
0058   reco::TrackRef gTrack = aMuon->globalTrack();
0059   const reco::HitPattern& q = gTrack->hitPattern();
0060 
0061   return (  //isMuonInAccept(aMuon) &&// no acceptance cuts!
0062       iTrack->found() > 11 && gTrack->chi2() / gTrack->ndof() < 20.0 && q.numberOfValidMuonHits() > 0 &&
0063       iTrack->chi2() / iTrack->ndof() < 4.0 && aMuon->muonID("TrackerMuonArbitrated") &&
0064       aMuon->muonID("TMLastStationAngTight") && p.pixelLayersWithMeasurement() > 1 &&
0065       std::abs(iTrack->dxy()) < 3.0 &&  //should be done w.r.t. PV!
0066       std::abs(iTrack->dz()) < 15.0     //should be done w.r.t. PV!
0067   );
0068 }
0069 
0070 bool MuScleFitMuonSelector::selTrackerMuon(const pat::Muon* aMuon) {
0071   reco::TrackRef iTrack = aMuon->innerTrack();
0072   const reco::HitPattern& p = iTrack->hitPattern();
0073 
0074   return (  //isMuonInAccept(aMuon) // no acceptance cuts!
0075       iTrack->found() > 11 && iTrack->chi2() / iTrack->ndof() < 4.0 && aMuon->muonID("TrackerMuonArbitrated") &&
0076       aMuon->muonID("TMLastStationAngTight") && p.pixelLayersWithMeasurement() > 1 &&
0077       std::abs(iTrack->dxy()) < 3.0 &&  //should be done w.r.t. PV!
0078       std::abs(iTrack->dz()) < 15.0     //should be done w.r.t. PV!
0079   );
0080 }
0081 
0082 // Note that at this level we save all the information. Events for which no suitable muon pair is found
0083 // are removed from the tree (together with the gen and sim information) at a later stage.
0084 // It would be better to remove them directly at this point.
0085 void MuScleFitMuonSelector::selectMuons(const edm::Event& event,
0086                                         std::vector<MuScleFitMuon>& muons,
0087                                         std::vector<GenMuonPair>& genPair,
0088                                         std::vector<std::pair<lorentzVector, lorentzVector> >& simPair,
0089                                         MuScleFitPlotter* plotter) {
0090   edm::Handle<pat::CompositeCandidateCollection> collAll = event.getHandle(onia2MuMuToken_);
0091   if (!collAll.isValid()) {
0092     edm::LogWarning("MuScleFitUtils") << "J/psi not present in event!";
0093   }
0094   std::vector<const pat::Muon*> collMuSel;
0095 
0096   //================onia cuts===========================/
0097 
0098   if (muonType_ <= -1 && PATmuons_) {
0099     std::vector<const pat::CompositeCandidate*> collSelGG;
0100     std::vector<const pat::CompositeCandidate*> collSelGT;
0101     std::vector<const pat::CompositeCandidate*> collSelTT;
0102     if (collAll.isValid()) {
0103       for (std::vector<pat::CompositeCandidate>::const_iterator it = collAll->begin(); it != collAll->end(); ++it) {
0104         const pat::CompositeCandidate* cand = &(*it);
0105         // cout << "Now checking candidate of type " << theJpsiCat << " with pt = " << cand->pt() << endl;
0106         const pat::Muon* muon1 = dynamic_cast<const pat::Muon*>(cand->daughter("muon1"));
0107         const pat::Muon* muon2 = dynamic_cast<const pat::Muon*>(cand->daughter("muon2"));
0108 
0109         if ((muon1->charge() * muon2->charge()) > 0)
0110           continue;
0111         // global + global?
0112         if (muon1->isGlobalMuon() && muon2->isGlobalMuon()) {
0113           if (selGlobalMuon(muon1) && selGlobalMuon(muon2) && cand->userFloat("vProb") > 0.001) {
0114             collSelGG.push_back(cand);
0115             continue;
0116           }
0117         }
0118         // global + tracker? (x2)
0119         if (muon1->isGlobalMuon() && muon2->isTrackerMuon()) {
0120           if (selGlobalMuon(muon1) && selTrackerMuon(muon2) && cand->userFloat("vProb") > 0.001) {
0121             collSelGT.push_back(cand);
0122             continue;
0123           }
0124         }
0125         if (muon2->isGlobalMuon() && muon1->isTrackerMuon()) {
0126           if (selGlobalMuon(muon2) && selTrackerMuon(muon1) && cand->userFloat("vProb") > 0.001) {
0127             collSelGT.push_back(cand);
0128             continue;
0129           }
0130         }
0131         // tracker + tracker?
0132         if (muon1->isTrackerMuon() && muon2->isTrackerMuon()) {
0133           if (selTrackerMuon(muon1) && selTrackerMuon(muon2) && cand->userFloat("vProb") > 0.001) {
0134             collSelTT.push_back(cand);
0135             continue;
0136           }
0137         }
0138       }
0139     }
0140     // Split them in independent collections if using muonType_ == -2, -3 or -4. Take them all if muonType_ == -1.
0141     std::vector<reco::Track> tracks;
0142     if (!collSelGG.empty()) {
0143       //CHECK THAT THEY ARE ORDERED BY PT !!!!!!!!!!!!!!!!!!!!!!!
0144       const pat::Muon* muon1 = dynamic_cast<const pat::Muon*>(collSelGG[0]->daughter("muon1"));
0145       const pat::Muon* muon2 = dynamic_cast<const pat::Muon*>(collSelGG[0]->daughter("muon2"));
0146       if (muonType_ == -1 || muonType_ == -2) {
0147         tracks.push_back(*(muon1->innerTrack()));
0148         tracks.push_back(*(muon2->innerTrack()));
0149         collMuSel.push_back(muon1);
0150         collMuSel.push_back(muon2);
0151       }
0152     } else if (collSelGG.empty() && !collSelGT.empty()) {
0153       //CHECK THAT THEY ARE ORDERED BY PT !!!!!!!!!!!!!!!!!!!!!!!
0154       const pat::Muon* muon1 = dynamic_cast<const pat::Muon*>(collSelGT[0]->daughter("muon1"));
0155       const pat::Muon* muon2 = dynamic_cast<const pat::Muon*>(collSelGT[0]->daughter("muon2"));
0156       if (muonType_ == -1 || muonType_ == -3) {
0157         tracks.push_back(*(muon1->innerTrack()));
0158         tracks.push_back(*(muon2->innerTrack()));
0159         collMuSel.push_back(muon1);
0160         collMuSel.push_back(muon2);
0161       }
0162     } else if (collSelGG.empty() && collSelGT.empty() && !collSelTT.empty()) {
0163       //CHECK THAT THEY ARE ORDERED BY PT !!!!!!!!!!!!!!!!!!!!!!!
0164       const pat::Muon* muon1 = dynamic_cast<const pat::Muon*>(collSelTT[0]->daughter("muon1"));
0165       const pat::Muon* muon2 = dynamic_cast<const pat::Muon*>(collSelTT[0]->daughter("muon2"));
0166       if (muonType_ == -1 || muonType_ == -4) {
0167         tracks.push_back(*(muon1->innerTrack()));
0168         tracks.push_back(*(muon2->innerTrack()));
0169         collMuSel.push_back(muon1);
0170         collMuSel.push_back(muon2);
0171       }
0172     }
0173     if (tracks.size() != 2 && !tracks.empty()) {
0174       std::cout << "ERROR strange number of muons selected by onia cuts!" << std::endl;
0175       abort();
0176     }
0177     muons = fillMuonCollection(tracks);
0178   } else if ((muonType_ < 4 && muonType_ >= 0) || muonType_ >= 10) {  // Muons (glb,sta,trk)
0179     std::vector<reco::Track> tracks;
0180     if (PATmuons_ == true) {
0181       edm::Handle<pat::MuonCollection> allMuons = event.getHandle(patMuonToken_);
0182       if (muonType_ == 0) {
0183         // Take directly the muon
0184         muons = fillMuonCollection(*allMuons);
0185       } else {
0186         for (std::vector<pat::Muon>::const_iterator muon = allMuons->begin(); muon != allMuons->end(); ++muon) {
0187           //std::cout<<"pat muon is global "<<muon->isGlobalMuon()<<std::endl;
0188           takeSelectedMuonType(muon, tracks);
0189         }
0190         muons = fillMuonCollection(tracks);
0191       }
0192     } else {
0193       edm::Handle<reco::MuonCollection> allMuons = event.getHandle(recoMuonToken_);
0194       if (muonType_ == 0) {
0195         // Take directly the muon
0196         muons = fillMuonCollection(*allMuons);
0197       } else {
0198         for (std::vector<reco::Muon>::const_iterator muon = allMuons->begin(); muon != allMuons->end(); ++muon) {
0199           takeSelectedMuonType(muon, tracks);
0200         }
0201         muons = fillMuonCollection(tracks);
0202       }
0203     }
0204   } else if (muonType_ == 4) {  //CaloMuons
0205     edm::Handle<reco::CaloMuonCollection> caloMuons = event.getHandle(caloMuonToken_);
0206     std::vector<reco::Track> tracks;
0207     for (std::vector<reco::CaloMuon>::const_iterator muon = caloMuons->begin(); muon != caloMuons->end(); ++muon) {
0208       tracks.push_back(*(muon->track()));
0209     }
0210     muons = fillMuonCollection(tracks);
0211   }
0212 
0213   else if (muonType_ == 5) {  // Inner tracker tracks
0214     edm::Handle<reco::TrackCollection> tracks = event.getHandle(trackCollectionToken_);
0215     muons = fillMuonCollection(*tracks);
0216   }
0217   plotter->fillRec(muons);
0218 
0219   // Generation and simulation part
0220   if (speedup_ == false) {
0221     if (PATmuons_) {
0222       // EM 2015.04.02 temporary fix to run on MINIAODSIM (which contains PAT muons) but not the "onia2MuMuPatTrkTrk" collection
0223       // selectGeneratedMuons(collAll, collMuSel, genPair, plotter);
0224       selectGenSimMuons(event, genPair, simPair, plotter);
0225     } else {
0226       selectGenSimMuons(event, genPair, simPair, plotter);
0227     }
0228   }
0229 }
0230 
0231 /// For PATmuons the generator information is read directly from the PAT object
0232 void MuScleFitMuonSelector::selectGeneratedMuons(const edm::Handle<pat::CompositeCandidateCollection>& collAll,
0233                                                  const std::vector<const pat::Muon*>& collMuSel,
0234                                                  std::vector<GenMuonPair>& genPair,
0235                                                  MuScleFitPlotter* plotter) {
0236   reco::GenParticleCollection genPatParticles{};
0237 
0238   //explicitly for JPsi but can be adapted!!!!!
0239   for (std::vector<pat::CompositeCandidate>::const_iterator it = collAll->begin(); it != collAll->end(); ++it) {
0240     reco::GenParticleRef genJpsi = it->genParticleRef();
0241     bool isMatched = (genJpsi.isAvailable() && genJpsi->pdgId() == 443);
0242     if (isMatched) {
0243       genPatParticles.push_back(*genJpsi.get());
0244     }
0245   }
0246 
0247   if (collMuSel.size() == 2) {
0248     reco::GenParticleRef genMu1 = collMuSel[0]->genParticleRef();
0249     reco::GenParticleRef genMu2 = collMuSel[1]->genParticleRef();
0250     bool isMuMatched = (genMu1.isAvailable() && genMu2.isAvailable() && genMu1->pdgId() * genMu2->pdgId() == -169);
0251     if (isMuMatched) {
0252       genPatParticles.push_back(*genMu1.get());
0253       genPatParticles.push_back(*genMu2.get());
0254 
0255       unsigned int motherId = 0;
0256       if (genMu1->mother() != nullptr) {
0257         motherId = genMu1->mother()->pdgId();
0258       }
0259       if (genMu1->pdgId() == 13)
0260         genPair.push_back(GenMuonPair(genMu1.get()->p4(), genMu2.get()->p4(), motherId));
0261       else
0262         // genPair.push_back(std::make_pair(genMu2.get()->p4(),genMu1.get()->p4()) );
0263         genPair.push_back(GenMuonPair(genMu2.get()->p4(), genMu1.get()->p4(), motherId));
0264 
0265       plotter->fillGen(genPatParticles, true);
0266 
0267       if (debug_ > 0)
0268         std::cout << "Found genParticles in PAT" << std::endl;
0269     } else {
0270       std::cout << "No recomuon selected so no access to generated info" << std::endl;
0271       // Fill it in any case, otherwise it will not be in sync with the event number
0272       // genPair.push_back( std::make_pair( lorentzVector(0.,0.,0.,0.), lorentzVector(0.,0.,0.,0.) ) );
0273       genPair.push_back(GenMuonPair(lorentzVector(0., 0., 0., 0.), lorentzVector(0., 0., 0., 0.), 0));
0274     }
0275   } else {
0276     std::cout << "No recomuon selected so no access to generated info" << std::endl;
0277     // Fill it in any case, otherwise it will not be in sync with the event number
0278     // genPair.push_back( std::make_pair( lorentzVector(0.,0.,0.,0.), lorentzVector(0.,0.,0.,0.) ) );
0279     genPair.push_back(GenMuonPair(lorentzVector(0., 0., 0., 0.), lorentzVector(0., 0., 0., 0.), 0));
0280   }
0281   if (debug_ > 0) {
0282     std::cout << "genParticles:" << std::endl;
0283     // std::cout << genPair.back().first << " , " << genPair.back().second << std::endl;
0284     std::cout << genPair.back().mu1 << " , " << genPair.back().mu2 << std::endl;
0285   }
0286 }
0287 
0288 void MuScleFitMuonSelector::selectGenSimMuons(const edm::Event& event,
0289                                               std::vector<GenMuonPair>& genPair,
0290                                               std::vector<std::pair<lorentzVector, lorentzVector> >& simPair,
0291                                               MuScleFitPlotter* plotter) {
0292   // Find and store in histograms the generated and simulated resonance and muons
0293   // ----------------------------------------------------------------------------
0294   edm::Handle<edm::HepMCProduct> evtMC = event.getHandle(evtMCToken_);
0295   edm::Handle<reco::GenParticleCollection> genParticles = event.getHandle(genParticleToken_);
0296 
0297   // Fill gen information only in the first loop
0298   bool ifHepMC = false;
0299   if (evtMC.isValid()) {
0300     genPair.push_back(findGenMuFromRes(evtMC.product()));
0301     plotter->fillGen(*evtMC, sherpa_);
0302     ifHepMC = true;
0303     if (debug_ > 0)
0304       std::cout << "Found hepMC" << std::endl;
0305   } else if (genParticles.isValid()) {
0306     genPair.push_back(findGenMuFromRes(genParticles.product()));
0307     plotter->fillGen(*genParticles);
0308     if (debug_ > 0)
0309       std::cout << "Found genParticles" << std::endl;
0310   } else {
0311     std::cout << "ERROR "
0312               << "non generation info and speedup true!!!!!!!!!!!!" << std::endl;
0313     // Fill it in any case, otherwise it will not be in sync with the event number
0314     // genPair.push_back( std::make_pair( lorentzVector(0.,0.,0.,0.), lorentzVector(0.,0.,0.,0.) ) );
0315     genPair.push_back(GenMuonPair(lorentzVector(0., 0., 0., 0.), lorentzVector(0., 0., 0., 0.), 0));
0316   }
0317   if (debug_ > 0) {
0318     std::cout << "genParticles:" << std::endl;
0319     std::cout << genPair.back().mu1 << " , " << genPair.back().mu2 << std::endl;
0320   }
0321   selectSimulatedMuons(event, ifHepMC, evtMC, simPair, plotter);
0322 }
0323 
0324 void MuScleFitMuonSelector::selectSimulatedMuons(const edm::Event& event,
0325                                                  const bool ifHepMC,
0326                                                  edm::Handle<edm::HepMCProduct> evtMC,
0327                                                  std::vector<std::pair<lorentzVector, lorentzVector> >& simPair,
0328                                                  MuScleFitPlotter* plotter) {
0329   edm::Handle<edm::SimTrackContainer> simTracks = event.getHandle(simTrackToken_);
0330   bool simTracksFound = false;
0331   if (simTracks.isValid()) {
0332     plotter->fillSim(simTracks);
0333     if (ifHepMC) {
0334       simPair.push_back(findSimMuFromRes(evtMC, simTracks));
0335       simTracksFound = true;
0336       plotter->fillGenSim(evtMC, simTracks);
0337     }
0338   } else {
0339     std::cout << "SimTracks not existent" << std::endl;
0340   }
0341   if (!simTracksFound) {
0342     simPair.push_back(std::make_pair(lorentzVector(0., 0., 0., 0.), lorentzVector(0., 0., 0., 0.)));
0343   }
0344 }
0345 
0346 GenMuonPair MuScleFitMuonSelector::findGenMuFromRes(const edm::HepMCProduct* evtMC) {
0347   const HepMC::GenEvent* Evt = evtMC->GetEvent();
0348   GenMuonPair muFromRes;
0349   //Loop on generated particles
0350   for (HepMC::GenEvent::particle_const_iterator part = Evt->particles_begin(); part != Evt->particles_end(); part++) {
0351     if (std::abs((*part)->pdg_id()) == 13 && (*part)->status() == 1) {
0352       bool fromRes = false;
0353       unsigned int motherPdgId = 0;
0354       for (HepMC::GenVertex::particle_iterator mother = (*part)->production_vertex()->particles_begin(HepMC::ancestors);
0355            mother != (*part)->production_vertex()->particles_end(HepMC::ancestors);
0356            ++mother) {
0357         motherPdgId = (*mother)->pdg_id();
0358 
0359         // For sherpa the resonance is not saved. The muons from the resonance can be identified
0360         // by having as mother a muon of status 3.
0361         if (sherpa_) {
0362           if (motherPdgId == 13 && (*mother)->status() == 3)
0363             fromRes = true;
0364         } else {
0365           for (int ires = 0; ires < 6; ++ires) {
0366             if (motherPdgId == motherPdgIdArray[ires] && resfind_[ires])
0367               fromRes = true;
0368           }
0369         }
0370       }
0371       if (fromRes) {
0372         if ((*part)->pdg_id() == 13) {
0373           //   muFromRes.first = (*part)->momentum();
0374           muFromRes.mu1 = (lorentzVector(
0375               (*part)->momentum().px(), (*part)->momentum().py(), (*part)->momentum().pz(), (*part)->momentum().e()));
0376         } else {
0377           muFromRes.mu2 = (lorentzVector(
0378               (*part)->momentum().px(), (*part)->momentum().py(), (*part)->momentum().pz(), (*part)->momentum().e()));
0379         }
0380         muFromRes.motherId = motherPdgId;
0381       }
0382     }
0383   }
0384   return muFromRes;
0385 }
0386 
0387 // std::pair<lorentzVector, lorentzVector> MuScleFitMuonSelector::findGenMuFromRes( const reco::GenParticleCollection* genParticles)
0388 GenMuonPair MuScleFitMuonSelector::findGenMuFromRes(const reco::GenParticleCollection* genParticles) {
0389   // std::pair<lorentzVector,lorentzVector> muFromRes;
0390   GenMuonPair muFromRes;
0391 
0392   //Loop on generated particles
0393   if (debug_ > 0)
0394     std::cout << "Starting loop on " << genParticles->size() << " genParticles" << std::endl;
0395   for (reco::GenParticleCollection::const_iterator part = genParticles->begin(); part != genParticles->end(); ++part) {
0396     if (debug_ > 0)
0397       std::cout << "genParticle has pdgId = " << std::abs(part->pdgId()) << " and status = " << part->status()
0398                 << std::endl;
0399     if (std::abs(part->pdgId()) == 13) {  // && part->status()==3) {
0400       bool fromRes = false;
0401       unsigned int motherPdgId = part->mother()->pdgId();
0402       if (debug_ > 0) {
0403         std::cout << "Found a muon with mother: " << motherPdgId << std::endl;
0404       }
0405       for (int ires = 0; ires < 6; ++ires) {
0406         // if( motherPdgId == motherPdgIdArray[ires] && resfind_[ires] ) fromRes = true; // changed by EM 2015.07.30
0407         // begin of comment
0408         // the list of resonances motherPdgIdArray does not contain the photon (PdgId = 21) while ~1% of the
0409         // mu+mu- events in the range [50,120] GeV has a photon as the mother.
0410         // It needs to be fixed without spoiling the logic of the selection of different resonances
0411         // e.g. mixing up onia etc.
0412         // end of comment
0413         if ((motherPdgId == motherPdgIdArray[ires] && resfind_[ires]) || motherPdgId == 21)
0414           fromRes = true;
0415       }
0416       if (fromRes) {
0417         if (debug_ > 0)
0418           std::cout << "fromRes = true, motherPdgId = " << motherPdgId << std::endl;
0419         const reco::Candidate* status3Muon = &(*part);
0420         const reco::Candidate* status1Muon = getStatus1Muon(status3Muon);
0421         if (part->pdgId() == 13) {
0422           if (status1Muon->p4().pt() != 0)
0423             muFromRes.mu1 = MuScleFitMuon(status1Muon->p4(), -1);
0424           else
0425             muFromRes.mu1 = MuScleFitMuon(status3Muon->p4(), -1);
0426           if (debug_ > 0)
0427             std::cout << "Found a genMuon - : " << muFromRes.mu1 << std::endl;
0428           //      muFromRes.first = (lorentzVector(status1Muon->p4().px(),status1Muon->p4().py(),
0429           //                       status1Muon->p4().pz(),status1Muon->p4().e()));
0430         } else {
0431           if (status1Muon->p4().pt() != 0)
0432             muFromRes.mu2 = MuScleFitMuon(status1Muon->p4(), +1);
0433           else
0434             muFromRes.mu2 = MuScleFitMuon(status3Muon->p4(), +1);
0435           if (debug_ > 0)
0436             std::cout << "Found a genMuon + : " << muFromRes.mu2 << std::endl;
0437           //      muFromRes.second = (lorentzVector(status1Muon->p4().px(),status1Muon->p4().py(),
0438           //                        status1Muon->p4().pz(),status1Muon->p4().e()));
0439         }
0440         muFromRes.motherId = motherPdgId;
0441       }
0442     }
0443   }
0444   return muFromRes;
0445 }
0446 
0447 std::pair<lorentzVector, lorentzVector> MuScleFitMuonSelector::findSimMuFromRes(
0448     const edm::Handle<edm::HepMCProduct>& evtMC, const edm::Handle<edm::SimTrackContainer>& simTracks) {
0449   //Loop on simulated tracks
0450   std::pair<lorentzVector, lorentzVector> simMuFromRes;
0451   for (edm::SimTrackContainer::const_iterator simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0452     //Chose muons
0453     if (std::abs((*simTrack).type()) == 13) {
0454       //If tracks from IP than find mother
0455       if ((*simTrack).genpartIndex() > 0) {
0456         HepMC::GenParticle* gp = evtMC->GetEvent()->barcode_to_particle((*simTrack).genpartIndex());
0457         if (gp != nullptr) {
0458           for (HepMC::GenVertex::particle_iterator mother = gp->production_vertex()->particles_begin(HepMC::ancestors);
0459                mother != gp->production_vertex()->particles_end(HepMC::ancestors);
0460                ++mother) {
0461             bool fromRes = false;
0462             unsigned int motherPdgId = (*mother)->pdg_id();
0463             for (int ires = 0; ires < 6; ++ires) {
0464               if ((motherPdgId == motherPdgIdArray[ires] && resfind_[ires]) || motherPdgId == 21)
0465                 fromRes = true;
0466             }
0467             if (fromRes) {
0468               if (gp->pdg_id() == 13)
0469                 simMuFromRes.first = lorentzVector(simTrack->momentum().px(),
0470                                                    simTrack->momentum().py(),
0471                                                    simTrack->momentum().pz(),
0472                                                    simTrack->momentum().e());
0473               else
0474                 simMuFromRes.second = lorentzVector(simTrack->momentum().px(),
0475                                                     simTrack->momentum().py(),
0476                                                     simTrack->momentum().pz(),
0477                                                     simTrack->momentum().e());
0478             }
0479           }
0480         }
0481         // else LogDebug("MuScleFitUtils") << "WARNING: no matching genParticle found for simTrack" << std::endl;
0482       }
0483     }
0484   }
0485   return simMuFromRes;
0486 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!