diff --git a/DQM/L1TMonitor/interface/L1TkMuonMonitor.h b/DQM/L1TMonitor/interface/L1TkMuonMonitor.h new file mode 100755 index 0000000000000..71f4ced903d61 --- /dev/null +++ b/DQM/L1TMonitor/interface/L1TkMuonMonitor.h @@ -0,0 +1,89 @@ +#ifndef DQM_L1TMonitor_L1TkMuonMonitor_h +#define DQM_L1TMonitor_L1TkMuonMonitor_h + +/* + * \file L1TkMuonMonitor.h + * + * Date: 2015/06/16 12:00:00 + * Revision: 0.03 + * \author S.Baranov + * +*/ + +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DQMServices/Core/interface/DQMStore.h" +#include "DQMServices/Core/interface/MonitorElement.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "DataFormats/Candidate/interface/LeafCandidate.h" + +#include +#include +#include + +// +// class decleration +// + +class L1TkMuonMonitor : public edm::EDAnalyzer { + +public: + + // Constructor + L1TkMuonMonitor(const edm::ParameterSet& ps); + + // Destructor + virtual ~L1TkMuonMonitor(); + +protected: + // Analyze + void analyze(const edm::Event& e, const edm::EventSetup& c); + + // BeginJob + void beginJob(void); + + // BeginRun + void beginRun(const edm::Run& r, const edm::EventSetup& c); + + // EndJob + void endJob(void); + +private: + // ----------member data --------------------------- + DQMStore * dbe; + + MonitorElement *l1tkmu_Pt, *l1tkmu_Pt_barrel, *l1tkmu_Pt_endcap; + MonitorElement *l1tkmu_Eta, *l1tkmu_Eta_barrel, *l1tkmu_Eta_endcap; + MonitorElement *l1tkmu_Phi, *l1tkmu_Phi_barrel, *l1tkmu_Phi_endcap; + MonitorElement *l1tkmu_Dxy, *l1tkmu_Dxy_barrel, *l1tkmu_Dxy_endcap,; + MonitorElement *l1tkmu_Z, *l1tkmu_Z_barrel, *l1tkmu_Z_endcap; + MonitorElement *l1tkmu_Q, *l1tkmu_Q_barrel, *l1tkmu_Q_endcap; + MonitorElement *l1tkmu_Pdgid, *l1tkmu_Pdgid_barrel, *l1tkmu_Pdgid_endcap; + + int nev_; // Number of events processed + std::string outputFile_; //file name for ROOT ouput + bool verbose_; + bool monitorDaemon_; + ofstream logFile_; + edm::InputTag gmtSource_ ; + + static const double piconv_; + double phiconv_(float phi); + void book_(const edm::EventSetup& c); + + edm::InputTag candInputTag_; + +}; +#endif diff --git a/DQM/L1TMonitor/python/L1TMonitor_cff.py b/DQM/L1TMonitor/python/L1TMonitor_cff.py index 2e3f915b91b14..1bd522212a721 100644 --- a/DQM/L1TMonitor/python/L1TMonitor_cff.py +++ b/DQM/L1TMonitor/python/L1TMonitor_cff.py @@ -52,6 +52,9 @@ # GMT DQM module from DQM.L1TMonitor.L1TGMT_cfi import * +# MUTK DQM module +from DQM.L1TMonitor.L1MUTK_cfi import * + # GT DQM module from DQM.L1TMonitor.L1TGT_cfi import * @@ -110,6 +113,7 @@ l1tCsctf + l1tRpctf + l1tGmt + + l1MuTk + l1tGt + l1ExtraDqmSeq + l1tBPTX + diff --git a/DQM/L1TMonitor/python/L1TkMuonMonitor_cfi.py b/DQM/L1TMonitor/python/L1TkMuonMonitor_cfi.py new file mode 100644 index 0000000000000..42f2eee358c11 --- /dev/null +++ b/DQM/L1TMonitor/python/L1TkMuonMonitor_cfi.py @@ -0,0 +1,9 @@ +import FWCore.ParameterSet.Config as cms + +L1TkMuonMonitor = cms.EDAnalyzer("L1TkMuonMonitor", + disableROOToutput = cms.untracked.bool(True), + verbose = cms.untracked.bool(False), + gmtSource = cms.InputTag("gtDigis"), + DQMStore = cms.untracked.bool(True), + candInputTag = cms.InputTag("L1TkMuonsMerge") +) diff --git a/DQM/L1TMonitor/src/L1TkMuonMonitor.cc b/DQM/L1TMonitor/src/L1TkMuonMonitor.cc new file mode 100755 index 0000000000000..fc1e90521c626 --- /dev/null +++ b/DQM/L1TMonitor/src/L1TkMuonMonitor.cc @@ -0,0 +1,242 @@ +/* + * \file L1TkMuonMonitor.cc + * + * Date: 2015/06/16 12:00:00 + * Revision: 0.03 + * \author S.Baranov + * + */ + +#include "DQM/L1TMonitor/interface/L1TkMuonMonitor.h" +#include "DQMServices/Core/interface/DQMStore.h" + +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/ESHandle.h" + +using namespace std; +using namespace edm; +using namespace reco; + +typedef edm::View CandView; + +const double L1TkMuonMonitor::piconv_ = 180./M_PI; + +L1TkMuonMonitor::L1TkMuonMonitor(const ParameterSet& ps) + : gmtSource_( ps.getParameter< InputTag >("gmtSource") ) + { + + // verbosity switch + verbose_ = ps.getUntrackedParameter("verbose", false); + + if(verbose_) cout << "L1TkMuonMonitor: constructor...." << endl; + + dbe = NULL; + if ( ps.getUntrackedParameter("DQMStore", false) ) + { + dbe = Service().operator->(); + dbe->setVerbose(0); + } + + outputFile_ = ps.getUntrackedParameter("outputFile", ""); + if ( outputFile_.size() != 0 ) { + LogInfo("L1TkMuonMonitor") << "L1T Monitoring histograms will be saved to " << outputFile_.c_str(); + } + + bool disable = ps.getUntrackedParameter("disableROOToutput", false); + if(disable){ + outputFile_=""; + } + + // input tags + candInputTag_ = ps.getParameter("candInputTag"); + LogInfo("L1TkMuonMonitor") << " L1TkMuonMonitor::L1TkMuonMonitor candInputTag " << candInputTag_; + + if ( dbe !=NULL ) { + dbe->setCurrentFolder("L1T/L1TkMuon"); + } + +} + +L1TkMuonMonitor::~L1TkMuonMonitor() +{ +} + +void L1TkMuonMonitor::beginJob() +{ + nev_ = 0; +} + +void L1TkMuonMonitor::beginRun(const edm::Run& r, const edm::EventSetup& c) +{ + if(nev_==0) { + book_(c); + } +} + +void L1TkMuonMonitor::endJob(void) +{ + if(verbose_) cout << "L1TkMuonMonitor: end job...." << endl; + LogInfo("EndJob") << "analyzed " << nev_ << " events"; + + if ( outputFile_.size() != 0 && dbe ) dbe->save(outputFile_); + + return; +} + +void L1TkMuonMonitor::analyze(const Event& e, const EventSetup& c) +{ + + nev_++; + + edm::Handle aH; + e.getByLabel(candInputTag_, aH); + + if(!aH.isValid()){ + LogDebug("L1TkMuonMonitor") << " ::analyze Can't find CandView with label **candInputTag**"; + } + + const CandView& cands(*aH.product()); + + int track = 0, trackPt=0; + double maxPt = 0.; + double dxy; + for (auto cand : cands){ + track++; + if (cand.pt() >= maxPt){ + trackPt = track; + maxPt = cand.pt(); + } + } + + track = 0; + for (auto cand : cands){ + track++; + if ( track == trackPt ){ + dxy = cand.vx()*sin(cand.phi())-cand.vy()*cos(cand.phi()); + + LogDebug("L1TkMuonMonitor") << " ::analyze track " + << track << " Pt " << cand.pt() << " Phi " << cand.phi() << " Eta " << cand.eta() << " Dxy " << dxy + << " Z " << cand.vz() << " Q " << cand.charge() << " cand.pdgId " << cand.pdgId(); + + l1tkmu_Pt->Fill(cand.pt()); + l1tkmu_Phi->Fill(cand.phi()); + l1tkmu_Eta->Fill(cand.eta()); + l1tkmu_Dxy->Fill(dxy); + l1tkmu_Z->Fill(cand.vz()); + l1tkmu_Q->Fill(float(cand.charge())); + l1tkmu_Pdgid->Fill(float(cand.pdgId())); + + if ( std::abs(cand.eta()) <= 1.1 ){ + // Barrel + l1tkmu_Pt_barrel->Fill(cand.pt()); + l1tkmu_Phi_barrel->Fill(cand.phi()); + l1tkmu_Eta_barrel->Fill(cand.eta()); + l1tkmu_Dxy_barrel->Fill(dxy); + l1tkmu_Z_barrel->Fill(cand.vz()); + l1tkmu_Q_barrel->Fill(float(cand.charge())); + l1tkmu_Pdgid_barrel->Fill(float(cand.pdgId())); + } else { + // Endcap + l1tkmu_Pt_endcap->Fill(cand.pt()); + l1tkmu_Phi_endcap->Fill(cand.phi()); + l1tkmu_Eta_endcap->Fill(cand.eta()); + l1tkmu_Dxy_endcap->Fill(dxy); + l1tkmu_Z_endcap->Fill(cand.vz()); + l1tkmu_Q_endcap->Fill(float(cand.charge())); + l1tkmu_Pdgid_endcap->Fill(float(cand.pdgId())); + } + } + } +} + +void L1TkMuonMonitor::book_(const EventSetup& c) +{ + // get hold of back-end interface + DQMStore* dbe = 0; + dbe = Service().operator->(); + + if ( dbe ) { + dbe->setCurrentFolder("L1T/L1TkMuon"); + dbe->rmdir("L1T/L1TkMuon"); + } + + if ( dbe ) + { + dbe->setCurrentFolder("L1T/L1TkMuon"); + + /* + * All + */ + l1tkmu_Pt = dbe->book1D("l1tkmu_Pt","L1 Muon Pt", 100, 0.,250. ); + l1tkmu_Pt->setAxisTitle("l1tkmu_Pt",1); + + l1tkmu_Phi = dbe->book1D("l1tkmu_Phi","L1 Muon Phi", 100, -3.,3. ); + l1tkmu_Phi->setAxisTitle("l1tkmu_Phi",1); + + l1tkmu_Eta = dbe->book1D("l1tkmu_Eta","L1 Muon Eta", 100, -3.,3. ); + l1tkmu_Eta->setAxisTitle("l1tkmu_Eta",1); + + l1tkmu_Dxy = dbe->book1D("l1tkmu_Dxy","L1 Muon Dxy", 100,-15.,15. ); + l1tkmu_Dxy->setAxisTitle("l1tkmu_Dxy",1); + + l1tkmu_Z = dbe->book1D("l1tkmu_Z","L1 Muon Z", 100,-20.,20. ); + l1tkmu_Z->setAxisTitle("l1tkmu_Z",1); + + l1tkmu_Q = dbe->book1D("l1tkmu_Q","L1 Muon Q", 4, -2.,2. ); + l1tkmu_Q->setAxisTitle("l1tkmu_Q",1); + + l1tkmu_Pdgid = dbe->book1D("l1tkmu_Pdgid","L1 Muon Pdgid", 100,-50.,50. ); + l1tkmu_Pdgid->setAxisTitle("l1tkmu_Pdgid",1); + + /* + * "barrel" (|eta|<=1.1) + */ + l1tkmu_Pt_barrel = dbe->book1D("l1tkmu_Pt_barrel","L1 Muon Pt barrel", 100, 0.,250. ); + l1tkmu_Pt_barrel->setAxisTitle("l1tkmu_Pt_barrel",1); + + l1tkmu_Phi_barrel = dbe->book1D("l1tkmu_Phi_barrel","L1 Muon Phi barrel", 100, -3.,3. ); + l1tkmu_Phi_barrel->setAxisTitle("l1tkmu_Phi_barrel",1); + + l1tkmu_Eta_barrel = dbe->book1D("l1tkmu_Eta_barrel","L1 Muon Eta barrel", 100, -3.,3. ); + l1tkmu_Eta_barrel->setAxisTitle("l1tkmu_Eta_barrel",1); + + l1tkmu_Dxy_barrel = dbe->book1D("l1tkmu_Dxy_barrel","L1 Muon Dxy barrel", 100,-15.,15. ); + l1tkmu_Dxy_barrel->setAxisTitle("l1tkmu_Dxy_barrel",1); + + l1tkmu_Z_barrel = dbe->book1D("l1tkmu_Z_barrel","L1 Muon Z barrel", 100,-20.,20. ); + l1tkmu_Z_barrel->setAxisTitle("l1tkmu_Z_barrel",1); + + l1tkmu_Q_barrel = dbe->book1D("l1tkmu_Q_barrel","L1 Muon Q barrel", 4, -2.,2. ); + l1tkmu_Q_barrel->setAxisTitle("l1tkmu_Q_barrel",1); + + l1tkmu_Pdgid_barrel = dbe->book1D("l1tkmu_Pdgid_barrel","L1 Muon Pdgid barrel", 100,-50.,50. ); + l1tkmu_Pdgid_barrel->setAxisTitle("l1tkmu_Pdgid_barrel",1); + + /* + * "endcap" (|eta|>1.1) + */ + l1tkmu_Pt_endcap = dbe->book1D("l1tkmu_Pt_endcap","L1 Muon Pt endcap", 100,0.,250. ); + l1tkmu_Pt_endcap->setAxisTitle("l1tkmu_Pt_endcap",1); + + l1tkmu_Phi_endcap = dbe->book1D("l1tkmu_Phi_endcap","L1 Muon Phi endcap", 100,-3.,3. ); + l1tkmu_Phi_endcap->setAxisTitle("l1tkmu_Phi_endcap",1); + + l1tkmu_Eta_endcap = dbe->book1D("l1tkmu_Eta_endcap","L1 Muon Eta endcap", 100,-3.,3. ); + l1tkmu_Eta_endcap->setAxisTitle("l1tkmu_Eta_endcap",1); + + l1tkmu_Dxy_endcap = dbe->book1D("l1tkmu_Dxy_endcap","L1 Muon Dxy endcap", 100,-15.,15. ); + l1tkmu_Dxy_endcap->setAxisTitle("l1tkmu_Dxy_endcap",1); + + l1tkmu_Z_endcap = dbe->book1D("l1tkmu_Z_endcap","L1 Muon Z endcap", 100,-20.,20. ); + l1tkmu_Z_endcap->setAxisTitle("l1tkmu_Z_endcap",1); + + l1tkmu_Q_endcap = dbe->book1D("l1tkmu_Q_endcap","L1 Muon Q endcap", 4, -2.,2. ); + l1tkmu_Q_endcap->setAxisTitle("l1tkmu_Q_endcap",1); + + l1tkmu_Pdgid_endcap = dbe->book1D("l1tkmu_Pdgid_endcap","L1 Muon Pdgid endcap", 100,-50.,50. ); + l1tkmu_Pdgid_endcap->setAxisTitle("l1tkmu_Pdgid_endcap",1); + + } +} + +DEFINE_FWK_MODULE(L1TkMuonMonitor); diff --git a/DQM/L1TMonitor/src/SealModule.cc b/DQM/L1TMonitor/src/SealModule.cc index 093e410586e8f..d02be3825890e 100644 --- a/DQM/L1TMonitor/src/SealModule.cc +++ b/DQM/L1TMonitor/src/SealModule.cc @@ -27,7 +27,6 @@ DEFINE_FWK_MODULE(L1TRPCTF); #include DEFINE_FWK_MODULE(L1TGMT); - #include DEFINE_FWK_MODULE(L1TGCT); diff --git a/GeneratorInterface/HiGenCommon/BuildFile.xml b/GeneratorInterface/HiGenCommon/BuildFile.xml index dee5e5e0a72ab..97d107d151eda 100644 --- a/GeneratorInterface/HiGenCommon/BuildFile.xml +++ b/GeneratorInterface/HiGenCommon/BuildFile.xml @@ -2,7 +2,7 @@ - + diff --git a/HLTrigger/HLTanalyzers/BuildFile.xml b/HLTrigger/HLTanalyzers/BuildFile.xml index bfc8267fa8bc4..4316c6292fdcf 100644 --- a/HLTrigger/HLTanalyzers/BuildFile.xml +++ b/HLTrigger/HLTanalyzers/BuildFile.xml @@ -19,7 +19,7 @@ - + diff --git a/SLHCUpgradeSimulations/L1TrackTrigger/plugins/L1TkMuonFromExtendedProducer.cc b/SLHCUpgradeSimulations/L1TrackTrigger/plugins/L1TkMuonFromExtendedProducer.cc index 39a63c76185f3..2439903c264c4 100644 --- a/SLHCUpgradeSimulations/L1TrackTrigger/plugins/L1TkMuonFromExtendedProducer.cc +++ b/SLHCUpgradeSimulations/L1TrackTrigger/plugins/L1TkMuonFromExtendedProducer.cc @@ -1,17 +1,15 @@ + // -*- 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 @@ -38,13 +36,28 @@ #include #include +#include +#include +#include -using namespace l1extra ; +using namespace l1extra; // // class declaration // +float minV=1000.; +float maxV=-1000.; + +// maximum Muons number +const int MuonMax = 20; + +// maximum Tracks number +const int TrackMax = 400; + +// maximum L1 Tracks Candidates number +const int TkMuonMax = 50; + class L1TkMuonFromExtendedProducer : public edm::EDProducer { public: @@ -53,9 +66,7 @@ class L1TkMuonFromExtendedProducer : public edm::EDProducer { struct PropState { //something simple, imagine it's hardware emulation PropState() : - pt(-99), eta(-99), phi(-99), - sigmaPt(-99), sigmaEta(-99), sigmaPhi(-99), - valid(false) {} + pt(-99), eta(-99), phi(-99), sigmaPt(-99), sigmaEta(-99), sigmaPhi(-99), valid(false) {} float pt; float eta; float phi; @@ -63,18 +74,92 @@ class L1TkMuonFromExtendedProducer : public edm::EDProducer { float sigmaEta; float sigmaPhi; bool valid; + }; + struct L1Muon { + float eta; + float phi; + float pt; + float feta; + float bx; + float sigmaEta; + float sigmaPhi; + unsigned int quality; + int idx; + bool gmtCand; + bool dtCand; + bool cscCand; + bool rpcCand; + int itk; + int Num; + }; + + struct L1Track { + float pt; + int nPars; + int q; + float z; + float chi2; + float eta; + float phi; + float dxy; + int nstubs; + int idx; + bool valid; + float propPt; + float propEta; + float propPhi; + float propSigmaEta; + float propSigmaPhi; + }; + + struct L1TrackProp { + float pt; + float p; + float eta; + float phi; + float q; + float z; + }; + + struct L1Cand { + float pt; + int nPars; + float q; + float z; + float eta; + float phi; + float dxy; + int l1TrackIndex; + int l1MuonIndex; + int quality; }; + L1Muon MuCn[MuonMax]; + L1Track L1Tk[TrackMax]; + L1TrackProp L1TkPr[TrackMax]; + L1Cand L1TkCn[TkMuonMax]; + explicit L1TkMuonFromExtendedProducer(const edm::ParameterSet&); ~L1TkMuonFromExtendedProducer(); static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); - + private: + virtual void produce(edm::Event&, const edm::EventSetup&); + int loadL1Muons(edm::Event&, L1Muon*); + int loadL1Tracks(edm::Event&, L1Track*, L1TrackProp*); + + int computeTkMuCandidates(const int, L1Muon*, const int, L1Track*, L1Cand*); + void loadTkMuCandidatesToEvent(edm::Event&, std::auto_ptr&, const int, L1Cand*); + PropState propagateToGMT(const L1TkTrackType& l1tk) const; + float Pt2Bin(float Pt); + float Phi2Bin(float Phi); + float Eta2Bin(float Eta); + float Z2Bin(float Eta); //configuration (preserving structure of L1TkMuonFromExtendedProducer edm::InputTag L1MuonsInputTag_; @@ -91,35 +176,176 @@ class L1TkMuonFromExtendedProducer : public edm::EDProducer { bool correctGMTPropForTkZ_; bool use5ParameterFit_; -} ; - +}; // -// constructors and destructor +// 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"); + use5ParameterFit_ = iConfig.getParameter("use5ParameterFit"); + produces(); +} + +L1TkMuonFromExtendedProducer::~L1TkMuonFromExtendedProducer() { +} +// ------------ Muon and Tracks converter ------------ +int +L1TkMuonFromExtendedProducer::loadL1Muons(edm::Event& iEvent, L1Muon* MuCn) +{ + using namespace edm; + using namespace std; - L1MuonsInputTag_ = iConfig.getParameter("L1MuonsInputTag"); - L1TrackInputTag_ = iConfig.getParameter("L1TrackInputTag"); + edm::Handle l1musH; + iEvent.getByLabel(L1MuonsInputTag_, l1musH); + const L1MuonParticleExtendedCollection& l1mus = *l1musH.product(); + + LogDebug("L1TkMuonFromExtendedProducer") << " l1mus.size= " << l1mus.size(); - 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"); + /* L1 Muons diapasons: + * min max + * Eta -2.425 2.425 + * Phi -3.14159 3.14159 + * pt 2 140 + * sigmaEta 0.0144337 0.0288675 + * sigmaPhi 0 0.0125959 + */ - correctGMTPropForTkZ_ = iConfig.getParameter("correctGMTPropForTkZ"); + // Muons filling - use5ParameterFit_ = iConfig.getParameter("use5ParameterFit"); - produces(); + int imu = -1; + int nL1Muons = -1; + for (auto l1mu : l1mus){ + imu++; + + if( imu l1tksH; + iEvent.getByLabel(L1TrackInputTag_, l1tksH); + const L1TkTrackCollectionType& l1tks = *l1tksH.product(); + + int nL1Tracks = -1; + int il1tk = -1; + for (auto l1tk : l1tks ){ + il1tk++; + + if( il1tk0? 1.:-1.; + L1TkPr[il1tk].z = l1tk.getPOCA().z(); + + /* + PropState pstate = propagateToGMT(l1tk); + if (!pstate.valid) { + + L1Tk[il1tk].valid = true; + L1Tk[il1tk].propPt = pstate.pt; + L1Tk[il1tk].propEta = pstate.eta; + L1Tk[il1tk].propPhi = pstate.phi; + L1Tk[il1tk].propSigmaEta = pstate.sigmaEta; + L1Tk[il1tk].propSigmaPhi = pstate.sigmaPhi; + + LogDebug("L1TkMuonFromExtendedProducer") << " L1Tk[" << il1tk << "].prop " + << " Pt " << L1Tk[il1tk].propPt << " eta " << L1Tk[il1tk].propEta << " phi " << L1Tk[il1tk].propPhi + << " SigmaEta " << L1Tk[il1tk].propSigmaEta << " SigmaPhi " << L1Tk[il1tk].propSigmaPhi; + } else { + L1Tk[il1tk].valid = false; + }; + */ + } + else { + LogWarning("L1TkMuonFromExtendedProducer") << " Maximum TrackMaxn number too small " << nL1Tracks; + } + }//over l1tks + LogDebug("L1TkMuonFromExtendedProducer") << " ::loadL1Muons nL1Tracks " << nL1Tracks; + return nL1Tracks; } // ------------ method called to produce the data ------------ @@ -128,166 +354,288 @@ L1TkMuonFromExtendedProducer::produce(edm::Event& iEvent, const edm::EventSetup& { using namespace edm; using namespace std; - - + std::auto_ptr tkMuons(new L1TkMuonParticleCollection); + + int nL1Muons = L1TkMuonFromExtendedProducer::loadL1Muons(iEvent, MuCn); + LogDebug("L1TkMuonFromExtendedProducer") << " nL1Muons " << nL1Muons; - edm::Handle l1musH; - iEvent.getByLabel(L1MuonsInputTag_, l1musH); - const L1MuonParticleExtendedCollection& l1mus = *l1musH.product(); + int nL1Tracks = L1TkMuonFromExtendedProducer::loadL1Tracks(iEvent, L1Tk, L1TkPr); + LogDebug("L1TkMuonFromExtendedProducer") << " nL1Tracks " << nL1Tracks; - edm::Handle l1tksH; - iEvent.getByLabel(L1TrackInputTag_, l1tksH); - const L1TkTrackCollectionType& l1tks = *l1tksH.product(); + int nL1TkCand = computeTkMuCandidates(nL1Muons, MuCn, nL1Tracks, L1Tk, L1TkCn); + if( nL1TkCand+1 > 0 ){ + LogDebug("L1TkMuonFromExtendedProducer") << " nL1TkCand " << nL1TkCand+1; + } + + L1TkMuonFromExtendedProducer::loadTkMuCandidatesToEvent(iEvent, tkMuons, nL1TkCand, L1TkCn); - L1TkMuonParticleCollection l1tkmuCands; - l1tkmuCands.reserve(l1mus.size()*4); //can do more if really needed + iEvent.put( tkMuons ); - 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; +// ------------ Tracks candidate computation ------------ +int +L1TkMuonFromExtendedProducer::computeTkMuCandidates(const int nL1Muons, L1Muon* MuCn, const int nL1Tracks, L1Track* L1Tk, + L1Cand* L1TkCn) +{ + using namespace edm; + using namespace std; + + // PropState calculation + for ( int itk=0; itk= 0){ + + L1TkPr[itk].pt = L1TkMuonFromExtendedProducer::Pt2Bin(L1TkPr[itk].pt); + L1TkPr[itk].phi = L1TkMuonFromExtendedProducer::Phi2Bin(L1TkPr[itk].phi); + L1TkPr[itk].eta = L1TkMuonFromExtendedProducer::Eta2Bin(L1TkPr[itk].eta); + L1TkPr[itk].z = L1TkMuonFromExtendedProducer::Z2Bin(L1TkPr[itk].z); + + if (!correctGMTPropForTkZ_) L1TkPr[itk].z = 0; + if ( (L1TkPr[itk].p<3.5) || (abs(L1TkPr[itk].eta) <1.1 && L1TkPr[itk].pt < 3.5) || (abs(L1TkPr[itk].eta) > 2.5) ){ + L1Tk[itk].valid = false; + } + else { + //0th order: + L1Tk[itk].valid = true; + + float dzCorrPhi = 1.; + float deta = 0; + + float etaProp = abs(L1TkPr[itk].eta); + + if (abs(L1TkPr[itk].eta) < 1.1){ + etaProp = 1.1; + deta = L1TkPr[itk].z/550./cosh(abs(L1TkPr[itk].eta)); + } else { + float delta = L1TkPr[itk].z/850.; //roughly scales as distance to 2nd station + if (L1TkPr[itk].eta > 0) delta *=-1; + dzCorrPhi = 1. + delta; + + float zOzs = L1TkPr[itk].z/850.; + if (L1TkPr[itk].eta > 0) deta = zOzs/(1. - zOzs); + else deta = zOzs/(1.+zOzs); + deta = deta*tanh(L1TkPr[itk].eta); + } + + float resPhi = L1TkPr[itk].phi - 1.464*L1TkPr[itk].q*cosh(1.7)/cosh(etaProp)/L1TkPr[itk].pt*dzCorrPhi - M_PI/144.; + if (resPhi > M_PI) resPhi -= 2.*M_PI; + if (resPhi < -M_PI) resPhi += 2.*M_PI; + + L1Tk[itk].propPt = L1TkPr[itk].pt; //not corrected for eloss + L1Tk[itk].propEta = L1TkPr[itk].eta + deta; + L1Tk[itk].propPhi = resPhi; + L1Tk[itk].propSigmaEta = 0.100/L1TkPr[itk].pt; //multiple scattering term + L1Tk[itk].propSigmaPhi = 0.106/L1TkPr[itk].pt; //need a better estimate for these + + LogDebug("L1TkMuonFromExtendedProducer") << " L1Tk[" << itk << "].prop " + << " Pt " << L1Tk[itk].propPt << " eta " << L1Tk[itk].propEta << " phi " << L1Tk[itk].propPhi + << " SigmaEta " << L1Tk[itk].propSigmaEta << " SigmaPhi " << L1Tk[itk].propSigmaPhi; + + }// else "true" + }// if idx + }// over itk - unsigned int l1mu_quality = l1mu.quality(); + // **** Calculation for Muons **** + int nL1TkCand = -1; + int il1tkcn = -1; - 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 + for (int imu=0; imu ETAMAX_) continue; + + // can skip quality cuts at the moment, keep bx=0 req + if (MuCn[imu].bx != 0) continue; + + if (!MuCn[imu].gmtCand){ + //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()){ + if (!MuCn[imu].dtCand){ // something can be called from here } - const auto& cscCand = l1mu.cscCand(); - if (!cscCand.empty()){ + if (!MuCn[imu].cscCand){ //apply something specific here } - const auto& rpcCand = l1mu.rpcCand(); - if (!rpcCand.empty()){ + if (!MuCn[imu].rpcCand){ //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; + // **** Calculation for Tracks **** + float dr2 = -1; + float dr2prop = 999.; + float drmin = 999.; + + for ( int itk=0; itk= 0){ + + L1Tk[itk].pt = L1TkMuonFromExtendedProducer::Pt2Bin(L1Tk[itk].pt); + L1Tk[itk].phi = L1TkMuonFromExtendedProducer::Phi2Bin(L1Tk[itk].phi); + L1Tk[itk].eta = L1TkMuonFromExtendedProducer::Eta2Bin(L1Tk[itk].eta); + L1Tk[itk].z = L1TkMuonFromExtendedProducer::Z2Bin(L1Tk[itk].z); + + if ( L1Tk[itk].pt < PTMINTRA_) continue; + if (fabs(L1Tk[itk].z) > ZMAX_) continue; + if ( L1Tk[itk].chi2 > CHI2MAX_) continue; + if ( L1Tk[itk].nstubs < nStubsmin_) continue; + // Do we need it !? // if ( L1Tk[itk].valid) continue; + + MuCn[imu].pt = L1TkMuonFromExtendedProducer::Pt2Bin(MuCn[imu].pt); + MuCn[imu].phi = L1TkMuonFromExtendedProducer::Phi2Bin(MuCn[imu].phi); + MuCn[imu].eta = L1TkMuonFromExtendedProducer::Eta2Bin(MuCn[imu].eta); + + dr2 = deltaR2(MuCn[imu].eta, MuCn[imu].phi, L1Tk[itk].eta, L1Tk[itk].phi); + if (dr2 > 0.3) continue; + + dr2prop = deltaR2(MuCn[imu].eta, MuCn[imu].phi, L1Tk[itk].propEta, L1Tk[itk].propPhi); + if (dr2prop < drmin) drmin = dr2prop; + if (dr2prop < drmin) continue; + + float etaCut = 3.*sqrt(MuCn[imu].sigmaEta*MuCn[imu].sigmaEta + L1Tk[itk].propSigmaEta*L1Tk[itk].propSigmaEta); + float phiCut = 4.*sqrt(MuCn[imu].sigmaPhi*MuCn[imu].sigmaPhi + L1Tk[itk].propSigmaPhi*L1Tk[itk].propSigmaPhi); + float dEta = std::abs(L1Tk[itk].propEta - MuCn[imu].eta); + float dPhi = std::abs(deltaPhi(L1Tk[itk].propPhi, MuCn[imu].phi)); + + LogDebug("L1TkMuonFromExtendedProducer") + << " match details: prop Pt "< ZMAX_) continue; + float MinPt = 2.; // GeV + float MaxPt = 200.; // GeV - float l1tk_chi2 = l1tk.getChi2(nPars); - if (l1tk_chi2 > CHI2MAX_) continue; - - int l1tk_nstubs = l1tk.getStubRefs().size(); - if ( l1tk_nstubs < nStubsmin_) continue; + // 2**6 = 64 + float PtBinStep = (MaxPt - MinPt)/exp2(5); - float l1tk_eta = l1tk.getMomentum(nPars).eta(); - float l1tk_phi = l1tk.getMomentum(nPars).phi(); + // float tmpPt = Pt; + Pt = max( (double)MinPt, min((double)Pt, (double)MaxPt)); + Pt = (MinPt + (float)( PtBinStep * (int)(Pt/PtBinStep))); + /* + cout << " Pt " << Pt << " PtBinStep " << PtBinStep << " sPt " << tmpPt - Pt << endl; + */ + return Pt; +} - 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 L1TkMuonFromExtendedProducer::Phi2Bin(float Phi){ + using namespace std; - 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 MinPhi = -3.14159; // Rad + float MaxPhi = 3.14159; // Rad - float dEta = std::abs(matchProp.eta - l1mu.eta()); - float dPhi = std::abs(deltaPhi(matchProp.phi, l1mu.phi())); + // 2**8 = 256 + float PhiBinStep = (MaxPhi - MinPhi)/exp2(8); - LogDebug("MYDEBUG")<<"match details: prop "< l1tkPtr(l1tksH, match_idx); + Phi = max( (double)MinPhi, min((double)Phi, (double)MaxPhi)); - unsigned int nPars = 4; - if (use5ParameterFit_) nPars = 5; - auto p3 = matchTk.getMomentum(nPars); - float p4e = sqrt(0.105658369*0.105658369 + p3.mag2() ); + //float tmpPhi = Phi; + Phi = (float)( PhiBinStep * (int)(Phi/PhiBinStep)); + if (Phi < 0.){ + Phi = -(float)( PhiBinStep * (int)(abs(Phi)/PhiBinStep)); + } - math::XYZTLorentzVector l1tkp4(p3.x(), p3.y(), p3.z(), p4e); + /* + cout << " tmpPhi " << tmpPhi << " Phi " << Phi << " dPhi " << (tmpPhi - Phi) << endl; + */ + return Phi; +} - auto tkv3=matchTk.getPOCA(nPars); - math::XYZPoint v3(tkv3.x(), tkv3.y(), tkv3.z()); - float trkisol = -999; - int l1tk_q = matchTk.getRInv(nPars)>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); +float L1TkMuonFromExtendedProducer::Eta2Bin(float Eta){ + using namespace std; - // EP: add the zvtx information - l1tkmu.setTrkzVtx( (float)tkv3.z() ); + float MinEta = -2.65; + float MaxEta = 2.65; - tkMuons->push_back(l1tkmu); - } - } + // 2**8 = 256 + float EtaBinStep = (MaxEta - MinEta)/exp2(8); - }//over l1mus - - - iEvent.put( tkMuons ); - + Eta = max( (double)MinEta, min((double)Eta, (double)MaxEta)); + + //float tmpEta = Eta; + Eta = (float)( EtaBinStep * (int)(Eta/EtaBinStep)); + if (Eta < 0.){ + Eta = -(float)( EtaBinStep * (int)(abs(Eta)/EtaBinStep)); + } + + /* + cout << " tmpEta " << tmpEta << " Eta " << Eta << " dEta " << (tmpEta - Eta) << endl; + */ + return Eta; } +float L1TkMuonFromExtendedProducer::Z2Bin(float Z){ + using namespace std; -// ------------ 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); + float MinZ = -17; + float MaxZ = 17; + + // 2**4 = 16 + float ZBinStep = (MaxZ - MinZ)/exp2(4); + + Z = max( (double)MinZ, min((double)Z, (double)MaxZ)); + + //float tmpZ = Z; + Z = (float)( ZBinStep * (int)(Z/ZBinStep)); + if (Z < 0.){ + Z = -(float)( ZBinStep * (int)(abs(Z)/ZBinStep)); + } + + /* + cout << " tmpZ " << tmpZ << " Z " << Z << " dZ " << (tmpZ - Z) << endl; + */ + return Z; } -L1TkMuonFromExtendedProducer::PropState L1TkMuonFromExtendedProducer::propagateToGMT(const L1TkMuonFromExtendedProducer::L1TkTrackType& tk) const { +L1TkMuonFromExtendedProducer::PropState L1TkMuonFromExtendedProducer::propagateToGMT(const L1TkMuonFromExtendedProducer::L1TkTrackType& tk) +const { + auto p3 = tk.getMomentum(); float tk_pt = p3.perp(); float tk_p = p3.mag(); @@ -296,6 +644,7 @@ L1TkMuonFromExtendedProducer::PropState L1TkMuonFromExtendedProducer::propagateT 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; @@ -330,14 +679,78 @@ L1TkMuonFromExtendedProducer::PropState L1TkMuonFromExtendedProducer::propagateT 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); + +// ------------ Tracks candidate ------------ +void +L1TkMuonFromExtendedProducer::loadTkMuCandidatesToEvent(edm::Event& iEvent, std::auto_ptr& tkMuons, +const int nL1TkCand, L1Cand* L1TkCn) +{ + using namespace edm; + using namespace std; + + edm::Handle l1tksH; + iEvent.getByLabel(L1TrackInputTag_, l1tksH); + const L1TkTrackCollectionType& l1tks = *l1tksH.product(); + + edm::Handle l1musH; + iEvent.getByLabel(L1MuonsInputTag_, l1musH); + const L1MuonParticleExtendedCollection& l1mus = *l1musH.product(); + + L1TkMuonParticleCollection l1tkmuCands; + l1tkmuCands.reserve(l1mus.size()*4); //can do more if really needed + LogDebug("L1TkMuonFromExtendedProducer") << " l1tks.size= " << l1tks.size(); + + for ( int il1tkcn=0; il1tkcn<=nL1TkCand; il1tkcn++ ){ + + Ptr< L1TkTrackType > l1tkPtr(l1tksH, L1TkCn[il1tkcn].l1TrackIndex); + const L1TkTrackType& matchTk = l1tks[L1TkCn[il1tkcn].l1TrackIndex]; + + auto p3 = matchTk.getMomentum(L1TkCn[il1tkcn].nPars); + float p4e = sqrt(0.105658369*0.105658369 + p3.mag2() ); + + math::XYZTLorentzVector l1tkp4(p3.x(), p3.y(), p3.z(), p4e); + + auto tkv3=matchTk.getPOCA(L1TkCn[il1tkcn].nPars); + math::XYZPoint v3(tkv3.x(), tkv3.y(), tkv3.z()); + + L1TkCn[il1tkcn].q = matchTk.getRInv(L1TkCn[il1tkcn].nPars)>0? 1: -1; + + L1TkMuonParticle l1tkmu(reco::LeafCandidate(L1TkCn[il1tkcn].q, l1tkp4, v3, -13*L1TkCn[il1tkcn].q )); + + l1tkmu.setTrkPtr(l1tkPtr); + + L1MuonParticleExtendedRef l1muRef(l1musH, L1TkCn[il1tkcn].l1MuonIndex); + l1tkmu.setMuExtendedRef(l1muRef); + + l1tkmu.setQuality(L1TkCn[il1tkcn].quality); + + float trkisol = -999; + l1tkmu.setTrkIsol(trkisol); + + // EP: add the zvtx information + l1tkmu.setTrkzVtx( (float)tkv3.z() ); + + tkMuons->push_back(l1tkmu); + + }//over il1tkcn +} +// ------------ 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); +} +//define this as a plug-in +DEFINE_FWK_MODULE(L1TkMuonFromExtendedProducer); diff --git a/SLHCUpgradeSimulations/L1TrackTrigger/test/BuildFile.xml b/SLHCUpgradeSimulations/L1TrackTrigger/test/BuildFile.xml index 70e1f0d1a5bd6..a5bf9e2cbf3ae 100644 --- a/SLHCUpgradeSimulations/L1TrackTrigger/test/BuildFile.xml +++ b/SLHCUpgradeSimulations/L1TrackTrigger/test/BuildFile.xml @@ -12,7 +12,9 @@ +