MuonRPCGeometry/plugins/RPCPhiEff.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RPCPhiEff
0004 // Class:      RPCPhiEff
0005 //
0006 /**\class RPCPhiEff RPCPhiEff.cc MyLib/RPCPhiEff/src/RPCPhiEff.cc
0007
0008 Description: <one line class summary>
0009
0010 Implementation:
0011 <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Tomasz Fruboes
0015 //         Created:  Wed Mar  7 08:31:57 CET 2007
0016 //
0017 //
0018
0019 #include "Validation/MuonRPCGeometry/plugins/RPCPhiEff.h"
0020 // system include files
0021 #include <memory>
0022
0023 // user include files
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029
0030 #include <L1Trigger/RPCTrigger/interface/RPCLogCone.h>
0031 #include <L1Trigger/RPCTrigger/interface/RPCConst.h>
0032
0033 #include "DataFormats/Math/interface/LorentzVector.h"
0034
0035 #include <iostream>
0036 #include <set>
0037 #include <fstream>
0038 #include <sstream>
0039
0040 #include "L1Trigger/RPCTrigger/interface/RPCConst.h"
0041
0042 #include "DataFormats/Math/interface/deltaR.h"
0043
0044 //
0045 // constants, enums and typedefs
0046 //
0047
0048 //
0049 // static data member definitions
0050 //
0051
0052 //
0053 // constructors and destructor
0054 //
0055 RPCPhiEff::RPCPhiEff(const edm::ParameterSet& ps)
0056     : m_rpcbToken(consumes<std::vector<L1MuRegionalCand> >(ps.getParameter<edm::InputTag>("rpcb"))),
0057       m_rpcfToken(consumes<std::vector<L1MuRegionalCand> >(ps.getParameter<edm::InputTag>("rpcf"))),
0058       m_g4Token(consumes<edm::SimTrackContainer>(ps.getParameter<edm::InputTag>("g4"))),
0059       m_rpcdigiToken(consumes<RPCDigiCollection>(ps.getParameter<edm::InputTag>("rpcdigi"))) {
0060   //m_outfileC.open("phieffC.txt");
0061   m_outfileR.open("muons.txt");
0062 }
0063
0064 RPCPhiEff::~RPCPhiEff() {
0065   //m_outfileC.close();
0066   m_outfileR.close();
0067 }
0068
0069 //
0070 // member functions
0071 //
0072
0073 // ------------ method called to for each event  ------------
0074 void RPCPhiEff::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0075   edm::Handle<std::vector<L1MuRegionalCand> > rpcBarrel;
0076   edm::Handle<std::vector<L1MuRegionalCand> > rpcForward;
0077
0078   iEvent.getByToken(m_rpcbToken, rpcBarrel);
0079   iEvent.getByToken(m_rpcfToken, rpcForward);
0080
0081   std::vector<edm::Handle<std::vector<L1MuRegionalCand> > > handleVec;
0082   handleVec.push_back(rpcBarrel);
0083   handleVec.push_back(rpcForward);
0084
0085   // Fetch MCdata
0086   edm::Handle<edm::SimTrackContainer> simTracks;
0087   iEvent.getByToken(m_g4Token, simTracks);
0088
0089   float etaGen = 0, phiGen = 0, ptGen = 0;
0090
0091   int noOfRecMuons = 0;
0092   int noOfMatchedRecMuons = 0;
0093
0094   bool firstRunForMuonMatchingCnt = true;
0095
0096   // ask MC muons to be separated
0097   double deltarMin = -1;
0098   edm::SimTrackContainer::const_iterator simTrk = simTracks->begin();
0099   for (; simTrk != simTracks->end(); ++simTrk) {
0100     if (simTrk->type() != -13 && simTrk->type() != 13)
0101       continue;
0102
0103     edm::SimTrackContainer::const_iterator simTrk2 = simTrk;
0104     ++simTrk2;
0105     for (; simTrk2 != simTracks->end(); ++simTrk2) {
0106       if (simTrk2->type() != -13 && simTrk2->type() != 13)
0107         continue;
0108       double drCand = reco::deltaR(simTrk2->momentum(), simTrk->momentum());
0109       if (drCand < deltarMin || deltarMin < 0)
0110         deltarMin = drCand;
0111     }
0112   }
0113
0114   //std::cout << deltarMin << std::endl;
0115   if (deltarMin < 0.7 && deltarMin > 0)
0116     return;
0117
0118   simTrk = simTracks->begin();
0119   for (; simTrk != simTracks->end(); simTrk++) {
0120     int type = simTrk->type();
0121     if (type == 13 || type == -13) {
0122       // Get the data
0123       const math::XYZTLorentzVectorD momentum = simTrk->momentum();
0124       etaGen = momentum.eta();
0125       ptGen = momentum.Pt();
0126       phiGen = momentum.phi();
0127
0128       bool matched = false;
0129       int ptCodeRec = 0;
0130       int towerRec = 0;
0131       int phiRec = 0;
0132       //int muonsFound=0;
0133       int qual = 0;
0134       int ghost = 0;  // number of ghost for montecarlo muon
0135
0136       // Iter rpc muon cands, perform delta R matching
0137       // todo perform matching also using eta...
0138       for (std::vector<edm::Handle<std::vector<L1MuRegionalCand> > >::iterator it = handleVec.begin();
0139            it != handleVec.end();
0140            ++it) {
0141         std::vector<L1MuRegionalCand>::const_iterator itRPC;
0142         for (itRPC = (*it)->begin(); itRPC != (*it)->end(); itRPC++) {
0143           int ptCode = itRPC->pt_packed();
0144           if (ptCode != 0) {
0145             if (firstRunForMuonMatchingCnt)
0146               ++noOfRecMuons;
0147             ptCodeRec = ptCode;
0148             phiRec = itRPC->phi_packed();
0149             qual = itRPC->quality();
0150             towerRec = itRPC->eta_packed();
0151             if (towerRec > 16) {
0152               towerRec = -((~towerRec & 63) + 1);
0153             }
0154             // match
0155             float pi = 3.14159265;
0156             //float phiPhys = 2*pi*float(phiRec)/144-pi;
0157             float phiScaled = phiGen;
0158             if (phiScaled < 0)
0159               phiScaled += 2 * pi;
0160             int phiHwGen = (phiScaled) / 2 / pi * 144;
0161
0162             //std::cout << "phi " << phiGen << " " << phiHwGen << " " << phiRec
0163             //          << " eta " << etaGen << " " << m_const.towerNumFromEta(etaGen) << " "<< towerRec << std::endl;
0164             if ((std::abs(phiHwGen - phiRec) < 10) && (std::abs(m_const.towerNumFromEta(etaGen) - towerRec) < 1)) {
0165               if (matched) {  // we have matched m.c. earlier, this is ghost
0166                 ++ghost;
0167               }
0168               matched = true;
0169               ++noOfMatchedRecMuons;
0170               m_outfileR << etaGen << " " << phiGen << " " << ptGen << " " << towerRec << " " << phiRec << " "
0171                          << ptCodeRec << " " << qual << " " << ghost << std::endl;
0172             }
0173           }  // (ptCode != 0)
0174         }    // muon cands iter ends
0175       }      // barrell/fwd iter ends
0176       firstRunForMuonMatchingCnt = false;
0177       if (!matched) {
0178         m_outfileR << etaGen << " " << phiGen << " " << ptGen << " " << 0 << " " << 0 << " " << 0 << " " << 0 << " "
0179                    << 0 << std::endl;
0180       }
0181     }
0182   }
0183
0184   edm::EventID id = iEvent.id();
0185   edm::EventNumber_t evNum = id.event();
0186   edm::RunNumber_t rnNum = id.run();
0187
0188   if (noOfMatchedRecMuons != noOfRecMuons) {
0189     edm::LogInfo("RPCEffWarn") << " MuonCands " << noOfRecMuons << " matched " << noOfMatchedRecMuons << " in run "
0190                                << rnNum << " event " << evNum;
0191   }
0192
0193   /*
0194     m_outfileC << etaGen << " " << phiGen << " " << ptGen << " "
0195     << phiRec << " " << towerRec << " " << muonsFound  << " "
0196     <<  fromCones(iEvent) <<std::endl;*/
0197   /*m_outfileR << etaGen << " " << phiGen << " " << ptGen << " "
0198     << phiRec << " " << towerRec << " " << muonsFound  << " "
0199     <<  fromRaw(iEvent) <<std::endl;*/
0200
0201   /*
0202     m_outfileR << etaGen << " " << phiGen << " " << ptGen << " "
0203     << phiRec << " " << towerRec << " " << muonsFound  << " "
0204     << std::endl;
0205
0206   */
0207 }
0208
0209 std::string RPCPhiEff::fromCones(const edm::Event& iEvent) { return ""; }
0210
0211 // ------------ Check hw planes fired using rpc digis
0212 std::string RPCPhiEff::fromRaw(const edm::Event& iEvent) {
0213   std::stringstream ss;
0214
0215   // Digi data.
0216   edm::Handle<RPCDigiCollection> rpcDigis;
0217   iEvent.getByToken(m_rpcdigiToken, rpcDigis);
0218
0219   std::set<int> hwPlanes;
0220
0221   RPCDigiCollection::DigiRangeIterator detUnitIt;
0222   for (detUnitIt = rpcDigis->begin(); detUnitIt != rpcDigis->end(); ++detUnitIt) {
0223     const RPCDigiCollection::Range& range = (*detUnitIt).second;
0224     bool hasBX0 = false;
0225     for (RPCDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
0226       if (digiIt->bx() == 0) {
0227         hasBX0 = true;
0228         break;
0229       }
0230     }
0231
0232     if (!hasBX0)
0233       continue;
0234
0235     const RPCDetId& id = (*detUnitIt).first;
0236     int station = id.station();
0237     int layer = id.layer();
0238     //int region = id.region();
0239
0240     if (station == 3)
0241       hwPlanes.insert(5);
0242
0243     else if (station == 4)
0244       hwPlanes.insert(6);
0245
0246     else if (station == 1 && layer == 1)
0247       hwPlanes.insert(1);
0248
0249     else if (station == 1 && layer == 2)
0250       hwPlanes.insert(2);
0251
0252     else if (station == 2 && layer == 1)
0253       hwPlanes.insert(3);
0254
0255     else if (station == 2 && layer == 2)
0256       hwPlanes.insert(4);
0257
0258     else
0259       std::cout << "??????????????" << std::endl;
0260   }
0261
0262   for (std::set<int>::iterator it = hwPlanes.begin(); it != hwPlanes.end(); ++it) {
0263     ss << " " << *it;
0264   }
0265
0266   return ss.str();
0267 }
0268
0269 // ------------ method called once each job just before starting event loop  ------------
0270 void RPCPhiEff::beginJob() {}
0271
0272 // ------------ method called once each job just after ending the event loop  ------------
0273 void RPCPhiEff::endJob() {}