L3MuonProducer/src/L3TkMuonProducer.cc

0001 /**  \class L3TkMuonProducer
0002  *
0003  *    This module creates a skimed list of reco::Track (pointing to the original TrackExtra and TrackingRecHitOwnedVector
0004  *    One highest pT track per L1/L2 is selected, requiring some quality.
0005  *
0006  *   \author  J-R Vlimant.
0007  */
0008
0009 // Framework
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "DataFormats/Common/interface/Handle.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include "RecoMuon/L3MuonProducer/src/L3TkMuonProducer.h"
0017 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0018
0019 #include <string>
0020
0021 using namespace edm;
0022 using namespace std;
0023 using namespace reco;
0024
0025 /// constructor with config
0026 L3TkMuonProducer::L3TkMuonProducer(const ParameterSet& parameterSet) {
0027   LogTrace("Muon|RecoMuon|L3TkMuonProducer") << " constructor called";
0028
0029   // StandAlone Collection Label
0030   theL3CollectionLabel = parameterSet.getParameter<InputTag>("InputObjects");
0031   trackToken_ = consumes<reco::TrackCollection>(theL3CollectionLabel);
0032   produces<TrackCollection>();
0033   produces<TrackExtraCollection>();
0034   produces<TrackingRecHitCollection>();
0035
0036   callWhenNewProductsRegistered([this](const edm::BranchDescription& iBD) {
0037     edm::TypeID id(typeid(L3MuonTrajectorySeedCollection));
0038     if (iBD.unwrappedTypeID() == id) {
0039       this->mayConsume<L3MuonTrajectorySeedCollection>(
0040           edm::InputTag{iBD.moduleLabel(), iBD.productInstanceName(), iBD.processName()});
0041     }
0042   });
0043 }
0044
0045 /// destructor
0046 L3TkMuonProducer::~L3TkMuonProducer() {
0047   LogTrace("Muon|RecoMuon|L3TkMuonProducer") << " L3TkMuonProducer destructor called";
0048 }
0049
0050 bool L3TkMuonProducer::sharedSeed(const L3MuonTrajectorySeed& s1, const L3MuonTrajectorySeed& s2) {
0051   //quit right away on nH=0
0052   if (s1.nHits() == 0 || s2.nHits() == 0)
0053     return false;
0054   //quit right away if not the same number of hits
0055   if (s1.nHits() != s2.nHits())
0056     return false;
0057   auto const& r1 = s1.recHits();
0058   auto const& r2 = s2.recHits();
0059   //quit right away if first detId does not match. front exist because of ==0 ->quit test
0060   if (r1.begin()->geographicalId() != r2.begin()->geographicalId())
0061     return false;
0062   //then check hit by hit if they are the same
0063   for (auto i1 = r1.begin(), i2 = r2.begin(); i1 != r1.end() && i2 != r2.end(); ++i1, ++i2) {
0064     if (!i1->sharesInput(&(*i2), TrackingRecHit::all))
0065       return false;
0066   }
0067   return true;
0068 }
0069
0070 string printvector(const vector<TrackRef>& v) {
0071   std::stringstream ss;
0072   for (unsigned int i = 0; i != v.size(); ++i) {
0073     if (i != 0)
0074       ss << "\n";
0075     ss << "track with ref: " << v[i].id().id() << ":" << v[i].key() << " and pT: " << v[i]->pt()
0076        << " with seedRef: " << v[i]->seedRef().id().id() << ":" << v[i]->seedRef().key();
0077   }
0078   return ss.str();
0079 }
0080
0081 string printvector(const vector<L3TkMuonProducer::SeedRef>& v) {
0082   std::stringstream ss;
0083   for (unsigned int i = 0; i != v.size(); ++i) {
0084     if (i != 0)
0085       ss << "\n";
0086     ss << "seed ref: " << v[i].id().id() << ":" << v[i].key();
0087     if (v[i]->l2Track().isNull())
0088       ss << " and pT: " << v[i]->l1Particle()->pt() << " of L1: " << v[i]->l1Particle().id().id() << ":"
0089          << v[i]->l1Particle().key();
0090     else
0091       ss << " and pT: " << v[i]->l2Track()->pt() << " of L2: " << v[i]->l2Track().id().id() << ":"
0092          << v[i]->l2Track().key();
0093   }
0094   return ss.str();
0095 }
0096
0097 string printseed(const L3TkMuonProducer::SeedRef& s) {
0098   std::stringstream ss;
0099   ss << " seed ref: " << s.id().id() << ":" << s.key() << " has " << s->nHits() << "rechits";
0100   for (auto const& hit : s->recHits()) {
0101     ss << "\n detId: " << std::hex << hit.geographicalId().rawId() << std::dec << " position: " << hit.localPosition()
0102        << " and error: " << hit.localPositionError();
0103   }
0104   return ss.str();
0105 }
0106
0107 /// reconstruct muons
0108 void L3TkMuonProducer::produce(Event& event, const EventSetup& eventSetup) {
0109   const string metname = "Muon|RecoMuon|L3TkMuonProducer";
0110
0111   // Take the L3 container
0112   LogDebug(metname) << " Taking the L3/GLB muons: " << theL3CollectionLabel.label();
0113   Handle<TrackCollection> tracks;
0114   event.getByToken(trackToken_, tracks);
0115
0116   //make the LX->L3s pools
0117   LXtoL3sMap LXtoL3s;
0118
0119   unsigned int maxI = tracks->size();
0120   bool gotL3seeds = false;
0121   edm::Handle<L3MuonTrajectorySeedCollection> l3seeds;
0122
0123   //make a list of reference to tracker tracks
0124   vector<TrackRef> orderedTrackTracks(maxI);
0125   for (unsigned int i = 0; i != maxI; i++)
0126     orderedTrackTracks[i] = TrackRef(tracks, i);
0127   LogDebug(metname) << "vector of L3 tracks before ordering:\n" << printvector(orderedTrackTracks);
0128   //order them in pT
0129   sort(orderedTrackTracks.begin(), orderedTrackTracks.end(), trackRefBypT);
0130   LogDebug(metname) << "vector of L3 tracks after ordering:\n" << printvector(orderedTrackTracks);
0131   //loop over then
0132   for (unsigned int i = 0; i != maxI; i++) {
0133     TrackRef& tk = orderedTrackTracks[i];
0134     SeedRef l3seedRef = tk->seedRef().castTo<SeedRef>();
0135
0136     vector<SeedRef> allPossibleOrderedLx;  // with identical hit-set
0137     //add the direct relation
0138     allPossibleOrderedLx.push_back(l3seedRef);
0139     LogDebug(metname) << "adding the seed ref: " << l3seedRef.id().id() << ":" << l3seedRef.key()
0140                       << " for this tracker track: " << tk.id().id() << ":" << tk.key();
0141
0142     //add the relations due to shared seeds
0143     //check whether there is a "shared" seed in addition
0144     if (!gotL3seeds) {
0145       //need to fetch the handle from the ref
0146       const edm::StableProvenance& seedsProv = event.getStableProvenance(l3seedRef.id());
0147       edm::InputTag l3seedsTag(seedsProv.moduleLabel(), seedsProv.productInstanceName(), seedsProv.processName());
0148       event.getByLabel(l3seedsTag, l3seeds);
0149       gotL3seeds = true;
0150       LogDebug(metname) << "got seeds handle from: " << l3seedsTag;
0151     }
0152     //loop the other seeds in the collection
0153     for (unsigned int iS = 0; iS != l3seeds->size(); ++iS) {
0154       const L3MuonTrajectorySeed& seed = (*l3seeds)[iS];
0155       const L3MuonTrajectorySeed& thisSeed = *l3seedRef;
0156       if (l3seedRef.key() == iS)
0157         continue;  //not care about this one
0158       //compare this seed with the seed in the collection
0159       if (sharedSeed(seed, thisSeed)) {
0160         SeedRef thisSharedSeedRef(l3seeds, iS);
0161         LogDebug(metname) << "shared seeds: \n"
0162                           << printseed(l3seedRef) << " and: \n"
0163                           << printseed(thisSharedSeedRef)
0164                           << "\nadding ANOTHER seed ref: " << thisSharedSeedRef.id().id() << ":"
0165                           << thisSharedSeedRef.key() << " for this tracker track: " << tk.id().id() << ":" << tk.key();
0166         //  edm::LogError(metname)<<" we have a shared seed right there.";
0167         allPossibleOrderedLx.push_back(thisSharedSeedRef);
0168       }  //seed is shared
0169     }    //loop all other existing seed for overlaps
0170
0171     //now you have the full list of Lx objects that have seed this tracker track.
0172     // order the list in pT of Lx objects
0173     LogDebug(metname) << "list of possible Lx objects for tracker track: " << tk.id().id() << ":" << tk.key()
0174                       << " before ordering\n"
0175                       << printvector(allPossibleOrderedLx);
0176     sort(allPossibleOrderedLx.begin(), allPossibleOrderedLx.end(), seedRefBypT);
0177     LogDebug(metname) << "list of possible Lx objects for tracker track: " << tk.id().id() << ":" << tk.key()
0178                       << " after ordering\n"
0179                       << printvector(allPossibleOrderedLx);
0180     // assign this tracker track to the highest pT Lx.
0181     for (unsigned int iL = 0; iL != allPossibleOrderedLx.size(); ++iL) {
0182       SeedRef thisRef = allPossibleOrderedLx[iL];
0183       pseudoRef ref = makePseudoRef(*thisRef);
0184       LogDebug(metname) << "seed ref: " << thisRef.id().id() << ":" << thisRef.key()
0185                         << " transcribe to pseudoref: " << ref.first << ":" << ref.second;
0186       LXtoL3sMap::iterator f = LXtoL3s.find(ref);
0187       if (f != LXtoL3s.end()) {
0188         //there's already an entry. because of the prior ordering in pT of the tracker track refs
0189         // the track ref already there *has* a higher pT: this one cannot compete and should be assigned to the next Lx;
0190         LogDebug(metname) << "this tracker track: " << tk.id().id() << ":" << tk.key() << " (" << tk->pt() << ")"
0191                           << "\n cannot compete in pT with track: " << f->second.first.id().id() << ":"
0192                           << f->second.first.key() << " (" << f->second.first->pt() << ")"
0193                           << "\n already assigned to pseudo ref: " << ref.first << ":" << ref.second
0194                           << " which corresponds to seedRef: " << f->second.second.id().id() << ":"
0195                           << f->second.second.key();
0196         continue;
0197       } else {
0198         //there was no entry yet. make the assignement
0199         LogDebug(metname) << "this tracker track: " << tk.id().id() << ":" << tk.key()
0200                           << " is assigned to pseudo ref: " << ref.first << ":" << ref.second
0201                           << " which corresponds to seedRef: " << thisRef.id().id() << ":" << thisRef.key();
0202         LXtoL3s[ref] = std::make_pair(tk, thisRef);
0203         //once assigned. break
0204         break;
0205       }
0206     }  //loop possible Lx for possible assignement
0207   }    //loop over ordered list of tracker track refs
0208
0209   //prepare the output
0210   auto outTracks = std::make_unique<TrackCollection>(LXtoL3s.size());
0211   auto outTrackExtras = std::make_unique<TrackExtraCollection>(LXtoL3s.size());
0212   reco::TrackExtraRefProd rTrackExtras = event.getRefBeforePut<TrackExtraCollection>();
0213   auto outRecHits = std::make_unique<TrackingRecHitCollection>();
0214   TrackingRecHitRefProd rHits = event.getRefBeforePut<TrackingRecHitCollection>();
0215
0216   LogDebug(metname) << "reading the map to make " << LXtoL3s.size() << "products.";
0217   //fill the collection from the map
0218   LXtoL3sMap::iterator f = LXtoL3s.begin();
0219   unsigned int i = 0;
0220   for (; f != LXtoL3s.end(); ++f, ++i) {
0221     LogDebug(metname) << "copy the track over, and make ref to extra";
0222     const Track& trk = *(f->second.first);
0223     (*outTracks)[i] = Track(trk);
0224     (*outTracks)[i].setExtra(TrackExtraRef(rTrackExtras, i));
0225
0226     LogDebug(metname) << "copy the trackExtra too, and change the seedref";
0227     edm::RefToBase<TrajectorySeed> seedRef(f->second.second);
0228     //do not use the copy constructor, otherwise the hit Ref are still the same
0229     (*outTrackExtras)[i] = TrackExtra(trk.outerPosition(),
0230                                       trk.outerMomentum(),
0231                                       trk.outerOk(),
0232                                       trk.innerPosition(),
0233                                       trk.innerMomentum(),
0234                                       trk.innerOk(),
0235                                       trk.outerStateCovariance(),
0236                                       trk.outerDetId(),
0237                                       trk.innerStateCovariance(),
0238                                       trk.innerDetId(),
0239                                       seedRef->direction(),
0240                                       seedRef);
0241
0242     LogDebug(metname) << "copy the hits too";
0243     unsigned int iRH = 0;
0244     for (trackingRecHit_iterator hit = trk.recHitsBegin(); hit != trk.recHitsEnd(); ++hit, ++iRH) {
0245       outRecHits->push_back((*hit)->clone());
0246     }
0247     (*outTrackExtras)[i].setHits(rHits, 0, iRH);
0248   }
0249
0250   LogDebug(metname) << "made: " << outTracks->size() << " tracks, " << outTrackExtras->size() << " extras and "
0251                     << outRecHits->size() << " rechits.";
0252
0253   //put the collection in the event
0254   LogDebug(metname) << "loading...";
0255   event.put(std::move(outTracks));
0256   event.put(std::move(outTrackExtras));
0257   event.put(std::move(outRecHits));
0258   LogDebug(metname) << " Event loaded"
0259                     << "================================";
0260 }