Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:02 

0001 // -*- C++ -*- // // Package: UserCode/L1Trigger // Class: L1MuonRecoTreeProducer // /**\class L1MuonRecoTreeProducer L1MuonRecoTreeProducer.cc
0002 /*
0003 UserCode/L1Trigger/src/L1MuonRecoTreeProducer.cc
0004 
0005  Description: Produce Muon Reco tree
0006 
0007  Implementation:
0008      
0009 */
0010 //
0011 // Original Author:  Luigi Guiducci
0012 //         Created:
0013 //
0014 //
0015 
0016 // system include files
0017 #include <memory>
0018 // framework
0019 #include "FWCore/Framework/interface/Frameworkfwd.h"
0020 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0026 #include "FWCore/Framework/interface/ESHandle.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include <FWCore/Framework/interface/LuminosityBlock.h>
0029 
0030 // Muons & Tracks Data Formats
0031 #include "DataFormats/MuonReco/interface/Muon.h"
0032 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0033 #include "DataFormats/TrackReco/interface/Track.h"
0034 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0035 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0036 #include "DataFormats/GeometrySurface/interface/Plane.h"
0037 #include "DataFormats/MuonReco/interface/MuonEnergy.h"
0038 #include "DataFormats/MuonReco/interface/MuonTime.h"
0039 #include "CondFormats/AlignmentRecord/interface/TrackerSurfaceDeformationRcd.h"
0040 
0041 // Transient tracks (for extrapolations)
0042 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0043 #include "TrackingTools/TransientTrack/interface/TrackTransientTrack.h"
0044 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0045 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0046 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0047 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0048 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0049 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0050 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0051 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0052 
0053 // B Field
0054 #include "MagneticField/Engine/interface/MagneticField.h"
0055 
0056 // Geometry
0057 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0058 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0059 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0060 #include <Geometry/RPCGeometry/interface/RPCGeometry.h>
0061 #include <Geometry/RPCGeometry/interface/RPCGeomServ.h>
0062 
0063 // ROOT output stuff
0064 #include "FWCore/ServiceRegistry/interface/Service.h"
0065 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0066 #include "TH1.h"
0067 #include "TTree.h"
0068 #include "TF1.h"
0069 #include "TMath.h"
0070 
0071 #include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoMuon.h"
0072 #include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoRpcHit.h"
0073 
0074 // GP
0075 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0076 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0077 #include "DataFormats/CSCRecHit/interface/CSCRecHit2D.h"
0078 
0079 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0080 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
0081 
0082 //vertex
0083 #include "DataFormats/VertexReco/interface/Vertex.h"
0084 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0085 
0086 #include <typeinfo>
0087 
0088 // RECO TRIGGER MATCHING:
0089 #include "DataFormats/Math/interface/deltaR.h"
0090 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0091 #include "DataFormats/Common/interface/TriggerResults.h"
0092 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0093 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0094 #include "TString.h"
0095 #include "TRegexp.h"
0096 #include <utility>
0097 
0098 //
0099 // class declaration
0100 //
0101 
0102 class L1MuonRecoTreeProducer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0103 public:
0104   explicit L1MuonRecoTreeProducer(const edm::ParameterSet &);
0105   ~L1MuonRecoTreeProducer() override;
0106   TrajectoryStateOnSurface cylExtrapTrkSam(reco::TrackRef track,
0107                                            double rho,
0108                                            const MagneticField *theMagneticField,
0109                                            const Propagator &propagatorAlong,
0110                                            const Propagator &propagatorOpposite);
0111   TrajectoryStateOnSurface surfExtrapTrkSam(reco::TrackRef track,
0112                                             double z,
0113                                             const MagneticField *theMagneticField,
0114                                             const Propagator &propagatorAlong,
0115                                             const Propagator &propagatorOpposite);
0116   void empty_global();
0117   void empty_tracker();
0118   void empty_standalone();
0119   void empty_hlt();
0120 
0121   double match_trigger(std::vector<int> &trigIndices,
0122                        const trigger::TriggerObjectCollection &trigObjs,
0123                        edm::Handle<trigger::TriggerEvent> &triggerEvent,
0124                        const reco::Muon &mu);
0125 
0126 private:
0127   void beginJob(void) override;
0128   void analyze(const edm::Event &, const edm::EventSetup &) override;
0129   void endJob() override;
0130   void beginRun(const edm::Run &, const edm::EventSetup &);
0131   void endRun(const edm::Run &, const edm::EventSetup &);
0132 
0133 public:
0134   L1Analysis::L1AnalysisRecoMuon *muon;
0135   L1Analysis::L1AnalysisRecoMuonDataFormat *muonData;
0136 
0137   L1Analysis::L1AnalysisRecoRpcHit *rpcHit;
0138   L1Analysis::L1AnalysisRecoRpcHitDataFormat *rpcHitData;
0139 
0140 private:
0141   unsigned maxMuon_;
0142   unsigned maxRpcHit_;
0143 
0144   //---------------------------------------------------------------------------
0145   // TRIGGER MATCHING
0146   // member variables needed for matching, using reco muons instead of PAT
0147   //---------------------------------------------------------------------------
0148   bool triggerMatching_;
0149   edm::InputTag triggerSummaryLabel_;
0150   std::string triggerProcessLabel_;
0151   std::vector<std::string> isoTriggerNames_;
0152   std::vector<std::string> triggerNames_;
0153 
0154   std::vector<int> isoTriggerIndices_;
0155   std::vector<int> triggerIndices_;
0156   double triggerMaxDeltaR_;
0157   HLTConfigProvider hltConfig_;
0158 
0159   enum { GL_MUON = 0, SA_MUON = 1, TR_MUON = 2, TRSA_MUON = 3 };
0160 
0161   edm::ESGetToken<CSCGeometry, MuonGeometryRecord> cscGeomToken_;
0162   edm::ESGetToken<RPCGeometry, MuonGeometryRecord> rpcGeomToken_;
0163 
0164   // The Magnetic field
0165   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> theBFieldToken_;
0166 
0167   // The GlobalTrackingGeometry
0168   edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> theTrackingGeometryToken_;
0169 
0170   // Extrapolator to cylinder
0171   edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorAlongToken_;
0172   edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorOppositeToken_;
0173 
0174   FreeTrajectoryState freeTrajStateMuon(reco::TrackRef track, const MagneticField *theMagneticField);
0175 
0176   // tree
0177   TTree *tree_;
0178 
0179   bool runOnPostLS1_;
0180 
0181   // EDM input tags
0182   edm::InputTag muonTag_;
0183   edm::InputTag rpcHitTag_;
0184 };
0185 
0186 L1MuonRecoTreeProducer::L1MuonRecoTreeProducer(const edm::ParameterSet &iConfig)
0187     : cscGeomToken_(esConsumes<CSCGeometry, MuonGeometryRecord>(edm::ESInputTag("", ""))),
0188       rpcGeomToken_(esConsumes<RPCGeometry, MuonGeometryRecord>(edm::ESInputTag("", ""))),
0189       theBFieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>(edm::ESInputTag("", ""))),
0190       theTrackingGeometryToken_(
0191           esConsumes<GlobalTrackingGeometry, GlobalTrackingGeometryRecord>(edm::ESInputTag("", ""))),
0192       propagatorAlongToken_(
0193           esConsumes<Propagator, TrackingComponentsRecord>(edm::ESInputTag("SmartPropagatorAny", ""))),
0194       propagatorOppositeToken_(
0195           esConsumes<Propagator, TrackingComponentsRecord>(edm::ESInputTag("SmartPropagatorAnyOpposite", ""))) {
0196   maxMuon_ = iConfig.getParameter<unsigned int>("maxMuon");
0197   maxRpcHit_ = iConfig.getParameter<unsigned int>("maxMuon");
0198 
0199   muonTag_ = iConfig.getParameter<edm::InputTag>("muonTag");
0200   rpcHitTag_ = iConfig.getParameter<edm::InputTag>("rpcHitTag");
0201 
0202   runOnPostLS1_ = iConfig.getParameter<bool>("runOnPostLS1");
0203 
0204   muon = new L1Analysis::L1AnalysisRecoMuon();
0205   muonData = muon->getData();
0206 
0207   rpcHit = new L1Analysis::L1AnalysisRecoRpcHit();
0208   rpcHitData = rpcHit->getData();
0209 
0210   usesResource(TFileService::kSharedResource);
0211 
0212   // set up output
0213   edm::Service<TFileService> fs;
0214   tree_ = fs->make<TTree>("MuonRecoTree", "MuonRecoTree");
0215   tree_->Branch("Muon", "L1Analysis::L1AnalysisRecoMuonDataFormat", &muonData, 32000, 3);
0216   tree_->Branch("RpcHit", "L1Analysis::L1AnalysisRecoRpcHitDataFormat", &rpcHitData, 32000, 3);
0217 
0218   //---------------------------------------------------------------------------
0219   // TRIGGER MATCHING
0220   // member variables needed for matching, if using reco muons instead of PAT
0221   //---------------------------------------------------------------------------
0222   triggerMatching_ = iConfig.getUntrackedParameter<bool>("triggerMatching");
0223   triggerSummaryLabel_ = iConfig.getParameter<edm::InputTag>("triggerSummaryLabel");
0224   triggerProcessLabel_ = iConfig.getUntrackedParameter<std::string>("triggerProcessLabel");
0225   triggerNames_ = iConfig.getParameter<std::vector<std::string> >("triggerNames");
0226   isoTriggerNames_ = iConfig.getParameter<std::vector<std::string> >("isoTriggerNames");
0227   triggerMaxDeltaR_ = iConfig.getParameter<double>("triggerMaxDeltaR");
0228 }
0229 
0230 L1MuonRecoTreeProducer::~L1MuonRecoTreeProducer() {}
0231 
0232 //
0233 // member functions
0234 //
0235 
0236 double L1MuonRecoTreeProducer::match_trigger(std::vector<int> &trigIndices,
0237                                              const trigger::TriggerObjectCollection &trigObjs,
0238                                              edm::Handle<trigger::TriggerEvent> &triggerEvent,
0239                                              const reco::Muon &mu) {
0240   double matchDeltaR = 9999;
0241 
0242   for (size_t iTrigIndex = 0; iTrigIndex < trigIndices.size(); ++iTrigIndex) {
0243     int triggerIndex = trigIndices[iTrigIndex];
0244     const std::vector<std::string> moduleLabels(hltConfig_.moduleLabels(triggerIndex));
0245     // find index of the last module:
0246     const unsigned moduleIndex = hltConfig_.size(triggerIndex) - 2;
0247     // find index of HLT trigger name:
0248     const unsigned hltFilterIndex =
0249         triggerEvent->filterIndex(edm::InputTag(moduleLabels[moduleIndex], "", triggerProcessLabel_));
0250 
0251     if (hltFilterIndex < triggerEvent->sizeFilters()) {
0252       const trigger::Keys triggerKeys(triggerEvent->filterKeys(hltFilterIndex));
0253       const trigger::Vids triggerVids(triggerEvent->filterIds(hltFilterIndex));
0254 
0255       const unsigned nTriggers = triggerVids.size();
0256       for (size_t iTrig = 0; iTrig < nTriggers; ++iTrig) {
0257         // loop over all trigger objects:
0258         const trigger::TriggerObject trigObject = trigObjs[triggerKeys[iTrig]];
0259 
0260         double dRtmp = deltaR(mu, trigObject);
0261 
0262         if (dRtmp < matchDeltaR) {
0263           matchDeltaR = dRtmp;
0264         }
0265 
0266       }  // loop over different trigger objects
0267     }  // if trigger is in event (should apply hltFilter with used trigger...)
0268   }  // loop over muon candidates
0269 
0270   return matchDeltaR;
0271 }
0272 
0273 void L1MuonRecoTreeProducer::empty_global() {
0274   muonData->ch.push_back(-999999);
0275   muonData->pt.push_back(-999999);
0276   muonData->p.push_back(-999999);
0277   muonData->eta.push_back(-999999);
0278   muonData->phi.push_back(-999999);
0279   muonData->normchi2.push_back(-999999);
0280   muonData->validhits.push_back(-999999);
0281   muonData->numberOfMatchedStations.push_back(-999999);
0282   muonData->numberOfValidMuonHits.push_back(-999999);
0283   muonData->imp_point_x.push_back(-999999);
0284   muonData->imp_point_y.push_back(-999999);
0285   muonData->imp_point_z.push_back(-999999);
0286   muonData->imp_point_p.push_back(-999999);
0287   muonData->imp_point_pt.push_back(-999999);
0288   muonData->phi_hb.push_back(-999999);
0289   muonData->z_hb.push_back(-999999);
0290   muonData->r_he_p.push_back(-999999);
0291   muonData->phi_he_p.push_back(-999999);
0292   muonData->r_he_n.push_back(-999999);
0293   muonData->phi_he_n.push_back(-999999);
0294   muonData->calo_energy.push_back(-999999);
0295   muonData->calo_energy3x3.push_back(-999999);
0296   muonData->ecal_time.push_back(-999999);
0297   muonData->ecal_terr.push_back(-999999);
0298   muonData->hcal_time.push_back(-999999);
0299   muonData->hcal_terr.push_back(-999999);
0300   muonData->time_dir.push_back(-999999);
0301   muonData->time_inout.push_back(-999999);
0302   muonData->time_inout_err.push_back(-999999);
0303   muonData->time_outin.push_back(-999999);
0304   muonData->time_outin_err.push_back(-999999);
0305 
0306   muonData->hlt_isomu.push_back(-999999);
0307   muonData->hlt_mu.push_back(-999999);
0308   muonData->hlt_isoDeltaR.push_back(-999999);
0309   muonData->hlt_deltaR.push_back(-999999);
0310 }
0311 
0312 void L1MuonRecoTreeProducer::empty_tracker() {
0313   muonData->tr_ch.push_back(-999999);
0314   muonData->tr_pt.push_back(-999999);
0315   muonData->tr_p.push_back(-999999);
0316   muonData->tr_eta.push_back(-999999);
0317   muonData->tr_phi.push_back(-999999);
0318   muonData->tr_normchi2.push_back(-999999);
0319   muonData->tr_validhits.push_back(-999999);
0320   muonData->tr_validpixhits.push_back(-999999);
0321   muonData->tr_d0.push_back(-999999);
0322   muonData->tr_imp_point_x.push_back(-999999);
0323   muonData->tr_imp_point_y.push_back(-999999);
0324   muonData->tr_imp_point_z.push_back(-999999);
0325   muonData->tr_imp_point_p.push_back(-999999);
0326   muonData->tr_imp_point_pt.push_back(-999999);
0327 
0328   muonData->tr_z_mb2.push_back(-999999);
0329   muonData->tr_phi_mb2.push_back(-999999);
0330   muonData->tr_r_me2_p.push_back(-999999);
0331   muonData->tr_phi_me2_p.push_back(-999999);
0332   muonData->tr_r_me2_n.push_back(-999999);
0333   muonData->tr_phi_me2_n.push_back(-999999);
0334 
0335   muonData->tr_z_mb1.push_back(-999999);
0336   muonData->tr_phi_mb1.push_back(-999999);
0337   muonData->tr_r_me1_p.push_back(-999999);
0338   muonData->tr_phi_me1_p.push_back(-999999);
0339   muonData->tr_r_me1_n.push_back(-999999);
0340   muonData->tr_phi_me1_n.push_back(-999999);
0341 }
0342 
0343 void L1MuonRecoTreeProducer::empty_standalone() {
0344   muonData->sa_imp_point_x.push_back(-999999);
0345   muonData->sa_imp_point_y.push_back(-999999);
0346   muonData->sa_imp_point_z.push_back(-999999);
0347   muonData->sa_imp_point_p.push_back(-999999);
0348   muonData->sa_imp_point_pt.push_back(-999999);
0349   muonData->sa_z_mb2.push_back(-999999);
0350   muonData->sa_phi_mb2.push_back(-999999);
0351   muonData->sa_pseta.push_back(-999999);
0352 
0353   muonData->sa_z_hb.push_back(-999999);
0354   muonData->sa_phi_hb.push_back(-999999);
0355 
0356   muonData->sa_r_he_p.push_back(-999999);
0357   muonData->sa_phi_he_p.push_back(-999999);
0358   muonData->sa_r_he_n.push_back(-999999);
0359   muonData->sa_phi_he_n.push_back(-999999);
0360 
0361   muonData->sa_normchi2.push_back(-999999);
0362   muonData->sa_validhits.push_back(-999999);
0363   muonData->sa_ch.push_back(-999999);
0364   muonData->sa_pt.push_back(-999999);
0365   muonData->sa_p.push_back(-999999);
0366   muonData->sa_eta.push_back(-999999);
0367   muonData->sa_phi.push_back(-999999);
0368   muonData->sa_outer_pt.push_back(-999999);
0369   muonData->sa_inner_pt.push_back(-999999);
0370   muonData->sa_outer_eta.push_back(-999999);
0371   muonData->sa_inner_eta.push_back(-999999);
0372   muonData->sa_outer_phi.push_back(-999999);
0373   muonData->sa_inner_phi.push_back(-999999);
0374   muonData->sa_outer_x.push_back(-999999);
0375   muonData->sa_outer_y.push_back(-999999);
0376   muonData->sa_outer_z.push_back(-999999);
0377   muonData->sa_inner_x.push_back(-999999);
0378   muonData->sa_inner_y.push_back(-999999);
0379   muonData->sa_inner_z.push_back(-999999);
0380 
0381   muonData->sa_r_me2_p.push_back(-999999);
0382   muonData->sa_phi_me2_p.push_back(-999999);
0383 
0384   muonData->sa_r_me2_n.push_back(-999999);
0385   muonData->sa_phi_me2_n.push_back(-999999);
0386 
0387   muonData->sa_z_mb1.push_back(-999999);
0388   muonData->sa_phi_mb1.push_back(-999999);
0389   muonData->sa_r_me1_p.push_back(-999999);
0390   muonData->sa_phi_me1_p.push_back(-999999);
0391   muonData->sa_r_me1_n.push_back(-999999);
0392   muonData->sa_phi_me1_n.push_back(-999999);
0393 
0394   muonData->sa_nChambers.push_back(-999);
0395   muonData->sa_nMatches.push_back(-999);
0396 }
0397 
0398 void L1MuonRecoTreeProducer::empty_hlt() {
0399   muonData->hlt_isomu.push_back(-999999);
0400   muonData->hlt_mu.push_back(-999999);
0401   muonData->hlt_isoDeltaR.push_back(-9999999);
0402   muonData->hlt_deltaR.push_back(-999999);
0403 }
0404 
0405 //
0406 // ------------ method called to for each event  ------------
0407 //
0408 void L1MuonRecoTreeProducer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0409   float pig = TMath::Pi();
0410 
0411   muon->Reset();
0412   rpcHit->Reset();
0413 
0414   const CSCGeometry &cscGeom = iSetup.getData(cscGeomToken_);
0415   const RPCGeometry &rpcGeom = iSetup.getData(rpcGeomToken_);
0416 
0417   //Get the Magnetic field from the setup
0418   const MagneticField &theBField = iSetup.getData(theBFieldToken_);
0419   const MagneticField *theMagneticField = &theBField;
0420 
0421   // Get the GlobalTrackingGeometry from the setup
0422   const GlobalTrackingGeometry &TrackingGeometry = iSetup.getData(theTrackingGeometryToken_);
0423   const GlobalTrackingGeometry *theTrackingGeometry = &TrackingGeometry;
0424 
0425   edm::Handle<RPCRecHitCollection> rpcRecHits;
0426 
0427   iEvent.getByLabel(rpcHitTag_, rpcRecHits);
0428   if (!rpcRecHits.isValid()) {
0429     edm::LogInfo("L1Prompt") << "can't find RPCRecHitCollection with label " << rpcHitTag_.label();
0430   } else {
0431     RPCRecHitCollection::const_iterator recHitIt = rpcRecHits->begin();
0432     RPCRecHitCollection::const_iterator recHitEnd = rpcRecHits->end();
0433 
0434     int iRpcRecHits = 0;
0435 
0436     for (; recHitIt != recHitEnd; ++recHitIt) {
0437       if ((unsigned int)iRpcRecHits > maxRpcHit_ - 1)
0438         continue;
0439 
0440       int cls = recHitIt->clusterSize();
0441       int firststrip = recHitIt->firstClusterStrip();
0442       int bx = recHitIt->BunchX();
0443 
0444       RPCDetId rpcId = recHitIt->rpcId();
0445       int region = rpcId.region();
0446       int stat = rpcId.station();
0447       int sect = rpcId.sector();
0448       int layer = rpcId.layer();
0449       int subsector = rpcId.subsector();
0450       int roll = rpcId.roll();
0451       int ring = rpcId.ring();
0452 
0453       LocalPoint recHitPosLoc = recHitIt->localPosition();
0454       const BoundPlane &RPCSurface = rpcGeom.roll(rpcId)->surface();
0455       GlobalPoint recHitPosGlob = RPCSurface.toGlobal(recHitPosLoc);
0456 
0457       float xLoc = recHitPosLoc.x();
0458       float phiGlob = recHitPosGlob.phi();
0459 
0460       rpcHitData->region.push_back(region);
0461       rpcHitData->clusterSize.push_back(cls);
0462       rpcHitData->strip.push_back(firststrip);
0463       rpcHitData->bx.push_back(bx);
0464       rpcHitData->xLoc.push_back(xLoc);
0465       rpcHitData->phiGlob.push_back(phiGlob);
0466       rpcHitData->station.push_back(stat);
0467       rpcHitData->sector.push_back(sect);
0468       rpcHitData->layer.push_back(layer);
0469       rpcHitData->subsector.push_back(subsector);
0470       rpcHitData->roll.push_back(roll);
0471       rpcHitData->ring.push_back(ring);
0472       rpcHitData->muonId.push_back(-999);  // CB set to invalid now, updated when looking at muon info
0473 
0474       iRpcRecHits++;
0475     }
0476 
0477     rpcHitData->nRpcHits = iRpcRecHits;
0478   }
0479 
0480   // Get the muon candidates
0481   edm::Handle<reco::MuonCollection> mucand;
0482   iEvent.getByLabel(muonTag_, mucand);
0483   if (!mucand.isValid()) {
0484     edm::LogInfo("L1Prompt") << "can't find Muon Collection with label " << muonTag_.label();
0485     return;
0486   }
0487 
0488   // Get the beamspot
0489   edm::Handle<reco::BeamSpot> beamSpot;
0490   iEvent.getByLabel("offlineBeamSpot", beamSpot);
0491 
0492   // Get the primary vertices
0493   edm::Handle<std::vector<reco::Vertex> > vertex;
0494   iEvent.getByLabel(edm::InputTag("offlinePrimaryVertices"), vertex);
0495 
0496   // Get the propagators
0497   const Propagator &propagatorAlong = iSetup.getData(propagatorAlongToken_);
0498   const Propagator &propagatorOpposite = iSetup.getData(propagatorOppositeToken_);
0499 
0500   for (reco::MuonCollection::const_iterator imu = mucand->begin();
0501        // for(pat::MuonCollection::const_iterator imu = mucand->begin();
0502        imu != mucand->end() && (unsigned)muonData->nMuons < maxMuon_;
0503        imu++) {
0504     int type = 0;
0505     if (imu->isGlobalMuon())
0506       type = type + 1;
0507     if (imu->isStandAloneMuon())
0508       type = type + 2;
0509     if (imu->isTrackerMuon())
0510       type = type + 4;
0511     if (imu->isCaloMuon())
0512       type = type + 8;
0513 
0514     bool isTIGHT =
0515         (!vertex->empty() && imu->isGlobalMuon() && imu->globalTrack()->normalizedChi2() < 10. &&
0516          imu->globalTrack()->hitPattern().numberOfValidMuonHits() > 0 && imu->numberOfMatchedStations() > 1 &&
0517          fabs(imu->innerTrack()->dxy(vertex->at(0).position())) < 0.2 &&
0518          fabs(imu->innerTrack()->dz(vertex->at(0).position())) < 0.5 &&
0519          imu->innerTrack()->hitPattern().numberOfValidPixelHits() > 0 &&
0520          imu->innerTrack()->hitPattern().trackerLayersWithMeasurement() > 5);
0521 
0522     if (isTIGHT)
0523       type = type + 16;
0524     if (imu->isPFMuon())
0525       type = type + 32;
0526 
0527     muonData->howmanytypes.push_back(type);  // muon type is counting calo muons and multiple assignments
0528 
0529     // bool type identifiers are exclusive. is this correct? CB to check
0530     bool isSA = (!imu->isGlobalMuon() && imu->isStandAloneMuon() && !imu->isTrackerMuon());
0531     bool isTR = (!imu->isGlobalMuon() && imu->isTrackerMuon() && !imu->isStandAloneMuon());
0532     bool isGL = (imu->isGlobalMuon());  //&&!(imu->isStandAloneMuon())&&!(imu->isTrackerMuon()));
0533     bool isTRSA = (!imu->isGlobalMuon() && imu->isStandAloneMuon() && imu->isTrackerMuon());
0534 
0535     // How we fill this. We have 3 blocks of variables: muons_; muons_sa_; muons_tr_
0536     //   GL   SA   TR      Description               muon_type     Bool id   Filled Vars               Variables to -99999
0537     //   0    0    0       Muon does not exist         /             /         /                             /
0538     //   0    0    1       It is a Tracker muon      TR_MUON       isTR      muons_tr_                  muons_sa_,muons_
0539     //   0    1    0       It is a SA muon           SA_MUON       isSA      muons_sa                   muons_tr,muons_
0540     //   0    1    1       It is a SA+Tracker muon   TRSA_MUON     isTRSA    muons_sa,muons_tr          muons_
0541     //   1    0    0       It is a Global only muon  GL_MUON       isGL      muons_,muons_sa,muons_tr        /
0542     //   1    0    1       Gl w/out SA cannot exist    /            /           /                            /
0543     //   1    1    0       Gl+SA (no trk-mu match)   GL_MUON       isGL      muons_,muons_sa,muons_tr     none
0544     //   1    1    1       GL+SA+Tr (all matched)    GL_MUON       isGL      muons_,muons_sa,muons_tr     none
0545 
0546     //---------------------------------------------------------------------
0547     // TRIGGER MATCHING:
0548     // if specified the reconstructed muons are matched to a trigger
0549     //---------------------------------------------------------------------
0550     if (triggerMatching_) {
0551       double isoMatchDeltaR = 9999.;
0552       double matchDeltaR = 9999.;
0553       int hasIsoTriggered = 0;
0554       int hasTriggered = 0;
0555 
0556       // first check if the trigger results are valid:
0557       edm::Handle<edm::TriggerResults> triggerResults;
0558       iEvent.getByLabel(edm::InputTag("TriggerResults", "", triggerProcessLabel_), triggerResults);
0559 
0560       if (triggerResults.isValid()) {
0561         edm::Handle<trigger::TriggerEvent> triggerEvent;
0562         iEvent.getByLabel(triggerSummaryLabel_, triggerEvent);
0563         if (triggerEvent.isValid()) {
0564           // get trigger objects:
0565           const trigger::TriggerObjectCollection triggerObjects = triggerEvent->getObjects();
0566 
0567           matchDeltaR = match_trigger(triggerIndices_, triggerObjects, triggerEvent, (*imu));
0568           if (matchDeltaR < triggerMaxDeltaR_)
0569             hasTriggered = 1;
0570 
0571           isoMatchDeltaR = match_trigger(isoTriggerIndices_, triggerObjects, triggerEvent, (*imu));
0572           if (isoMatchDeltaR < triggerMaxDeltaR_)
0573             hasIsoTriggered = 1;
0574         }  // end if (triggerEvent.isValid())
0575       }  // end if (triggerResults.isValid())
0576 
0577       // fill trigger matching variables:
0578       muonData->hlt_isomu.push_back(hasIsoTriggered);
0579       muonData->hlt_mu.push_back(hasTriggered);
0580       muonData->hlt_isoDeltaR.push_back(isoMatchDeltaR);
0581       muonData->hlt_deltaR.push_back(matchDeltaR);
0582     } else {
0583       empty_hlt();
0584     }  // end if (triggerMatching_)
0585 
0586     if (isGL || isTR || isSA || isTRSA) {
0587       muonData->nMuons = muonData->nMuons + 1;
0588       if (isTR)
0589         muonData->type.push_back(TR_MUON);
0590       if (isGL)
0591         muonData->type.push_back(GL_MUON);
0592       if (isSA)
0593         muonData->type.push_back(SA_MUON);
0594       if (isTRSA)
0595         muonData->type.push_back(TRSA_MUON);
0596       if (!isGL)
0597         empty_global();
0598       if (!isTR && !isGL && !isTRSA)
0599         empty_tracker();
0600       if (!isSA && !isGL && !isTRSA)
0601         empty_standalone();
0602 
0603       // begin GP
0604       //---------------------------------------------------------------------
0605       // RECHIT information in CSC: only for standalone/global muons!
0606       //---------------------------------------------------------------------
0607       // An artificial rank to sort RecHits
0608       // the closer RecHit to key station/layer -> the smaller the rank
0609       // KK's original idea
0610       const int lutCSC[4][6] = {
0611           {26, 24, 22, 21, 23, 25}, {6, 4, 2, 1, 3, 5}, {16, 14, 12, 11, 13, 15}, {36, 34, 32, 31, 33, 35}};
0612 
0613       int globalTypeRCH = -999999;
0614       // float  localEtaRCH = -999999;
0615       // float  localPhiRCH = -999999;
0616       float globalEtaRCH = -999999;
0617       float globalPhiRCH = -999999;
0618 
0619       if (isSA || isGL) {
0620         trackingRecHit_iterator hit = imu->outerTrack()->recHitsBegin();
0621         trackingRecHit_iterator hitEnd = imu->outerTrack()->recHitsEnd();
0622 
0623         for (; hit != hitEnd; ++hit) {
0624           if (!((*hit)->isValid()))
0625             continue;
0626 
0627           // Hardware ID of the RecHit (in terms of wire/strips/chambers)
0628           DetId detid = (*hit)->geographicalId();
0629 
0630           // Interested in muon systems only
0631           // Look only at CSC Hits (CSC id is 2)
0632           if (detid.det() != DetId::Muon)
0633             continue;
0634           if (detid.subdetId() != MuonSubdetId::CSC)
0635             continue;
0636 
0637           CSCDetId id(detid.rawId());
0638           //std::cout << "before Lut id.station = " <<id.station()
0639           //          << " id.layer() = " << id.layer() << " globalTypeRCH = "
0640           //          << globalTypeRCH  << std::endl;
0641           // another sanity check
0642           if (id.station() < 1)
0643             continue;
0644 
0645           // Look up some stuff specific to CSCRecHit2D
0646           // std::cout << " typeid().name() " << typeid(**hit).name() << std::endl;
0647           // const CSCRecHit2D* CSChit =dynamic_cast<const CSCRecHit2D*>(&**hit);
0648 
0649           const CSCSegment *cscSegment = dynamic_cast<const CSCSegment *>(&**hit);
0650           //std::cout << "cscSegment = " << cscSegment << std::endl;
0651           if (cscSegment == nullptr)
0652             continue;
0653           // const CSCRecHit2D* CSChit =(CSCRecHit2D*)(&**hit);
0654 
0655           // std::cout << " after CSCRecHit2D, CSChit = "  << CSChit << std::endl;
0656           // LocalPoint rhitlocal = CSChit->localPosition();
0657           LocalPoint rhitlocal = cscSegment->localPosition();
0658 
0659           // for debugging purpouses
0660           //if (printLevel > 0) {
0661           //std::cout << "!!!!rhitlocal.phi = "<<rhitlocal.phi() << std::endl;
0662           //std::cout << "rhitlocal.y="<<rhitlocal.y();
0663           //std::cout << "rhitlocal.z="<<rhitlocal.z();
0664           //}
0665 
0666           GlobalPoint gp = GlobalPoint(0.0, 0.0, 0.0);
0667 
0668           const CSCChamber *cscchamber = cscGeom.chamber(id);
0669 
0670           if (!cscchamber)
0671             continue;
0672 
0673           gp = cscchamber->toGlobal(rhitlocal);
0674 
0675           // identify the rechit position
0676           //int pos = ( ((id.station()-1)*6) + (id.layer()-1) ) + (MAX_CSC_RECHIT*whichMuon);
0677 
0678           // --------------------------------------------------
0679           // this part has to be deprecated once we are sure of
0680           // the TMatrixF usage ;)
0681           // fill the rechits array
0682           //rchEtaList[pos]      = gp.eta();
0683           //rchPhiList[pos]      = gp.phi();
0684           //float phi02PI = gp.phi();
0685           //if (gp.phi() < 0) phi02PI += (2*PI);
0686 
0687           //rchPhiList_02PI[pos] = phi02PI;
0688           // --------------------------------------------------
0689 
0690           // --------------------------------------------------
0691           // See if this hit is closer to the "key" position
0692           //if( lutCSC[id.station()-1][id.layer()-1]<globalTypeRCH || globalTypeRCH<0 ){
0693           if (lutCSC[id.station() - 1][3 - 1] < globalTypeRCH || globalTypeRCH < 0) {
0694             //globalTypeRCH    = lutCSC[id.station()-1][id.layer()-1];
0695             globalTypeRCH = lutCSC[id.station() - 1][3 - 1];
0696             // localEtaRCH      = rhitlocal.eta();
0697             // localPhiRCH      = rhitlocal.phi();
0698             globalEtaRCH = gp.eta();
0699             globalPhiRCH = gp.phi();  // phi from -pi to pi
0700             //std::cout << "globalEtaRCH =  "  <<globalEtaRCH
0701             //          << " globalPhiRCH = " << globalPhiRCH << std::endl;
0702             if (globalPhiRCH < 0)
0703               globalPhiRCH = globalPhiRCH + 2 * pig;  // convert to [0; 2pi]
0704           }
0705           // --------------------------------------------------
0706         }
0707 
0708         hit = imu->outerTrack()->recHitsBegin();
0709 
0710         for (; hit != hitEnd; hit++) {
0711           if (!((*hit)->isValid()))
0712             continue;
0713 
0714           DetId detId = (*hit)->geographicalId();
0715 
0716           if (detId.det() != DetId::Muon)
0717             continue;
0718           if (detId.subdetId() != MuonSubdetId::RPC)
0719             continue;
0720 
0721           RPCDetId rpcId = (RPCDetId)(*hit)->geographicalId();
0722 
0723           int region = rpcId.region();
0724           int stat = rpcId.station();
0725           int sect = rpcId.sector();
0726           int layer = rpcId.layer();
0727           int subsector = rpcId.subsector();
0728           int roll = rpcId.roll();
0729           int ring = rpcId.ring();
0730 
0731           float xLoc = (*hit)->localPosition().x();
0732 
0733           for (int iRpcHit = 0; iRpcHit < rpcHitData->nRpcHits; ++iRpcHit) {
0734             if (region == rpcHitData->region.at(iRpcHit) && stat == rpcHitData->station.at(iRpcHit) &&
0735                 sect == rpcHitData->sector.at(iRpcHit) && layer == rpcHitData->layer.at(iRpcHit) &&
0736                 subsector == rpcHitData->subsector.at(iRpcHit) && roll == rpcHitData->roll.at(iRpcHit) &&
0737                 ring == rpcHitData->ring.at(iRpcHit) && fabs(xLoc - rpcHitData->xLoc.at(iRpcHit)) < 0.01)
0738 
0739               rpcHitData->muonId.at(iRpcHit) = muonData->nMuons;  // CB rpc hit belongs to mu imu
0740                   // due to cleaning as in 2012 1 rpc hit belogns to 1 mu
0741           }
0742         }
0743       }
0744       // at the end of the loop, write only the best rechit
0745       // if(globalPhiRCH > -100.) std::cout << "globalPhiRCH = " << globalPhiRCH
0746       // << "globalEtaRCH = " << globalEtaRCH << std::endl;
0747       muonData->rchCSCtype.push_back(globalTypeRCH);
0748       muonData->rchPhi.push_back(globalPhiRCH);
0749       muonData->rchEta.push_back(globalEtaRCH);
0750       // end GP
0751 
0752       if (isGL) {
0753         // Filling calo energies and calo times
0754         if (imu->isEnergyValid()) {
0755           reco::MuonEnergy muon_energy;
0756           muon_energy = imu->calEnergy();
0757           muonData->calo_energy.push_back(muon_energy.tower);
0758           muonData->calo_energy3x3.push_back(muon_energy.towerS9);
0759           muonData->ecal_time.push_back(muon_energy.ecal_time);
0760           muonData->ecal_terr.push_back(muon_energy.ecal_timeError);
0761           muonData->hcal_time.push_back(muon_energy.hcal_time);
0762           muonData->hcal_terr.push_back(muon_energy.hcal_timeError);
0763         } else {
0764           muonData->calo_energy.push_back(-999999);
0765           muonData->calo_energy3x3.push_back(-999999);
0766           muonData->ecal_time.push_back(-999999);
0767           muonData->ecal_terr.push_back(-999999);
0768           muonData->hcal_time.push_back(-999999);
0769           muonData->hcal_terr.push_back(-999999);
0770         }
0771 
0772         // Filling muon time data
0773         if (imu->isTimeValid()) {
0774           reco::MuonTime muon_time;
0775           muon_time = imu->time();
0776           muonData->time_dir.push_back(muon_time.direction());
0777           muonData->time_inout.push_back(muon_time.timeAtIpInOut);
0778           muonData->time_inout_err.push_back(muon_time.timeAtIpInOutErr);
0779           muonData->time_outin.push_back(muon_time.timeAtIpOutIn);
0780           muonData->time_outin_err.push_back(muon_time.timeAtIpOutInErr);
0781         } else {
0782           muonData->time_dir.push_back(-999999);
0783           muonData->time_inout.push_back(-999999);
0784           muonData->time_inout_err.push_back(-999999);
0785           muonData->time_outin.push_back(-999999);
0786           muonData->time_outin_err.push_back(-999999);
0787         }
0788 
0789         // Use the track now, and make a transient track out of it (for extrapolations)
0790         reco::TrackRef glb_mu = imu->globalTrack();
0791         reco::TransientTrack ttrack(*glb_mu, theMagneticField, theTrackingGeometry);
0792 
0793         // Track quantities
0794         muonData->ch.push_back(glb_mu->charge());
0795         muonData->pt.push_back(glb_mu->pt());
0796         muonData->p.push_back(glb_mu->p());
0797         muonData->eta.push_back(glb_mu->eta());
0798         muonData->phi.push_back(glb_mu->phi());
0799         muonData->normchi2.push_back(glb_mu->normalizedChi2());
0800         muonData->validhits.push_back(glb_mu->numberOfValidHits());
0801         muonData->numberOfMatchedStations.push_back(imu->numberOfMatchedStations());
0802         muonData->numberOfValidMuonHits.push_back(glb_mu->hitPattern().numberOfValidMuonHits());
0803 
0804         // Extrapolation to IP
0805         if (ttrack.impactPointTSCP().isValid()) {
0806           muonData->imp_point_x.push_back(ttrack.impactPointTSCP().position().x());
0807           muonData->imp_point_y.push_back(ttrack.impactPointTSCP().position().y());
0808           muonData->imp_point_z.push_back(ttrack.impactPointTSCP().position().z());
0809           muonData->imp_point_p.push_back(sqrt(ttrack.impactPointTSCP().position().mag()));
0810           muonData->imp_point_pt.push_back(sqrt(ttrack.impactPointTSCP().position().perp2()));
0811         } else {
0812           muonData->imp_point_x.push_back(-999999);
0813           muonData->imp_point_y.push_back(-999999);
0814           muonData->imp_point_z.push_back(-999999);
0815           muonData->imp_point_p.push_back(-999999);
0816           muonData->imp_point_pt.push_back(-999999);
0817         }
0818 
0819         if (!runOnPostLS1_) {
0820           //Extrapolation to HB
0821           TrajectoryStateOnSurface tsos;
0822           tsos = cylExtrapTrkSam(glb_mu, 235, theMagneticField, propagatorAlong, propagatorOpposite);
0823           if (tsos.isValid()) {
0824             double xx = tsos.globalPosition().x();
0825             double yy = tsos.globalPosition().y();
0826             double zz = tsos.globalPosition().z();
0827             muonData->z_hb.push_back(zz);
0828             double rr = sqrt(xx * xx + yy * yy);
0829             double cosphi = xx / rr;
0830             if (yy >= 0)
0831               muonData->phi_hb.push_back(acos(cosphi));
0832             else
0833               muonData->phi_hb.push_back(2 * pig - acos(cosphi));
0834           } else {
0835             muonData->phi_hb.push_back(-999999);
0836             muonData->z_hb.push_back(-999999);
0837           }
0838 
0839           //Extrapolation to HE+
0840           tsos = surfExtrapTrkSam(glb_mu, 479, theMagneticField, propagatorAlong, propagatorOpposite);
0841           if (tsos.isValid()) {
0842             double xx = tsos.globalPosition().x();
0843             double yy = tsos.globalPosition().y();
0844             double rr = sqrt(xx * xx + yy * yy);
0845             muonData->r_he_p.push_back(rr);
0846             double cosphi = xx / rr;
0847             if (yy >= 0)
0848               muonData->phi_he_p.push_back(acos(cosphi));
0849             else
0850               muonData->phi_he_p.push_back(2 * pig - acos(cosphi));
0851           } else {
0852             muonData->r_he_p.push_back(-999999);
0853             muonData->phi_he_p.push_back(-999999);
0854           }
0855 
0856           //Extrapolation to HE-
0857           tsos = surfExtrapTrkSam(glb_mu, -479, theMagneticField, propagatorAlong, propagatorOpposite);
0858           if (tsos.isValid()) {
0859             double xx = tsos.globalPosition().x();
0860             double yy = tsos.globalPosition().y();
0861             double rr = sqrt(xx * xx + yy * yy);
0862             muonData->r_he_n.push_back(rr);
0863             double cosphi = xx / rr;
0864             if (yy >= 0)
0865               muonData->phi_he_n.push_back(acos(cosphi));
0866             else
0867               muonData->phi_he_n.push_back(2 * pig - acos(cosphi));
0868           } else {
0869             muonData->r_he_n.push_back(-999999);
0870             muonData->phi_he_n.push_back(-999999);
0871           }
0872         } else {
0873           muonData->phi_hb.push_back(-999999);
0874           muonData->z_hb.push_back(-999999);
0875           muonData->r_he_p.push_back(-999999);
0876           muonData->phi_he_p.push_back(-999999);
0877           muonData->r_he_n.push_back(-999999);
0878           muonData->phi_he_n.push_back(-999999);
0879         }
0880       }  // end of IF IS GLOBAL
0881 
0882       // If global muon or tracker muon, fill the tracker track quantities
0883       if (isTR || isGL || isTRSA) {
0884         // Take the tracker track and build a transient track out of it
0885         reco::TrackRef tr_mu = imu->innerTrack();
0886         reco::TransientTrack ttrack(*tr_mu, theMagneticField, theTrackingGeometry);
0887         // Fill track quantities
0888         muonData->tr_ch.push_back(tr_mu->charge());
0889         muonData->tr_pt.push_back(tr_mu->pt());
0890         muonData->tr_p.push_back(tr_mu->p());
0891         muonData->tr_eta.push_back(tr_mu->eta());
0892         muonData->tr_phi.push_back(tr_mu->phi());
0893         muonData->tr_normchi2.push_back(tr_mu->normalizedChi2());
0894         muonData->tr_validhits.push_back(tr_mu->numberOfValidHits());
0895         muonData->tr_validpixhits.push_back(tr_mu->hitPattern().numberOfValidPixelHits());
0896         // find d0 from vertex position
0897         //////////// Beam spot //////////////
0898         edm::Handle<reco::BeamSpot> beamSpot;
0899         iEvent.getByLabel("offlineBeamSpot", beamSpot);
0900         muonData->tr_d0.push_back(tr_mu->dxy(beamSpot->position()));
0901 
0902         // Extrapolation to the IP
0903         if (ttrack.impactPointTSCP().isValid()) {
0904           muonData->tr_imp_point_x.push_back(ttrack.impactPointTSCP().position().x());
0905           muonData->tr_imp_point_y.push_back(ttrack.impactPointTSCP().position().y());
0906           muonData->tr_imp_point_z.push_back(ttrack.impactPointTSCP().position().z());
0907           muonData->tr_imp_point_p.push_back(sqrt(ttrack.impactPointTSCP().position().mag()));
0908           muonData->tr_imp_point_pt.push_back(sqrt(ttrack.impactPointTSCP().position().perp2()));
0909         } else {
0910           muonData->tr_imp_point_x.push_back(-999999);
0911           muonData->tr_imp_point_y.push_back(-999999);
0912           muonData->tr_imp_point_z.push_back(-999999);
0913           muonData->tr_imp_point_p.push_back(-999999);
0914           muonData->tr_imp_point_pt.push_back(-999999);
0915         }
0916 
0917         if (!runOnPostLS1_) {
0918           TrajectoryStateOnSurface tsos;
0919           tsos = cylExtrapTrkSam(
0920               tr_mu, 410, theMagneticField, propagatorAlong, propagatorOpposite);  // track at MB1 radius - extrapolation
0921           if (tsos.isValid()) {
0922             double xx = tsos.globalPosition().x();
0923             double yy = tsos.globalPosition().y();
0924             double zz = tsos.globalPosition().z();
0925             muonData->tr_z_mb1.push_back(zz);
0926             double rr = sqrt(xx * xx + yy * yy);
0927             double cosphi = xx / rr;
0928             if (yy >= 0)
0929               muonData->tr_phi_mb1.push_back(acos(cosphi));
0930             else
0931               muonData->tr_phi_mb1.push_back(2 * pig - acos(cosphi));
0932           } else {
0933             muonData->tr_z_mb1.push_back(-999999);
0934             muonData->tr_phi_mb1.push_back(-999999);
0935           }
0936 
0937           tsos = cylExtrapTrkSam(
0938               tr_mu, 500, theMagneticField, propagatorAlong, propagatorOpposite);  // track at MB2 radius - extrapolation
0939           if (tsos.isValid()) {
0940             double xx = tsos.globalPosition().x();
0941             double yy = tsos.globalPosition().y();
0942             double zz = tsos.globalPosition().z();
0943             muonData->tr_z_mb2.push_back(zz);
0944             double rr = sqrt(xx * xx + yy * yy);
0945             double cosphi = xx / rr;
0946             if (yy >= 0)
0947               muonData->tr_phi_mb2.push_back(acos(cosphi));
0948             else
0949               muonData->tr_phi_mb2.push_back(2 * pig - acos(cosphi));
0950           } else {
0951             muonData->tr_z_mb2.push_back(-999999);
0952             muonData->tr_phi_mb2.push_back(-999999);
0953           }
0954 
0955           tsos = surfExtrapTrkSam(
0956               tr_mu, 630, theMagneticField, propagatorAlong, propagatorOpposite);  // track at ME1+ plane - extrapolation
0957           if (tsos.isValid()) {
0958             double xx = tsos.globalPosition().x();
0959             double yy = tsos.globalPosition().y();
0960             double rr = sqrt(xx * xx + yy * yy);
0961             muonData->tr_r_me1_p.push_back(rr);
0962             double cosphi = xx / rr;
0963             if (yy >= 0)
0964               muonData->tr_phi_me1_p.push_back(acos(cosphi));
0965             else
0966               muonData->tr_phi_me1_p.push_back(2 * pig - acos(cosphi));
0967           } else {
0968             muonData->tr_r_me1_p.push_back(-999999);
0969             muonData->tr_phi_me1_p.push_back(-999999);
0970           }
0971 
0972           tsos = surfExtrapTrkSam(
0973               tr_mu, 790, theMagneticField, propagatorAlong, propagatorOpposite);  // track at ME2+ plane - extrapolation
0974           if (tsos.isValid()) {
0975             double xx = tsos.globalPosition().x();
0976             double yy = tsos.globalPosition().y();
0977             double rr = sqrt(xx * xx + yy * yy);
0978             muonData->tr_r_me2_p.push_back(rr);
0979             double cosphi = xx / rr;
0980             if (yy >= 0)
0981               muonData->tr_phi_me2_p.push_back(acos(cosphi));
0982             else
0983               muonData->tr_phi_me2_p.push_back(2 * pig - acos(cosphi));
0984           } else {
0985             muonData->tr_r_me2_p.push_back(-999999);
0986             muonData->tr_phi_me2_p.push_back(-999999);
0987           }
0988 
0989           tsos = surfExtrapTrkSam(tr_mu,
0990                                   -630,
0991                                   theMagneticField,
0992                                   propagatorAlong,
0993                                   propagatorOpposite);  // track at ME1- plane - extrapolation
0994           if (tsos.isValid()) {
0995             double xx = tsos.globalPosition().x();
0996             double yy = tsos.globalPosition().y();
0997             double rr = sqrt(xx * xx + yy * yy);
0998             muonData->tr_r_me1_n.push_back(rr);
0999             double cosphi = xx / rr;
1000             if (yy >= 0)
1001               muonData->tr_phi_me1_n.push_back(acos(cosphi));
1002             else
1003               muonData->tr_phi_me1_n.push_back(2 * pig - acos(cosphi));
1004           } else {
1005             muonData->tr_r_me1_n.push_back(-999999);
1006             muonData->tr_phi_me1_n.push_back(-999999);
1007           }
1008 
1009           tsos = surfExtrapTrkSam(tr_mu,
1010                                   -790,
1011                                   theMagneticField,
1012                                   propagatorAlong,
1013                                   propagatorOpposite);  // track at ME2- plane - extrapolation
1014           if (tsos.isValid()) {
1015             double xx = tsos.globalPosition().x();
1016             double yy = tsos.globalPosition().y();
1017             double rr = sqrt(xx * xx + yy * yy);
1018             muonData->tr_r_me2_n.push_back(rr);
1019             double cosphi = xx / rr;
1020             if (yy >= 0)
1021               muonData->tr_phi_me2_n.push_back(acos(cosphi));
1022             else
1023               muonData->tr_phi_me2_n.push_back(2 * pig - acos(cosphi));
1024           } else {
1025             muonData->tr_r_me2_n.push_back(-999999);
1026             muonData->tr_phi_me2_n.push_back(-999999);
1027           }
1028         } else {
1029           muonData->tr_z_mb1.push_back(-999999);
1030           muonData->tr_phi_mb1.push_back(-999999);
1031           muonData->tr_z_mb2.push_back(-999999);
1032           muonData->tr_phi_mb2.push_back(-999999);
1033           muonData->tr_r_me1_p.push_back(-999999);
1034           muonData->tr_phi_me1_p.push_back(-999999);
1035           muonData->tr_r_me2_p.push_back(-999999);
1036           muonData->tr_phi_me2_p.push_back(-999999);
1037           muonData->tr_r_me1_n.push_back(-999999);
1038           muonData->tr_phi_me1_n.push_back(-999999);
1039           muonData->tr_r_me2_n.push_back(-999999);
1040           muonData->tr_phi_me2_n.push_back(-999999);
1041         }
1042 
1043       }  // end of IF IS TRACKER
1044 
1045       // If global muon or sa muon, fill the sa track quantities
1046       if (isGL || isSA || isTRSA) {
1047         muonData->sa_nChambers.push_back(imu->numberOfChambers());
1048         muonData->sa_nMatches.push_back(imu->numberOfMatches());
1049 
1050         // Take the SA track and build a transient track out of it
1051         reco::TrackRef sa_mu = imu->outerTrack();
1052         reco::TransientTrack ttrack(*sa_mu, theMagneticField, theTrackingGeometry);
1053 
1054         // Extrapolation to IP
1055         if (ttrack.impactPointTSCP().isValid()) {
1056           muonData->sa_imp_point_x.push_back(ttrack.impactPointTSCP().position().x());
1057           muonData->sa_imp_point_y.push_back(ttrack.impactPointTSCP().position().y());
1058           muonData->sa_imp_point_z.push_back(ttrack.impactPointTSCP().position().z());
1059           muonData->sa_imp_point_p.push_back(sqrt(ttrack.impactPointTSCP().position().mag()));
1060           muonData->sa_imp_point_pt.push_back(sqrt(ttrack.impactPointTSCP().position().perp2()));
1061         } else {
1062           muonData->sa_imp_point_x.push_back(-999999);
1063           muonData->sa_imp_point_y.push_back(-999999);
1064           muonData->sa_imp_point_z.push_back(-999999);
1065           muonData->sa_imp_point_p.push_back(-999999);
1066           muonData->sa_imp_point_pt.push_back(-999999);
1067         }
1068 
1069         // Extrapolation to MB2
1070 
1071         TrajectoryStateOnSurface tsos;
1072         tsos = cylExtrapTrkSam(
1073             sa_mu, 410, theMagneticField, propagatorAlong, propagatorOpposite);  // track at MB1 radius - extrapolation
1074         if (tsos.isValid()) {
1075           double xx = tsos.globalPosition().x();
1076           double yy = tsos.globalPosition().y();
1077           double zz = tsos.globalPosition().z();
1078           muonData->sa_z_mb1.push_back(zz);
1079           double rr = sqrt(xx * xx + yy * yy);
1080           double cosphi = xx / rr;
1081           if (yy >= 0)
1082             muonData->sa_phi_mb1.push_back(acos(cosphi));
1083           else
1084             muonData->sa_phi_mb1.push_back(2 * pig - acos(cosphi));
1085         } else {
1086           muonData->sa_z_mb1.push_back(-999999);
1087           muonData->sa_phi_mb1.push_back(-999999);
1088         }
1089 
1090         tsos = cylExtrapTrkSam(
1091             sa_mu, 500, theMagneticField, propagatorAlong, propagatorOpposite);  // track at MB2 radius - extrapolation
1092         if (tsos.isValid()) {
1093           double xx = tsos.globalPosition().x();
1094           double yy = tsos.globalPosition().y();
1095           double zz = tsos.globalPosition().z();
1096           muonData->sa_z_mb2.push_back(zz);
1097           double rr = sqrt(xx * xx + yy * yy);
1098           double cosphi = xx / rr;
1099           if (yy >= 0)
1100             muonData->sa_phi_mb2.push_back(acos(cosphi));
1101           else
1102             muonData->sa_phi_mb2.push_back(2 * pig - acos(cosphi));
1103           double abspseta = -log(tan(atan(fabs(rr / zz)) / 2.0));
1104           if (zz >= 0)
1105             muonData->sa_pseta.push_back(abspseta);
1106           else
1107             muonData->sa_pseta.push_back(-abspseta);
1108         } else {
1109           muonData->sa_z_mb2.push_back(-999999);
1110           muonData->sa_phi_mb2.push_back(-999999);
1111           muonData->sa_pseta.push_back(-999999);
1112         }
1113 
1114         tsos = surfExtrapTrkSam(
1115             sa_mu, 630, theMagneticField, propagatorAlong, propagatorOpposite);  // track at ME1+ plane - extrapolation
1116         if (tsos.isValid()) {
1117           double xx = tsos.globalPosition().x();
1118           double yy = tsos.globalPosition().y();
1119           double rr = sqrt(xx * xx + yy * yy);
1120           muonData->sa_r_me1_p.push_back(rr);
1121           double cosphi = xx / rr;
1122           if (yy >= 0)
1123             muonData->sa_phi_me1_p.push_back(acos(cosphi));
1124           else
1125             muonData->sa_phi_me1_p.push_back(2 * pig - acos(cosphi));
1126         } else {
1127           muonData->sa_r_me1_p.push_back(-999999);
1128           muonData->sa_phi_me1_p.push_back(-999999);
1129         }
1130 
1131         // Extrapolation to ME2+
1132         tsos = surfExtrapTrkSam(sa_mu, 790, theMagneticField, propagatorAlong, propagatorOpposite);
1133         if (tsos.isValid()) {
1134           double xx = tsos.globalPosition().x();
1135           double yy = tsos.globalPosition().y();
1136           double rr = sqrt(xx * xx + yy * yy);
1137           muonData->sa_r_me2_p.push_back(rr);
1138           double cosphi = xx / rr;
1139           if (yy >= 0)
1140             muonData->sa_phi_me2_p.push_back(acos(cosphi));
1141           else
1142             muonData->sa_phi_me2_p.push_back(2 * pig - acos(cosphi));
1143         } else {
1144           muonData->sa_r_me2_p.push_back(-999999);
1145           muonData->sa_phi_me2_p.push_back(-999999);
1146         }
1147 
1148         tsos = surfExtrapTrkSam(
1149             sa_mu, -630, theMagneticField, propagatorAlong, propagatorOpposite);  // track at ME1- plane - extrapolation
1150         if (tsos.isValid()) {
1151           double xx = tsos.globalPosition().x();
1152           double yy = tsos.globalPosition().y();
1153           double rr = sqrt(xx * xx + yy * yy);
1154           muonData->sa_r_me1_n.push_back(rr);
1155           double cosphi = xx / rr;
1156           if (yy >= 0)
1157             muonData->sa_phi_me1_n.push_back(acos(cosphi));
1158           else
1159             muonData->sa_phi_me1_n.push_back(2 * pig - acos(cosphi));
1160         } else {
1161           muonData->sa_r_me1_n.push_back(-999999);
1162           muonData->sa_phi_me1_n.push_back(-999999);
1163         }
1164 
1165         // Extrapolation to ME2-
1166         tsos = surfExtrapTrkSam(
1167             sa_mu, -790, theMagneticField, propagatorAlong, propagatorOpposite);  // track at ME2- disk - extrapolation
1168         if (tsos.isValid()) {
1169           double xx = tsos.globalPosition().x();
1170           double yy = tsos.globalPosition().y();
1171           double rr = sqrt(xx * xx + yy * yy);
1172           muonData->sa_r_me2_n.push_back(rr);
1173           double cosphi = xx / rr;
1174           if (yy >= 0)
1175             muonData->sa_phi_me2_n.push_back(acos(cosphi));
1176           else
1177             muonData->sa_phi_me2_n.push_back(2 * pig - acos(cosphi));
1178         } else {
1179           muonData->sa_r_me2_n.push_back(-999999);
1180           muonData->sa_phi_me2_n.push_back(-999999);
1181         }
1182 
1183         // Extrapolation to HB
1184         tsos = cylExtrapTrkSam(
1185             sa_mu, 235, theMagneticField, propagatorAlong, propagatorOpposite);  // track at HB radius - extrapolation
1186         if (tsos.isValid()) {
1187           double xx = tsos.globalPosition().x();
1188           double yy = tsos.globalPosition().y();
1189           double zz = tsos.globalPosition().z();
1190           muonData->sa_z_hb.push_back(zz);
1191           double rr = sqrt(xx * xx + yy * yy);
1192           double cosphi = xx / rr;
1193           if (yy >= 0)
1194             muonData->sa_phi_hb.push_back(acos(cosphi));
1195           else
1196             muonData->sa_phi_hb.push_back(2 * pig - acos(cosphi));
1197         } else {
1198           muonData->sa_z_hb.push_back(-999999);
1199           muonData->sa_phi_hb.push_back(-999999);
1200         }
1201 
1202         // Extrapolation to HE+
1203         tsos = surfExtrapTrkSam(sa_mu, 479, theMagneticField, propagatorAlong, propagatorOpposite);
1204         if (tsos.isValid()) {
1205           double xx = tsos.globalPosition().x();
1206           double yy = tsos.globalPosition().y();
1207           double rr = sqrt(xx * xx + yy * yy);
1208           muonData->sa_r_he_p.push_back(rr);
1209           double cosphi = xx / rr;
1210           if (yy >= 0)
1211             muonData->sa_phi_he_p.push_back(acos(cosphi));
1212           else
1213             muonData->sa_phi_he_p.push_back(2 * pig - acos(cosphi));
1214         } else {
1215           muonData->sa_r_he_p.push_back(-999999);
1216           muonData->sa_phi_he_p.push_back(-999999);
1217         }
1218 
1219         // Extrapolation to HE-
1220         tsos = surfExtrapTrkSam(sa_mu, -479, theMagneticField, propagatorAlong, propagatorOpposite);
1221         if (tsos.isValid()) {
1222           double xx = tsos.globalPosition().x();
1223           double yy = tsos.globalPosition().y();
1224           double rr = sqrt(xx * xx + yy * yy);
1225           muonData->sa_r_he_n.push_back(rr);
1226           double cosphi = xx / rr;
1227           if (yy >= 0)
1228             muonData->sa_phi_he_n.push_back(acos(cosphi));
1229           else
1230             muonData->sa_phi_he_n.push_back(2 * pig - acos(cosphi));
1231         } else {
1232           muonData->sa_r_he_n.push_back(-999999);
1233           muonData->sa_phi_he_n.push_back(-999999);
1234         }
1235 
1236         //  Several track quantities
1237         muonData->sa_normchi2.push_back(sa_mu->normalizedChi2());
1238         muonData->sa_validhits.push_back(sa_mu->numberOfValidHits());
1239         muonData->sa_ch.push_back(sa_mu->charge());
1240         muonData->sa_pt.push_back(sa_mu->pt());
1241         muonData->sa_p.push_back(sa_mu->p());
1242         muonData->sa_eta.push_back(sa_mu->eta());
1243         muonData->sa_phi.push_back(sa_mu->phi());
1244         muonData->sa_outer_pt.push_back(sqrt(sa_mu->outerMomentum().Perp2()));
1245         muonData->sa_inner_pt.push_back(sqrt(sa_mu->innerMomentum().Perp2()));
1246         muonData->sa_outer_eta.push_back(sa_mu->outerMomentum().Eta());
1247         muonData->sa_inner_eta.push_back(sa_mu->innerMomentum().Eta());
1248         muonData->sa_outer_phi.push_back(sa_mu->outerMomentum().Phi());
1249         muonData->sa_inner_phi.push_back(sa_mu->innerMomentum().Phi());
1250         muonData->sa_outer_x.push_back(sa_mu->outerPosition().x());
1251         muonData->sa_outer_y.push_back(sa_mu->outerPosition().y());
1252         muonData->sa_outer_z.push_back(sa_mu->outerPosition().z());
1253         muonData->sa_inner_x.push_back(sa_mu->innerPosition().x());
1254         muonData->sa_inner_y.push_back(sa_mu->innerPosition().y());
1255         muonData->sa_inner_z.push_back(sa_mu->innerPosition().z());
1256 
1257       }  // end of IF IS STANDALONE
1258 
1259     }  // end of if type 012 found
1260   }  // end of muon loop
1261 
1262   tree_->Fill();
1263 }
1264 
1265 // to get the track position info at a particular rho
1266 TrajectoryStateOnSurface L1MuonRecoTreeProducer::cylExtrapTrkSam(reco::TrackRef track,
1267                                                                  double rho,
1268                                                                  const MagneticField *theMagneticField,
1269                                                                  const Propagator &propagatorAlong,
1270                                                                  const Propagator &propagatorOpposite) {
1271   Cylinder::PositionType pos(0, 0, 0);
1272   Cylinder::RotationType rot;
1273   Cylinder::CylinderPointer myCylinder = Cylinder::build(pos, rot, rho);
1274 
1275   FreeTrajectoryState recoStart = freeTrajStateMuon(track, theMagneticField);
1276   TrajectoryStateOnSurface recoProp;
1277   recoProp = propagatorAlong.propagate(recoStart, *myCylinder);
1278   if (!recoProp.isValid()) {
1279     recoProp = propagatorOpposite.propagate(recoStart, *myCylinder);
1280   }
1281   return recoProp;
1282 }
1283 
1284 // to get track position at a particular (xy) plane given its z
1285 TrajectoryStateOnSurface L1MuonRecoTreeProducer::surfExtrapTrkSam(reco::TrackRef track,
1286                                                                   double z,
1287                                                                   const MagneticField *theMagneticField,
1288                                                                   const Propagator &propagatorAlong,
1289                                                                   const Propagator &propagatorOpposite) {
1290   Plane::PositionType pos(0, 0, z);
1291   Plane::RotationType rot;
1292   Plane::PlanePointer myPlane = Plane::build(pos, rot);
1293 
1294   FreeTrajectoryState recoStart = freeTrajStateMuon(track, theMagneticField);
1295   TrajectoryStateOnSurface recoProp;
1296   recoProp = propagatorAlong.propagate(recoStart, *myPlane);
1297   if (!recoProp.isValid()) {
1298     recoProp = propagatorOpposite.propagate(recoStart, *myPlane);
1299   }
1300   return recoProp;
1301 }
1302 
1303 FreeTrajectoryState L1MuonRecoTreeProducer::freeTrajStateMuon(reco::TrackRef track,
1304                                                               const MagneticField *theMagneticField) {
1305   GlobalPoint innerPoint(track->innerPosition().x(), track->innerPosition().y(), track->innerPosition().z());
1306   GlobalVector innerVec(track->innerMomentum().x(), track->innerMomentum().y(), track->innerMomentum().z());
1307 
1308   FreeTrajectoryState recoStart(innerPoint, innerVec, track->charge(), theMagneticField);
1309 
1310   return recoStart;
1311 }
1312 
1313 void L1MuonRecoTreeProducer::beginRun(const edm::Run &run, const edm::EventSetup &eventSetup) {
1314   // Prepare for trigger matching for each new run:
1315   // Look up triggetIndices in the HLT config for the different paths
1316   if (triggerMatching_) {
1317     bool changed = true;
1318     if (!hltConfig_.init(run, eventSetup, triggerProcessLabel_, changed)) {
1319       // if you can't initialize hlt configuration, crash!
1320       std::cout << "Error: didn't find process" << triggerProcessLabel_ << std::endl;
1321       assert(false);
1322     }
1323 
1324     bool enableWildcard = true;
1325     for (size_t iTrig = 0; iTrig < triggerNames_.size(); ++iTrig) {
1326       // prepare for regular expression (with wildcards) functionality:
1327       TString tNameTmp = TString(triggerNames_[iTrig]);
1328       TRegexp tNamePattern = TRegexp(tNameTmp, enableWildcard);
1329       int tIndex = -1;
1330       // find the trigger index:
1331       for (unsigned ipath = 0; ipath < hltConfig_.size(); ++ipath) {
1332         // use TString since it provides reg exp functionality:
1333         TString tmpName = TString(hltConfig_.triggerName(ipath));
1334         if (tmpName.Contains(tNamePattern)) {
1335           tIndex = int(ipath);
1336           triggerIndices_.push_back(tIndex);
1337         }
1338       }
1339       if (tIndex < 0) {  // if can't find trigger path at all, give warning:
1340         std::cout << "Warning: Could not find trigger" << triggerNames_[iTrig] << std::endl;
1341         //assert(false);
1342       }
1343     }  // end for triggerNames
1344     for (size_t iTrig = 0; iTrig < isoTriggerNames_.size(); ++iTrig) {
1345       // prepare for regular expression functionality:
1346       TString tNameTmp = TString(isoTriggerNames_[iTrig]);
1347       TRegexp tNamePattern = TRegexp(tNameTmp, enableWildcard);
1348       int tIndex = -1;
1349       // find the trigger index:
1350       for (unsigned ipath = 0; ipath < hltConfig_.size(); ++ipath) {
1351         // use TString since it provides reg exp functionality:
1352         TString tmpName = TString(hltConfig_.triggerName(ipath));
1353         if (tmpName.Contains(tNamePattern)) {
1354           tIndex = int(ipath);
1355           isoTriggerIndices_.push_back(tIndex);
1356         }
1357       }
1358       if (tIndex < 0) {  // if can't find trigger path at all, give warning:
1359         std::cout << "Warning: Could not find trigger" << isoTriggerNames_[iTrig] << std::endl;
1360         //assert(false);
1361       }
1362     }  // end for isoTriggerNames
1363   }  // end if (triggerMatching_)
1364 }
1365 
1366 void L1MuonRecoTreeProducer::endRun(const edm::Run &run, const edm::EventSetup &eventSetup) {}
1367 // ------------ method called once each job just before starting event loop  ------------
1368 void L1MuonRecoTreeProducer::beginJob(void) {}
1369 
1370 // ------------ method called once each job just after ending the event loop  ------------
1371 void L1MuonRecoTreeProducer::endJob() { delete muon; }
1372 
1373 //define this as a plug-in
1374 DEFINE_FWK_MODULE(L1MuonRecoTreeProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!