00001
00002
00003
00004
00005
00013
00014
00015
00016
00017
00018
00019
00020 #include <memory>
00021 #include <iostream>
00022 #include <vector>
00023 #include <map>
00024 #include <string>
00025 #include <iomanip>
00026 #include <fstream>
00027 #include <ctime>
00028
00029
00030 #include "FWCore/Framework/interface/Frameworkfwd.h"
00031 #include "FWCore/Framework/interface/EDFilter.h"
00032
00033 #include "FWCore/Framework/interface/Event.h"
00034 #include "FWCore/Framework/interface/MakerMacros.h"
00035
00036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00037
00038 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
00039 #include <FWCore/ServiceRegistry/interface/Service.h>
00040 #include <FWCore/MessageLogger/interface/MessageLogger.h>
00041 #include <DataFormats/RPCRecHit/interface/RPCRecHit.h>
00042 #include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h"
00043 #include <DataFormats/RPCDigi/interface/RPCDigiCollection.h>
00044 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
00045
00046 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
00047 #include <Geometry/RPCGeometry/interface/RPCRoll.h>
00048 #include <Geometry/Records/interface/MuonGeometryRecord.h>
00049
00050 #include "EventFilter/RPCRawToDigi/interface/RPCRawDataCounts.h"
00051 #include "EventFilter/RPCRawToDigi/interface/RPCRecordFormatter.h"
00052
00053
00054 #include "CondFormats/RPCObjects/interface/RPCReadOutMapping.h"
00055 #include "CondFormats/RPCObjects/interface/RPCEMap.h"
00056 #include "CondFormats/DataRecord/interface/RPCEMapRcd.h"
00057
00058 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
00059 #include "DataFormats/CSCDigi/interface/CSCWireDigi.h"
00060 #include "DataFormats/CSCDigi/interface/CSCWireDigiCollection.h"
00061 #include "DataFormats/CSCDigi/interface/CSCStripDigi.h"
00062 #include "DataFormats/CSCDigi/interface/CSCStripDigiCollection.h"
00063
00064 #include "DataFormats/DTDigi/interface/DTDigiCollection.h"
00065 #include "CondFormats/DTObjects/interface/DTT0.h"
00066
00067 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00068 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00069 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
00070 #include "DataFormats/GeometryVector/interface/LocalVector.h"
00071
00072 #include "TVector3.h"
00073 #include "TH1F.h"
00074 #include "TH2F.h"
00075 #include "TFile.h"
00076 #include "TString.h"
00077 #include "TTree.h"
00078 #include "TProfile.h"
00079
00080 using namespace std;
00081 using namespace edm;
00082
00083
00084
00085
00086
00087
00088 class RPCNoise : public edm::EDFilter {
00089 public:
00090 explicit RPCNoise(const edm::ParameterSet&);
00091 ~RPCNoise();
00092
00093 private:
00094 virtual void beginJob() ;
00095 virtual bool filter(edm::Event&, const edm::EventSetup&) override;
00096 virtual void endJob() ;
00097
00098
00099 int nEventsAnalyzed;
00100 int nEventsSelected;
00101 int iRun;
00102 int iEvent;
00103 int firstOrbit;
00104 int lastOrbit;
00105 int thisOrbit;
00106
00107 bool fillHistograms;
00108 int nRPCHitsCut;
00109 int nCSCStripsCut;
00110 int nCSCWiresCut;
00111 int nDTDigisCut;
00112
00113 std::string histogramFileName;
00114
00115 TFile *theHistogramFile;
00116
00117 TH1F *nWires;
00118 TH1F *nStrips;
00119 TH1F *nWiresH;
00120 TH1F *nStripsH;
00121 TH1F *nDTDigis;
00122 TH1F *nDTDigisH;
00123 TH1F *t0All;
00124 TH1F *t0AllH;
00125 TH1F *nDTDigisIn;
00126 TH1F *nDTDigisInH;
00127 TH1F *nDTDigisOut;
00128 TH1F *nDTDigisOutH;
00129 TH1F *fDTDigisOut;
00130 TH1F *fDTDigisOutH;
00131 TH1F *nRPCRecHits;
00132 TH1F *nRPCRecHitsLong;
00133 TProfile *hitsVsSerial;
00134 TProfile *orbitVsSerial;
00135 TProfile *hitsVsOrbit;
00136 TH1F *dOrbit;
00137 TH1F *RPCBX;
00138 TH1F *RPCClSize;
00139 TH1F *RPCBXH;
00140 TH1F *RPCClSizeH;
00141 TH1F *rpcStation;
00142 TH1F *rpcStationH;
00143 TH1F *rpcRing;
00144 TH1F *rpcRingH;
00145 TH1F *rpcSector;
00146 TH1F *rpcSectorH;
00147 TH1F *rpcLayer;
00148 TH1F *rpcLayerH;
00149 TProfile *rpcStationVsOrbit;
00150 TProfile *rpcSectorVsOrbit;
00151 TProfile *rpcRingVsOrbit;
00152 TH1F *rpcCorner;
00153 TH1F *rpcCornerH;
00154 TProfile *rpcCornerVsOrbit;
00155
00156 };
00157
00158 RPCNoise::RPCNoise(const edm::ParameterSet& pset)
00159 {
00160 histogramFileName = pset.getUntrackedParameter<std::string>("histogramFileName","histos.root");
00161 fillHistograms = pset.getUntrackedParameter<bool>("fillHistograms",true);
00162 nRPCHitsCut = pset.getUntrackedParameter<int>("nRPCHitsCut",40);
00163 nCSCStripsCut = pset.getUntrackedParameter<int>("nCSCStripsCut",50);
00164 nCSCWiresCut = pset.getUntrackedParameter<int>("nCSCWiresCut",10);
00165 nDTDigisCut = pset.getUntrackedParameter<int>("nDTDigisCut",10);
00166 }
00167 RPCNoise::~RPCNoise()
00168 {
00169 }
00170
00171
00172 void
00173 RPCNoise::beginJob()
00174 {
00175
00176 nEventsAnalyzed = 0;
00177 nEventsSelected = 0;
00178 iRun = 0;
00179 iEvent = 0;
00180 firstOrbit = lastOrbit = thisOrbit = 0;
00181
00182 if (fillHistograms) {
00183
00184 theHistogramFile = new TFile(histogramFileName.c_str(), "RECREATE");
00185 theHistogramFile->cd();
00186
00187 nWires = new TH1F("nWires","number of wire digis", 121, -0.5, 120.5);
00188 nStrips = new TH1F("nStrips","number of strip digis", 201, -0.5, 200.5);
00189 nWiresH = new TH1F("nWiresH","number of wire digis HIGH", 121, -0.5, 120.5);
00190 nStripsH = new TH1F("nStripsH","number of strip digis HIGH", 201, -0.5, 200.5);
00191 nDTDigis = new TH1F("nDTDigis","number of DT digis",201,-0.5,200.5);
00192 nDTDigisH = new TH1F("nDTDigisH","number of DT digis HIGH",201,-0.5,200.5);
00193 nDTDigisIn = new TH1F("nDTDigisIn","N DT digis in window",75,0.,150.);
00194 nDTDigisInH = new TH1F("nDTDigisInH","N DT digis in window HIGH",75,0.,150.);
00195 nDTDigisOut = new TH1F("nDTDigisOut","N DT digis out window",75,0.,150.);
00196 nDTDigisOutH = new TH1F("nDTDigisOutH","N DT digis out window HIGH",75,0.,150.);
00197 fDTDigisOut = new TH1F("fDTDigisOut", "fraction DT digis outside window",55,0.,1.1);
00198 fDTDigisOutH = new TH1F("fDTDigisOutH","fraction DT digis outside window HIGH",55,0.,1.1);
00199
00200 t0All = new TH1F("t0All","t0",700,0.,7000.);
00201 t0AllH = new TH1F("t0AllH","t0 HIGH",700,0.,7000.);
00202 RPCBX = new TH1F("RPCBX","RPC BX",21,-10.5,10.5);
00203 RPCBXH = new TH1F("RPCBXH","RPC BX HIGH",21,-10.5,10.5);
00204 RPCClSize = new TH1F("RPCClSize","RPC cluster size",61,-0.5,60.5);
00205 RPCClSizeH = new TH1F("RPCClSizeH","RPC cluster size HIGH",61,-0.5,60.5);
00206
00207 nRPCRecHits = new TH1F("nRPCRecHits","number of RPC RecHits",101,-0.5,100.5);
00208 nRPCRecHitsLong = new TH1F("nRPCRecHitsLong","number of RPC RecHits",601,-0.5,600.5);
00209 hitsVsSerial = new TProfile("hitsVsSerial","mean RPC hits vs serial event number",4000,0.,40000.,0.,1000.);
00210 orbitVsSerial = new TProfile("orbitVsSerial","relative orbit number vs serial event number",4000,0.,40000.,0.,1.e10);
00211 hitsVsOrbit = new TProfile("hitsVsOrbit","mean RPC hits vs orbit number",3000,0.,1200000.,0.,1000.);
00212 dOrbit = new TH1F("dOrbit","difference in orbit number",121,-0.5,120.5);
00213
00214 rpcStation = new TH1F("rpcStation", "RPC station", 6,-0.5,5.5);
00215 rpcStationH = new TH1F("rpcStationH", "RPC station HIGH", 6,-0.5,5.5);
00216 rpcRing = new TH1F("rpcRing", "RPC ring", 9,-4.5,4.5);
00217 rpcRingH = new TH1F("rpcRingH", "RPC ring HIGH", 9,-4.5,4.5);
00218 rpcSector = new TH1F("rpcSector", "RPC sector", 15,-0.5,14.5);
00219 rpcSectorH = new TH1F("rpcSectorH", "RPC sector HIGH", 15,-0.5,14.5);
00220 rpcLayer = new TH1F("rpcLayer", "RPC layer", 4,-0.5,3.5);
00221 rpcLayerH = new TH1F("rpcLayerH", "RPC layer HIGH", 4,-0.5,3.5);
00222 rpcStationVsOrbit = new TProfile("rpcStationVsOrbit","mean RPC station vs. Orbit",3000,0.,1200000.,0.,20.);
00223 rpcSectorVsOrbit = new TProfile("rpcSectorVsOrbit","mean RPC sector vs. Orbit", 3000,0.,1200000.,0.,20.);
00224 rpcRingVsOrbit = new TProfile("rpcRingVsOrbit","mean RPC ring vs. Orbit", 3000,0.,1200000.,-20.,20.);
00225 rpcCorner = new TH1F("rpcCorner", "special corner designation", 4,-0.5,3.5);
00226 rpcCornerH = new TH1F("rpcCornerH","special corner designation HIGH", 4,-0.5,3.5);
00227 rpcCornerVsOrbit = new TProfile("rpcCornerVsOrbit","special corner vs. Orbit",3000,0.,1200000.,-20.,20.);
00228 }
00229
00230 }
00231
00232
00233 void
00234 RPCNoise::endJob() {
00235 std::cout << "\n\t===============================================================\n"
00236 << "\tnumber of events analyzed = " << nEventsAnalyzed << std::endl
00237 << "\tnumber of events selected = " << nEventsSelected << std::endl
00238 << "\tfirst and last orbit number : " << firstOrbit << ", " << lastOrbit << ", " << lastOrbit-firstOrbit << std::endl
00239 << "\t===============================================================\n\n";
00240
00241
00242 if (fillHistograms) {
00243
00244 printf("\n\n======= write out my histograms ====\n\n");
00245 theHistogramFile->cd();
00246 nWires->Write();
00247 nStrips->Write();
00248 nWiresH->Write();
00249 nStripsH->Write();
00250 nDTDigis->Write();
00251 nDTDigisH->Write();
00252 nDTDigisIn->Write();
00253 nDTDigisInH->Write();
00254 nDTDigisOut->Write();
00255 nDTDigisOutH->Write();
00256 fDTDigisOut->Write();
00257 fDTDigisOutH->Write();
00258 nRPCRecHits->Write();
00259 nRPCRecHitsLong->Write();
00260 hitsVsSerial->Write();
00261 hitsVsOrbit->Write();
00262 orbitVsSerial->Write();
00263 t0All->Write();
00264 t0AllH->Write();
00265 RPCBX->Write();
00266 RPCClSize->Write();
00267 RPCBXH->Write();
00268 RPCClSizeH->Write();
00269 rpcStation->Write();
00270 rpcStationH->Write();
00271 rpcRing->Write();
00272 rpcRingH->Write();
00273 rpcSector->Write();
00274 rpcSectorH->Write();
00275 rpcLayer->Write();
00276 rpcLayerH->Write();
00277 dOrbit->Write();
00278 rpcStationVsOrbit->Write();
00279 rpcSectorVsOrbit->Write();
00280 rpcRingVsOrbit->Write();
00281 rpcCorner->Write();
00282 rpcCornerH->Write();
00283 rpcCornerVsOrbit->Write();
00284 theHistogramFile->Close();
00285 }
00286 }
00287
00288
00289
00290
00291 bool
00292 RPCNoise::filter(edm::Event& event, const edm::EventSetup& eventSetup){
00293
00294 bool selectThisEvent = false;
00295
00296
00297 nEventsAnalyzed++;
00298
00299 iRun = event.id().run();
00300 iEvent = event.id().event();
00301
00302 bool printThisLine = (nEventsAnalyzed%100 == 0);
00303 if (printThisLine) {
00304 std::cout << "======================================"
00305 << " analyzed= " << nEventsAnalyzed
00306 << ", selected= " << nEventsSelected
00307 << "\trun,event: " << iRun << ", " << iEvent << std::endl;
00308 }
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322 int bx = event.bunchCrossing();
00323 int thisOrbit = event.orbitNumber();
00324 long mTime = 3564*thisOrbit + bx;
00325 if (firstOrbit == 0) {
00326 firstOrbit = thisOrbit;
00327 lastOrbit = thisOrbit;
00328 }
00329 int deltaOrbit = thisOrbit - lastOrbit;
00330 lastOrbit = thisOrbit;
00331 int relativeOrbit = thisOrbit - firstOrbit;
00332
00333 if (fillHistograms) {dOrbit->Fill(deltaOrbit);}
00334
00335 if (nEventsAnalyzed < 200) {
00336 std::cout << iEvent
00337
00338
00339
00340 << "\torbit,bx,mTime: " << thisOrbit << "," << bx << "," << mTime
00341 << "\tdelta= " << deltaOrbit
00342 << std::endl;
00343 }
00344
00345
00346
00347
00348
00349 edm::Handle<RPCRecHitCollection> rpcRecHits;
00350 event.getByLabel("rpcRecHits","",rpcRecHits);
00351
00352
00353 int nRPC = 0;
00354 RPCRecHitCollection::const_iterator rpcIt;
00355 for (rpcIt = rpcRecHits->begin(); rpcIt != rpcRecHits->end(); rpcIt++) {
00356
00357
00358 nRPC++;
00359 }
00360
00361
00362 for (rpcIt = rpcRecHits->begin(); rpcIt != rpcRecHits->end(); rpcIt++) {
00363 RPCDetId id = (RPCDetId)(*rpcIt).rpcId();
00364 int kRegion = id.region();
00365 int kStation = id.station();
00366 int kRing = id.ring();
00367 int kSector = id.sector();
00368 int kLayer = id.layer();
00369 int bx = (*rpcIt).BunchX();
00370 int clSize = (*rpcIt).clusterSize();
00371 int cornerFlag = 0;
00372 if ( (kStation>3) && (kSector<3) ) {
00373 cornerFlag = 1;
00374 if (kRing < 0) cornerFlag = 2;
00375 }
00376 if (nEventsAnalyzed < 100) {
00377 std::cout << "Region/Station/Ring/Sector/Layer: "
00378 << kRegion << " / "
00379 << kStation << " / "
00380 << kRing << " / "
00381 << kSector << " / "
00382 << kLayer
00383 << "\tbx,clSize: " << bx << ", " << clSize
00384 << std::endl;
00385 }
00386 if (fillHistograms) {
00387 RPCBX->Fill(bx);
00388 RPCClSize->Fill(min((float)clSize,(float)60.));
00389 rpcStation->Fill(kStation);
00390 rpcRing->Fill(kRing);
00391 rpcSector->Fill(kSector);
00392 rpcLayer->Fill(kLayer);
00393 rpcStationVsOrbit->Fill(relativeOrbit,kStation);
00394 rpcSectorVsOrbit->Fill(relativeOrbit,kSector);
00395 rpcRingVsOrbit->Fill(relativeOrbit,kRing);
00396 rpcCorner->Fill(cornerFlag);
00397 rpcCornerVsOrbit->Fill(relativeOrbit,cornerFlag);
00398 if (nRPC > nRPCHitsCut) {
00399 RPCBXH->Fill(bx);
00400 RPCClSizeH->Fill(min((float)clSize,(float)60.));
00401 rpcStationH->Fill(kStation);
00402 rpcRingH->Fill(kRing);
00403 rpcSectorH->Fill(kSector);
00404 rpcLayerH->Fill(kLayer);
00405 rpcCornerH->Fill(cornerFlag);
00406 }
00407 }
00408
00409 }
00410
00411
00412
00413
00414
00415 edm::Handle<CSCWireDigiCollection> wires;
00416 edm::Handle<CSCStripDigiCollection> strips;
00417 event.getByLabel("muonCSCDigis","MuonCSCWireDigi",wires);
00418 event.getByLabel("muonCSCDigis","MuonCSCStripDigi",strips);
00419
00420
00421 int nW = 0;
00422 for (CSCWireDigiCollection::DigiRangeIterator jW=wires->begin(); jW!=wires->end(); jW++) {
00423 std::vector<CSCWireDigi>::const_iterator wireIterA = (*jW).second.first;
00424 std::vector<CSCWireDigi>::const_iterator lWireA = (*jW).second.second;
00425 for( ; wireIterA != lWireA; ++wireIterA) {
00426 nW++;
00427 }
00428 }
00429
00430
00431
00432
00433 int nS = 0;
00434 for (CSCStripDigiCollection::DigiRangeIterator jS=strips->begin(); jS!=strips->end(); jS++) {
00435 std::vector<CSCStripDigi>::const_iterator stripItA = (*jS).second.first;
00436 std::vector<CSCStripDigi>::const_iterator lastStripA = (*jS).second.second;
00437 for( ; stripItA != lastStripA; ++stripItA) {
00438 std::vector<int> myADCVals = stripItA->getADCCounts();
00439 int iDiff = myADCVals[4]+myADCVals[5]-myADCVals[0]-myADCVals[1];
00440 if (iDiff > 30) {
00441 nS++;
00442 }
00443 }
00444 }
00445
00446
00447
00448
00449
00450
00451 edm::Handle<DTDigiCollection> dtDIGIs;
00452 event.getByLabel("muonDTDigis",dtDIGIs);
00453
00454
00455 int nDT = 0;
00456 int nDTin = 0;
00457 int nDTout = 0;
00458 for (DTDigiCollection::DigiRangeIterator jDT=dtDIGIs->begin(); jDT!=dtDIGIs->end(); ++jDT) {
00459 const DTDigiCollection::Range& digiRange = (*jDT).second;
00460 for (DTDigiCollection::const_iterator digi = digiRange.first;
00461 digi != digiRange.second;
00462 digi++) {
00463 double t0 = (*digi).countsTDC();
00464 nDT++;
00465 if ((t0>3050) && (t0<3700)) {
00466 nDTin++;
00467 } else {
00468 nDTout++;
00469 }
00470 if (fillHistograms) {
00471 t0All->Fill(t0);
00472 if (nRPC > nRPCHitsCut) {t0AllH->Fill(t0);}
00473 }
00474 }
00475 }
00476
00477
00478
00479
00480
00481 if (nEventsAnalyzed < 1000) {std::cout << "\tnumber of CSC DIGIS = " << nW << ", " << nS
00482 << "\tDT DIGIS = " << nDT
00483 << "\tRPC Rechits = " << nRPC << std::endl;}
00484
00485 if (fillHistograms) {
00486
00487 nWires->Fill(min((float)nW,(float)120.));
00488 nStrips->Fill(min((float)nS,(float)200.));
00489
00490 nDTDigis->Fill(min((float)nDT,(float)200.));
00491 nDTDigisIn->Fill(min((float)nDTin,(float)200.));
00492 nDTDigisOut->Fill(min((float)nDTout,(float)200.));
00493 if (nDT > 0) {
00494 float fracOut = float(nDTout)/float(nDT);
00495 fDTDigisOut->Fill(fracOut);
00496 }
00497 nRPCRecHits->Fill(min((float)nRPC,(float)100.));
00498 nRPCRecHitsLong->Fill(min((float)nRPC,(float)1000.));
00499 hitsVsSerial->Fill(nEventsAnalyzed,nRPC);
00500 hitsVsOrbit->Fill(relativeOrbit,nRPC);
00501 orbitVsSerial->Fill(nEventsAnalyzed,relativeOrbit);
00502
00503 if (nRPC > nRPCHitsCut) {
00504 nWiresH->Fill(min((float)nW,(float)120.));
00505 nStripsH->Fill(min((float)nS,(float)200.));
00506 nDTDigisH->Fill(min((float)nDT,(float)200.));
00507 nDTDigisInH->Fill(min((float)nDTin,(float)200.));
00508 nDTDigisOutH->Fill(min((float)nDTout,(float)200.));
00509 if (nDT > 0) {
00510 float fracOut = float(nDTout)/float(nDT);
00511 fDTDigisOutH->Fill(fracOut);
00512 }
00513 }
00514 }
00515
00516
00517
00518 selectThisEvent = (nRPC > nRPCHitsCut) && (nW > nCSCWiresCut || nS > nCSCStripsCut) && (nDT > nDTDigisCut);
00519 if (selectThisEvent) {nEventsSelected++;}
00520
00521 return selectThisEvent;
00522 }
00523
00524
00525 DEFINE_FWK_MODULE(RPCNoise);
1.7.3