diff --git a/DataFormats/L1TrackTrigger/interface/L1TkMuonParticle.h b/DataFormats/L1TrackTrigger/interface/L1TkMuonParticle.h index 00381625cf427..39cc7e970d1e1 100644 --- a/DataFormats/L1TrackTrigger/interface/L1TkMuonParticle.h +++ b/DataFormats/L1TrackTrigger/interface/L1TkMuonParticle.h @@ -22,6 +22,8 @@ #include "DataFormats/L1Trigger/interface/L1MuonParticle.h" #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h" +#include "DataFormats/L1Trigger/interface/L1MuonParticleExtended.h" +#include "DataFormats/L1Trigger/interface/L1MuonParticleExtendedFwd.h" namespace l1extra { @@ -32,7 +34,7 @@ namespace l1extra typedef TTTrack< Ref_PixelDigi_ > L1TkTrackType; typedef std::vector< L1TkTrackType > L1TkTrackCollectionType; - L1TkMuonParticle(); + L1TkMuonParticle() : theIsolation(-999.), TrkzVtx_(999.), quality_(999) {} L1TkMuonParticle( const LorentzVector& p4, const edm::Ptr< DTMatch >& muRef, float tkisol = -999. ); @@ -42,6 +44,9 @@ namespace l1extra const edm::Ptr< L1TkTrackType >& trkPtr, float tkisol = -999. ); + //! more basic constructor, in case refs/ptrs can't be set or to be set separately + L1TkMuonParticle(const reco::LeafCandidate& cand) : reco::LeafCandidate(cand), theIsolation(-999.), TrkzVtx_(999.), quality_(999) {} + virtual ~L1TkMuonParticle() {} // const member functions @@ -50,16 +55,27 @@ namespace l1extra const edm::Ptr< L1TkTrackType >& getTrkPtr() const { return trkPtr_ ; } + const L1MuonParticleRef& getMuRef() const + { return muRef_ ; } + + const L1MuonParticleExtendedRef& getMuExtendedRef() const + { return muExtendedRef_ ; } + float getTrkIsol() const { return theIsolation; } float getTrkzVtx() const { return TrkzVtx_ ; } int bx() const ; - void setTrkzVtx(float TrkzVtx) { TrkzVtx_ = TrkzVtx ; } - void setTrkIsol(float TrkIsol) { theIsolation = TrkIsol ; } + unsigned int quality() const {return quality_;} - unsigned int quality() const; + void setTrkPtr(const edm::Ptr< L1TkTrackType >& p) {trkPtr_ = p;} + void setMuRef(const L1MuonParticleRef& r){muRef_ = r;} + void setMuExtendedRef(const L1MuonParticleExtendedRef& r){muExtendedRef_ = r;} + + void setTrkzVtx(float TrkzVtx) { TrkzVtx_ = TrkzVtx ; } + void setTrkIsol(float TrkIsol) { theIsolation = TrkIsol ; } + void setQuality(unsigned int q){ quality_ = q;} private: @@ -69,11 +85,12 @@ namespace l1extra // used for the Naive producer edm::Ref< L1MuonParticleCollection > muRef_ ; + L1MuonParticleExtendedRef muExtendedRef_ ; edm::Ptr< L1TkTrackType > trkPtr_ ; float theIsolation; float TrkzVtx_ ; - + unsigned int quality_; }; } diff --git a/DataFormats/L1TrackTrigger/src/L1TkMuonParticle.cc b/DataFormats/L1TrackTrigger/src/L1TkMuonParticle.cc index 65d6fa2183940..abb2b7bd86427 100644 --- a/DataFormats/L1TrackTrigger/src/L1TkMuonParticle.cc +++ b/DataFormats/L1TrackTrigger/src/L1TkMuonParticle.cc @@ -7,15 +7,15 @@ using namespace l1extra; -L1TkMuonParticle::L1TkMuonParticle() {} - // Padova's TkMuons L1TkMuonParticle::L1TkMuonParticle( const LorentzVector& p4, const edm::Ptr< DTMatch > &muRef, float tkisol ) : LeafCandidate( ( char ) 0, p4 ), theDTMatch ( muRef ) , - theIsolation ( tkisol ) + theIsolation ( tkisol ), + TrkzVtx_(999.), + quality_(999) { // need to set the z of the track @@ -30,14 +30,18 @@ L1TkMuonParticle::L1TkMuonParticle( const LorentzVector& p4, : LeafCandidate( ( char ) 0, p4 ), muRef_ ( muRef ) , trkPtr_ ( trkPtr ) , - theIsolation ( tkisol ) - + theIsolation ( tkisol ), + TrkzVtx_(999), + quality_(999) { if ( trkPtr_.isNonnull() ) { float z = getTrkPtr() -> getPOCA().z(); setTrkzVtx( z ); } + if (muRef_.isNonnull()) { + quality_ = muRef_->gmtMuonCand().quality(); + } } @@ -48,26 +52,19 @@ int L1TkMuonParticle::bx() const { return theDTMatch->getDTBX(); } - if (muRef_.isNonnull() ) { - return dummy; + // PL. In case Pierluigi's objects have a bx + if ( muRef_.isNonnull() ) { + return (getMuRef() -> bx()) ; + } + else if (muExtendedRef_.isNonnull()) { + return getMuExtendedRef()->bx(); + } return dummy; } -unsigned int L1TkMuonParticle::quality() const { - - int dummy = 999; - - if ( muRef_.isNonnull() ) { - return (muRef_ -> gmtMuonCand().quality() ); - } - else { - return dummy; - } -} - diff --git a/DataFormats/L1TrackTrigger/src/classes_def.xml b/DataFormats/L1TrackTrigger/src/classes_def.xml index ea6a351143da2..12e7b9ef72a2d 100644 --- a/DataFormats/L1TrackTrigger/src/classes_def.xml +++ b/DataFormats/L1TrackTrigger/src/classes_def.xml @@ -84,7 +84,8 @@ - + + diff --git a/DataFormats/L1Trigger/interface/L1MuonParticleExtended.h b/DataFormats/L1Trigger/interface/L1MuonParticleExtended.h new file mode 100644 index 0000000000000..cba01c898844e --- /dev/null +++ b/DataFormats/L1Trigger/interface/L1MuonParticleExtended.h @@ -0,0 +1,107 @@ +#ifndef L1Trigger_L1MuonParticleExtended_h +#define L1Trigger_L1MuonParticleExtended_h +// -*- C++ -*- +// +// Package: L1Trigger +// Class : L1MuonParticleExtended +// +// Description: L1MuonParticle with detector layer information. +// Should be useful from the start for detailed analysis of efficiency etc +// Expect this to become feasible for future upgrade hardware. + +// user include files +#include "DataFormats/L1Trigger/interface/L1MuonParticle.h" +#include "DataFormats/Candidate/interface/LeafCandidate.h" + +#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuRegionalCand.h" +#include "DataFormats/DetId/interface/DetId.h" + +// forward declarations + +namespace l1extra { + + class L1MuonParticleExtended : public L1MuonParticle { + public: + L1MuonParticleExtended() : sigmaEta_(-99), sigmaPhi_(-99), quality_(0) {} + + L1MuonParticleExtended( const L1MuonParticle& l1mu) : L1MuonParticle(l1mu), sigmaEta_(-99), sigmaPhi_(-99), quality_(0) {} + + virtual ~L1MuonParticleExtended() {} + + + struct StationData { + StationData() : station(-1), ringOrWheel(-99), bx(-99), quality(-99), + phi(-99), sigmaPhi(-99), + eta(-99), sigmaEta(-99), + bendPhi(-99), bendEta(-99), + bendPhiInt(-99), bendEtaInt(-99), + valid(false) {} + DetId id; + int station; //a bit wasteful + int ringOrWheel; + int bx; + int quality; + float phi; + float sigmaPhi; + float eta; + float sigmaEta; + float bendPhi; //! direction.phi - position.phi + float bendEta; //! direction.eta - position.eta + + int bendPhiInt; //magic word from primitives about phi direction + int bendEtaInt; //magic word from primitives about eta or theta direction + + bool valid; + }; + + // ---------- const member functions --------------------- + virtual L1MuonParticleExtended* clone() const + { return new L1MuonParticleExtended( *this ) ; } + + // ---------- member functions --------------------------- + const L1MuRegionalCand& cscCand() const {return cscCand_;} + void setCscCand(const L1MuRegionalCand& cand) {cscCand_ = cand;} + + const L1MuRegionalCand& rpcCand() const {return rpcCand_;} + void setRpcCand(const L1MuRegionalCand& cand) {rpcCand_ = cand;} + + const L1MuRegionalCand& dtCand() const {return dtCand_;} + void setDtCand(const L1MuRegionalCand& cand) {dtCand_ = cand;} + + float sigmaEta() const { return sigmaEta_; } + void setSigmaEta(float val) {sigmaEta_ = val;} + + float sigmaPhi() const { return sigmaPhi_; } + void setSigmaPhi(float val) {sigmaPhi_ = val;} + + unsigned int quality() const {return quality_;} + void setQuality(unsigned int q) {quality_ = q;} + + const StationData& cscData(unsigned int s ) const {return (s>0 && s<= 4) ? cscData_[s-1] : dummyData_; } + void setCscData(const StationData& sd, unsigned int s) {if (s>0 && s<= 4) cscData_[s-1] = sd; } + + const StationData& dtData(unsigned int s ) const {return (s>0 && s<= 4) ? dtData_[s-1] : dummyData_; } + void setDtData(const StationData& sd, unsigned int s) {if (s>0 && s<= 4) dtData_[s-1] = sd; } + + const StationData& rpcData(unsigned int s ) const {return (s>0 && s<= 4) ? rpcData_[s-1] : dummyData_; } + void setRpcData(const StationData& sd, unsigned int s) {if (s>0 && s<= 4) rpcData_[s-1] = sd; } + private: + // ---------- member data -------------------------------- + L1MuRegionalCand cscCand_; + L1MuRegionalCand rpcCand_; + L1MuRegionalCand dtCand_; + + float sigmaEta_; + float sigmaPhi_; + + unsigned int quality_; + + StationData cscData_[4]; + StationData dtData_[4]; + StationData rpcData_[4]; + + StationData dummyData_; + }; +} + +#endif diff --git a/DataFormats/L1Trigger/interface/L1MuonParticleExtendedFwd.h b/DataFormats/L1Trigger/interface/L1MuonParticleExtendedFwd.h new file mode 100644 index 0000000000000..5605c7b039beb --- /dev/null +++ b/DataFormats/L1Trigger/interface/L1MuonParticleExtendedFwd.h @@ -0,0 +1,25 @@ +#ifndef L1Trigger_L1MuonParticleExtendedFwd_h +#define L1Trigger_L1MuonParticleExtendedFwd_h +// -*- C++ -*- +// +// Package: L1Trigger +// Class : L1MuonParticleExtendedFwd +// + +#include +#include "DataFormats/Common/interface/Ref.h" +#include "DataFormats/Common/interface/RefVector.h" + + +namespace l1extra { + + class L1MuonParticleExtended ; + + typedef std::vector< L1MuonParticleExtended > L1MuonParticleExtendedCollection ; + + typedef edm::Ref< L1MuonParticleExtendedCollection > L1MuonParticleExtendedRef ; + typedef edm::RefVector< L1MuonParticleExtendedCollection > L1MuonParticleExtendedRefVector ; + typedef std::vector< L1MuonParticleExtendedRef > L1MuonParticleExtendedVectorRef ; +} + +#endif diff --git a/DataFormats/L1Trigger/src/classes.h b/DataFormats/L1Trigger/src/classes.h index fac481f156227..9dfa6b55a4182 100644 --- a/DataFormats/L1Trigger/src/classes.h +++ b/DataFormats/L1Trigger/src/classes.h @@ -11,6 +11,8 @@ #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h" #include "DataFormats/L1Trigger/interface/L1MuonParticle.h" #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h" +#include "DataFormats/L1Trigger/interface/L1MuonParticleExtended.h" +#include "DataFormats/L1Trigger/interface/L1MuonParticleExtendedFwd.h" #include "DataFormats/L1Trigger/interface/L1EtMissParticle.h" #include "DataFormats/L1Trigger/interface/L1EtMissParticleFwd.h" #include "DataFormats/L1Trigger/interface/L1ParticleMap.h" @@ -29,6 +31,7 @@ namespace { l1extra::L1EmParticleCollection emColl ; l1extra::L1JetParticleCollection jetColl ; l1extra::L1MuonParticleCollection muonColl ; + l1extra::L1MuonParticleExtendedCollection muonCollE ; l1extra::L1EtMissParticle etMiss ; l1extra::L1EtMissParticleCollection etMissColl ; l1extra::L1ParticleMapCollection mapColl ; @@ -37,6 +40,7 @@ namespace { edm::Wrapper w_emColl; edm::Wrapper w_jetColl; edm::Wrapper w_muonColl; + edm::Wrapper w_muonCollE; edm::Wrapper w_etMiss; edm::Wrapper w_etMissColl; edm::Wrapper w_mapColl; @@ -45,18 +49,21 @@ namespace { l1extra::L1EmParticleRef refEm ; l1extra::L1JetParticleRef refJet ; l1extra::L1MuonParticleRef refMuon ; + l1extra::L1MuonParticleExtendedRef refMuonE ; l1extra::L1EtMissParticleRef refEtMiss ; l1extra::L1HFRingsRef refHFRings ; l1extra::L1EmParticleRefVector refVecEmColl ; l1extra::L1JetParticleRefVector refVecJetColl ; l1extra::L1MuonParticleRefVector refVecMuonColl ; + l1extra::L1MuonParticleExtendedRefVector refVecMuonCollE ; l1extra::L1EtMissParticleRefVector refVecEtMiss ; l1extra::L1HFRingsRefVector refVecHFRings ; l1extra::L1EmParticleVectorRef vecRefEmColl ; l1extra::L1JetParticleVectorRef vecRefJetColl ; l1extra::L1MuonParticleVectorRef vecRefMuonColl ; + l1extra::L1MuonParticleExtendedVectorRef vecRefMuonCollE ; l1extra::L1EtMissParticleVectorRef vecRefEtMissColl ; l1extra::L1HFRingsVectorRef vecRefHFRingsColl ; @@ -64,6 +71,7 @@ namespace { edm::reftobase::Holder rtbe; edm::reftobase::Holder rtbm; + edm::reftobase::Holder rtbmE; edm::reftobase::Holder rtbj; edm::reftobase::Holder rtbm1; edm::reftobase::Holder rtbm2; diff --git a/DataFormats/L1Trigger/src/classes_def.xml b/DataFormats/L1Trigger/src/classes_def.xml index 522304c763cc6..f8eff615bf793 100644 --- a/DataFormats/L1Trigger/src/classes_def.xml +++ b/DataFormats/L1Trigger/src/classes_def.xml @@ -8,6 +8,12 @@ + + + + + + @@ -20,12 +26,14 @@ + + @@ -34,18 +42,21 @@ + + + @@ -55,6 +66,8 @@ + + diff --git a/L1Trigger/L1ExtraFromDigis/BuildFile.xml b/L1Trigger/L1ExtraFromDigis/BuildFile.xml index f74cab6b5bf89..95c472fedfbd8 100644 --- a/L1Trigger/L1ExtraFromDigis/BuildFile.xml +++ b/L1Trigger/L1ExtraFromDigis/BuildFile.xml @@ -4,7 +4,10 @@ + + + diff --git a/L1Trigger/L1ExtraFromDigis/python/l1extraMuExtended_cfi.py b/L1Trigger/L1ExtraFromDigis/python/l1extraMuExtended_cfi.py new file mode 100644 index 0000000000000..bac9c68fc222c --- /dev/null +++ b/L1Trigger/L1ExtraFromDigis/python/l1extraMuExtended_cfi.py @@ -0,0 +1,8 @@ +import FWCore.ParameterSet.Config as cms + +l1extraMuExtended = cms.EDProducer( + "L1MuonParticleExtendedProducer", + gmtROSource = cms.InputTag("simGmtDigis"), + csctfSource = cms.InputTag("simCsctfTrackDigis"), + writeAllCSCTFs = cms.bool(True) + ) diff --git a/L1Trigger/L1ExtraFromDigis/src/L1MuonParticleExtendedProducer.cc b/L1Trigger/L1ExtraFromDigis/src/L1MuonParticleExtendedProducer.cc new file mode 100644 index 0000000000000..6ebdf1933ac6f --- /dev/null +++ b/L1Trigger/L1ExtraFromDigis/src/L1MuonParticleExtendedProducer.cc @@ -0,0 +1,317 @@ +// -*- C++ -*- +// +// Class: L1MuonParticleExtendedProducer +// +// + + +// system include files +#include + +// user include files +#include "Geometry/Records/interface/MuonGeometryRecord.h" +#include "Geometry/CSCGeometry/interface/CSCGeometry.h" + +#include "CondFormats/L1TObjects/interface/L1MuTriggerScales.h" +#include "CondFormats/DataRecord/interface/L1MuTriggerScalesRcd.h" +#include "CondFormats/L1TObjects/interface/L1MuTriggerPtScale.h" +#include "CondFormats/DataRecord/interface/L1MuTriggerPtScaleRcd.h" +#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutCollection.h" +#include "DataFormats/L1CSCTrackFinder/interface/L1CSCTrackCollection.h" + +#include "DataFormats/L1Trigger/interface/L1MuonParticleExtended.h" +#include "DataFormats/L1Trigger/interface/L1MuonParticleExtendedFwd.h" + +#include "L1Trigger/CSCCommonTrigger/interface/CSCConstants.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "DataFormats/Common/interface/OrphanHandle.h" + +#include "DataFormats/Math/interface/deltaPhi.h" + +// +// class declaration +// +class L1MuonParticleExtendedProducer : public edm::EDProducer { +public: + explicit L1MuonParticleExtendedProducer(const edm::ParameterSet&); + ~L1MuonParticleExtendedProducer() {} + +private: + virtual void produce(edm::Event&, const edm::EventSetup&); + + edm::InputTag gmtROSource_ ; + edm::InputTag csctfSource_ ; + + bool writeAllCSCTFs_; + + static const double muonMassGeV_ ; + +}; + + +// +// static data member definitions +// + +const double L1MuonParticleExtendedProducer::muonMassGeV_ = 0.105658369 ; // PDG06 + +// +// constructors and destructor +// +L1MuonParticleExtendedProducer::L1MuonParticleExtendedProducer(const edm::ParameterSet& iConfig) + : gmtROSource_( iConfig.getParameter< edm::InputTag >("gmtROSource" ) ), + csctfSource_(iConfig.getParameter< edm::InputTag >("csctfSource")), + writeAllCSCTFs_(iConfig.getParameter< bool>("writeAllCSCTFs")) +{ + using namespace l1extra ; + + //register your products + produces< L1MuonParticleExtendedCollection >() ; + if (writeAllCSCTFs_) produces< L1MuonParticleExtendedCollection >("csc") ; +} + + +// ------------ method called to produce the data ------------ +void +L1MuonParticleExtendedProducer::produce( edm::Event& iEvent, + const edm::EventSetup& iSetup) +{ + using namespace edm ; + using namespace l1extra ; + using namespace std ; + using namespace reco ; + + auto_ptr< L1MuonParticleExtendedCollection > muColl( new L1MuonParticleExtendedCollection ); + auto_ptr< L1MuonParticleExtendedCollection > cscColl( new L1MuonParticleExtendedCollection ); + + ESHandle< L1MuTriggerScales > muScales ; + iSetup.get< L1MuTriggerScalesRcd >().get( muScales ) ; + + ESHandle< L1MuTriggerPtScale > muPtScale ; + iSetup.get< L1MuTriggerPtScaleRcd >().get( muPtScale ) ; + + Handle< L1MuGMTReadoutCollection > gmtRO_h ; + iEvent.getByLabel( gmtROSource_, gmtRO_h ) ; + const L1MuGMTReadoutCollection& gmtROs(*gmtRO_h.product()); + + Handle csctftH; + iEvent.getByLabel(csctfSource_, csctftH); + const L1CSCTrackCollection& csctfts(*csctftH.product()); + + ESHandle cscGH; + iSetup.get().get(cscGH); + const CSCGeometry* csc_geometry = &*cscGH; + + + int curMuCollIdx = -1; //index of the GMT-based + for (auto gmtRO : gmtROs.getRecords()){ + vector cscGmtIdxs(4,-1); //index into muColl for a given CSC cand + for (auto const gmtCand : gmtRO.getGMTCands() ){ + if (gmtCand.empty()) continue; + + float pt = muPtScale->getPtScale()->getLowEdge( gmtCand.ptIndex() ) + 1.e-6 ; + float eta = muScales->getGMTEtaScale()->getCenter( gmtCand.etaIndex() ) ; + float phi = muScales->getPhiScale()->getLowEdge( gmtCand.phiIndex() ) ; + + + float etaP = muScales->getGMTEtaScale()->getCenter( gmtCand.etaIndex() +1 ) ; + float phiP = muScales->getPhiScale()->getLowEdge( gmtCand.phiIndex() +1 ) ; + + float etaM = muScales->getGMTEtaScale()->getCenter( gmtCand.etaIndex() -1) ; + float phiM = muScales->getPhiScale()->getLowEdge( gmtCand.phiIndex() -1) ; + + math::PtEtaPhiMLorentzVector p4( pt, eta, phi, muonMassGeV_); + L1MuonParticle l1muP(gmtCand.charge(), p4, gmtCand, gmtCand.bx()); + int cscGmtIdx = -1; + if (!gmtCand.isRPC() && gmtCand.isFwd()) cscGmtIdx = gmtCand.getDTCSCIndex(); + + muColl->push_back(L1MuonParticleExtended(l1muP)); + curMuCollIdx++; + if (cscGmtIdx>=0) cscGmtIdxs[cscGmtIdx] = curMuCollIdx; + L1MuonParticleExtended* aGmtCand = &muColl->back(); + + float sigmaEta=99; + float dEtaP = std::abs(etaP - eta); + float dEtaM = std::abs(etaM - eta); + if (dEtaP < 1 && dEtaM < 1 && dEtaP > 0 && dEtaM > 0){ + sigmaEta = sqrt(dEtaP*dEtaM); //take geom mean for no particular reason + } else if (dEtaP < 1 || dEtaM < 1) { + sigmaEta = dEtaP < dEtaM && dEtaP >0 ? dEtaP : dEtaM; + } + sigmaEta = sigmaEta/sqrt(12.); + + float sigmaPhi=99; + float dPhiP = std::abs(deltaPhi(phiP, phi)); + float dPhiM = std::abs(deltaPhi(phiM, phi)); + if (dPhiP < 1 && dPhiM < 1 && dPhiP > 0 && dPhiM > 0){ + sigmaPhi = sqrt(dPhiP*dPhiM); //take geom mean for no particular reason + } else if (dPhiP < 1 || dPhiM < 1) { + sigmaPhi = dPhiP < dPhiM && dPhiP > 0 ? dPhiP : dPhiM; + } + sigmaPhi = sigmaPhi/sqrt(12.); + + aGmtCand->setSigmaEta(sigmaEta); + aGmtCand->setSigmaPhi(sigmaPhi); + aGmtCand->setQuality(gmtCand.quality()); + + } + if (1< 2){ //at the moment this is not really "ALL", just up to 4 from CSCTF + unsigned cscInd = -1; + for (auto const csctfReg : gmtRO.getCSCCands() ){ + cscInd++; + if (csctfReg.empty()) continue; + + float pt = muPtScale->getPtScale()->getLowEdge( csctfReg.pt_packed() ) + 1.e-6 ; + float eta = muScales->getRegionalEtaScale(csctfReg.type_idx())->getCenter(csctfReg.eta_packed()); + float phi = muScales->getPhiScale()->getLowEdge( csctfReg.phi_packed() ) ; + + float etaP = muScales->getRegionalEtaScale(csctfReg.type_idx())->getCenter(csctfReg.eta_packed() + 1); + float phiP = muScales->getPhiScale()->getLowEdge( csctfReg.phi_packed() +1 ) ; + float etaM = muScales->getRegionalEtaScale(csctfReg.type_idx())->getCenter(csctfReg.eta_packed() -1); + float phiM = muScales->getPhiScale()->getLowEdge( csctfReg.phi_packed() -1 ) ; + + int gmtIdx = cscGmtIdxs[cscInd]; + L1MuonParticleExtended* gmtParticle = 0; + if (gmtIdx >=0) gmtParticle = &(*muColl)[gmtIdx]; + L1MuGMTExtendedCand gmtCand = gmtParticle != 0 ? gmtParticle->gmtMuonCand() : L1MuGMTExtendedCand(); + + math::PtEtaPhiMLorentzVector p4( pt, eta, phi, muonMassGeV_); + int qCand = csctfReg.chargeValid() ? csctfReg.chargeValue() : 0; + L1MuonParticle l1muP(qCand, p4, gmtCand, csctfReg.bx()); + + L1MuonParticleExtended cscParticleTmp(l1muP); + cscParticleTmp.setCscCand(csctfReg); + if (gmtParticle) gmtParticle->setCscCand(csctfReg); + cscColl->push_back(cscParticleTmp); + L1MuonParticleExtended* cscParticle = &cscColl->back(); + + float sigmaEta=99; + float dEtaP = std::abs(etaP - eta); + float dEtaM = std::abs(etaM - eta); + if (dEtaP < 1 && dEtaM < 1 && dEtaP > 0 && dEtaM > 0){ + sigmaEta = sqrt(dEtaP*dEtaM); //take geom mean for no particular reason + } else if (dEtaP < 1 || dEtaM < 1) { + sigmaEta = dEtaP < dEtaM && dEtaP >0 ? dEtaP : dEtaM; + } + sigmaEta = sigmaEta/sqrt(12.); + + float sigmaPhi=99; + float dPhiP = std::abs(deltaPhi(phiP, phi)); + float dPhiM = std::abs(deltaPhi(phiM, phi)); + if (dPhiP < 1 && dPhiM < 1 && dPhiP > 0 && dPhiM > 0){ + sigmaPhi = sqrt(dPhiP*dPhiM); //take geom mean for no particular reason + } else if (dPhiP < 1 || dPhiM < 1) { + sigmaPhi = dPhiP < dPhiM && dPhiP > 0 ? dPhiP : dPhiM; + } + sigmaPhi = sigmaPhi/sqrt(12.); + + cscParticle->setSigmaEta(sigmaEta); + cscParticle->setSigmaPhi(sigmaPhi); + + unsigned int qual = 0; + if (gmtParticle) qual = gmtParticle->quality(); + else { + unsigned int qualC = csctfReg.quality(); + if (qualC ==2) qual = 3; + if (qualC ==3) qual = 6; + } + cscParticle->setQuality(qual); + + //now get the details, not the most optimal way .. not a big deal here + unsigned int roWord = csctfReg.getDataWord(); + for (auto csctfFull : csctfts){ + auto cscL1Tk = csctfFull.first; + unsigned int fWord = cscL1Tk.getDataWord(); + //there's probably a more obvious way to match (just by local phi,sector,eta,etc) + if (cscL1Tk.endcap() == 2) fWord |= (1<< (L1MuRegionalCand::ETA_START+L1MuRegionalCand::ETA_LENGTH-1)); + unsigned int fQuality = 0; + unsigned int fPt = 0; + cscL1Tk.decodeRank(cscL1Tk.rank(), fPt, fQuality); + fWord |= (fQuality<< L1MuRegionalCand::QUAL_START); + fWord |= (fPt << L1MuRegionalCand::PT_START); + //from L1Trigger/CSCTrackFinder/src/CSCTFMuonSorter.cc + unsigned int fgPhi = cscL1Tk.localPhi() + (cscL1Tk.sector() -1)*24 + 6; + if(fgPhi > 143) fgPhi -= 143; + fWord |= (fgPhi << L1MuRegionalCand::PHI_START); + // edm::LogWarning("MYDEBUG")<<"Checking or vs other "<4 || aStation<1) { + edm::LogError("BADLCTID")<<"got station "<chamber(anID)->layer(CSCConstants::KEY_CLCT_LAYER)->geometry(); + + typedef L1MuonParticleExtended::StationData StationData; + StationData sData; + for (; aDigi != anLCT.second.second; ++aDigi){ + sData.valid = true; + sData.id = anID.rawId(); + sData.station = aStation; + sData.ringOrWheel = aRing; + sData.bx = aDigi->getBX()-6; //offset to be centered at zero as all higher level objects have it + sData.quality = aDigi->getQuality(); + + int aWG = aDigi->getKeyWG(); + int aHS = aDigi->getStrip(); + float fStrip = 0.5*(aHS+1) - 0.25 ; //mind the offset + float fWire = layer_geometry->middleWireOfGroup(aWG+1); + LocalPoint wireStripPoint = layer_geometry->intersectionOfStripAndWire(fStrip, fWire); + CSCDetId key_id(anID.endcap(), anID.station(), anID.ring(), anID.chamber(), CSCConstants::KEY_CLCT_LAYER); + GlobalPoint gp = csc_geometry->idToDet(key_id)->surface().toGlobal(wireStripPoint); + sData.phi = gp.phi(); + sData.eta = gp.eta(); + + int deltaWire = aWG +1 < layer_geometry->numberOfWireGroups() ? 1 : -1; + int deltaStrip = fStrip +1 <= layer_geometry->numberOfStrips() ? 1 : -1; + + float fStripD1 = fStrip+0.5*deltaStrip; + float fWireD1 = layer_geometry->middleWireOfGroup(aWG+1 + deltaWire); + LocalPoint wireStripPlus1 = layer_geometry->intersectionOfStripAndWire(fStripD1, fWireD1); + GlobalPoint gpD1 = csc_geometry->idToDet(key_id)->surface().toGlobal(wireStripPlus1); + float dPhi = std::abs(deltaPhi(gp.phi(), gpD1.phi())); + + float dEta = fabs(gp.eta() - gpD1.eta()); + + if (dEta<1E-3){ + edm::LogWarning("SmallDeltaEta")<<" st "<< aStation<<" r "<getCLCTPattern(); + if (aDigi->getBend()) sData.bendPhiInt*=-1; + } //LCTs for a given anID + + cscParticle->setCscData(sData, aStation); + if (gmtParticle) gmtParticle->setCscData(sData, aStation); + } //LCT range for anID + }//matched full CSCTF track collection with the current L1MuRegionalCand + }//loop over csctfts + } + } + } + + iEvent.put( muColl ); + if (writeAllCSCTFs_) iEvent.put( cscColl, "csc" ); +} + + +//define this as a plug-in +DEFINE_FWK_MODULE(L1MuonParticleExtendedProducer); diff --git a/SLHCUpgradeSimulations/L1TrackTrigger/plugins/L1TkMuonFromExtendedProducer.cc b/SLHCUpgradeSimulations/L1TrackTrigger/plugins/L1TkMuonFromExtendedProducer.cc new file mode 100644 index 0000000000000..80a2276df9778 --- /dev/null +++ b/SLHCUpgradeSimulations/L1TrackTrigger/plugins/L1TkMuonFromExtendedProducer.cc @@ -0,0 +1,325 @@ +// -*- C++ -*- +// +// input: L1 TkTracks and L1MuonParticleExtended (standalone with component details) +// match the two and produce a collection of L1TkMuonParticle +// eventually, this should be made modular and allow to swap out different algorithms + + +#include "DataFormats/L1TrackTrigger/interface/L1TkMuonParticle.h" +#include "DataFormats/L1TrackTrigger/interface/L1TkMuonParticleFwd.h" + + +// for L1Tracks: +#include "SimDataFormats/SLHC/interface/StackedTrackerTypes.h" + +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/EventSetup.h" + +#include "DataFormats/Common/interface/Handle.h" +#include "FWCore/Utilities/interface/InputTag.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/Math/interface/LorentzVector.h" +#include "DataFormats/Math/interface/deltaR.h" +#include "DataFormats/Math/interface/deltaPhi.h" + +#include "TMath.h" + +// system include files +#include +#include + + +using namespace l1extra ; + +// +// class declaration +// + +class L1TkMuonFromExtendedProducer : public edm::EDProducer { +public: + + typedef TTTrack< Ref_PixelDigi_ > L1TkTrackType; + typedef std::vector< L1TkTrackType > L1TkTrackCollectionType; + + struct PropState { //something simple, imagine it's hardware emulation + PropState() : + pt(-99), eta(-99), phi(-99), + sigmaPt(-99), sigmaEta(-99), sigmaPhi(-99), + valid(false) {} + float pt; + float eta; + float phi; + float sigmaPt; + float sigmaEta; + float sigmaPhi; + bool valid; + + }; + + explicit L1TkMuonFromExtendedProducer(const edm::ParameterSet&); + ~L1TkMuonFromExtendedProducer(); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + virtual void produce(edm::Event&, const edm::EventSetup&); + + PropState propagateToGMT(const L1TkTrackType& l1tk) const; + + //configuration (preserving structure of L1TkMuonFromExtendedProducer + edm::InputTag L1MuonsInputTag_; + edm::InputTag L1TrackInputTag_; + + float ETAMIN_; + float ETAMAX_; + float ZMAX_; // |z_track| < ZMAX in cm + float CHI2MAX_; + float PTMINTRA_; + // float DRmax_; + int nStubsmin_ ; // minimum number of stubs + // bool closest_ ; + bool correctGMTPropForTkZ_; + +} ; + + +// +// constructors and destructor +// +L1TkMuonFromExtendedProducer::L1TkMuonFromExtendedProducer(const edm::ParameterSet& iConfig) +{ + + + L1MuonsInputTag_ = iConfig.getParameter("L1MuonsInputTag"); + L1TrackInputTag_ = iConfig.getParameter("L1TrackInputTag"); + + ETAMIN_ = (float)iConfig.getParameter("ETAMIN"); + ETAMAX_ = (float)iConfig.getParameter("ETAMAX"); + ZMAX_ = (float)iConfig.getParameter("ZMAX"); + CHI2MAX_ = (float)iConfig.getParameter("CHI2MAX"); + PTMINTRA_ = (float)iConfig.getParameter("PTMINTRA"); + // DRmax_ = (float)iConfig.getParameter("DRmax"); + nStubsmin_ = iConfig.getParameter("nStubsmin"); + // closest_ = iConfig.getParameter("closest"); + + correctGMTPropForTkZ_ = iConfig.getParameter("correctGMTPropForTkZ"); + + produces(); +} + +L1TkMuonFromExtendedProducer::~L1TkMuonFromExtendedProducer() { +} + +// ------------ method called to produce the data ------------ +void +L1TkMuonFromExtendedProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + using namespace edm; + using namespace std; + + + std::auto_ptr tkMuons(new L1TkMuonParticleCollection); + + edm::Handle l1musH; + iEvent.getByLabel(L1MuonsInputTag_, l1musH); + const L1MuonParticleExtendedCollection& l1mus = *l1musH.product(); + + edm::Handle l1tksH; + iEvent.getByLabel(L1TrackInputTag_, l1tksH); + const L1TkTrackCollectionType& l1tks = *l1tksH.product(); + + L1TkMuonParticleCollection l1tkmuCands; + l1tkmuCands.reserve(l1mus.size()*4); //can do more if really needed + + int imu = -1; + for (auto l1mu : l1mus){ + imu++; + L1MuonParticleExtendedRef l1muRef(l1musH, imu); + + float l1mu_eta = l1mu.eta(); + float l1mu_phi = l1mu.phi(); + + float l1mu_feta = fabs( l1mu_eta ); + if (l1mu_feta < ETAMIN_) continue; + if (l1mu_feta > ETAMAX_) continue; + + // can skip quality cuts at the moment, keep bx=0 req + if (l1mu.bx() != 0) continue; + + unsigned int l1mu_quality = l1mu.quality(); + + const auto& gmtCand = l1mu.gmtMuonCand(); + if (!gmtCand.empty()){ + //some selections here + //currently the input can be a merge from different track finders + //so, the GMT cand may be missing + } + + const auto& dtCand = l1mu.dtCand(); + if (!dtCand.empty()){ + // something can be called from here + } + + const auto& cscCand = l1mu.cscCand(); + if (!cscCand.empty()){ + //apply something specific here + } + + float drmin = 999; + float ptmax = -1; + if (ptmax < 0) ptmax = -1; // dummy + + PropState matchProp; + int match_idx = -1; + int il1tk = -1; + + LogDebug("MDEBUG")<<"have a gmt, look for a match "; + for (auto l1tk : l1tks ){ + il1tk++; + float l1tk_pt = l1tk.getMomentum().perp(); + if (l1tk_pt < PTMINTRA_) continue; + + float l1tk_z = l1tk.getPOCA().z(); + if (fabs(l1tk_z) > ZMAX_) continue; + + float l1tk_chi2 = l1tk.getChi2(); + if (l1tk_chi2 > CHI2MAX_) continue; + + int l1tk_nstubs = l1tk.getStubRefs().size(); + if ( l1tk_nstubs < nStubsmin_) continue; + + float l1tk_eta = l1tk.getMomentum().eta(); + float l1tk_phi = l1tk.getMomentum().phi(); + + + float dr2 = deltaR2(l1mu_eta, l1mu_phi, l1tk_eta, l1tk_phi); + + if (dr2 > 0.3) continue; + + PropState pstate = propagateToGMT(l1tk); + if (!pstate.valid) continue; + + float dr2prop = deltaR2(l1mu_eta, l1mu_phi, pstate.eta, pstate.phi); + + if (dr2prop < drmin){ + drmin = dr2prop; + match_idx = il1tk; + matchProp = pstate; + } + }// over l1tks + + LogDebug("MYDEBUG")<<"matching index is "<= 0){ + const L1TkTrackType& matchTk = l1tks[match_idx]; + + float etaCut = 3.*sqrt(l1mu.sigmaEta()*l1mu.sigmaEta() + matchProp.sigmaEta*matchProp.sigmaEta); + float phiCut = 4.*sqrt(l1mu.sigmaPhi()*l1mu.sigmaPhi() + matchProp.sigmaPhi*matchProp.sigmaPhi); + + float dEta = std::abs(matchProp.eta - l1mu.eta()); + float dPhi = std::abs(deltaPhi(matchProp.phi, l1mu.phi())); + + LogDebug("MYDEBUG")<<"match details: prop "< l1tkPtr(l1tksH, match_idx); + auto p3 = matchTk.getMomentum(); + float p4e = sqrt(0.105658369*0.105658369 + p3.mag()); + math::XYZTLorentzVector l1tkp4(p3.x(), p3.y(), p3.z(), p4e); + + auto tkv3=matchTk.getPOCA(); + math::XYZPoint v3(tkv3.x(), tkv3.y(), tkv3.z()); + float trkisol = -999; + int l1tk_q = matchTk.getRInv()>0? 1: -1; + + L1TkMuonParticle l1tkmu(reco::LeafCandidate(l1tk_q, l1tkp4, v3, -13*l1tk_q )); + l1tkmu.setTrkPtr(l1tkPtr); + l1tkmu.setMuExtendedRef(l1muRef); + l1tkmu.setQuality(l1mu_quality); + l1tkmu.setTrkIsol(trkisol); + tkMuons->push_back(l1tkmu); + } + } + + }//over l1mus + + + iEvent.put( tkMuons ); + +} + + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void +L1TkMuonFromExtendedProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + //The following says we do not know what parameters are allowed so do no validation + // Please change this to state exactly what you do use, even if it is no parameters + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + + +L1TkMuonFromExtendedProducer::PropState L1TkMuonFromExtendedProducer::propagateToGMT(const L1TkMuonFromExtendedProducer::L1TkTrackType& tk) const { + auto p3 = tk.getMomentum(); + float tk_pt = p3.perp(); + float tk_p = p3.mag(); + float tk_eta = p3.eta(); + float tk_aeta = std::abs(tk_eta); + float tk_phi = p3.phi(); + float tk_q = tk.getRInv()>0? 1.: -1.; + float tk_z = tk.getPOCA().z(); + if (!correctGMTPropForTkZ_) tk_z = 0; + + L1TkMuonFromExtendedProducer::PropState dest; + if (tk_p<3.5 ) return dest; + if (tk_aeta <1.1 && tk_pt < 3.5) return dest; + if (tk_aeta > 2.5) return dest; + + //0th order: + dest.valid = true; + + float dzCorrPhi = 1.; + float deta = 0; + float etaProp = tk_aeta; + + if (tk_aeta < 1.1){ + etaProp = 1.1; + deta = tk_z/550./cosh(tk_aeta); + } else { + float delta = tk_z/850.; //roughly scales as distance to 2nd station + if (tk_eta > 0) delta *=-1; + dzCorrPhi = 1. + delta; + + float zOzs = tk_z/850.; + if (tk_eta > 0) deta = zOzs/(1. - zOzs); + else deta = zOzs/(1.+zOzs); + deta = deta*tanh(tk_eta); + } + float resPhi = tk_phi - 1.464*tk_q*cosh(1.7)/cosh(etaProp)/tk_pt*dzCorrPhi - M_PI/144.; + if (resPhi > M_PI) resPhi -= 2.*M_PI; + if (resPhi < -M_PI) resPhi += 2.*M_PI; + + dest.eta = tk_eta + deta; + dest.phi = resPhi; + dest.pt = tk_pt; //not corrected for eloss + + dest.sigmaEta = 0.100/tk_pt; //multiple scattering term + dest.sigmaPhi = 0.106/tk_pt; //need a better estimate for these + return dest; +} + +//define this as a plug-in +DEFINE_FWK_MODULE(L1TkMuonFromExtendedProducer); + + + diff --git a/SLHCUpgradeSimulations/L1TrackTrigger/python/l1TkMuonsExt_cff.py b/SLHCUpgradeSimulations/L1TrackTrigger/python/l1TkMuonsExt_cff.py new file mode 100644 index 0000000000000..56e44adbfa292 --- /dev/null +++ b/SLHCUpgradeSimulations/L1TrackTrigger/python/l1TkMuonsExt_cff.py @@ -0,0 +1,14 @@ +import FWCore.ParameterSet.Config as cms + +from SLHCUpgradeSimulations.L1TrackTrigger.l1TkMuonsExt_cfi import l1TkMuonsExt +l1TkMuonsExtCSC = l1TkMuonsExt.clone( + L1MuonsInputTag = cms.InputTag("l1extraMuExtended", "csc"), + ETAMIN = cms.double(1.1) + ) + +l1TkMuonsExtNoZCor = l1TkMuonsExt.clone( correctGMTPropForTkZ = cms.bool(False) ) +l1TkMuonsExtCSCNoZCor = l1TkMuonsExtCSC.clone( correctGMTPropForTkZ = cms.bool(False) ) + +l1TkMuonsExtSequence = cms.Sequence (l1TkMuonsExt * l1TkMuonsExtCSC + * l1TkMuonsExtNoZCor * l1TkMuonsExtCSCNoZCor) + diff --git a/SLHCUpgradeSimulations/L1TrackTrigger/python/l1TkMuonsExt_cfi.py b/SLHCUpgradeSimulations/L1TrackTrigger/python/l1TkMuonsExt_cfi.py new file mode 100644 index 0000000000000..7b1baad9c89d2 --- /dev/null +++ b/SLHCUpgradeSimulations/L1TrackTrigger/python/l1TkMuonsExt_cfi.py @@ -0,0 +1,19 @@ +import FWCore.ParameterSet.Config as cms + +l1TkMuonsExt = cms.EDProducer( + "L1TkMuonFromExtendedProducer", + L1MuonsInputTag = cms.InputTag("l1extraMuExtended"), + L1TrackInputTag = cms.InputTag("TTTracksFromPixelDigis", "Level1TTTracks"), + ETAMIN = cms.double(0), + ETAMAX = cms.double(5.), # no cut + ZMAX = cms.double( 25. ), # in cm + CHI2MAX = cms.double( 100. ), + PTMINTRA = cms.double( 2. ), # in GeV +# DRmax = cms.double( 0.5 ), + nStubsmin = cms.int32( 3 ), # minimum number of stubs +# closest = cms.bool( True ), + correctGMTPropForTkZ = cms.bool(True) + +) + +