L2MuonSeedGenerator/src/L2MuonSeedGenerator.cc

0001 //-------------------------------------------------
0002 //
0003 /**  \class L2MuonSeedGenerator
0004  *
0005  *   L2 muon seed generator:
0006  *   Transform the L1 informations in seeds for the
0007  *   L2 muon reconstruction
0008  *
0009  *
0010  *
0011  *   \author  A.Everett, R.Bellan, J. Alcaraz
0012  *
0013  *    ORCA's author: N. Neumeister
0014  */
0015 //
0016 //--------------------------------------------------
0017
0018 // Class Header
0019 #include "RecoMuon/L2MuonSeedGenerator/src/L2MuonSeedGenerator.h"
0020
0021 // Framework
0022 #include "FWCore/Framework/interface/ConsumesCollector.h"
0023 #include "FWCore/Framework/interface/EventSetup.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0027
0028 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0029 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
0030 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0031 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0032 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0033 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0034 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0035
0036 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0037 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0038
0039 using namespace std;
0040 using namespace edm;
0041 using namespace l1extra;
0042
0043 // constructors
0044 L2MuonSeedGenerator::L2MuonSeedGenerator(const edm::ParameterSet& iConfig)
0045     : theSource(iConfig.getParameter<InputTag>("InputObjects")),
0046       theL1GMTReadoutCollection(iConfig.getParameter<InputTag>("GMTReadoutCollection")),
0047       thePropagatorName(iConfig.getParameter<string>("Propagator")),
0048       theL1MinPt(iConfig.getParameter<double>("L1MinPt")),
0049       theL1MaxEta(iConfig.getParameter<double>("L1MaxEta")),
0050       theL1MinQuality(iConfig.getParameter<unsigned int>("L1MinQuality")),
0051       useOfflineSeed(iConfig.getUntrackedParameter<bool>("UseOfflineSeed", false)),
0052       useUnassociatedL1(iConfig.existsAs<bool>("UseUnassociatedL1") ? iConfig.getParameter<bool>("UseUnassociatedL1")
0053                                                                     : true) {
0054   gmtToken_ = consumes<L1MuGMTReadoutCollection>(theL1GMTReadoutCollection);
0055   muCollToken_ = consumes<L1MuonParticleCollection>(theSource);
0056
0057   if (useOfflineSeed) {
0058     theOfflineSeedLabel = iConfig.getUntrackedParameter<InputTag>("OfflineSeedLabel");
0059     offlineSeedToken_ = consumes<edm::View<TrajectorySeed> >(theOfflineSeedLabel);
0060   }
0061
0062   // service parameters
0063   ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
0064
0065   // the services
0066   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0067
0068   // the estimator
0069   theEstimator = new Chi2MeasurementEstimator(10000.);
0070
0071   produces<L2MuonTrajectorySeedCollection>();
0072 }
0073
0074 // destructor
0075 L2MuonSeedGenerator::~L2MuonSeedGenerator() {
0076   if (theService)
0077     delete theService;
0078   if (theEstimator)
0079     delete theEstimator;
0080 }
0081
0082 void L2MuonSeedGenerator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0083   const std::string metname = "Muon|RecoMuon|L2MuonSeedGenerator";
0084   MuonPatternRecoDumper debug;
0085
0086   auto output = std::make_unique<L2MuonTrajectorySeedCollection>();
0087
0088   // Muon particles and GMT readout collection
0089   edm::Handle<L1MuGMTReadoutCollection> gmtrc_handle;
0090   iEvent.getByToken(gmtToken_, gmtrc_handle);
0091   L1MuGMTReadoutRecord const& gmtrr = gmtrc_handle.product()->getRecord(0);
0092
0093   edm::Handle<L1MuonParticleCollection> muColl;
0094   iEvent.getByToken(muCollToken_, muColl);
0095   LogTrace(metname) << "Number of muons " << muColl->size() << endl;
0096
0097   edm::Handle<edm::View<TrajectorySeed> > offlineSeedHandle;
0098   vector<int> offlineSeedMap;
0099   if (useOfflineSeed) {
0100     iEvent.getByToken(offlineSeedToken_, offlineSeedHandle);
0101     LogTrace(metname) << "Number of offline seeds " << offlineSeedHandle->size() << endl;
0102     offlineSeedMap = vector<int>(offlineSeedHandle->size(), 0);
0103   }
0104
0105   L1MuonParticleCollection::const_iterator it;
0106   L1MuonParticleRef::key_type l1ParticleIndex = 0;
0107
0108   for (it = muColl->begin(); it != muColl->end(); ++it, ++l1ParticleIndex) {
0109     const L1MuGMTExtendedCand muonCand = (*it).gmtMuonCand();
0110     unsigned int quality = 0;
0111     bool valid_charge = false;
0112     ;
0113
0114     if (muonCand.empty()) {
0115       LogWarning(metname) << "L2MuonSeedGenerator: WARNING, no L1MuGMTCand! " << endl;
0116       LogWarning(metname) << "L2MuonSeedGenerator:   this should make sense only within MC tests" << endl;
0117       // FIXME! Temporary to handle the MC input
0118       quality = 7;
0119       valid_charge = true;
0120     } else {
0121       quality = muonCand.quality();
0122       valid_charge = muonCand.charge_valid();
0123     }
0124
0125     float pt = (*it).pt();
0126     float eta = (*it).eta();
0127     float theta = 2 * atan(exp(-eta));
0128     float phi = (*it).phi();
0129     int charge = (*it).charge();
0130     // Set charge=0 for the time being if the valid charge bit is zero
0131     if (!valid_charge)
0132       charge = 0;
0133     bool barrel = !(*it).isForward();
0134
0135     // Get a better eta and charge from regional information
0136     // Phi has the same resolution in GMT than regionally, is not it?
0137     if (!(muonCand.empty())) {
0138       int idx = -1;
0139       vector<L1MuRegionalCand> rmc;
0140       if (!muonCand.isRPC()) {
0141         idx = muonCand.getDTCSCIndex();
0142         if (muonCand.isFwd())
0143           rmc = gmtrr.getCSCCands();
0144         else
0145           rmc = gmtrr.getDTBXCands();
0146       } else {
0147         idx = muonCand.getRPCIndex();
0148         if (muonCand.isFwd())
0149           rmc = gmtrr.getFwdRPCCands();
0150         else
0151           rmc = gmtrr.getBrlRPCCands();
0152       }
0153       if (idx >= 0) {
0154         eta = rmc[idx].etaValue();
0155         //phi = rmc[idx].phiValue();
0156         // Use this charge if the valid charge bit is zero
0157         if (!valid_charge)
0158           charge = rmc[idx].chargeValue();
0159       }
0160     }
0161
0162     if (pt < theL1MinPt || fabs(eta) > theL1MaxEta)
0163       continue;
0164
0165     LogTrace(metname) << "New L2 Muon Seed";
0166     LogTrace(metname) << "Pt = " << pt << " GeV/c";
0167     LogTrace(metname) << "eta = " << eta;
0168     LogTrace(metname) << "theta = " << theta << " rad";
0169     LogTrace(metname) << "phi = " << phi << " rad";
0170     LogTrace(metname) << "charge = " << charge;
0171     LogTrace(metname) << "In Barrel? = " << barrel;
0172
0173     if (quality <= theL1MinQuality)
0174       continue;
0175     LogTrace(metname) << "quality = " << quality;
0176
0177     // Update the services
0178     theService->update(iSetup);
0179
0180     const DetLayer* detLayer = nullptr;
0181     float radius = 0.;
0182
0183     CLHEP::Hep3Vector vec(0., 1., 0.);
0184     vec.setTheta(theta);
0185     vec.setPhi(phi);
0186
0187     // Get the det layer on which the state should be put
0188     if (barrel) {
0189       LogTrace(metname) << "The seed is in the barrel";
0190
0191       // MB2
0192       DetId id = DTChamberId(0, 2, 0);
0193       detLayer = theService->detLayerGeometry()->idToLayer(id);
0194       LogTrace(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0195
0196       const BoundSurface* sur = &(detLayer->surface());
0197       const BoundCylinder* bc = dynamic_cast<const BoundCylinder*>(sur);
0198
0199       radius = fabs(bc->radius() / sin(theta));
0200
0201       LogTrace(metname) << "radius " << radius;
0202
0203       if (pt < 3.5)
0204         pt = 3.5;
0205     } else {
0206       LogTrace(metname) << "The seed is in the endcap";
0207
0208       DetId id;
0209       // ME2
0210       if (theta < Geom::pi() / 2.)
0211         id = CSCDetId(1, 2, 0, 0, 0);
0212       else
0213         id = CSCDetId(2, 2, 0, 0, 0);
0214
0215       detLayer = theService->detLayerGeometry()->idToLayer(id);
0216       LogTrace(metname) << "L2 Layer: " << debug.dumpLayer(detLayer);
0217
0218       radius = fabs(detLayer->position().z() / cos(theta));
0219
0220       if (pt < 1.0)
0221         pt = 1.0;
0222     }
0223
0224     vec.setMag(radius);
0225
0226     GlobalPoint pos(vec.x(), vec.y(), vec.z());
0227
0228     GlobalVector mom(pt * cos(phi), pt * sin(phi), pt * cos(theta) / sin(theta));
0229
0230     GlobalTrajectoryParameters param(pos, mom, charge, &*theService->magneticField());
0231     AlgebraicSymMatrix55 mat;
0232
0233     mat[0][0] = (0.25 / pt) * (0.25 / pt);  // sigma^2(charge/abs_momentum)
0234     if (!barrel)
0235       mat[0][0] = (0.4 / pt) * (0.4 / pt);
0236
0237     //Assign q/pt = 0 +- 1/pt if charge has been declared invalid
0238     if (!valid_charge)
0239       mat[0][0] = (1. / pt) * (1. / pt);
0240
0241     mat[1][1] = 0.05 * 0.05;  // sigma^2(lambda)
0242     mat[2][2] = 0.2 * 0.2;    // sigma^2(phi)
0243     mat[3][3] = 20. * 20.;    // sigma^2(x_transverse))
0244     mat[4][4] = 20. * 20.;    // sigma^2(y_transverse))
0245
0246     CurvilinearTrajectoryError error(mat);
0247
0248     const FreeTrajectoryState state(param, error);
0249
0250     LogTrace(metname) << "Free trajectory State from the parameters";
0251     LogTrace(metname) << debug.dumpFTS(state);
0252
0253     // Propagate the state on the MB2/ME2 surface
0254     TrajectoryStateOnSurface tsos = theService->propagator(thePropagatorName)->propagate(state, detLayer->surface());
0255
0256     LogTrace(metname) << "State after the propagation on the layer";
0257     LogTrace(metname) << debug.dumpLayer(detLayer);
0258     LogTrace(metname) << debug.dumpFTS(state);
0259
0260     if (tsos.isValid()) {
0261       // Get the compatible dets on the layer
0262       std::vector<pair<const GeomDet*, TrajectoryStateOnSurface> > detsWithStates =
0263           detLayer->compatibleDets(tsos, *theService->propagator(thePropagatorName), *theEstimator);
0264       if (!detsWithStates.empty()) {
0265         TrajectoryStateOnSurface newTSOS = detsWithStates.front().second;
0266         const GeomDet* newTSOSDet = detsWithStates.front().first;
0267
0268         LogTrace(metname) << "Most compatible det";
0269         LogTrace(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
0270
0271         LogDebug(metname) << "L1 info: Det and State:";
0272         LogDebug(metname) << debug.dumpMuonId(newTSOSDet->geographicalId());
0273
0274         if (newTSOS.isValid()) {
0275           //LogDebug(metname) << "(x, y, z) = (" << newTSOS.globalPosition().x() << ", "
0276           //                  << newTSOS.globalPosition().y() << ", " << newTSOS.globalPosition().z() << ")";
0277           LogDebug(metname) << "pos: (r=" << newTSOS.globalPosition().mag()
0278                             << ", phi=" << newTSOS.globalPosition().phi() << ", eta=" << newTSOS.globalPosition().eta()
0279                             << ")";
0280           LogDebug(metname) << "mom: (q*pt=" << newTSOS.charge() * newTSOS.globalMomentum().perp()
0281                             << ", phi=" << newTSOS.globalMomentum().phi() << ", eta=" << newTSOS.globalMomentum().eta()
0282                             << ")";
0283
0284           //LogDebug(metname) << "State on it";
0285           //LogDebug(metname) << debug.dumpTSOS(newTSOS);
0286
0287           //PTrajectoryStateOnDet seedTSOS;
0288           edm::OwnVector<TrackingRecHit> container;
0289
0290           if (useOfflineSeed) {
0291             const TrajectorySeed* assoOffseed = associateOfflineSeedToL1(offlineSeedHandle, offlineSeedMap, newTSOS);
0292
0293             if (assoOffseed != nullptr) {
0294               PTrajectoryStateOnDet const& seedTSOS = assoOffseed->startingState();
0295               for (auto const& tsci : assoOffseed->recHits()) {
0296                 container.push_back(tsci);
0297               }
0298               output->push_back(
0299                   L2MuonTrajectorySeed(seedTSOS, container, alongMomentum, L1MuonParticleRef(muColl, l1ParticleIndex)));
0300             } else {
0301               if (useUnassociatedL1) {
0302                 // convert the TSOS into a PTSOD
0303                 PTrajectoryStateOnDet const& seedTSOS =
0304                     trajectoryStateTransform::persistentState(newTSOS, newTSOSDet->geographicalId().rawId());
0305                 output->push_back(L2MuonTrajectorySeed(
0306                     seedTSOS, container, alongMomentum, L1MuonParticleRef(muColl, l1ParticleIndex)));
0307               }
0308             }
0309           } else {
0310             // convert the TSOS into a PTSOD
0311             PTrajectoryStateOnDet const& seedTSOS =
0312                 trajectoryStateTransform::persistentState(newTSOS, newTSOSDet->geographicalId().rawId());
0313             output->push_back(
0314                 L2MuonTrajectorySeed(seedTSOS, container, alongMomentum, L1MuonParticleRef(muColl, l1ParticleIndex)));
0315           }
0316         }
0317       }
0318     }
0319   }
0320
0321   iEvent.put(std::move(output));
0322 }
0323
0324 // FIXME: does not resolve ambiguities yet!
0325 const TrajectorySeed* L2MuonSeedGenerator::associateOfflineSeedToL1(edm::Handle<edm::View<TrajectorySeed> >& offseeds,
0326                                                                     std::vector<int>& offseedMap,
0327                                                                     TrajectoryStateOnSurface& newTsos) {
0328   const std::string metlabel = "Muon|RecoMuon|L2MuonSeedGenerator";
0329   MuonPatternRecoDumper debugtmp;
0330
0331   edm::View<TrajectorySeed>::const_iterator offseed, endOffseed = offseeds->end();
0332   const TrajectorySeed* selOffseed = nullptr;
0333   double bestDr = 99999.;
0334   unsigned int nOffseed(0);
0335   int lastOffseed(-1);
0336
0337   for (offseed = offseeds->begin(); offseed != endOffseed; ++offseed, ++nOffseed) {
0338     if (offseedMap[nOffseed] != 0)
0339       continue;
0340     GlobalPoint glbPos = theService->trackingGeometry()
0341                              ->idToDet(offseed->startingState().detId())
0342                              ->surface()
0343                              .toGlobal(offseed->startingState().parameters().position());
0344     GlobalVector glbMom = theService->trackingGeometry()
0345                               ->idToDet(offseed->startingState().detId())
0346                               ->surface()
0347                               .toGlobal(offseed->startingState().parameters().momentum());
0348
0349     // Preliminary check
0350     double preDr = deltaR(newTsos.globalPosition().eta(), newTsos.globalPosition().phi(), glbPos.eta(), glbPos.phi());
0351     if (preDr > 1.0)
0352       continue;
0353
0354     const FreeTrajectoryState offseedFTS(
0355         glbPos, glbMom, offseed->startingState().parameters().charge(), &*theService->magneticField());
0356     TrajectoryStateOnSurface offseedTsos =
0357         theService->propagator(thePropagatorName)->propagate(offseedFTS, newTsos.surface());
0358     LogDebug(metlabel) << "Offline seed info: Det and State" << std::endl;
0359     LogDebug(metlabel) << debugtmp.dumpMuonId(offseed->startingState().detId()) << std::endl;
0360     //LogDebug(metlabel) << "(x, y, z) = (" << newTSOS.globalPosition().x() << ", "
0361     //                   << newTSOS.globalPosition().y() << ", " << newTSOS.globalPosition().z() << ")" << std::endl;
0362     LogDebug(metlabel) << "pos: (r=" << offseedFTS.position().mag() << ", phi=" << offseedFTS.position().phi()
0363                        << ", eta=" << offseedFTS.position().eta() << ")" << std::endl;
0364     LogDebug(metlabel) << "mom: (q*pt=" << offseedFTS.charge() * offseedFTS.momentum().perp()
0365                        << ", phi=" << offseedFTS.momentum().phi() << ", eta=" << offseedFTS.momentum().eta() << ")"
0366                        << std::endl
0367                        << std::endl;
0368     //LogDebug(metlabel) << debugtmp.dumpFTS(offseedFTS) << std::endl;
0369
0370     if (offseedTsos.isValid()) {
0371       LogDebug(metlabel) << "Offline seed info after propagation to L1 layer:" << std::endl;
0372       //LogDebug(metlabel) << "(x, y, z) = (" << offseedTsos.globalPosition().x() << ", "
0373       //                   << offseedTsos.globalPosition().y() << ", " << offseedTsos.globalPosition().z() << ")" << std::endl;
0374       LogDebug(metlabel) << "pos: (r=" << offseedTsos.globalPosition().mag()
0375                          << ", phi=" << offseedTsos.globalPosition().phi()
0376                          << ", eta=" << offseedTsos.globalPosition().eta() << ")" << std::endl;
0377       LogDebug(metlabel) << "mom: (q*pt=" << offseedTsos.charge() * offseedTsos.globalMomentum().perp()
0378                          << ", phi=" << offseedTsos.globalMomentum().phi()
0379                          << ", eta=" << offseedTsos.globalMomentum().eta() << ")" << std::endl
0380                          << std::endl;
0381       //LogDebug(metlabel) << debugtmp.dumpTSOS(offseedTsos) << std::endl;
0382       double newDr = deltaR(newTsos.globalPosition().eta(),
0383                             newTsos.globalPosition().phi(),
0384                             offseedTsos.globalPosition().eta(),
0385                             offseedTsos.globalPosition().phi());
0386       LogDebug(metlabel) << "   -- DR = " << newDr << std::endl;
0387       if (newDr < 0.3 && newDr < bestDr) {
0388         LogDebug(metlabel) << "          --> OK! " << newDr << std::endl << std::endl;
0389         selOffseed = &*offseed;
0390         bestDr = newDr;
0391         offseedMap[nOffseed] = 1;
0392         if (lastOffseed > -1)
0393           offseedMap[lastOffseed] = 0;
0394         lastOffseed = nOffseed;
0395       } else {
0396         LogDebug(metlabel) << "          --> Rejected. " << newDr << std::endl << std::endl;
0397       }
0398     } else {
0399       LogDebug(metlabel) << "Invalid offline seed TSOS after propagation!" << std::endl << std::endl;
0400     }
0401   }
0402
0403   return selOffseed;
0404 }