Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-28 05:58:10

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 noOfMuons = 0;
0092   int noOfRecMuons = 0;
0093   int noOfMatchedRecMuons = 0;
0094 
0095   bool firstRunForMuonMatchingCnt = true;
0096 
0097   // ask MC muons to be separated
0098   double deltarMin = -1;
0099   edm::SimTrackContainer::const_iterator simTrk = simTracks->begin();
0100   for (; simTrk != simTracks->end(); ++simTrk) {
0101     if (simTrk->type() != -13 && simTrk->type() != 13)
0102       continue;
0103 
0104     edm::SimTrackContainer::const_iterator simTrk2 = simTrk;
0105     ++simTrk2;
0106     for (; simTrk2 != simTracks->end(); ++simTrk2) {
0107       if (simTrk2->type() != -13 && simTrk2->type() != 13)
0108         continue;
0109       double drCand = reco::deltaR(simTrk2->momentum(), simTrk->momentum());
0110       if (drCand < deltarMin || deltarMin < 0)
0111         deltarMin = drCand;
0112     }
0113   }
0114 
0115   //std::cout << deltarMin << std::endl;
0116   if (deltarMin < 0.7 && deltarMin > 0)
0117     return;
0118 
0119   simTrk = simTracks->begin();
0120   for (; simTrk != simTracks->end(); simTrk++) {
0121     int type = simTrk->type();
0122     if (type == 13 || type == -13) {
0123       // Get the data
0124       const math::XYZTLorentzVectorD momentum = simTrk->momentum();
0125       etaGen = momentum.eta();
0126       ptGen = momentum.Pt();
0127       phiGen = momentum.phi();
0128       noOfMuons++;
0129 
0130       bool matched = false;
0131       int ptCodeRec = 0;
0132       int towerRec = 0;
0133       int phiRec = 0;
0134       //int muonsFound=0;
0135       int qual = 0;
0136       int ghost = 0;  // number of ghost for montecarlo muon
0137 
0138       // Iter rpc muon cands, perform delta R matching
0139       // todo perform matching also using eta...
0140       for (std::vector<edm::Handle<std::vector<L1MuRegionalCand> > >::iterator it = handleVec.begin();
0141            it != handleVec.end();
0142            ++it) {
0143         std::vector<L1MuRegionalCand>::const_iterator itRPC;
0144         for (itRPC = (*it)->begin(); itRPC != (*it)->end(); itRPC++) {
0145           int ptCode = itRPC->pt_packed();
0146           if (ptCode != 0) {
0147             if (firstRunForMuonMatchingCnt)
0148               ++noOfRecMuons;
0149             ptCodeRec = ptCode;
0150             phiRec = itRPC->phi_packed();
0151             qual = itRPC->quality();
0152             towerRec = itRPC->eta_packed();
0153             if (towerRec > 16) {
0154               towerRec = -((~towerRec & 63) + 1);
0155             }
0156             // match
0157             float pi = 3.14159265;
0158             //float phiPhys = 2*pi*float(phiRec)/144-pi;
0159             float phiScaled = phiGen;
0160             if (phiScaled < 0)
0161               phiScaled += 2 * pi;
0162             int phiHwGen = (phiScaled) / 2 / pi * 144;
0163 
0164             //std::cout << "phi " << phiGen << " " << phiHwGen << " " << phiRec
0165             //          << " eta " << etaGen << " " << m_const.towerNumFromEta(etaGen) << " "<< towerRec << std::endl;
0166             if ((std::abs(phiHwGen - phiRec) < 10) && (std::abs(m_const.towerNumFromEta(etaGen) - towerRec) < 1)) {
0167               if (matched) {  // we have matched m.c. earlier, this is ghost
0168                 ++ghost;
0169               }
0170               matched = true;
0171               ++noOfMatchedRecMuons;
0172               m_outfileR << etaGen << " " << phiGen << " " << ptGen << " " << towerRec << " " << phiRec << " "
0173                          << ptCodeRec << " " << qual << " " << ghost << std::endl;
0174             }
0175           }  // (ptCode != 0)
0176         }    // muon cands iter ends
0177       }      // barrell/fwd iter ends
0178       firstRunForMuonMatchingCnt = false;
0179       if (!matched) {
0180         m_outfileR << etaGen << " " << phiGen << " " << ptGen << " " << 0 << " " << 0 << " " << 0 << " " << 0 << " "
0181                    << 0 << std::endl;
0182       }
0183     }
0184   }
0185 
0186   edm::EventID id = iEvent.id();
0187   edm::EventNumber_t evNum = id.event();
0188   edm::RunNumber_t rnNum = id.run();
0189 
0190   if (noOfMatchedRecMuons != noOfRecMuons) {
0191     edm::LogInfo("RPCEffWarn") << " MuonCands " << noOfRecMuons << " matched " << noOfMatchedRecMuons << " in run "
0192                                << rnNum << " event " << evNum;
0193   }
0194 
0195   /*
0196     m_outfileC << etaGen << " " << phiGen << " " << ptGen << " "
0197     << phiRec << " " << towerRec << " " << muonsFound  << " "
0198     <<  fromCones(iEvent) <<std::endl;*/
0199   /*m_outfileR << etaGen << " " << phiGen << " " << ptGen << " "
0200     << phiRec << " " << towerRec << " " << muonsFound  << " "
0201     <<  fromRaw(iEvent) <<std::endl;*/
0202 
0203   /*
0204     m_outfileR << etaGen << " " << phiGen << " " << ptGen << " "
0205     << phiRec << " " << towerRec << " " << muonsFound  << " "
0206     << std::endl;
0207 
0208   */
0209 }
0210 
0211 std::string RPCPhiEff::fromCones(const edm::Event& iEvent) { return ""; }
0212 
0213 // ------------ Check hw planes fired using rpc digis
0214 std::string RPCPhiEff::fromRaw(const edm::Event& iEvent) {
0215   std::stringstream ss;
0216 
0217   // Digi data.
0218   edm::Handle<RPCDigiCollection> rpcDigis;
0219   iEvent.getByToken(m_rpcdigiToken, rpcDigis);
0220 
0221   std::set<int> hwPlanes;
0222 
0223   RPCDigiCollection::DigiRangeIterator detUnitIt;
0224   for (detUnitIt = rpcDigis->begin(); detUnitIt != rpcDigis->end(); ++detUnitIt) {
0225     const RPCDigiCollection::Range& range = (*detUnitIt).second;
0226     bool hasBX0 = false;
0227     for (RPCDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
0228       if (digiIt->bx() == 0) {
0229         hasBX0 = true;
0230         break;
0231       }
0232     }
0233 
0234     if (!hasBX0)
0235       continue;
0236 
0237     const RPCDetId& id = (*detUnitIt).first;
0238     int station = id.station();
0239     int layer = id.layer();
0240     //int region = id.region();
0241 
0242     if (station == 3)
0243       hwPlanes.insert(5);
0244 
0245     else if (station == 4)
0246       hwPlanes.insert(6);
0247 
0248     else if (station == 1 && layer == 1)
0249       hwPlanes.insert(1);
0250 
0251     else if (station == 1 && layer == 2)
0252       hwPlanes.insert(2);
0253 
0254     else if (station == 2 && layer == 1)
0255       hwPlanes.insert(3);
0256 
0257     else if (station == 2 && layer == 2)
0258       hwPlanes.insert(4);
0259 
0260     else
0261       std::cout << "??????????????" << std::endl;
0262   }
0263 
0264   for (std::set<int>::iterator it = hwPlanes.begin(); it != hwPlanes.end(); ++it) {
0265     ss << " " << *it;
0266   }
0267 
0268   return ss.str();
0269 }
0270 
0271 // ------------ method called once each job just before starting event loop  ------------
0272 void RPCPhiEff::beginJob() {}
0273 
0274 // ------------ method called once each job just after ending the event loop  ------------
0275 void RPCPhiEff::endJob() {}