CSCValidation/src/CSCValidation.cc

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