Skims/src/RPCNoise.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RPCNoise
0004 // Class:      RPCNoise
0005 //
0006 /**\class RPCNoise RPCNoise.cc RecoLocalMuon/RPCNoise/src/RPCNoise.cc
0007
0008  Description: <simple analyis of RPC noise, and event filter>
0009
0010  Implementation:
0011      <simple application of EDFilter>
0012 */
0013 //
0014 // Original Author:  Michael Henry Schmitt
0015 //         Created:  Thu Oct 30 21:31:44 CET 2008
0016 // $Id: RPCNoise.cc,v 1.3 2010/08/07 14:55:55 wmtan Exp $
0017 //
0018 //
0019 // system include files
0020 #include <memory>
0021 #include <iostream>
0022 #include <vector>
0023 #include <map>
0024 #include <string>
0025 #include <iomanip>
0026 #include <fstream>
0027 #include <ctime>
0028
0029 // user include files
0030 #include "FWCore/Framework/interface/Frameworkfwd.h"
0031 #include "FWCore/Framework/interface/one/EDFilter.h"
0032
0033 #include "FWCore/Framework/interface/Event.h"
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037
0038 #include <Geometry/CommonDetUnit/interface/GeomDet.h>  //
0039 #include <FWCore/ServiceRegistry/interface/Service.h>
0040 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0041 #include <DataFormats/RPCRecHit/interface/RPCRecHit.h>
0042 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
0043 #include <DataFormats/RPCDigi/interface/RPCDigiCollection.h>
0044 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0045
0046 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0047 #include <Geometry/RPCGeometry/interface/RPCRoll.h>
0048 #include <Geometry/Records/interface/MuonGeometryRecord.h>
0049
0050 #include "EventFilter/RPCRawToDigi/interface/RPCRawDataCounts.h"
0051 #include "EventFilter/RPCRawToDigi/interface/RPCRecordFormatter.h"
0052 //#include "EventFilter/RPCRawToDigi/interface/RPCRawSynchro.h"
0053
0054 #include "CondFormats/RPCObjects/interface/RPCReadOutMapping.h"
0055 #include "CondFormats/RPCObjects/interface/RPCEMap.h"
0056 #include "CondFormats/DataRecord/interface/RPCEMapRcd.h"
0057
0058 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0059 #include "DataFormats/CSCDigi/interface/CSCWireDigi.h"
0060 #include "DataFormats/CSCDigi/interface/CSCWireDigiCollection.h"
0061 #include "DataFormats/CSCDigi/interface/CSCStripDigi.h"
0062 #include "DataFormats/CSCDigi/interface/CSCStripDigiCollection.h"
0063
0064 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
0065 #include "CondFormats/DTObjects/interface/DTT0.h"
0066
0067 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0068 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0069 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0070 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0071
0072 #include "TVector3.h"
0073 #include "TH1F.h"
0074 #include "TH2F.h"
0075 #include "TFile.h"
0076 #include "TString.h"
0077 #include "TTree.h"
0078 #include "TProfile.h"
0079
0080 using namespace std;
0081 using namespace edm;
0082
0083 //
0084 // class declaration
0085 //
0086
0087 class RPCNoise : public edm::one::EDFilter<> {
0088 public:
0089   explicit RPCNoise(const edm::ParameterSet &);
0090   ~RPCNoise() override;
0091
0092 private:
0093   void beginJob() override;
0094   bool filter(edm::Event &, const edm::EventSetup &) override;
0095   void endJob() override;
0096
0097   // counters
0098   int nEventsAnalyzed;
0099   int nEventsSelected;
0100   int iRun;
0101   int iEvent;
0102   int firstOrbit;
0103   int lastOrbit;
0104   int thisOrbit;
0105   // control parameters
0106   bool fillHistograms;
0107   int nRPCHitsCut;
0108   int nCSCStripsCut;
0109   int nCSCWiresCut;
0110   int nDTDigisCut;
0111   // histogram
0112   std::string histogramFileName;
0113   // The root file for the histograms.
0114   TFile *theHistogramFile;
0115   // histograms
0116   TH1F *nWires;
0117   TH1F *nStrips;
0118   TH1F *nWiresH;
0119   TH1F *nStripsH;
0120   TH1F *nDTDigis;
0121   TH1F *nDTDigisH;
0122   TH1F *t0All;
0123   TH1F *t0AllH;
0124   TH1F *nDTDigisIn;
0125   TH1F *nDTDigisInH;
0126   TH1F *nDTDigisOut;
0127   TH1F *nDTDigisOutH;
0128   TH1F *fDTDigisOut;
0129   TH1F *fDTDigisOutH;
0130   TH1F *nRPCRecHits;
0131   TH1F *nRPCRecHitsLong;
0132   TProfile *hitsVsSerial;
0133   TProfile *orbitVsSerial;
0134   TProfile *hitsVsOrbit;
0135   TH1F *dOrbit;
0136   TH1F *RPCBX;
0137   TH1F *RPCClSize;
0138   TH1F *RPCBXH;
0139   TH1F *RPCClSizeH;
0140   TH1F *rpcStation;
0141   TH1F *rpcStationH;
0142   TH1F *rpcRing;
0143   TH1F *rpcRingH;
0144   TH1F *rpcSector;
0145   TH1F *rpcSectorH;
0146   TH1F *rpcLayer;
0147   TH1F *rpcLayerH;
0148   TProfile *rpcStationVsOrbit;
0149   TProfile *rpcSectorVsOrbit;
0150   TProfile *rpcRingVsOrbit;
0151   TH1F *rpcCorner;
0152   TH1F *rpcCornerH;
0153   TProfile *rpcCornerVsOrbit;
0154 };
0155
0156 RPCNoise::RPCNoise(const edm::ParameterSet &pset) {
0157   histogramFileName = pset.getUntrackedParameter<std::string>("histogramFileName", "histos.root");
0158   fillHistograms = pset.getUntrackedParameter<bool>("fillHistograms", true);
0159   nRPCHitsCut = pset.getUntrackedParameter<int>("nRPCHitsCut", 40);
0160   nCSCStripsCut = pset.getUntrackedParameter<int>("nCSCStripsCut", 50);
0161   nCSCWiresCut = pset.getUntrackedParameter<int>("nCSCWiresCut", 10);
0162   nDTDigisCut = pset.getUntrackedParameter<int>("nDTDigisCut", 10);
0163 }
0164 RPCNoise::~RPCNoise() {}
0165
0166 void RPCNoise::beginJob() {
0167   // initialize variables
0168   nEventsAnalyzed = 0;
0169   nEventsSelected = 0;
0170   iRun = 0;
0171   iEvent = 0;
0172   firstOrbit = lastOrbit = thisOrbit = 0;
0173
0174   if (fillHistograms) {
0175     // Create the root file for the histograms
0176     theHistogramFile = new TFile(histogramFileName.c_str(), "RECREATE");
0177     theHistogramFile->cd();
0178     // book histograms
0179     nWires = new TH1F("nWires", "number of wire digis", 121, -0.5, 120.5);
0180     nStrips = new TH1F("nStrips", "number of strip digis", 201, -0.5, 200.5);
0181     nWiresH = new TH1F("nWiresH", "number of wire digis HIGH", 121, -0.5, 120.5);
0182     nStripsH = new TH1F("nStripsH", "number of strip digis HIGH", 201, -0.5, 200.5);
0183     nDTDigis = new TH1F("nDTDigis", "number of DT digis", 201, -0.5, 200.5);
0184     nDTDigisH = new TH1F("nDTDigisH", "number of DT digis HIGH", 201, -0.5, 200.5);
0185     nDTDigisIn = new TH1F("nDTDigisIn", "N DT digis in window", 75, 0., 150.);
0186     nDTDigisInH = new TH1F("nDTDigisInH", "N DT digis in window HIGH", 75, 0., 150.);
0187     nDTDigisOut = new TH1F("nDTDigisOut", "N DT digis out window", 75, 0., 150.);
0188     nDTDigisOutH = new TH1F("nDTDigisOutH", "N DT digis out window HIGH", 75, 0., 150.);
0189     fDTDigisOut = new TH1F("fDTDigisOut", "fraction DT digis outside window", 55, 0., 1.1);
0190     fDTDigisOutH = new TH1F("fDTDigisOutH", "fraction DT digis outside window HIGH", 55, 0., 1.1);
0191
0192     t0All = new TH1F("t0All", "t0", 700, 0., 7000.);
0193     t0AllH = new TH1F("t0AllH", "t0 HIGH", 700, 0., 7000.);
0194     RPCBX = new TH1F("RPCBX", "RPC BX", 21, -10.5, 10.5);
0195     RPCBXH = new TH1F("RPCBXH", "RPC BX HIGH", 21, -10.5, 10.5);
0196     RPCClSize = new TH1F("RPCClSize", "RPC cluster size", 61, -0.5, 60.5);
0197     RPCClSizeH = new TH1F("RPCClSizeH", "RPC cluster size HIGH", 61, -0.5, 60.5);
0198
0199     nRPCRecHits = new TH1F("nRPCRecHits", "number of RPC RecHits", 101, -0.5, 100.5);
0200     nRPCRecHitsLong = new TH1F("nRPCRecHitsLong", "number of RPC RecHits", 601, -0.5, 600.5);
0201     hitsVsSerial = new TProfile("hitsVsSerial", "mean RPC hits vs serial event number", 4000, 0., 40000., 0., 1000.);
0202     orbitVsSerial =
0203         new TProfile("orbitVsSerial", "relative orbit number vs serial event number", 4000, 0., 40000., 0., 1.e10);
0204     hitsVsOrbit = new TProfile("hitsVsOrbit", "mean RPC hits vs orbit number", 3000, 0., 1200000., 0., 1000.);
0205     dOrbit = new TH1F("dOrbit", "difference in orbit number", 121, -0.5, 120.5);
0206
0207     rpcStation = new TH1F("rpcStation", "RPC station", 6, -0.5, 5.5);
0208     rpcStationH = new TH1F("rpcStationH", "RPC station HIGH", 6, -0.5, 5.5);
0209     rpcRing = new TH1F("rpcRing", "RPC ring", 9, -4.5, 4.5);
0210     rpcRingH = new TH1F("rpcRingH", "RPC ring HIGH", 9, -4.5, 4.5);
0211     rpcSector = new TH1F("rpcSector", "RPC sector", 15, -0.5, 14.5);
0212     rpcSectorH = new TH1F("rpcSectorH", "RPC sector HIGH", 15, -0.5, 14.5);
0213     rpcLayer = new TH1F("rpcLayer", "RPC layer", 4, -0.5, 3.5);
0214     rpcLayerH = new TH1F("rpcLayerH", "RPC layer HIGH", 4, -0.5, 3.5);
0215     rpcStationVsOrbit = new TProfile("rpcStationVsOrbit", "mean RPC station vs. Orbit", 3000, 0., 1200000., 0., 20.);
0216     rpcSectorVsOrbit = new TProfile("rpcSectorVsOrbit", "mean RPC sector vs. Orbit", 3000, 0., 1200000., 0., 20.);
0217     rpcRingVsOrbit = new TProfile("rpcRingVsOrbit", "mean RPC ring vs. Orbit", 3000, 0., 1200000., -20., 20.);
0218     rpcCorner = new TH1F("rpcCorner", "special corner designation", 4, -0.5, 3.5);
0219     rpcCornerH = new TH1F("rpcCornerH", "special corner designation HIGH", 4, -0.5, 3.5);
0220     rpcCornerVsOrbit = new TProfile("rpcCornerVsOrbit", "special corner vs. Orbit", 3000, 0., 1200000., -20., 20.);
0221   }
0222 }
0223
0224 void RPCNoise::endJob() {
0225   std::cout << "\n\t===============================================================\n"
0226             << "\tnumber of events analyzed      = " << nEventsAnalyzed << std::endl
0227             << "\tnumber of events selected      = " << nEventsSelected << std::endl
0228             << "\tfirst and last orbit number    : " << firstOrbit << ", " << lastOrbit << ", "
0229             << lastOrbit - firstOrbit << std::endl
0230             << "\t===============================================================\n\n";
0231
0232   if (fillHistograms) {
0233     // Write the histos to file
0234     printf("\n\n======= write out my histograms ====\n\n");
0235     theHistogramFile->cd();
0236     nWires->Write();
0237     nStrips->Write();
0238     nWiresH->Write();
0239     nStripsH->Write();
0240     nDTDigis->Write();
0241     nDTDigisH->Write();
0242     nDTDigisIn->Write();
0243     nDTDigisInH->Write();
0244     nDTDigisOut->Write();
0245     nDTDigisOutH->Write();
0246     fDTDigisOut->Write();
0247     fDTDigisOutH->Write();
0248     nRPCRecHits->Write();
0249     nRPCRecHitsLong->Write();
0250     hitsVsSerial->Write();
0251     hitsVsOrbit->Write();
0252     orbitVsSerial->Write();
0253     t0All->Write();
0254     t0AllH->Write();
0255     RPCBX->Write();
0256     RPCClSize->Write();
0257     RPCBXH->Write();
0258     RPCClSizeH->Write();
0259     rpcStation->Write();
0260     rpcStationH->Write();
0261     rpcRing->Write();
0262     rpcRingH->Write();
0263     rpcSector->Write();
0264     rpcSectorH->Write();
0265     rpcLayer->Write();
0266     rpcLayerH->Write();
0267     dOrbit->Write();
0268     rpcStationVsOrbit->Write();
0269     rpcSectorVsOrbit->Write();
0270     rpcRingVsOrbit->Write();
0271     rpcCorner->Write();
0272     rpcCornerH->Write();
0273     rpcCornerVsOrbit->Write();
0274     theHistogramFile->Close();
0275   }
0276 }
0277
0278 bool RPCNoise::filter(edm::Event &event, const edm::EventSetup &eventSetup) {
0279   bool selectThisEvent = false;
0280
0281   // increment counter
0282   nEventsAnalyzed++;
0283
0284   iRun = event.id().run();
0285   iEvent = event.id().event();
0286
0287   bool printThisLine = (nEventsAnalyzed % 100 == 0);
0288   if (printThisLine) {
0289     std::cout << "======================================"
0290               << " analyzed= " << nEventsAnalyzed << ", selected= " << nEventsSelected << "\trun,event: " << iRun
0291               << ", " << iEvent << std::endl;
0292   }
0293
0294   /*
0295   const edm::Timestamp jTime = event.time();
0296   unsigned int sec  = jTime.value() >> 32;
0297   unsigned int usec = 0xFFFFFFFF & jTime.value() ;
0298   double floatTime = sec + usec/(float)1000000.;
0299   */
0300
0301   // first event gives
0302   // sec = 1225315493
0303   //    orbit = 202375185
0304   //    bx = 764
0305   //    mtime = 205094517
0306   int bx = event.bunchCrossing();
0307   int thisOrbit = event.orbitNumber();
0308   long mTime = 3564 * thisOrbit + bx;
0309   if (firstOrbit == 0) {
0310     firstOrbit = thisOrbit;
0311     lastOrbit = thisOrbit;
0312   }
0313   int deltaOrbit = thisOrbit - lastOrbit;
0314   lastOrbit = thisOrbit;
0315   int relativeOrbit = thisOrbit - firstOrbit;
0316
0317   if (fillHistograms) {
0318     dOrbit->Fill(deltaOrbit);
0319   }
0320
0321   if (nEventsAnalyzed < 200) {
0322     std::cout << iEvent
0323               //     << "\tsec,usec: " << sec << ", " << usec
0324               //     << "\tfloatTime= " << std::setprecision(16) << floatTime
0325               //     << "\tctime: " << ctime(sec)
0326               << "\torbit,bx,mTime: " << thisOrbit << "," << bx << "," << mTime << "\tdelta= " << deltaOrbit
0327               << std::endl;
0328   }
0329
0330   // ================
0331   // RPC recHits
0332   // ================
0333   edm::Handle<RPCRecHitCollection> rpcRecHits;
0334   event.getByLabel("rpcRecHits", "", rpcRecHits);
0335
0336   // count the number of RPC rechits
0337   int nRPC = 0;
0338   RPCRecHitCollection::const_iterator rpcIt;
0339   for (rpcIt = rpcRecHits->begin(); rpcIt != rpcRecHits->end(); rpcIt++) {
0340     //    RPCDetId id = (RPCDetId)(*rpcIt).rpcId();
0341     //    LocalPoint rhitlocal = (*rpcIt).localPosition();
0342     nRPC++;
0343   }
0344
0345   // loop again, this time fill histograms
0346   for (rpcIt = rpcRecHits->begin(); rpcIt != rpcRecHits->end(); rpcIt++) {
0347     RPCDetId id = (RPCDetId)(*rpcIt).rpcId();
0348     int kRegion = id.region();
0349     int kStation = id.station();
0350     int kRing = id.ring();
0351     int kSector = id.sector();
0352     int kLayer = id.layer();
0353     int bx = (*rpcIt).BunchX();
0354     int clSize = (*rpcIt).clusterSize();
0355     int cornerFlag = 0;
0356     if ((kStation > 3) && (kSector < 3)) {
0357       cornerFlag = 1;
0358       if (kRing < 0)
0359         cornerFlag = 2;
0360     }
0361     if (nEventsAnalyzed < 100) {
0362       std::cout << "Region/Station/Ring/Sector/Layer: " << kRegion << " / " << kStation << " / " << kRing << " / "
0363                 << kSector << " / " << kLayer << "\tbx,clSize: " << bx << ", " << clSize << std::endl;
0364     }
0365     if (fillHistograms) {
0366       RPCBX->Fill(bx);
0367       RPCClSize->Fill(min((float)clSize, (float)60.));
0368       rpcStation->Fill(kStation);
0369       rpcRing->Fill(kRing);
0370       rpcSector->Fill(kSector);
0371       rpcLayer->Fill(kLayer);
0372       rpcStationVsOrbit->Fill(relativeOrbit, kStation);
0373       rpcSectorVsOrbit->Fill(relativeOrbit, kSector);
0374       rpcRingVsOrbit->Fill(relativeOrbit, kRing);
0375       rpcCorner->Fill(cornerFlag);
0376       rpcCornerVsOrbit->Fill(relativeOrbit, cornerFlag);
0377       if (nRPC > nRPCHitsCut) {
0378         RPCBXH->Fill(bx);
0379         RPCClSizeH->Fill(min((float)clSize, (float)60.));
0380         rpcStationH->Fill(kStation);
0381         rpcRingH->Fill(kRing);
0382         rpcSectorH->Fill(kSector);
0383         rpcLayerH->Fill(kLayer);
0384         rpcCornerH->Fill(cornerFlag);
0385       }
0386     }
0387   }
0388
0389   // ===============
0390   // CSC DIGIs
0391   // ===============
0392   edm::Handle<CSCWireDigiCollection> wires;
0393   edm::Handle<CSCStripDigiCollection> strips;
0394   event.getByLabel("muonCSCDigis", "MuonCSCWireDigi", wires);
0395   event.getByLabel("muonCSCDigis", "MuonCSCStripDigi", strips);
0396
0397   // count the number of wire digis.
0398   int nW = 0;
0399   for (CSCWireDigiCollection::DigiRangeIterator jW = wires->begin(); jW != wires->end(); jW++) {
0400     std::vector<CSCWireDigi>::const_iterator wireIterA = (*jW).second.first;
0401     std::vector<CSCWireDigi>::const_iterator lWireA = (*jW).second.second;
0402     for (; wireIterA != lWireA; ++wireIterA) {
0403       nW++;
0404     }
0405   }
0406
0407   // count the number of fired strips.
0408   // I am using a crude indicator of signal - this is fast and adequate for
0409   // this purpose, but it would be poor for actual CSC studies.
0410   int nS = 0;
0411   for (CSCStripDigiCollection::DigiRangeIterator jS = strips->begin(); jS != strips->end(); jS++) {
0412     std::vector<CSCStripDigi>::const_iterator stripItA = (*jS).second.first;
0413     std::vector<CSCStripDigi>::const_iterator lastStripA = (*jS).second.second;
0414     for (; stripItA != lastStripA; ++stripItA) {
0415       std::vector<int> myADCVals = stripItA->getADCCounts();
0416       int iDiff = myADCVals[4] + myADCVals[5] - myADCVals[0] - myADCVals[1];
0417       if (iDiff > 30) {
0418         nS++;
0419       }
0420     }
0421   }
0422
0423   // ===============
0424   // DT DIGIs
0425   // ===============
0426   // see: CalibMuon/DTCalibration/plugins/DTT0Calibration.cc
0427   edm::Handle<DTDigiCollection> dtDIGIs;
0428   event.getByLabel("muonDTDigis", dtDIGIs);
0429
0430   // count the number of digis.
0431   int nDT = 0;
0432   int nDTin = 0;
0433   int nDTout = 0;
0434   for (DTDigiCollection::DigiRangeIterator jDT = dtDIGIs->begin(); jDT != dtDIGIs->end(); ++jDT) {
0435     const DTDigiCollection::Range &digiRange = (*jDT).second;
0436     for (DTDigiCollection::const_iterator digi = digiRange.first; digi != digiRange.second; digi++) {
0437       double t0 = (*digi).countsTDC();
0438       nDT++;
0439       if ((t0 > 3050) && (t0 < 3700)) {
0440         nDTin++;
0441       } else {
0442         nDTout++;
0443       }
0444       if (fillHistograms) {
0445         t0All->Fill(t0);
0446         if (nRPC > nRPCHitsCut) {
0447           t0AllH->Fill(t0);
0448         }
0449       }
0450     }
0451   }
0452
0453   //==============
0454   // Analysis
0455   //==============
0456
0457   if (nEventsAnalyzed < 1000) {
0458     std::cout << "\tnumber of CSC DIGIS = " << nW << ", " << nS << "\tDT DIGIS = " << nDT << "\tRPC Rechits = " << nRPC
0459               << std::endl;
0460   }
0461
0462   if (fillHistograms) {
0463     nWires->Fill(min((float)nW, (float)120.));
0464     nStrips->Fill(min((float)nS, (float)200.));
0465
0466     nDTDigis->Fill(min((float)nDT, (float)200.));
0467     nDTDigisIn->Fill(min((float)nDTin, (float)200.));
0468     nDTDigisOut->Fill(min((float)nDTout, (float)200.));
0469     if (nDT > 0) {
0470       float fracOut = float(nDTout) / float(nDT);
0471       fDTDigisOut->Fill(fracOut);
0472     }
0473     nRPCRecHits->Fill(min((float)nRPC, (float)100.));
0474     nRPCRecHitsLong->Fill(min((float)nRPC, (float)1000.));
0475     hitsVsSerial->Fill(nEventsAnalyzed, nRPC);
0476     hitsVsOrbit->Fill(relativeOrbit, nRPC);
0477     orbitVsSerial->Fill(nEventsAnalyzed, relativeOrbit);
0478
0479     if (nRPC > nRPCHitsCut) {
0480       nWiresH->Fill(min((float)nW, (float)120.));
0481       nStripsH->Fill(min((float)nS, (float)200.));
0482       nDTDigisH->Fill(min((float)nDT, (float)200.));
0483       nDTDigisInH->Fill(min((float)nDTin, (float)200.));
0484       nDTDigisOutH->Fill(min((float)nDTout, (float)200.));
0485       if (nDT > 0) {
0486         float fracOut = float(nDTout) / float(nDT);
0487         fDTDigisOutH->Fill(fracOut);
0488       }
0489     }
0490   }
0491
0492   // select this event for output?
0493
0494   selectThisEvent = (nRPC > nRPCHitsCut) && (nW > nCSCWiresCut || nS > nCSCStripsCut) && (nDT > nDTDigisCut);
0495   if (selectThisEvent) {
0496     nEventsSelected++;
0497   }
0498
0499   return selectThisEvent;
0500 }
0501
0502 //define this as a plug-in
0503 DEFINE_FWK_MODULE(RPCNoise);