Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-01 22:41:01

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     int nmuons = 0;
0099     for (reco::MuonCollection::const_iterator staMuon = muonsHandle->begin(); staMuon != muonsHandle->end();
0100          ++staMuon) {
0101       //select sta muons
0102       if (!staMuon->isStandAloneMuon()) {
0103         LogDebug("CosmicRegionalSeedGenerator") << "This muon is not a stand alone muon";
0104         continue;
0105       }
0106 
0107       //bit 25 as a coverage -1.4 < eta < 1.4
0108       if (abs(staMuon->standAloneMuon()->eta()) > 1.5)
0109         continue;
0110 
0111       //debug
0112       nmuons++;
0113       LogDebug("CosmicRegionalSeedGenerator") << "Muon stand alone found in the collection - in muons chambers: \n "
0114                                               << "Position = " << staMuon->standAloneMuon()->outerPosition() << "\n "
0115                                               << "Momentum = " << staMuon->standAloneMuon()->outerMomentum() << "\n "
0116                                               << "Eta = " << staMuon->standAloneMuon()->eta() << "\n "
0117                                               << "Phi = " << staMuon->standAloneMuon()->phi();
0118 
0119       //initial position, momentum, charge
0120 
0121       GlobalPoint initialRegionPosition(staMuon->standAloneMuon()->referencePoint().x(),
0122                                         staMuon->standAloneMuon()->referencePoint().y(),
0123                                         staMuon->standAloneMuon()->referencePoint().z());
0124       GlobalVector initialRegionMomentum(staMuon->standAloneMuon()->momentum().x(),
0125                                          staMuon->standAloneMuon()->momentum().y(),
0126                                          staMuon->standAloneMuon()->momentum().z());
0127       int charge = (int)staMuon->standAloneMuon()->charge();
0128 
0129       LogDebug("CosmicRegionalSeedGenerator") << "Initial region - Reference point of the sta muon: \n "
0130                                               << "Position = " << initialRegionPosition << "\n "
0131                                               << "Momentum = " << initialRegionMomentum << "\n "
0132                                               << "Eta = " << initialRegionPosition.eta() << "\n "
0133                                               << "Phi = " << initialRegionPosition.phi() << "\n "
0134                                               << "Charge = " << charge;
0135 
0136       //propagation on the last layers of TOB
0137       if (staMuon->standAloneMuon()->outerPosition().y() > 0)
0138         initialRegionMomentum *= -1;
0139       GlobalTrajectoryParameters glb_parameters(
0140           initialRegionPosition, initialRegionMomentum, charge, propagator.magneticField());
0141       FreeTrajectoryState fts(glb_parameters);
0142       StateOnTrackerBound onBounds(&propagator);
0143       TrajectoryStateOnSurface outer = onBounds(fts);
0144 
0145       if (!outer.isValid()) {
0146         //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
0147         LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid";
0148         continue;
0149       }
0150 
0151       //final position & momentum
0152       GlobalPoint regionPosition = outer.globalPosition();
0153       GlobalVector regionMom = outer.globalMomentum();
0154 
0155       LogDebug("CosmicRegionalSeedGenerator")
0156           << "Region after propagation: \n "
0157           << "Position = " << outer.globalPosition() << "\n "
0158           << "Momentum = " << outer.globalMomentum() << "\n "
0159           << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
0160           << "Eta = " << outer.globalPosition().eta() << "\n "
0161           << "Phi = " << outer.globalPosition().phi();
0162 
0163       //step back
0164       double stepBack = 1;
0165       GlobalPoint center = regionPosition + stepBack * regionMom.unit();
0166       GlobalVector v = stepBack * regionMom.unit();
0167       LogDebug("CosmicRegionalSeedGenerator") << "Step back vector =  " << v << "\n";
0168 
0169       //exclude region built in jets
0170       if (doJetsExclusionCheck_) {
0171         double delta_R_min = 1000.;
0172         for (CaloJetCollection::const_iterator jet = caloJetsHandle->begin(); jet != caloJetsHandle->end(); jet++) {
0173           if (jet->pt() < jetsPtMin_)
0174             continue;
0175 
0176           double deta = center.eta() - jet->eta();
0177           double dphi = fabs(center.phi() - jet->phi());
0178           if (dphi > TMath::Pi())
0179             dphi = 2 * TMath::Pi() - dphi;
0180 
0181           double delta_R = sqrt(deta * deta + dphi * dphi);
0182           if (delta_R < delta_R_min)
0183             delta_R_min = delta_R;
0184 
0185         }  //end loop on jets
0186 
0187         if (delta_R_min < deltaRExclusionSize_) {
0188           LogDebug("CosmicRegionalSeedGenerator") << "Region built too close from a jet";
0189           continue;
0190         }
0191       }  //end if doJetsExclusionCheck
0192 
0193       //definition of the region
0194 
0195       result.push_back(std::make_unique<CosmicTrackingRegion>(
0196           (-1) * regionMom, center, ptMin_, rVertex_, zVertex_, deltaEta_, deltaPhi_, magfield, measurementTracker));
0197 
0198       LogDebug("CosmicRegionalSeedGenerator")
0199           << "Final CosmicTrackingRegion \n "
0200           << "Position = " << center << "\n "
0201           << "Direction = " << result.back()->direction() << "\n "
0202           << "Distance from the region on the layer = " << (regionPosition - center).mag() << "\n "
0203           << "Eta = " << center.eta() << "\n "
0204           << "Phi = " << center.phi();
0205 
0206     }  //end loop on muons
0207 
0208   }  //end if SeedOnStaMuon
0209 
0210   //________________________________________
0211   //
0212   //Seeding on cosmic muons (MC && Datas)
0213   //________________________________________
0214 
0215   if (regionBase_ == "seedOnCosmicMuon") {
0216     LogDebug("CosmicRegionalSeedGenerator") << "Seeding on cosmic muons tracks";
0217 
0218     //get collections
0219     //+++++++++++++++
0220 
0221     //get the muon collection
0222     edm::Handle<reco::TrackCollection> cosmicMuonsHandle;
0223     event.getByToken(recoTrackMuonsToken_, cosmicMuonsHandle);
0224     if (!cosmicMuonsHandle.isValid()) {
0225       edm::LogError("CollectionNotFound")
0226           << "Error::No cosmic muons collection (" << recoTrackMuonsCollection_
0227           << ") in the event - Please verify the name of the muon reco track collection";
0228       return result;
0229     }
0230 
0231     LogDebug("CosmicRegionalSeedGenerator") << "Cosmic muons tracks collection size = " << cosmicMuonsHandle->size();
0232 
0233     //get the jet collection
0234     edm::Handle<CaloJetCollection> caloJetsHandle;
0235     event.getByToken(recoCaloJetsToken_, caloJetsHandle);
0236 
0237     //definition of the region
0238     //+++++++++++++++++++++++++
0239 
0240     int nmuons = 0;
0241     for (reco::TrackCollection::const_iterator cosmicMuon = cosmicMuonsHandle->begin();
0242          cosmicMuon != cosmicMuonsHandle->end();
0243          ++cosmicMuon) {
0244       //bit 25 as a coverage -1.4 < eta < 1.4
0245       if (abs(cosmicMuon->eta()) > 1.5)
0246         continue;
0247 
0248       nmuons++;
0249 
0250       //initial position, momentum, charge
0251       GlobalPoint initialRegionPosition(
0252           cosmicMuon->referencePoint().x(), cosmicMuon->referencePoint().y(), cosmicMuon->referencePoint().z());
0253       GlobalVector initialRegionMomentum(
0254           cosmicMuon->momentum().x(), cosmicMuon->momentum().y(), cosmicMuon->momentum().z());
0255       int charge = (int)cosmicMuon->charge();
0256 
0257       LogDebug("CosmicRegionalSeedGenerator") << "Position and momentum of the muon track in the muon chambers: \n "
0258                                               << "x = " << cosmicMuon->outerPosition().x() << "\n "
0259                                               << "y = " << cosmicMuon->outerPosition().y() << "\n "
0260                                               << "y = " << cosmicMuon->pt() << "\n "
0261                                               << "Initial region - Reference point of the cosmic muon track: \n "
0262                                               << "Position = " << initialRegionPosition << "\n "
0263                                               << "Momentum = " << initialRegionMomentum << "\n "
0264                                               << "Eta = " << initialRegionPosition.eta() << "\n "
0265                                               << "Phi = " << initialRegionPosition.phi() << "\n "
0266                                               << "Charge = " << charge;
0267 
0268       //propagation on the last layers of TOB
0269       if (cosmicMuon->outerPosition().y() > 0 && cosmicMuon->momentum().y() < 0)
0270         initialRegionMomentum *= -1;
0271       GlobalTrajectoryParameters glb_parameters(
0272           initialRegionPosition, initialRegionMomentum, charge, propagator.magneticField());
0273       FreeTrajectoryState fts(glb_parameters);
0274       StateOnTrackerBound onBounds(&propagator);
0275       TrajectoryStateOnSurface outer = onBounds(fts);
0276 
0277       if (!outer.isValid()) {
0278         //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
0279         LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid";
0280         continue;
0281       }
0282 
0283       //final position & momentum
0284       GlobalPoint regionPosition = outer.globalPosition();
0285       GlobalVector regionMom = outer.globalMomentum();
0286 
0287       LogDebug("CosmicRegionalSeedGenerator")
0288           << "Region after propagation: \n "
0289           << "Position = " << outer.globalPosition() << "\n "
0290           << "Momentum = " << outer.globalMomentum() << "\n "
0291           << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
0292           << "Eta = " << outer.globalPosition().eta() << "\n "
0293           << "Phi = " << outer.globalPosition().phi();
0294 
0295       //step back
0296       double stepBack = 1;
0297       GlobalPoint center = regionPosition + stepBack * regionMom.unit();
0298       GlobalVector v = stepBack * regionMom.unit();
0299       LogDebug("CosmicRegionalSeedGenerator") << "Step back vector =  " << v << "\n";
0300 
0301       //exclude region built in jets
0302       if (doJetsExclusionCheck_) {
0303         double delta_R_min = 1000.;
0304         for (CaloJetCollection::const_iterator jet = caloJetsHandle->begin(); jet != caloJetsHandle->end(); jet++) {
0305           if (jet->pt() < jetsPtMin_)
0306             continue;
0307 
0308           double deta = center.eta() - jet->eta();
0309           double dphi = fabs(center.phi() - jet->phi());
0310           if (dphi > TMath::Pi())
0311             dphi = 2 * TMath::Pi() - dphi;
0312 
0313           double delta_R = sqrt(deta * deta + dphi * dphi);
0314           if (delta_R < delta_R_min)
0315             delta_R_min = delta_R;
0316 
0317         }  //end loop on jets
0318 
0319         if (delta_R_min < deltaRExclusionSize_) {
0320           LogDebug("CosmicRegionalSeedGenerator") << "Region built too close from a jet";
0321           continue;
0322         }
0323       }  // end if doJetsExclusionCheck
0324 
0325       //definition of the region
0326       result.push_back(std::make_unique<CosmicTrackingRegion>(
0327           (-1) * regionMom, center, ptMin_, rVertex_, zVertex_, deltaEta_, deltaPhi_, magfield, measurementTracker));
0328 
0329       LogDebug("CosmicRegionalSeedGenerator")
0330           << "Final CosmicTrackingRegion \n "
0331           << "Position = " << center << "\n "
0332           << "Direction = " << result.back()->direction() << "\n "
0333           << "Distance from the region on the layer = " << (regionPosition - center).mag() << "\n "
0334           << "Eta = " << center.eta() << "\n "
0335           << "Phi = " << center.phi();
0336 
0337     }  //end loop on muons
0338 
0339   }  //end if SeedOnCosmicMuon
0340 
0341   //________________________________________
0342   //
0343   //Seeding on L2 muons (MC && Datas)
0344   //________________________________________
0345 
0346   if (regionBase_ == "seedOnL2Muon") {
0347     LogDebug("CosmicRegionalSeedGenerator") << "Seeding on L2 muons";
0348 
0349     //get collections
0350     //+++++++++++++++
0351 
0352     //get the muon collection
0353     edm::Handle<reco::RecoChargedCandidateCollection> L2MuonsHandle;
0354     event.getByToken(recoL2MuonsToken_, L2MuonsHandle);
0355 
0356     if (!L2MuonsHandle.isValid()) {
0357       edm::LogError("CollectionNotFound") << "Error::No L2 muons collection (" << recoL2MuonsCollection_
0358                                           << ") in the event - Please verify the name of the L2 muon collection";
0359       return result;
0360     }
0361 
0362     LogDebug("CosmicRegionalSeedGenerator") << "L2 muons collection size = " << L2MuonsHandle->size();
0363 
0364     //definition of the region
0365     //+++++++++++++++++++++++++
0366 
0367     int nmuons = 0;
0368     for (reco::RecoChargedCandidateCollection::const_iterator L2Muon = L2MuonsHandle->begin();
0369          L2Muon != L2MuonsHandle->end();
0370          ++L2Muon) {
0371       reco::TrackRef tkL2Muon = L2Muon->get<reco::TrackRef>();
0372 
0373       //bit 25 as a coverage -1.4 < eta < 1.4
0374       if (abs(tkL2Muon->eta()) > 1.5)
0375         continue;
0376 
0377       nmuons++;
0378 
0379       //initial position, momentum, charge
0380       GlobalPoint initialRegionPosition(
0381           tkL2Muon->referencePoint().x(), tkL2Muon->referencePoint().y(), tkL2Muon->referencePoint().z());
0382       GlobalVector initialRegionMomentum(tkL2Muon->momentum().x(), tkL2Muon->momentum().y(), tkL2Muon->momentum().z());
0383       int charge = (int)tkL2Muon->charge();
0384 
0385       LogDebug("CosmicRegionalSeedGenerator") << "Position and momentum of the L2 muon track in the muon chambers: \n "
0386                                               << "x = " << tkL2Muon->outerPosition().x() << "\n "
0387                                               << "y = " << tkL2Muon->outerPosition().y() << "\n "
0388                                               << "y = " << tkL2Muon->pt() << "\n "
0389                                               << "Initial region - Reference point of the L2 muon track: \n "
0390                                               << "Position = " << initialRegionPosition << "\n "
0391                                               << "Momentum = " << initialRegionMomentum << "\n "
0392                                               << "Eta = " << initialRegionPosition.eta() << "\n "
0393                                               << "Phi = " << initialRegionPosition.phi() << "\n "
0394                                               << "Charge = " << charge;
0395 
0396       //seeding only in the bottom
0397       if (tkL2Muon->outerPosition().y() > 0) {
0398         LogDebug("CosmicRegionalSeedGenerator") << "L2 muon in the TOP --- Region not created";
0399         return result;
0400       }
0401 
0402       GlobalTrajectoryParameters glb_parameters(
0403           initialRegionPosition, initialRegionMomentum, charge, propagator.magneticField());
0404       FreeTrajectoryState fts(glb_parameters);
0405       StateOnTrackerBound onBounds(&propagator);
0406       TrajectoryStateOnSurface outer = onBounds(fts);
0407 
0408       if (!outer.isValid()) {
0409         //edm::LogError("FailedPropagation") << "Trajectory state on surface not valid" ;
0410         LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid";
0411         continue;
0412       }
0413 
0414       //final position & momentum
0415       GlobalPoint regionPosition = outer.globalPosition();
0416       GlobalVector regionMom = outer.globalMomentum();
0417 
0418       LogDebug("CosmicRegionalSeedGenerator")
0419           << "Region after propagation: \n "
0420           << "Position = " << outer.globalPosition() << "\n "
0421           << "Momentum = " << outer.globalMomentum() << "\n "
0422           << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
0423           << "Eta = " << outer.globalPosition().eta() << "\n "
0424           << "Phi = " << outer.globalPosition().phi();
0425 
0426       //step back
0427       double stepBack = 1;
0428       GlobalPoint center = regionPosition + stepBack * regionMom.unit();
0429       GlobalVector v = stepBack * regionMom.unit();
0430       LogDebug("CosmicRegionalSeedGenerator") << "Step back vector =  " << v << "\n";
0431 
0432       //definition of the region
0433       result.push_back(std::make_unique<CosmicTrackingRegion>(
0434           (-1) * regionMom, center, ptMin_, rVertex_, zVertex_, deltaEta_, deltaPhi_, magfield, measurementTracker));
0435 
0436       LogDebug("CosmicRegionalSeedGenerator")
0437           << "Final L2TrackingRegion \n "
0438           << "Position = " << center << "\n "
0439           << "Direction = " << result.back()->direction() << "\n "
0440           << "Distance from the region on the layer = " << (regionPosition - center).mag() << "\n "
0441           << "Eta = " << center.eta() << "\n "
0442           << "Phi = " << center.phi();
0443 
0444     }  //end loop on muons
0445 
0446   }  //end if SeedOnL2Muon
0447 
0448   return result;
0449 }