Skims/src/RPCRecHitFilter.cc

0001 // Orso Iorio, INFN Napoli
0002 //
0003 //This module is a filter based on the number of RPC Hits
0004 //
0005
0006 #include "DPGAnalysis/Skims/src/RPCRecHitFilter.h"
0007
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "DataFormats/Common/interface/Handle.h"
0013 #include <DataFormats/RPCRecHit/interface/RPCRecHit.h>
0014 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
0015
0016 #include <Geometry/DTGeometry/interface/DTGeometry.h>
0017 #include "Geometry/RPCGeometry/interface/RPCGeomServ.h"
0018 #include "MagneticField/Engine/interface/MagneticField.h"
0019 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0020 #include <Geometry/Records/interface/MuonGeometryRecord.h>
0021 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0022 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0023 #include "DataFormats/DetId/interface/DetId.h"
0024 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0025
0026 #include <FWCore/Utilities/interface/Exception.h>
0027 #include "DataFormats/Provenance/interface/Timestamp.h"
0028
0029 /////Trigger
0030
0031 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
0032 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuRegionalCand.h"
0033 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutCollection.h"
0034
0035 #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTCand.h"
0036 #include "CondFormats/L1TObjects/interface/L1GtParameters.h"
0037 #include "CondFormats/DataRecord/interface/L1GtParametersRcd.h"
0038
0039 using namespace reco;
0040
0041 typedef std::vector<L1MuRegionalCand> L1MuRegionalCandCollection;
0042
0043 #include <sys/time.h>
0044 #include <algorithm>
0045 #include <memory>
0046 #include <cmath>
0047 #include <cmath>
0048 #include "TFile.h"
0049 #include "TMath.h"
0050 #include "TTree.h"
0051
0052 #include "TDirectory.h"
0053 #include "TFile.h"
0054 #include "TTree.h"
0055 #include <cstdlib>
0056 #include <cstdio>
0057 #include <cstdlib>
0058 #include <string>
0059 #include <memory>
0060 #include <cmath>
0061
0062 using namespace edm;
0063 using namespace reco;
0064 using namespace std;
0065
0066 RPCRecHitFilter::RPCRecHitFilter(const edm::ParameterSet& iConfig)
0067     : rpcGeomToken_(esConsumes()), trackingGeoToken_(esConsumes()) {
0068   LogTrace("RPCEffTrackExtrapolation") << "inside constructor" << std::endl;
0069
0070   RPCDataLabel_ = iConfig.getUntrackedParameter<std::string>("rpcRecHitLabel");
0071   rpcRecHitToken_ = consumes<RPCRecHitCollection>(edm::InputTag(RPCDataLabel_));
0072
0073   centralBX_ = iConfig.getUntrackedParameter<int>("CentralBunchCrossing", 0);
0074   BXWindow_ = iConfig.getUntrackedParameter<int>("BunchCrossingWindow", 9999);
0075
0076   minHits_ = iConfig.getUntrackedParameter<int>("minimumNumberOfHits", 0);
0077   hitsInStations_ = iConfig.getUntrackedParameter<int>("HitsInStation", 0);
0078
0079   Debug_ = iConfig.getUntrackedParameter<bool>("Debug", false);
0080   Verbose_ = iConfig.getUntrackedParameter<bool>("Verbose", false);
0081
0082   Barrel_ = iConfig.getUntrackedParameter<bool>("UseBarrel", true);
0083   EndcapPositive_ = iConfig.getUntrackedParameter<bool>("UseEndcapPositive", true);
0084   EndcapNegative_ = iConfig.getUntrackedParameter<bool>("UseEndcapNegative", true);
0085
0086   cosmicsVeto_ = iConfig.getUntrackedParameter<bool>("CosmicsVeto", false);
0087 }
0088
0089 bool RPCRecHitFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0090   edm::ESHandle<RPCGeometry> rpcGeo = iSetup.getHandle(rpcGeomToken_);
0091   edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry = iSetup.getHandle(trackingGeoToken_);
0092
0093   edm::Handle<RPCRecHitCollection> rpcHits = iEvent.getHandle(rpcRecHitToken_);
0094
0095   std::map<int, int> numberOfRecHitsBarrel;
0096   std::map<int, int> numberOfDigisBarrel;
0097   std::map<int, int> numberOfRecHitsEndcap;
0098   std::map<int, int> numberOfDigisEndcap;
0099
0100   std::map<pair<int, int>, std::vector<RPCDetId> > numberOfRecHitsSameWheelSameSector;
0101   std::map<pair<int, int>, std::vector<RPCDetId> > numberOfDigisSameWheelSameSector;
0102   std::map<pair<int, int>, std::vector<RPCDetId> > numberOfHitsSameDiskSectorPositive;
0103   std::map<pair<int, int>, std::vector<RPCDetId> > numberOfHitsSameDiskSectorNegative;
0104
0105   const std::vector<const RPCRoll*> rls = rpcGeo->rolls();
0106
0107   bool condition = true;
0108
0109   int nBarrel = 0;
0110   int nEndcap = 0;
0111
0112   /////BEGIN LOOP ON THE ROLLS
0113   for (int i = 0; i < (int)rls.size(); ++i) {
0114     RPCDetId did = rls[i]->id();
0115
0116     /// LOOP OVER THE RECHITS
0117     RPCRecHitCollection::range rpcRecHitRange = rpcHits->get(did);
0118     RPCRecHitCollection::const_iterator RecHitsIt;
0119
0120     for (RecHitsIt = rpcRecHitRange.first; RecHitsIt != rpcRecHitRange.second; ++RecHitsIt) {
0121       //std::cout<< " roll is "<< did << " bx recHit is " << RecHitsIt->BunchX()<< " event number " <<eventNumber  <<std::endl;
0122
0123       if (did.region() == 0) {
0124         if (cosmicsVeto_) {
0125           for (int u = -1; u <= 1; ++u) {
0126             for (int t = -1; t <= 1; ++t) {
0127               numberOfRecHitsSameWheelSameSector[pair<int, int>(did.ring() + u, did.sector() + t)].push_back(did);
0128             }
0129           }
0130         }
0131
0132         else {
0133           if (did.station() == 1) {
0134             for (int u = -1; u <= 1; ++u) {
0135               for (int t = -1; t <= 1; ++t) {
0136                 numberOfRecHitsSameWheelSameSector[pair<int, int>(did.ring() + u, did.sector() + t)].push_back(did);
0137               }
0138             }
0139           }
0140         }
0141
0142         ++numberOfRecHitsBarrel[did.ring()];
0143         ++nBarrel;
0144
0145       }
0146
0147       else {
0148         if (did.region() == -1) {
0149           for (int t = -1; t <= 1; ++t) {
0150             numberOfHitsSameDiskSectorNegative[pair<int, int>(did.ring(), did.sector() + t)].push_back(did);
0151           }
0152         }
0153
0154         if (did.region() == 1) {
0155           for (int t = -1; t <= 1; ++t) {
0156             numberOfHitsSameDiskSectorPositive[pair<int, int>(did.ring(), did.sector() + t)].push_back(did);
0157           }
0158         }
0159         ++numberOfRecHitsEndcap[did.station()];
0160
0161         ++nEndcap;
0162       }
0163
0164       break;
0165     }
0166   }
0167   /////END OF THE LOOP ON THE ROLLS
0168
0169   ////CONDITION IS THAT THERE ARE TWO HITS IN RB1in and RB1out OF TWO ADJACENT SECTORS AND WHEELS
0170   ////OR THAT IN TWO DISKS SAME OF THE ENDCAP THERE ARE TWO HITS IN ADJACENT SECTORS IN THE SAME RING
0171   bool cond1 = false;
0172   bool cond2 = false;
0173   bool cond3 = false;
0174
0175   std::map<int, bool> vectorBarrelCands;
0176   std::map<int, bool> vectorEndcapCandsPositive;
0177   std::map<int, bool> vectorEndcapCandsNegative;
0178
0179   // barrel
0180   for (std::map<pair<int, int>, std::vector<RPCDetId> >::const_iterator iter =
0181            numberOfRecHitsSameWheelSameSector.begin();
0182        iter != numberOfRecHitsSameWheelSameSector.end();
0183        ++iter) {
0184     vectorBarrelCands[1] = false;
0185     vectorBarrelCands[2] = false;
0186
0187     if (iter->second.size() > 1) {
0188       for (size_t i = 0; i < iter->second.size(); ++i) {
0189         if (iter->second[i].layer() == 1 && iter->second[i].station() == 1)
0190           vectorBarrelCands[0] = true;
0191         if (iter->second[i].layer() == 2 && iter->second[i].station() == 1)
0192           vectorBarrelCands[1] = true;
0193         if (cosmicsVeto_)
0194           if (iter->second[i].station() > 2) {
0195             vectorBarrelCands[1] = false;
0196             vectorBarrelCands[2] = false;
0197             break;
0198           }
0199       }
0200     }
0201
0202     if ((vectorBarrelCands[0] && vectorBarrelCands[1])) {
0203       cond1 = true;
0204       break;
0205     }
0206   }
0207
0208   // endcap positive
0209   for (std::map<pair<int, int>, std::vector<RPCDetId> >::const_iterator iter =
0210            numberOfHitsSameDiskSectorPositive.begin();
0211        iter != numberOfHitsSameDiskSectorPositive.end();
0212        ++iter) {
0213     vectorEndcapCandsPositive[1] = false;
0214     vectorEndcapCandsPositive[2] = false;
0215     vectorEndcapCandsPositive[3] = false;
0216
0217     if (iter->second.size() > 1) {
0218       for (size_t i = 0; i < iter->second.size(); ++i) {
0219         if (iter->second[i].station() == 1)
0220           vectorEndcapCandsPositive[1] = true;
0221         if (iter->second[i].station() == 2)
0222           vectorEndcapCandsPositive[2] = true;
0223         if (iter->second[i].station() == 3)
0224           vectorEndcapCandsPositive[3] = true;
0225       }
0226     }
0227
0228     if (((vectorEndcapCandsPositive[1] && vectorEndcapCandsPositive[2]) ||
0229          (vectorEndcapCandsPositive[1] && vectorEndcapCandsPositive[3]) ||
0230          (vectorEndcapCandsPositive[2] && vectorEndcapCandsPositive[3]))) {
0231       cond2 = true;
0232       break;
0233     }
0234   }
0235
0236   // endcap negative
0237   for (std::map<pair<int, int>, std::vector<RPCDetId> >::const_iterator iter =
0238            numberOfHitsSameDiskSectorNegative.begin();
0239        iter != numberOfHitsSameDiskSectorNegative.end();
0240        ++iter) {
0241     vectorEndcapCandsNegative[1] = false;
0242     vectorEndcapCandsNegative[2] = false;
0243     vectorEndcapCandsNegative[3] = false;
0244
0245     if (iter->second.size() > 1) {
0246       for (size_t i = 0; i < iter->second.size(); ++i) {
0247         if (iter->second[i].station() == 1)
0248           vectorEndcapCandsNegative[1] = true;
0249         if (iter->second[i].station() == 2)
0250           vectorEndcapCandsNegative[2] = true;
0251         if (iter->second[i].station() == 3)
0252           vectorEndcapCandsNegative[3] = true;
0253       }
0254     }
0255
0256     if (((vectorEndcapCandsNegative[1] && vectorEndcapCandsNegative[2]) ||
0257          (vectorEndcapCandsNegative[1] && vectorEndcapCandsNegative[3]) ||
0258          (vectorEndcapCandsNegative[2] && vectorEndcapCandsNegative[3]))) {
0259       cond3 = true;
0260       break;
0261     }
0262   }
0263
0264   condition = condition && (nBarrel + nEndcap >= minHits_);
0265
0266   cond1 = Barrel_ && cond1;
0267   cond2 = EndcapPositive_ && cond2;
0268   cond3 = EndcapNegative_ && cond3;
0269
0270   bool condition2 = (cond1 || cond2 || cond3);
0271   if (Barrel_ || EndcapPositive_ || EndcapNegative_)
0272     condition = condition && condition2;
0273
0274   return condition;
0275 }
0276
0277 DEFINE_FWK_MODULE(RPCRecHitFilter);