Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:24:47

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