Reference Manual

Go to the documentation of this file.
00001 #include "Validation/RPCRecHits/interface/RPCRecHitValid.h"
00002 #include "FWCore/Framework/interface/MakerMacros.h"
00003 
00004 #include "FWCore/Framework/interface/ESHandle.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 
00007 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
00008 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
00009 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
00010 #include "DataFormats/TrackReco/interface/Track.h"
00011 #include "DataFormats/MuonReco/interface/Muon.h"
00012 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00013 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
00014 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00015 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
00016 #include "Geometry/CommonTopologies/interface/StripTopology.h"
00017 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
00018 #include "Geometry/RPCGeometry/interface/RPCRollSpecs.h"
00019 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
00020 #include "Geometry/RPCGeometry/interface/RPCGeomServ.h"
00021 #include "Geometry/Records/interface/MuonGeometryRecord.h"
00022 #include "SimGeneral/TrackingAnalysis/interface/SimHitTPAssociationProducer.h"
00023 
00024 #include <algorithm>
00025 
00026 using namespace std;
00027 
00028 typedef MonitorElement* MEP;
00029 
00030 RPCRecHitValid::RPCRecHitValid(const edm::ParameterSet& pset)
00031 {
00032   simHitLabel_ = pset.getParameter<edm::InputTag>("simHit");
00033   recHitLabel_ = pset.getParameter<edm::InputTag>("recHit");
00034   simParticleLabel_ = pset.getParameter<edm::InputTag>("simTrack");
00035   simHitAssocLabel_ = pset.getParameter<edm::InputTag>("simHitAssoc");
00036   muonLabel_ = pset.getParameter<edm::InputTag>("muon");
00037   dbe_ = edm::Service<DQMStore>().operator->();
00038   if ( !dbe_ )
00039   {
00040     edm::LogError("RPCRecHitValid") << "No DQMStore instance\n";
00041     return;
00042   }
00043 
00044   // Book MonitorElements
00045   subDir_ = pset.getParameter<std::string>("subDir");
00046   h_.bookHistograms(dbe_, subDir_);
00047 
00048   // SimHit plots, not compatible to RPCPoint-RPCRecHit comparison
00049   dbe_->setCurrentFolder(subDir_+"/HitProperty");
00050   h_simParticleType = dbe_->book1D("SimHitPType", "SimHit particle type", 11, 0, 11);
00051   h_simParticleType->getTH1()->SetMinimum(0);
00052   if ( TH1* h = h_simParticleType->getTH1() )
00053   {
00054     h->GetXaxis()->SetBinLabel(1 , "#mu^{-}");
00055     h->GetXaxis()->SetBinLabel(2 , "#mu^{+}");
00056     h->GetXaxis()->SetBinLabel(3 , "e^{-}"  );
00057     h->GetXaxis()->SetBinLabel(4 , "e^{+}"  );
00058     h->GetXaxis()->SetBinLabel(5 , "#pi^{+}");
00059     h->GetXaxis()->SetBinLabel(6 , "#pi^{-}");
00060     h->GetXaxis()->SetBinLabel(7 , "K^{+}"  );
00061     h->GetXaxis()->SetBinLabel(8 , "K^{-}"  );
00062     h->GetXaxis()->SetBinLabel(9 , "p^{+}"  );
00063     h->GetXaxis()->SetBinLabel(10, "p^{-}"  );
00064     h->GetXaxis()->SetBinLabel(11, "Other"  );
00065   }
00066 
00067   dbe_->setCurrentFolder(subDir_+"/Track");
00068 
00069   h_nRPCHitPerSimMuon        = dbe_->book1D("NRPCHitPerSimMuon"       , "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
00070   h_nRPCHitPerSimMuonBarrel  = dbe_->book1D("NRPCHitPerSimMuonBarrel" , "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
00071   h_nRPCHitPerSimMuonOverlap = dbe_->book1D("NRPCHitPerSimMuonOverlap", "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
00072   h_nRPCHitPerSimMuonEndcap  = dbe_->book1D("NRPCHitPerSimMuonEndcap" , "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
00073 
00074   h_nRPCHitPerRecoMuon        = dbe_->book1D("NRPCHitPerRecoMuon"       , "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
00075   h_nRPCHitPerRecoMuonBarrel  = dbe_->book1D("NRPCHitPerRecoMuonBarrel" , "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
00076   h_nRPCHitPerRecoMuonOverlap = dbe_->book1D("NRPCHitPerRecoMuonOverlap", "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
00077   h_nRPCHitPerRecoMuonEndcap  = dbe_->book1D("NRPCHitPerRecoMuonEndcap" , "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
00078 
00079   h_nRPCHitPerSimMuon        ->getTH1()->SetMinimum(0);
00080   h_nRPCHitPerSimMuonBarrel  ->getTH1()->SetMinimum(0);
00081   h_nRPCHitPerSimMuonOverlap ->getTH1()->SetMinimum(0);
00082   h_nRPCHitPerSimMuonEndcap  ->getTH1()->SetMinimum(0);
00083 
00084   h_nRPCHitPerRecoMuon       ->getTH1()->SetMinimum(0);
00085   h_nRPCHitPerRecoMuonBarrel ->getTH1()->SetMinimum(0);
00086   h_nRPCHitPerRecoMuonOverlap->getTH1()->SetMinimum(0);
00087   h_nRPCHitPerRecoMuonEndcap ->getTH1()->SetMinimum(0);
00088 
00089   float ptBins[] = {0, 1, 2, 5, 10, 20, 30, 50, 100, 200, 300, 500};
00090   const int nPtBins = sizeof(ptBins)/sizeof(float)-1;
00091   h_simMuonBarrel_pt   = dbe_->book1D("SimMuonBarrel_pt"  , "SimMuon RPCHit in Barrel  p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
00092   h_simMuonOverlap_pt  = dbe_->book1D("SimMuonOverlap_pt" , "SimMuon RPCHit in Overlap p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
00093   h_simMuonEndcap_pt   = dbe_->book1D("SimMuonEndcap_pt"  , "SimMuon RPCHit in Endcap  p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
00094   h_simMuonNoRPC_pt    = dbe_->book1D("SimMuonNoRPC_pt" , "SimMuon without RPCHit p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
00095   h_simMuonBarrel_eta  = dbe_->book1D("SimMuonBarrel_eta" , "SimMuon RPCHit in Barrel  #eta;#eta", 50, -2.5, 2.5);
00096   h_simMuonOverlap_eta = dbe_->book1D("SimMuonOverlap_eta", "SimMuon RPCHit in Overlap #eta;#eta", 50, -2.5, 2.5);
00097   h_simMuonEndcap_eta  = dbe_->book1D("SimMuonEndcap_eta" , "SimMuon RPCHit in Endcap  #eta;#eta", 50, -2.5, 2.5);
00098   h_simMuonNoRPC_eta   = dbe_->book1D("SimMuonNoRPC_eta", "SimMuon without RPCHit #eta;#eta", 50, -2.5, 2.5);
00099   h_simMuonBarrel_phi  = dbe_->book1D("SimMuonBarrel_phi" , "SimMuon RPCHit in Barrel  #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
00100   h_simMuonOverlap_phi = dbe_->book1D("SimMuonOverlap_phi", "SimMuon RPCHit in Overlap #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
00101   h_simMuonEndcap_phi  = dbe_->book1D("SimMuonEndcap_phi" , "SimMuon RPCHit in Endcap  #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
00102   h_simMuonNoRPC_phi   = dbe_->book1D("SimMuonNoRPC_phi", "SimMuon without RPCHit #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
00103 
00104   h_recoMuonBarrel_pt   = dbe_->book1D("RecoMuonBarrel_pt"  , "RecoMuon RPCHit in Barrel  p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
00105   h_recoMuonOverlap_pt  = dbe_->book1D("RecoMuonOverlap_pt" , "RecoMuon RPCHit in Overlap p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
00106   h_recoMuonEndcap_pt   = dbe_->book1D("RecoMuonEndcap_pt"  , "RecoMuon RPCHit in Endcap  p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
00107   h_recoMuonNoRPC_pt    = dbe_->book1D("RecoMuonNoRPC_pt" , "RecoMuon without RPCHit p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
00108   h_recoMuonBarrel_eta  = dbe_->book1D("RecoMuonBarrel_eta" , "RecoMuon RPCHit in Barrel  #eta;#eta", 50, -2.5, 2.5);
00109   h_recoMuonOverlap_eta = dbe_->book1D("RecoMuonOverlap_eta", "RecoMuon RPCHit in Overlap #eta;#eta", 50, -2.5, 2.5);
00110   h_recoMuonEndcap_eta  = dbe_->book1D("RecoMuonEndcap_eta" , "RecoMuon RPCHit in Endcap  #eta;#eta", 50, -2.5, 2.5);
00111   h_recoMuonNoRPC_eta   = dbe_->book1D("RecoMuonNoRPC_eta", "RecoMuon without RPCHit #eta;#eta", 50, -2.5, 2.5);
00112   h_recoMuonBarrel_phi  = dbe_->book1D("RecoMuonBarrel_phi" , "RecoMuon RPCHit in Barrel  #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
00113   h_recoMuonOverlap_phi = dbe_->book1D("RecoMuonOverlap_phi", "RecoMuon RPCHit in Overlap #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
00114   h_recoMuonEndcap_phi  = dbe_->book1D("RecoMuonEndcap_phi" , "RecoMuon RPCHit in Endcap  #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
00115   h_recoMuonNoRPC_phi   = dbe_->book1D("RecoMuonNoRPC_phi", "RecoMuon without RPCHit #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
00116 
00117   h_simMuonBarrel_pt   ->getTH1()->SetMinimum(0);
00118   h_simMuonOverlap_pt  ->getTH1()->SetMinimum(0);
00119   h_simMuonEndcap_pt   ->getTH1()->SetMinimum(0);
00120   h_simMuonNoRPC_pt    ->getTH1()->SetMinimum(0);
00121   h_simMuonBarrel_eta  ->getTH1()->SetMinimum(0);
00122   h_simMuonOverlap_eta ->getTH1()->SetMinimum(0);
00123   h_simMuonEndcap_eta  ->getTH1()->SetMinimum(0);
00124   h_simMuonNoRPC_eta   ->getTH1()->SetMinimum(0);
00125   h_simMuonBarrel_phi  ->getTH1()->SetMinimum(0);
00126   h_simMuonOverlap_phi ->getTH1()->SetMinimum(0);
00127   h_simMuonEndcap_phi  ->getTH1()->SetMinimum(0);
00128   h_simMuonNoRPC_phi   ->getTH1()->SetMinimum(0);
00129 
00130   h_recoMuonBarrel_pt  ->getTH1()->SetMinimum(0);
00131   h_recoMuonOverlap_pt ->getTH1()->SetMinimum(0);
00132   h_recoMuonEndcap_pt  ->getTH1()->SetMinimum(0);
00133   h_recoMuonNoRPC_pt   ->getTH1()->SetMinimum(0);
00134   h_recoMuonBarrel_eta ->getTH1()->SetMinimum(0);
00135   h_recoMuonOverlap_eta->getTH1()->SetMinimum(0);
00136   h_recoMuonEndcap_eta ->getTH1()->SetMinimum(0);
00137   h_recoMuonNoRPC_eta  ->getTH1()->SetMinimum(0);
00138   h_recoMuonBarrel_phi ->getTH1()->SetMinimum(0);
00139   h_recoMuonOverlap_phi->getTH1()->SetMinimum(0);
00140   h_recoMuonEndcap_phi ->getTH1()->SetMinimum(0);
00141   h_recoMuonNoRPC_phi  ->getTH1()->SetMinimum(0);
00142 
00143   dbe_->setCurrentFolder(subDir_+"/Occupancy");
00144 
00145   h_eventCount = dbe_->book1D("EventCount", "Event count", 3, 1, 4);
00146   h_eventCount->getTH1()->SetMinimum(0);
00147   if ( h_eventCount )
00148   {
00149     TH1* h = h_eventCount->getTH1();
00150     h->GetXaxis()->SetBinLabel(1, "eventBegin");
00151     h->GetXaxis()->SetBinLabel(2, "eventEnd");
00152     h->GetXaxis()->SetBinLabel(3, "run");
00153   }
00154 
00155   h_refPunchOccupancyBarrel_wheel   = dbe_->book1D("RefPunchOccupancyBarrel_wheel"  , "RefPunchthrough occupancy", 5, -2.5, 2.5);
00156   h_refPunchOccupancyEndcap_disk    = dbe_->book1D("RefPunchOccupancyEndcap_disk"   , "RefPunchthrough occupancy", 9, -4.5, 4.5);
00157   h_refPunchOccupancyBarrel_station = dbe_->book1D("RefPunchOccupancyBarrel_station", "RefPunchthrough occupancy", 4,  0.5, 4.5);
00158   h_recPunchOccupancyBarrel_wheel   = dbe_->book1D("RecPunchOccupancyBarrel_wheel"  , "Punchthrough recHit occupancy", 5, -2.5, 2.5);
00159   h_recPunchOccupancyEndcap_disk    = dbe_->book1D("RecPunchOccupancyEndcap_disk"   , "Punchthrough recHit occupancy", 9, -4.5, 4.5);
00160   h_recPunchOccupancyBarrel_station = dbe_->book1D("RecPunchOccupancyBarrel_station", "Punchthrough recHit occupancy", 4,  0.5, 4.5);
00161 
00162   h_refPunchOccupancyBarrel_wheel   ->getTH1()->SetMinimum(0);
00163   h_refPunchOccupancyEndcap_disk    ->getTH1()->SetMinimum(0);
00164   h_refPunchOccupancyBarrel_station ->getTH1()->SetMinimum(0);
00165   h_recPunchOccupancyBarrel_wheel   ->getTH1()->SetMinimum(0);
00166   h_recPunchOccupancyEndcap_disk    ->getTH1()->SetMinimum(0);
00167   h_recPunchOccupancyBarrel_station ->getTH1()->SetMinimum(0);
00168 
00169   h_refPunchOccupancyBarrel_wheel_station = dbe_->book2D("RefPunchOccupancyBarrel_wheel_station", "RefPunchthrough occupancy", 5, -2.5, 2.5, 4, 0.5, 4.5);
00170   h_refPunchOccupancyEndcap_disk_ring     = dbe_->book2D("RefPunchOccupancyEndcap_disk_ring"    , "RefPunchthrough occupancy", 9, -4.5, 4.5, 4, 0.5, 4.5);
00171   h_recPunchOccupancyBarrel_wheel_station = dbe_->book2D("RecPunchOccupancyBarrel_wheel_station", "Punchthrough recHit occupancy", 5, -2.5, 2.5, 4, 0.5, 4.5);
00172   h_recPunchOccupancyEndcap_disk_ring     = dbe_->book2D("RecPunchOccupancyEndcap_disk_ring"    , "Punchthrough recHit occupancy", 9, -4.5, 4.5, 4, 0.5, 4.5);
00173 
00174   h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetOption("COLZ");
00175   h_refPunchOccupancyEndcap_disk_ring    ->getTH2F()->SetOption("COLZ");
00176   h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetOption("COLZ");
00177   h_recPunchOccupancyEndcap_disk_ring    ->getTH2F()->SetOption("COLZ");
00178 
00179   h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetContour(10);
00180   h_refPunchOccupancyEndcap_disk_ring    ->getTH2F()->SetContour(10);
00181   h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetContour(10);
00182   h_recPunchOccupancyEndcap_disk_ring    ->getTH2F()->SetContour(10);
00183 
00184   h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetStats(0);
00185   h_refPunchOccupancyEndcap_disk_ring    ->getTH2F()->SetStats(0);
00186   h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetStats(0);
00187   h_recPunchOccupancyEndcap_disk_ring    ->getTH2F()->SetStats(0);
00188 
00189   h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetMinimum(0);
00190   h_refPunchOccupancyEndcap_disk_ring    ->getTH2F()->SetMinimum(0);
00191   h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetMinimum(0);
00192   h_recPunchOccupancyEndcap_disk_ring    ->getTH2F()->SetMinimum(0);
00193 
00194   for ( int i=1; i<=5; ++i )
00195   {
00196     TString binLabel = Form("Wheel %d", i-3);
00197     h_refPunchOccupancyBarrel_wheel->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
00198     h_refPunchOccupancyBarrel_wheel_station->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
00199     h_recPunchOccupancyBarrel_wheel->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
00200     h_recPunchOccupancyBarrel_wheel_station->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
00201   }
00202 
00203   for ( int i=1; i<=9; ++i )
00204   {
00205     TString binLabel = Form("Disk %d", i-5);
00206     h_refPunchOccupancyEndcap_disk  ->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
00207     h_refPunchOccupancyEndcap_disk_ring  ->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
00208     h_recPunchOccupancyEndcap_disk  ->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
00209     h_recPunchOccupancyEndcap_disk_ring  ->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
00210   }
00211 
00212   for ( int i=1; i<=4; ++i )
00213   {
00214     TString binLabel = Form("Station %d", i);
00215     h_refPunchOccupancyBarrel_station  ->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
00216     h_refPunchOccupancyBarrel_wheel_station  ->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
00217     h_recPunchOccupancyBarrel_station  ->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
00218     h_recPunchOccupancyBarrel_wheel_station  ->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
00219   }
00220 
00221   for ( int i=1; i<=4; ++i )
00222   {
00223     TString binLabel = Form("Ring %d", i);
00224     h_refPunchOccupancyEndcap_disk_ring  ->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
00225     h_recPunchOccupancyEndcap_disk_ring  ->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
00226   }
00227 }
00228 
00229 RPCRecHitValid::~RPCRecHitValid()
00230 {
00231 }
00232 
00233 void RPCRecHitValid::beginJob()
00234 {
00235 }
00236 
00237 void RPCRecHitValid::endJob()
00238 {
00239 }
00240 
00241 void RPCRecHitValid::beginRun(const edm::Run& run, const edm::EventSetup& eventSetup)
00242 {
00243   if ( !dbe_ ) return;
00244   h_eventCount->Fill(3);
00245 
00246   // Book roll-by-roll histograms
00247   edm::ESHandle<RPCGeometry> rpcGeom;
00248   eventSetup.get<MuonGeometryRecord>().get(rpcGeom);
00249 
00250   int nRPCRollBarrel = 0, nRPCRollEndcap = 0;
00251 
00252   TrackingGeometry::DetContainer rpcDets = rpcGeom->dets();
00253   for ( auto det : rpcDets )
00254   {
00255     RPCChamber* rpcCh = dynamic_cast<RPCChamber*>(det);
00256     if ( !rpcCh ) continue;
00257 
00258     std::vector<const RPCRoll*> rolls = rpcCh->rolls();
00259     for ( auto roll : rolls )
00260     {
00261       if ( !roll ) continue;
00262 
00263       //RPCGeomServ rpcSrv(roll->id());
00264       const int rawId = roll->geographicalId().rawId();
00265       //if ( !roll->specs()->isRPC() ) { cout << "\nNoRPC : " << rpcSrv.name() << ' ' << rawId << endl; continue; }
00266 
00267       if ( roll->isBarrel() )
00268       {
00269         detIdToIndexMapBarrel_[rawId] = nRPCRollBarrel;
00270         //rollIdToNameMapBarrel_[rawId] = rpcSrv.name();
00271         ++nRPCRollBarrel;
00272       }
00273       else
00274       {
00275         detIdToIndexMapEndcap_[rawId] = nRPCRollEndcap;
00276         //rollIdToNameMapEndcap_[rawId] = rpcSrv.name();
00277         ++nRPCRollEndcap;
00278       }
00279     }
00280   }
00281 
00282   dbe_->setCurrentFolder(subDir_+"/Occupancy");
00283   h_matchOccupancyBarrel_detId = dbe_->book1D("MatchOccupancyBarrel_detId", "Matched hit occupancy;roll index (can be arbitrary)", nRPCRollBarrel, 0, nRPCRollBarrel);
00284   h_matchOccupancyEndcap_detId = dbe_->book1D("MatchOccupancyEndcap_detId", "Matched hit occupancy;roll index (can be arbitrary)", nRPCRollEndcap, 0, nRPCRollEndcap);
00285   h_refOccupancyBarrel_detId = dbe_->book1D("RefOccupancyBarrel_detId", "Reference hit occupancy;roll index (can be arbitrary)", nRPCRollBarrel, 0, nRPCRollBarrel);
00286   h_refOccupancyEndcap_detId = dbe_->book1D("RefOccupancyEndcap_detId", "Reference hit occupancy;roll index (can be arbitrary)", nRPCRollEndcap, 0, nRPCRollEndcap);
00287   h_noiseOccupancyBarrel_detId = dbe_->book1D("NoiseOccupancyBarrel_detId", "Noise occupancy;roll index (can be arbitrary)", nRPCRollBarrel, 0, nRPCRollBarrel);
00288   h_noiseOccupancyEndcap_detId = dbe_->book1D("NoiseOccupancyEndcap_detId", "Noise occupancy;roll index (can be arbitrary)", nRPCRollEndcap, 0, nRPCRollEndcap);
00289 
00290   h_matchOccupancyBarrel_detId->getTH1()->SetMinimum(0);
00291   h_matchOccupancyEndcap_detId->getTH1()->SetMinimum(0);
00292   h_refOccupancyBarrel_detId  ->getTH1()->SetMinimum(0);
00293   h_refOccupancyEndcap_detId  ->getTH1()->SetMinimum(0);
00294   h_noiseOccupancyBarrel_detId->getTH1()->SetMinimum(0);
00295   h_noiseOccupancyEndcap_detId->getTH1()->SetMinimum(0);
00296 
00297   h_rollAreaBarrel_detId = dbe_->bookProfile("RollAreaBarrel_detId", "Roll area;roll index;Area", nRPCRollBarrel, 0., 1.*nRPCRollBarrel, 0., 1e5);
00298   h_rollAreaEndcap_detId = dbe_->bookProfile("RollAreaEndcap_detId", "Roll area;roll index;Area", nRPCRollEndcap, 0., 1.*nRPCRollEndcap, 0., 1e5);
00299 
00300   for ( auto detIdToIndex : detIdToIndexMapBarrel_ )
00301   {
00302     const int rawId = detIdToIndex.first;
00303     const int index = detIdToIndex.second;
00304 
00305     const RPCDetId rpcDetId = static_cast<const RPCDetId>(rawId);
00306     const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(rpcDetId));
00307 
00308     //RPCGeomServ rpcSrv(roll->id());
00309     //if ( !roll->specs()->isRPC() ) { cout << "\nNoRPC : " << rpcSrv.name() << ' ' << rawId << endl; continue; }
00310 
00311     const StripTopology& topol = roll->specificTopology();
00312     const double area = topol.stripLength()*topol.nstrips()*topol.pitch();
00313 
00314     h_rollAreaBarrel_detId->Fill(index, area);
00315   }
00316 
00317   for ( auto detIdToIndex : detIdToIndexMapEndcap_ )
00318   {
00319     const int rawId = detIdToIndex.first;
00320     const int index = detIdToIndex.second;
00321 
00322     const RPCDetId rpcDetId = static_cast<const RPCDetId>(rawId);
00323     const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(rpcDetId));
00324 
00325     //RPCGeomServ rpcSrv(roll->id());
00326     //if ( !roll->specs()->isRPC() ) { cout << "\nNoRPC : " << rpcSrv.name() << ' ' << rawId << endl; continue; }
00327 
00328     const StripTopology& topol = roll->specificTopology();
00329     const double area = topol.stripLength()*topol.nstrips()*topol.pitch();
00330 
00331     h_rollAreaEndcap_detId->Fill(index, area);
00332   }
00333 }
00334 
00335 void RPCRecHitValid::endRun(const edm::Run& run, const edm::EventSetup& eventSetup)
00336 {
00337   if ( !dbe_ ) return;
00338 
00339 }
00340 
00341 void RPCRecHitValid::analyze(const edm::Event& event, const edm::EventSetup& eventSetup)
00342 {
00343   if ( !dbe_ ) return;
00344   h_eventCount->Fill(1);
00345 
00346   // Get the RPC Geometry
00347   edm::ESHandle<RPCGeometry> rpcGeom;
00348   eventSetup.get<MuonGeometryRecord>().get(rpcGeom);
00349 
00350   // Retrieve SimHits from the event
00351   edm::Handle<edm::PSimHitContainer> simHitHandle;
00352   if ( !event.getByLabel(simHitLabel_, simHitHandle) )
00353   {
00354     edm::LogInfo("RPCRecHitValid") << "Cannot find simHit collection\n";
00355     return;
00356   }
00357 
00358   // Retrieve RecHits from the event
00359   edm::Handle<RPCRecHitCollection> recHitHandle;
00360   if ( !event.getByLabel(recHitLabel_, recHitHandle) )
00361   {
00362     edm::LogInfo("RPCRecHitValid") << "Cannot find recHit collection\n";
00363     return;
00364   }
00365 
00366   // Get SimParticles
00367   edm::Handle<TrackingParticleCollection> simParticleHandle;
00368   if ( !event.getByLabel(simParticleLabel_, simParticleHandle) )
00369   {
00370     edm::LogInfo("RPCRecHitValid") << "Cannot find TrackingParticle collection\n";
00371     return;
00372   }
00373 
00374   typedef std::pair<TrackingParticleRef, TrackPSimHitRef> SimHitTPPair;
00375   typedef std::vector<SimHitTPPair> SimHitTPAssociationList;
00376   // Get SimParticle to SimHit association map
00377   edm::Handle<SimHitTPAssociationList> simHitsTPAssoc;
00378   if ( !event.getByLabel(simHitAssocLabel_, simHitsTPAssoc) )
00379   {
00380     edm::LogInfo("RPCRecHitValid") << "Cannot find TrackingParticle to SimHit association map\n";
00381     return;
00382   }
00383 
00384   // Get RecoMuons
00385   edm::Handle<edm::View<reco::Muon> > muonHandle;
00386   if ( !event.getByLabel(muonLabel_, muonHandle) )
00387   {
00388     edm::LogInfo("RPCRecHitValid") << "Cannot find muon collection\n";
00389     return;
00390   }
00391 
00392   typedef edm::PSimHitContainer::const_iterator SimHitIter;
00393   typedef RPCRecHitCollection::const_iterator RecHitIter;
00394   typedef std::vector<TrackPSimHitRef> SimHitRefs;
00395 
00396   // TrackingParticles with (and without) RPC simHits
00397   SimHitRefs muonSimHits, pthrSimHits;
00398 
00399   for ( int i=0, n=simParticleHandle->size(); i<n; ++i )
00400   {
00401     TrackingParticleRef simParticle(simParticleHandle, i);
00402     if ( simParticle->pt() < 1.0 or simParticle->p() < 2.5 ) continue; // globalMuon acceptance
00403 
00404     // Collect SimHits from this Tracking Particle
00405     SimHitRefs simHitsFromParticle;
00406     auto range = std::equal_range(simHitsTPAssoc->begin(), simHitsTPAssoc->end(),
00407                                   std::make_pair(simParticle, TrackPSimHitRef()),
00408                                   SimHitTPAssociationProducer::simHitTPAssociationListGreater);
00409     for ( auto simParticleToHit = range.first; simParticleToHit != range.second; ++simParticleToHit )
00410     {
00411       auto simHit = simParticleToHit->second;
00412       const DetId detId(simHit->detUnitId());
00413       if ( detId.det() != DetId::Muon or detId.subdetId() != MuonSubdetId::RPC ) continue;
00414 
00415       simHitsFromParticle.push_back(simParticleToHit->second);
00416     }
00417     const int nRPCHit = simHitsFromParticle.size();
00418     const bool hasRPCHit = nRPCHit > 0;
00419 
00420     if ( abs(simParticle->pdgId()) == 13 )
00421     {
00422       muonSimHits.insert(muonSimHits.end(), simHitsFromParticle.begin(), simHitsFromParticle.end());
00423 
00424       // Count number of Barrel hits and Endcap hits
00425       int nRPCHitBarrel = 0;
00426       int nRPCHitEndcap = 0;
00427       for ( auto simHit : simHitsFromParticle )
00428       {
00429         const RPCDetId rpcDetId = static_cast<const RPCDetId>(simHit->detUnitId());
00430         const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(rpcDetId));
00431         if ( !roll ) continue;
00432 
00433         if ( rpcDetId.region() == 0 ) ++nRPCHitBarrel;
00434         else ++nRPCHitEndcap;
00435       }
00436 
00437       // Fill TrackingParticle related histograms
00438       h_nRPCHitPerSimMuon->Fill(nRPCHit);
00439       if ( nRPCHitBarrel and nRPCHitEndcap )
00440       {
00441         h_nRPCHitPerSimMuonOverlap->Fill(nRPCHit);
00442         h_simMuonOverlap_pt->Fill(simParticle->pt());
00443         h_simMuonOverlap_eta->Fill(simParticle->eta());
00444         h_simMuonOverlap_phi->Fill(simParticle->phi());
00445       }
00446       else if ( nRPCHitBarrel )
00447       {
00448         h_nRPCHitPerSimMuonBarrel->Fill(nRPCHit);
00449         h_simMuonBarrel_pt->Fill(simParticle->pt());
00450         h_simMuonBarrel_eta->Fill(simParticle->eta());
00451         h_simMuonBarrel_phi->Fill(simParticle->phi());
00452       }
00453       else if ( nRPCHitEndcap )
00454       {
00455         h_nRPCHitPerSimMuonEndcap->Fill(nRPCHit);
00456         h_simMuonEndcap_pt->Fill(simParticle->pt());
00457         h_simMuonEndcap_eta->Fill(simParticle->eta());
00458         h_simMuonEndcap_phi->Fill(simParticle->phi());
00459       }
00460       else
00461       {
00462         h_simMuonNoRPC_pt->Fill(simParticle->pt());
00463         h_simMuonNoRPC_eta->Fill(simParticle->eta());
00464         h_simMuonNoRPC_phi->Fill(simParticle->phi());
00465       }
00466     }
00467     else
00468     {
00469       pthrSimHits.insert(pthrSimHits.end(), simHitsFromParticle.begin(), simHitsFromParticle.end());
00470     }
00471 
00472     if ( hasRPCHit )
00473     {
00474       switch ( simParticle->pdgId() )
00475       {
00476         case    13: h_simParticleType->Fill( 0); break;
00477         case   -13: h_simParticleType->Fill( 1); break;
00478         case    11: h_simParticleType->Fill( 2); break;
00479         case   -11: h_simParticleType->Fill( 3); break;
00480         case   211: h_simParticleType->Fill( 4); break;
00481         case  -211: h_simParticleType->Fill( 5); break;
00482         case   321: h_simParticleType->Fill( 6); break;
00483         case  -321: h_simParticleType->Fill( 7); break;
00484         case  2212: h_simParticleType->Fill( 8); break;
00485         case -2212: h_simParticleType->Fill( 9); break;
00486         default:    h_simParticleType->Fill(10); break;
00487       }
00488     }
00489   }
00490 
00491   // Loop over muon simHits, fill histograms which does not need associations
00492   int nRefHitBarrel = 0, nRefHitEndcap = 0;
00493   for ( auto simHit : muonSimHits )
00494   {
00495     const RPCDetId detId = static_cast<const RPCDetId>(simHit->detUnitId());
00496     const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(detId));
00497 
00498     const int region = roll->id().region();
00499     const int ring = roll->id().ring();
00500     //const int sector = roll->id().sector();
00501     const int station = roll->id().station();
00502     //const int layer = roll->id().layer();
00503     //const int subSector = roll->id().subsector();
00504 
00505     if ( region == 0 )
00506     {
00507       ++nRefHitBarrel;
00508       h_.refHitOccupancyBarrel_wheel->Fill(ring);
00509       h_.refHitOccupancyBarrel_station->Fill(station);
00510       h_.refHitOccupancyBarrel_wheel_station->Fill(ring, station);
00511 
00512       h_refOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[simHit->detUnitId()]);
00513     }
00514     else
00515     {
00516       ++nRefHitEndcap;
00517       h_.refHitOccupancyEndcap_disk->Fill(region*station);
00518       h_.refHitOccupancyEndcap_disk_ring->Fill(region*station, ring);
00519 
00520       h_refOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[simHit->detUnitId()]);
00521     }
00522   }
00523 
00524   // Loop over punch-through simHits, fill histograms which does not need associations
00525   for ( auto simHit : pthrSimHits )
00526   {
00527     const RPCDetId detId = static_cast<const RPCDetId>(simHit->detUnitId());
00528     const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(detId()));
00529 
00530     const int region = roll->id().region();
00531     const int ring = roll->id().ring();
00532     //const int sector = roll->id().sector();
00533     const int station = roll->id().station();
00534     //const int layer = roll->id().layer();
00535     //const int subSector = roll->id().subsector();
00536 
00537     if ( region == 0 )
00538     {
00539       ++nRefHitBarrel;
00540       h_refPunchOccupancyBarrel_wheel->Fill(ring);
00541       h_refPunchOccupancyBarrel_station->Fill(station);
00542       h_refPunchOccupancyBarrel_wheel_station->Fill(ring, station);
00543 
00544       h_refOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[simHit->detUnitId()]);
00545     }
00546     else
00547     {
00548       ++nRefHitEndcap;
00549       h_refPunchOccupancyEndcap_disk->Fill(region*station);
00550       h_refPunchOccupancyEndcap_disk_ring->Fill(region*station, ring);
00551 
00552       h_refOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[simHit->detUnitId()]);
00553     }
00554   }
00555   h_.nRefHitBarrel->Fill(nRefHitBarrel);
00556   h_.nRefHitEndcap->Fill(nRefHitEndcap);
00557 
00558   // Loop over recHits, fill histograms which does not need associations
00559   int sumClusterSizeBarrel = 0, sumClusterSizeEndcap = 0;
00560   int nRecHitBarrel = 0, nRecHitEndcap = 0;
00561   for ( RecHitIter recHitIter = recHitHandle->begin();
00562         recHitIter != recHitHandle->end(); ++recHitIter )
00563   {
00564     const RPCDetId detId = static_cast<const RPCDetId>(recHitIter->rpcId());
00565     const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(detId()));
00566     if ( !roll ) continue;
00567 
00568     const int region = roll->id().region();
00569     const int ring = roll->id().ring();
00570     //const int sector = roll->id().sector();
00571     const int station = roll->id().station();
00572     //const int layer = roll->id().layer();
00573     //const int subSector = roll->id().subsector();
00574 
00575     h_.clusterSize->Fill(recHitIter->clusterSize());
00576 
00577     if ( region == 0 )
00578     {
00579       ++nRecHitBarrel;
00580       sumClusterSizeBarrel += recHitIter->clusterSize();
00581       h_.clusterSizeBarrel->Fill(recHitIter->clusterSize());
00582       h_.recHitOccupancyBarrel_wheel->Fill(ring);
00583       h_.recHitOccupancyBarrel_station->Fill(station);
00584       h_.recHitOccupancyBarrel_wheel_station->Fill(ring, station);
00585     }
00586     else
00587     {
00588       ++nRecHitEndcap;
00589       sumClusterSizeEndcap += recHitIter->clusterSize();
00590       h_.clusterSizeEndcap->Fill(recHitIter->clusterSize());
00591       h_.recHitOccupancyEndcap_disk->Fill(region*station);
00592       h_.recHitOccupancyEndcap_disk_ring->Fill(region*station, ring);
00593     }
00594 
00595   }
00596   const double nRecHit = nRecHitBarrel+nRecHitEndcap;
00597   h_.nRecHitBarrel->Fill(nRecHitBarrel);
00598   h_.nRecHitEndcap->Fill(nRecHitEndcap);
00599   if ( nRecHit > 0 )
00600   {
00601     const int sumClusterSize = sumClusterSizeBarrel+sumClusterSizeEndcap;
00602     h_.avgClusterSize->Fill(double(sumClusterSize)/nRecHit);
00603 
00604     if ( nRecHitBarrel > 0 )
00605     {
00606       h_.avgClusterSizeBarrel->Fill(double(sumClusterSizeBarrel)/nRecHitBarrel);
00607     }
00608     if ( nRecHitEndcap > 0 )
00609     {
00610       h_.avgClusterSizeEndcap->Fill(double(sumClusterSizeEndcap)/nRecHitEndcap);
00611     }
00612   }
00613 
00614   // Start matching SimHits to RecHits
00615   typedef std::map<TrackPSimHitRef, RecHitIter> SimToRecHitMap;
00616   SimToRecHitMap simToRecHitMap;
00617 
00618   for ( auto simHit : muonSimHits )
00619   {
00620     const RPCDetId simDetId = static_cast<const RPCDetId>(simHit->detUnitId());
00621     //const RPCRoll* simRoll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(simDetId));
00622 
00623     const double simX = simHit->localPosition().x();
00624 
00625     for ( RecHitIter recHitIter = recHitHandle->begin();
00626           recHitIter != recHitHandle->end(); ++recHitIter )
00627     {
00628       const RPCDetId recDetId = static_cast<const RPCDetId>(recHitIter->rpcId());
00629       const RPCRoll* recRoll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(recDetId));
00630       if ( !recRoll ) continue;
00631 
00632       if ( simDetId != recDetId ) continue;
00633 
00634       const double recX = recHitIter->localPosition().x();
00635       const double newDx = fabs(recX - simX);
00636 
00637       // Associate SimHit to RecHit
00638       SimToRecHitMap::const_iterator prevSimToReco = simToRecHitMap.find(simHit);
00639       if ( prevSimToReco == simToRecHitMap.end() )
00640       {
00641         simToRecHitMap.insert(std::make_pair(simHit, recHitIter));
00642       }
00643       else
00644       {
00645         const double oldDx = fabs(prevSimToReco->second->localPosition().x() - simX);
00646 
00647         if ( newDx < oldDx )
00648         {
00649           simToRecHitMap[simHit] = recHitIter;
00650         }
00651       }
00652     }
00653   }
00654 
00655   // Now we have simHit-recHit mapping
00656   // So we can fill up relavant histograms
00657   int nMatchHitBarrel = 0, nMatchHitEndcap = 0;
00658   for ( auto match : simToRecHitMap )
00659   {
00660     TrackPSimHitRef simHit = match.first;
00661     RecHitIter recHitIter = match.second;
00662 
00663     const RPCDetId detId = static_cast<const RPCDetId>(simHit->detUnitId());
00664     const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(detId));
00665 
00666     const int region = roll->id().region();
00667     const int ring = roll->id().ring();
00668     //const int sector = roll->id().sector();
00669     const int station = roll->id().station();
00670     //const int layer = roll->id().layer();
00671     //const int subsector = roll->id().subsector();
00672 
00673     const double simX = simHit->localPosition().x();
00674     const double recX = recHitIter->localPosition().x();
00675     const double errX = sqrt(recHitIter->localPositionError().xx());
00676     const double dX = recX - simX;
00677     const double pull = errX == 0 ? -999 : dX/errX;
00678 
00679     //const GlobalPoint simPos = roll->toGlobal(simHitIter->localPosition());
00680     //const GlobalPoint recPos = roll->toGlobal(recHitIter->localPosition());
00681 
00682     if ( region == 0 )
00683     {
00684       ++nMatchHitBarrel;
00685       h_.resBarrel->Fill(dX);
00686       h_.pullBarrel->Fill(pull);
00687       h_.matchOccupancyBarrel_wheel->Fill(ring);
00688       h_.matchOccupancyBarrel_station->Fill(station);
00689       h_.matchOccupancyBarrel_wheel_station->Fill(ring, station);
00690 
00691       h_.res_wheel_res->Fill(ring, dX);
00692       h_.res_station_res->Fill(station, dX);
00693       h_.pull_wheel_pull->Fill(ring, pull);
00694       h_.pull_station_pull->Fill(station, pull);
00695 
00696       h_matchOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[detId.rawId()]);
00697     }
00698     else
00699     {
00700       ++nMatchHitEndcap;
00701       h_.resEndcap->Fill(dX);
00702       h_.pullEndcap->Fill(pull);
00703       h_.matchOccupancyEndcap_disk->Fill(region*station);
00704       h_.matchOccupancyEndcap_disk_ring->Fill(region*station, ring);
00705 
00706       h_.res_disk_res->Fill(region*station, dX);
00707       h_.res_ring_res->Fill(ring, dX);
00708       h_.pull_disk_pull->Fill(region*station, pull);
00709       h_.pull_ring_pull->Fill(ring, pull);
00710 
00711       h_matchOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[detId.rawId()]);
00712     }
00713 
00714   }
00715   h_.nMatchHitBarrel->Fill(nMatchHitBarrel);
00716   h_.nMatchHitEndcap->Fill(nMatchHitEndcap);
00717 
00718   // Reco Muon hits
00719   for ( edm::View<reco::Muon>::const_iterator muon = muonHandle->begin();
00720         muon != muonHandle->end(); ++muon )
00721   {
00722     if ( !muon->isGlobalMuon() ) continue;
00723 
00724     int nRPCHitBarrel = 0;
00725     int nRPCHitEndcap = 0;
00726 
00727     const reco::TrackRef glbTrack = muon->globalTrack();
00728     for ( trackingRecHit_iterator recHit = glbTrack->recHitsBegin();
00729           recHit != glbTrack->recHitsEnd(); ++recHit )
00730     {
00731       if ( !(*recHit)->isValid() ) continue;
00732       const DetId detId = (*recHit)->geographicalId();
00733       if ( detId.det() != DetId::Muon or detId.subdetId() != MuonSubdetId::RPC ) continue;
00734       const RPCDetId rpcDetId = static_cast<const RPCDetId>(detId);
00735 
00736       if ( rpcDetId.region() == 0 ) ++nRPCHitBarrel;
00737       else ++nRPCHitEndcap;
00738     }
00739 
00740     const int nRPCHit = nRPCHitBarrel + nRPCHitEndcap;
00741     h_nRPCHitPerRecoMuon->Fill(nRPCHit);
00742     if ( nRPCHitBarrel and nRPCHitEndcap )
00743     {
00744       h_nRPCHitPerRecoMuonOverlap->Fill(nRPCHit);
00745       h_recoMuonOverlap_pt->Fill(muon->pt());
00746       h_recoMuonOverlap_eta->Fill(muon->eta());
00747       h_recoMuonOverlap_phi->Fill(muon->phi());
00748     }
00749     else if ( nRPCHitBarrel )
00750     {
00751       h_nRPCHitPerRecoMuonBarrel->Fill(nRPCHit);
00752       h_recoMuonBarrel_pt->Fill(muon->pt());
00753       h_recoMuonBarrel_eta->Fill(muon->eta());
00754       h_recoMuonBarrel_phi->Fill(muon->phi());
00755     }
00756     else if ( nRPCHitEndcap )
00757     {
00758       h_nRPCHitPerRecoMuonEndcap->Fill(nRPCHit);
00759       h_recoMuonEndcap_pt->Fill(muon->pt());
00760       h_recoMuonEndcap_eta->Fill(muon->eta());
00761       h_recoMuonEndcap_phi->Fill(muon->phi());
00762     }
00763     else
00764     {
00765       h_recoMuonNoRPC_pt->Fill(muon->pt());
00766       h_recoMuonNoRPC_eta->Fill(muon->eta());
00767       h_recoMuonNoRPC_phi->Fill(muon->phi());
00768     }
00769   }
00770 
00771   // Find Non-muon hits
00772   for ( RecHitIter recHitIter = recHitHandle->begin();
00773         recHitIter != recHitHandle->end(); ++recHitIter )
00774   {
00775     const RPCDetId detId = static_cast<const RPCDetId>(recHitIter->rpcId());
00776     const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(detId));
00777 
00778     const int region = roll->id().region();
00779     const int ring = roll->id().ring();
00780     //const int sector = roll->id().sector();
00781     const int station = roll->id().station();
00782     //const int layer = roll->id().layer();
00783     //const int subsector = roll->id().subsector();
00784 
00785     bool matched = false;
00786     for ( auto match : simToRecHitMap )
00787     {
00788       if ( recHitIter == match.second )
00789       {
00790         matched = true;
00791         break;
00792       }
00793     }
00794 
00795     if ( !matched )
00796     {
00797       // FIXME : kept for backward compatibility //
00798       if ( region == 0 )
00799       {
00800         h_.umOccupancyBarrel_wheel->Fill(ring);
00801         h_.umOccupancyBarrel_station->Fill(station);
00802         h_.umOccupancyBarrel_wheel_station->Fill(ring, station);
00803       }
00804       else
00805       {
00806         h_.umOccupancyEndcap_disk->Fill(region*station);
00807         h_.umOccupancyEndcap_disk_ring->Fill(region*station, ring);
00808       }
00809       // End of FIXME //
00810 
00811       int nPunchMatched = 0;
00812       // Check if this recHit came from non-muon simHit
00813       for ( auto simHit : pthrSimHits )
00814       {
00815         const int absSimHitPType = abs(simHit->particleType());
00816         if ( absSimHitPType == 13 ) continue;
00817 
00818         const RPCDetId simDetId = static_cast<const RPCDetId>(simHit->detUnitId());
00819         if ( simDetId == detId ) ++nPunchMatched;
00820       }
00821 
00822       if ( nPunchMatched > 0 )
00823       {
00824         if ( region == 0 )
00825         {
00826           h_recPunchOccupancyBarrel_wheel->Fill(ring);
00827           h_recPunchOccupancyBarrel_station->Fill(station);
00828           h_recPunchOccupancyBarrel_wheel_station->Fill(ring, station);
00829         }
00830         else
00831         {
00832           h_recPunchOccupancyEndcap_disk->Fill(region*station);
00833           h_recPunchOccupancyEndcap_disk_ring->Fill(region*station, ring);
00834         }
00835       }
00836     }
00837   }
00838 
00839   // Find noise recHits : RecHits without SimHit match
00840   for ( RecHitIter recHitIter = recHitHandle->begin();
00841         recHitIter != recHitHandle->end(); ++recHitIter )
00842   {
00843     const RPCDetId recDetId = static_cast<const RPCDetId>(recHitIter->rpcId());
00844     const RPCRoll* roll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(recDetId));
00845 
00846     const int region = roll->id().region();
00847     //    const int ring = roll->id().ring(); // UNUSED VARIABLE
00848     //const int sector = roll->id().sector();
00849     //    const int station = roll->id().station(); // UNUSED VARIABLE
00850     //const int layer = roll->id().layer();
00851     //const int subsector = roll->id().subsector();
00852 
00853     const double recX = recHitIter->localPosition().x();
00854     const double recErrX = sqrt(recHitIter->localPositionError().xx());
00855 
00856     bool matched = false;
00857     for ( SimHitIter simHitIter = simHitHandle->begin();
00858           simHitIter != simHitHandle->end(); ++simHitIter )
00859     {
00860       const RPCDetId simDetId = static_cast<const RPCDetId>(simHitIter->detUnitId());
00861       const RPCRoll* simRoll = dynamic_cast<const RPCRoll*>(rpcGeom->roll(simDetId));
00862       if ( !simRoll ) continue;
00863 
00864       if ( simDetId != recDetId ) continue;
00865 
00866       const double simX = simHitIter->localPosition().x();
00867       const double dX = fabs(recX-simX);
00868 
00869       if ( dX/recErrX < 5 )
00870       {
00871         matched = true;
00872         break;
00873       }
00874     }
00875 
00876     if ( !matched )
00877     {
00878       if ( region == 0 )
00879       {
00880         h_noiseOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[recDetId.rawId()]);
00881       }
00882       else
00883       {
00884         h_noiseOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[recDetId.rawId()]);
00885       }
00886     }
00887   }
00888 
00889   h_eventCount->Fill(2);
00890 }
00891 
00892 DEFINE_FWK_MODULE(RPCRecHitValid);
00893
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/Validation/RPCRecHits/src/RPCRecHitValid.cc
