00001 #include <iostream>
00002 #include <memory>
00003 #include <string>
00004
00005
00006 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00007
00008
00009 #include "DataFormats/MuonReco/interface/Muon.h"
00010 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00011
00012 #include "DataFormats/JetReco/interface/CaloJet.h"
00013 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00014
00015
00016 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
00017 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
00018
00019 #include "RecoTracker/SpecialSeedGenerators/interface/CosmicRegionalSeedGenerator.h"
00020
00021 #include <TMath.h>
00022
00023 using namespace std;
00024 using namespace trigger;
00025 using namespace reco;
00026 using namespace edm;
00027
00028 CosmicRegionalSeedGenerator::CosmicRegionalSeedGenerator(edm::ParameterSet const& conf) :
00029 conf_(conf)
00030 {
00031 edm::LogInfo ("CosmicRegionalSeedGenerator") << "Begin Run:: Constructing CosmicRegionalSeedGenerator";
00032
00033 regionPSet = conf_.getParameter<edm::ParameterSet>("RegionPSet");
00034 ptMin_ = regionPSet.getParameter<double>("ptMin");
00035 rVertex_ = regionPSet.getParameter<double>("rVertex");
00036 zVertex_ = regionPSet.getParameter<double>("zVertex");
00037 deltaEta_ = regionPSet.getParameter<double>("deltaEtaRegion");
00038 deltaPhi_ = regionPSet.getParameter<double>("deltaPhiRegion");
00039
00040 edm::ParameterSet toolsPSet = conf_.getParameter<edm::ParameterSet>("ToolsPSet");
00041 thePropagatorName_ = toolsPSet.getParameter<std::string>("thePropagatorName");
00042 regionBase_ = toolsPSet.getParameter<std::string>("regionBase");
00043
00044 edm::ParameterSet collectionsPSet = conf_.getParameter<edm::ParameterSet>("CollectionsPSet");
00045 recoMuonsCollection_ = collectionsPSet.getParameter<edm::InputTag>("recoMuonsCollection");
00046 recoTrackMuonsCollection_ = collectionsPSet.getParameter<edm::InputTag>("recoTrackMuonsCollection");
00047 recoL2MuonsCollection_ = collectionsPSet.getParameter<edm::InputTag>("recoL2MuonsCollection");
00048
00049 edm::ParameterSet regionInJetsCheckPSet = conf_.getParameter<edm::ParameterSet>("RegionInJetsCheckPSet");
00050 doJetsExclusionCheck_ = regionInJetsCheckPSet.getParameter<bool>("doJetsExclusionCheck");
00051 deltaRExclusionSize_ = regionInJetsCheckPSet.getParameter<double>("deltaRExclusionSize");
00052 jetsPtMin_ = regionInJetsCheckPSet.getParameter<double>("jetsPtMin");
00053 recoCaloJetsCollection_ = regionInJetsCheckPSet.getParameter<edm::InputTag>("recoCaloJetsCollection");
00054
00055
00056 edm::LogInfo ("CosmicRegionalSeedGenerator") << "Reco muons collection: " << recoMuonsCollection_ << "\n"
00057 << "Reco tracks muons collection: " << recoTrackMuonsCollection_<< "\n"
00058 << "Reco L2 muons collection: " << recoL2MuonsCollection_;
00059 }
00060
00061 std::vector<TrackingRegion*, std::allocator<TrackingRegion*> > CosmicRegionalSeedGenerator::regions(const edm::Event& event, const edm::EventSetup& es) const
00062 {
00063
00064 std::vector<TrackingRegion* > result;
00065
00066
00067
00068
00069
00070
00071
00072
00073 if(regionBase_=="seedOnStaMuon"||regionBase_=="") {
00074
00075 LogDebug("CosmicRegionalSeedGenerator") << "Seeding on stand alone muons ";
00076
00077
00078
00079
00080
00081 edm::Handle<reco::MuonCollection> muonsHandle;
00082 event.getByLabel(recoMuonsCollection_,muonsHandle);
00083 if (!muonsHandle.isValid())
00084 {
00085 edm::LogError("CollectionNotFound") << "Error::No reco muons collection (" << recoMuonsCollection_ << ") in the event - Please verify the name of the muon collection";
00086 return result;
00087 }
00088
00089 LogDebug("CosmicRegionalSeedGenerator") << "Muons collection size = " << muonsHandle->size();
00090
00091
00092 edm::Handle<CaloJetCollection> caloJetsHandle;
00093 event.getByLabel(recoCaloJetsCollection_,caloJetsHandle);
00094
00095
00096 edm::ESHandle<Propagator> thePropagator;
00097 es.get<TrackingComponentsRecord>().get(thePropagatorName_, thePropagator);
00098
00099
00100 edm::ESHandle<TrackerGeometry> theTrackerGeometry;
00101 es.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
00102
00103 DetId outerid;
00104
00105
00106
00107
00108
00109 int nmuons = 0;
00110 for (reco::MuonCollection::const_iterator staMuon = muonsHandle->begin(); staMuon != muonsHandle->end(); ++staMuon) {
00111
00112
00113 if (!staMuon->isStandAloneMuon()) {
00114 LogDebug("CosmicRegionalSeedGenerator") << "This muon is not a stand alone muon";
00115 continue;
00116 }
00117
00118
00119 if ( abs( staMuon->standAloneMuon()->eta() ) > 1.5 ) continue;
00120
00121
00122 nmuons++;
00123 LogDebug("CosmicRegionalSeedGenerator") << "Muon stand alone found in the collection - in muons chambers: \n "
00124 << "Position = " << staMuon->standAloneMuon()->outerPosition() << "\n "
00125 << "Momentum = " << staMuon->standAloneMuon()->outerMomentum() << "\n "
00126 << "Eta = " << staMuon->standAloneMuon()->eta() << "\n "
00127 << "Phi = " << staMuon->standAloneMuon()->phi();
00128
00129
00130
00131 GlobalPoint initialRegionPosition(staMuon->standAloneMuon()->referencePoint().x(), staMuon->standAloneMuon()->referencePoint().y(), staMuon->standAloneMuon()->referencePoint().z());
00132 GlobalVector initialRegionMomentum(staMuon->standAloneMuon()->momentum().x(), staMuon->standAloneMuon()->momentum().y(), staMuon->standAloneMuon()->momentum().z());
00133 int charge = (int) staMuon->standAloneMuon()->charge();
00134
00135 LogDebug("CosmicRegionalSeedGenerator") << "Initial region - Reference point of the sta muon: \n "
00136 << "Position = " << initialRegionPosition << "\n "
00137 << "Momentum = " << initialRegionMomentum << "\n "
00138 << "Eta = " << initialRegionPosition.eta() << "\n "
00139 << "Phi = " << initialRegionPosition.phi() << "\n "
00140 << "Charge = " << charge;
00141
00142
00143 if ( staMuon->standAloneMuon()->outerPosition().y()>0 ) initialRegionMomentum *=-1;
00144 GlobalTrajectoryParameters glb_parameters(initialRegionPosition,
00145 initialRegionMomentum,
00146 charge,
00147 thePropagator->magneticField());
00148 FreeTrajectoryState fts(glb_parameters);
00149 StateOnTrackerBound onBounds(thePropagator.product());
00150 TrajectoryStateOnSurface outer = onBounds(fts);
00151
00152 if (!outer.isValid())
00153 {
00154
00155 LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid" ;
00156 continue;
00157 }
00158
00159
00160
00161 GlobalPoint regionPosition = outer.globalPosition();
00162 GlobalVector regionMom = outer.globalMomentum();
00163
00164 LogDebug("CosmicRegionalSeedGenerator") << "Region after propagation: \n "
00165 << "Position = " << outer.globalPosition() << "\n "
00166 << "Momentum = " << outer.globalMomentum() << "\n "
00167 << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
00168 << "Eta = " << outer.globalPosition().eta() << "\n "
00169 << "Phi = " << outer.globalPosition().phi();
00170
00171
00172
00173 double stepBack = 1;
00174 GlobalPoint center = regionPosition + stepBack * regionMom.unit();
00175 GlobalVector v = stepBack * regionMom.unit();
00176 LogDebug("CosmicRegionalSeedGenerator") << "Step back vector = " << v << "\n";
00177
00178
00179 if ( doJetsExclusionCheck_ ) {
00180 double delta_R_min = 1000.;
00181 for ( CaloJetCollection::const_iterator jet = caloJetsHandle->begin (); jet != caloJetsHandle->end(); jet++ ) {
00182 if ( jet->pt() < jetsPtMin_ ) continue;
00183
00184 double deta = center.eta() - jet->eta();
00185 double dphi = fabs( center.phi() - jet->phi() );
00186 if ( dphi > TMath::Pi() ) dphi = 2*TMath::Pi() - dphi;
00187
00188 double delta_R = sqrt(deta*deta + dphi*dphi);
00189 if ( delta_R < delta_R_min ) delta_R_min = delta_R;
00190
00191 }
00192
00193 if ( delta_R_min < deltaRExclusionSize_ ) {
00194 LogDebug("CosmicRegionalSeedGenerator") << "Region built too close from a jet";
00195 continue;
00196 }
00197 }
00198
00199
00200
00201 CosmicTrackingRegion *etaphiRegion = new CosmicTrackingRegion((-1)*regionMom,
00202 center,
00203 ptMin_,
00204 rVertex_,
00205 zVertex_,
00206 deltaEta_,
00207 deltaPhi_,
00208 regionPSet
00209 );
00210
00211
00212
00213
00214 result.push_back(etaphiRegion);
00215
00216 LogDebug("CosmicRegionalSeedGenerator") << "Final CosmicTrackingRegion \n "
00217 << "Position = "<< center << "\n "
00218 << "Direction = "<< etaphiRegion->direction() << "\n "
00219 << "Distance from the region on the layer = " << (regionPosition -center).mag() << "\n "
00220 << "Eta = " << center.eta() << "\n "
00221 << "Phi = " << center.phi();
00222
00223
00224 }
00225
00226 }
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 if(regionBase_=="seedOnCosmicMuon") {
00239
00240 LogDebug("CosmicRegionalSeedGenerator") << "Seeding on cosmic muons tracks";
00241
00242
00243
00244
00245
00246 edm::Handle<reco::TrackCollection> cosmicMuonsHandle;
00247 event.getByLabel(recoTrackMuonsCollection_,cosmicMuonsHandle);
00248 if (!cosmicMuonsHandle.isValid())
00249 {
00250 edm::LogError("CollectionNotFound") << "Error::No cosmic muons collection (" << recoTrackMuonsCollection_ << ") in the event - Please verify the name of the muon reco track collection";
00251 return result;
00252 }
00253
00254 LogDebug("CosmicRegionalSeedGenerator") << "Cosmic muons tracks collection size = " << cosmicMuonsHandle->size();
00255
00256
00257 edm::Handle<CaloJetCollection> caloJetsHandle;
00258 event.getByLabel(recoCaloJetsCollection_,caloJetsHandle);
00259
00260
00261 edm::ESHandle<Propagator> thePropagator;
00262 es.get<TrackingComponentsRecord>().get(thePropagatorName_, thePropagator);
00263
00264
00265 edm::ESHandle<TrackerGeometry> theTrackerGeometry;
00266 es.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
00267 DetId outerid;
00268
00269
00270
00271
00272
00273 int nmuons = 0;
00274 for (reco::TrackCollection::const_iterator cosmicMuon = cosmicMuonsHandle->begin(); cosmicMuon != cosmicMuonsHandle->end(); ++cosmicMuon) {
00275
00276
00277 if ( abs( cosmicMuon->eta() ) > 1.5 ) continue;
00278
00279 nmuons++;
00280
00281
00282 GlobalPoint initialRegionPosition(cosmicMuon->referencePoint().x(), cosmicMuon->referencePoint().y(), cosmicMuon->referencePoint().z());
00283 GlobalVector initialRegionMomentum(cosmicMuon->momentum().x(), cosmicMuon->momentum().y(), cosmicMuon->momentum().z());
00284 int charge = (int) cosmicMuon->charge();
00285
00286 LogDebug("CosmicRegionalSeedGenerator") << "Position and momentum of the muon track in the muon chambers: \n "
00287 << "x = " << cosmicMuon->outerPosition().x() << "\n "
00288 << "y = " << cosmicMuon->outerPosition().y() << "\n "
00289 << "y = " << cosmicMuon->pt() << "\n "
00290 << "Initial region - Reference point of the cosmic muon track: \n "
00291 << "Position = " << initialRegionPosition << "\n "
00292 << "Momentum = " << initialRegionMomentum << "\n "
00293 << "Eta = " << initialRegionPosition.eta() << "\n "
00294 << "Phi = " << initialRegionPosition.phi() << "\n "
00295 << "Charge = " << charge;
00296
00297
00298 if ( cosmicMuon->outerPosition().y()>0 && cosmicMuon->momentum().y()<0 ) initialRegionMomentum *=-1;
00299 GlobalTrajectoryParameters glb_parameters(initialRegionPosition,
00300 initialRegionMomentum,
00301 charge,
00302 thePropagator->magneticField());
00303 FreeTrajectoryState fts(glb_parameters);
00304 StateOnTrackerBound onBounds(thePropagator.product());
00305 TrajectoryStateOnSurface outer = onBounds(fts);
00306
00307 if (!outer.isValid())
00308 {
00309
00310 LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid" ;
00311 continue;
00312 }
00313
00314
00315
00316 GlobalPoint regionPosition = outer.globalPosition();
00317 GlobalVector regionMom = outer.globalMomentum();
00318
00319 LogDebug("CosmicRegionalSeedGenerator") << "Region after propagation: \n "
00320 << "Position = " << outer.globalPosition() << "\n "
00321 << "Momentum = " << outer.globalMomentum() << "\n "
00322 << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
00323 << "Eta = " << outer.globalPosition().eta() << "\n "
00324 << "Phi = " << outer.globalPosition().phi();
00325
00326
00327
00328 double stepBack = 1;
00329 GlobalPoint center = regionPosition + stepBack * regionMom.unit();
00330 GlobalVector v = stepBack * regionMom.unit();
00331 LogDebug("CosmicRegionalSeedGenerator") << "Step back vector = " << v << "\n";
00332
00333
00334 if ( doJetsExclusionCheck_ ) {
00335 double delta_R_min = 1000.;
00336 for ( CaloJetCollection::const_iterator jet = caloJetsHandle->begin (); jet != caloJetsHandle->end(); jet++ ) {
00337 if ( jet->pt() < jetsPtMin_ ) continue;
00338
00339 double deta = center.eta() - jet->eta();
00340 double dphi = fabs( center.phi() - jet->phi() );
00341 if ( dphi > TMath::Pi() ) dphi = 2*TMath::Pi() - dphi;
00342
00343 double delta_R = sqrt(deta*deta + dphi*dphi);
00344 if ( delta_R < delta_R_min ) delta_R_min = delta_R;
00345
00346 }
00347
00348 if ( delta_R_min < deltaRExclusionSize_ ) {
00349 LogDebug("CosmicRegionalSeedGenerator") << "Region built too close from a jet";
00350 continue;
00351 }
00352 }
00353
00354
00355 CosmicTrackingRegion *etaphiRegion = new CosmicTrackingRegion((-1)*regionMom,
00356 center,
00357 ptMin_,
00358 rVertex_,
00359 zVertex_,
00360 deltaEta_,
00361 deltaPhi_,
00362 regionPSet
00363 );
00364
00365
00366
00367 result.push_back(etaphiRegion);
00368
00369 LogDebug("CosmicRegionalSeedGenerator") << "Final CosmicTrackingRegion \n "
00370 << "Position = "<< center << "\n "
00371 << "Direction = "<< etaphiRegion->direction() << "\n "
00372 << "Distance from the region on the layer = " << (regionPosition -center).mag() << "\n "
00373 << "Eta = " << center.eta() << "\n "
00374 << "Phi = " << center.phi();
00375
00376 }
00377
00378 }
00379
00380
00381
00382
00383
00384
00385
00386 if(regionBase_=="seedOnL2Muon") {
00387
00388 LogDebug("CosmicRegionalSeedGenerator") << "Seeding on L2 muons";
00389
00390
00391
00392
00393
00394 edm::Handle<reco::RecoChargedCandidateCollection> L2MuonsHandle;
00395 event.getByLabel(recoL2MuonsCollection_,L2MuonsHandle);
00396
00397 if (!L2MuonsHandle.isValid())
00398 {
00399 edm::LogError("CollectionNotFound") << "Error::No L2 muons collection (" << recoL2MuonsCollection_ <<") in the event - Please verify the name of the L2 muon collection";
00400 return result;
00401 }
00402
00403 LogDebug("CosmicRegionalSeedGenerator") << "L2 muons collection size = " << L2MuonsHandle->size();
00404
00405
00406 edm::ESHandle<Propagator> thePropagator;
00407 es.get<TrackingComponentsRecord>().get(thePropagatorName_, thePropagator);
00408
00409
00410 edm::ESHandle<TrackerGeometry> theTrackerGeometry;
00411 es.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
00412 DetId outerid;
00413
00414
00415
00416
00417
00418 int nmuons = 0;
00419 for (reco::RecoChargedCandidateCollection::const_iterator L2Muon = L2MuonsHandle->begin(); L2Muon != L2MuonsHandle->end(); ++L2Muon) {
00420 reco::TrackRef tkL2Muon = L2Muon->get<reco::TrackRef>();
00421
00422
00423 if ( abs( tkL2Muon->eta() ) > 1.5 ) continue;
00424
00425 nmuons++;
00426
00427
00428 GlobalPoint initialRegionPosition(tkL2Muon->referencePoint().x(), tkL2Muon->referencePoint().y(), tkL2Muon->referencePoint().z());
00429 GlobalVector initialRegionMomentum(tkL2Muon->momentum().x(), tkL2Muon->momentum().y(), tkL2Muon->momentum().z());
00430 int charge = (int) tkL2Muon->charge();
00431
00432 LogDebug("CosmicRegionalSeedGenerator") << "Position and momentum of the L2 muon track in the muon chambers: \n "
00433 << "x = " << tkL2Muon->outerPosition().x() << "\n "
00434 << "y = " << tkL2Muon->outerPosition().y() << "\n "
00435 << "y = " << tkL2Muon->pt() << "\n "
00436 << "Initial region - Reference point of the L2 muon track: \n "
00437 << "Position = " << initialRegionPosition << "\n "
00438 << "Momentum = " << initialRegionMomentum << "\n "
00439 << "Eta = " << initialRegionPosition.eta() << "\n "
00440 << "Phi = " << initialRegionPosition.phi() << "\n "
00441 << "Charge = " << charge;
00442
00443
00444
00445 if ( tkL2Muon->outerPosition().y() > 0 )
00446 {
00447 LogDebug("CosmicRegionalSeedGenerator") << "L2 muon in the TOP --- Region not created";
00448 return result;
00449 }
00450
00451 GlobalTrajectoryParameters glb_parameters(initialRegionPosition,
00452 initialRegionMomentum,
00453 charge,
00454 thePropagator->magneticField());
00455 FreeTrajectoryState fts(glb_parameters);
00456 StateOnTrackerBound onBounds(thePropagator.product());
00457 TrajectoryStateOnSurface outer = onBounds(fts);
00458
00459 if (!outer.isValid())
00460 {
00461
00462 LogDebug("CosmicRegionalSeedGenerator") << "Trajectory state on surface not valid" ;
00463 continue;
00464 }
00465
00466
00467
00468 GlobalPoint regionPosition = outer.globalPosition();
00469 GlobalVector regionMom = outer.globalMomentum();
00470
00471 LogDebug("CosmicRegionalSeedGenerator") << "Region after propagation: \n "
00472 << "Position = " << outer.globalPosition() << "\n "
00473 << "Momentum = " << outer.globalMomentum() << "\n "
00474 << "R = " << regionPosition.perp() << " ---- z = " << regionPosition.z() << "\n "
00475 << "Eta = " << outer.globalPosition().eta() << "\n "
00476 << "Phi = " << outer.globalPosition().phi();
00477
00478
00479
00480 double stepBack = 1;
00481 GlobalPoint center = regionPosition + stepBack * regionMom.unit();
00482 GlobalVector v = stepBack * regionMom.unit();
00483 LogDebug("CosmicRegionalSeedGenerator") << "Step back vector = " << v << "\n";
00484
00485
00486
00487 CosmicTrackingRegion *etaphiRegion = new CosmicTrackingRegion((-1)*regionMom,
00488 center,
00489 ptMin_,
00490 rVertex_,
00491 zVertex_,
00492 deltaEta_,
00493 deltaPhi_,
00494 regionPSet
00495 );
00496
00497 result.push_back(etaphiRegion);
00498
00499 LogDebug("CosmicRegionalSeedGenerator") << "Final L2TrackingRegion \n "
00500 << "Position = "<< center << "\n "
00501 << "Direction = "<< etaphiRegion->direction() << "\n "
00502 << "Distance from the region on the layer = " << (regionPosition -center).mag() << "\n "
00503 << "Eta = " << center.eta() << "\n "
00504 << "Phi = " << center.phi();
00505
00506
00507 }
00508
00509 }
00510
00511 return result;
00512
00513 }
00514