SpecialSeedGenerators/src/CosmicRegionalSeedGenerator.cc

0001 #include <iostream>
0002 #include <memory>
0003 #include <string>
0004
0005 #include "DataFormats/MuonReco/interface/Muon.h"
0006 #include "DataFormats/JetReco/interface/CaloJet.h"
0007 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0008
0009 #include "RecoTracker/SpecialSeedGenerators/interface/CosmicRegionalSeedGenerator.h"
0010
0011 #include <TMath.h>
0012
0013 using namespace std;
0014 using namespace trigger;
0015 using namespace reco;
0016 using namespace edm;
0017
0018 CosmicRegionalSeedGenerator::CosmicRegionalSeedGenerator(edm::ParameterSet const& conf, edm::ConsumesCollector&& iC)
0019     : magfieldToken_(iC.esConsumes()) {
0020   edm::LogInfo("CosmicRegionalSeedGenerator") << "Begin Run:: Constructing  CosmicRegionalSeedGenerator";
0021
0022   edm::ParameterSet regionPSet = conf.getParameter<edm::ParameterSet>("RegionPSet");
0023   ptMin_ = regionPSet.getParameter<double>("ptMin");
0024   rVertex_ = regionPSet.getParameter<double>("rVertex");
0025   zVertex_ = regionPSet.getParameter<double>("zVertex");
0026   deltaEta_ = regionPSet.getParameter<double>("deltaEtaRegion");
0027   deltaPhi_ = regionPSet.getParameter<double>("deltaPhiRegion");
0028
0029   edm::ParameterSet toolsPSet = conf.getParameter<edm::ParameterSet>("ToolsPSet");
0030   propagatorToken_ = iC.esConsumes(edm::ESInputTag("", toolsPSet.getParameter<std::string>("thePropagatorName")));
0031   regionBase_ = toolsPSet.getParameter<std::string>("regionBase");
0032
0033   edm::ParameterSet collectionsPSet = conf.getParameter<edm::ParameterSet>("CollectionsPSet");
0034   recoMuonsCollection_ = collectionsPSet.getParameter<edm::InputTag>("recoMuonsCollection");
0035   recoTrackMuonsCollection_ = collectionsPSet.getParameter<edm::InputTag>("recoTrackMuonsCollection");
0036   recoL2MuonsCollection_ = collectionsPSet.getParameter<edm::InputTag>("recoL2MuonsCollection");
0037
0038   edm::ParameterSet regionInJetsCheckPSet = conf.getParameter<edm::ParameterSet>("RegionInJetsCheckPSet");
0039   doJetsExclusionCheck_ = regionInJetsCheckPSet.getParameter<bool>("doJetsExclusionCheck");
0040   deltaRExclusionSize_ = regionInJetsCheckPSet.getParameter<double>("deltaRExclusionSize");
0041   jetsPtMin_ = regionInJetsCheckPSet.getParameter<double>("jetsPtMin");
0042   recoCaloJetsCollection_ = regionInJetsCheckPSet.getParameter<edm::InputTag>("recoCaloJetsCollection");
0043   recoCaloJetsToken_ = iC.consumes<reco::CaloJetCollection>(recoCaloJetsCollection_);
0044   recoMuonsToken_ = iC.consumes<reco::MuonCollection>(recoMuonsCollection_);
0045   recoTrackMuonsToken_ = iC.consumes<reco::TrackCollection>(recoTrackMuonsCollection_);
0046   recoL2MuonsToken_ = iC.consumes<reco::RecoChargedCandidateCollection>(recoL2MuonsCollection_);
0047   measurementTrackerEventToken_ = iC.consumes<MeasurementTrackerEvent>(edm::InputTag("MeasurementTrackerEvent"));
0048
0049   edm::LogInfo("CosmicRegionalSeedGenerator") << "Reco muons collection: " << recoMuonsCollection_ << "\n"
0050                                               << "Reco tracks muons collection: " << recoTrackMuonsCollection_ << "\n"
0051                                               << "Reco L2 muons collection: " << recoL2MuonsCollection_;
0052 }
0053
0054 std::vector<std::unique_ptr<TrackingRegion>> CosmicRegionalSeedGenerator::regions(const edm::Event& event,
0055                                                                                   const edm::EventSetup& es) const {
0056   std::vector<std::unique_ptr<TrackingRegion>> result;
0057
0058   const MeasurementTrackerEvent* measurementTracker = nullptr;
0059   if (!measurementTrackerEventToken_.isUninitialized()) {
0060     edm::Handle<MeasurementTrackerEvent> hmte;
0061     event.getByToken(measurementTrackerEventToken_, hmte);
0062     measurementTracker = hmte.product();
0063   }
0064
0065   //get the propagator
0066   auto const& propagator = es.getData(propagatorToken_);
0067   auto const& magfield = es.getData(magfieldToken_);
0068
0069   //________________________________________
0070   //
0071   //Seeding on Sta muon (MC && Datas)
0072   //________________________________________
0073
0074   if (regionBase_ == "seedOnStaMuon" || regionBase_.empty()) {
0075     LogDebug("CosmicRegionalSeedGenerator") << "Seeding on stand alone muons ";
0076
0077     //get collections
0078     //+++++++++++++++
0079
0080     //get the muon collection
0081     edm::Handle<reco::MuonCollection> muonsHandle;
0082     event.getByToken(recoMuonsToken_, muonsHandle);
0083     if (!muonsHandle.isValid()) {
0084       edm::LogError("CollectionNotFound") << "Error::No reco muons collection (" << recoMuonsCollection_
0085                                           << ") in the event - Please verify the name of the muon collection";
0086       return result;
0087     }
0088
0089     LogDebug("CosmicRegionalSeedGenerator") << "Muons collection size = " << muonsHandle->size();
0090
0091     //get the jet collection
0092     edm::Handle<CaloJetCollection> caloJetsHandle;
0093     event.getByToken(recoCaloJetsToken_, caloJetsHandle);
0094
0095     //definition of the region
0096     //+++++++++++++++++++++++++
0097
0098     for (reco::MuonCollection::const_iterator staMuon = muonsHandle->begin(); staMuon != muonsHandle->end();
0099          ++staMuon) {
0100       //select sta muons
0101       if (!staMuon->isStandAloneMuon()) {
0102         LogDebug("CosmicRegionalSeedGenerator") << "This muon is not a stand alone muon";
0103         continue;
0104       }
0105
0106       //bit 25 as a coverage -1.4 < eta < 1.4
0107       if (abs(staMuon->standAloneMuon()->eta()) > 1.5)
0108         continue;
0109
0110       //debug
0111       LogDebug("CosmicRegionalSeedGenerator") << "Muon stand alone found in the collection - in muons chambers: \n "
0112                                               << "Position = " << staMuon->standAloneMuon()->outerPosition() << "\n "
0113                                               << "Momentum = " << staMuon->standAloneMuon()->outerMomentum() << "\n "
0114                                               << "Eta = " << staMuon->standAloneMuon()->eta() << "\n "
0115                                               << "Phi = " << staMuon->standAloneMuon()->phi();
0116
0117       //initial position, momentum, charge
0118
0119       GlobalPoint initialRegionPosition(staMuon->standAloneMuon()->referencePoint().x(),
0120                                         staMuon->standAloneMuon()->referencePoint().y(),
0121                                         staMuon->standAloneMuon()->referencePoint().z());
0122       GlobalVector initialRegionMomentum(staMuon->standAloneMuon()->momentum().x(),
0123                                          staMuon->standAloneMuon()->momentum().y(),
0124                                          staMuon->standAloneMuon()->momentum().z());
0125       int charge = (int)staMuon->standAloneMuon()->charge();
0126
0127       LogDebug("CosmicRegionalSeedGenerator") << "Initial region - Reference point of the sta muon: \n "
0128                                               << "Position = " << initialRegionPosition << "\n "
0129                                               << "Momentum = " << initialRegionMomentum << "\n "
0130                                               << "Eta = " << initialRegionPosition.eta() << "\n "
0131                                               << "Phi = " << initialRegionPosition.phi() << "\n "
0132                                               << "Charge = " << charge;
0133
0134       //propagation on the last layers of TOB
0135       if (staMuon->standAloneMuon()->outerPosition().y() > 0)
0136         initialRegionMomentum *= -1;
0137       GlobalTrajectoryParameters glb_parameters(
0138           initialRegionPosition, initialRegionMomentum, charge, propagator.magneticField());
0139       FreeTrajectoryState fts(glb_parameters);
0140       StateOnTrackerBound onBounds(&propagator);
0141       TrajectoryStateOnSurface outer = onBounds(fts);
0142
0143       if (!outer.isValid()) {
0144         //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
0145         LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid";
0146         continue;
0147       }
0148
0149       //final position & momentum
0150       GlobalPoint regionPosition = outer.globalPosition();
0151       GlobalVector regionMom = outer.globalMomentum();
0152
0153       LogDebug("CosmicRegionalSeedGenerator")
0154           << "Region after propagation: \n "
0155           << "Position = " << outer.globalPosition() << "\n "
0156           << "Momentum = " << outer.globalMomentum() << "\n "
0157           << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
0158           << "Eta = " << outer.globalPosition().eta() << "\n "
0159           << "Phi = " << outer.globalPosition().phi();
0160
0161       //step back
0162       double stepBack = 1;
0163       GlobalPoint center = regionPosition + stepBack * regionMom.unit();
0164       GlobalVector v = stepBack * regionMom.unit();
0165       LogDebug("CosmicRegionalSeedGenerator") << "Step back vector =  " << v << "\n";
0166
0167       //exclude region built in jets
0168       if (doJetsExclusionCheck_) {
0169         double delta_R_min = 1000.;
0170         for (CaloJetCollection::const_iterator jet = caloJetsHandle->begin(); jet != caloJetsHandle->end(); jet++) {
0171           if (jet->pt() < jetsPtMin_)
0172             continue;
0173
0174           double deta = center.eta() - jet->eta();
0175           double dphi = fabs(center.phi() - jet->phi());
0176           if (dphi > TMath::Pi())
0177             dphi = 2 * TMath::Pi() - dphi;
0178
0179           double delta_R = sqrt(deta * deta + dphi * dphi);
0180           if (delta_R < delta_R_min)
0181             delta_R_min = delta_R;
0182
0183         }  //end loop on jets
0184
0185         if (delta_R_min < deltaRExclusionSize_) {
0186           LogDebug("CosmicRegionalSeedGenerator") << "Region built too close from a jet";
0187           continue;
0188         }
0189       }  //end if doJetsExclusionCheck
0190
0191       //definition of the region
0192
0193       result.push_back(std::make_unique<CosmicTrackingRegion>(
0194           (-1) * regionMom, center, ptMin_, rVertex_, zVertex_, deltaEta_, deltaPhi_, magfield, measurementTracker));
0195
0196       LogDebug("CosmicRegionalSeedGenerator")
0197           << "Final CosmicTrackingRegion \n "
0198           << "Position = " << center << "\n "
0199           << "Direction = " << result.back()->direction() << "\n "
0200           << "Distance from the region on the layer = " << (regionPosition - center).mag() << "\n "
0201           << "Eta = " << center.eta() << "\n "
0202           << "Phi = " << center.phi();
0203
0204     }  //end loop on muons
0205
0206   }  //end if SeedOnStaMuon
0207
0208   //________________________________________
0209   //
0210   //Seeding on cosmic muons (MC && Datas)
0211   //________________________________________
0212
0213   if (regionBase_ == "seedOnCosmicMuon") {
0214     LogDebug("CosmicRegionalSeedGenerator") << "Seeding on cosmic muons tracks";
0215
0216     //get collections
0217     //+++++++++++++++
0218
0219     //get the muon collection
0220     edm::Handle<reco::TrackCollection> cosmicMuonsHandle;
0221     event.getByToken(recoTrackMuonsToken_, cosmicMuonsHandle);
0222     if (!cosmicMuonsHandle.isValid()) {
0223       edm::LogError("CollectionNotFound")
0224           << "Error::No cosmic muons collection (" << recoTrackMuonsCollection_
0225           << ") in the event - Please verify the name of the muon reco track collection";
0226       return result;
0227     }
0228
0229     LogDebug("CosmicRegionalSeedGenerator") << "Cosmic muons tracks collection size = " << cosmicMuonsHandle->size();
0230
0231     //get the jet collection
0232     edm::Handle<CaloJetCollection> caloJetsHandle;
0233     event.getByToken(recoCaloJetsToken_, caloJetsHandle);
0234
0235     //definition of the region
0236     //+++++++++++++++++++++++++
0237
0238     for (reco::TrackCollection::const_iterator cosmicMuon = cosmicMuonsHandle->begin();
0239          cosmicMuon != cosmicMuonsHandle->end();
0240          ++cosmicMuon) {
0241       //bit 25 as a coverage -1.4 < eta < 1.4
0242       if (abs(cosmicMuon->eta()) > 1.5)
0243         continue;
0244
0245       //initial position, momentum, charge
0246       GlobalPoint initialRegionPosition(
0247           cosmicMuon->referencePoint().x(), cosmicMuon->referencePoint().y(), cosmicMuon->referencePoint().z());
0248       GlobalVector initialRegionMomentum(
0249           cosmicMuon->momentum().x(), cosmicMuon->momentum().y(), cosmicMuon->momentum().z());
0250       int charge = (int)cosmicMuon->charge();
0251
0252       LogDebug("CosmicRegionalSeedGenerator") << "Position and momentum of the muon track in the muon chambers: \n "
0253                                               << "x = " << cosmicMuon->outerPosition().x() << "\n "
0254                                               << "y = " << cosmicMuon->outerPosition().y() << "\n "
0255                                               << "y = " << cosmicMuon->pt() << "\n "
0256                                               << "Initial region - Reference point of the cosmic muon track: \n "
0257                                               << "Position = " << initialRegionPosition << "\n "
0258                                               << "Momentum = " << initialRegionMomentum << "\n "
0259                                               << "Eta = " << initialRegionPosition.eta() << "\n "
0260                                               << "Phi = " << initialRegionPosition.phi() << "\n "
0261                                               << "Charge = " << charge;
0262
0263       //propagation on the last layers of TOB
0264       if (cosmicMuon->outerPosition().y() > 0 && cosmicMuon->momentum().y() < 0)
0265         initialRegionMomentum *= -1;
0266       GlobalTrajectoryParameters glb_parameters(
0267           initialRegionPosition, initialRegionMomentum, charge, propagator.magneticField());
0268       FreeTrajectoryState fts(glb_parameters);
0269       StateOnTrackerBound onBounds(&propagator);
0270       TrajectoryStateOnSurface outer = onBounds(fts);
0271
0272       if (!outer.isValid()) {
0273         //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
0274         LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid";
0275         continue;
0276       }
0277
0278       //final position & momentum
0279       GlobalPoint regionPosition = outer.globalPosition();
0280       GlobalVector regionMom = outer.globalMomentum();
0281
0282       LogDebug("CosmicRegionalSeedGenerator")
0283           << "Region after propagation: \n "
0284           << "Position = " << outer.globalPosition() << "\n "
0285           << "Momentum = " << outer.globalMomentum() << "\n "
0286           << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
0287           << "Eta = " << outer.globalPosition().eta() << "\n "
0288           << "Phi = " << outer.globalPosition().phi();
0289
0290       //step back
0291       double stepBack = 1;
0292       GlobalPoint center = regionPosition + stepBack * regionMom.unit();
0293       GlobalVector v = stepBack * regionMom.unit();
0294       LogDebug("CosmicRegionalSeedGenerator") << "Step back vector =  " << v << "\n";
0295
0296       //exclude region built in jets
0297       if (doJetsExclusionCheck_) {
0298         double delta_R_min = 1000.;
0299         for (CaloJetCollection::const_iterator jet = caloJetsHandle->begin(); jet != caloJetsHandle->end(); jet++) {
0300           if (jet->pt() < jetsPtMin_)
0301             continue;
0302
0303           double deta = center.eta() - jet->eta();
0304           double dphi = fabs(center.phi() - jet->phi());
0305           if (dphi > TMath::Pi())
0306             dphi = 2 * TMath::Pi() - dphi;
0307
0308           double delta_R = sqrt(deta * deta + dphi * dphi);
0309           if (delta_R < delta_R_min)
0310             delta_R_min = delta_R;
0311
0312         }  //end loop on jets
0313
0314         if (delta_R_min < deltaRExclusionSize_) {
0315           LogDebug("CosmicRegionalSeedGenerator") << "Region built too close from a jet";
0316           continue;
0317         }
0318       }  // end if doJetsExclusionCheck
0319
0320       //definition of the region
0321       result.push_back(std::make_unique<CosmicTrackingRegion>(
0322           (-1) * regionMom, center, ptMin_, rVertex_, zVertex_, deltaEta_, deltaPhi_, magfield, measurementTracker));
0323
0324       LogDebug("CosmicRegionalSeedGenerator")
0325           << "Final CosmicTrackingRegion \n "
0326           << "Position = " << center << "\n "
0327           << "Direction = " << result.back()->direction() << "\n "
0328           << "Distance from the region on the layer = " << (regionPosition - center).mag() << "\n "
0329           << "Eta = " << center.eta() << "\n "
0330           << "Phi = " << center.phi();
0331
0332     }  //end loop on muons
0333
0334   }  //end if SeedOnCosmicMuon
0335
0336   //________________________________________
0337   //
0338   //Seeding on L2 muons (MC && Datas)
0339   //________________________________________
0340
0341   if (regionBase_ == "seedOnL2Muon") {
0342     LogDebug("CosmicRegionalSeedGenerator") << "Seeding on L2 muons";
0343
0344     //get collections
0345     //+++++++++++++++
0346
0347     //get the muon collection
0348     edm::Handle<reco::RecoChargedCandidateCollection> L2MuonsHandle;
0349     event.getByToken(recoL2MuonsToken_, L2MuonsHandle);
0350
0351     if (!L2MuonsHandle.isValid()) {
0352       edm::LogError("CollectionNotFound") << "Error::No L2 muons collection (" << recoL2MuonsCollection_
0353                                           << ") in the event - Please verify the name of the L2 muon collection";
0354       return result;
0355     }
0356
0357     LogDebug("CosmicRegionalSeedGenerator") << "L2 muons collection size = " << L2MuonsHandle->size();
0358
0359     //definition of the region
0360     //+++++++++++++++++++++++++
0361
0362     for (reco::RecoChargedCandidateCollection::const_iterator L2Muon = L2MuonsHandle->begin();
0363          L2Muon != L2MuonsHandle->end();
0364          ++L2Muon) {
0365       reco::TrackRef tkL2Muon = L2Muon->get<reco::TrackRef>();
0366
0367       //bit 25 as a coverage -1.4 < eta < 1.4
0368       if (abs(tkL2Muon->eta()) > 1.5)
0369         continue;
0370
0371       //initial position, momentum, charge
0372       GlobalPoint initialRegionPosition(
0373           tkL2Muon->referencePoint().x(), tkL2Muon->referencePoint().y(), tkL2Muon->referencePoint().z());
0374       GlobalVector initialRegionMomentum(tkL2Muon->momentum().x(), tkL2Muon->momentum().y(), tkL2Muon->momentum().z());
0375       int charge = (int)tkL2Muon->charge();
0376
0377       LogDebug("CosmicRegionalSeedGenerator") << "Position and momentum of the L2 muon track in the muon chambers: \n "
0378                                               << "x = " << tkL2Muon->outerPosition().x() << "\n "
0379                                               << "y = " << tkL2Muon->outerPosition().y() << "\n "
0380                                               << "y = " << tkL2Muon->pt() << "\n "
0381                                               << "Initial region - Reference point of the L2 muon track: \n "
0382                                               << "Position = " << initialRegionPosition << "\n "
0383                                               << "Momentum = " << initialRegionMomentum << "\n "
0384                                               << "Eta = " << initialRegionPosition.eta() << "\n "
0385                                               << "Phi = " << initialRegionPosition.phi() << "\n "
0386                                               << "Charge = " << charge;
0387
0388       //seeding only in the bottom
0389       if (tkL2Muon->outerPosition().y() > 0) {
0390         LogDebug("CosmicRegionalSeedGenerator") << "L2 muon in the TOP --- Region not created";
0391         return result;
0392       }
0393
0394       GlobalTrajectoryParameters glb_parameters(
0395           initialRegionPosition, initialRegionMomentum, charge, propagator.magneticField());
0396       FreeTrajectoryState fts(glb_parameters);
0397       StateOnTrackerBound onBounds(&propagator);
0398       TrajectoryStateOnSurface outer = onBounds(fts);
0399
0400       if (!outer.isValid()) {
0401         //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
0402         LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid";
0403         continue;
0404       }
0405
0406       //final position & momentum
0407       GlobalPoint regionPosition = outer.globalPosition();
0408       GlobalVector regionMom = outer.globalMomentum();
0409
0410       LogDebug("CosmicRegionalSeedGenerator")
0411           << "Region after propagation: \n "
0412           << "Position = " << outer.globalPosition() << "\n "
0413           << "Momentum = " << outer.globalMomentum() << "\n "
0414           << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
0415           << "Eta = " << outer.globalPosition().eta() << "\n "
0416           << "Phi = " << outer.globalPosition().phi();
0417
0418       //step back
0419       double stepBack = 1;
0420       GlobalPoint center = regionPosition + stepBack * regionMom.unit();
0421       GlobalVector v = stepBack * regionMom.unit();
0422       LogDebug("CosmicRegionalSeedGenerator") << "Step back vector =  " << v << "\n";
0423
0424       //definition of the region
0425       result.push_back(std::make_unique<CosmicTrackingRegion>(
0426           (-1) * regionMom, center, ptMin_, rVertex_, zVertex_, deltaEta_, deltaPhi_, magfield, measurementTracker));
0427
0428       LogDebug("CosmicRegionalSeedGenerator")
0429           << "Final L2TrackingRegion \n "
0430           << "Position = " << center << "\n "
0431           << "Direction = " << result.back()->direction() << "\n "
0432           << "Distance from the region on the layer = " << (regionPosition - center).mag() << "\n "
0433           << "Eta = " << center.eta() << "\n "
0434           << "Phi = " << center.phi();
0435
0436     }  //end loop on muons
0437
0438   }  //end if SeedOnL2Muon
0439
0440   return result;
0441 }