File indexing completed on 2023-05-08 23:19:18
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include "RecoLocalMuon/CSCValidation/src/CSCValidation.h"
0015
0016 using namespace std;
0017 using namespace edm;
0018
0019
0020
0021
0022 CSCValidation::CSCValidation(const ParameterSet& pset) {
0023
0024 rootFileName = pset.getUntrackedParameter<std::string>("rootFileName", "valHists.root");
0025 isSimulation = pset.getUntrackedParameter<bool>("isSimulation", false);
0026 writeTreeToFile = pset.getUntrackedParameter<bool>("writeTreeToFile", true);
0027 detailedAnalysis = pset.getUntrackedParameter<bool>("detailedAnalysis", false);
0028 useDigis = pset.getUntrackedParameter<bool>("useDigis", true);
0029
0030
0031 useQualityFilter = pset.getUntrackedParameter<bool>("useQualityFilter", false);
0032 pMin = pset.getUntrackedParameter<double>("pMin", 4.);
0033 chisqMax = pset.getUntrackedParameter<double>("chisqMax", 20.);
0034 nCSCHitsMin = pset.getUntrackedParameter<int>("nCSCHitsMin", 10);
0035 nCSCHitsMax = pset.getUntrackedParameter<int>("nCSCHitsMax", 25);
0036 lengthMin = pset.getUntrackedParameter<double>("lengthMin", 140.);
0037 lengthMax = pset.getUntrackedParameter<double>("lengthMax", 600.);
0038 deltaPhiMax = pset.getUntrackedParameter<double>("deltaPhiMax", 0.2);
0039 polarMax = pset.getUntrackedParameter<double>("polarMax", 0.7);
0040 polarMin = pset.getUntrackedParameter<double>("polarMin", 0.3);
0041
0042
0043 useTriggerFilter = pset.getUntrackedParameter<bool>("useTriggerFilter", false);
0044
0045
0046 rd_token = consumes<FEDRawDataCollection>(pset.getParameter<edm::InputTag>("rawDataTag"));
0047 sd_token = consumes<CSCStripDigiCollection>(pset.getParameter<edm::InputTag>("stripDigiTag"));
0048 wd_token = consumes<CSCWireDigiCollection>(pset.getParameter<edm::InputTag>("wireDigiTag"));
0049 cd_token = consumes<CSCComparatorDigiCollection>(pset.getParameter<edm::InputTag>("compDigiTag"));
0050 al_token = consumes<CSCALCTDigiCollection>(pset.getParameter<edm::InputTag>("alctDigiTag"));
0051 cl_token = consumes<CSCCLCTDigiCollection>(pset.getParameter<edm::InputTag>("clctDigiTag"));
0052 co_token = consumes<CSCCorrelatedLCTDigiCollection>(pset.getParameter<edm::InputTag>("corrlctDigiTag"));
0053 rh_token = consumes<CSCRecHit2DCollection>(pset.getParameter<edm::InputTag>("cscRecHitTag"));
0054 se_token = consumes<CSCSegmentCollection>(pset.getParameter<edm::InputTag>("cscSegTag"));
0055 sa_token = consumes<reco::TrackCollection>(pset.getParameter<edm::InputTag>("saMuonTag"));
0056 l1_token = consumes<L1MuGMTReadoutCollection>(pset.getParameter<edm::InputTag>("l1aTag"));
0057 tr_token = consumes<TriggerResults>(pset.getParameter<edm::InputTag>("hltTag"));
0058 sh_token = consumes<PSimHitContainer>(pset.getParameter<edm::InputTag>("simHitTag"));
0059
0060 pedestalsToken_ = esConsumes<CSCDBPedestals, CSCDBPedestalsRcd>();
0061 gainsToken_ = esConsumes<CSCDBGains, CSCDBGainsRcd>();
0062 noiseToken_ = esConsumes<CSCDBNoiseMatrix, CSCDBNoiseMatrixRcd>();
0063 crosstalkToken_ = esConsumes<CSCDBCrosstalk, CSCDBCrosstalkRcd>();
0064
0065 crateToken_ = esConsumes<CSCCrateMap, CSCCrateMapRcd>();
0066
0067
0068 makeOccupancyPlots = pset.getUntrackedParameter<bool>("makeOccupancyPlots", true);
0069 makeTriggerPlots = pset.getUntrackedParameter<bool>("makeTriggerPlots", false);
0070 makeStripPlots = pset.getUntrackedParameter<bool>("makeStripPlots", true);
0071 makeWirePlots = pset.getUntrackedParameter<bool>("makeWirePlots", true);
0072 makeRecHitPlots = pset.getUntrackedParameter<bool>("makeRecHitPlots", true);
0073 makeSimHitPlots = pset.getUntrackedParameter<bool>("makeSimHitPlots", true);
0074 makeSegmentPlots = pset.getUntrackedParameter<bool>("makeSegmentPlots", true);
0075 makeResolutionPlots = pset.getUntrackedParameter<bool>("makeResolutionPlots", true);
0076 makePedNoisePlots = pset.getUntrackedParameter<bool>("makePedNoisePlots", true);
0077 makeEfficiencyPlots = pset.getUntrackedParameter<bool>("makeEfficiencyPlots", true);
0078 makeGasGainPlots = pset.getUntrackedParameter<bool>("makeGasGainPlots", true);
0079 makeAFEBTimingPlots = pset.getUntrackedParameter<bool>("makeAFEBTimingPlots", true);
0080 makeCompTimingPlots = pset.getUntrackedParameter<bool>("makeCompTimingPlots", true);
0081 makeADCTimingPlots = pset.getUntrackedParameter<bool>("makeADCTimingPlots", true);
0082 makeRHNoisePlots = pset.getUntrackedParameter<bool>("makeRHNoisePlots", false);
0083 makeCalibPlots = pset.getUntrackedParameter<bool>("makeCalibPlots", false);
0084 makeStandalonePlots = pset.getUntrackedParameter<bool>("makeStandalonePlots", false);
0085 makeTimeMonitorPlots = pset.getUntrackedParameter<bool>("makeTimeMonitorPlots", false);
0086 makeHLTPlots = pset.getUntrackedParameter<bool>("makeHLTPlots", false);
0087
0088
0089 nEventsAnalyzed = 0;
0090 rhTreeCount = 0;
0091 segTreeCount = 0;
0092 firstEvent = true;
0093
0094
0095 theFile = new TFile(rootFileName.c_str(), "RECREATE");
0096 theFile->cd();
0097
0098
0099 histos = new CSCValHists();
0100
0101
0102 hOWires = new TH2I("hOWires", "Wire Digi Occupancy", 36, 0.5, 36.5, 20, 0.5, 20.5);
0103 hOStrips = new TH2I("hOStrips", "Strip Digi Occupancy", 36, 0.5, 36.5, 20, 0.5, 20.5);
0104 hORecHits = new TH2I("hORecHits", "RecHit Occupancy", 36, 0.5, 36.5, 20, 0.5, 20.5);
0105 hOSegments = new TH2I("hOSegments", "Segments Occupancy", 36, 0.5, 36.5, 20, 0.5, 20.5);
0106
0107
0108 hSSTE = new TH1F("hSSTE", "hSSTE", 40, 0, 40);
0109 hRHSTE = new TH1F("hRHSTE", "hRHSTE", 40, 0, 40);
0110 hSEff = new TH1F("hSEff", "Segment Efficiency", 20, 0.5, 20.5);
0111 hRHEff = new TH1F("hRHEff", "recHit Efficiency", 20, 0.5, 20.5);
0112
0113 const int nChambers = 36;
0114 const int nTypes = 20;
0115 float nCH_min = 0.5;
0116 float nCh_max = float(nChambers) + 0.5;
0117 float nT_min = 0.5;
0118 float nT_max = float(nTypes) + 0.5;
0119
0120 hSSTE2 = new TH2F("hSSTE2", "hSSTE2", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0121 hRHSTE2 = new TH2F("hRHSTE2", "hRHSTE2", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0122 hStripSTE2 = new TH2F("hStripSTE2", "hStripSTE2", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0123 hWireSTE2 = new TH2F("hWireSTE2", "hWireSTE2", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0124
0125 hEffDenominator = new TH2F("hEffDenominator", "hEffDenominator", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0126 hSEff2 = new TH2F("hSEff2", "Segment Efficiency 2D", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0127 hRHEff2 = new TH2F("hRHEff2", "recHit Efficiency 2D", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0128
0129 hStripEff2 = new TH2F("hStripEff2", "strip Efficiency 2D", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0130 hWireEff2 = new TH2F("hWireEff2", "wire Efficiency 2D", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0131
0132 hSensitiveAreaEvt =
0133 new TH2F("hSensitiveAreaEvt", "events in sensitive area", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0134
0135 hSSTETight = new TH1F("hSSTETight", "hSSTE Tight", 40, 0, 40);
0136 hRHSTETight = new TH1F("hRHSTETight", "hRHSTE Tight", 40, 0, 40);
0137
0138 hSSTE2Tight = new TH2F("hSSTE2Tight", "hSSTE2 Tight", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0139 hRHSTE2Tight = new TH2F("hRHSTE2Tight", "hRHSTE2 Tight", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0140 hStripSTE2Tight =
0141 new TH2F("hStripSTE2Tight", "hStripSTE2 Tight", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0142 hWireSTE2Tight = new TH2F("hWireSTE2Tight", "hWireSTE2 Tight", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0143 hEffDenominatorTight =
0144 new TH2F("hEffDenominatorTight", "hEffDenominator Tight", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
0145
0146
0147 if (writeTreeToFile)
0148 histos->setupTrees();
0149
0150 geomToken_ = esConsumes<CSCGeometry, MuonGeometryRecord>();
0151 }
0152
0153
0154
0155
0156 CSCValidation::~CSCValidation() {
0157
0158 histoEfficiency(hRHSTE, hRHEff);
0159 histoEfficiency(hSSTE, hSEff);
0160 hSEff2->Divide(hSSTE2, hEffDenominator, 1., 1., "B");
0161 hRHEff2->Divide(hRHSTE2, hEffDenominator, 1., 1., "B");
0162 hStripEff2->Divide(hStripSTE2, hEffDenominator, 1., 1., "B");
0163 hWireEff2->Divide(hWireSTE2, hEffDenominator, 1., 1., "B");
0164
0165 histos->insertPlot(hRHSTETight, "hRHSTETight", "Efficiency");
0166 histos->insertPlot(hSSTETight, "hSSTETight", "Efficiency");
0167 histos->insertPlot(hStripSTE2Tight, "hStripSTE2Tight", "Efficiency");
0168 histos->insertPlot(hWireSTE2Tight, "hWireSTE2Tight", "Efficiency");
0169 histos->insertPlot(hRHSTE2Tight, "hRHSTE2Tight", "Efficiency");
0170 histos->insertPlot(hSSTE2Tight, "hSSTE2Tight", "Efficiency");
0171 histos->insertPlot(hEffDenominatorTight, "hEffDenominatorTight", "Efficiency");
0172
0173 histos->insertPlot(hRHSTE, "hRHSTE", "Efficiency");
0174 histos->insertPlot(hSSTE, "hSSTE", "Efficiency");
0175 histos->insertPlot(hSSTE2, "hSSTE2", "Efficiency");
0176 histos->insertPlot(hEffDenominator, "hEffDenominator", "Efficiency");
0177 histos->insertPlot(hRHSTE2, "hRHSTE2", "Efficiency");
0178 histos->insertPlot(hStripSTE2, "hStripSTE2", "Efficiency");
0179 histos->insertPlot(hWireSTE2, "hWireSTE2", "Efficiency");
0180
0181
0182 histos->insertPlot(hSEff, "hSEff", "Efficiency");
0183 histos->insertPlot(hRHEff, "hRHEff", "Efficiency");
0184
0185 histos->insertPlot(hSEff2, "hSEff2", "Efficiency");
0186 histos->insertPlot(hRHEff2, "hRHEff2", "Efficiency");
0187 histos->insertPlot(hStripEff2, "hStripff2", "Efficiency");
0188 histos->insertPlot(hWireEff2, "hWireff2", "Efficiency");
0189
0190 histos->insertPlot(hSensitiveAreaEvt, "", "Efficiency");
0191
0192
0193 histos->insertPlot(hOWires, "hOWires", "Digis");
0194 histos->insertPlot(hOStrips, "hOStrips", "Digis");
0195 histos->insertPlot(hORecHits, "hORecHits", "recHits");
0196 histos->insertPlot(hOSegments, "hOSegments", "Segments");
0197
0198
0199 histos->writeHists(theFile);
0200 if (writeTreeToFile)
0201 histos->writeTrees(theFile);
0202 theFile->Close();
0203 }
0204
0205
0206
0207
0208 void CSCValidation::analyze(edm::Event const& event, edm::EventSetup const& eventSetup) {
0209
0210 nEventsAnalyzed++;
0211
0212
0213
0214
0215
0216 edm::Handle<CSCWireDigiCollection> wires;
0217 edm::Handle<CSCStripDigiCollection> strips;
0218 edm::Handle<CSCComparatorDigiCollection> compars;
0219 edm::Handle<CSCALCTDigiCollection> alcts;
0220 edm::Handle<CSCCLCTDigiCollection> clcts;
0221 edm::Handle<CSCCorrelatedLCTDigiCollection> correlatedlcts;
0222 if (useDigis) {
0223 event.getByToken(sd_token, strips);
0224 event.getByToken(wd_token, wires);
0225 event.getByToken(cd_token, compars);
0226 event.getByToken(al_token, alcts);
0227 event.getByToken(cl_token, clcts);
0228 event.getByToken(co_token, correlatedlcts);
0229 }
0230
0231
0232 edm::ESHandle<CSCGeometry> cscGeom = eventSetup.getHandle(geomToken_);
0233
0234
0235 edm::Handle<CSCRecHit2DCollection> recHits;
0236 event.getByToken(rh_token, recHits);
0237
0238
0239
0240
0241
0242
0243
0244 edm::Handle<PSimHitContainer> simHits;
0245 if (isSimulation)
0246 event.getByToken(sh_token, simHits);
0247
0248
0249 edm::Handle<CSCSegmentCollection> cscSegments;
0250 event.getByToken(se_token, cscSegments);
0251
0252
0253 edm::Handle<L1MuGMTReadoutCollection> pCollection;
0254 if (makeTriggerPlots || useTriggerFilter || (useDigis && makeTimeMonitorPlots)) {
0255 event.getByToken(l1_token, pCollection);
0256 }
0257 edm::Handle<TriggerResults> hlt;
0258 if (makeHLTPlots) {
0259 event.getByToken(tr_token, hlt);
0260 }
0261
0262
0263 edm::Handle<reco::TrackCollection> saMuons;
0264 if (makeStandalonePlots || useQualityFilter) {
0265 event.getByToken(sa_token, saMuons);
0266 }
0267
0268
0269
0270
0271
0272
0273
0274 if (nEventsAnalyzed == 1 && makeCalibPlots)
0275 doCalibrations(eventSetup);
0276
0277
0278 bool CSCL1A = false;
0279 if (makeTriggerPlots || useTriggerFilter)
0280 CSCL1A = doTrigger(pCollection);
0281 if (!useTriggerFilter)
0282 CSCL1A = true;
0283
0284 cleanEvent = false;
0285 if (makeStandalonePlots || useQualityFilter)
0286 cleanEvent = filterEvents(recHits, cscSegments, saMuons);
0287 if (!useQualityFilter)
0288 cleanEvent = true;
0289
0290
0291
0292 if (makeOccupancyPlots && CSCL1A)
0293 doOccupancies(strips, wires, recHits, cscSegments);
0294
0295 if (makeHLTPlots)
0296 doHLT(hlt);
0297
0298 if (cleanEvent && CSCL1A) {
0299
0300 if (makeStripPlots && useDigis)
0301 doStripDigis(strips);
0302
0303
0304 if (makeWirePlots && useDigis)
0305 doWireDigis(wires);
0306
0307
0308 if (makeRecHitPlots)
0309 doRecHits(recHits, cscGeom);
0310
0311
0312 if (isSimulation && makeSimHitPlots)
0313 doSimHits(recHits, simHits);
0314
0315
0316 if (makeSegmentPlots)
0317 doSegments(cscSegments, cscGeom);
0318
0319
0320 if (makeResolutionPlots)
0321 doResolution(cscSegments, cscGeom);
0322
0323
0324 if (makePedNoisePlots && useDigis)
0325 doPedestalNoise(strips);
0326
0327
0328 if (makeEfficiencyPlots)
0329 doEfficiencies(wires, strips, recHits, cscSegments, cscGeom);
0330
0331
0332 if (makeGasGainPlots && useDigis)
0333 doGasGain(*wires, *strips, *recHits);
0334
0335
0336 if (makeAFEBTimingPlots && useDigis)
0337 doAFEBTiming(*wires);
0338
0339
0340 if (makeCompTimingPlots && useDigis)
0341 doCompTiming(*compars);
0342
0343
0344 if (makeADCTimingPlots)
0345 doADCTiming(*recHits);
0346
0347
0348 if (makeRHNoisePlots && useDigis)
0349 doNoiseHits(recHits, cscSegments, cscGeom, strips);
0350
0351
0352 if (makeStandalonePlots)
0353 doStandalone(saMuons);
0354
0355
0356 if (makeTimeMonitorPlots)
0357 doTimeMonitoring(recHits, cscSegments, alcts, clcts, correlatedlcts, pCollection, cscGeom, eventSetup, event);
0358
0359 firstEvent = false;
0360 }
0361 }
0362
0363
0364
0365
0366
0367
0368
0369 bool CSCValidation::filterEvents(edm::Handle<CSCRecHit2DCollection> recHits,
0370 edm::Handle<CSCSegmentCollection> cscSegments,
0371 edm::Handle<reco::TrackCollection> saMuons) {
0372
0373
0374 if (recHits->size() < 4 || recHits->size() > 100)
0375 return false;
0376 if (cscSegments->size() < 1 || cscSegments->size() > 15)
0377 return false;
0378 return true;
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410
0411
0412
0413
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442 }
0443
0444
0445
0446
0447
0448
0449
0450 void CSCValidation::doOccupancies(edm::Handle<CSCStripDigiCollection> strips,
0451 edm::Handle<CSCWireDigiCollection> wires,
0452 edm::Handle<CSCRecHit2DCollection> recHits,
0453 edm::Handle<CSCSegmentCollection> cscSegments) {
0454 bool wireo[2][4][4][36];
0455 bool stripo[2][4][4][36];
0456 bool rechito[2][4][4][36];
0457 bool segmento[2][4][4][36];
0458
0459 bool hasWires = false;
0460 bool hasStrips = false;
0461 bool hasRecHits = false;
0462 bool hasSegments = false;
0463
0464 for (int e = 0; e < 2; e++) {
0465 for (int s = 0; s < 4; s++) {
0466 for (int r = 0; r < 4; r++) {
0467 for (int c = 0; c < 36; c++) {
0468 wireo[e][s][r][c] = false;
0469 stripo[e][s][r][c] = false;
0470 rechito[e][s][r][c] = false;
0471 segmento[e][s][r][c] = false;
0472 }
0473 }
0474 }
0475 }
0476
0477 if (useDigis) {
0478
0479 for (CSCWireDigiCollection::DigiRangeIterator wi = wires->begin(); wi != wires->end(); wi++) {
0480 CSCDetId id = (CSCDetId)(*wi).first;
0481 int kEndcap = id.endcap();
0482 int kRing = id.ring();
0483 int kStation = id.station();
0484 int kChamber = id.chamber();
0485 std::vector<CSCWireDigi>::const_iterator wireIt = (*wi).second.first;
0486 std::vector<CSCWireDigi>::const_iterator lastWire = (*wi).second.second;
0487 for (; wireIt != lastWire; ++wireIt) {
0488 if (!wireo[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1]) {
0489 wireo[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1] = true;
0490 hOWires->Fill(kChamber, typeIndex(id));
0491 histos->fill1DHist(
0492 chamberSerial(id), "hOWireSerial", "Wire Occupancy by Chamber Serial", 601, -0.5, 600.5, "Digis");
0493 hasWires = true;
0494 }
0495 }
0496 }
0497
0498
0499 for (CSCStripDigiCollection::DigiRangeIterator si = strips->begin(); si != strips->end(); si++) {
0500 CSCDetId id = (CSCDetId)(*si).first;
0501 int kEndcap = id.endcap();
0502 int kRing = id.ring();
0503 int kStation = id.station();
0504 int kChamber = id.chamber();
0505 std::vector<CSCStripDigi>::const_iterator stripIt = (*si).second.first;
0506 std::vector<CSCStripDigi>::const_iterator lastStrip = (*si).second.second;
0507 for (; stripIt != lastStrip; ++stripIt) {
0508 std::vector<int> myADCVals = stripIt->getADCCounts();
0509 bool thisStripFired = false;
0510 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
0511 float threshold = 13.3;
0512 float diff = 0.;
0513 for (unsigned int iCount = 0; iCount < myADCVals.size(); iCount++) {
0514 diff = (float)myADCVals[iCount] - thisPedestal;
0515 if (diff > threshold) {
0516 thisStripFired = true;
0517 }
0518 }
0519 if (thisStripFired) {
0520 if (!stripo[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1]) {
0521 stripo[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1] = true;
0522 hOStrips->Fill(kChamber, typeIndex(id));
0523 histos->fill1DHist(
0524 chamberSerial(id), "hOStripSerial", "Strip Occupancy by Chamber Serial", 601, -0.5, 600.5, "Digis");
0525 hasStrips = true;
0526 }
0527 }
0528 }
0529 }
0530 }
0531
0532
0533 CSCRecHit2DCollection::const_iterator recIt;
0534 for (recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
0535 CSCDetId idrec = (CSCDetId)(*recIt).cscDetId();
0536 int kEndcap = idrec.endcap();
0537 int kRing = idrec.ring();
0538 int kStation = idrec.station();
0539 int kChamber = idrec.chamber();
0540 if (!rechito[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1]) {
0541 rechito[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1] = true;
0542 histos->fill1DHist(
0543 chamberSerial(idrec), "hORecHitsSerial", "RecHit Occupancy by Chamber Serial", 601, -0.5, 600.5, "recHits");
0544 hORecHits->Fill(kChamber, typeIndex(idrec));
0545 hasRecHits = true;
0546 }
0547 }
0548
0549
0550 for (CSCSegmentCollection::const_iterator segIt = cscSegments->begin(); segIt != cscSegments->end(); segIt++) {
0551 CSCDetId id = (CSCDetId)(*segIt).cscDetId();
0552 int kEndcap = id.endcap();
0553 int kRing = id.ring();
0554 int kStation = id.station();
0555 int kChamber = id.chamber();
0556 if (!segmento[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1]) {
0557 segmento[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1] = true;
0558 histos->fill1DHist(
0559 chamberSerial(id), "hOSegmentsSerial", "Segment Occupancy by Chamber Serial", 601, -0.5, 600.5, "Segments");
0560 hOSegments->Fill(kChamber, typeIndex(id));
0561 hasSegments = true;
0562 }
0563 }
0564
0565
0566 histos->fill1DHist(1, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
0567 if (hasWires)
0568 histos->fill1DHist(3, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
0569 if (hasStrips)
0570 histos->fill1DHist(5, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
0571 if (hasWires && hasStrips)
0572 histos->fill1DHist(7, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
0573 if (hasRecHits)
0574 histos->fill1DHist(9, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
0575 if (hasSegments)
0576 histos->fill1DHist(11, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
0577 if (!cleanEvent)
0578 histos->fill1DHist(13, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
0579 }
0580
0581
0582
0583
0584
0585
0586
0587 bool CSCValidation::doTrigger(edm::Handle<L1MuGMTReadoutCollection> pCollection) {
0588 std::vector<L1MuGMTReadoutRecord> L1Mrec = pCollection->getRecords();
0589 std::vector<L1MuGMTReadoutRecord>::const_iterator igmtrr;
0590
0591 bool csc_l1a = false;
0592 bool dt_l1a = false;
0593 bool rpcf_l1a = false;
0594 bool rpcb_l1a = false;
0595 bool beamHaloTrigger = false;
0596
0597 int myBXNumber = -1000;
0598
0599 for (igmtrr = L1Mrec.begin(); igmtrr != L1Mrec.end(); igmtrr++) {
0600 std::vector<L1MuRegionalCand>::const_iterator iter1;
0601 std::vector<L1MuRegionalCand> rmc;
0602
0603
0604 int icsc = 0;
0605 rmc = igmtrr->getCSCCands();
0606 for (iter1 = rmc.begin(); iter1 != rmc.end(); iter1++) {
0607 if (!(*iter1).empty()) {
0608 icsc++;
0609 int kQuality = (*iter1).quality();
0610 if (kQuality == 1)
0611 beamHaloTrigger = true;
0612 }
0613 }
0614 if (igmtrr->getBxInEvent() == 0 && icsc > 0)
0615 csc_l1a = true;
0616 if (igmtrr->getBxInEvent() == 0) {
0617 myBXNumber = igmtrr->getBxNr();
0618 }
0619
0620
0621 int idt = 0;
0622 rmc = igmtrr->getDTBXCands();
0623 for (iter1 = rmc.begin(); iter1 != rmc.end(); iter1++) {
0624 if (!(*iter1).empty()) {
0625 idt++;
0626 }
0627 }
0628 if (igmtrr->getBxInEvent() == 0 && idt > 0)
0629 dt_l1a = true;
0630
0631
0632 int irpcb = 0;
0633 rmc = igmtrr->getBrlRPCCands();
0634 for (iter1 = rmc.begin(); iter1 != rmc.end(); iter1++) {
0635 if (!(*iter1).empty()) {
0636 irpcb++;
0637 }
0638 }
0639 if (igmtrr->getBxInEvent() == 0 && irpcb > 0)
0640 rpcb_l1a = true;
0641
0642
0643 int irpcf = 0;
0644 rmc = igmtrr->getFwdRPCCands();
0645 for (iter1 = rmc.begin(); iter1 != rmc.end(); iter1++) {
0646 if (!(*iter1).empty()) {
0647 irpcf++;
0648 }
0649 }
0650 if (igmtrr->getBxInEvent() == 0 && irpcf > 0)
0651 rpcf_l1a = true;
0652 }
0653
0654
0655 if (csc_l1a)
0656 histos->fill1DHist(myBXNumber, "vtBXNumber", "BX Number", 4001, -0.5, 4000.5, "Trigger");
0657 if (csc_l1a)
0658 histos->fill1DHist(1, "vtBits", "trigger bits", 11, -0.5, 10.5, "Trigger");
0659 if (dt_l1a)
0660 histos->fill1DHist(2, "vtBits", "trigger bits", 11, -0.5, 10.5, "Trigger");
0661 if (rpcb_l1a)
0662 histos->fill1DHist(3, "vtBits", "trigger bits", 11, -0.5, 10.5, "Trigger");
0663 if (rpcf_l1a)
0664 histos->fill1DHist(4, "vtBits", "trigger bits", 11, -0.5, 10.5, "Trigger");
0665 if (beamHaloTrigger)
0666 histos->fill1DHist(8, "vtBits", "trigger bits", 11, -0.5, 10.5, "Trigger");
0667
0668 if (csc_l1a) {
0669 histos->fill1DHist(1, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
0670 if (dt_l1a)
0671 histos->fill1DHist(2, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
0672 if (rpcb_l1a)
0673 histos->fill1DHist(3, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
0674 if (rpcf_l1a)
0675 histos->fill1DHist(4, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
0676 if (dt_l1a || rpcb_l1a || rpcf_l1a)
0677 histos->fill1DHist(5, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
0678 if (!(dt_l1a || rpcb_l1a || rpcf_l1a))
0679 histos->fill1DHist(6, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
0680 } else {
0681 histos->fill1DHist(1, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
0682 if (dt_l1a)
0683 histos->fill1DHist(2, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
0684 if (rpcb_l1a)
0685 histos->fill1DHist(3, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
0686 if (rpcf_l1a)
0687 histos->fill1DHist(4, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
0688 if (dt_l1a || rpcb_l1a || rpcf_l1a)
0689 histos->fill1DHist(5, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
0690 if (!(dt_l1a || rpcb_l1a || rpcf_l1a))
0691 histos->fill1DHist(6, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
0692 }
0693
0694
0695
0696 if (csc_l1a)
0697 return true;
0698
0699 return false;
0700 }
0701
0702
0703
0704
0705
0706
0707
0708 bool CSCValidation::doHLT(Handle<TriggerResults> hlt) {
0709
0710 int hltSize = hlt->size();
0711 for (int i = 0; i < hltSize; ++i) {
0712 if (hlt->accept(i))
0713 histos->fill1DHist(i, "hltBits", "HLT Trigger Bits", hltSize + 1, -0.5, (float)hltSize + 0.5, "Trigger");
0714 }
0715
0716 return true;
0717 }
0718
0719
0720
0721
0722
0723
0724
0725 void CSCValidation::doCalibrations(const edm::EventSetup& eventSetup) {
0726
0727 if (nEventsAnalyzed == 1) {
0728 LogDebug("Calibrations") << "Loading Calibrations...";
0729
0730
0731 edm::ESHandle<CSCDBGains> hGains = eventSetup.getHandle(gainsToken_);
0732 const CSCDBGains* pGains = hGains.product();
0733
0734 edm::ESHandle<CSCDBCrosstalk> hCrosstalk = eventSetup.getHandle(crosstalkToken_);
0735 const CSCDBCrosstalk* pCrosstalk = hCrosstalk.product();
0736
0737 edm::ESHandle<CSCDBNoiseMatrix> hNoiseMatrix = eventSetup.getHandle(noiseToken_);
0738 const CSCDBNoiseMatrix* pNoiseMatrix = hNoiseMatrix.product();
0739
0740 edm::ESHandle<CSCDBPedestals> hPedestals = eventSetup.getHandle(pedestalsToken_);
0741 const CSCDBPedestals* pPedestals = hPedestals.product();
0742
0743 LogDebug("Calibrations") << "Calibrations Loaded!";
0744
0745 for (int i = 0; i < 400; i++) {
0746 int bin = i + 1;
0747 histos->fillCalibHist(pGains->gains[i].gain_slope, "hCalibGainsS", "Gains Slope", 400, 0, 400, bin, "Calib");
0748 histos->fillCalibHist(
0749 pCrosstalk->crosstalk[i].xtalk_slope_left, "hCalibXtalkSL", "Xtalk Slope Left", 400, 0, 400, bin, "Calib");
0750 histos->fillCalibHist(
0751 pCrosstalk->crosstalk[i].xtalk_slope_right, "hCalibXtalkSR", "Xtalk Slope Right", 400, 0, 400, bin, "Calib");
0752 histos->fillCalibHist(pCrosstalk->crosstalk[i].xtalk_intercept_left,
0753 "hCalibXtalkIL",
0754 "Xtalk Intercept Left",
0755 400,
0756 0,
0757 400,
0758 bin,
0759 "Calib");
0760 histos->fillCalibHist(pCrosstalk->crosstalk[i].xtalk_intercept_right,
0761 "hCalibXtalkIR",
0762 "Xtalk Intercept Right",
0763 400,
0764 0,
0765 400,
0766 bin,
0767 "Calib");
0768 histos->fillCalibHist(pPedestals->pedestals[i].ped, "hCalibPedsP", "Peds", 400, 0, 400, bin, "Calib");
0769 histos->fillCalibHist(pPedestals->pedestals[i].rms, "hCalibPedsR", "Peds RMS", 400, 0, 400, bin, "Calib");
0770 histos->fillCalibHist(
0771 pNoiseMatrix->matrix[i].elem33, "hCalibNoise33", "Noise Matrix 33", 400, 0, 400, bin, "Calib");
0772 histos->fillCalibHist(
0773 pNoiseMatrix->matrix[i].elem34, "hCalibNoise34", "Noise Matrix 34", 400, 0, 400, bin, "Calib");
0774 histos->fillCalibHist(
0775 pNoiseMatrix->matrix[i].elem35, "hCalibNoise35", "Noise Matrix 35", 400, 0, 400, bin, "Calib");
0776 histos->fillCalibHist(
0777 pNoiseMatrix->matrix[i].elem44, "hCalibNoise44", "Noise Matrix 44", 400, 0, 400, bin, "Calib");
0778 histos->fillCalibHist(
0779 pNoiseMatrix->matrix[i].elem45, "hCalibNoise45", "Noise Matrix 45", 400, 0, 400, bin, "Calib");
0780 histos->fillCalibHist(
0781 pNoiseMatrix->matrix[i].elem46, "hCalibNoise46", "Noise Matrix 46", 400, 0, 400, bin, "Calib");
0782 histos->fillCalibHist(
0783 pNoiseMatrix->matrix[i].elem55, "hCalibNoise55", "Noise Matrix 55", 400, 0, 400, bin, "Calib");
0784 histos->fillCalibHist(
0785 pNoiseMatrix->matrix[i].elem56, "hCalibNoise56", "Noise Matrix 56", 400, 0, 400, bin, "Calib");
0786 histos->fillCalibHist(
0787 pNoiseMatrix->matrix[i].elem57, "hCalibNoise57", "Noise Matrix 57", 400, 0, 400, bin, "Calib");
0788 histos->fillCalibHist(
0789 pNoiseMatrix->matrix[i].elem66, "hCalibNoise66", "Noise Matrix 66", 400, 0, 400, bin, "Calib");
0790 histos->fillCalibHist(
0791 pNoiseMatrix->matrix[i].elem67, "hCalibNoise67", "Noise Matrix 67", 400, 0, 400, bin, "Calib");
0792 histos->fillCalibHist(
0793 pNoiseMatrix->matrix[i].elem77, "hCalibNoise77", "Noise Matrix 77", 400, 0, 400, bin, "Calib");
0794 }
0795 }
0796 }
0797
0798
0799
0800
0801
0802
0803
0804 void CSCValidation::doWireDigis(edm::Handle<CSCWireDigiCollection> wires) {
0805 int nWireGroupsTotal = 0;
0806 for (CSCWireDigiCollection::DigiRangeIterator dWDiter = wires->begin(); dWDiter != wires->end(); dWDiter++) {
0807 CSCDetId id = (CSCDetId)(*dWDiter).first;
0808 std::vector<CSCWireDigi>::const_iterator wireIter = (*dWDiter).second.first;
0809 std::vector<CSCWireDigi>::const_iterator lWire = (*dWDiter).second.second;
0810 for (; wireIter != lWire; ++wireIter) {
0811 int myWire = wireIter->getWireGroup();
0812 int myTBin = wireIter->getTimeBin();
0813 nWireGroupsTotal++;
0814 histos->fill1DHistByType(myWire, "hWireWire", "Wiregroup Numbers Fired", id, 113, -0.5, 112.5, "Digis");
0815 histos->fill1DHistByType(myTBin, "hWireTBin", "Wire TimeBin Fired", id, 17, -0.5, 16.5, "Digis");
0816 histos->fillProfile(
0817 chamberSerial(id), myTBin, "hWireTBinProfile", "Wire TimeBin Fired", 601, -0.5, 600.5, -0.5, 16.5, "Digis");
0818 if (detailedAnalysis) {
0819 histos->fill1DHistByLayer(
0820 myWire, "hWireWire", "Wiregroup Numbers Fired", id, 113, -0.5, 112.5, "WireNumberByLayer");
0821 histos->fill1DHistByLayer(myTBin, "hWireTBin", "Wire TimeBin Fired", id, 17, -0.5, 16.5, "WireTimeByLayer");
0822 }
0823 }
0824 }
0825
0826
0827 if (nWireGroupsTotal == 0)
0828 nWireGroupsTotal = -1;
0829
0830 histos->fill1DHist(nWireGroupsTotal, "hWirenGroupsTotal", "Wires Fired Per Event", 251, -0.5, 250.5, "Digis");
0831 }
0832
0833
0834
0835
0836
0837
0838
0839 void CSCValidation::doStripDigis(edm::Handle<CSCStripDigiCollection> strips) {
0840 int nStripsFired = 0;
0841 for (CSCStripDigiCollection::DigiRangeIterator dSDiter = strips->begin(); dSDiter != strips->end(); dSDiter++) {
0842 CSCDetId id = (CSCDetId)(*dSDiter).first;
0843 std::vector<CSCStripDigi>::const_iterator stripIter = (*dSDiter).second.first;
0844 std::vector<CSCStripDigi>::const_iterator lStrip = (*dSDiter).second.second;
0845 for (; stripIter != lStrip; ++stripIter) {
0846 int myStrip = stripIter->getStrip();
0847 std::vector<int> myADCVals = stripIter->getADCCounts();
0848 bool thisStripFired = false;
0849 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
0850 float threshold = 13.3;
0851 float diff = 0.;
0852 for (unsigned int iCount = 0; iCount < myADCVals.size(); iCount++) {
0853 diff = (float)myADCVals[iCount] - thisPedestal;
0854 if (diff > threshold) {
0855 thisStripFired = true;
0856 }
0857 }
0858 if (thisStripFired) {
0859 nStripsFired++;
0860
0861 histos->fill1DHistByType(myStrip, "hStripStrip", "Strip Number", id, 81, -0.5, 80.5, "Digis");
0862 if (detailedAnalysis) {
0863 histos->fill1DHistByLayer(myStrip, "hStripStrip", "Strip Number", id, 81, -0.5, 80.5, "StripNumberByLayer");
0864 }
0865 }
0866 }
0867 }
0868
0869 if (nStripsFired == 0)
0870 nStripsFired = -1;
0871
0872 histos->fill1DHist(nStripsFired, "hStripNFired", "Fired Strips per Event", 351, -0.5, 350.5, "Digis");
0873 }
0874
0875
0876
0877
0878
0879
0880
0881 void CSCValidation::doPedestalNoise(edm::Handle<CSCStripDigiCollection> strips) {
0882 constexpr float threshold = 13.3;
0883 for (CSCStripDigiCollection::DigiRangeIterator dPNiter = strips->begin(); dPNiter != strips->end(); dPNiter++) {
0884 CSCDetId id = (CSCDetId)(*dPNiter).first;
0885 std::vector<CSCStripDigi>::const_iterator pedIt = (*dPNiter).second.first;
0886 std::vector<CSCStripDigi>::const_iterator lStrip = (*dPNiter).second.second;
0887 for (; pedIt != lStrip; ++pedIt) {
0888 int myStrip = pedIt->getStrip();
0889 std::vector<int> myADCVals = pedIt->getADCCounts();
0890 float TotalADC = getSignal(*strips, id, myStrip);
0891 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
0892 float thisSignal =
0893 (1. / 6) * (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
0894 bool thisStripFired = TotalADC > threshold;
0895 if (!thisStripFired) {
0896 float ADC = thisSignal - thisPedestal;
0897 histos->fill1DHist(ADC, "hStripPed", "Pedestal Noise Distribution", 50, -25., 25., "PedestalNoise");
0898 histos->fill1DHistByType(ADC, "hStripPedME", "Pedestal Noise Distribution", id, 50, -25., 25., "PedestalNoise");
0899 histos->fillProfile(chamberSerial(id),
0900 ADC,
0901 "hStripPedMEProfile",
0902 "Wire TimeBin Fired",
0903 601,
0904 -0.5,
0905 600.5,
0906 -25,
0907 25,
0908 "PedestalNoise");
0909 if (detailedAnalysis) {
0910 histos->fill1DHistByLayer(
0911 ADC, "hStripPedME", "Pedestal Noise Distribution", id, 50, -25., 25., "PedestalNoiseByLayer");
0912 }
0913 }
0914 }
0915 }
0916 }
0917
0918
0919
0920
0921
0922
0923
0924 void CSCValidation::doRecHits(edm::Handle<CSCRecHit2DCollection> recHits, edm::ESHandle<CSCGeometry> cscGeom) {
0925
0926 int nRecHits = recHits->size();
0927
0928
0929
0930
0931
0932
0933 CSCRecHit2DCollection::const_iterator dRHIter;
0934 for (dRHIter = recHits->begin(); dRHIter != recHits->end(); dRHIter++) {
0935
0936 CSCDetId idrec = (CSCDetId)(*dRHIter).cscDetId();
0937 int kEndcap = idrec.endcap();
0938 int kRing = idrec.ring();
0939 int kStation = idrec.station();
0940 int kChamber = idrec.chamber();
0941 int kLayer = idrec.layer();
0942
0943
0944 LocalPoint rhitlocal = (*dRHIter).localPosition();
0945 float xreco = rhitlocal.x();
0946 float yreco = rhitlocal.y();
0947 LocalError rerrlocal = (*dRHIter).localPositionError();
0948
0949 float xxerr = rerrlocal.xx();
0950 float yyerr = rerrlocal.yy();
0951 float xyerr = rerrlocal.xy();
0952
0953 float stpos = (*dRHIter).positionWithinStrip();
0954 float sterr = (*dRHIter).errorWithinStrip();
0955
0956
0957 float rHSumQ = 0;
0958 float sumsides = 0.;
0959 int adcsize = dRHIter->nStrips() * dRHIter->nTimeBins();
0960 for (unsigned int i = 0; i < dRHIter->nStrips(); i++) {
0961 for (unsigned int j = 0; j < dRHIter->nTimeBins() - 1; j++) {
0962 rHSumQ += dRHIter->adcs(i, j);
0963 if (i != 1)
0964 sumsides += dRHIter->adcs(i, j);
0965 }
0966 }
0967
0968 float rHratioQ = sumsides / rHSumQ;
0969 if (adcsize != 12)
0970 rHratioQ = -99;
0971
0972
0973 float rHtime = 0;
0974 rHtime = (*dRHIter).tpeak() / 50.;
0975
0976
0977 const CSCLayer* csclayer = cscGeom->layer(idrec);
0978
0979
0980 GlobalPoint rhitglobal = csclayer->toGlobal(rhitlocal);
0981 float grecx = rhitglobal.x();
0982 float grecy = rhitglobal.y();
0983
0984
0985 if (writeTreeToFile && rhTreeCount < 1500000) {
0986 histos->fillRechitTree(xreco, yreco, grecx, grecy, kEndcap, kStation, kRing, kChamber, kLayer);
0987 rhTreeCount++;
0988 }
0989
0990
0991
0992 histos->fill2DHistByStation(
0993 grecx, grecy, "hRHGlobal", "recHit Global Position", idrec, 100, -800., 800., 100, -800., 800., "recHits");
0994 if (kStation == 1 && (kRing == 1 || kRing == 4))
0995 histos->fill1DHistByType(rHSumQ, "hRHSumQ", "Sum 3x3 recHit Charge", idrec, 125, 0, 4000, "recHits");
0996 else
0997 histos->fill1DHistByType(rHSumQ, "hRHSumQ", "Sum 3x3 recHit Charge", idrec, 125, 0, 2000, "recHits");
0998 histos->fill1DHistByType(rHratioQ, "hRHRatioQ", "Charge Ratio (Ql+Qr)/Qt", idrec, 120, -0.1, 1.1, "recHits");
0999 histos->fill1DHistByType(rHtime, "hRHTiming", "recHit Timing", idrec, 200, -10, 10, "recHits");
1000 histos->fill1DHistByType(sqrt(xxerr), "hRHxerr", "RecHit Error on Local X", idrec, 100, -0.1, 2, "recHits");
1001 histos->fill1DHistByType(sqrt(yyerr), "hRHyerr", "RecHit Error on Local Y", idrec, 100, -0.1, 2, "recHits");
1002 histos->fill1DHistByType(xyerr, "hRHxyerr", "Corr. RecHit XY Error", idrec, 100, -1, 2, "recHits");
1003 if (adcsize == 12)
1004 histos->fill1DHistByType(stpos, "hRHstpos", "Reconstructed Position on Strip", idrec, 120, -0.6, 0.6, "recHits");
1005 histos->fill1DHistByType(
1006 sterr, "hRHsterr", "Estimated Error on Strip Measurement", idrec, 120, -0.05, 0.25, "recHits");
1007 histos->fillProfile(
1008 chamberSerial(idrec), rHSumQ, "hRHSumQProfile", "Sum 3x3 recHit Charge", 601, -0.5, 600.5, 0, 4000, "recHits");
1009 histos->fillProfile(
1010 chamberSerial(idrec), rHtime, "hRHTimingProfile", "recHit Timing", 601, -0.5, 600.5, -11, 11, "recHits");
1011 if (detailedAnalysis) {
1012 if (kStation == 1 && (kRing == 1 || kRing == 4))
1013 histos->fill1DHistByLayer(rHSumQ, "hRHSumQ", "Sum 3x3 recHit Charge", idrec, 125, 0, 4000, "RHQByLayer");
1014 else
1015 histos->fill1DHistByLayer(rHSumQ, "hRHSumQ", "Sum 3x3 recHit Charge", idrec, 125, 0, 2000, "RHQByLayer");
1016 histos->fill1DHistByLayer(rHratioQ, "hRHRatioQ", "Charge Ratio (Ql+Qr)/Qt", idrec, 120, -0.1, 1.1, "RHQByLayer");
1017 histos->fill1DHistByLayer(rHtime, "hRHTiming", "recHit Timing", idrec, 200, -10, 10, "RHTimingByLayer");
1018 histos->fill2DHistByLayer(xreco,
1019 yreco,
1020 "hRHLocalXY",
1021 "recHit Local Position",
1022 idrec,
1023 50,
1024 -100.,
1025 100.,
1026 75,
1027 -150.,
1028 150.,
1029 "RHLocalXYByLayer");
1030 histos->fill1DHistByLayer(
1031 sqrt(xxerr), "hRHxerr", "RecHit Error on Local X", idrec, 100, -0.1, 2, "RHErrorsByLayer");
1032 histos->fill1DHistByLayer(
1033 sqrt(yyerr), "hRHyerr", "RecHit Error on Local Y", idrec, 100, -0.1, 2, "RHErrorsByLayer");
1034 histos->fill1DHistByType(
1035 stpos, "hRHstpos", "Reconstructed Position on Strip", idrec, 120, -0.6, 0.6, "RHStripPosByLayer");
1036 histos->fill1DHistByType(
1037 sterr, "hRHsterr", "Estimated Error on Strip Measurement", idrec, 120, -0.05, 0.25, "RHStripPosByLayer");
1038 }
1039
1040 }
1041
1042 if (nRecHits == 0)
1043 nRecHits = -1;
1044
1045 histos->fill1DHist(nRecHits, "hRHnrechits", "recHits per Event (all chambers)", 201, -0.5, 200.5, "recHits");
1046 }
1047
1048
1049
1050
1051
1052
1053
1054 void CSCValidation::doSimHits(edm::Handle<CSCRecHit2DCollection> recHits, edm::Handle<PSimHitContainer> simHits) {
1055 CSCRecHit2DCollection::const_iterator dSHrecIter;
1056 for (dSHrecIter = recHits->begin(); dSHrecIter != recHits->end(); dSHrecIter++) {
1057 CSCDetId idrec = (CSCDetId)(*dSHrecIter).cscDetId();
1058 LocalPoint rhitlocal = (*dSHrecIter).localPosition();
1059 float xreco = rhitlocal.x();
1060 float yreco = rhitlocal.y();
1061 float xError = sqrt((*dSHrecIter).localPositionError().xx());
1062 float yError = sqrt((*dSHrecIter).localPositionError().yy());
1063 float simHitXres = -99;
1064 float simHitYres = -99;
1065 float xPull = -99;
1066 float yPull = -99;
1067 float mindiffX = 99;
1068 float mindiffY = 10;
1069
1070 PSimHitContainer::const_iterator dSHsimIter;
1071 for (dSHsimIter = simHits->begin(); dSHsimIter != simHits->end(); dSHsimIter++) {
1072
1073 CSCDetId sId = (CSCDetId)(*dSHsimIter).detUnitId();
1074
1075
1076 if (sId == idrec && abs((*dSHsimIter).particleType()) == 13) {
1077
1078 LocalPoint sHitlocal = (*dSHsimIter).localPosition();
1079
1080
1081 if ((sHitlocal.x() - xreco) < mindiffX && (sHitlocal.y() - yreco) < mindiffY) {
1082 simHitXres = (sHitlocal.x() - xreco);
1083 simHitYres = (sHitlocal.y() - yreco);
1084 mindiffX = (sHitlocal.x() - xreco);
1085 xPull = simHitXres / xError;
1086 yPull = simHitYres / yError;
1087 }
1088 }
1089 }
1090
1091 histos->fill1DHistByType(
1092 simHitXres, "hSimXResid", "SimHitX - Reconstructed X", idrec, 100, -1.0, 1.0, "Resolution");
1093 histos->fill1DHistByType(
1094 simHitYres, "hSimYResid", "SimHitY - Reconstructed Y", idrec, 100, -5.0, 5.0, "Resolution");
1095 histos->fill1DHistByType(xPull, "hSimXPull", "Local X Pulls", idrec, 100, -5.0, 5.0, "Resolution");
1096 histos->fill1DHistByType(yPull, "hSimYPull", "Local Y Pulls", idrec, 100, -5.0, 5.0, "Resolution");
1097 }
1098 }
1099
1100
1101
1102
1103
1104
1105
1106 void CSCValidation::doSegments(edm::Handle<CSCSegmentCollection> cscSegments, edm::ESHandle<CSCGeometry> cscGeom) {
1107
1108 int nSegments = cscSegments->size();
1109
1110
1111
1112
1113 for (CSCSegmentCollection::const_iterator dSiter = cscSegments->begin(); dSiter != cscSegments->end(); dSiter++) {
1114
1115 CSCDetId id = (CSCDetId)(*dSiter).cscDetId();
1116 int kEndcap = id.endcap();
1117 int kRing = id.ring();
1118 int kStation = id.station();
1119 int kChamber = id.chamber();
1120
1121
1122 float chisq = (*dSiter).chi2();
1123 int nhits = (*dSiter).nRecHits();
1124 int nDOF = 2 * nhits - 4;
1125 double chisqProb = ChiSquaredProbability((double)chisq, nDOF);
1126 LocalPoint localPos = (*dSiter).localPosition();
1127 float segX = localPos.x();
1128 float segY = localPos.y();
1129 LocalVector segDir = (*dSiter).localDirection();
1130 double theta = segDir.theta();
1131
1132
1133 float globX = 0.;
1134 float globY = 0.;
1135 float globTheta = 0.;
1136 float globPhi = 0.;
1137 const CSCChamber* cscchamber = cscGeom->chamber(id);
1138 if (cscchamber) {
1139 GlobalPoint globalPosition = cscchamber->toGlobal(localPos);
1140 globX = globalPosition.x();
1141 globY = globalPosition.y();
1142 GlobalVector globalDirection = cscchamber->toGlobal(segDir);
1143 globTheta = globalDirection.theta();
1144 globPhi = globalDirection.phi();
1145 }
1146
1147
1148 if (writeTreeToFile && segTreeCount < 1500000) {
1149 histos->fillSegmentTree(segX, segY, globX, globY, kEndcap, kStation, kRing, kChamber);
1150 segTreeCount++;
1151 }
1152
1153
1154 histos->fill2DHistByStation(globX,
1155 globY,
1156 "hSGlobal",
1157 "Segment Global Positions;global x (cm)",
1158 id,
1159 100,
1160 -800.,
1161 800.,
1162 100,
1163 -800.,
1164 800.,
1165 "Segments");
1166 histos->fill1DHistByType(nhits, "hSnHits", "N hits on Segments", id, 8, -0.5, 7.5, "Segments");
1167 histos->fill1DHistByType(theta, "hSTheta", "local theta segments", id, 128, -3.2, 3.2, "Segments");
1168 histos->fill1DHistByType((chisq / nDOF), "hSChiSq", "segments chi-squared/ndof", id, 110, -0.05, 10.5, "Segments");
1169 histos->fill1DHistByType(
1170 chisqProb, "hSChiSqProb", "segments chi-squared probability", id, 110, -0.05, 1.05, "Segments");
1171 histos->fill1DHist(globTheta, "hSGlobalTheta", "segment global theta", 128, 0, 3.2, "Segments");
1172 histos->fill1DHist(globPhi, "hSGlobalPhi", "segment global phi", 128, -3.2, 3.2, "Segments");
1173 histos->fillProfile(
1174 chamberSerial(id), nhits, "hSnHitsProfile", "N hits on Segments", 601, -0.5, 600.5, -0.5, 7.5, "Segments");
1175 if (detailedAnalysis) {
1176 histos->fill1DHistByChamber(nhits, "hSnHits", "N hits on Segments", id, 8, -0.5, 7.5, "HitsOnSegmentByChamber");
1177 histos->fill1DHistByChamber(theta, "hSTheta", "local theta segments", id, 128, -3.2, 3.2, "DetailedSegments");
1178 histos->fill1DHistByChamber(
1179 (chisq / nDOF), "hSChiSq", "segments chi-squared/ndof", id, 110, -0.05, 10.5, "SegChi2ByChamber");
1180 histos->fill1DHistByChamber(
1181 chisqProb, "hSChiSqProb", "segments chi-squared probability", id, 110, -0.05, 1.05, "SegChi2ByChamber");
1182 }
1183
1184 }
1185
1186 if (nSegments == 0)
1187 nSegments = -1;
1188
1189 histos->fill1DHist(nSegments, "hSnSegments", "Segments per Event", 31, -0.5, 30.5, "Segments");
1190 }
1191
1192
1193
1194
1195
1196
1197
1198 void CSCValidation::doResolution(edm::Handle<CSCSegmentCollection> cscSegments, edm::ESHandle<CSCGeometry> cscGeom) {
1199 for (CSCSegmentCollection::const_iterator dSiter = cscSegments->begin(); dSiter != cscSegments->end(); dSiter++) {
1200 CSCDetId id = (CSCDetId)(*dSiter).cscDetId();
1201
1202
1203
1204 std::vector<CSCRecHit2D> theseRecHits = (*dSiter).specificRecHits();
1205 int nRH = (*dSiter).nRecHits();
1206 CLHEP::HepMatrix sp(6, 1);
1207 CLHEP::HepMatrix se(6, 1);
1208 for (std::vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
1209 CSCDetId idRH = (CSCDetId)(*iRH).cscDetId();
1210 int kRing = idRH.ring();
1211 int kStation = idRH.station();
1212 int kLayer = idRH.layer();
1213
1214
1215 int centerid = iRH->nStrips() / 2;
1216 int centerStrip = iRH->channels(centerid);
1217
1218
1219 if (nRH == 6) {
1220 float stpos = (*iRH).positionWithinStrip();
1221 se(kLayer, 1) = (*iRH).errorWithinStrip();
1222
1223 if (kStation == 1 && (kRing == 1 || kRing == 4))
1224 sp(kLayer, 1) = stpos + centerStrip;
1225 else {
1226 if (kLayer == 1 || kLayer == 3 || kLayer == 5)
1227 sp(kLayer, 1) = stpos + centerStrip;
1228 if (kLayer == 2 || kLayer == 4 || kLayer == 6)
1229 sp(kLayer, 1) = stpos - 0.5 + centerStrip;
1230 }
1231 }
1232 }
1233
1234 float residual = -99;
1235 float pull = -99;
1236
1237 if (nRH == 6) {
1238 float expected = fitX(sp, se);
1239 residual = expected - sp(3, 1);
1240 pull = residual / se(3, 1);
1241 }
1242
1243
1244 histos->fill1DHistByType(
1245 residual, "hSResid", "Fitted Position on Strip - Reconstructed for Layer 3", id, 100, -0.5, 0.5, "Resolution");
1246 histos->fill1DHistByType(pull, "hSStripPosPull", "Strip Measurement Pulls", id, 100, -5.0, 5.0, "Resolution");
1247 histos->fillProfile(chamberSerial(id),
1248 residual,
1249 "hSResidProfile",
1250 "Fitted Position on Strip - Reconstructed for Layer 3",
1251 601,
1252 -0.5,
1253 600.5,
1254 -0.5,
1255 0.5,
1256 "Resolution");
1257 if (detailedAnalysis) {
1258 histos->fill1DHistByChamber(residual,
1259 "hSResid",
1260 "Fitted Position on Strip - Reconstructed for Layer 3",
1261 id,
1262 100,
1263 -0.5,
1264 0.5,
1265 "DetailedResolution");
1266 histos->fill1DHistByChamber(pull, "hSStripPosPull", "Strip Measurement Pulls", id, 100, -5.0, 5.0, "Resolution");
1267 }
1268 }
1269 }
1270
1271
1272
1273
1274
1275
1276
1277 void CSCValidation::doStandalone(Handle<reco::TrackCollection> saMuons) {
1278 int nSAMuons = saMuons->size();
1279 histos->fill1DHist(nSAMuons, "trNSAMuons", "N Standalone Muons per Event", 6, -0.5, 5.5, "STAMuons");
1280
1281 for (reco::TrackCollection::const_iterator muon = saMuons->begin(); muon != saMuons->end(); ++muon) {
1282 float preco = muon->p();
1283 float ptreco = muon->pt();
1284 int n = muon->recHitsSize();
1285 float chi2 = muon->chi2();
1286 float normchi2 = muon->normalizedChi2();
1287
1288
1289 int nDTHits = 0;
1290 int nCSCHits = 0;
1291 int nCSCHitsp = 0;
1292 int nCSCHitsm = 0;
1293 int nRPCHits = 0;
1294 int nRPCHitsp = 0;
1295 int nRPCHitsm = 0;
1296 int np = 0;
1297 int nm = 0;
1298 std::vector<CSCDetId> staChambers;
1299 for (trackingRecHit_iterator hit = muon->recHitsBegin(); hit != muon->recHitsEnd(); ++hit) {
1300 const DetId detId((*hit)->geographicalId());
1301 if (detId.det() == DetId::Muon) {
1302 if (detId.subdetId() == MuonSubdetId::RPC) {
1303 RPCDetId rpcId(detId.rawId());
1304 nRPCHits++;
1305 if (rpcId.region() == 1) {
1306 nRPCHitsp++;
1307 np++;
1308 }
1309 if (rpcId.region() == -1) {
1310 nRPCHitsm++;
1311 nm++;
1312 }
1313 }
1314 if (detId.subdetId() == MuonSubdetId::DT) {
1315 nDTHits++;
1316 } else if (detId.subdetId() == MuonSubdetId::CSC) {
1317 CSCDetId cscId(detId.rawId());
1318 staChambers.push_back(detId.rawId());
1319 nCSCHits++;
1320 if (cscId.endcap() == 1) {
1321 nCSCHitsp++;
1322 np++;
1323 }
1324 if (cscId.endcap() == 2) {
1325 nCSCHitsm++;
1326 nm++;
1327 }
1328 }
1329 }
1330 }
1331
1332 GlobalPoint innerPnt(muon->innerPosition().x(), muon->innerPosition().y(), muon->innerPosition().z());
1333 GlobalPoint outerPnt(muon->outerPosition().x(), muon->outerPosition().y(), muon->outerPosition().z());
1334 GlobalVector innerKin(muon->innerMomentum().x(), muon->innerMomentum().y(), muon->innerMomentum().z());
1335 GlobalVector outerKin(muon->outerMomentum().x(), muon->outerMomentum().y(), muon->outerMomentum().z());
1336 GlobalVector deltaPnt = innerPnt - outerPnt;
1337 double crudeLength = deltaPnt.mag();
1338 double deltaPhi = innerPnt.phi() - outerPnt.phi();
1339 double innerGlobalPolarAngle = innerKin.theta();
1340 double outerGlobalPolarAngle = outerKin.theta();
1341
1342
1343 histos->fill1DHist(n, "trN", "N hits on a STA Muon Track", 51, -0.5, 50.5, "STAMuons");
1344 if (np != 0)
1345 histos->fill1DHist(np, "trNp", "N hits on a STA Muon Track (plus endcap)", 51, -0.5, 50.5, "STAMuons");
1346 if (nm != 0)
1347 histos->fill1DHist(nm, "trNm", "N hits on a STA Muon Track (minus endcap)", 51, -0.5, 50.5, "STAMuons");
1348 histos->fill1DHist(nDTHits, "trNDT", "N DT hits on a STA Muon Track", 51, -0.5, 50.5, "STAMuons");
1349 histos->fill1DHist(nCSCHits, "trNCSC", "N CSC hits on a STA Muon Track", 51, -0.5, 50.5, "STAMuons");
1350 if (nCSCHitsp != 0)
1351 histos->fill1DHist(nCSCHitsp, "trNCSCp", "N CSC hits on a STA Muon Track (+ endcap)", 51, -0.5, 50.5, "STAMuons");
1352 if (nCSCHitsm != 0)
1353 histos->fill1DHist(nCSCHitsm, "trNCSCm", "N CSC hits on a STA Muon Track (- endcap)", 51, -0.5, 50.5, "STAMuons");
1354 histos->fill1DHist(nRPCHits, "trNRPC", "N RPC hits on a STA Muon Track", 51, -0.5, 50.5, "STAMuons");
1355 if (nRPCHitsp != 0)
1356 histos->fill1DHist(nRPCHitsp, "trNRPCp", "N RPC hits on a STA Muon Track (+ endcap)", 51, -0.5, 50.5, "STAMuons");
1357 if (nRPCHitsm != 0)
1358 histos->fill1DHist(nRPCHitsm, "trNRPCm", "N RPC hits on a STA Muon Track (- endcap)", 51, -0.5, 50.5, "STAMuons");
1359 histos->fill1DHist(preco, "trP", "STA Muon Momentum", 100, 0, 300, "STAMuons");
1360 histos->fill1DHist(ptreco, "trPT", "STA Muon pT", 100, 0, 40, "STAMuons");
1361 histos->fill1DHist(chi2, "trChi2", "STA Muon Chi2", 100, 0, 200, "STAMuons");
1362 histos->fill1DHist(normchi2, "trNormChi2", "STA Muon Normalized Chi2", 100, 0, 10, "STAMuons");
1363 histos->fill1DHist(crudeLength, "trLength", "Straight Line Length of STA Muon", 120, 0., 2400., "STAMuons");
1364 histos->fill1DHist(
1365 deltaPhi, "trDeltaPhi", "Delta-Phi Between Inner and Outer STA Muon Pos.", 100, -0.5, 0.5, "STAMuons");
1366 histos->fill1DHist(
1367 innerGlobalPolarAngle, "trInnerPolar", "Polar Angle of Inner P Vector (STA muons)", 128, 0, 3.2, "STAMuons");
1368 histos->fill1DHist(
1369 outerGlobalPolarAngle, "trOuterPolar", "Polar Angle of Outer P Vector (STA muons)", 128, 0, 3.2, "STAMuons");
1370 histos->fill1DHist(innerPnt.phi(), "trInnerPhi", "Phi of Inner Position (STA muons)", 256, -3.2, 3.2, "STAMuons");
1371 histos->fill1DHist(outerPnt.phi(), "trOuterPhi", "Phi of Outer Position (STA muons)", 256, -3.2, 3.2, "STAMuons");
1372 }
1373 }
1374
1375
1376
1377
1378
1379 int CSCValidation::chamberSerial(CSCDetId id) {
1380 int st = id.station();
1381 int ri = id.ring();
1382 int ch = id.chamber();
1383 int ec = id.endcap();
1384 int kSerial = ch;
1385 if (st == 1 && ri == 1)
1386 kSerial = ch;
1387 if (st == 1 && ri == 2)
1388 kSerial = ch + 36;
1389 if (st == 1 && ri == 3)
1390 kSerial = ch + 72;
1391 if (st == 1 && ri == 4)
1392 kSerial = ch;
1393 if (st == 2 && ri == 1)
1394 kSerial = ch + 108;
1395 if (st == 2 && ri == 2)
1396 kSerial = ch + 126;
1397 if (st == 3 && ri == 1)
1398 kSerial = ch + 162;
1399 if (st == 3 && ri == 2)
1400 kSerial = ch + 180;
1401 if (st == 4 && ri == 1)
1402 kSerial = ch + 216;
1403 if (st == 4 && ri == 2)
1404 kSerial = ch + 234;
1405 if (ec == 2)
1406 kSerial = kSerial + 300;
1407 return kSerial;
1408 }
1409
1410
1411
1412
1413
1414 int CSCValidation::ringSerial(CSCDetId id) {
1415 int st = id.station();
1416 int ri = id.ring();
1417 int ec = id.endcap();
1418 int kSerial = 0;
1419 if (st == 1 && ri == 1)
1420 kSerial = ri;
1421 if (st == 1 && ri == 2)
1422 kSerial = ri;
1423 if (st == 1 && ri == 3)
1424 kSerial = ri;
1425 if (st == 1 && ri == 4)
1426 kSerial = 1;
1427 if (st == 2)
1428 kSerial = ri + 3;
1429 if (st == 3)
1430 kSerial = ri + 5;
1431 if (st == 4)
1432 kSerial = ri + 7;
1433 if (ec == 2)
1434 kSerial = kSerial * (-1);
1435 return kSerial;
1436 }
1437
1438
1439
1440
1441
1442
1443 float CSCValidation::fitX(const CLHEP::HepMatrix& points, const CLHEP::HepMatrix& errors) {
1444 float S = 0;
1445 float Sx = 0;
1446 float Sy = 0;
1447 float Sxx = 0;
1448 float Sxy = 0;
1449 float sigma2 = 0;
1450
1451 for (int i = 1; i < 7; i++) {
1452 if (i != 3) {
1453 sigma2 = errors(i, 1) * errors(i, 1);
1454 S = S + (1 / sigma2);
1455 Sy = Sy + (points(i, 1) / sigma2);
1456 Sx = Sx + ((i) / sigma2);
1457 Sxx = Sxx + (i * i) / sigma2;
1458 Sxy = Sxy + (((i)*points(i, 1)) / sigma2);
1459 }
1460 }
1461
1462 float delta = S * Sxx - Sx * Sx;
1463 float intercept = (Sxx * Sy - Sx * Sxy) / delta;
1464 float slope = (S * Sxy - Sx * Sy) / delta;
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475 return (intercept + slope * 3);
1476 }
1477
1478
1479
1480
1481
1482
1483 void CSCValidation::doEfficiencies(edm::Handle<CSCWireDigiCollection> wires,
1484 edm::Handle<CSCStripDigiCollection> strips,
1485 edm::Handle<CSCRecHit2DCollection> recHits,
1486 edm::Handle<CSCSegmentCollection> cscSegments,
1487 edm::ESHandle<CSCGeometry> cscGeom) {
1488 bool allWires[2][4][4][36][6];
1489 bool allStrips[2][4][4][36][6];
1490 bool AllRecHits[2][4][4][36][6];
1491 bool AllSegments[2][4][4][36];
1492
1493
1494 for (int iE = 0; iE < 2; iE++) {
1495 for (int iS = 0; iS < 4; iS++) {
1496 for (int iR = 0; iR < 4; iR++) {
1497 for (int iC = 0; iC < 36; iC++) {
1498 AllSegments[iE][iS][iR][iC] = false;
1499
1500 for (int iL = 0; iL < 6; iL++) {
1501 allWires[iE][iS][iR][iC][iL] = false;
1502 allStrips[iE][iS][iR][iC][iL] = false;
1503 AllRecHits[iE][iS][iR][iC][iL] = false;
1504 }
1505 }
1506 }
1507 }
1508 }
1509
1510 if (useDigis) {
1511
1512 for (CSCWireDigiCollection::DigiRangeIterator dWDiter = wires->begin(); dWDiter != wires->end(); dWDiter++) {
1513 CSCDetId idrec = (CSCDetId)(*dWDiter).first;
1514 std::vector<CSCWireDigi>::const_iterator wireIter = (*dWDiter).second.first;
1515 std::vector<CSCWireDigi>::const_iterator lWire = (*dWDiter).second.second;
1516 for (; wireIter != lWire; ++wireIter) {
1517 allWires[idrec.endcap() - 1][idrec.station() - 1][idrec.ring() - 1][idrec.chamber() - 1][idrec.layer() - 1] =
1518 true;
1519 break;
1520 }
1521 }
1522
1523
1524 for (CSCStripDigiCollection::DigiRangeIterator dSDiter = strips->begin(); dSDiter != strips->end(); dSDiter++) {
1525 CSCDetId idrec = (CSCDetId)(*dSDiter).first;
1526 std::vector<CSCStripDigi>::const_iterator stripIter = (*dSDiter).second.first;
1527 std::vector<CSCStripDigi>::const_iterator lStrip = (*dSDiter).second.second;
1528 for (; stripIter != lStrip; ++stripIter) {
1529 std::vector<int> myADCVals = stripIter->getADCCounts();
1530 bool thisStripFired = false;
1531 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
1532 float threshold = 13.3;
1533 float diff = 0.;
1534 for (unsigned int iCount = 0; iCount < myADCVals.size(); iCount++) {
1535 diff = (float)myADCVals[iCount] - thisPedestal;
1536 if (diff > threshold) {
1537 thisStripFired = true;
1538 break;
1539 }
1540 }
1541 if (thisStripFired) {
1542 allStrips[idrec.endcap() - 1][idrec.station() - 1][idrec.ring() - 1][idrec.chamber() - 1][idrec.layer() - 1] =
1543 true;
1544 break;
1545 }
1546 }
1547 }
1548 }
1549
1550
1551 for (CSCRecHit2DCollection::const_iterator recEffIt = recHits->begin(); recEffIt != recHits->end(); recEffIt++) {
1552
1553 CSCDetId idrec = (CSCDetId)(*recEffIt).cscDetId();
1554 AllRecHits[idrec.endcap() - 1][idrec.station() - 1][idrec.ring() - 1][idrec.chamber() - 1][idrec.layer() - 1] =
1555 true;
1556 }
1557
1558 std::vector<unsigned int> seg_ME2(2, 0);
1559 std::vector<unsigned int> seg_ME3(2, 0);
1560 std::vector<std::pair<CSCDetId, CSCSegment> > theSegments(4);
1561
1562 for (CSCSegmentCollection::const_iterator segEffIt = cscSegments->begin(); segEffIt != cscSegments->end();
1563 segEffIt++) {
1564 CSCDetId idseg = (CSCDetId)(*segEffIt).cscDetId();
1565
1566
1567
1568 AllSegments[idseg.endcap() - 1][idseg.station() - 1][idseg.ring() - 1][idseg.chamber() - 1] = true;
1569
1570
1571
1572 if (2 == idseg.station() || 3 == idseg.station()) {
1573 unsigned int seg_tmp;
1574 if (2 == idseg.station()) {
1575 ++seg_ME2[idseg.endcap() - 1];
1576 seg_tmp = seg_ME2[idseg.endcap() - 1];
1577 } else {
1578 ++seg_ME3[idseg.endcap() - 1];
1579 seg_tmp = seg_ME3[idseg.endcap() - 1];
1580 }
1581
1582 if (1 == seg_tmp && 6 == (*segEffIt).nRecHits() && (*segEffIt).chi2() / (*segEffIt).degreesOfFreedom() < 3.) {
1583 std::pair<CSCDetId, CSCSegment> specSeg = make_pair((CSCDetId)(*segEffIt).cscDetId(), *segEffIt);
1584 theSegments[2 * (idseg.endcap() - 1) + (idseg.station() - 2)] = specSeg;
1585 }
1586 }
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603 }
1604
1605 for (int iE = 0; iE < 2; iE++) {
1606 for (int iS = 0; iS < 4; iS++) {
1607 for (int iR = 0; iR < 4; iR++) {
1608 for (int iC = 0; iC < 36; iC++) {
1609 int NumberOfLayers = 0;
1610 for (int iL = 0; iL < 6; iL++) {
1611 if (AllRecHits[iE][iS][iR][iC][iL]) {
1612 NumberOfLayers++;
1613 }
1614 }
1615 int bin = 0;
1616 if (iS == 0)
1617 bin = iR + 1 + (iE * 10);
1618 else
1619 bin = (iS + 1) * 2 + (iR + 1) + (iE * 10);
1620 if (NumberOfLayers > 1) {
1621
1622 if (AllSegments[iE][iS][iR][iC]) {
1623
1624
1625 hSSTE->Fill(bin - 0.5);
1626 if (NumberOfLayers > 3)
1627 hSSTETight->Fill(bin - 0.5);
1628 }
1629
1630
1631 hSSTE->Fill(20 + bin - 0.5);
1632 if (NumberOfLayers > 3)
1633 hSSTETight->Fill(20 + bin - 0.5);
1634
1635 }
1636 if (AllSegments[iE][iS][iR][iC]) {
1637 if (NumberOfLayers == 6) {
1638
1639
1640 hRHSTE->Fill(bin - 0.5);
1641 hRHSTETight->Fill(bin - 0.5);
1642 ;
1643 }
1644
1645
1646 hRHSTE->Fill(20 + bin - 0.5);
1647 if (NumberOfLayers > 3)
1648 hRHSTETight->Fill(20 + bin - 0.5);
1649 ;
1650 }
1651 }
1652 }
1653 }
1654 }
1655
1656
1657 std::vector<std::pair<CSCDetId, CSCSegment>*> theSeg;
1658 if (1 == seg_ME2[0])
1659 theSeg.push_back(&theSegments[0]);
1660 if (1 == seg_ME3[0])
1661 theSeg.push_back(&theSegments[1]);
1662 if (1 == seg_ME2[1])
1663 theSeg.push_back(&theSegments[2]);
1664 if (1 == seg_ME3[1])
1665 theSeg.push_back(&theSegments[3]);
1666
1667
1668
1669
1670 std::map<std::string, float> chamberTypes;
1671 chamberTypes["ME1/a"] = 0.5;
1672 chamberTypes["ME1/b"] = 1.5;
1673 chamberTypes["ME1/2"] = 2.5;
1674 chamberTypes["ME1/3"] = 3.5;
1675 chamberTypes["ME2/1"] = 4.5;
1676 chamberTypes["ME2/2"] = 5.5;
1677 chamberTypes["ME3/1"] = 6.5;
1678 chamberTypes["ME3/2"] = 7.5;
1679 chamberTypes["ME4/1"] = 8.5;
1680 chamberTypes["ME4/2"] = 9.5;
1681
1682 if (!theSeg.empty()) {
1683 std::map<int, GlobalPoint> extrapolatedPoint;
1684 std::map<int, GlobalPoint>::iterator it;
1685 const CSCGeometry::ChamberContainer& ChamberContainer = cscGeom->chambers();
1686
1687 for (size_t nCh = 0; nCh < ChamberContainer.size(); nCh++) {
1688 const CSCChamber* cscchamber = ChamberContainer[nCh];
1689 std::pair<CSCDetId, CSCSegment>* thisSegment = nullptr;
1690 for (size_t iSeg = 0; iSeg < theSeg.size(); ++iSeg) {
1691 if (cscchamber->id().endcap() == theSeg[iSeg]->first.endcap()) {
1692 if (1 == cscchamber->id().station() || 3 == cscchamber->id().station()) {
1693 if (2 == theSeg[iSeg]->first.station()) {
1694 thisSegment = theSeg[iSeg];
1695 }
1696 } else if (2 == cscchamber->id().station() || 4 == cscchamber->id().station()) {
1697 if (3 == theSeg[iSeg]->first.station()) {
1698 thisSegment = theSeg[iSeg];
1699 }
1700 }
1701 }
1702 }
1703
1704 if (thisSegment) {
1705 CSCSegment* seg = &(thisSegment->second);
1706 const CSCChamber* segChamber = cscGeom->chamber(thisSegment->first);
1707 LocalPoint localCenter(0., 0., 0);
1708 GlobalPoint cscchamberCenter = cscchamber->toGlobal(localCenter);
1709
1710 it = extrapolatedPoint.find(int(cscchamberCenter.z()));
1711 if (it == extrapolatedPoint.end()) {
1712 GlobalPoint segPos = segChamber->toGlobal(seg->localPosition());
1713 GlobalVector segDir = segChamber->toGlobal(seg->localDirection());
1714 double paramaterLine = lineParametrization(segPos.z(), cscchamberCenter.z(), segDir.z());
1715 double xExtrapolated = extrapolate1D(segPos.x(), segDir.x(), paramaterLine);
1716 double yExtrapolated = extrapolate1D(segPos.y(), segDir.y(), paramaterLine);
1717 GlobalPoint globP(xExtrapolated, yExtrapolated, cscchamberCenter.z());
1718 extrapolatedPoint[int(cscchamberCenter.z())] = globP;
1719 }
1720
1721 LocalPoint extrapolatedPointLocal = cscchamber->toLocal(extrapolatedPoint[int(cscchamberCenter.z())]);
1722 const CSCLayer* layer_p = cscchamber->layer(1);
1723 const CSCLayerGeometry* layerGeom = layer_p->geometry();
1724 const std::array<const float, 4>& layerBounds = layerGeom->parameters();
1725 float shiftFromEdge = 15.;
1726 float shiftFromDeadZone = 10.;
1727
1728 bool pass = withinSensitiveRegion(extrapolatedPointLocal,
1729 layerBounds,
1730 cscchamber->id().station(),
1731 cscchamber->id().ring(),
1732 shiftFromEdge,
1733 shiftFromDeadZone);
1734 if (pass) {
1735
1736
1737
1738
1739
1740 int nRHLayers = 0;
1741 for (int iL = 0; iL < 6; ++iL) {
1742 if (AllRecHits[cscchamber->id().endcap() - 1][cscchamber->id().station() - 1][cscchamber->id().ring() - 1]
1743 [cscchamber->id().chamber() - 1][iL]) {
1744 ++nRHLayers;
1745 }
1746 }
1747
1748 float verticalScale = chamberTypes[cscchamber->specs()->chamberTypeName()];
1749 if (cscchamberCenter.z() < 0) {
1750 verticalScale = -verticalScale;
1751 }
1752 verticalScale += 10.5;
1753 hSensitiveAreaEvt->Fill(float(cscchamber->id().chamber()), verticalScale);
1754 if (nRHLayers > 1) {
1755
1756
1757
1758
1759 hEffDenominator->Fill(float(cscchamber->id().chamber()), verticalScale);
1760 if (nRHLayers > 3)
1761 hEffDenominatorTight->Fill(float(cscchamber->id().chamber()), verticalScale);
1762
1763 if (AllSegments[cscchamber->id().endcap() - 1][cscchamber->id().station() - 1][cscchamber->id().ring() - 1]
1764 [cscchamber->id().chamber() - 1]) {
1765 hSSTE2->Fill(float(cscchamber->id().chamber()), float(verticalScale));
1766 if (nRHLayers > 3)
1767 hSSTE2Tight->Fill(float(cscchamber->id().chamber()), float(verticalScale));
1768 }
1769
1770 for (int iL = 0; iL < 6; ++iL) {
1771 float weight = 1. / 6.;
1772
1773
1774 if (AllRecHits[cscchamber->id().endcap() - 1][cscchamber->id().station() - 1][cscchamber->id().ring() - 1]
1775 [cscchamber->id().chamber() - 1][iL]) {
1776 hRHSTE2->Fill(float(cscchamber->id().chamber()), float(verticalScale), weight);
1777 if (nRHLayers > 3)
1778 hRHSTE2Tight->Fill(float(cscchamber->id().chamber()), float(verticalScale), weight);
1779 }
1780 if (useDigis) {
1781
1782 if (allWires[cscchamber->id().endcap() - 1][cscchamber->id().station() - 1][cscchamber->id().ring() - 1]
1783 [cscchamber->id().chamber() - 1][iL]) {
1784
1785 hWireSTE2->Fill(float(cscchamber->id().chamber()), float(verticalScale), weight);
1786 if (nRHLayers > 3)
1787 hWireSTE2Tight->Fill(float(cscchamber->id().chamber()), float(verticalScale), weight);
1788 }
1789
1790 if (allStrips[cscchamber->id().endcap() - 1][cscchamber->id().station() - 1]
1791 [cscchamber->id().ring() - 1][cscchamber->id().chamber() - 1][iL]) {
1792
1793 hStripSTE2->Fill(float(cscchamber->id().chamber()), float(verticalScale), weight);
1794 if (nRHLayers > 3)
1795 hStripSTE2Tight->Fill(float(cscchamber->id().chamber()), float(verticalScale), weight);
1796 }
1797 }
1798 }
1799 }
1800 }
1801 }
1802 }
1803 }
1804
1805 }
1806
1807 void CSCValidation::getEfficiency(float bin, float Norm, std::vector<float>& eff) {
1808
1809 float Efficiency = 0.;
1810 float EffError = 0.;
1811 if (fabs(Norm) > 0.000000001) {
1812 Efficiency = bin / Norm;
1813 if (bin < Norm) {
1814 EffError = sqrt((1. - Efficiency) * Efficiency / Norm);
1815 }
1816 }
1817 eff[0] = Efficiency;
1818 eff[1] = EffError;
1819 }
1820
1821 void CSCValidation::histoEfficiency(TH1F* readHisto, TH1F* writeHisto) {
1822 std::vector<float> eff(2);
1823 int Nbins = readHisto->GetSize() - 2;
1824 std::vector<float> bins(Nbins);
1825 std::vector<float> Efficiency(Nbins);
1826 std::vector<float> EffError(Nbins);
1827 float Num = 1;
1828 float Den = 1;
1829 for (int i = 0; i < 20; i++) {
1830 Num = readHisto->GetBinContent(i + 1);
1831 Den = readHisto->GetBinContent(i + 21);
1832 getEfficiency(Num, Den, eff);
1833 Efficiency[i] = eff[0];
1834 EffError[i] = eff[1];
1835 writeHisto->SetBinContent(i + 1, Efficiency[i]);
1836 writeHisto->SetBinError(i + 1, EffError[i]);
1837 }
1838 }
1839
1840 bool CSCValidation::withinSensitiveRegion(LocalPoint localPos,
1841 const std::array<const float, 4>& layerBounds,
1842 int station,
1843 int ring,
1844 float shiftFromEdge,
1845 float shiftFromDeadZone) {
1846
1847 bool pass = false;
1848
1849 float y_center = 0.;
1850 double yUp = layerBounds[3] + y_center;
1851 double yDown = -layerBounds[3] + y_center;
1852 double xBound1Shifted = layerBounds[0] - shiftFromEdge;
1853 double xBound2Shifted = layerBounds[1] - shiftFromEdge;
1854 double lineSlope = (yUp - yDown) / (xBound2Shifted - xBound1Shifted);
1855 double lineConst = yUp - lineSlope * xBound2Shifted;
1856 double yBorder = lineSlope * abs(localPos.x()) + lineConst;
1857
1858
1859 std::vector<float> deadZoneCenter(6);
1860 float cutZone = shiftFromDeadZone;
1861
1862 if (station > 1 && station < 5) {
1863 if (2 == ring) {
1864 deadZoneCenter[0] = -162.48;
1865 deadZoneCenter[1] = -81.8744;
1866 deadZoneCenter[2] = -21.18165;
1867 deadZoneCenter[3] = 39.51105;
1868 deadZoneCenter[4] = 100.2939;
1869 deadZoneCenter[5] = 160.58;
1870
1871 if (localPos.y() > yBorder &&
1872 ((localPos.y() > deadZoneCenter[0] + cutZone && localPos.y() < deadZoneCenter[1] - cutZone) ||
1873 (localPos.y() > deadZoneCenter[1] + cutZone && localPos.y() < deadZoneCenter[2] - cutZone) ||
1874 (localPos.y() > deadZoneCenter[2] + cutZone && localPos.y() < deadZoneCenter[3] - cutZone) ||
1875 (localPos.y() > deadZoneCenter[3] + cutZone && localPos.y() < deadZoneCenter[4] - cutZone) ||
1876 (localPos.y() > deadZoneCenter[4] + cutZone && localPos.y() < deadZoneCenter[5] - cutZone))) {
1877 pass = true;
1878 }
1879 } else if (1 == ring) {
1880 if (2 == station) {
1881 deadZoneCenter[0] = -95.80;
1882 deadZoneCenter[1] = -27.47;
1883 deadZoneCenter[2] = 33.67;
1884 deadZoneCenter[3] = 90.85;
1885 } else if (3 == station) {
1886 deadZoneCenter[0] = -89.305;
1887 deadZoneCenter[1] = -39.705;
1888 deadZoneCenter[2] = 20.195;
1889 deadZoneCenter[3] = 77.395;
1890 } else if (4 == station) {
1891 deadZoneCenter[0] = -75.645;
1892 deadZoneCenter[1] = -26.055;
1893 deadZoneCenter[2] = 23.855;
1894 deadZoneCenter[3] = 70.575;
1895 }
1896 if (localPos.y() > yBorder &&
1897 ((localPos.y() > deadZoneCenter[0] + cutZone && localPos.y() < deadZoneCenter[1] - cutZone) ||
1898 (localPos.y() > deadZoneCenter[1] + cutZone && localPos.y() < deadZoneCenter[2] - cutZone) ||
1899 (localPos.y() > deadZoneCenter[2] + cutZone && localPos.y() < deadZoneCenter[3] - cutZone))) {
1900 pass = true;
1901 }
1902 }
1903 } else if (1 == station) {
1904 if (3 == ring) {
1905 deadZoneCenter[0] = -83.155;
1906 deadZoneCenter[1] = -22.7401;
1907 deadZoneCenter[2] = 27.86665;
1908 deadZoneCenter[3] = 81.005;
1909 if (localPos.y() > yBorder &&
1910 ((localPos.y() > deadZoneCenter[0] + cutZone && localPos.y() < deadZoneCenter[1] - cutZone) ||
1911 (localPos.y() > deadZoneCenter[1] + cutZone && localPos.y() < deadZoneCenter[2] - cutZone) ||
1912 (localPos.y() > deadZoneCenter[2] + cutZone && localPos.y() < deadZoneCenter[3] - cutZone))) {
1913 pass = true;
1914 }
1915 } else if (2 == ring) {
1916 deadZoneCenter[0] = -86.285;
1917 deadZoneCenter[1] = -32.88305;
1918 deadZoneCenter[2] = 32.867423;
1919 deadZoneCenter[3] = 88.205;
1920 if (localPos.y() > (yBorder) &&
1921 ((localPos.y() > deadZoneCenter[0] + cutZone && localPos.y() < deadZoneCenter[1] - cutZone) ||
1922 (localPos.y() > deadZoneCenter[1] + cutZone && localPos.y() < deadZoneCenter[2] - cutZone) ||
1923 (localPos.y() > deadZoneCenter[2] + cutZone && localPos.y() < deadZoneCenter[3] - cutZone))) {
1924 pass = true;
1925 }
1926 } else if (1 == ring) {
1927 deadZoneCenter[0] = -31.5;
1928 deadZoneCenter[1] = 86.0;
1929 if (localPos.y() > (yBorder) &&
1930 (localPos.y() > deadZoneCenter[0] && localPos.y() < deadZoneCenter[1] - cutZone)) {
1931 pass = true;
1932 }
1933 } else if (4 == ring) {
1934 deadZoneCenter[0] = -86.0;
1935 deadZoneCenter[1] = -31.5;
1936 if (localPos.y() > (yBorder) &&
1937 (localPos.y() > deadZoneCenter[0] + cutZone && localPos.y() < deadZoneCenter[1])) {
1938 pass = true;
1939 }
1940 }
1941 }
1942 return pass;
1943 }
1944
1945
1946
1947
1948
1949
1950
1951
1952 float CSCValidation::getSignal(const CSCStripDigiCollection& stripdigis, CSCDetId idCS, int centerStrip) {
1953 float SigADC[5];
1954 float TotalADC = 0;
1955 SigADC[0] = 0;
1956 SigADC[1] = 0;
1957 SigADC[2] = 0;
1958 SigADC[3] = 0;
1959 SigADC[4] = 0;
1960
1961
1962 CSCStripDigiCollection::DigiRangeIterator sIt;
1963
1964 for (sIt = stripdigis.begin(); sIt != stripdigis.end(); sIt++) {
1965 CSCDetId id = (CSCDetId)(*sIt).first;
1966 if (id == idCS) {
1967
1968 std::vector<CSCStripDigi>::const_iterator digiItr = (*sIt).second.first;
1969 std::vector<CSCStripDigi>::const_iterator last = (*sIt).second.second;
1970 for (; digiItr != last; ++digiItr) {
1971 int thisStrip = digiItr->getStrip();
1972 if (thisStrip == (centerStrip)) {
1973 std::vector<int> myADCVals = digiItr->getADCCounts();
1974 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
1975 float thisSignal = (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
1976 SigADC[0] = thisSignal - 6 * thisPedestal;
1977 }
1978
1979 if (thisStrip == (centerStrip + 1)) {
1980 std::vector<int> myADCVals = digiItr->getADCCounts();
1981 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
1982 float thisSignal = (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
1983 SigADC[1] = thisSignal - 6 * thisPedestal;
1984 }
1985 if (thisStrip == (centerStrip + 2)) {
1986 std::vector<int> myADCVals = digiItr->getADCCounts();
1987 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
1988 float thisSignal = (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
1989 SigADC[2] = thisSignal - 6 * thisPedestal;
1990 }
1991 if (thisStrip == (centerStrip - 1)) {
1992 std::vector<int> myADCVals = digiItr->getADCCounts();
1993 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
1994 float thisSignal = (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
1995 SigADC[3] = thisSignal - 6 * thisPedestal;
1996 }
1997 if (thisStrip == (centerStrip - 2)) {
1998 std::vector<int> myADCVals = digiItr->getADCCounts();
1999 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
2000 float thisSignal = (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
2001 SigADC[4] = thisSignal - 6 * thisPedestal;
2002 }
2003 }
2004 TotalADC = 0.2 * (SigADC[0] + SigADC[1] + SigADC[2] + SigADC[3] + SigADC[4]);
2005 }
2006 }
2007 return TotalADC;
2008 }
2009
2010
2011
2012
2013
2014
2015 void CSCValidation::doNoiseHits(edm::Handle<CSCRecHit2DCollection> recHits,
2016 edm::Handle<CSCSegmentCollection> cscSegments,
2017 edm::ESHandle<CSCGeometry> cscGeom,
2018 edm::Handle<CSCStripDigiCollection> strips) {
2019 CSCRecHit2DCollection::const_iterator recIt;
2020 for (recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
2021 CSCDetId idrec = (CSCDetId)(*recIt).cscDetId();
2022
2023
2024 AllRechits.insert(std::pair<CSCDetId, CSCRecHit2D>(idrec, *recIt));
2025
2026
2027 int centerid = recIt->nStrips() / 2;
2028 int centerStrip = recIt->channels(centerid);
2029
2030 float rHsignal = getthisSignal(*strips, idrec, centerStrip);
2031 histos->fill1DHist(
2032 rHsignal, "hrHSignal", "Signal in the 4th time bin for centre strip", 1100, -99, 1000, "recHits");
2033 }
2034
2035 for (CSCSegmentCollection::const_iterator it = cscSegments->begin(); it != cscSegments->end(); it++) {
2036 std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits();
2037 for (std::vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
2038 CSCDetId idRH = (CSCDetId)(*iRH).cscDetId();
2039 LocalPoint lpRH = (*iRH).localPosition();
2040 float xrec = lpRH.x();
2041 float yrec = lpRH.y();
2042 float zrec = lpRH.z();
2043 bool RHalreadyinMap = false;
2044
2045 multimap<CSCDetId, CSCRecHit2D>::iterator segRHit;
2046 segRHit = SegRechits.find(idRH);
2047 if (segRHit != SegRechits.end()) {
2048 for (; segRHit != SegRechits.upper_bound(idRH); ++segRHit) {
2049
2050 LocalPoint lposRH = (segRHit->second).localPosition();
2051 float xpos = lposRH.x();
2052 float ypos = lposRH.y();
2053 float zpos = lposRH.z();
2054 if (xrec == xpos && yrec == ypos && zrec == zpos) {
2055 RHalreadyinMap = true;
2056
2057 break;
2058 }
2059 }
2060 }
2061 if (!RHalreadyinMap) {
2062 SegRechits.insert(std::pair<CSCDetId, CSCRecHit2D>(idRH, *iRH));
2063 }
2064 }
2065 }
2066
2067 findNonAssociatedRecHits(cscGeom, strips);
2068 }
2069
2070
2071
2072
2073
2074
2075
2076 void CSCValidation::findNonAssociatedRecHits(edm::ESHandle<CSCGeometry> cscGeom,
2077 edm::Handle<CSCStripDigiCollection> strips) {
2078 for (std::multimap<CSCDetId, CSCRecHit2D>::iterator allRHiter = AllRechits.begin(); allRHiter != AllRechits.end();
2079 ++allRHiter) {
2080 CSCDetId idRH = allRHiter->first;
2081 LocalPoint lpRH = (allRHiter->second).localPosition();
2082 float xrec = lpRH.x();
2083 float yrec = lpRH.y();
2084 float zrec = lpRH.z();
2085
2086 bool foundmatch = false;
2087 multimap<CSCDetId, CSCRecHit2D>::iterator segRHit;
2088 segRHit = SegRechits.find(idRH);
2089 if (segRHit != SegRechits.end()) {
2090 for (; segRHit != SegRechits.upper_bound(idRH); ++segRHit) {
2091 LocalPoint lposRH = (segRHit->second).localPosition();
2092 float xpos = lposRH.x();
2093 float ypos = lposRH.y();
2094 float zpos = lposRH.z();
2095
2096 if (xrec == xpos && yrec == ypos && zrec == zpos) {
2097 foundmatch = true;
2098 }
2099
2100 float d = 0.;
2101 float dclose = 1000.;
2102
2103 if (!foundmatch) {
2104 d = sqrt(pow(xrec - xpos, 2) + pow(yrec - ypos, 2) + pow(zrec - zpos, 2));
2105 if (d < dclose) {
2106 dclose = d;
2107 if (distRHmap.find((allRHiter->second)) ==
2108 distRHmap.end()) {
2109 distRHmap.insert(make_pair(allRHiter->second, dclose));
2110 } else {
2111
2112 distRHmap.erase(allRHiter->second);
2113 distRHmap.insert(
2114 make_pair(allRHiter->second, dclose));
2115 }
2116 }
2117 }
2118 }
2119 }
2120 if (!foundmatch) {
2121 NonAssociatedRechits.insert(std::pair<CSCDetId, CSCRecHit2D>(idRH, allRHiter->second));
2122 }
2123 }
2124
2125 for (std::map<CSCRecHit2D, float, ltrh>::iterator iter = distRHmap.begin(); iter != distRHmap.end(); ++iter) {
2126 histos->fill1DHist(iter->second,
2127 "hdistRH",
2128 "Distance of Non Associated RecHit from closest Segment RecHit",
2129 500,
2130 0.,
2131 100.,
2132 "NonAssociatedRechits");
2133 }
2134
2135 for (std::multimap<CSCDetId, CSCRecHit2D>::iterator iter = NonAssociatedRechits.begin();
2136 iter != NonAssociatedRechits.end();
2137 ++iter) {
2138 CSCDetId idrec = iter->first;
2139 int kEndcap = idrec.endcap();
2140 int cEndcap = idrec.endcap();
2141 if (kEndcap == 2)
2142 cEndcap = -1;
2143 int kRing = idrec.ring();
2144 int kStation = idrec.station();
2145 int kChamber = idrec.chamber();
2146 int kLayer = idrec.layer();
2147
2148
2149 LocalPoint rhitlocal = (iter->second).localPosition();
2150 float xreco = rhitlocal.x();
2151 float yreco = rhitlocal.y();
2152
2153
2154 int centerid = (iter->second).nStrips() / 2;
2155 int centerStrip = (iter->second).channels(centerid);
2156
2157
2158 float rHSumQ = 0;
2159 float sumsides = 0.;
2160 int adcsize = (iter->second).nStrips() * (iter->second).nTimeBins();
2161 for (unsigned int i = 0; i < (iter->second).nStrips(); i++) {
2162 for (unsigned int j = 0; j < (iter->second).nTimeBins() - 1; j++) {
2163 rHSumQ += (iter->second).adcs(i, j);
2164 if (i != 1)
2165 sumsides += (iter->second).adcs(i, j);
2166 }
2167 }
2168
2169 float rHratioQ = sumsides / rHSumQ;
2170 if (adcsize != 12)
2171 rHratioQ = -99;
2172
2173
2174 float rHtime = (iter->second).tpeak() / 50;
2175
2176
2177 int rHwidth = getWidth(*strips, idrec, centerStrip);
2178
2179
2180 const CSCLayer* csclayer = cscGeom->layer(idrec);
2181
2182
2183 GlobalPoint rhitglobal = csclayer->toGlobal(rhitlocal);
2184 float grecx = rhitglobal.x();
2185 float grecy = rhitglobal.y();
2186
2187
2188 int kCodeBroad = cEndcap * (4 * (kStation - 1) + kRing);
2189 int kCodeNarrow = cEndcap * (100 * (kRing - 1) + kChamber);
2190
2191
2192 histos->fill1DHist(
2193 kCodeBroad, "hNARHCodeBroad", "broad scope code for recHits", 33, -16.5, 16.5, "NonAssociatedRechits");
2194 if (kStation == 1)
2195 histos->fill1DHist(kCodeNarrow,
2196 "hNARHCodeNarrow1",
2197 "narrow scope recHit code station 1",
2198 801,
2199 -400.5,
2200 400.5,
2201 "NonAssociatedRechits");
2202 if (kStation == 2)
2203 histos->fill1DHist(kCodeNarrow,
2204 "hNARHCodeNarrow2",
2205 "narrow scope recHit code station 2",
2206 801,
2207 -400.5,
2208 400.5,
2209 "NonAssociatedRechits");
2210 if (kStation == 3)
2211 histos->fill1DHist(kCodeNarrow,
2212 "hNARHCodeNarrow3",
2213 "narrow scope recHit code station 3",
2214 801,
2215 -400.5,
2216 400.5,
2217 "NonAssociatedRechits");
2218 if (kStation == 4)
2219 histos->fill1DHist(kCodeNarrow,
2220 "hNARHCodeNarrow4",
2221 "narrow scope recHit code station 4",
2222 801,
2223 -400.5,
2224 400.5,
2225 "NonAssociatedRechits");
2226 histos->fill1DHistByType(kLayer, "hNARHLayer", "RecHits per Layer", idrec, 8, -0.5, 7.5, "NonAssociatedRechits");
2227 histos->fill1DHistByType(xreco, "hNARHX", "Local X of recHit", idrec, 160, -80., 80., "NonAssociatedRechits");
2228 histos->fill1DHistByType(yreco, "hNARHY", "Local Y of recHit", idrec, 60, -180., 180., "NonAssociatedRechits");
2229 if (kStation == 1 && (kRing == 1 || kRing == 4))
2230 histos->fill1DHistByType(
2231 rHSumQ, "hNARHSumQ", "Sum 3x3 recHit Charge", idrec, 250, 0, 4000, "NonAssociatedRechits");
2232 else
2233 histos->fill1DHistByType(
2234 rHSumQ, "hNARHSumQ", "Sum 3x3 recHit Charge", idrec, 250, 0, 2000, "NonAssociatedRechits");
2235 histos->fill1DHistByType(
2236 rHratioQ, "hNARHRatioQ", "Ratio (Ql+Qr)/Qt)", idrec, 120, -0.1, 1.1, "NonAssociatedRechits");
2237 histos->fill1DHistByType(rHtime, "hNARHTiming", "recHit Timing", idrec, 200, -10, 10, "NonAssociatedRechits");
2238 histos->fill2DHistByStation(grecx,
2239 grecy,
2240 "hNARHGlobal",
2241 "recHit Global Position",
2242 idrec,
2243 400,
2244 -800.,
2245 800.,
2246 400,
2247 -800.,
2248 800.,
2249 "NonAssociatedRechits");
2250 histos->fill1DHistByType(
2251 rHwidth, "hNARHwidth", "width for Non associated recHit", idrec, 21, -0.5, 20.5, "NonAssociatedRechits");
2252 }
2253
2254 for (std::multimap<CSCDetId, CSCRecHit2D>::iterator iter = SegRechits.begin(); iter != SegRechits.end(); ++iter) {
2255 CSCDetId idrec = iter->first;
2256 int kEndcap = idrec.endcap();
2257 int cEndcap = idrec.endcap();
2258 if (kEndcap == 2)
2259 cEndcap = -1;
2260 int kRing = idrec.ring();
2261 int kStation = idrec.station();
2262 int kChamber = idrec.chamber();
2263 int kLayer = idrec.layer();
2264
2265
2266 LocalPoint rhitlocal = (iter->second).localPosition();
2267 float xreco = rhitlocal.x();
2268 float yreco = rhitlocal.y();
2269
2270
2271 int centerid = (iter->second).nStrips() / 2;
2272 int centerStrip = (iter->second).channels(centerid);
2273
2274
2275
2276 float rHSumQ = 0;
2277 float sumsides = 0.;
2278 int adcsize = (iter->second).nStrips() * (iter->second).nTimeBins();
2279 for (unsigned int i = 0; i < (iter->second).nStrips(); i++) {
2280 for (unsigned int j = 0; j < (iter->second).nTimeBins() - 1; j++) {
2281 rHSumQ += (iter->second).adcs(i, j);
2282 if (i != 1)
2283 sumsides += (iter->second).adcs(i, j);
2284 }
2285 }
2286
2287 float rHratioQ = sumsides / rHSumQ;
2288 if (adcsize != 12)
2289 rHratioQ = -99;
2290
2291
2292 float rHtime = (iter->second).tpeak() / 50;
2293
2294
2295 int rHwidth = getWidth(*strips, idrec, centerStrip);
2296
2297
2298 const CSCLayer* csclayer = cscGeom->layer(idrec);
2299
2300
2301 GlobalPoint rhitglobal = csclayer->toGlobal(rhitlocal);
2302 float grecx = rhitglobal.x();
2303 float grecy = rhitglobal.y();
2304
2305
2306 int kCodeBroad = cEndcap * (4 * (kStation - 1) + kRing);
2307 int kCodeNarrow = cEndcap * (100 * (kRing - 1) + kChamber);
2308
2309
2310 histos->fill1DHist(
2311 kCodeBroad, "hSegRHCodeBroad", "broad scope code for recHits", 33, -16.5, 16.5, "AssociatedRechits");
2312 if (kStation == 1)
2313 histos->fill1DHist(kCodeNarrow,
2314 "hSegRHCodeNarrow1",
2315 "narrow scope recHit code station 1",
2316 801,
2317 -400.5,
2318 400.5,
2319 "AssociatedRechits");
2320 if (kStation == 2)
2321 histos->fill1DHist(kCodeNarrow,
2322 "hSegRHCodeNarrow2",
2323 "narrow scope recHit code station 2",
2324 801,
2325 -400.5,
2326 400.5,
2327 "AssociatedRechits");
2328 if (kStation == 3)
2329 histos->fill1DHist(kCodeNarrow,
2330 "hSegRHCodeNarrow3",
2331 "narrow scope recHit code station 3",
2332 801,
2333 -400.5,
2334 400.5,
2335 "AssociatedRechits");
2336 if (kStation == 4)
2337 histos->fill1DHist(kCodeNarrow,
2338 "hSegRHCodeNarrow4",
2339 "narrow scope recHit code station 4",
2340 801,
2341 -400.5,
2342 400.5,
2343 "AssociatedRechits");
2344 histos->fill1DHistByType(kLayer, "hSegRHLayer", "RecHits per Layer", idrec, 8, -0.5, 7.5, "AssociatedRechits");
2345 histos->fill1DHistByType(xreco, "hSegRHX", "Local X of recHit", idrec, 160, -80., 80., "AssociatedRechits");
2346 histos->fill1DHistByType(yreco, "hSegRHY", "Local Y of recHit", idrec, 60, -180., 180., "AssociatedRechits");
2347 if (kStation == 1 && (kRing == 1 || kRing == 4))
2348 histos->fill1DHistByType(rHSumQ, "hSegRHSumQ", "Sum 3x3 recHit Charge", idrec, 250, 0, 4000, "AssociatedRechits");
2349 else
2350 histos->fill1DHistByType(rHSumQ, "hSegRHSumQ", "Sum 3x3 recHit Charge", idrec, 250, 0, 2000, "AssociatedRechits");
2351 histos->fill1DHistByType(rHratioQ, "hSegRHRatioQ", "Ratio (Ql+Qr)/Qt)", idrec, 120, -0.1, 1.1, "AssociatedRechits");
2352 histos->fill1DHistByType(rHtime, "hSegRHTiming", "recHit Timing", idrec, 200, -10, 10, "AssociatedRechits");
2353 histos->fill2DHistByStation(grecx,
2354 grecy,
2355 "hSegRHGlobal",
2356 "recHit Global Position",
2357 idrec,
2358 400,
2359 -800.,
2360 800.,
2361 400,
2362 -800.,
2363 800.,
2364 "AssociatedRechits");
2365 histos->fill1DHistByType(
2366 rHwidth, "hSegRHwidth", "width for Non associated recHit", idrec, 21, -0.5, 20.5, "AssociatedRechits");
2367 }
2368
2369 distRHmap.clear();
2370 AllRechits.clear();
2371 SegRechits.clear();
2372 NonAssociatedRechits.clear();
2373 }
2374
2375 float CSCValidation::getthisSignal(const CSCStripDigiCollection& stripdigis, CSCDetId idRH, int centerStrip) {
2376
2377 CSCStripDigiCollection::DigiRangeIterator sIt;
2378 float thisADC = 0.;
2379
2380
2381 for (sIt = stripdigis.begin(); sIt != stripdigis.end(); sIt++) {
2382 CSCDetId id = (CSCDetId)(*sIt).first;
2383
2384 if (id == idRH) {
2385
2386 vector<CSCStripDigi>::const_iterator digiItr = (*sIt).second.first;
2387 vector<CSCStripDigi>::const_iterator last = (*sIt).second.second;
2388
2389 int St = idRH.station();
2390 int Rg = idRH.ring();
2391 if (St == 1 && Rg == 4) {
2392 while (centerStrip > 16)
2393 centerStrip -= 16;
2394 }
2395 for (; digiItr != last; ++digiItr) {
2396 int thisStrip = digiItr->getStrip();
2397
2398 std::vector<int> myADCVals = digiItr->getADCCounts();
2399 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
2400 float Signal = (float)myADCVals[3];
2401 if (thisStrip == (centerStrip)) {
2402 thisADC = Signal - thisPedestal;
2403
2404
2405
2406
2407 }
2408 if (thisStrip == (centerStrip + 1)) {
2409 std::vector<int> myADCVals = digiItr->getADCCounts();
2410 }
2411 if (thisStrip == (centerStrip - 1)) {
2412 std::vector<int> myADCVals = digiItr->getADCCounts();
2413 }
2414 }
2415 }
2416 }
2417
2418 return thisADC;
2419 }
2420
2421
2422
2423
2424
2425
2426
2427 int CSCValidation::getWidth(const CSCStripDigiCollection& stripdigis, CSCDetId idRH, int centerStrip) {
2428 int width = 1;
2429 int widthpos = 0;
2430 int widthneg = 0;
2431
2432
2433 CSCStripDigiCollection::DigiRangeIterator sIt;
2434
2435 for (sIt = stripdigis.begin(); sIt != stripdigis.end(); sIt++) {
2436 CSCDetId id = (CSCDetId)(*sIt).first;
2437 if (id == idRH) {
2438 std::vector<CSCStripDigi>::const_iterator digiItr = (*sIt).second.first;
2439 std::vector<CSCStripDigi>::const_iterator first = (*sIt).second.first;
2440 std::vector<CSCStripDigi>::const_iterator last = (*sIt).second.second;
2441 std::vector<CSCStripDigi>::const_iterator it = (*sIt).second.first;
2442 std::vector<CSCStripDigi>::const_iterator itr = (*sIt).second.first;
2443
2444 int St = idRH.station();
2445 int Rg = idRH.ring();
2446 if (St == 1 && Rg == 4) {
2447 while (centerStrip > 16)
2448 centerStrip -= 16;
2449 }
2450 for (; digiItr != last; ++digiItr) {
2451 int thisStrip = digiItr->getStrip();
2452 if (thisStrip == (centerStrip)) {
2453 it = digiItr;
2454 for (; it != last; ++it) {
2455 int strip = it->getStrip();
2456 std::vector<int> myADCVals = it->getADCCounts();
2457 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
2458 if (((float)myADCVals[3] - thisPedestal) < 6 || widthpos == 10 || it == last) {
2459 break;
2460 }
2461 if (strip != centerStrip) {
2462 widthpos += 1;
2463 }
2464 }
2465 itr = digiItr;
2466 for (; itr != first; --itr) {
2467 int strip = itr->getStrip();
2468 std::vector<int> myADCVals = itr->getADCCounts();
2469 float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
2470 if (((float)myADCVals[3] - thisPedestal) < 6 || widthneg == 10 || itr == first) {
2471 break;
2472 }
2473 if (strip != centerStrip) {
2474 widthneg += 1;
2475 }
2476 }
2477 }
2478 }
2479 }
2480 }
2481
2482 width = width + widthneg + widthpos;
2483
2484 return width;
2485 }
2486
2487
2488
2489
2490
2491
2492 void CSCValidation::doGasGain(const CSCWireDigiCollection& wirecltn,
2493 const CSCStripDigiCollection& strpcltn,
2494 const CSCRecHit2DCollection& rechitcltn) {
2495 int channel = 0, mult, wire, layer, idlayer, idchamber, ring;
2496 int wire_strip_rechit_present;
2497 std::string name, title, endcapstr;
2498 ostringstream ss;
2499 CSCIndexer indexer;
2500 std::map<int, int>::iterator intIt;
2501
2502 m_single_wire_layer.clear();
2503
2504 if (firstEvent) {
2505
2506
2507 m_wire_hvsegm.clear();
2508 std::map<int, std::vector<int> >::iterator intvecIt;
2509
2510 int csctype[10] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
2511 int hvsegm_layer[10] = {1, 1, 3, 3, 3, 5, 3, 5, 3, 5};
2512 int id;
2513 nmbhvsegm.clear();
2514 for (int i = 0; i < 10; i++)
2515 nmbhvsegm.push_back(hvsegm_layer[i]);
2516
2517 std::vector<int> zer_1_1a(49, 0);
2518 id = csctype[0];
2519 if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2520 m_wire_hvsegm[id] = zer_1_1a;
2521 intvecIt = m_wire_hvsegm.find(id);
2522 for (int wire = 1; wire <= 48; wire++)
2523 intvecIt->second[wire] = 1;
2524
2525
2526 std::vector<int> zer_1_1b(49, 0);
2527 id = csctype[1];
2528 if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2529 m_wire_hvsegm[id] = zer_1_1b;
2530 intvecIt = m_wire_hvsegm.find(id);
2531 for (int wire = 1; wire <= 48; wire++)
2532 intvecIt->second[wire] = 1;
2533
2534
2535 std::vector<int> zer_1_2(65, 0);
2536 id = csctype[2];
2537 if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2538 m_wire_hvsegm[id] = zer_1_2;
2539 intvecIt = m_wire_hvsegm.find(id);
2540 for (int wire = 1; wire <= 24; wire++)
2541 intvecIt->second[wire] = 1;
2542 for (int wire = 25; wire <= 48; wire++)
2543 intvecIt->second[wire] = 2;
2544 for (int wire = 49; wire <= 64; wire++)
2545 intvecIt->second[wire] = 3;
2546
2547
2548 std::vector<int> zer_1_3(33, 0);
2549 id = csctype[3];
2550 if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2551 m_wire_hvsegm[id] = zer_1_3;
2552 intvecIt = m_wire_hvsegm.find(id);
2553 for (int wire = 1; wire <= 12; wire++)
2554 intvecIt->second[wire] = 1;
2555 for (int wire = 13; wire <= 22; wire++)
2556 intvecIt->second[wire] = 2;
2557 for (int wire = 23; wire <= 32; wire++)
2558 intvecIt->second[wire] = 3;
2559
2560
2561 std::vector<int> zer_2_1(113, 0);
2562 id = csctype[4];
2563 if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2564 m_wire_hvsegm[id] = zer_2_1;
2565 intvecIt = m_wire_hvsegm.find(id);
2566 for (int wire = 1; wire <= 44; wire++)
2567 intvecIt->second[wire] = 1;
2568 for (int wire = 45; wire <= 80; wire++)
2569 intvecIt->second[wire] = 2;
2570 for (int wire = 81; wire <= 112; wire++)
2571 intvecIt->second[wire] = 3;
2572
2573
2574 std::vector<int> zer_2_2(65, 0);
2575 id = csctype[5];
2576 if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2577 m_wire_hvsegm[id] = zer_2_2;
2578 intvecIt = m_wire_hvsegm.find(id);
2579 for (int wire = 1; wire <= 16; wire++)
2580 intvecIt->second[wire] = 1;
2581 for (int wire = 17; wire <= 28; wire++)
2582 intvecIt->second[wire] = 2;
2583 for (int wire = 29; wire <= 40; wire++)
2584 intvecIt->second[wire] = 3;
2585 for (int wire = 41; wire <= 52; wire++)
2586 intvecIt->second[wire] = 4;
2587 for (int wire = 53; wire <= 64; wire++)
2588 intvecIt->second[wire] = 5;
2589
2590
2591 std::vector<int> zer_3_1(97, 0);
2592 id = csctype[6];
2593 if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2594 m_wire_hvsegm[id] = zer_3_1;
2595 intvecIt = m_wire_hvsegm.find(id);
2596 for (int wire = 1; wire <= 32; wire++)
2597 intvecIt->second[wire] = 1;
2598 for (int wire = 33; wire <= 64; wire++)
2599 intvecIt->second[wire] = 2;
2600 for (int wire = 65; wire <= 96; wire++)
2601 intvecIt->second[wire] = 3;
2602
2603
2604 std::vector<int> zer_3_2(65, 0);
2605 id = csctype[7];
2606 if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2607 m_wire_hvsegm[id] = zer_3_2;
2608 intvecIt = m_wire_hvsegm.find(id);
2609 for (int wire = 1; wire <= 16; wire++)
2610 intvecIt->second[wire] = 1;
2611 for (int wire = 17; wire <= 28; wire++)
2612 intvecIt->second[wire] = 2;
2613 for (int wire = 29; wire <= 40; wire++)
2614 intvecIt->second[wire] = 3;
2615 for (int wire = 41; wire <= 52; wire++)
2616 intvecIt->second[wire] = 4;
2617 for (int wire = 53; wire <= 64; wire++)
2618 intvecIt->second[wire] = 5;
2619
2620
2621 std::vector<int> zer_4_1(97, 0);
2622 id = csctype[8];
2623 if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2624 m_wire_hvsegm[id] = zer_4_1;
2625 intvecIt = m_wire_hvsegm.find(id);
2626 for (int wire = 1; wire <= 32; wire++)
2627 intvecIt->second[wire] = 1;
2628 for (int wire = 33; wire <= 64; wire++)
2629 intvecIt->second[wire] = 2;
2630 for (int wire = 65; wire <= 96; wire++)
2631 intvecIt->second[wire] = 3;
2632
2633
2634 std::vector<int> zer_4_2(65, 0);
2635 id = csctype[9];
2636 if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2637 m_wire_hvsegm[id] = zer_4_2;
2638 intvecIt = m_wire_hvsegm.find(id);
2639 for (int wire = 1; wire <= 16; wire++)
2640 intvecIt->second[wire] = 1;
2641 for (int wire = 17; wire <= 28; wire++)
2642 intvecIt->second[wire] = 2;
2643 for (int wire = 29; wire <= 40; wire++)
2644 intvecIt->second[wire] = 3;
2645 for (int wire = 41; wire <= 52; wire++)
2646 intvecIt->second[wire] = 4;
2647 for (int wire = 53; wire <= 64; wire++)
2648 intvecIt->second[wire] = 5;
2649
2650 }
2651
2652
2653 wire_strip_rechit_present = 0;
2654 if (wirecltn.begin() != wirecltn.end())
2655 wire_strip_rechit_present = wire_strip_rechit_present + 1;
2656 if (strpcltn.begin() != strpcltn.end())
2657 wire_strip_rechit_present = wire_strip_rechit_present + 2;
2658 if (rechitcltn.begin() != rechitcltn.end())
2659 wire_strip_rechit_present = wire_strip_rechit_present + 4;
2660
2661 if (wire_strip_rechit_present == 7) {
2662
2663
2664
2665
2666 CSCWireDigiCollection::DigiRangeIterator wiredetUnitIt;
2667
2668 for (wiredetUnitIt = wirecltn.begin(); wiredetUnitIt != wirecltn.end(); ++wiredetUnitIt) {
2669 const CSCDetId id = (*wiredetUnitIt).first;
2670 idlayer = indexer.dbIndex(id, channel);
2671
2672 mult = 0;
2673 wire = 0;
2674 const CSCWireDigiCollection::Range& range = (*wiredetUnitIt).second;
2675 for (CSCWireDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
2676 wire = (*digiIt).getWireGroup();
2677 mult++;
2678 }
2679
2680
2681 if (mult == 1) {
2682 if (m_single_wire_layer.find(idlayer) == m_single_wire_layer.end())
2683 m_single_wire_layer[idlayer] = wire;
2684 }
2685 }
2686
2687
2688 CSCRecHit2DCollection::const_iterator recIt;
2689 CSCRecHit2D::ADCContainer m_adc;
2690
2691 for (recIt = rechitcltn.begin(); recIt != rechitcltn.end(); ++recIt) {
2692 CSCDetId id = (CSCDetId)(*recIt).cscDetId();
2693 idlayer = indexer.dbIndex(id, channel);
2694 idchamber = idlayer / 10;
2695 layer = id.layer();
2696
2697 if (m_single_wire_layer.find(idlayer) != m_single_wire_layer.end()) {
2698 if (recIt->nStrips() == 3) {
2699
2700 unsigned int binmx = 0;
2701 float adcmax = 0.0;
2702
2703 for (unsigned int i = 0; i < recIt->nStrips(); i++)
2704 for (unsigned int j = 0; j < recIt->nTimeBins(); j++)
2705 if (recIt->adcs(i, j) > adcmax) {
2706 adcmax = recIt->adcs(i, j);
2707 binmx = j;
2708 }
2709
2710 float adc_3_3_sum = 0.0;
2711
2712 for (unsigned int i = 0; i < recIt->nStrips(); i++)
2713 for (unsigned int j = binmx - 1; j <= binmx + 1; j++)
2714 adc_3_3_sum += recIt->adcs(i, j);
2715
2716 if (adc_3_3_sum > 0.0 && adc_3_3_sum < 2000.0) {
2717
2718 int flag = 0;
2719 if (id.station() == 1 && id.ring() == 4 && recIt->channels(1) > 16)
2720 flag = 1;
2721
2722 if (flag == 0) {
2723 wire = m_single_wire_layer[idlayer];
2724 int chambertype = id.iChamberType(id.station(), id.ring());
2725 int hvsgmtnmb = m_wire_hvsegm[chambertype][wire];
2726 int nmbofhvsegm = nmbhvsegm[chambertype - 1];
2727 int location = (layer - 1) * nmbofhvsegm + hvsgmtnmb;
2728
2729 ss << "gas_gain_rechit_adc_3_3_sum_location_ME_" << idchamber;
2730 name = ss.str();
2731 ss.str("");
2732 if (id.endcap() == 1)
2733 endcapstr = "+";
2734 ring = id.ring();
2735 if (id.station() == 1 && id.ring() == 4)
2736 ring = 1;
2737 if (id.endcap() == 2)
2738 endcapstr = "-";
2739 ss << "Gas Gain Rechit ADC3X3 Sum ME" << endcapstr << id.station() << "/" << ring << "/" << id.chamber();
2740 title = ss.str();
2741 ss.str("");
2742 float x = location;
2743 float y = adc_3_3_sum;
2744 histos->fill2DHist(x, y, name, title, 30, 1.0, 31.0, 50, 0.0, 2000.0, "GasGain");
2745
2746
2747
2748
2749
2750
2751
2752 }
2753 }
2754 }
2755 }
2756 }
2757 }
2758 }
2759
2760
2761
2762
2763
2764
2765 void CSCValidation::doAFEBTiming(const CSCWireDigiCollection& wirecltn) {
2766 ostringstream ss;
2767 std::string name, title, endcapstr;
2768 float x, y;
2769 int wire, wiretbin, nmbwiretbin, layer, afeb, idlayer, idchamber;
2770 int channel = 0;
2771 CSCIndexer indexer;
2772
2773 if (wirecltn.begin() != wirecltn.end()) {
2774
2775
2776
2777
2778
2779 CSCWireDigiCollection::DigiRangeIterator wiredetUnitIt;
2780 for (wiredetUnitIt = wirecltn.begin(); wiredetUnitIt != wirecltn.end(); ++wiredetUnitIt) {
2781 const CSCDetId id = (*wiredetUnitIt).first;
2782 idlayer = indexer.dbIndex(id, channel);
2783 idchamber = idlayer / 10;
2784 layer = id.layer();
2785
2786 if (id.endcap() == 1)
2787 endcapstr = "+";
2788 if (id.endcap() == 2)
2789 endcapstr = "-";
2790
2791
2792
2793 const CSCWireDigiCollection::Range& range = (*wiredetUnitIt).second;
2794 for (CSCWireDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
2795 wire = (*digiIt).getWireGroup();
2796 wiretbin = (*digiIt).getTimeBin();
2797 nmbwiretbin = (*digiIt).getTimeBinsOn().size();
2798 afeb = 3 * ((wire - 1) / 8) + (layer + 1) / 2;
2799
2800
2801 x = afeb;
2802 y = wiretbin;
2803 ss << "afeb_time_bin_vs_afeb_occupancy_ME_" << idchamber;
2804 name = ss.str();
2805 ss.str("");
2806 ss << "Time Bin vs AFEB Occupancy ME" << endcapstr << id.station() << "/" << id.ring() << "/" << id.chamber();
2807 title = ss.str();
2808 ss.str("");
2809 histos->fill2DHist(x, y, name, title, 42, 1., 43., 16, 0., 16., "AFEBTiming");
2810
2811
2812 x = afeb;
2813 y = nmbwiretbin;
2814 ss << "nmb_afeb_time_bins_vs_afeb_ME_" << idchamber;
2815 name = ss.str();
2816 ss.str("");
2817 ss << "Number of Time Bins vs AFEB ME" << endcapstr << id.station() << "/" << id.ring() << "/" << id.chamber();
2818 title = ss.str();
2819 ss.str("");
2820 histos->fill2DHist(x, y, name, title, 42, 1., 43., 16, 0., 16., "AFEBTiming");
2821
2822 }
2823 }
2824 }
2825 }
2826
2827
2828
2829
2830
2831
2832 void CSCValidation::doCompTiming(const CSCComparatorDigiCollection& compars) {
2833 ostringstream ss;
2834 std::string name, title, endcap;
2835 float x, y;
2836 int strip, tbin, cfeb, idlayer, idchamber;
2837 int channel = 0;
2838 CSCIndexer indexer;
2839
2840 if (compars.begin() != compars.end()) {
2841
2842
2843
2844
2845
2846 CSCComparatorDigiCollection::DigiRangeIterator compdetUnitIt;
2847 for (compdetUnitIt = compars.begin(); compdetUnitIt != compars.end(); ++compdetUnitIt) {
2848 const CSCDetId id = (*compdetUnitIt).first;
2849 idlayer = indexer.dbIndex(id, channel);
2850 idchamber = idlayer / 10;
2851
2852 if (id.endcap() == 1)
2853 endcap = "+";
2854 if (id.endcap() == 2)
2855 endcap = "-";
2856
2857 const CSCComparatorDigiCollection::Range& range = (*compdetUnitIt).second;
2858 for (CSCComparatorDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
2859 strip = (*digiIt).getStrip();
2860
2861
2862
2863
2864
2865 indexer.dbIndex(id, strip);
2866
2867 tbin = (*digiIt).getTimeBin();
2868 cfeb = (strip - 1) / 16 + 1;
2869
2870
2871
2872 x = cfeb;
2873 y = tbin;
2874 ss << "comp_time_bin_vs_cfeb_occupancy_ME_" << idchamber;
2875 name = ss.str();
2876 ss.str("");
2877 ss << "Comparator Time Bin vs CFEB Occupancy ME" << endcap << id.station() << "/" << id.ring() << "/"
2878 << id.chamber();
2879 title = ss.str();
2880 ss.str("");
2881 histos->fill2DHist(x, y, name, title, 5, 1., 6., 16, 0., 16., "CompTiming");
2882
2883 }
2884 }
2885 }
2886 }
2887
2888
2889
2890
2891
2892
2893 void CSCValidation::doADCTiming(const CSCRecHit2DCollection& rechitcltn) {
2894 float adc_3_3_sum, adc_3_3_wtbin, x, y;
2895 int cfeb, idchamber, ring;
2896
2897 std::string name, title, endcapstr;
2898 ostringstream ss;
2899 std::vector<float> zer(6, 0.0);
2900
2901 CSCIndexer indexer;
2902 std::map<int, int>::iterator intIt;
2903
2904 if (rechitcltn.begin() != rechitcltn.end()) {
2905
2906
2907
2908 CSCRecHit2DCollection::const_iterator recIt;
2909 CSCRecHit2D::ADCContainer m_adc;
2910 for (recIt = rechitcltn.begin(); recIt != rechitcltn.end(); ++recIt) {
2911 CSCDetId id = (CSCDetId)(*recIt).cscDetId();
2912
2913 if (recIt->nStrips() == 3) {
2914
2915
2916 unsigned int binmx = 0;
2917 float adcmax = 0.0;
2918
2919 for (unsigned int i = 0; i < recIt->nStrips(); i++)
2920 for (unsigned int j = 0; j < recIt->nTimeBins(); j++)
2921 if (recIt->adcs(i, j) > adcmax) {
2922 adcmax = recIt->adcs(i, j);
2923 binmx = j;
2924 }
2925
2926 adc_3_3_sum = 0.0;
2927
2928 for (unsigned int i = 0; i < recIt->nStrips(); i++)
2929 for (unsigned int j = binmx - 1; j <= binmx + 1; j++)
2930 adc_3_3_sum += recIt->adcs(i, j);
2931
2932
2933 if (adc_3_3_sum > 100.0) {
2934 int centerStrip = recIt->channels(1);
2935
2936 int flag = 0;
2937 if (id.station() == 1 && id.ring() == 4 && centerStrip > 16)
2938 flag = 1;
2939
2940 if (flag == 0) {
2941 adc_3_3_wtbin = (*recIt).tpeak() / 50;
2942 idchamber = indexer.dbIndex(id, centerStrip) / 10;
2943
2944
2945
2946
2947
2948
2949
2950 ss << "adc_3_3_weight_time_bin_vs_cfeb_occupancy_ME_" << idchamber;
2951 name = ss.str();
2952 ss.str("");
2953
2954 std::string endcapstr;
2955 if (id.endcap() == 1)
2956 endcapstr = "+";
2957 if (id.endcap() == 2)
2958 endcapstr = "-";
2959 ring = id.ring();
2960 if (id.ring() == 4)
2961 ring = 1;
2962 ss << "ADC 3X3 Weighted Time Bin vs CFEB Occupancy ME" << endcapstr << id.station() << "/" << ring << "/"
2963 << id.chamber();
2964 title = ss.str();
2965 ss.str("");
2966
2967 cfeb = (centerStrip - 1) / 16 + 1;
2968 x = cfeb;
2969 y = adc_3_3_wtbin;
2970 histos->fill2DHist(x, y, name, title, 5, 1., 6., 80, -8., 8., "ADCTiming");
2971 }
2972 }
2973 }
2974 }
2975 }
2976 }
2977
2978
2979
2980
2981
2982
2983 void CSCValidation::doTimeMonitoring(edm::Handle<CSCRecHit2DCollection> recHits,
2984 edm::Handle<CSCSegmentCollection> cscSegments,
2985 edm::Handle<CSCALCTDigiCollection> alcts,
2986 edm::Handle<CSCCLCTDigiCollection> clcts,
2987 edm::Handle<CSCCorrelatedLCTDigiCollection> correlatedlcts,
2988 edm::Handle<L1MuGMTReadoutCollection> pCollection,
2989 edm::ESHandle<CSCGeometry> cscGeom,
2990 const edm::EventSetup& eventSetup,
2991 const edm::Event& event) {
2992 map<CSCDetId, float> segment_median_map;
2993 map<CSCDetId, GlobalPoint> segment_position_map;
2994
2995
2996
2997
2998 for (CSCSegmentCollection::const_iterator dSiter = cscSegments->begin(); dSiter != cscSegments->end(); dSiter++) {
2999 CSCDetId id = (CSCDetId)(*dSiter).cscDetId();
3000 LocalPoint localPos = (*dSiter).localPosition();
3001 GlobalPoint globalPosition = GlobalPoint(0.0, 0.0, 0.0);
3002 const CSCChamber* cscchamber = cscGeom->chamber(id);
3003 if (cscchamber) {
3004 globalPosition = cscchamber->toGlobal(localPos);
3005 }
3006
3007
3008 std::vector<CSCRecHit2D> theseRecHits = (*dSiter).specificRecHits();
3009 int nRH = (*dSiter).nRecHits();
3010 if (nRH < 4)
3011 continue;
3012
3013
3014 vector<float> non_zero;
3015
3016 for (vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
3017 non_zero.push_back(iRH->tpeak());
3018
3019 }
3020
3021
3022 sort(non_zero.begin(), non_zero.end());
3023 int middle_index = non_zero.size() / 2;
3024 float average_two = (non_zero.at(middle_index - 1) + non_zero.at(middle_index)) / 2.;
3025 if (non_zero.size() % 2)
3026 average_two = non_zero.at(middle_index);
3027
3028
3029 segment_median_map[id] = average_two;
3030 segment_position_map[id] = globalPosition;
3031
3032 double distToIP = sqrt(globalPosition.x() * globalPosition.x() + globalPosition.y() * globalPosition.y() +
3033 globalPosition.z() * globalPosition.z());
3034
3035 histos->fillProfile(chamberSerial(id),
3036 average_two,
3037 "timeChamber",
3038 "Segment mean time",
3039 601,
3040 -0.5,
3041 600.5,
3042 -400.,
3043 400.,
3044 "TimeMonitoring");
3045 histos->fillProfileByType(id.chamber(),
3046 average_two,
3047 "timeChamberByType",
3048 "Segment mean time by chamber",
3049 id,
3050 36,
3051 0.5,
3052 36.5,
3053 -400,
3054 400.,
3055 "TimeMonitoring");
3056 histos->fill2DHist(distToIP,
3057 average_two,
3058 "seg_time_vs_distToIP",
3059 "Segment time vs. Distance to IP",
3060 80,
3061 600.,
3062 1400.,
3063 800,
3064 -400,
3065 400.,
3066 "TimeMonitoring");
3067 histos->fill2DHist(globalPosition.z(),
3068 average_two,
3069 "seg_time_vs_globZ",
3070 "Segment time vs. z position",
3071 240,
3072 -1200,
3073 1200,
3074 800,
3075 -400.,
3076 400.,
3077 "TimeMonitoring");
3078 histos->fill2DHist(fabs(globalPosition.z()),
3079 average_two,
3080 "seg_time_vs_absglobZ",
3081 "Segment time vs. abs(z position)",
3082 120,
3083 0.,
3084 1200.,
3085 800,
3086 -400.,
3087 400.,
3088 "TimeMonitoring");
3089
3090 }
3091
3092
3093 map<CSCDetId, float>::const_iterator it_outer;
3094 map<CSCDetId, float>::const_iterator it_inner;
3095 for (it_outer = segment_median_map.begin(); it_outer != segment_median_map.end(); it_outer++) {
3096 CSCDetId id_outer = it_outer->first;
3097 float t_outer = it_outer->second;
3098
3099
3100 for (it_inner = segment_median_map.begin(); it_inner != segment_median_map.end(); it_inner++) {
3101 CSCDetId id_inner = it_inner->first;
3102 float t_inner = it_inner->second;
3103
3104
3105
3106
3107 if (chamberSerial(id_outer) == chamberSerial(id_inner))
3108 continue;
3109
3110
3111
3112
3113
3114
3115
3116
3117 if (id_outer.endcap() == 1 && id_inner.endcap() == 2 && id_outer.station() == id_inner.station() &&
3118 id_outer.ring() == id_inner.ring()) {
3119 histos->fill1DHist(t_outer - t_inner,
3120 "diff_opposite_endcaps",
3121 "#Delta t [ME+]-[ME-] for chambers in same station and ring",
3122 800,
3123 -400.,
3124 400.,
3125 "TimeMonitoring");
3126 histos->fill1DHistByType(t_outer - t_inner,
3127 "diff_opposite_endcaps_byType",
3128 "#Delta t [ME+]-[ME-] for chambers in same station and ring",
3129 id_outer,
3130 800,
3131 -400.,
3132 400.,
3133 "TimeMonitoring");
3134 }
3135
3136 }
3137 }
3138
3139
3140 if (!useDigis)
3141 return;
3142
3143
3144 vector<L1MuGMTReadoutRecord> L1Mrec = pCollection->getRecords();
3145 vector<L1MuGMTReadoutRecord>::const_iterator igmtrr;
3146 int L1GMT_BXN = -100;
3147 bool has_CSCTrigger = false;
3148 bool has_beamHaloTrigger = false;
3149 for (igmtrr = L1Mrec.begin(); igmtrr != L1Mrec.end(); igmtrr++) {
3150 std::vector<L1MuRegionalCand>::const_iterator iter1;
3151 std::vector<L1MuRegionalCand> rmc;
3152
3153 int icsc = 0;
3154 rmc = igmtrr->getCSCCands();
3155 for (iter1 = rmc.begin(); iter1 != rmc.end(); iter1++) {
3156 if (!(*iter1).empty()) {
3157 icsc++;
3158 int kQuality = (*iter1).quality();
3159 if (kQuality == 1)
3160 has_beamHaloTrigger = true;
3161 }
3162 }
3163 if (igmtrr->getBxInEvent() == 0 && icsc > 0) {
3164
3165 L1GMT_BXN = igmtrr->getBxNr();
3166 has_CSCTrigger = true;
3167 } else if (igmtrr->getBxInEvent() == 0) {
3168
3169 L1GMT_BXN = igmtrr->getBxNr();
3170 }
3171 }
3172
3173
3174
3175
3176
3177 int n_alcts = 0;
3178 map<CSCDetId, int> ALCT_KeyWG_map;
3179 for (CSCALCTDigiCollection::DigiRangeIterator j = alcts->begin(); j != alcts->end(); j++) {
3180 const CSCALCTDigiCollection::Range& range = (*j).second;
3181 const CSCDetId& idALCT = (*j).first;
3182 for (CSCALCTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
3183
3184 if ((*digiIt).isValid()) {
3185 n_alcts++;
3186 histos->fill1DHist((*digiIt).getBX(), "ALCT_getBX", "ALCT.getBX()", 11, -0.5, 10.5, "TimeMonitoring");
3187 histos->fill1DHist(
3188 (*digiIt).getFullBX(), "ALCT_getFullBX", "ALCT.getFullBX()", 3601, -0.5, 3600.5, "TimeMonitoring");
3189
3190 if (ALCT_KeyWG_map.find(idALCT.chamberId()) == ALCT_KeyWG_map.end()) {
3191 ALCT_KeyWG_map[idALCT.chamberId()] = (*digiIt).getKeyWG();
3192
3193 }
3194 }
3195 }
3196 }
3197
3198
3199
3200
3201 int n_clcts = 0;
3202 map<CSCDetId, int> CLCT_getFullBx_map;
3203 for (CSCCLCTDigiCollection::DigiRangeIterator j = clcts->begin(); j != clcts->end(); j++) {
3204 const CSCCLCTDigiCollection::Range& range = (*j).second;
3205 const CSCDetId& idCLCT = (*j).first;
3206 for (CSCCLCTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
3207
3208 if ((*digiIt).isValid()) {
3209 n_clcts++;
3210 histos->fill1DHist((*digiIt).getBX(), "CLCT_getBX", "CLCT.getBX()", 11, -0.5, 10.5, "TimeMonitoring");
3211 histos->fill1DHist(
3212 (*digiIt).getFullBX(), "CLCT_getFullBX", "CLCT.getFullBX()", 3601, -0.5, 3600.5, "TimeMonitoring");
3213
3214 if (CLCT_getFullBx_map.find(idCLCT.chamberId()) == CLCT_getFullBx_map.end()) {
3215 CLCT_getFullBx_map[idCLCT.chamberId()] = (*digiIt).getFullBX();
3216
3217 }
3218 }
3219 }
3220 }
3221
3222
3223
3224
3225 int n_correlatedlcts = 0;
3226 for (CSCCorrelatedLCTDigiCollection::DigiRangeIterator j = correlatedlcts->begin(); j != correlatedlcts->end(); j++) {
3227 const CSCCorrelatedLCTDigiCollection::Range& range = (*j).second;
3228 for (CSCCorrelatedLCTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
3229 if ((*digiIt).isValid()) {
3230 n_correlatedlcts++;
3231 histos->fill1DHist(
3232 (*digiIt).getBX(), "CorrelatedLCTS_getBX", "CorrelatedLCT.getBX()", 11, -0.5, 10.5, "TimeMonitoring");
3233 }
3234 }
3235 }
3236
3237 int nRecHits = recHits->size();
3238 int nSegments = cscSegments->size();
3239 if (has_CSCTrigger) {
3240 histos->fill1DHist(L1GMT_BXN, "BX_L1CSCCand", "BX of L1 CSC Cand", 4001, -0.5, 4000.5, "TimeMonitoring");
3241 histos->fill2DHist(L1GMT_BXN,
3242 n_alcts,
3243 "n_ALCTs_v_BX_L1CSCCand",
3244 "Number of ALCTs vs. BX of L1 CSC Cand",
3245 4001,
3246 -0.5,
3247 4000.5,
3248 51,
3249 -0.5,
3250 50.5,
3251 "TimeMonitoring");
3252 histos->fill2DHist(L1GMT_BXN,
3253 n_clcts,
3254 "n_CLCTs_v_BX_L1CSCCand",
3255 "Number of CLCTs vs. BX of L1 CSC Cand",
3256 4001,
3257 -0.5,
3258 4000.5,
3259 51,
3260 -0.5,
3261 50.5,
3262 "TimeMonitoring");
3263 histos->fill2DHist(L1GMT_BXN,
3264 n_correlatedlcts,
3265 "n_CorrelatedLCTs_v_BX_L1CSCCand",
3266 "Number of CorrelatedLCTs vs. BX of L1 CSC Cand",
3267 4001,
3268 -0.5,
3269 4000.5,
3270 51,
3271 -0.5,
3272 50.5,
3273 "TimeMonitoring");
3274 histos->fill2DHist(L1GMT_BXN,
3275 nRecHits,
3276 "n_RecHits_v_BX_L1CSCCand",
3277 "Number of RecHits vs. BX of L1 CSC Cand",
3278 4001,
3279 -0.5,
3280 4000.5,
3281 101,
3282 -0.5,
3283 100.5,
3284 "TimeMonitoring");
3285 histos->fill2DHist(L1GMT_BXN,
3286 nSegments,
3287 "n_Segments_v_BX_L1CSCCand",
3288 "Number of Segments vs. BX of L1 CSC Cand",
3289 4001,
3290 -0.5,
3291 4000.5,
3292 51,
3293 -0.5,
3294 50.5,
3295 "TimeMonitoring");
3296 }
3297 if (has_CSCTrigger && has_beamHaloTrigger) {
3298 histos->fill1DHist(
3299 L1GMT_BXN, "BX_L1CSCCand_w_beamHalo", "BX of L1 CSC (w beamHalo bit)", 4001, -0.5, 4000.5, "TimeMonitoring");
3300 histos->fill2DHist(L1GMT_BXN,
3301 n_alcts,
3302 "n_ALCTs_v_BX_L1CSCCand_w_beamHalo",
3303 "Number of ALCTs vs. BX of L1 CSC Cand (w beamHalo bit)",
3304 4001,
3305 -0.5,
3306 4000.5,
3307 51,
3308 -0.5,
3309 50.5,
3310 "TimeMonitoring");
3311 histos->fill2DHist(L1GMT_BXN,
3312 n_clcts,
3313 "n_CLCTs_v_BX_L1CSCCand_w_beamHalo",
3314 "Number of CLCTs vs. BX of L1 CSC Cand (w beamHalo bit)",
3315 4001,
3316 -0.5,
3317 4000.5,
3318 51,
3319 -0.5,
3320 50.5,
3321 "TimeMonitoring");
3322 histos->fill2DHist(L1GMT_BXN,
3323 n_correlatedlcts,
3324 "n_CorrelatedLCTs_v_BX_L1CSCCand_w_beamHalo",
3325 "Number of CorrelatedLCTs vs. BX of L1 CSC Cand (w beamHalo bit)",
3326 4001,
3327 -0.5,
3328 4000.5,
3329 51,
3330 -0.5,
3331 50.5,
3332 "TimeMonitoring");
3333 histos->fill2DHist(L1GMT_BXN,
3334 nRecHits,
3335 "n_RecHits_v_BX_L1CSCCand_w_beamHalo",
3336 "Number of RecHits vs. BX of L1 CSC Cand (w beamHalo bit)",
3337 4001,
3338 -0.5,
3339 4000.5,
3340 101,
3341 -0.5,
3342 100.5,
3343 "TimeMonitoring");
3344 histos->fill2DHist(L1GMT_BXN,
3345 nSegments,
3346 "n_Segments_v_BX_L1CSCCand_w_beamHalo",
3347 "Number of Segments vs. BX of L1 CSC Cand (w beamHalo bit)",
3348 4001,
3349 -0.5,
3350 4000.5,
3351 51,
3352 -0.5,
3353 50.5,
3354 "TimeMonitoring");
3355 }
3356
3357
3358
3359
3360
3361
3362
3363 edm::ESHandle<CSCCrateMap> hcrate = eventSetup.getHandle(crateToken_);
3364 const CSCCrateMap* pcrate = hcrate.product();
3365
3366
3367 edm::Handle<FEDRawDataCollection> rawdata;
3368 event.getByToken(rd_token, rawdata);
3369
3370
3371
3372 unsigned long dccBinCheckMask = 0x06080016;
3373 unsigned int examinerMask = 0x1FEBF3F6;
3374 unsigned int errorMask = 0x0;
3375
3376 for (int id = FEDNumbering::MINCSCFEDID; id <= FEDNumbering::MAXCSCFEDID; ++id) {
3377
3378
3379
3380
3381
3382 const FEDRawData& fedData = rawdata->FEDData(id);
3383 unsigned long length = fedData.size();
3384
3385 if (length >= 32) {
3386 std::stringstream examiner_out, examiner_err;
3387
3388 CSCDCCExaminer* examiner = new CSCDCCExaminer();
3389 if (examinerMask & 0x40000)
3390 examiner->crcCFEB(true);
3391 if (examinerMask & 0x8000)
3392 examiner->crcTMB(true);
3393 if (examinerMask & 0x0400)
3394 examiner->crcALCT(true);
3395 examiner->setMask(examinerMask);
3396 const short unsigned int* data = (short unsigned int*)fedData.data();
3397
3398 bool goodEvent;
3399 if (examiner->check(data, long(fedData.size() / 2)) < 0) {
3400 goodEvent = false;
3401 } else {
3402 goodEvent = !(examiner->errors() & dccBinCheckMask);
3403 }
3404
3405 if (goodEvent) {
3406
3407 CSCDCCExaminer* ptrExaminer = examiner;
3408 CSCDCCEventData dccData((short unsigned int*)fedData.data(), ptrExaminer);
3409
3410
3411 const std::vector<CSCDDUEventData>& dduData = dccData.dduData();
3412
3413
3414 CSCDetId layer(1, 1, 1, 1, 1);
3415
3416 for (unsigned int iDDU = 0; iDDU < dduData.size(); ++iDDU) {
3417
3418
3419 if (dduData[iDDU].trailer().errorstat() & errorMask) {
3420 LogTrace("CSCDCCUnpacker|CSCRawToDigi") << "DDU# " << iDDU << " has serious error - no digis unpacked! "
3421 << std::hex << dduData[iDDU].trailer().errorstat();
3422 continue;
3423 }
3424
3425
3426 const std::vector<CSCEventData>& cscData = dduData[iDDU].cscData();
3427 for (unsigned int iCSC = 0; iCSC < cscData.size(); ++iCSC) {
3428
3429 int vmecrate = cscData[iCSC].dmbHeader()->crateID();
3430 int dmb = cscData[iCSC].dmbHeader()->dmbID();
3431
3432
3433
3434
3435 int icfeb = 0;
3436 int ilayer = 0;
3437
3438 if ((vmecrate >= 1) && (vmecrate <= 60) && (dmb >= 1) && (dmb <= 10) && (dmb != 6)) {
3439 layer = pcrate->detId(vmecrate, dmb, icfeb, ilayer);
3440 } else {
3441 LogTrace("CSCTimingAlignment|CSCDCCUnpacker|CSCRawToDigi") << " detID input out of range!!! ";
3442 LogTrace("CSCTimingAlignment|CSCDCCUnpacker|CSCRawToDigi")
3443 << " skipping chamber vme= " << vmecrate << " dmb= " << dmb;
3444 continue;
3445 }
3446
3447
3448 int nalct = cscData[iCSC].dmbHeader()->nalct();
3449 bool goodALCT = false;
3450
3451 if (nalct && cscData[iCSC].alctHeader()) {
3452 if (cscData[iCSC].alctHeader()->check()) {
3453 goodALCT = true;
3454 }
3455 }
3456
3457
3458 int nclct = cscData[iCSC].dmbHeader()->nclct();
3459 bool goodTMB = false;
3460 if (nclct && cscData[iCSC].tmbData()) {
3461 if (cscData[iCSC].tmbHeader()->check()) {
3462 if (cscData[iCSC].comparatorData()->check())
3463 goodTMB = true;
3464 }
3465 }
3466
3467 if (goodTMB && goodALCT) {
3468 if (ALCT_KeyWG_map.find(layer) == ALCT_KeyWG_map.end()) {
3469 printf("no ALCT info for Chamber %d %d %d %d \n",
3470 layer.endcap(),
3471 layer.station(),
3472 layer.ring(),
3473 layer.chamber());
3474 continue;
3475 }
3476 if (CLCT_getFullBx_map.find(layer) == CLCT_getFullBx_map.end()) {
3477 printf("no CLCT info for Chamber %d %d %d %d \n",
3478 layer.endcap(),
3479 layer.station(),
3480 layer.ring(),
3481 layer.chamber());
3482 continue;
3483 }
3484 int ALCT0Key = ALCT_KeyWG_map.find(layer)->second;
3485 int CLCTPretrigger = CLCT_getFullBx_map.find(layer)->second;
3486
3487 const CSCTMBHeader* tmbHead = cscData[iCSC].tmbHeader();
3488
3489 histos->fill1DHistByStation(tmbHead->BXNCount(),
3490 "TMB_BXNCount",
3491 "TMB_BXNCount",
3492 layer.chamberId(),
3493 3601,
3494 -0.5,
3495 3600.5,
3496 "TimeMonitoring");
3497 histos->fill1DHistByStation(tmbHead->ALCTMatchTime(),
3498 "TMB_ALCTMatchTime",
3499 "TMB_ALCTMatchTime",
3500 layer.chamberId(),
3501 7,
3502 -0.5,
3503 6.5,
3504 "TimeMonitoring");
3505
3506 histos->fill1DHist(
3507 tmbHead->BXNCount(), "TMB_BXNCount", "TMB_BXNCount", 3601, -0.5, 3600.5, "TimeMonitoring");
3508 histos->fill1DHist(
3509 tmbHead->ALCTMatchTime(), "TMB_ALCTMatchTime", "TMB_ALCTMatchTime", 7, -0.5, 6.5, "TimeMonitoring");
3510
3511 histos->fill1DHistByType(tmbHead->ALCTMatchTime(),
3512 "TMB_ALCTMatchTime",
3513 "TMB_ALCTMatchTime",
3514 layer.chamberId(),
3515 7,
3516 -0.5,
3517 6.5,
3518 "TimeMonitoring");
3519
3520 histos->fillProfile(chamberSerial(layer.chamberId()),
3521 tmbHead->ALCTMatchTime(),
3522 "prof_TMB_ALCTMatchTime",
3523 "prof_TMB_ALCTMatchTime",
3524 601,
3525 -0.5,
3526 600.5,
3527 -0.5,
3528 7.5,
3529 "TimeMonitoring");
3530 histos->fillProfile(ALCT0Key,
3531 tmbHead->ALCTMatchTime(),
3532 "prof_TMB_ALCTMatchTime_v_ALCT0KeyWG",
3533 "prof_TMB_ALCTMatchTime_v_ALCT0KeyWG",
3534 128,
3535 -0.5,
3536 127.5,
3537 0,
3538 7,
3539 "TimeMonitoring");
3540 histos->fillProfileByType(ALCT0Key,
3541 tmbHead->ALCTMatchTime(),
3542 "prf_TMB_ALCTMatchTime_v_ALCT0KeyWG",
3543 "prf_TMB_ALCTMatchTime_v_ALCT0KeyWG",
3544 layer.chamberId(),
3545 128,
3546 -0.5,
3547 127.5,
3548 0,
3549 7,
3550 "TimeMonitoring");
3551
3552
3553
3554 int TMB_ALCT_rel_L1A = tmbHead->BXNCount() - (CLCTPretrigger + 2 + tmbHead->ALCTMatchTime());
3555 if (TMB_ALCT_rel_L1A > 3563)
3556 TMB_ALCT_rel_L1A = TMB_ALCT_rel_L1A - 3564;
3557 if (TMB_ALCT_rel_L1A < 0)
3558 TMB_ALCT_rel_L1A = TMB_ALCT_rel_L1A + 3564;
3559
3560
3561 histos->fill1DHist(
3562 TMB_ALCT_rel_L1A, "h1D_TMB_ALCT_rel_L1A", "h1D_TMB_ALCT_rel_L1A", 11, 144.5, 155.5, "TimeMonitoring");
3563 histos->fill2DHist(chamberSerial(layer.chamberId()),
3564 TMB_ALCT_rel_L1A,
3565 "h2D_TMB_ALCT_rel_L1A",
3566 "h2D_TMB_ALCT_rel_L1A",
3567 601,
3568 -0.5,
3569 600.5,
3570 11,
3571 144.5,
3572 155.5,
3573 "TimeMonitoring");
3574 histos->fill2DHist(ringSerial(layer.chamberId()),
3575 TMB_ALCT_rel_L1A,
3576 "h2D_TMB_ALCT_rel_L1A_by_ring",
3577 "h2D_TMB_ALCT_rel_L1A_by_ring",
3578 19,
3579 -9.5,
3580 9.5,
3581 11,
3582 144.5,
3583 155.5,
3584 "TimeMonitoring");
3585 histos->fillProfile(chamberSerial(layer.chamberId()),
3586 TMB_ALCT_rel_L1A,
3587 "prof_TMB_ALCT_rel_L1A",
3588 "prof_TMB_ALCT_rel_L1A",
3589 601,
3590 -0.5,
3591 600.5,
3592 145,
3593 155,
3594 "TimeMonitoring");
3595 histos->fillProfile(ringSerial(layer.chamberId()),
3596 TMB_ALCT_rel_L1A,
3597 "prof_TMB_ALCT_rel_L1A_by_ring",
3598 "prof_TMB_ALCT_rel_L1A_by_ring",
3599 19,
3600 -9.5,
3601 9.5,
3602 145,
3603 155,
3604 "TimeMonitoring");
3605
3606 histos->fill2DHist(ALCT0Key,
3607 TMB_ALCT_rel_L1A,
3608 "h2D_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3609 "h2D_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3610 128,
3611 -0.5,
3612 127.5,
3613 11,
3614 144.5,
3615 155.5,
3616 "TimeMonitoring");
3617 histos->fillProfile(ALCT0Key,
3618 TMB_ALCT_rel_L1A,
3619 "prof_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3620 "prof_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3621 128,
3622 -0.5,
3623 127.5,
3624 145,
3625 155,
3626 "TimeMonitoring");
3627 histos->fillProfileByType(ALCT0Key,
3628 TMB_ALCT_rel_L1A,
3629 "prf_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3630 "prf_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3631 layer.chamberId(),
3632 128,
3633 -0.5,
3634 127.5,
3635 145,
3636 155,
3637 "TimeMonitoring");
3638 }
3639
3640 }
3641 }
3642 }
3643 if (examiner != nullptr)
3644 delete examiner;
3645 }
3646 }
3647 }
3648
3649 void CSCValidation::beginJob() { std::cout << "CSCValidation starting..." << std::endl; }
3650
3651 void CSCValidation::endJob() { std::cout << "CSCValidation: Events analyzed " << nEventsAnalyzed << std::endl; }
3652
3653 DEFINE_FWK_MODULE(CSCValidation);