CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
CSCValidation.cc
Go to the documentation of this file.
1 /*
2  * validation package for CSC DIGIs, RECHITs and SEGMENTs + more.
3  *
4  * Michael Schmitt
5  * Andy Kubik
6  * Northwestern University
7  */
8 #include "RecoLocalMuon/CSCValidation/src/CSCValidation.h"
9 
10 using namespace std;
11 using namespace edm;
12 
14 // CONSTRUCTOR //
16 CSCValidation::CSCValidation(const ParameterSet& pset) {
17  // Get the various input parameters
18  rootFileName = pset.getUntrackedParameter<std::string>("rootFileName", "valHists.root");
19  isSimulation = pset.getUntrackedParameter<bool>("isSimulation", false);
20  writeTreeToFile = pset.getUntrackedParameter<bool>("writeTreeToFile", true);
21  detailedAnalysis = pset.getUntrackedParameter<bool>("detailedAnalysis", false);
22  useDigis = pset.getUntrackedParameter<bool>("useDigis", true);
23 
24  // event quality filter
25  useQualityFilter = pset.getUntrackedParameter<bool>("useQualityFilter", false);
26  pMin = pset.getUntrackedParameter<double>("pMin", 4.);
27  chisqMax = pset.getUntrackedParameter<double>("chisqMax", 20.);
28  nCSCHitsMin = pset.getUntrackedParameter<int>("nCSCHitsMin", 10);
29  nCSCHitsMax = pset.getUntrackedParameter<int>("nCSCHitsMax", 25);
30  lengthMin = pset.getUntrackedParameter<double>("lengthMin", 140.);
31  lengthMax = pset.getUntrackedParameter<double>("lengthMax", 600.);
32  deltaPhiMax = pset.getUntrackedParameter<double>("deltaPhiMax", 0.2);
33  polarMax = pset.getUntrackedParameter<double>("polarMax", 0.7);
34  polarMin = pset.getUntrackedParameter<double>("polarMin", 0.3);
35 
36  // trigger filter
37  useTriggerFilter = pset.getUntrackedParameter<bool>("useTriggerFilter", false);
38 
39  // input tags for collections
40  rd_token = consumes<FEDRawDataCollection>(pset.getParameter<edm::InputTag>("rawDataTag"));
41  sd_token = consumes<CSCStripDigiCollection>(pset.getParameter<edm::InputTag>("stripDigiTag"));
42  wd_token = consumes<CSCWireDigiCollection>(pset.getParameter<edm::InputTag>("wireDigiTag"));
43  cd_token = consumes<CSCComparatorDigiCollection>(pset.getParameter<edm::InputTag>("compDigiTag"));
44  al_token = consumes<CSCALCTDigiCollection>(pset.getParameter<edm::InputTag>("alctDigiTag"));
45  cl_token = consumes<CSCCLCTDigiCollection>(pset.getParameter<edm::InputTag>("clctDigiTag"));
46  co_token = consumes<CSCCorrelatedLCTDigiCollection>(pset.getParameter<edm::InputTag>("corrlctDigiTag"));
47  rh_token = consumes<CSCRecHit2DCollection>(pset.getParameter<edm::InputTag>("cscRecHitTag"));
48  se_token = consumes<CSCSegmentCollection>(pset.getParameter<edm::InputTag>("cscSegTag"));
49  sa_token = consumes<reco::TrackCollection>(pset.getParameter<edm::InputTag>("saMuonTag"));
50  l1_token = consumes<L1MuGMTReadoutCollection>(pset.getParameter<edm::InputTag>("l1aTag"));
51  tr_token = consumes<TriggerResults>(pset.getParameter<edm::InputTag>("hltTag"));
52  sh_token = consumes<PSimHitContainer>(pset.getParameter<edm::InputTag>("simHitTag"));
53 
54  // flags to switch on/off individual modules
55  makeOccupancyPlots = pset.getUntrackedParameter<bool>("makeOccupancyPlots", true);
56  makeTriggerPlots = pset.getUntrackedParameter<bool>("makeTriggerPlots", false);
57  makeStripPlots = pset.getUntrackedParameter<bool>("makeStripPlots", true);
58  makeWirePlots = pset.getUntrackedParameter<bool>("makeWirePlots", true);
59  makeRecHitPlots = pset.getUntrackedParameter<bool>("makeRecHitPlots", true);
60  makeSimHitPlots = pset.getUntrackedParameter<bool>("makeSimHitPlots", true);
61  makeSegmentPlots = pset.getUntrackedParameter<bool>("makeSegmentPlots", true);
62  makeResolutionPlots = pset.getUntrackedParameter<bool>("makeResolutionPlots", true);
63  makePedNoisePlots = pset.getUntrackedParameter<bool>("makePedNoisePlots", true);
64  makeEfficiencyPlots = pset.getUntrackedParameter<bool>("makeEfficiencyPlots", true);
65  makeGasGainPlots = pset.getUntrackedParameter<bool>("makeGasGainPlots", true);
66  makeAFEBTimingPlots = pset.getUntrackedParameter<bool>("makeAFEBTimingPlots", true);
67  makeCompTimingPlots = pset.getUntrackedParameter<bool>("makeCompTimingPlots", true);
68  makeADCTimingPlots = pset.getUntrackedParameter<bool>("makeADCTimingPlots", true);
69  makeRHNoisePlots = pset.getUntrackedParameter<bool>("makeRHNoisePlots", false);
70  makeCalibPlots = pset.getUntrackedParameter<bool>("makeCalibPlots", false);
71  makeStandalonePlots = pset.getUntrackedParameter<bool>("makeStandalonePlots", false);
72  makeTimeMonitorPlots = pset.getUntrackedParameter<bool>("makeTimeMonitorPlots", false);
73  makeHLTPlots = pset.getUntrackedParameter<bool>("makeHLTPlots", false);
74 
75  // set counters to zero
76  nEventsAnalyzed = 0;
77  rhTreeCount = 0;
78  segTreeCount = 0;
79  firstEvent = true;
80 
81  // Create the root file for the histograms
82  theFile = new TFile(rootFileName.c_str(), "RECREATE");
83  theFile->cd();
84 
85  // Create object of class CSCValHists to manage histograms
86  histos = new CSCValHists();
87 
88  // book Occupancy Histos
89  hOWires = new TH2I("hOWires", "Wire Digi Occupancy", 36, 0.5, 36.5, 20, 0.5, 20.5);
90  hOStrips = new TH2I("hOStrips", "Strip Digi Occupancy", 36, 0.5, 36.5, 20, 0.5, 20.5);
91  hORecHits = new TH2I("hORecHits", "RecHit Occupancy", 36, 0.5, 36.5, 20, 0.5, 20.5);
92  hOSegments = new TH2I("hOSegments", "Segments Occupancy", 36, 0.5, 36.5, 20, 0.5, 20.5);
93 
94  // book Eff histos
95  hSSTE = new TH1F("hSSTE", "hSSTE", 40, 0, 40);
96  hRHSTE = new TH1F("hRHSTE", "hRHSTE", 40, 0, 40);
97  hSEff = new TH1F("hSEff", "Segment Efficiency", 20, 0.5, 20.5);
98  hRHEff = new TH1F("hRHEff", "recHit Efficiency", 20, 0.5, 20.5);
99 
100  const int nChambers = 36;
101  const int nTypes = 18;
102  float nCH_min = 0.5;
103  float nCh_max = float(nChambers) + 0.5;
104  float nT_min = 0.5;
105  float nT_max = float(nTypes) + 0.5;
106 
107  hSSTE2 = new TH2F("hSSTE2", "hSSTE2", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
108  hRHSTE2 = new TH2F("hRHSTE2", "hRHSTE2", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
109  hStripSTE2 = new TH2F("hStripSTE2", "hStripSTE2", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
110  hWireSTE2 = new TH2F("hWireSTE2", "hWireSTE2", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
111 
112  hEffDenominator = new TH2F("hEffDenominator", "hEffDenominator", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
113  hSEff2 = new TH2F("hSEff2", "Segment Efficiency 2D", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
114  hRHEff2 = new TH2F("hRHEff2", "recHit Efficiency 2D", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
115 
116  hStripEff2 = new TH2F("hStripEff2", "strip Efficiency 2D", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
117  hWireEff2 = new TH2F("hWireEff2", "wire Efficiency 2D", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
118 
119  hSensitiveAreaEvt =
120  new TH2F("hSensitiveAreaEvt", "events in sensitive area", nChambers, nCH_min, nCh_max, nTypes, nT_min, nT_max);
121 
122  // setup trees to hold global position data for rechits and segments
123  if (writeTreeToFile)
124  histos->setupTrees();
125 
126  geomToken_ = esConsumes<CSCGeometry, MuonGeometryRecord>();
127 }
128 
130 // DESTRUCTOR //
132 CSCValidation::~CSCValidation() {
133  // produce final efficiency histograms
134  histoEfficiency(hRHSTE, hRHEff);
135  histoEfficiency(hSSTE, hSEff);
136  hSEff2->Divide(hSSTE2, hEffDenominator, 1., 1., "B");
137  hRHEff2->Divide(hRHSTE2, hEffDenominator, 1., 1., "B");
138  hStripEff2->Divide(hStripSTE2, hEffDenominator, 1., 1., "B");
139  hWireEff2->Divide(hWireSTE2, hEffDenominator, 1., 1., "B");
140 
141  histos->insertPlot(hRHSTE, "hRHSTE", "Efficiency");
142  histos->insertPlot(hSSTE, "hSSTE", "Efficiency");
143  histos->insertPlot(hSSTE2, "hSSTE2", "Efficiency");
144  histos->insertPlot(hEffDenominator, "hEffDenominator", "Efficiency");
145  histos->insertPlot(hRHSTE2, "hRHSTE2", "Efficiency");
146  histos->insertPlot(hStripSTE2, "hStripSTE2", "Efficiency");
147  histos->insertPlot(hWireSTE2, "hWireSTE2", "Efficiency");
148 
149  //moving this to post job macros
150  histos->insertPlot(hSEff, "hSEff", "Efficiency");
151  histos->insertPlot(hRHEff, "hRHEff", "Efficiency");
152 
153  histos->insertPlot(hSEff2, "hSEff2", "Efficiency");
154  histos->insertPlot(hRHEff2, "hRHEff2", "Efficiency");
155  histos->insertPlot(hStripEff2, "hStripff2", "Efficiency");
156  histos->insertPlot(hWireEff2, "hWireff2", "Efficiency");
157 
158  histos->insertPlot(hSensitiveAreaEvt, "", "Efficiency");
159 
160  // throw in occupancy plots so they're saved
161  histos->insertPlot(hOWires, "hOWires", "Digis");
162  histos->insertPlot(hOStrips, "hOStrips", "Digis");
163  histos->insertPlot(hORecHits, "hORecHits", "recHits");
164  histos->insertPlot(hOSegments, "hOSegments", "Segments");
165 
166  // write histos to the specified file
167  histos->writeHists(theFile);
168  if (writeTreeToFile)
169  histos->writeTrees(theFile);
170  theFile->Close();
171 }
172 
174 // Analysis //
176 void CSCValidation::analyze(const Event& event, const EventSetup& eventSetup) {
177  // increment counter
178  nEventsAnalyzed++;
179 
180  //int iRun = event.id().run();
181  //int iEvent = event.id().event();
182 
183  // Get the Digis
184  edm::Handle<CSCWireDigiCollection> wires;
185  edm::Handle<CSCStripDigiCollection> strips;
186  edm::Handle<CSCComparatorDigiCollection> compars;
187  edm::Handle<CSCALCTDigiCollection> alcts;
188  edm::Handle<CSCCLCTDigiCollection> clcts;
189  edm::Handle<CSCCorrelatedLCTDigiCollection> correlatedlcts;
190  if (useDigis) {
191  event.getByToken(sd_token, strips);
192  event.getByToken(wd_token, wires);
193  event.getByToken(cd_token, compars);
194  event.getByToken(al_token, alcts);
195  event.getByToken(cl_token, clcts);
196  event.getByToken(co_token, correlatedlcts);
197  }
198 
199  // Get the CSC Geometry :
200  edm::ESHandle<CSCGeometry> cscGeom = eventSetup.getHandle(geomToken_);
201 
202  // Get the RecHits collection :
203  edm::Handle<CSCRecHit2DCollection> recHits;
204  event.getByToken(rh_token, recHits);
205 
206  //CSCRecHit2DCollection::const_iterator recIt;
207  //for (recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
208  // recIt->print();
209  // }
210 
211  // Get the SimHits (if applicable)
212  edm::Handle<PSimHitContainer> simHits;
213  if (isSimulation)
214  event.getByToken(sh_token, simHits);
215 
216  // get CSC segment collection
217  edm::Handle<CSCSegmentCollection> cscSegments;
218  event.getByToken(se_token, cscSegments);
219 
220  // get the trigger collection
221  edm::Handle<L1MuGMTReadoutCollection> pCollection;
222  if (makeTriggerPlots || useTriggerFilter || (useDigis && makeTimeMonitorPlots)) {
223  event.getByToken(l1_token, pCollection);
224  }
225  edm::Handle<TriggerResults> hlt;
226  if (makeHLTPlots) {
227  event.getByToken(tr_token, hlt);
228  }
229 
230  // get the standalone muon collection
231  edm::Handle<reco::TrackCollection> saMuons;
232  if (makeStandalonePlots || useQualityFilter) {
233  event.getByToken(sa_token, saMuons);
234  }
235 
237  // Run the modules //
239 
240  // Only do this for the first event
241  // this is probably outdated and needs to be looked at
242  if (nEventsAnalyzed == 1 && makeCalibPlots)
243  doCalibrations(eventSetup);
244 
245  // Look at the l1a trigger info (returns true if csc L1A present)
246  bool CSCL1A = false;
247  if (makeTriggerPlots || useTriggerFilter)
248  CSCL1A = doTrigger(pCollection);
249  if (!useTriggerFilter)
250  CSCL1A = true; // always true if not filtering on trigger
251 
252  cleanEvent = false;
253  if (makeStandalonePlots || useQualityFilter)
254  cleanEvent = filterEvents(recHits, cscSegments, saMuons);
255  if (!useQualityFilter)
256  cleanEvent = true; // always true if not filtering on event quality
257 
258  // look at various chamber occupancies
259  // keep this outside of filter for diagnostics???
260  if (makeOccupancyPlots && CSCL1A)
261  doOccupancies(strips, wires, recHits, cscSegments);
262 
263  if (makeHLTPlots)
264  doHLT(hlt);
265 
266  if (cleanEvent && CSCL1A) {
267  // general look at strip digis
268  if (makeStripPlots && useDigis)
269  doStripDigis(strips);
270 
271  // general look at wire digis
272  if (makeWirePlots && useDigis)
273  doWireDigis(wires);
274 
275  // general look at rechits
276  if (makeRecHitPlots)
277  doRecHits(recHits, cscGeom);
278 
279  // look at simHits
280  if (isSimulation && makeSimHitPlots)
281  doSimHits(recHits, simHits);
282 
283  // general look at Segments
284  if (makeSegmentPlots)
285  doSegments(cscSegments, cscGeom);
286 
287  // look at hit resolution
288  if (makeResolutionPlots)
289  doResolution(cscSegments, cscGeom);
290 
291  // look at Pedestal Noise
292  if (makePedNoisePlots && useDigis)
293  doPedestalNoise(strips);
294 
295  // look at recHit and segment efficiencies
296  if (makeEfficiencyPlots)
297  doEfficiencies(wires, strips, recHits, cscSegments, cscGeom);
298 
299  // gas gain
300  if (makeGasGainPlots && useDigis)
301  doGasGain(*wires, *strips, *recHits);
302 
303  // AFEB timing
304  if (makeAFEBTimingPlots && useDigis)
305  doAFEBTiming(*wires);
306 
307  // Comparators timing
308  if (makeCompTimingPlots && useDigis)
309  doCompTiming(*compars);
310 
311  // strip ADC timing
312  if (makeADCTimingPlots)
313  doADCTiming(*recHits);
314 
315  // recHit Noise
316  if (makeRHNoisePlots && useDigis)
317  doNoiseHits(recHits, cscSegments, cscGeom, strips);
318 
319  // look at standalone muons (not implemented yet)
320  if (makeStandalonePlots)
321  doStandalone(saMuons);
322 
323  // make plots for monitoring the trigger and offline timing
324  if (makeTimeMonitorPlots)
325  doTimeMonitoring(recHits, cscSegments, alcts, clcts, correlatedlcts, pCollection, cscGeom, eventSetup, event);
326 
327  firstEvent = false;
328  }
329 }
330 
331 // ==============================================
332 //
333 // event filter, returns true only for events with "good" standalone muon
334 //
335 // ==============================================
336 
337 bool CSCValidation::filterEvents(edm::Handle<CSCRecHit2DCollection> recHits,
338  edm::Handle<CSCSegmentCollection> cscSegments,
339  edm::Handle<reco::TrackCollection> saMuons) {
340  //int nGoodSAMuons = 0;
341 
342  if (recHits->size() < 4 || recHits->size() > 100)
343  return false;
344  if (cscSegments->size() < 1 || cscSegments->size() > 15)
345  return false;
346  return true;
347  //if (saMuons->size() != 1) return false;
348  /*
349  for(reco::TrackCollection::const_iterator muon = saMuons->begin(); muon != saMuons->end(); ++ muon ) {
350  double p = muon->p();
351  double reducedChisq = muon->normalizedChi2();
352 
353  GlobalPoint innerPnt(muon->innerPosition().x(),muon->innerPosition().y(),muon->innerPosition().z());
354  GlobalPoint outerPnt(muon->outerPosition().x(),muon->outerPosition().y(),muon->outerPosition().z());
355  GlobalVector innerKin(muon->innerMomentum().x(),muon->innerMomentum().y(),muon->innerMomentum().z());
356  GlobalVector outerKin(muon->outerMomentum().x(),muon->outerMomentum().y(),muon->outerMomentum().z());
357  GlobalVector deltaPnt = innerPnt - outerPnt;
358  double crudeLength = deltaPnt.mag();
359  double deltaPhi = innerPnt.phi() - outerPnt.phi();
360  double innerGlobalPolarAngle = innerKin.theta();
361  double outerGlobalPolarAngle = outerKin.theta();
362 
363  int nCSCHits = 0;
364  for (trackingRecHit_iterator hit = muon->recHitsBegin(); hit != muon->recHitsEnd(); ++hit ) {
365  if ( (*hit)->isValid() ) {
366  const DetId detId( (*hit)->geographicalId() );
367  if (detId.det() == DetId::Muon) {
368  if (detId.subdetId() == MuonSubdetId::CSC) {
369  nCSCHits++;
370  } // this is a CSC hit
371  } // this is a muon hit
372  } // hit is valid
373  } // end loop over rechits
374 
375  bool goodSAMuon = (p > pMin)
376  && ( reducedChisq < chisqMax )
377  && ( nCSCHits >= nCSCHitsMin )
378  && ( nCSCHits <= nCSCHitsMax )
379  && ( crudeLength > lengthMin )
380  && ( crudeLength < lengthMax );
381 
382 
383  goodSAMuon = goodSAMuon && ( fabs(deltaPhi) < deltaPhiMax );
384  goodSAMuon = goodSAMuon &&
385  (
386  ( ( innerGlobalPolarAngle > polarMin) && ( innerGlobalPolarAngle < polarMax) ) ||
387  ( (M_PI-innerGlobalPolarAngle > polarMin) && (M_PI-innerGlobalPolarAngle < polarMax) )
388  );
389  goodSAMuon = goodSAMuon &&
390  (
391  ( ( outerGlobalPolarAngle > polarMin) && ( outerGlobalPolarAngle < polarMax) ) ||
392  ( (M_PI-outerGlobalPolarAngle > polarMin) && (M_PI-outerGlobalPolarAngle < polarMax) )
393  );
394 
395  //goodSAMuon = goodSAMuon && (nCSCHits > nCSCHitsMin) && (nCSCHits < 13);
396  //goodSAMuon = goodSAMuon && (nCSCHits > 13) && (nCSCHits < 19);
397  //goodSAMuon = goodSAMuon && (nCSCHits > 19) && (nCSCHits < nCSCHitsMax);
398 
399 
400  if (goodSAMuon) nGoodSAMuons++;
401 
402  } // end loop over stand-alone muon collection
403 
404 
405  histos->fill1DHist(nGoodSAMuons,"hNGoodMuons", "Number of Good STA Muons per Event",11,-0.5,10.5,"STAMuons");
406 
407  if (nGoodSAMuons == 1) return true;
408  return false;
409  */
410 }
411 
412 // ==============================================
413 //
414 // look at Occupancies
415 //
416 // ==============================================
417 
418 void CSCValidation::doOccupancies(edm::Handle<CSCStripDigiCollection> strips,
419  edm::Handle<CSCWireDigiCollection> wires,
420  edm::Handle<CSCRecHit2DCollection> recHits,
421  edm::Handle<CSCSegmentCollection> cscSegments) {
422  bool wireo[2][4][4][36];
423  bool stripo[2][4][4][36];
424  bool rechito[2][4][4][36];
425  bool segmento[2][4][4][36];
426 
427  bool hasWires = false;
428  bool hasStrips = false;
429  bool hasRecHits = false;
430  bool hasSegments = false;
431 
432  for (int e = 0; e < 2; e++) {
433  for (int s = 0; s < 4; s++) {
434  for (int r = 0; r < 4; r++) {
435  for (int c = 0; c < 36; c++) {
436  wireo[e][s][r][c] = false;
437  stripo[e][s][r][c] = false;
438  rechito[e][s][r][c] = false;
439  segmento[e][s][r][c] = false;
440  }
441  }
442  }
443  }
444 
445  if (useDigis) {
446  //wires
447  for (CSCWireDigiCollection::DigiRangeIterator wi = wires->begin(); wi != wires->end(); wi++) {
448  CSCDetId id = (CSCDetId)(*wi).first;
449  int kEndcap = id.endcap();
450  int kRing = id.ring();
451  int kStation = id.station();
452  int kChamber = id.chamber();
453  std::vector<CSCWireDigi>::const_iterator wireIt = (*wi).second.first;
454  std::vector<CSCWireDigi>::const_iterator lastWire = (*wi).second.second;
455  for (; wireIt != lastWire; ++wireIt) {
456  if (!wireo[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1]) {
457  wireo[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1] = true;
458  hOWires->Fill(kChamber, typeIndex(id));
459  histos->fill1DHist(
460  chamberSerial(id), "hOWireSerial", "Wire Occupancy by Chamber Serial", 601, -0.5, 600.5, "Digis");
461  hasWires = true;
462  }
463  }
464  }
465 
466  //strips
467  for (CSCStripDigiCollection::DigiRangeIterator si = strips->begin(); si != strips->end(); si++) {
468  CSCDetId id = (CSCDetId)(*si).first;
469  int kEndcap = id.endcap();
470  int kRing = id.ring();
471  int kStation = id.station();
472  int kChamber = id.chamber();
473  std::vector<CSCStripDigi>::const_iterator stripIt = (*si).second.first;
474  std::vector<CSCStripDigi>::const_iterator lastStrip = (*si).second.second;
475  for (; stripIt != lastStrip; ++stripIt) {
476  std::vector<int> myADCVals = stripIt->getADCCounts();
477  bool thisStripFired = false;
478  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
479  float threshold = 13.3;
480  float diff = 0.;
481  for (unsigned int iCount = 0; iCount < myADCVals.size(); iCount++) {
482  diff = (float)myADCVals[iCount] - thisPedestal;
483  if (diff > threshold) {
484  thisStripFired = true;
485  }
486  }
487  if (thisStripFired) {
488  if (!stripo[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1]) {
489  stripo[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1] = true;
490  hOStrips->Fill(kChamber, typeIndex(id));
491  histos->fill1DHist(
492  chamberSerial(id), "hOStripSerial", "Strip Occupancy by Chamber Serial", 601, -0.5, 600.5, "Digis");
493  hasStrips = true;
494  }
495  }
496  }
497  }
498  }
499 
500  //rechits
501  CSCRecHit2DCollection::const_iterator recIt;
502  for (recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
503  CSCDetId idrec = (CSCDetId)(*recIt).cscDetId();
504  int kEndcap = idrec.endcap();
505  int kRing = idrec.ring();
506  int kStation = idrec.station();
507  int kChamber = idrec.chamber();
508  if (!rechito[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1]) {
509  rechito[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1] = true;
510  histos->fill1DHist(
511  chamberSerial(idrec), "hORecHitsSerial", "RecHit Occupancy by Chamber Serial", 601, -0.5, 600.5, "recHits");
512  hORecHits->Fill(kChamber, typeIndex(idrec));
513  hasRecHits = true;
514  }
515  }
516 
517  //segments
518  for (CSCSegmentCollection::const_iterator segIt = cscSegments->begin(); segIt != cscSegments->end(); segIt++) {
519  CSCDetId id = (CSCDetId)(*segIt).cscDetId();
520  int kEndcap = id.endcap();
521  int kRing = id.ring();
522  int kStation = id.station();
523  int kChamber = id.chamber();
524  if (!segmento[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1]) {
525  segmento[kEndcap - 1][kStation - 1][kRing - 1][kChamber - 1] = true;
526  histos->fill1DHist(
527  chamberSerial(id), "hOSegmentsSerial", "Segment Occupancy by Chamber Serial", 601, -0.5, 600.5, "Segments");
528  hOSegments->Fill(kChamber, typeIndex(id));
529  hasSegments = true;
530  }
531  }
532 
533  // overall CSC occupancy (events with CSC data compared to total)
534  histos->fill1DHist(1, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
535  if (hasWires)
536  histos->fill1DHist(3, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
537  if (hasStrips)
538  histos->fill1DHist(5, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
539  if (hasWires && hasStrips)
540  histos->fill1DHist(7, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
541  if (hasRecHits)
542  histos->fill1DHist(9, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
543  if (hasSegments)
544  histos->fill1DHist(11, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
545  if (!cleanEvent)
546  histos->fill1DHist(13, "hCSCOccupancy", "overall CSC occupancy", 15, -0.5, 14.5, "GeneralHists");
547 }
548 
549 // ==============================================
550 //
551 // look at Trigger info
552 //
553 // ==============================================
554 
555 bool CSCValidation::doTrigger(edm::Handle<L1MuGMTReadoutCollection> pCollection) {
556  std::vector<L1MuGMTReadoutRecord> L1Mrec = pCollection->getRecords();
557  std::vector<L1MuGMTReadoutRecord>::const_iterator igmtrr;
558 
559  bool csc_l1a = false;
560  bool dt_l1a = false;
561  bool rpcf_l1a = false;
562  bool rpcb_l1a = false;
563  bool beamHaloTrigger = false;
564 
565  int myBXNumber = -1000;
566 
567  for (igmtrr = L1Mrec.begin(); igmtrr != L1Mrec.end(); igmtrr++) {
568  std::vector<L1MuRegionalCand>::const_iterator iter1;
569  std::vector<L1MuRegionalCand> rmc;
570 
571  // CSC
572  int icsc = 0;
573  rmc = igmtrr->getCSCCands();
574  for (iter1 = rmc.begin(); iter1 != rmc.end(); iter1++) {
575  if (!(*iter1).empty()) {
576  icsc++;
577  int kQuality = (*iter1).quality(); // kQuality = 1 means beam halo
578  if (kQuality == 1)
579  beamHaloTrigger = true;
580  }
581  }
582  if (igmtrr->getBxInEvent() == 0 && icsc > 0)
583  csc_l1a = true;
584  if (igmtrr->getBxInEvent() == 0) {
585  myBXNumber = igmtrr->getBxNr();
586  }
587 
588  // DT
589  int idt = 0;
590  rmc = igmtrr->getDTBXCands();
591  for (iter1 = rmc.begin(); iter1 != rmc.end(); iter1++) {
592  if (!(*iter1).empty()) {
593  idt++;
594  }
595  }
596  if (igmtrr->getBxInEvent() == 0 && idt > 0)
597  dt_l1a = true;
598 
599  // RPC Barrel
600  int irpcb = 0;
601  rmc = igmtrr->getBrlRPCCands();
602  for (iter1 = rmc.begin(); iter1 != rmc.end(); iter1++) {
603  if (!(*iter1).empty()) {
604  irpcb++;
605  }
606  }
607  if (igmtrr->getBxInEvent() == 0 && irpcb > 0)
608  rpcb_l1a = true;
609 
610  // RPC Forward
611  int irpcf = 0;
612  rmc = igmtrr->getFwdRPCCands();
613  for (iter1 = rmc.begin(); iter1 != rmc.end(); iter1++) {
614  if (!(*iter1).empty()) {
615  irpcf++;
616  }
617  }
618  if (igmtrr->getBxInEvent() == 0 && irpcf > 0)
619  rpcf_l1a = true;
620  }
621 
622  // Fill some histograms with L1A info
623  if (csc_l1a)
624  histos->fill1DHist(myBXNumber, "vtBXNumber", "BX Number", 4001, -0.5, 4000.5, "Trigger");
625  if (csc_l1a)
626  histos->fill1DHist(1, "vtBits", "trigger bits", 11, -0.5, 10.5, "Trigger");
627  if (dt_l1a)
628  histos->fill1DHist(2, "vtBits", "trigger bits", 11, -0.5, 10.5, "Trigger");
629  if (rpcb_l1a)
630  histos->fill1DHist(3, "vtBits", "trigger bits", 11, -0.5, 10.5, "Trigger");
631  if (rpcf_l1a)
632  histos->fill1DHist(4, "vtBits", "trigger bits", 11, -0.5, 10.5, "Trigger");
633  if (beamHaloTrigger)
634  histos->fill1DHist(8, "vtBits", "trigger bits", 11, -0.5, 10.5, "Trigger");
635 
636  if (csc_l1a) {
637  histos->fill1DHist(1, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
638  if (dt_l1a)
639  histos->fill1DHist(2, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
640  if (rpcb_l1a)
641  histos->fill1DHist(3, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
642  if (rpcf_l1a)
643  histos->fill1DHist(4, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
644  if (dt_l1a || rpcb_l1a || rpcf_l1a)
645  histos->fill1DHist(5, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
646  if (!(dt_l1a || rpcb_l1a || rpcf_l1a))
647  histos->fill1DHist(6, "vtCSCY", "trigger bits", 11, -0.5, 10.5, "Trigger");
648  } else {
649  histos->fill1DHist(1, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
650  if (dt_l1a)
651  histos->fill1DHist(2, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
652  if (rpcb_l1a)
653  histos->fill1DHist(3, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
654  if (rpcf_l1a)
655  histos->fill1DHist(4, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
656  if (dt_l1a || rpcb_l1a || rpcf_l1a)
657  histos->fill1DHist(5, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
658  if (!(dt_l1a || rpcb_l1a || rpcf_l1a))
659  histos->fill1DHist(6, "vtCSCN", "trigger bits", 11, -0.5, 10.5, "Trigger");
660  }
661 
662  // if valid CSC L1A then return true for possible use elsewhere
663 
664  if (csc_l1a)
665  return true;
666 
667  return false;
668 }
669 
670 // ==============================================
671 //
672 // look at HLT Trigger
673 //
674 // ==============================================
675 
676 bool CSCValidation::doHLT(Handle<TriggerResults> hlt) {
677  // HLT stuff
678  int hltSize = hlt->size();
679  for (int i = 0; i < hltSize; ++i) {
680  if (hlt->accept(i))
681  histos->fill1DHist(i, "hltBits", "HLT Trigger Bits", hltSize + 1, -0.5, (float)hltSize + 0.5, "Trigger");
682  }
683 
684  return true;
685 }
686 
687 // ==============================================
688 //
689 // look at Calibrations
690 //
691 // ==============================================
692 
693 void CSCValidation::doCalibrations(const edm::EventSetup& eventSetup) {
694  // Only do this for the first event
695  if (nEventsAnalyzed == 1) {
696  LogDebug("Calibrations") << "Loading Calibrations...";
697 
698  // get the gains
699  edm::ESHandle<CSCDBGains> hGains;
700  eventSetup.get<CSCDBGainsRcd>().get(hGains);
701  const CSCDBGains* pGains = hGains.product();
702  // get the crosstalks
703  edm::ESHandle<CSCDBCrosstalk> hCrosstalk;
704  eventSetup.get<CSCDBCrosstalkRcd>().get(hCrosstalk);
705  const CSCDBCrosstalk* pCrosstalk = hCrosstalk.product();
706  // get the noise matrix
707  edm::ESHandle<CSCDBNoiseMatrix> hNoiseMatrix;
708  eventSetup.get<CSCDBNoiseMatrixRcd>().get(hNoiseMatrix);
709  const CSCDBNoiseMatrix* pNoiseMatrix = hNoiseMatrix.product();
710  // get pedestals
711  edm::ESHandle<CSCDBPedestals> hPedestals;
712  eventSetup.get<CSCDBPedestalsRcd>().get(hPedestals);
713  const CSCDBPedestals* pPedestals = hPedestals.product();
714 
715  LogDebug("Calibrations") << "Calibrations Loaded!";
716 
717  for (int i = 0; i < 400; i++) {
718  int bin = i + 1;
719  histos->fillCalibHist(pGains->gains[i].gain_slope, "hCalibGainsS", "Gains Slope", 400, 0, 400, bin, "Calib");
720  histos->fillCalibHist(
721  pCrosstalk->crosstalk[i].xtalk_slope_left, "hCalibXtalkSL", "Xtalk Slope Left", 400, 0, 400, bin, "Calib");
722  histos->fillCalibHist(
723  pCrosstalk->crosstalk[i].xtalk_slope_right, "hCalibXtalkSR", "Xtalk Slope Right", 400, 0, 400, bin, "Calib");
724  histos->fillCalibHist(pCrosstalk->crosstalk[i].xtalk_intercept_left,
725  "hCalibXtalkIL",
726  "Xtalk Intercept Left",
727  400,
728  0,
729  400,
730  bin,
731  "Calib");
732  histos->fillCalibHist(pCrosstalk->crosstalk[i].xtalk_intercept_right,
733  "hCalibXtalkIR",
734  "Xtalk Intercept Right",
735  400,
736  0,
737  400,
738  bin,
739  "Calib");
740  histos->fillCalibHist(pPedestals->pedestals[i].ped, "hCalibPedsP", "Peds", 400, 0, 400, bin, "Calib");
741  histos->fillCalibHist(pPedestals->pedestals[i].rms, "hCalibPedsR", "Peds RMS", 400, 0, 400, bin, "Calib");
742  histos->fillCalibHist(
743  pNoiseMatrix->matrix[i].elem33, "hCalibNoise33", "Noise Matrix 33", 400, 0, 400, bin, "Calib");
744  histos->fillCalibHist(
745  pNoiseMatrix->matrix[i].elem34, "hCalibNoise34", "Noise Matrix 34", 400, 0, 400, bin, "Calib");
746  histos->fillCalibHist(
747  pNoiseMatrix->matrix[i].elem35, "hCalibNoise35", "Noise Matrix 35", 400, 0, 400, bin, "Calib");
748  histos->fillCalibHist(
749  pNoiseMatrix->matrix[i].elem44, "hCalibNoise44", "Noise Matrix 44", 400, 0, 400, bin, "Calib");
750  histos->fillCalibHist(
751  pNoiseMatrix->matrix[i].elem45, "hCalibNoise45", "Noise Matrix 45", 400, 0, 400, bin, "Calib");
752  histos->fillCalibHist(
753  pNoiseMatrix->matrix[i].elem46, "hCalibNoise46", "Noise Matrix 46", 400, 0, 400, bin, "Calib");
754  histos->fillCalibHist(
755  pNoiseMatrix->matrix[i].elem55, "hCalibNoise55", "Noise Matrix 55", 400, 0, 400, bin, "Calib");
756  histos->fillCalibHist(
757  pNoiseMatrix->matrix[i].elem56, "hCalibNoise56", "Noise Matrix 56", 400, 0, 400, bin, "Calib");
758  histos->fillCalibHist(
759  pNoiseMatrix->matrix[i].elem57, "hCalibNoise57", "Noise Matrix 57", 400, 0, 400, bin, "Calib");
760  histos->fillCalibHist(
761  pNoiseMatrix->matrix[i].elem66, "hCalibNoise66", "Noise Matrix 66", 400, 0, 400, bin, "Calib");
762  histos->fillCalibHist(
763  pNoiseMatrix->matrix[i].elem67, "hCalibNoise67", "Noise Matrix 67", 400, 0, 400, bin, "Calib");
764  histos->fillCalibHist(
765  pNoiseMatrix->matrix[i].elem77, "hCalibNoise77", "Noise Matrix 77", 400, 0, 400, bin, "Calib");
766  }
767  }
768 }
769 
770 // ==============================================
771 //
772 // look at WIRE DIGIs
773 //
774 // ==============================================
775 
776 void CSCValidation::doWireDigis(edm::Handle<CSCWireDigiCollection> wires) {
777  int nWireGroupsTotal = 0;
778  for (CSCWireDigiCollection::DigiRangeIterator dWDiter = wires->begin(); dWDiter != wires->end(); dWDiter++) {
779  CSCDetId id = (CSCDetId)(*dWDiter).first;
780  std::vector<CSCWireDigi>::const_iterator wireIter = (*dWDiter).second.first;
781  std::vector<CSCWireDigi>::const_iterator lWire = (*dWDiter).second.second;
782  for (; wireIter != lWire; ++wireIter) {
783  int myWire = wireIter->getWireGroup();
784  int myTBin = wireIter->getTimeBin();
785  nWireGroupsTotal++;
786  histos->fill1DHistByType(myWire, "hWireWire", "Wiregroup Numbers Fired", id, 113, -0.5, 112.5, "Digis");
787  histos->fill1DHistByType(myTBin, "hWireTBin", "Wire TimeBin Fired", id, 17, -0.5, 16.5, "Digis");
788  histos->fillProfile(
789  chamberSerial(id), myTBin, "hWireTBinProfile", "Wire TimeBin Fired", 601, -0.5, 600.5, -0.5, 16.5, "Digis");
790  if (detailedAnalysis) {
791  histos->fill1DHistByLayer(
792  myWire, "hWireWire", "Wiregroup Numbers Fired", id, 113, -0.5, 112.5, "WireNumberByLayer");
793  histos->fill1DHistByLayer(myTBin, "hWireTBin", "Wire TimeBin Fired", id, 17, -0.5, 16.5, "WireTimeByLayer");
794  }
795  }
796  } // end wire loop
797 
798  // this way you can zero suppress but still store info on # events with no digis
799  if (nWireGroupsTotal == 0)
800  nWireGroupsTotal = -1;
801 
802  histos->fill1DHist(nWireGroupsTotal, "hWirenGroupsTotal", "Wires Fired Per Event", 151, -0.5, 150.5, "Digis");
803 }
804 
805 // ==============================================
806 //
807 // look at STRIP DIGIs
808 //
809 // ==============================================
810 
811 void CSCValidation::doStripDigis(edm::Handle<CSCStripDigiCollection> strips) {
812  int nStripsFired = 0;
813  for (CSCStripDigiCollection::DigiRangeIterator dSDiter = strips->begin(); dSDiter != strips->end(); dSDiter++) {
814  CSCDetId id = (CSCDetId)(*dSDiter).first;
815  std::vector<CSCStripDigi>::const_iterator stripIter = (*dSDiter).second.first;
816  std::vector<CSCStripDigi>::const_iterator lStrip = (*dSDiter).second.second;
817  for (; stripIter != lStrip; ++stripIter) {
818  int myStrip = stripIter->getStrip();
819  std::vector<int> myADCVals = stripIter->getADCCounts();
820  bool thisStripFired = false;
821  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
822  float threshold = 13.3;
823  float diff = 0.;
824  for (unsigned int iCount = 0; iCount < myADCVals.size(); iCount++) {
825  diff = (float)myADCVals[iCount] - thisPedestal;
826  if (diff > threshold) {
827  thisStripFired = true;
828  }
829  }
830  if (thisStripFired) {
831  nStripsFired++;
832  // fill strip histos
833  histos->fill1DHistByType(myStrip, "hStripStrip", "Strip Number", id, 81, -0.5, 80.5, "Digis");
834  if (detailedAnalysis) {
835  histos->fill1DHistByLayer(myStrip, "hStripStrip", "Strip Number", id, 81, -0.5, 80.5, "StripNumberByLayer");
836  }
837  }
838  }
839  } // end strip loop
840 
841  if (nStripsFired == 0)
842  nStripsFired = -1;
843 
844  histos->fill1DHist(nStripsFired, "hStripNFired", "Fired Strips per Event", 251, -0.5, 250.5, "Digis");
845 }
846 
847 //=======================================================
848 //
849 // Look at the Pedestal Noise Distributions
850 //
851 //=======================================================
852 
853 void CSCValidation::doPedestalNoise(edm::Handle<CSCStripDigiCollection> strips) {
854  constexpr float threshold = 13.3;
855  for (CSCStripDigiCollection::DigiRangeIterator dPNiter = strips->begin(); dPNiter != strips->end(); dPNiter++) {
856  CSCDetId id = (CSCDetId)(*dPNiter).first;
857  std::vector<CSCStripDigi>::const_iterator pedIt = (*dPNiter).second.first;
858  std::vector<CSCStripDigi>::const_iterator lStrip = (*dPNiter).second.second;
859  for (; pedIt != lStrip; ++pedIt) {
860  int myStrip = pedIt->getStrip();
861  std::vector<int> myADCVals = pedIt->getADCCounts();
862  float TotalADC = getSignal(*strips, id, myStrip);
863  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
864  float thisSignal =
865  (1. / 6) * (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
866  bool thisStripFired = TotalADC > threshold;
867  if (!thisStripFired) {
868  float ADC = thisSignal - thisPedestal;
869  histos->fill1DHist(ADC, "hStripPed", "Pedestal Noise Distribution", 50, -25., 25., "PedestalNoise");
870  histos->fill1DHistByType(ADC, "hStripPedME", "Pedestal Noise Distribution", id, 50, -25., 25., "PedestalNoise");
871  histos->fillProfile(chamberSerial(id),
872  ADC,
873  "hStripPedMEProfile",
874  "Wire TimeBin Fired",
875  601,
876  -0.5,
877  600.5,
878  -25,
879  25,
880  "PedestalNoise");
881  if (detailedAnalysis) {
882  histos->fill1DHistByLayer(
883  ADC, "hStripPedME", "Pedestal Noise Distribution", id, 50, -25., 25., "PedestalNoiseByLayer");
884  }
885  }
886  }
887  }
888 }
889 
890 // ==============================================
891 //
892 // look at RECHITs
893 //
894 // ==============================================
895 
896 void CSCValidation::doRecHits(edm::Handle<CSCRecHit2DCollection> recHits, edm::ESHandle<CSCGeometry> cscGeom) {
897  // Get the RecHits collection :
898  int nRecHits = recHits->size();
899 
900  // ---------------------
901  // Loop over rechits
902  // ---------------------
903  int iHit = 0;
904 
905  // Build iterator for rechits and loop :
906  CSCRecHit2DCollection::const_iterator dRHIter;
907  for (dRHIter = recHits->begin(); dRHIter != recHits->end(); dRHIter++) {
908  iHit++;
909 
910  // Find chamber with rechits in CSC
911  CSCDetId idrec = (CSCDetId)(*dRHIter).cscDetId();
912  int kEndcap = idrec.endcap();
913  int kRing = idrec.ring();
914  int kStation = idrec.station();
915  int kChamber = idrec.chamber();
916  int kLayer = idrec.layer();
917 
918  // Store rechit as a Local Point:
919  LocalPoint rhitlocal = (*dRHIter).localPosition();
920  float xreco = rhitlocal.x();
921  float yreco = rhitlocal.y();
922  LocalError rerrlocal = (*dRHIter).localPositionError();
923  //these errors are squared!
924  float xxerr = rerrlocal.xx();
925  float yyerr = rerrlocal.yy();
926  float xyerr = rerrlocal.xy();
927  // errors in strip units
928  float stpos = (*dRHIter).positionWithinStrip();
929  float sterr = (*dRHIter).errorWithinStrip();
930 
931  // Find the charge associated with this hit
932  float rHSumQ = 0;
933  float sumsides = 0.;
934  int adcsize = dRHIter->nStrips() * dRHIter->nTimeBins();
935  for (unsigned int i = 0; i < dRHIter->nStrips(); i++) {
936  for (unsigned int j = 0; j < dRHIter->nTimeBins() - 1; j++) {
937  rHSumQ += dRHIter->adcs(i, j);
938  if (i != 1)
939  sumsides += dRHIter->adcs(i, j);
940  }
941  }
942 
943  float rHratioQ = sumsides / rHSumQ;
944  if (adcsize != 12)
945  rHratioQ = -99;
946 
947  // Get the signal timing of this hit
948  float rHtime = 0;
949  rHtime = (*dRHIter).tpeak() / 50.;
950 
951  // Get pointer to the layer:
952  const CSCLayer* csclayer = cscGeom->layer(idrec);
953 
954  // Transform hit position from local chamber geometry to global CMS geom
955  GlobalPoint rhitglobal = csclayer->toGlobal(rhitlocal);
956  float grecx = rhitglobal.x();
957  float grecy = rhitglobal.y();
958 
959  // Fill the rechit position branch
960  if (writeTreeToFile && rhTreeCount < 1500000) {
961  histos->fillRechitTree(xreco, yreco, grecx, grecy, kEndcap, kStation, kRing, kChamber, kLayer);
962  rhTreeCount++;
963  }
964 
965  // Fill some histograms
966  // only fill if 3 strips were used in the hit
967  histos->fill2DHistByStation(
968  grecx, grecy, "hRHGlobal", "recHit Global Position", idrec, 100, -800., 800., 100, -800., 800., "recHits");
969  if (kStation == 1 && (kRing == 1 || kRing == 4))
970  histos->fill1DHistByType(rHSumQ, "hRHSumQ", "Sum 3x3 recHit Charge", idrec, 125, 0, 4000, "recHits");
971  else
972  histos->fill1DHistByType(rHSumQ, "hRHSumQ", "Sum 3x3 recHit Charge", idrec, 125, 0, 2000, "recHits");
973  histos->fill1DHistByType(rHratioQ, "hRHRatioQ", "Charge Ratio (Ql+Qr)/Qt", idrec, 120, -0.1, 1.1, "recHits");
974  histos->fill1DHistByType(rHtime, "hRHTiming", "recHit Timing", idrec, 200, -10, 10, "recHits");
975  histos->fill1DHistByType(sqrt(xxerr), "hRHxerr", "RecHit Error on Local X", idrec, 100, -0.1, 2, "recHits");
976  histos->fill1DHistByType(sqrt(yyerr), "hRHyerr", "RecHit Error on Local Y", idrec, 100, -0.1, 2, "recHits");
977  histos->fill1DHistByType(xyerr, "hRHxyerr", "Corr. RecHit XY Error", idrec, 100, -1, 2, "recHits");
978  if (adcsize == 12)
979  histos->fill1DHistByType(stpos, "hRHstpos", "Reconstructed Position on Strip", idrec, 120, -0.6, 0.6, "recHits");
980  histos->fill1DHistByType(
981  sterr, "hRHsterr", "Estimated Error on Strip Measurement", idrec, 120, -0.05, 0.25, "recHits");
982  histos->fillProfile(
983  chamberSerial(idrec), rHSumQ, "hRHSumQProfile", "Sum 3x3 recHit Charge", 601, -0.5, 600.5, 0, 4000, "recHits");
984  histos->fillProfile(
985  chamberSerial(idrec), rHtime, "hRHTimingProfile", "recHit Timing", 601, -0.5, 600.5, -11, 11, "recHits");
986  if (detailedAnalysis) {
987  if (kStation == 1 && (kRing == 1 || kRing == 4))
988  histos->fill1DHistByLayer(rHSumQ, "hRHSumQ", "Sum 3x3 recHit Charge", idrec, 125, 0, 4000, "RHQByLayer");
989  else
990  histos->fill1DHistByLayer(rHSumQ, "hRHSumQ", "Sum 3x3 recHit Charge", idrec, 125, 0, 2000, "RHQByLayer");
991  histos->fill1DHistByLayer(rHratioQ, "hRHRatioQ", "Charge Ratio (Ql+Qr)/Qt", idrec, 120, -0.1, 1.1, "RHQByLayer");
992  histos->fill1DHistByLayer(rHtime, "hRHTiming", "recHit Timing", idrec, 200, -10, 10, "RHTimingByLayer");
993  histos->fill2DHistByLayer(xreco,
994  yreco,
995  "hRHLocalXY",
996  "recHit Local Position",
997  idrec,
998  50,
999  -100.,
1000  100.,
1001  75,
1002  -150.,
1003  150.,
1004  "RHLocalXYByLayer");
1005  histos->fill1DHistByLayer(
1006  sqrt(xxerr), "hRHxerr", "RecHit Error on Local X", idrec, 100, -0.1, 2, "RHErrorsByLayer");
1007  histos->fill1DHistByLayer(
1008  sqrt(yyerr), "hRHyerr", "RecHit Error on Local Y", idrec, 100, -0.1, 2, "RHErrorsByLayer");
1009  histos->fill1DHistByType(
1010  stpos, "hRHstpos", "Reconstructed Position on Strip", idrec, 120, -0.6, 0.6, "RHStripPosByLayer");
1011  histos->fill1DHistByType(
1012  sterr, "hRHsterr", "Estimated Error on Strip Measurement", idrec, 120, -0.05, 0.25, "RHStripPosByLayer");
1013  }
1014 
1015  } //end rechit loop
1016 
1017  if (nRecHits == 0)
1018  nRecHits = -1;
1019 
1020  histos->fill1DHist(nRecHits, "hRHnrechits", "recHits per Event (all chambers)", 151, -0.5, 150.5, "recHits");
1021 }
1022 
1023 // ==============================================
1024 //
1025 // look at SIMHITS
1026 //
1027 // ==============================================
1028 
1029 void CSCValidation::doSimHits(edm::Handle<CSCRecHit2DCollection> recHits, edm::Handle<PSimHitContainer> simHits) {
1030  CSCRecHit2DCollection::const_iterator dSHrecIter;
1031  for (dSHrecIter = recHits->begin(); dSHrecIter != recHits->end(); dSHrecIter++) {
1032  CSCDetId idrec = (CSCDetId)(*dSHrecIter).cscDetId();
1033  LocalPoint rhitlocal = (*dSHrecIter).localPosition();
1034  float xreco = rhitlocal.x();
1035  float yreco = rhitlocal.y();
1036  float xError = sqrt((*dSHrecIter).localPositionError().xx());
1037  float yError = sqrt((*dSHrecIter).localPositionError().yy());
1038  float simHitXres = -99;
1039  float simHitYres = -99;
1040  float xPull = -99;
1041  float yPull = -99;
1042  float mindiffX = 99;
1043  float mindiffY = 10;
1044  // If MC, find closest muon simHit to check resolution:
1045  PSimHitContainer::const_iterator dSHsimIter;
1046  for (dSHsimIter = simHits->begin(); dSHsimIter != simHits->end(); dSHsimIter++) {
1047  // Get DetID for this simHit:
1048  CSCDetId sId = (CSCDetId)(*dSHsimIter).detUnitId();
1049  // Check if the simHit detID matches that of current recHit
1050  // and make sure it is a muon hit:
1051  if (sId == idrec && abs((*dSHsimIter).particleType()) == 13) {
1052  // Get the position of this simHit in local coordinate system:
1053  LocalPoint sHitlocal = (*dSHsimIter).localPosition();
1054  // Now we need to make reasonably sure that this simHit is
1055  // responsible for this recHit:
1056  if ((sHitlocal.x() - xreco) < mindiffX && (sHitlocal.y() - yreco) < mindiffY) {
1057  simHitXres = (sHitlocal.x() - xreco);
1058  simHitYres = (sHitlocal.y() - yreco);
1059  mindiffX = (sHitlocal.x() - xreco);
1060  xPull = simHitXres / xError;
1061  yPull = simHitYres / yError;
1062  }
1063  }
1064  }
1065 
1066  histos->fill1DHistByType(
1067  simHitXres, "hSimXResid", "SimHitX - Reconstructed X", idrec, 100, -1.0, 1.0, "Resolution");
1068  histos->fill1DHistByType(
1069  simHitYres, "hSimYResid", "SimHitY - Reconstructed Y", idrec, 100, -5.0, 5.0, "Resolution");
1070  histos->fill1DHistByType(xPull, "hSimXPull", "Local X Pulls", idrec, 100, -5.0, 5.0, "Resolution");
1071  histos->fill1DHistByType(yPull, "hSimYPull", "Local Y Pulls", idrec, 100, -5.0, 5.0, "Resolution");
1072  }
1073 }
1074 
1075 // ==============================================
1076 //
1077 // look at SEGMENTs
1078 //
1079 // ===============================================
1080 
1081 void CSCValidation::doSegments(edm::Handle<CSCSegmentCollection> cscSegments, edm::ESHandle<CSCGeometry> cscGeom) {
1082  // get CSC segment collection
1083  int nSegments = cscSegments->size();
1084 
1085  // -----------------------
1086  // loop over segments
1087  // -----------------------
1088  int iSegment = 0;
1089  for (CSCSegmentCollection::const_iterator dSiter = cscSegments->begin(); dSiter != cscSegments->end(); dSiter++) {
1090  iSegment++;
1091  //
1092  CSCDetId id = (CSCDetId)(*dSiter).cscDetId();
1093  int kEndcap = id.endcap();
1094  int kRing = id.ring();
1095  int kStation = id.station();
1096  int kChamber = id.chamber();
1097 
1098  //
1099  float chisq = (*dSiter).chi2();
1100  int nhits = (*dSiter).nRecHits();
1101  int nDOF = 2 * nhits - 4;
1102  double chisqProb = ChiSquaredProbability((double)chisq, nDOF);
1103  LocalPoint localPos = (*dSiter).localPosition();
1104  float segX = localPos.x();
1105  float segY = localPos.y();
1106  LocalVector segDir = (*dSiter).localDirection();
1107  double theta = segDir.theta();
1108 
1109  // global transformation
1110  float globX = 0.;
1111  float globY = 0.;
1112  float globTheta = 0.;
1113  float globPhi = 0.;
1114  const CSCChamber* cscchamber = cscGeom->chamber(id);
1115  if (cscchamber) {
1116  GlobalPoint globalPosition = cscchamber->toGlobal(localPos);
1117  globX = globalPosition.x();
1118  globY = globalPosition.y();
1119  GlobalVector globalDirection = cscchamber->toGlobal(segDir);
1120  globTheta = globalDirection.theta();
1121  globPhi = globalDirection.phi();
1122  }
1123 
1124  // Fill segment position branch
1125  if (writeTreeToFile && segTreeCount < 1500000) {
1126  histos->fillSegmentTree(segX, segY, globX, globY, kEndcap, kStation, kRing, kChamber);
1127  segTreeCount++;
1128  }
1129 
1130  // Fill histos
1131  histos->fill2DHistByStation(globX,
1132  globY,
1133  "hSGlobal",
1134  "Segment Global Positions;global x (cm)",
1135  id,
1136  100,
1137  -800.,
1138  800.,
1139  100,
1140  -800.,
1141  800.,
1142  "Segments");
1143  histos->fill1DHistByType(nhits, "hSnHits", "N hits on Segments", id, 8, -0.5, 7.5, "Segments");
1144  histos->fill1DHistByType(theta, "hSTheta", "local theta segments", id, 128, -3.2, 3.2, "Segments");
1145  histos->fill1DHistByType((chisq / nDOF), "hSChiSq", "segments chi-squared/ndof", id, 110, -0.05, 10.5, "Segments");
1146  histos->fill1DHistByType(
1147  chisqProb, "hSChiSqProb", "segments chi-squared probability", id, 110, -0.05, 1.05, "Segments");
1148  histos->fill1DHist(globTheta, "hSGlobalTheta", "segment global theta", 128, 0, 3.2, "Segments");
1149  histos->fill1DHist(globPhi, "hSGlobalPhi", "segment global phi", 128, -3.2, 3.2, "Segments");
1150  histos->fillProfile(
1151  chamberSerial(id), nhits, "hSnHitsProfile", "N hits on Segments", 601, -0.5, 600.5, -0.5, 7.5, "Segments");
1152  if (detailedAnalysis) {
1153  histos->fill1DHistByChamber(nhits, "hSnHits", "N hits on Segments", id, 8, -0.5, 7.5, "HitsOnSegmentByChamber");
1154  histos->fill1DHistByChamber(theta, "hSTheta", "local theta segments", id, 128, -3.2, 3.2, "DetailedSegments");
1155  histos->fill1DHistByChamber(
1156  (chisq / nDOF), "hSChiSq", "segments chi-squared/ndof", id, 110, -0.05, 10.5, "SegChi2ByChamber");
1157  histos->fill1DHistByChamber(
1158  chisqProb, "hSChiSqProb", "segments chi-squared probability", id, 110, -0.05, 1.05, "SegChi2ByChamber");
1159  }
1160 
1161  } // end segment loop
1162 
1163  if (nSegments == 0)
1164  nSegments = -1;
1165 
1166  histos->fill1DHist(nSegments, "hSnSegments", "Segments per Event", 31, -0.5, 30.5, "Segments");
1167 }
1168 
1169 // ==============================================
1170 //
1171 // look at hit Resolution
1172 //
1173 // ===============================================
1174 
1175 void CSCValidation::doResolution(edm::Handle<CSCSegmentCollection> cscSegments, edm::ESHandle<CSCGeometry> cscGeom) {
1176  for (CSCSegmentCollection::const_iterator dSiter = cscSegments->begin(); dSiter != cscSegments->end(); dSiter++) {
1177  CSCDetId id = (CSCDetId)(*dSiter).cscDetId();
1178 
1179  //
1180  // try to get the CSC recHits that contribute to this segment.
1181  std::vector<CSCRecHit2D> theseRecHits = (*dSiter).specificRecHits();
1182  int nRH = (*dSiter).nRecHits();
1183  int jRH = 0;
1184  CLHEP::HepMatrix sp(6, 1);
1185  CLHEP::HepMatrix se(6, 1);
1186  for (std::vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
1187  jRH++;
1188  CSCDetId idRH = (CSCDetId)(*iRH).cscDetId();
1189  int kRing = idRH.ring();
1190  int kStation = idRH.station();
1191  int kLayer = idRH.layer();
1192 
1193  // Find the strip containing this hit
1194  int centerid = iRH->nStrips() / 2;
1195  int centerStrip = iRH->channels(centerid);
1196 
1197  // If this segment has 6 hits, find the position of each hit on the strip in units of stripwidth and store values
1198  if (nRH == 6) {
1199  float stpos = (*iRH).positionWithinStrip();
1200  se(kLayer, 1) = (*iRH).errorWithinStrip();
1201  // Take into account half-strip staggering of layers (ME1/1 has no staggering)
1202  if (kStation == 1 && (kRing == 1 || kRing == 4))
1203  sp(kLayer, 1) = stpos + centerStrip;
1204  else {
1205  if (kLayer == 1 || kLayer == 3 || kLayer == 5)
1206  sp(kLayer, 1) = stpos + centerStrip;
1207  if (kLayer == 2 || kLayer == 4 || kLayer == 6)
1208  sp(kLayer, 1) = stpos - 0.5 + centerStrip;
1209  }
1210  }
1211  }
1212 
1213  float residual = -99;
1214  float pull = -99;
1215  // Fit all points except layer 3, then compare expected value for layer 3 to reconstructed value
1216  if (nRH == 6) {
1217  float expected = fitX(sp, se);
1218  residual = expected - sp(3, 1);
1219  pull = residual / se(3, 1);
1220  }
1221 
1222  // Fill histos
1223  histos->fill1DHistByType(
1224  residual, "hSResid", "Fitted Position on Strip - Reconstructed for Layer 3", id, 100, -0.5, 0.5, "Resolution");
1225  histos->fill1DHistByType(pull, "hSStripPosPull", "Strip Measurement Pulls", id, 100, -5.0, 5.0, "Resolution");
1226  histos->fillProfile(chamberSerial(id),
1227  residual,
1228  "hSResidProfile",
1229  "Fitted Position on Strip - Reconstructed for Layer 3",
1230  601,
1231  -0.5,
1232  600.5,
1233  -0.5,
1234  0.5,
1235  "Resolution");
1236  if (detailedAnalysis) {
1237  histos->fill1DHistByChamber(residual,
1238  "hSResid",
1239  "Fitted Position on Strip - Reconstructed for Layer 3",
1240  id,
1241  100,
1242  -0.5,
1243  0.5,
1244  "DetailedResolution");
1245  histos->fill1DHistByChamber(pull, "hSStripPosPull", "Strip Measurement Pulls", id, 100, -5.0, 5.0, "Resolution");
1246  }
1247  }
1248 }
1249 
1250 // ==============================================
1251 //
1252 // look at Standalone Muons
1253 //
1254 // ==============================================
1255 
1256 void CSCValidation::doStandalone(Handle<reco::TrackCollection> saMuons) {
1257  int nSAMuons = saMuons->size();
1258  histos->fill1DHist(nSAMuons, "trNSAMuons", "N Standalone Muons per Event", 6, -0.5, 5.5, "STAMuons");
1259 
1260  for (reco::TrackCollection::const_iterator muon = saMuons->begin(); muon != saMuons->end(); ++muon) {
1261  float preco = muon->p();
1262  float ptreco = muon->pt();
1263  int n = muon->recHitsSize();
1264  float chi2 = muon->chi2();
1265  float normchi2 = muon->normalizedChi2();
1266 
1267  // loop over hits
1268  int nDTHits = 0;
1269  int nCSCHits = 0;
1270  int nCSCHitsp = 0;
1271  int nCSCHitsm = 0;
1272  int nRPCHits = 0;
1273  int nRPCHitsp = 0;
1274  int nRPCHitsm = 0;
1275  int np = 0;
1276  int nm = 0;
1277  std::vector<CSCDetId> staChambers;
1278  for (trackingRecHit_iterator hit = muon->recHitsBegin(); hit != muon->recHitsEnd(); ++hit) {
1279  const DetId detId((*hit)->geographicalId());
1280  if (detId.det() == DetId::Muon) {
1281  if (detId.subdetId() == MuonSubdetId::RPC) {
1282  RPCDetId rpcId(detId.rawId());
1283  nRPCHits++;
1284  if (rpcId.region() == 1) {
1285  nRPCHitsp++;
1286  np++;
1287  }
1288  if (rpcId.region() == -1) {
1289  nRPCHitsm++;
1290  nm++;
1291  }
1292  }
1293  if (detId.subdetId() == MuonSubdetId::DT) {
1294  nDTHits++;
1295  } else if (detId.subdetId() == MuonSubdetId::CSC) {
1296  CSCDetId cscId(detId.rawId());
1297  staChambers.push_back(detId.rawId());
1298  nCSCHits++;
1299  if (cscId.endcap() == 1) {
1300  nCSCHitsp++;
1301  np++;
1302  }
1303  if (cscId.endcap() == 2) {
1304  nCSCHitsm++;
1305  nm++;
1306  }
1307  }
1308  }
1309  }
1310 
1311  GlobalPoint innerPnt(muon->innerPosition().x(), muon->innerPosition().y(), muon->innerPosition().z());
1312  GlobalPoint outerPnt(muon->outerPosition().x(), muon->outerPosition().y(), muon->outerPosition().z());
1313  GlobalVector innerKin(muon->innerMomentum().x(), muon->innerMomentum().y(), muon->innerMomentum().z());
1314  GlobalVector outerKin(muon->outerMomentum().x(), muon->outerMomentum().y(), muon->outerMomentum().z());
1315  GlobalVector deltaPnt = innerPnt - outerPnt;
1316  double crudeLength = deltaPnt.mag();
1317  double deltaPhi = innerPnt.phi() - outerPnt.phi();
1318  double innerGlobalPolarAngle = innerKin.theta();
1319  double outerGlobalPolarAngle = outerKin.theta();
1320 
1321  // fill histograms
1322  histos->fill1DHist(n, "trN", "N hits on a STA Muon Track", 51, -0.5, 50.5, "STAMuons");
1323  if (np != 0)
1324  histos->fill1DHist(np, "trNp", "N hits on a STA Muon Track (plus endcap)", 51, -0.5, 50.5, "STAMuons");
1325  if (nm != 0)
1326  histos->fill1DHist(nm, "trNm", "N hits on a STA Muon Track (minus endcap)", 51, -0.5, 50.5, "STAMuons");
1327  histos->fill1DHist(nDTHits, "trNDT", "N DT hits on a STA Muon Track", 51, -0.5, 50.5, "STAMuons");
1328  histos->fill1DHist(nCSCHits, "trNCSC", "N CSC hits on a STA Muon Track", 51, -0.5, 50.5, "STAMuons");
1329  if (nCSCHitsp != 0)
1330  histos->fill1DHist(nCSCHitsp, "trNCSCp", "N CSC hits on a STA Muon Track (+ endcap)", 51, -0.5, 50.5, "STAMuons");
1331  if (nCSCHitsm != 0)
1332  histos->fill1DHist(nCSCHitsm, "trNCSCm", "N CSC hits on a STA Muon Track (- endcap)", 51, -0.5, 50.5, "STAMuons");
1333  histos->fill1DHist(nRPCHits, "trNRPC", "N RPC hits on a STA Muon Track", 51, -0.5, 50.5, "STAMuons");
1334  if (nRPCHitsp != 0)
1335  histos->fill1DHist(nRPCHitsp, "trNRPCp", "N RPC hits on a STA Muon Track (+ endcap)", 51, -0.5, 50.5, "STAMuons");
1336  if (nRPCHitsm != 0)
1337  histos->fill1DHist(nRPCHitsm, "trNRPCm", "N RPC hits on a STA Muon Track (- endcap)", 51, -0.5, 50.5, "STAMuons");
1338  histos->fill1DHist(preco, "trP", "STA Muon Momentum", 100, 0, 300, "STAMuons");
1339  histos->fill1DHist(ptreco, "trPT", "STA Muon pT", 100, 0, 40, "STAMuons");
1340  histos->fill1DHist(chi2, "trChi2", "STA Muon Chi2", 100, 0, 200, "STAMuons");
1341  histos->fill1DHist(normchi2, "trNormChi2", "STA Muon Normalized Chi2", 100, 0, 10, "STAMuons");
1342  histos->fill1DHist(crudeLength, "trLength", "Straight Line Length of STA Muon", 120, 0., 2400., "STAMuons");
1343  histos->fill1DHist(
1344  deltaPhi, "trDeltaPhi", "Delta-Phi Between Inner and Outer STA Muon Pos.", 100, -0.5, 0.5, "STAMuons");
1345  histos->fill1DHist(
1346  innerGlobalPolarAngle, "trInnerPolar", "Polar Angle of Inner P Vector (STA muons)", 128, 0, 3.2, "STAMuons");
1347  histos->fill1DHist(
1348  outerGlobalPolarAngle, "trOuterPolar", "Polar Angle of Outer P Vector (STA muons)", 128, 0, 3.2, "STAMuons");
1349  histos->fill1DHist(innerPnt.phi(), "trInnerPhi", "Phi of Inner Position (STA muons)", 256, -3.2, 3.2, "STAMuons");
1350  histos->fill1DHist(outerPnt.phi(), "trOuterPhi", "Phi of Outer Position (STA muons)", 256, -3.2, 3.2, "STAMuons");
1351  }
1352 }
1353 
1354 //--------------------------------------------------------------
1355 // Compute a serial number for the chamber.
1356 // This is useful when filling histograms and working with arrays.
1357 //--------------------------------------------------------------
1358 int CSCValidation::chamberSerial(CSCDetId id) {
1359  int st = id.station();
1360  int ri = id.ring();
1361  int ch = id.chamber();
1362  int ec = id.endcap();
1363  int kSerial = ch;
1364  if (st == 1 && ri == 1)
1365  kSerial = ch;
1366  if (st == 1 && ri == 2)
1367  kSerial = ch + 36;
1368  if (st == 1 && ri == 3)
1369  kSerial = ch + 72;
1370  if (st == 1 && ri == 4)
1371  kSerial = ch;
1372  if (st == 2 && ri == 1)
1373  kSerial = ch + 108;
1374  if (st == 2 && ri == 2)
1375  kSerial = ch + 126;
1376  if (st == 3 && ri == 1)
1377  kSerial = ch + 162;
1378  if (st == 3 && ri == 2)
1379  kSerial = ch + 180;
1380  if (st == 4 && ri == 1)
1381  kSerial = ch + 216;
1382  if (st == 4 && ri == 2)
1383  kSerial = ch + 234; // one day...
1384  if (ec == 2)
1385  kSerial = kSerial + 300;
1386  return kSerial;
1387 }
1388 
1389 //--------------------------------------------------------------
1390 // Compute a serial number for the ring.
1391 // This is useful when filling histograms and working with arrays.
1392 //--------------------------------------------------------------
1393 int CSCValidation::ringSerial(CSCDetId id) {
1394  int st = id.station();
1395  int ri = id.ring();
1396  int ec = id.endcap();
1397  int kSerial = 0;
1398  if (st == 1 && ri == 1)
1399  kSerial = ri;
1400  if (st == 1 && ri == 2)
1401  kSerial = ri;
1402  if (st == 1 && ri == 3)
1403  kSerial = ri;
1404  if (st == 1 && ri == 4)
1405  kSerial = 1;
1406  if (st == 2)
1407  kSerial = ri + 3;
1408  if (st == 3)
1409  kSerial = ri + 5;
1410  if (st == 4)
1411  kSerial = ri + 7;
1412  if (ec == 2)
1413  kSerial = kSerial * (-1);
1414  return kSerial;
1415 }
1416 
1417 //-------------------------------------------------------------------------------------
1418 // Fits a straight line to a set of 5 points with errors. Functions assumes 6 points
1419 // and removes hit in layer 3. It then returns the expected position value in layer 3
1420 // based on the fit.
1421 //-------------------------------------------------------------------------------------
1422 float CSCValidation::fitX(const CLHEP::HepMatrix& points, const CLHEP::HepMatrix& errors) {
1423  float S = 0;
1424  float Sx = 0;
1425  float Sy = 0;
1426  float Sxx = 0;
1427  float Sxy = 0;
1428  float sigma2 = 0;
1429 
1430  for (int i = 1; i < 7; i++) {
1431  if (i != 3) {
1432  sigma2 = errors(i, 1) * errors(i, 1);
1433  S = S + (1 / sigma2);
1434  Sy = Sy + (points(i, 1) / sigma2);
1435  Sx = Sx + ((i) / sigma2);
1436  Sxx = Sxx + (i * i) / sigma2;
1437  Sxy = Sxy + (((i)*points(i, 1)) / sigma2);
1438  }
1439  }
1440 
1441  float delta = S * Sxx - Sx * Sx;
1442  float intercept = (Sxx * Sy - Sx * Sxy) / delta;
1443  float slope = (S * Sxy - Sx * Sy) / delta;
1444 
1445  //float chi = 0;
1446  //float chi2 = 0;
1447 
1448  // calculate chi2 (not currently used)
1449  //for (int i=1;i<7;i++){
1450  // chi = (points(i,1) - intercept - slope*i)/(errors(i,1));
1451  // chi2 = chi2 + chi*chi;
1452  //}
1453 
1454  return (intercept + slope * 3);
1455 }
1456 
1457 //----------------------------------------------------------------------------
1458 // Calculate basic efficiencies for recHits and Segments
1459 // Author: S. Stoynev
1460 //----------------------------------------------------------------------------
1461 
1462 void CSCValidation::doEfficiencies(edm::Handle<CSCWireDigiCollection> wires,
1463  edm::Handle<CSCStripDigiCollection> strips,
1464  edm::Handle<CSCRecHit2DCollection> recHits,
1465  edm::Handle<CSCSegmentCollection> cscSegments,
1466  edm::ESHandle<CSCGeometry> cscGeom) {
1467  bool allWires[2][4][4][36][6];
1468  bool allStrips[2][4][4][36][6];
1469  bool AllRecHits[2][4][4][36][6];
1470  bool AllSegments[2][4][4][36];
1471 
1472  //bool MultiSegments[2][4][4][36];
1473  for (int iE = 0; iE < 2; iE++) {
1474  for (int iS = 0; iS < 4; iS++) {
1475  for (int iR = 0; iR < 4; iR++) {
1476  for (int iC = 0; iC < 36; iC++) {
1477  AllSegments[iE][iS][iR][iC] = false;
1478  //MultiSegments[iE][iS][iR][iC] = false;
1479  for (int iL = 0; iL < 6; iL++) {
1480  allWires[iE][iS][iR][iC][iL] = false;
1481  allStrips[iE][iS][iR][iC][iL] = false;
1482  AllRecHits[iE][iS][iR][iC][iL] = false;
1483  }
1484  }
1485  }
1486  }
1487  }
1488 
1489  if (useDigis) {
1490  // Wires
1491  for (CSCWireDigiCollection::DigiRangeIterator dWDiter = wires->begin(); dWDiter != wires->end(); dWDiter++) {
1492  CSCDetId idrec = (CSCDetId)(*dWDiter).first;
1493  std::vector<CSCWireDigi>::const_iterator wireIter = (*dWDiter).second.first;
1494  std::vector<CSCWireDigi>::const_iterator lWire = (*dWDiter).second.second;
1495  for (; wireIter != lWire; ++wireIter) {
1496  allWires[idrec.endcap() - 1][idrec.station() - 1][idrec.ring() - 1][idrec.chamber() - 1][idrec.layer() - 1] =
1497  true;
1498  break;
1499  }
1500  }
1501 
1502  //---- STRIPS
1503  for (CSCStripDigiCollection::DigiRangeIterator dSDiter = strips->begin(); dSDiter != strips->end(); dSDiter++) {
1504  CSCDetId idrec = (CSCDetId)(*dSDiter).first;
1505  std::vector<CSCStripDigi>::const_iterator stripIter = (*dSDiter).second.first;
1506  std::vector<CSCStripDigi>::const_iterator lStrip = (*dSDiter).second.second;
1507  for (; stripIter != lStrip; ++stripIter) {
1508  std::vector<int> myADCVals = stripIter->getADCCounts();
1509  bool thisStripFired = false;
1510  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
1511  float threshold = 13.3;
1512  float diff = 0.;
1513  for (unsigned int iCount = 0; iCount < myADCVals.size(); iCount++) {
1514  diff = (float)myADCVals[iCount] - thisPedestal;
1515  if (diff > threshold) {
1516  thisStripFired = true;
1517  break;
1518  }
1519  }
1520  if (thisStripFired) {
1521  allStrips[idrec.endcap() - 1][idrec.station() - 1][idrec.ring() - 1][idrec.chamber() - 1][idrec.layer() - 1] =
1522  true;
1523  break;
1524  }
1525  }
1526  }
1527  }
1528 
1529  // Rechits
1530  for (CSCRecHit2DCollection::const_iterator recEffIt = recHits->begin(); recEffIt != recHits->end(); recEffIt++) {
1531  //CSCDetId idrec = (CSCDetId)(*recIt).cscDetId();
1532  CSCDetId idrec = (CSCDetId)(*recEffIt).cscDetId();
1533  AllRecHits[idrec.endcap() - 1][idrec.station() - 1][idrec.ring() - 1][idrec.chamber() - 1][idrec.layer() - 1] =
1534  true;
1535  }
1536 
1537  std::vector<unsigned int> seg_ME2(2, 0);
1538  std::vector<unsigned int> seg_ME3(2, 0);
1539  std::vector<std::pair<CSCDetId, CSCSegment> > theSegments(4);
1540  // Segments
1541  for (CSCSegmentCollection::const_iterator segEffIt = cscSegments->begin(); segEffIt != cscSegments->end();
1542  segEffIt++) {
1543  CSCDetId idseg = (CSCDetId)(*segEffIt).cscDetId();
1544  //if(AllSegments[idrec.endcap() -1][idrec.station() -1][idrec.ring() -1][idrec.chamber()]){
1545  //MultiSegments[idrec.endcap() -1][idrec.station() -1][idrec.ring() -1][idrec.chamber()] = true;
1546  //}
1547  AllSegments[idseg.endcap() - 1][idseg.station() - 1][idseg.ring() - 1][idseg.chamber() - 1] = true;
1548  // "Intrinsic" efficiency measurement relies on "good" segment extrapolation - we need the pre-selection below
1549  // station 2 "good" segment will be used for testing efficiencies in ME1 and ME3
1550  // station 3 "good" segment will be used for testing efficiencies in ME2 and ME4
1551  if (2 == idseg.station() || 3 == idseg.station()) {
1552  unsigned int seg_tmp;
1553  if (2 == idseg.station()) {
1554  ++seg_ME2[idseg.endcap() - 1];
1555  seg_tmp = seg_ME2[idseg.endcap() - 1];
1556  } else {
1557  ++seg_ME3[idseg.endcap() - 1];
1558  seg_tmp = seg_ME3[idseg.endcap() - 1];
1559  }
1560  // is the segment good
1561  if (1 == seg_tmp && 6 == (*segEffIt).nRecHits() && (*segEffIt).chi2() / (*segEffIt).degreesOfFreedom() < 3.) {
1562  std::pair<CSCDetId, CSCSegment> specSeg = make_pair((CSCDetId)(*segEffIt).cscDetId(), *segEffIt);
1563  theSegments[2 * (idseg.endcap() - 1) + (idseg.station() - 2)] = specSeg;
1564  }
1565  }
1566  /*
1567  if(2==idseg.station()){
1568  ++seg_ME2[idseg.endcap() -1];
1569  if(1==seg_ME2[idseg.endcap() -1] && 6==(*segEffIt).nRecHits() && (*segEffIt).chi2()/(*segEffIt).degreesOfFreedom()<3.){
1570  std::pair <CSCDetId, CSCSegment> specSeg = make_pair( (CSCDetId)(*segEffIt).cscDetId(),*segEffIt);
1571  theSegments[2*(idseg.endcap()-1)+(idseg.station() -2)] = specSeg;
1572  }
1573  }
1574  else if(3==idseg.station()){
1575  ++seg_ME3[idseg.endcap() -1];
1576  if(1==seg_ME3[idseg.endcap() -1] && 6==(*segEffIt).nRecHits() && (*segEffIt).chi2()/(*segEffIt).degreesOfFreedom()<3.){
1577  std::pair <CSCDetId, CSCSegment> specSeg = make_pair( (CSCDetId)(*segEffIt).cscDetId(),*segEffIt);
1578  theSegments[2*(idseg.endcap()-1)+(idseg.station() -2)] = specSeg;
1579  }
1580  }
1581  */
1582  }
1583  // Simple efficiency calculations
1584  for (int iE = 0; iE < 2; iE++) {
1585  for (int iS = 0; iS < 4; iS++) {
1586  for (int iR = 0; iR < 4; iR++) {
1587  for (int iC = 0; iC < 36; iC++) {
1588  int NumberOfLayers = 0;
1589  for (int iL = 0; iL < 6; iL++) {
1590  if (AllRecHits[iE][iS][iR][iC][iL]) {
1591  NumberOfLayers++;
1592  }
1593  }
1594  int bin = 0;
1595  if (iS == 0)
1596  bin = iR + 1 + (iE * 10);
1597  else
1598  bin = (iS + 1) * 2 + (iR + 1) + (iE * 10);
1599  if (NumberOfLayers > 1) {
1600  //if(!(MultiSegments[iE][iS][iR][iC])){
1601  if (AllSegments[iE][iS][iR][iC]) {
1602  //---- Efficient segment evenents
1603  hSSTE->AddBinContent(bin);
1604  }
1605  //---- All segment events (normalization)
1606  hSSTE->AddBinContent(20 + bin);
1607  //}
1608  }
1609  if (AllSegments[iE][iS][iR][iC]) {
1610  if (NumberOfLayers == 6) {
1611  //---- Efficient rechit events
1612  hRHSTE->AddBinContent(bin);
1613  ;
1614  }
1615  //---- All rechit events (normalization)
1616  hRHSTE->AddBinContent(20 + bin);
1617  ;
1618  }
1619  }
1620  }
1621  }
1622  }
1623 
1624  // pick a segment only if there are no others in the station
1625  std::vector<std::pair<CSCDetId, CSCSegment>*> theSeg;
1626  if (1 == seg_ME2[0])
1627  theSeg.push_back(&theSegments[0]);
1628  if (1 == seg_ME3[0])
1629  theSeg.push_back(&theSegments[1]);
1630  if (1 == seg_ME2[1])
1631  theSeg.push_back(&theSegments[2]);
1632  if (1 == seg_ME3[1])
1633  theSeg.push_back(&theSegments[3]);
1634 
1635  // Needed for plots
1636  // at the end the chamber types will be numbered as 1 to 18
1637  // (ME-4/1, -ME3/2, -ME3/1, ..., +ME3/1, +ME3/2, ME+4/1 )
1638  std::map<std::string, float> chamberTypes;
1639  chamberTypes["ME1/a"] = 0.5;
1640  chamberTypes["ME1/b"] = 1.5;
1641  chamberTypes["ME1/2"] = 2.5;
1642  chamberTypes["ME1/3"] = 3.5;
1643  chamberTypes["ME2/1"] = 4.5;
1644  chamberTypes["ME2/2"] = 5.5;
1645  chamberTypes["ME3/1"] = 6.5;
1646  chamberTypes["ME3/2"] = 7.5;
1647  chamberTypes["ME4/1"] = 8.5;
1648 
1649  if (!theSeg.empty()) {
1650  std::map<int, GlobalPoint> extrapolatedPoint;
1651  std::map<int, GlobalPoint>::iterator it;
1652  const CSCGeometry::ChamberContainer& ChamberContainer = cscGeom->chambers();
1653  // Pick which chamber with which segment to test
1654  for (size_t nCh = 0; nCh < ChamberContainer.size(); nCh++) {
1655  const CSCChamber* cscchamber = ChamberContainer[nCh];
1656  std::pair<CSCDetId, CSCSegment>* thisSegment = nullptr;
1657  for (size_t iSeg = 0; iSeg < theSeg.size(); ++iSeg) {
1658  if (cscchamber->id().endcap() == theSeg[iSeg]->first.endcap()) {
1659  if (1 == cscchamber->id().station() || 3 == cscchamber->id().station()) {
1660  if (2 == theSeg[iSeg]->first.station()) {
1661  thisSegment = theSeg[iSeg];
1662  }
1663  } else if (2 == cscchamber->id().station() || 4 == cscchamber->id().station()) {
1664  if (3 == theSeg[iSeg]->first.station()) {
1665  thisSegment = theSeg[iSeg];
1666  }
1667  }
1668  }
1669  }
1670  // this chamber is to be tested with thisSegment
1671  if (thisSegment) {
1672  CSCSegment* seg = &(thisSegment->second);
1673  const CSCChamber* segChamber = cscGeom->chamber(thisSegment->first);
1674  LocalPoint localCenter(0., 0., 0);
1675  GlobalPoint cscchamberCenter = cscchamber->toGlobal(localCenter);
1676  // try to save some time (extrapolate a segment to a certain position only once)
1677  it = extrapolatedPoint.find(int(cscchamberCenter.z()));
1678  if (it == extrapolatedPoint.end()) {
1679  GlobalPoint segPos = segChamber->toGlobal(seg->localPosition());
1680  GlobalVector segDir = segChamber->toGlobal(seg->localDirection());
1681  double paramaterLine = lineParametrization(segPos.z(), cscchamberCenter.z(), segDir.z());
1682  double xExtrapolated = extrapolate1D(segPos.x(), segDir.x(), paramaterLine);
1683  double yExtrapolated = extrapolate1D(segPos.y(), segDir.y(), paramaterLine);
1684  GlobalPoint globP(xExtrapolated, yExtrapolated, cscchamberCenter.z());
1685  extrapolatedPoint[int(cscchamberCenter.z())] = globP;
1686  }
1687  // Where does the extrapolated point lie in the (tested) chamber local frame? Here:
1688  LocalPoint extrapolatedPointLocal = cscchamber->toLocal(extrapolatedPoint[int(cscchamberCenter.z())]);
1689  const CSCLayer* layer_p = cscchamber->layer(1); //layer 1
1690  const CSCLayerGeometry* layerGeom = layer_p->geometry();
1691  const std::array<const float, 4>& layerBounds = layerGeom->parameters();
1692  float shiftFromEdge = 15.; //cm
1693  float shiftFromDeadZone = 10.;
1694  // is the extrapolated point within a sensitive region
1695  bool pass = withinSensitiveRegion(extrapolatedPointLocal,
1696  layerBounds,
1697  cscchamber->id().station(),
1698  cscchamber->id().ring(),
1699  shiftFromEdge,
1700  shiftFromDeadZone);
1701  if (pass) { // the extrapolation point of the segment lies within sensitive region of that chamber
1702  // how many rechit layers are there in the chamber?
1703  // 0 - maybe the muon died or is deflected at large angle? do not use that case
1704  // 1 - could be noise...
1705  // 2 or more - this is promissing; this is our definition of a reliable signal; use it below
1706  // is other definition better?
1707  int nRHLayers = 0;
1708  for (int iL = 0; iL < 6; ++iL) {
1709  if (AllRecHits[cscchamber->id().endcap() - 1][cscchamber->id().station() - 1][cscchamber->id().ring() - 1]
1710  [cscchamber->id().chamber() - 1][iL]) {
1711  ++nRHLayers;
1712  }
1713  }
1714  //std::cout<<" nRHLayers = "<<nRHLayers<<std::endl;
1715  float verticalScale = chamberTypes[cscchamber->specs()->chamberTypeName()];
1716  if (cscchamberCenter.z() < 0) {
1717  verticalScale = -verticalScale;
1718  }
1719  verticalScale += 9.5;
1720  hSensitiveAreaEvt->Fill(float(cscchamber->id().chamber()), verticalScale);
1721  if (nRHLayers > 1) { // this chamber contains a reliable signal
1722  //chamberTypes[cscchamber->specs()->chamberTypeName()];
1723  // "intrinsic" efficiencies
1724  //std::cout<<" verticalScale = "<<verticalScale<<" chType = "<<cscchamber->specs()->chamberTypeName()<<std::endl;
1725  // this is the denominator forr all efficiencies
1726  hEffDenominator->Fill(float(cscchamber->id().chamber()), verticalScale);
1727  // Segment efficiency
1728  if (AllSegments[cscchamber->id().endcap() - 1][cscchamber->id().station() - 1][cscchamber->id().ring() - 1]
1729  [cscchamber->id().chamber() - 1]) {
1730  hSSTE2->Fill(float(cscchamber->id().chamber()), float(verticalScale));
1731  }
1732 
1733  for (int iL = 0; iL < 6; ++iL) {
1734  float weight = 1. / 6.;
1735  // one shold account for the weight in the efficiency...
1736  // Rechit efficiency
1737  if (AllRecHits[cscchamber->id().endcap() - 1][cscchamber->id().station() - 1][cscchamber->id().ring() - 1]
1738  [cscchamber->id().chamber() - 1][iL]) {
1739  hRHSTE2->Fill(float(cscchamber->id().chamber()), float(verticalScale), weight);
1740  }
1741  if (useDigis) {
1742  // Wire efficiency
1743  if (allWires[cscchamber->id().endcap() - 1][cscchamber->id().station() - 1][cscchamber->id().ring() - 1]
1744  [cscchamber->id().chamber() - 1][iL]) {
1745  // one shold account for the weight in the efficiency...
1746  hWireSTE2->Fill(float(cscchamber->id().chamber()), float(verticalScale), weight);
1747  }
1748  // Strip efficiency
1749  if (allStrips[cscchamber->id().endcap() - 1][cscchamber->id().station() - 1]
1750  [cscchamber->id().ring() - 1][cscchamber->id().chamber() - 1][iL]) {
1751  // one shold account for the weight in the efficiency...
1752  hStripSTE2->Fill(float(cscchamber->id().chamber()), float(verticalScale), weight);
1753  }
1754  }
1755  }
1756  }
1757  }
1758  }
1759  }
1760  }
1761  //
1762 }
1763 
1764 void CSCValidation::getEfficiency(float bin, float Norm, std::vector<float>& eff) {
1765  //---- Efficiency with binomial error
1766  float Efficiency = 0.;
1767  float EffError = 0.;
1768  if (fabs(Norm) > 0.000000001) {
1769  Efficiency = bin / Norm;
1770  if (bin < Norm) {
1771  EffError = sqrt((1. - Efficiency) * Efficiency / Norm);
1772  }
1773  }
1774  eff[0] = Efficiency;
1775  eff[1] = EffError;
1776 }
1777 
1778 void CSCValidation::histoEfficiency(TH1F* readHisto, TH1F* writeHisto) {
1779  std::vector<float> eff(2);
1780  int Nbins = readHisto->GetSize() - 2; //without underflows and overflows
1781  std::vector<float> bins(Nbins);
1782  std::vector<float> Efficiency(Nbins);
1783  std::vector<float> EffError(Nbins);
1784  float Num = 1;
1785  float Den = 1;
1786  for (int i = 0; i < 20; i++) {
1787  Num = readHisto->GetBinContent(i + 1);
1788  Den = readHisto->GetBinContent(i + 21);
1789  getEfficiency(Num, Den, eff);
1790  Efficiency[i] = eff[0];
1791  EffError[i] = eff[1];
1792  writeHisto->SetBinContent(i + 1, Efficiency[i]);
1793  writeHisto->SetBinError(i + 1, EffError[i]);
1794  }
1795 }
1796 
1797 bool CSCValidation::withinSensitiveRegion(LocalPoint localPos,
1798  const std::array<const float, 4>& layerBounds,
1799  int station,
1800  int ring,
1801  float shiftFromEdge,
1802  float shiftFromDeadZone) {
1803  //---- check if it is in a good local region (sensitive area - geometrical and HV boundaries excluded)
1804  bool pass = false;
1805 
1806  float y_center = 0.;
1807  double yUp = layerBounds[3] + y_center;
1808  double yDown = -layerBounds[3] + y_center;
1809  double xBound1Shifted = layerBounds[0] - shiftFromEdge; //
1810  double xBound2Shifted = layerBounds[1] - shiftFromEdge; //
1811  double lineSlope = (yUp - yDown) / (xBound2Shifted - xBound1Shifted);
1812  double lineConst = yUp - lineSlope * xBound2Shifted;
1813  double yBorder = lineSlope * abs(localPos.x()) + lineConst;
1814 
1815  //bool withinChamberOnly = false;// false = "good region"; true - boundaries only
1816  std::vector<float> deadZoneCenter(6);
1817  float cutZone = shiftFromDeadZone; //cm
1818  //---- hardcoded... not good
1819  if (station > 1 && station < 5) {
1820  if (2 == ring) {
1821  deadZoneCenter[0] = -162.48;
1822  deadZoneCenter[1] = -81.8744;
1823  deadZoneCenter[2] = -21.18165;
1824  deadZoneCenter[3] = 39.51105;
1825  deadZoneCenter[4] = 100.2939;
1826  deadZoneCenter[5] = 160.58;
1827 
1828  if (localPos.y() > yBorder &&
1829  ((localPos.y() > deadZoneCenter[0] + cutZone && localPos.y() < deadZoneCenter[1] - cutZone) ||
1830  (localPos.y() > deadZoneCenter[1] + cutZone && localPos.y() < deadZoneCenter[2] - cutZone) ||
1831  (localPos.y() > deadZoneCenter[2] + cutZone && localPos.y() < deadZoneCenter[3] - cutZone) ||
1832  (localPos.y() > deadZoneCenter[3] + cutZone && localPos.y() < deadZoneCenter[4] - cutZone) ||
1833  (localPos.y() > deadZoneCenter[4] + cutZone && localPos.y() < deadZoneCenter[5] - cutZone))) {
1834  pass = true;
1835  }
1836  } else if (1 == ring) {
1837  if (2 == station) {
1838  deadZoneCenter[0] = -95.80;
1839  deadZoneCenter[1] = -27.47;
1840  deadZoneCenter[2] = 33.67;
1841  deadZoneCenter[3] = 90.85;
1842  } else if (3 == station) {
1843  deadZoneCenter[0] = -89.305;
1844  deadZoneCenter[1] = -39.705;
1845  deadZoneCenter[2] = 20.195;
1846  deadZoneCenter[3] = 77.395;
1847  } else if (4 == station) {
1848  deadZoneCenter[0] = -75.645;
1849  deadZoneCenter[1] = -26.055;
1850  deadZoneCenter[2] = 23.855;
1851  deadZoneCenter[3] = 70.575;
1852  }
1853  if (localPos.y() > yBorder &&
1854  ((localPos.y() > deadZoneCenter[0] + cutZone && localPos.y() < deadZoneCenter[1] - cutZone) ||
1855  (localPos.y() > deadZoneCenter[1] + cutZone && localPos.y() < deadZoneCenter[2] - cutZone) ||
1856  (localPos.y() > deadZoneCenter[2] + cutZone && localPos.y() < deadZoneCenter[3] - cutZone))) {
1857  pass = true;
1858  }
1859  }
1860  } else if (1 == station) {
1861  if (3 == ring) {
1862  deadZoneCenter[0] = -83.155;
1863  deadZoneCenter[1] = -22.7401;
1864  deadZoneCenter[2] = 27.86665;
1865  deadZoneCenter[3] = 81.005;
1866  if (localPos.y() > yBorder &&
1867  ((localPos.y() > deadZoneCenter[0] + cutZone && localPos.y() < deadZoneCenter[1] - cutZone) ||
1868  (localPos.y() > deadZoneCenter[1] + cutZone && localPos.y() < deadZoneCenter[2] - cutZone) ||
1869  (localPos.y() > deadZoneCenter[2] + cutZone && localPos.y() < deadZoneCenter[3] - cutZone))) {
1870  pass = true;
1871  }
1872  } else if (2 == ring) {
1873  deadZoneCenter[0] = -86.285;
1874  deadZoneCenter[1] = -32.88305;
1875  deadZoneCenter[2] = 32.867423;
1876  deadZoneCenter[3] = 88.205;
1877  if (localPos.y() > (yBorder) &&
1878  ((localPos.y() > deadZoneCenter[0] + cutZone && localPos.y() < deadZoneCenter[1] - cutZone) ||
1879  (localPos.y() > deadZoneCenter[1] + cutZone && localPos.y() < deadZoneCenter[2] - cutZone) ||
1880  (localPos.y() > deadZoneCenter[2] + cutZone && localPos.y() < deadZoneCenter[3] - cutZone))) {
1881  pass = true;
1882  }
1883  } else {
1884  deadZoneCenter[0] = -81.0;
1885  deadZoneCenter[1] = 81.0;
1886  if (localPos.y() > (yBorder) &&
1887  (localPos.y() > deadZoneCenter[0] + cutZone && localPos.y() < deadZoneCenter[1] - cutZone)) {
1888  pass = true;
1889  }
1890  }
1891  }
1892  return pass;
1893 }
1894 
1895 //---------------------------------------------------------------------------------------
1896 // Given a set of digis, the CSCDetId, and the central strip of your choosing, returns
1897 // the avg. Signal-Pedestal for 6 time bin x 5 strip .
1898 //
1899 // Author: P. Jindal
1900 //---------------------------------------------------------------------------------------
1901 
1902 float CSCValidation::getSignal(const CSCStripDigiCollection& stripdigis, CSCDetId idCS, int centerStrip) {
1903  float SigADC[5];
1904  float TotalADC = 0;
1905  SigADC[0] = 0;
1906  SigADC[1] = 0;
1907  SigADC[2] = 0;
1908  SigADC[3] = 0;
1909  SigADC[4] = 0;
1910 
1911  // Loop over strip digis
1912  CSCStripDigiCollection::DigiRangeIterator sIt;
1913 
1914  for (sIt = stripdigis.begin(); sIt != stripdigis.end(); sIt++) {
1915  CSCDetId id = (CSCDetId)(*sIt).first;
1916  if (id == idCS) {
1917  // First, find the Signal-Pedestal for center strip
1918  std::vector<CSCStripDigi>::const_iterator digiItr = (*sIt).second.first;
1919  std::vector<CSCStripDigi>::const_iterator last = (*sIt).second.second;
1920  for (; digiItr != last; ++digiItr) {
1921  int thisStrip = digiItr->getStrip();
1922  if (thisStrip == (centerStrip)) {
1923  std::vector<int> myADCVals = digiItr->getADCCounts();
1924  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
1925  float thisSignal = (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
1926  SigADC[0] = thisSignal - 6 * thisPedestal;
1927  }
1928  // Now,find the Signal-Pedestal for neighbouring 4 strips
1929  if (thisStrip == (centerStrip + 1)) {
1930  std::vector<int> myADCVals = digiItr->getADCCounts();
1931  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
1932  float thisSignal = (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
1933  SigADC[1] = thisSignal - 6 * thisPedestal;
1934  }
1935  if (thisStrip == (centerStrip + 2)) {
1936  std::vector<int> myADCVals = digiItr->getADCCounts();
1937  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
1938  float thisSignal = (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
1939  SigADC[2] = thisSignal - 6 * thisPedestal;
1940  }
1941  if (thisStrip == (centerStrip - 1)) {
1942  std::vector<int> myADCVals = digiItr->getADCCounts();
1943  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
1944  float thisSignal = (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
1945  SigADC[3] = thisSignal - 6 * thisPedestal;
1946  }
1947  if (thisStrip == (centerStrip - 2)) {
1948  std::vector<int> myADCVals = digiItr->getADCCounts();
1949  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
1950  float thisSignal = (myADCVals[2] + myADCVals[3] + myADCVals[4] + myADCVals[5] + myADCVals[6] + myADCVals[7]);
1951  SigADC[4] = thisSignal - 6 * thisPedestal;
1952  }
1953  }
1954  TotalADC = 0.2 * (SigADC[0] + SigADC[1] + SigADC[2] + SigADC[3] + SigADC[4]);
1955  }
1956  }
1957  return TotalADC;
1958 }
1959 
1960 //---------------------------------------------------------------------------------------
1961 // Look at non-associated recHits
1962 // Author: P. Jindal
1963 //---------------------------------------------------------------------------------------
1964 
1965 void CSCValidation::doNoiseHits(edm::Handle<CSCRecHit2DCollection> recHits,
1966  edm::Handle<CSCSegmentCollection> cscSegments,
1967  edm::ESHandle<CSCGeometry> cscGeom,
1968  edm::Handle<CSCStripDigiCollection> strips) {
1969  CSCRecHit2DCollection::const_iterator recIt;
1970  for (recIt = recHits->begin(); recIt != recHits->end(); recIt++) {
1971  CSCDetId idrec = (CSCDetId)(*recIt).cscDetId();
1972 
1973  //Store the Rechits into a Map
1974  AllRechits.insert(std::pair<CSCDetId, CSCRecHit2D>(idrec, *recIt));
1975 
1976  // Find the strip containing this hit
1977  int centerid = recIt->nStrips() / 2;
1978  int centerStrip = recIt->channels(centerid);
1979 
1980  float rHsignal = getthisSignal(*strips, idrec, centerStrip);
1981  histos->fill1DHist(
1982  rHsignal, "hrHSignal", "Signal in the 4th time bin for centre strip", 1100, -99, 1000, "recHits");
1983  }
1984 
1985  for (CSCSegmentCollection::const_iterator it = cscSegments->begin(); it != cscSegments->end(); it++) {
1986  std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits();
1987  for (std::vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
1988  CSCDetId idRH = (CSCDetId)(*iRH).cscDetId();
1989  LocalPoint lpRH = (*iRH).localPosition();
1990  float xrec = lpRH.x();
1991  float yrec = lpRH.y();
1992  float zrec = lpRH.z();
1993  bool RHalreadyinMap = false;
1994  //Store the rechits associated with segments into a Map
1995  multimap<CSCDetId, CSCRecHit2D>::iterator segRHit;
1996  segRHit = SegRechits.find(idRH);
1997  if (segRHit != SegRechits.end()) {
1998  for (; segRHit != SegRechits.upper_bound(idRH); ++segRHit) {
1999  //for( segRHit = SegRechits.begin(); segRHit != SegRechits.end() ;++segRHit){
2000  LocalPoint lposRH = (segRHit->second).localPosition();
2001  float xpos = lposRH.x();
2002  float ypos = lposRH.y();
2003  float zpos = lposRH.z();
2004  if (xrec == xpos && yrec == ypos && zrec == zpos) {
2005  RHalreadyinMap = true;
2006  //std::cout << " Already exists " <<std ::endl;
2007  break;
2008  }
2009  }
2010  }
2011  if (!RHalreadyinMap) {
2012  SegRechits.insert(std::pair<CSCDetId, CSCRecHit2D>(idRH, *iRH));
2013  }
2014  }
2015  }
2016 
2017  findNonAssociatedRecHits(cscGeom, strips);
2018 }
2019 
2020 //---------------------------------------------------------------------------------------
2021 // Given the list of all rechits and the rechits on a segment finds the rechits
2022 // not associated to a segment and stores in a list
2023 //
2024 //---------------------------------------------------------------------------------------
2025 
2026 void CSCValidation::findNonAssociatedRecHits(edm::ESHandle<CSCGeometry> cscGeom,
2027  edm::Handle<CSCStripDigiCollection> strips) {
2028  for (std::multimap<CSCDetId, CSCRecHit2D>::iterator allRHiter = AllRechits.begin(); allRHiter != AllRechits.end();
2029  ++allRHiter) {
2030  CSCDetId idRH = allRHiter->first;
2031  LocalPoint lpRH = (allRHiter->second).localPosition();
2032  float xrec = lpRH.x();
2033  float yrec = lpRH.y();
2034  float zrec = lpRH.z();
2035 
2036  bool foundmatch = false;
2037  multimap<CSCDetId, CSCRecHit2D>::iterator segRHit;
2038  segRHit = SegRechits.find(idRH);
2039  if (segRHit != SegRechits.end()) {
2040  for (; segRHit != SegRechits.upper_bound(idRH); ++segRHit) {
2041  LocalPoint lposRH = (segRHit->second).localPosition();
2042  float xpos = lposRH.x();
2043  float ypos = lposRH.y();
2044  float zpos = lposRH.z();
2045 
2046  if (xrec == xpos && yrec == ypos && zrec == zpos) {
2047  foundmatch = true;
2048  }
2049 
2050  float d = 0.;
2051  float dclose = 1000.;
2052 
2053  if (!foundmatch) {
2054  d = sqrt(pow(xrec - xpos, 2) + pow(yrec - ypos, 2) + pow(zrec - zpos, 2));
2055  if (d < dclose) {
2056  dclose = d;
2057  if (distRHmap.find((allRHiter->second)) ==
2058  distRHmap.end()) { // entry for rechit does not yet exist, create one
2059  distRHmap.insert(make_pair(allRHiter->second, dclose));
2060  } else {
2061  // we already have an entry for the detid.
2062  distRHmap.erase(allRHiter->second);
2063  distRHmap.insert(
2064  make_pair(allRHiter->second, dclose)); // fill rechits for the segment with the given detid
2065  }
2066  }
2067  }
2068  }
2069  }
2070  if (!foundmatch) {
2071  NonAssociatedRechits.insert(std::pair<CSCDetId, CSCRecHit2D>(idRH, allRHiter->second));
2072  }
2073  }
2074 
2075  for (std::map<CSCRecHit2D, float, ltrh>::iterator iter = distRHmap.begin(); iter != distRHmap.end(); ++iter) {
2076  histos->fill1DHist(iter->second,
2077  "hdistRH",
2078  "Distance of Non Associated RecHit from closest Segment RecHit",
2079  500,
2080  0.,
2081  100.,
2082  "NonAssociatedRechits");
2083  }
2084 
2085  for (std::multimap<CSCDetId, CSCRecHit2D>::iterator iter = NonAssociatedRechits.begin();
2086  iter != NonAssociatedRechits.end();
2087  ++iter) {
2088  CSCDetId idrec = iter->first;
2089  int kEndcap = idrec.endcap();
2090  int cEndcap = idrec.endcap();
2091  if (kEndcap == 2)
2092  cEndcap = -1;
2093  int kRing = idrec.ring();
2094  int kStation = idrec.station();
2095  int kChamber = idrec.chamber();
2096  int kLayer = idrec.layer();
2097 
2098  // Store rechit as a Local Point:
2099  LocalPoint rhitlocal = (iter->second).localPosition();
2100  float xreco = rhitlocal.x();
2101  float yreco = rhitlocal.y();
2102 
2103  // Find the strip containing this hit
2104  int centerid = (iter->second).nStrips() / 2;
2105  int centerStrip = (iter->second).channels(centerid);
2106 
2107  // Find the charge associated with this hit
2108  float rHSumQ = 0;
2109  float sumsides = 0.;
2110  int adcsize = (iter->second).nStrips() * (iter->second).nTimeBins();
2111  for (unsigned int i = 0; i < (iter->second).nStrips(); i++) {
2112  for (unsigned int j = 0; j < (iter->second).nTimeBins() - 1; j++) {
2113  rHSumQ += (iter->second).adcs(i, j);
2114  if (i != 1)
2115  sumsides += (iter->second).adcs(i, j);
2116  }
2117  }
2118 
2119  float rHratioQ = sumsides / rHSumQ;
2120  if (adcsize != 12)
2121  rHratioQ = -99;
2122 
2123  // Get the signal timing of this hit
2124  float rHtime = (iter->second).tpeak() / 50;
2125 
2126  // Get the width of this hit
2127  int rHwidth = getWidth(*strips, idrec, centerStrip);
2128 
2129  // Get pointer to the layer:
2130  const CSCLayer* csclayer = cscGeom->layer(idrec);
2131 
2132  // Transform hit position from local chamber geometry to global CMS geom
2133  GlobalPoint rhitglobal = csclayer->toGlobal(rhitlocal);
2134  float grecx = rhitglobal.x();
2135  float grecy = rhitglobal.y();
2136 
2137  // Simple occupancy variables
2138  int kCodeBroad = cEndcap * (4 * (kStation - 1) + kRing);
2139  int kCodeNarrow = cEndcap * (100 * (kRing - 1) + kChamber);
2140 
2141  //Fill the non-associated rechits parameters in histogram
2142  histos->fill1DHist(
2143  kCodeBroad, "hNARHCodeBroad", "broad scope code for recHits", 33, -16.5, 16.5, "NonAssociatedRechits");
2144  if (kStation == 1)
2145  histos->fill1DHist(kCodeNarrow,
2146  "hNARHCodeNarrow1",
2147  "narrow scope recHit code station 1",
2148  801,
2149  -400.5,
2150  400.5,
2151  "NonAssociatedRechits");
2152  if (kStation == 2)
2153  histos->fill1DHist(kCodeNarrow,
2154  "hNARHCodeNarrow2",
2155  "narrow scope recHit code station 2",
2156  801,
2157  -400.5,
2158  400.5,
2159  "NonAssociatedRechits");
2160  if (kStation == 3)
2161  histos->fill1DHist(kCodeNarrow,
2162  "hNARHCodeNarrow3",
2163  "narrow scope recHit code station 3",
2164  801,
2165  -400.5,
2166  400.5,
2167  "NonAssociatedRechits");
2168  if (kStation == 4)
2169  histos->fill1DHist(kCodeNarrow,
2170  "hNARHCodeNarrow4",
2171  "narrow scope recHit code station 4",
2172  801,
2173  -400.5,
2174  400.5,
2175  "NonAssociatedRechits");
2176  histos->fill1DHistByType(kLayer, "hNARHLayer", "RecHits per Layer", idrec, 8, -0.5, 7.5, "NonAssociatedRechits");
2177  histos->fill1DHistByType(xreco, "hNARHX", "Local X of recHit", idrec, 160, -80., 80., "NonAssociatedRechits");
2178  histos->fill1DHistByType(yreco, "hNARHY", "Local Y of recHit", idrec, 60, -180., 180., "NonAssociatedRechits");
2179  if (kStation == 1 && (kRing == 1 || kRing == 4))
2180  histos->fill1DHistByType(
2181  rHSumQ, "hNARHSumQ", "Sum 3x3 recHit Charge", idrec, 250, 0, 4000, "NonAssociatedRechits");
2182  else
2183  histos->fill1DHistByType(
2184  rHSumQ, "hNARHSumQ", "Sum 3x3 recHit Charge", idrec, 250, 0, 2000, "NonAssociatedRechits");
2185  histos->fill1DHistByType(
2186  rHratioQ, "hNARHRatioQ", "Ratio (Ql+Qr)/Qt)", idrec, 120, -0.1, 1.1, "NonAssociatedRechits");
2187  histos->fill1DHistByType(rHtime, "hNARHTiming", "recHit Timing", idrec, 200, -10, 10, "NonAssociatedRechits");
2188  histos->fill2DHistByStation(grecx,
2189  grecy,
2190  "hNARHGlobal",
2191  "recHit Global Position",
2192  idrec,
2193  400,
2194  -800.,
2195  800.,
2196  400,
2197  -800.,
2198  800.,
2199  "NonAssociatedRechits");
2200  histos->fill1DHistByType(
2201  rHwidth, "hNARHwidth", "width for Non associated recHit", idrec, 21, -0.5, 20.5, "NonAssociatedRechits");
2202  }
2203 
2204  for (std::multimap<CSCDetId, CSCRecHit2D>::iterator iter = SegRechits.begin(); iter != SegRechits.end(); ++iter) {
2205  CSCDetId idrec = iter->first;
2206  int kEndcap = idrec.endcap();
2207  int cEndcap = idrec.endcap();
2208  if (kEndcap == 2)
2209  cEndcap = -1;
2210  int kRing = idrec.ring();
2211  int kStation = idrec.station();
2212  int kChamber = idrec.chamber();
2213  int kLayer = idrec.layer();
2214 
2215  // Store rechit as a Local Point:
2216  LocalPoint rhitlocal = (iter->second).localPosition();
2217  float xreco = rhitlocal.x();
2218  float yreco = rhitlocal.y();
2219 
2220  // Find the strip containing this hit
2221  int centerid = (iter->second).nStrips() / 2;
2222  int centerStrip = (iter->second).channels(centerid);
2223 
2224  // Find the charge associated with this hit
2225 
2226  float rHSumQ = 0;
2227  float sumsides = 0.;
2228  int adcsize = (iter->second).nStrips() * (iter->second).nTimeBins();
2229  for (unsigned int i = 0; i < (iter->second).nStrips(); i++) {
2230  for (unsigned int j = 0; j < (iter->second).nTimeBins() - 1; j++) {
2231  rHSumQ += (iter->second).adcs(i, j);
2232  if (i != 1)
2233  sumsides += (iter->second).adcs(i, j);
2234  }
2235  }
2236 
2237  float rHratioQ = sumsides / rHSumQ;
2238  if (adcsize != 12)
2239  rHratioQ = -99;
2240 
2241  // Get the signal timing of this hit
2242  float rHtime = (iter->second).tpeak() / 50;
2243 
2244  // Get the width of this hit
2245  int rHwidth = getWidth(*strips, idrec, centerStrip);
2246 
2247  // Get pointer to the layer:
2248  const CSCLayer* csclayer = cscGeom->layer(idrec);
2249 
2250  // Transform hit position from local chamber geometry to global CMS geom
2251  GlobalPoint rhitglobal = csclayer->toGlobal(rhitlocal);
2252  float grecx = rhitglobal.x();
2253  float grecy = rhitglobal.y();
2254 
2255  // Simple occupancy variables
2256  int kCodeBroad = cEndcap * (4 * (kStation - 1) + kRing);
2257  int kCodeNarrow = cEndcap * (100 * (kRing - 1) + kChamber);
2258 
2259  //Fill the non-associated rechits global position in histogram
2260  histos->fill1DHist(
2261  kCodeBroad, "hSegRHCodeBroad", "broad scope code for recHits", 33, -16.5, 16.5, "AssociatedRechits");
2262  if (kStation == 1)
2263  histos->fill1DHist(kCodeNarrow,
2264  "hSegRHCodeNarrow1",
2265  "narrow scope recHit code station 1",
2266  801,
2267  -400.5,
2268  400.5,
2269  "AssociatedRechits");
2270  if (kStation == 2)
2271  histos->fill1DHist(kCodeNarrow,
2272  "hSegRHCodeNarrow2",
2273  "narrow scope recHit code station 2",
2274  801,
2275  -400.5,
2276  400.5,
2277  "AssociatedRechits");
2278  if (kStation == 3)
2279  histos->fill1DHist(kCodeNarrow,
2280  "hSegRHCodeNarrow3",
2281  "narrow scope recHit code station 3",
2282  801,
2283  -400.5,
2284  400.5,
2285  "AssociatedRechits");
2286  if (kStation == 4)
2287  histos->fill1DHist(kCodeNarrow,
2288  "hSegRHCodeNarrow4",
2289  "narrow scope recHit code station 4",
2290  801,
2291  -400.5,
2292  400.5,
2293  "AssociatedRechits");
2294  histos->fill1DHistByType(kLayer, "hSegRHLayer", "RecHits per Layer", idrec, 8, -0.5, 7.5, "AssociatedRechits");
2295  histos->fill1DHistByType(xreco, "hSegRHX", "Local X of recHit", idrec, 160, -80., 80., "AssociatedRechits");
2296  histos->fill1DHistByType(yreco, "hSegRHY", "Local Y of recHit", idrec, 60, -180., 180., "AssociatedRechits");
2297  if (kStation == 1 && (kRing == 1 || kRing == 4))
2298  histos->fill1DHistByType(rHSumQ, "hSegRHSumQ", "Sum 3x3 recHit Charge", idrec, 250, 0, 4000, "AssociatedRechits");
2299  else
2300  histos->fill1DHistByType(rHSumQ, "hSegRHSumQ", "Sum 3x3 recHit Charge", idrec, 250, 0, 2000, "AssociatedRechits");
2301  histos->fill1DHistByType(rHratioQ, "hSegRHRatioQ", "Ratio (Ql+Qr)/Qt)", idrec, 120, -0.1, 1.1, "AssociatedRechits");
2302  histos->fill1DHistByType(rHtime, "hSegRHTiming", "recHit Timing", idrec, 200, -10, 10, "AssociatedRechits");
2303  histos->fill2DHistByStation(grecx,
2304  grecy,
2305  "hSegRHGlobal",
2306  "recHit Global Position",
2307  idrec,
2308  400,
2309  -800.,
2310  800.,
2311  400,
2312  -800.,
2313  800.,
2314  "AssociatedRechits");
2315  histos->fill1DHistByType(
2316  rHwidth, "hSegRHwidth", "width for Non associated recHit", idrec, 21, -0.5, 20.5, "AssociatedRechits");
2317  }
2318 
2319  distRHmap.clear();
2320  AllRechits.clear();
2321  SegRechits.clear();
2322  NonAssociatedRechits.clear();
2323 }
2324 
2325 float CSCValidation::getthisSignal(const CSCStripDigiCollection& stripdigis, CSCDetId idRH, int centerStrip) {
2326  // Loop over strip digis responsible for this recHit
2327  CSCStripDigiCollection::DigiRangeIterator sIt;
2328  float thisADC = 0.;
2329  //bool foundRHid = false;
2330  // std::cout<<"iD S/R/C/L = "<<idRH<<" "<<idRH.station()<<"/"<<idRH.ring()<<"/"<<idRH.chamber()<<"/"<<idRH.layer()<<std::endl;
2331  for (sIt = stripdigis.begin(); sIt != stripdigis.end(); sIt++) {
2332  CSCDetId id = (CSCDetId)(*sIt).first;
2333  //std::cout<<"STRIPS: id S/R/C/L = "<<id<<" "<<id.station()<<"/"<<id.ring()<<"/"<<id.chamber()<<"/"<<id.layer()<<std::endl;
2334  if (id == idRH) {
2335  //foundRHid = true;
2336  vector<CSCStripDigi>::const_iterator digiItr = (*sIt).second.first;
2337  vector<CSCStripDigi>::const_iterator last = (*sIt).second.second;
2338  //if(digiItr == last ) {std::cout << " Attention1 :: Size of digi collection is zero " << std::endl;}
2339  int St = idRH.station();
2340  int Rg = idRH.ring();
2341  if (St == 1 && Rg == 4) {
2342  while (centerStrip > 16)
2343  centerStrip -= 16;
2344  }
2345  for (; digiItr != last; ++digiItr) {
2346  int thisStrip = digiItr->getStrip();
2347  //std::cout<<" thisStrip = "<<thisStrip<<" centerStrip = "<<centerStrip<<std::endl;
2348  std::vector<int> myADCVals = digiItr->getADCCounts();
2349  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
2350  float Signal = (float)myADCVals[3];
2351  if (thisStrip == (centerStrip)) {
2352  thisADC = Signal - thisPedestal;
2353  //if(thisADC >= 0. && thisADC <2.) {std::cout << " Attention2 :: The Signal is equal to the pedestal " << std::endl;
2354  //}
2355  //if(thisADC < 0.) {std::cout << " Attention3 :: The Signal is less than the pedestal " << std::endl;
2356  //}
2357  }
2358  if (thisStrip == (centerStrip + 1)) {
2359  std::vector<int> myADCVals = digiItr->getADCCounts();
2360  }
2361  if (thisStrip == (centerStrip - 1)) {
2362  std::vector<int> myADCVals = digiItr->getADCCounts();
2363  }
2364  }
2365  }
2366  }
2367  //if(!foundRHid){std::cout << " Attention4 :: Did not find a matching RH id in the Strip Digi collection " << std::endl;}
2368  return thisADC;
2369 }
2370 
2371 //---------------------------------------------------------------------------------------
2372 //
2373 // Function is meant to take the DetId and center strip number of a recHit and return
2374 // the width in terms of strips
2375 //---------------------------------------------------------------------------------------
2376 
2377 int CSCValidation::getWidth(const CSCStripDigiCollection& stripdigis, CSCDetId idRH, int centerStrip) {
2378  int width = 1;
2379  int widthpos = 0;
2380  int widthneg = 0;
2381 
2382  // Loop over strip digis responsible for this recHit and sum charge
2383  CSCStripDigiCollection::DigiRangeIterator sIt;
2384 
2385  for (sIt = stripdigis.begin(); sIt != stripdigis.end(); sIt++) {
2386  CSCDetId id = (CSCDetId)(*sIt).first;
2387  if (id == idRH) {
2388  std::vector<CSCStripDigi>::const_iterator digiItr = (*sIt).second.first;
2389  std::vector<CSCStripDigi>::const_iterator first = (*sIt).second.first;
2390  std::vector<CSCStripDigi>::const_iterator last = (*sIt).second.second;
2391  std::vector<CSCStripDigi>::const_iterator it = (*sIt).second.first;
2392  std::vector<CSCStripDigi>::const_iterator itr = (*sIt).second.first;
2393  //std::cout << " IDRH " << id <<std::endl;
2394  int St = idRH.station();
2395  int Rg = idRH.ring();
2396  if (St == 1 && Rg == 4) {
2397  while (centerStrip > 16)
2398  centerStrip -= 16;
2399  }
2400  for (; digiItr != last; ++digiItr) {
2401  int thisStrip = digiItr->getStrip();
2402  if (thisStrip == (centerStrip)) {
2403  it = digiItr;
2404  for (; it != last; ++it) {
2405  int strip = it->getStrip();
2406  std::vector<int> myADCVals = it->getADCCounts();
2407  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
2408  if (((float)myADCVals[3] - thisPedestal) < 6 || widthpos == 10 || it == last) {
2409  break;
2410  }
2411  if (strip != centerStrip) {
2412  widthpos += 1;
2413  }
2414  }
2415  itr = digiItr;
2416  for (; itr != first; --itr) {
2417  int strip = itr->getStrip();
2418  std::vector<int> myADCVals = itr->getADCCounts();
2419  float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
2420  if (((float)myADCVals[3] - thisPedestal) < 6 || widthneg == 10 || itr == first) {
2421  break;
2422  }
2423  if (strip != centerStrip) {
2424  widthneg += 1;
2425  }
2426  }
2427  }
2428  }
2429  }
2430  }
2431  //std::cout << "Widthneg - " << widthneg << "Widthpos + " << widthpos << std::endl;
2432  width = width + widthneg + widthpos;
2433  //std::cout << "Width " << width << std::endl;
2434  return width;
2435 }
2436 
2437 //---------------------------------------------------------------------------
2438 // Module for looking at gas gains
2439 // Author N. Terentiev
2440 //---------------------------------------------------------------------------
2441 
2442 void CSCValidation::doGasGain(const CSCWireDigiCollection& wirecltn,
2443  const CSCStripDigiCollection& strpcltn,
2444  const CSCRecHit2DCollection& rechitcltn) {
2445  int channel = 0, mult, wire, layer, idlayer, idchamber, ring;
2446  int wire_strip_rechit_present;
2447  std::string name, title, endcapstr;
2448  ostringstream ss;
2449  CSCIndexer indexer;
2450  std::map<int, int>::iterator intIt;
2451 
2452  m_single_wire_layer.clear();
2453 
2454  if (firstEvent) {
2455  // HV segments, their # and location in terms of wire groups
2456 
2457  m_wire_hvsegm.clear();
2458  std::map<int, std::vector<int> >::iterator intvecIt;
2459  // ME1a ME1b ME1/2 ME1/3 ME2/1 ME2/2 ME3/1 ME3/2 ME4/1 ME4/2
2460  int csctype[10] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
2461  int hvsegm_layer[10] = {1, 1, 3, 3, 3, 5, 3, 5, 3, 5};
2462  int id;
2463  nmbhvsegm.clear();
2464  for (int i = 0; i < 10; i++)
2465  nmbhvsegm.push_back(hvsegm_layer[i]);
2466  // For ME1/1a
2467  std::vector<int> zer_1_1a(49, 0);
2468  id = csctype[0];
2469  if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2470  m_wire_hvsegm[id] = zer_1_1a;
2471  intvecIt = m_wire_hvsegm.find(id);
2472  for (int wire = 1; wire <= 48; wire++)
2473  intvecIt->second[wire] = 1; // Segment 1
2474 
2475  // For ME1/1b
2476  std::vector<int> zer_1_1b(49, 0);
2477  id = csctype[1];
2478  if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2479  m_wire_hvsegm[id] = zer_1_1b;
2480  intvecIt = m_wire_hvsegm.find(id);
2481  for (int wire = 1; wire <= 48; wire++)
2482  intvecIt->second[wire] = 1; // Segment 1
2483 
2484  // For ME1/2
2485  std::vector<int> zer_1_2(65, 0);
2486  id = csctype[2];
2487  if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2488  m_wire_hvsegm[id] = zer_1_2;
2489  intvecIt = m_wire_hvsegm.find(id);
2490  for (int wire = 1; wire <= 24; wire++)
2491  intvecIt->second[wire] = 1; // Segment 1
2492  for (int wire = 25; wire <= 48; wire++)
2493  intvecIt->second[wire] = 2; // Segment 2
2494  for (int wire = 49; wire <= 64; wire++)
2495  intvecIt->second[wire] = 3; // Segment 3
2496 
2497  // For ME1/3
2498  std::vector<int> zer_1_3(33, 0);
2499  id = csctype[3];
2500  if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2501  m_wire_hvsegm[id] = zer_1_3;
2502  intvecIt = m_wire_hvsegm.find(id);
2503  for (int wire = 1; wire <= 12; wire++)
2504  intvecIt->second[wire] = 1; // Segment 1
2505  for (int wire = 13; wire <= 22; wire++)
2506  intvecIt->second[wire] = 2; // Segment 2
2507  for (int wire = 23; wire <= 32; wire++)
2508  intvecIt->second[wire] = 3; // Segment 3
2509 
2510  // For ME2/1
2511  std::vector<int> zer_2_1(113, 0);
2512  id = csctype[4];
2513  if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2514  m_wire_hvsegm[id] = zer_2_1;
2515  intvecIt = m_wire_hvsegm.find(id);
2516  for (int wire = 1; wire <= 44; wire++)
2517  intvecIt->second[wire] = 1; // Segment 1
2518  for (int wire = 45; wire <= 80; wire++)
2519  intvecIt->second[wire] = 2; // Segment 2
2520  for (int wire = 81; wire <= 112; wire++)
2521  intvecIt->second[wire] = 3; // Segment 3
2522 
2523  // For ME2/2
2524  std::vector<int> zer_2_2(65, 0);
2525  id = csctype[5];
2526  if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2527  m_wire_hvsegm[id] = zer_2_2;
2528  intvecIt = m_wire_hvsegm.find(id);
2529  for (int wire = 1; wire <= 16; wire++)
2530  intvecIt->second[wire] = 1; // Segment 1
2531  for (int wire = 17; wire <= 28; wire++)
2532  intvecIt->second[wire] = 2; // Segment 2
2533  for (int wire = 29; wire <= 40; wire++)
2534  intvecIt->second[wire] = 3; // Segment 3
2535  for (int wire = 41; wire <= 52; wire++)
2536  intvecIt->second[wire] = 4; // Segment 4
2537  for (int wire = 53; wire <= 64; wire++)
2538  intvecIt->second[wire] = 5; // Segment 5
2539 
2540  // For ME3/1
2541  std::vector<int> zer_3_1(97, 0);
2542  id = csctype[6];
2543  if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2544  m_wire_hvsegm[id] = zer_3_1;
2545  intvecIt = m_wire_hvsegm.find(id);
2546  for (int wire = 1; wire <= 32; wire++)
2547  intvecIt->second[wire] = 1; // Segment 1
2548  for (int wire = 33; wire <= 64; wire++)
2549  intvecIt->second[wire] = 2; // Segment 2
2550  for (int wire = 65; wire <= 96; wire++)
2551  intvecIt->second[wire] = 3; // Segment 3
2552 
2553  // For ME3/2
2554  std::vector<int> zer_3_2(65, 0);
2555  id = csctype[7];
2556  if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2557  m_wire_hvsegm[id] = zer_3_2;
2558  intvecIt = m_wire_hvsegm.find(id);
2559  for (int wire = 1; wire <= 16; wire++)
2560  intvecIt->second[wire] = 1; // Segment 1
2561  for (int wire = 17; wire <= 28; wire++)
2562  intvecIt->second[wire] = 2; // Segment 2
2563  for (int wire = 29; wire <= 40; wire++)
2564  intvecIt->second[wire] = 3; // Segment 3
2565  for (int wire = 41; wire <= 52; wire++)
2566  intvecIt->second[wire] = 4; // Segment 4
2567  for (int wire = 53; wire <= 64; wire++)
2568  intvecIt->second[wire] = 5; // Segment 5
2569 
2570  // For ME4/1
2571  std::vector<int> zer_4_1(97, 0);
2572  id = csctype[8];
2573  if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2574  m_wire_hvsegm[id] = zer_4_1;
2575  intvecIt = m_wire_hvsegm.find(id);
2576  for (int wire = 1; wire <= 32; wire++)
2577  intvecIt->second[wire] = 1; // Segment 1
2578  for (int wire = 33; wire <= 64; wire++)
2579  intvecIt->second[wire] = 2; // Segment 2
2580  for (int wire = 65; wire <= 96; wire++)
2581  intvecIt->second[wire] = 3; // Segment 3
2582 
2583  // For ME4/2
2584  std::vector<int> zer_4_2(65, 0);
2585  id = csctype[9];
2586  if (m_wire_hvsegm.find(id) == m_wire_hvsegm.end())
2587  m_wire_hvsegm[id] = zer_4_2;
2588  intvecIt = m_wire_hvsegm.find(id);
2589  for (int wire = 1; wire <= 16; wire++)
2590  intvecIt->second[wire] = 1; // Segment 1
2591  for (int wire = 17; wire <= 28; wire++)
2592  intvecIt->second[wire] = 2; // Segment 2
2593  for (int wire = 29; wire <= 40; wire++)
2594  intvecIt->second[wire] = 3; // Segment 3
2595  for (int wire = 41; wire <= 52; wire++)
2596  intvecIt->second[wire] = 4; // Segment 4
2597  for (int wire = 53; wire <= 64; wire++)
2598  intvecIt->second[wire] = 5; // Segment 5
2599 
2600  } // end of if(nEventsAnalyzed==1)
2601 
2602  // do wires, strips and rechits present?
2603  wire_strip_rechit_present = 0;
2604  if (wirecltn.begin() != wirecltn.end())
2605  wire_strip_rechit_present = wire_strip_rechit_present + 1;
2606  if (strpcltn.begin() != strpcltn.end())
2607  wire_strip_rechit_present = wire_strip_rechit_present + 2;
2608  if (rechitcltn.begin() != rechitcltn.end())
2609  wire_strip_rechit_present = wire_strip_rechit_present + 4;
2610 
2611  if (wire_strip_rechit_present == 7) {
2612  // std::cout<<"Event "<<nEventsAnalyzed<<std::endl;
2613  // std::cout<<std::endl;
2614 
2615  // cycle on wire collection for all CSC to select single wire hit layers
2616  CSCWireDigiCollection::DigiRangeIterator wiredetUnitIt;
2617 
2618  for (wiredetUnitIt = wirecltn.begin(); wiredetUnitIt != wirecltn.end(); ++wiredetUnitIt) {
2619  const CSCDetId id = (*wiredetUnitIt).first;
2620  idlayer = indexer.dbIndex(id, channel);
2621  // looping in the layer of given CSC
2622  mult = 0;
2623  wire = 0;
2624  const CSCWireDigiCollection::Range& range = (*wiredetUnitIt).second;
2625  for (CSCWireDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
2626  wire = (*digiIt).getWireGroup();
2627  mult++;
2628  } // end of digis loop in layer
2629 
2630  // select layers with single wire hit
2631  if (mult == 1) {
2632  if (m_single_wire_layer.find(idlayer) == m_single_wire_layer.end())
2633  m_single_wire_layer[idlayer] = wire;
2634  } // end of if(mult==1)
2635  } // end of cycle on detUnit
2636 
2637  // Looping thru rechit collection
2638  CSCRecHit2DCollection::const_iterator recIt;
2639  CSCRecHit2D::ADCContainer m_adc;
2640 
2641  for (recIt = rechitcltn.begin(); recIt != rechitcltn.end(); ++recIt) {
2642  CSCDetId id = (CSCDetId)(*recIt).cscDetId();
2643  idlayer = indexer.dbIndex(id, channel);
2644  idchamber = idlayer / 10;
2645  layer = id.layer();
2646  // select layer with single wire rechit
2647  if (m_single_wire_layer.find(idlayer) != m_single_wire_layer.end()) {
2648  if (recIt->nStrips() == 3) {
2649  // get 3X3 ADC Sum
2650  unsigned int binmx = 0;
2651  float adcmax = 0.0;
2652 
2653  for (unsigned int i = 0; i < recIt->nStrips(); i++)
2654  for (unsigned int j = 0; j < recIt->nTimeBins(); j++)
2655  if (recIt->adcs(i, j) > adcmax) {
2656  adcmax = recIt->adcs(i, j);
2657  binmx = j;
2658  }
2659 
2660  float adc_3_3_sum = 0.0;
2661  //well, this really only works for 3 strips in readout - not sure the right fix for general case
2662  for (unsigned int i = 0; i < recIt->nStrips(); i++)
2663  for (unsigned int j = binmx - 1; j <= binmx + 1; j++)
2664  adc_3_3_sum += recIt->adcs(i, j);
2665 
2666  if (adc_3_3_sum > 0.0 && adc_3_3_sum < 2000.0) {
2667  // temporary fix for ME1/1a to avoid triple entries
2668  int flag = 0;
2669  if (id.station() == 1 && id.ring() == 4 && recIt->channels(1) > 16)
2670  flag = 1;
2671  // end of temporary fix
2672  if (flag == 0) {
2673  wire = m_single_wire_layer[idlayer];
2674  int chambertype = id.iChamberType(id.station(), id.ring());
2675  int hvsgmtnmb = m_wire_hvsegm[chambertype][wire];
2676  int nmbofhvsegm = nmbhvsegm[chambertype - 1];
2677  int location = (layer - 1) * nmbofhvsegm + hvsgmtnmb;
2678 
2679  ss << "gas_gain_rechit_adc_3_3_sum_location_ME_" << idchamber;
2680  name = ss.str();
2681  ss.str("");
2682  if (id.endcap() == 1)
2683  endcapstr = "+";
2684  ring = id.ring();
2685  if (id.station() == 1 && id.ring() == 4)
2686  ring = 1;
2687  if (id.endcap() == 2)
2688  endcapstr = "-";
2689  ss << "Gas Gain Rechit ADC3X3 Sum ME" << endcapstr << id.station() << "/" << ring << "/" << id.chamber();
2690  title = ss.str();
2691  ss.str("");
2692  float x = location;
2693  float y = adc_3_3_sum;
2694  histos->fill2DHist(x, y, name, title, 30, 1.0, 31.0, 50, 0.0, 2000.0, "GasGain");
2695 
2696  /*
2697  std::cout<<idchamber<<" "<<id.station()<<" "<<id.ring()<<" "
2698  <<id.chamber()<<" "<<layer<<" "<< wire<<" "<<m_strip[1]<<" "<<
2699  chambertype<<" "<< hvsgmtnmb<<" "<< nmbofhvsegm<<" "<<
2700  location<<" "<<adc_3_3_sum<<std::endl;
2701  */
2702  } // end of if flag==0
2703  } // end if(adcsum>0.0 && adcsum<2000.0)
2704  } // end of if if(m_strip.size()==3
2705  } // end of if single wire
2706  } // end of looping thru rechit collection
2707  } // end of if wire and strip and rechit present
2708 }
2709 
2710 //---------------------------------------------------------------------------
2711 // Module for looking at AFEB Timing
2712 // Author N. Terentiev
2713 //---------------------------------------------------------------------------
2714 
2715 void CSCValidation::doAFEBTiming(const CSCWireDigiCollection& wirecltn) {
2716  ostringstream ss;
2717  std::string name, title, endcapstr;
2718  float x, y;
2719  int wire, wiretbin, nmbwiretbin, layer, afeb, idlayer, idchamber;
2720  int channel = 0; // for CSCIndexer::dbIndex(id, channel); irrelevant here
2721  CSCIndexer indexer;
2722 
2723  if (wirecltn.begin() != wirecltn.end()) {
2724  //std::cout<<std::endl;
2725  //std::cout<<"Event "<<nEventsAnalyzed<<std::endl;
2726  //std::cout<<std::endl;
2727 
2728  // cycle on wire collection for all CSC
2729  CSCWireDigiCollection::DigiRangeIterator wiredetUnitIt;
2730  for (wiredetUnitIt = wirecltn.begin(); wiredetUnitIt != wirecltn.end(); ++wiredetUnitIt) {
2731  const CSCDetId id = (*wiredetUnitIt).first;
2732  idlayer = indexer.dbIndex(id, channel);
2733  idchamber = idlayer / 10;
2734  layer = id.layer();
2735 
2736  if (id.endcap() == 1)
2737  endcapstr = "+";
2738  if (id.endcap() == 2)
2739  endcapstr = "-";
2740 
2741  // looping in the layer of given CSC
2742 
2743  const CSCWireDigiCollection::Range& range = (*wiredetUnitIt).second;
2744  for (CSCWireDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
2745  wire = (*digiIt).getWireGroup();
2746  wiretbin = (*digiIt).getTimeBin();
2747  nmbwiretbin = (*digiIt).getTimeBinsOn().size();
2748  afeb = 3 * ((wire - 1) / 8) + (layer + 1) / 2;
2749 
2750  // Anode wire group time bin vs afeb for each CSC
2751  x = afeb;
2752  y = wiretbin;
2753  ss << "afeb_time_bin_vs_afeb_occupancy_ME_" << idchamber;
2754  name = ss.str();
2755  ss.str("");
2756  ss << "Time Bin vs AFEB Occupancy ME" << endcapstr << id.station() << "/" << id.ring() << "/" << id.chamber();
2757  title = ss.str();
2758  ss.str("");
2759  histos->fill2DHist(x, y, name, title, 42, 1., 43., 16, 0., 16., "AFEBTiming");
2760 
2761  // Number of anode wire group time bin vs afeb for each CSC
2762  x = afeb;
2763  y = nmbwiretbin;
2764  ss << "nmb_afeb_time_bins_vs_afeb_ME_" << idchamber;
2765  name = ss.str();
2766  ss.str("");
2767  ss << "Number of Time Bins vs AFEB ME" << endcapstr << id.station() << "/" << id.ring() << "/" << id.chamber();
2768  title = ss.str();
2769  ss.str("");
2770  histos->fill2DHist(x, y, name, title, 42, 1., 43., 16, 0., 16., "AFEBTiming");
2771 
2772  } // end of digis loop in layer
2773  } // end of wire collection loop
2774  } // end of if(wirecltn.begin() != wirecltn.end())
2775 }
2776 
2777 //---------------------------------------------------------------------------
2778 // Module for looking at Comparitor Timing
2779 // Author N. Terentiev
2780 //---------------------------------------------------------------------------
2781 
2782 void CSCValidation::doCompTiming(const CSCComparatorDigiCollection& compars) {
2783  ostringstream ss;
2784  std::string name, title, endcap;
2785  float x, y;
2786  int strip, tbin, cfeb, idlayer, idchamber;
2787  int channel = 0; // for CSCIndexer::dbIndex(id, channel); irrelevant here
2788  CSCIndexer indexer;
2789 
2790  if (compars.begin() != compars.end()) {
2791  //std::cout<<std::endl;
2792  //std::cout<<"Event "<<nEventsAnalyzed<<std::endl;
2793  //std::cout<<std::endl;
2794 
2795  // cycle on comparators collection for all CSC
2796  CSCComparatorDigiCollection::DigiRangeIterator compdetUnitIt;
2797  for (compdetUnitIt = compars.begin(); compdetUnitIt != compars.end(); ++compdetUnitIt) {
2798  const CSCDetId id = (*compdetUnitIt).first;
2799  idlayer = indexer.dbIndex(id, channel); // channel irrelevant here
2800  idchamber = idlayer / 10;
2801 
2802  if (id.endcap() == 1)
2803  endcap = "+";
2804  if (id.endcap() == 2)
2805  endcap = "-";
2806  // looping in the layer of given CSC
2807  const CSCComparatorDigiCollection::Range& range = (*compdetUnitIt).second;
2808  for (CSCComparatorDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
2809  strip = (*digiIt).getStrip();
2810  /*
2811  if(id.station()==1 && (id.ring()==1 || id.ring()==4))
2812  std::cout<<idchamber<<" "<<id.station()<<" "<<id.ring()<<" "
2813  <<strip <<std::endl;
2814  */
2815  indexer.dbIndex(id, strip); // strips 1-16 of ME1/1a
2816  // become strips 65-80 of ME1/1
2817  tbin = (*digiIt).getTimeBin();
2818  cfeb = (strip - 1) / 16 + 1;
2819 
2820  // time bin vs cfeb for each CSC
2821 
2822  x = cfeb;
2823  y = tbin;
2824  ss << "comp_time_bin_vs_cfeb_occupancy_ME_" << idchamber;
2825  name = ss.str();
2826  ss.str("");
2827  ss << "Comparator Time Bin vs CFEB Occupancy ME" << endcap << id.station() << "/" << id.ring() << "/"
2828  << id.chamber();
2829  title = ss.str();
2830  ss.str("");
2831  histos->fill2DHist(x, y, name, title, 5, 1., 6., 16, 0., 16., "CompTiming");
2832 
2833  } // end of digis loop in layer
2834  } // end of collection loop
2835  } // end of if(compars.begin() !=compars.end())
2836 }
2837 
2838 //---------------------------------------------------------------------------
2839 // Module for looking at Strip Timing
2840 // Author N. Terentiev
2841 //---------------------------------------------------------------------------
2842 
2843 void CSCValidation::doADCTiming(const CSCRecHit2DCollection& rechitcltn) {
2844  float adc_3_3_sum, adc_3_3_wtbin, x, y;
2845  int cfeb, idchamber, ring;
2846 
2847  std::string name, title, endcapstr;
2848  ostringstream ss;
2849  std::vector<float> zer(6, 0.0);
2850 
2851  CSCIndexer indexer;
2852  std::map<int, int>::iterator intIt;
2853 
2854  if (rechitcltn.begin() != rechitcltn.end()) {
2855  // std::cout<<"Event "<<nEventsAnalyzed <<std::endl;
2856 
2857  // Looping thru rechit collection
2858  CSCRecHit2DCollection::const_iterator recIt;
2859  CSCRecHit2D::ADCContainer m_adc;
2860  for (recIt = rechitcltn.begin(); recIt != rechitcltn.end(); ++recIt) {
2861  CSCDetId id = (CSCDetId)(*recIt).cscDetId();
2862  // getting strips comprising rechit
2863  if (recIt->nStrips() == 3) {
2864  // get 3X3 ADC Sum
2865  // get 3X3 ADC Sum
2866  unsigned int binmx = 0;
2867  float adcmax = 0.0;
2868 
2869  for (unsigned int i = 0; i < recIt->nStrips(); i++)
2870  for (unsigned int j = 0; j < recIt->nTimeBins(); j++)
2871  if (recIt->adcs(i, j) > adcmax) {
2872  adcmax = recIt->adcs(i, j);
2873  binmx = j;
2874  }
2875 
2876  adc_3_3_sum = 0.0;
2877  //well, this really only works for 3 strips in readout - not sure the right fix for general case
2878  for (unsigned int i = 0; i < recIt->nStrips(); i++)
2879  for (unsigned int j = binmx - 1; j <= binmx + 1; j++)
2880  adc_3_3_sum += recIt->adcs(i, j);
2881 
2882  // ADC weighted time bin
2883  if (adc_3_3_sum > 100.0) {
2884  int centerStrip = recIt->channels(1); //take central from 3 strips;
2885  // temporary fix
2886  int flag = 0;
2887  if (id.station() == 1 && id.ring() == 4 && centerStrip > 16)
2888  flag = 1;
2889  // end of temporary fix
2890  if (flag == 0) {
2891  adc_3_3_wtbin = (*recIt).tpeak() / 50; //getTiming(strpcltn, id, centerStrip);
2892  idchamber = indexer.dbIndex(id, centerStrip) / 10; //strips 1-16 ME1/1a
2893  // become strips 65-80 ME1/1 !!!
2894  /*
2895  if(id.station()==1 && (id.ring()==1 || id.ring()==4))
2896  std::cout<<idchamber<<" "<<id.station()<<" "<<id.ring()<<" "<<m_strip[1]<<" "<<
2897  " "<<centerStrip<<
2898  " "<<adc_3_3_wtbin<<" "<<adc_3_3_sum<<std::endl;
2899  */
2900  ss << "adc_3_3_weight_time_bin_vs_cfeb_occupancy_ME_" << idchamber;
2901  name = ss.str();
2902  ss.str("");
2903 
2904  std::string endcapstr;
2905  if (id.endcap() == 1)
2906  endcapstr = "+";
2907  if (id.endcap() == 2)
2908  endcapstr = "-";
2909  ring = id.ring();
2910  if (id.ring() == 4)
2911  ring = 1;
2912  ss << "ADC 3X3 Weighted Time Bin vs CFEB Occupancy ME" << endcapstr << id.station() << "/" << ring << "/"
2913  << id.chamber();
2914  title = ss.str();
2915  ss.str("");
2916 
2917  cfeb = (centerStrip - 1) / 16 + 1;
2918  x = cfeb;
2919  y = adc_3_3_wtbin;
2920  histos->fill2DHist(x, y, name, title, 5, 1., 6., 80, -8., 8., "ADCTiming");
2921  } // end of if flag==0
2922  } // end of if (adc_3_3_sum > 100.0)
2923  } // end of if if(m_strip.size()==3
2924  } // end of the pass thru CSCRecHit2DCollection
2925  } // end of if (rechitcltn.begin() != rechitcltn.end())
2926 }
2927 
2928 //---------------------------------------------------------------------------------------
2929 // Construct histograms for monitoring the trigger and offline timing
2930 // Author: A. Deisher
2931 //---------------------------------------------------------------------------------------
2932 
2933 void CSCValidation::doTimeMonitoring(edm::Handle<CSCRecHit2DCollection> recHits,
2934  edm::Handle<CSCSegmentCollection> cscSegments,
2935  edm::Handle<CSCALCTDigiCollection> alcts,
2936  edm::Handle<CSCCLCTDigiCollection> clcts,
2937  edm::Handle<CSCCorrelatedLCTDigiCollection> correlatedlcts,
2938  edm::Handle<L1MuGMTReadoutCollection> pCollection,
2939  edm::ESHandle<CSCGeometry> cscGeom,
2940  const edm::EventSetup& eventSetup,
2941  const edm::Event& event) {
2942  map<CSCDetId, float> segment_median_map; //structure for storing the median time for segments in a chamber
2943  map<CSCDetId, GlobalPoint> segment_position_map; //structure for storing the global position for segments in a chamber
2944 
2945  // -----------------------
2946  // loop over segments
2947  // -----------------------
2948  int iSegment = 0;
2949  for (CSCSegmentCollection::const_iterator dSiter = cscSegments->begin(); dSiter != cscSegments->end(); dSiter++) {
2950  iSegment++;
2951 
2952  CSCDetId id = (CSCDetId)(*dSiter).cscDetId();
2953  LocalPoint localPos = (*dSiter).localPosition();
2954  GlobalPoint globalPosition = GlobalPoint(0.0, 0.0, 0.0);
2955  const CSCChamber* cscchamber = cscGeom->chamber(id);
2956  if (cscchamber) {
2957  globalPosition = cscchamber->toGlobal(localPos);
2958  }
2959 
2960  // try to get the CSC recHits that contribute to this segment.
2961  std::vector<CSCRecHit2D> theseRecHits = (*dSiter).specificRecHits();
2962  int nRH = (*dSiter).nRecHits();
2963  if (nRH < 4)
2964  continue;
2965 
2966  //Store the recHit times of a segment in a vector for later sorting
2967  vector<float> non_zero;
2968 
2969  for (vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
2970  non_zero.push_back(iRH->tpeak());
2971 
2972  } // end rechit loop
2973 
2974  //Sort the vector of hit times for this segment and average the center two
2975  sort(non_zero.begin(), non_zero.end());
2976  int middle_index = non_zero.size() / 2;
2977  float average_two = (non_zero.at(middle_index - 1) + non_zero.at(middle_index)) / 2.;
2978  if (non_zero.size() % 2)
2979  average_two = non_zero.at(middle_index);
2980 
2981  //If we've vetoed events with multiple segments per chamber, this should never overwrite informations
2982  segment_median_map[id] = average_two;
2983  segment_position_map[id] = globalPosition;
2984 
2985  double distToIP = sqrt(globalPosition.x() * globalPosition.x() + globalPosition.y() * globalPosition.y() +
2986  globalPosition.z() * globalPosition.z());
2987 
2988  histos->fillProfile(chamberSerial(id),
2989  average_two,
2990  "timeChamber",
2991  "Segment mean time",
2992  601,
2993  -0.5,
2994  600.5,
2995  -400.,
2996  400.,
2997  "TimeMonitoring");
2998  histos->fillProfileByType(id.chamber(),
2999  average_two,
3000  "timeChamberByType",
3001  "Segment mean time by chamber",
3002  id,
3003  36,
3004  0.5,
3005  36.5,
3006  -400,
3007  400.,
3008  "TimeMonitoring");
3009  histos->fill2DHist(distToIP,
3010  average_two,
3011  "seg_time_vs_distToIP",
3012  "Segment time vs. Distance to IP",
3013  80,
3014  600.,
3015  1400.,
3016  800,
3017  -400,
3018  400.,
3019  "TimeMonitoring");
3020  histos->fill2DHist(globalPosition.z(),
3021  average_two,
3022  "seg_time_vs_globZ",
3023  "Segment time vs. z position",
3024  240,
3025  -1200,
3026  1200,
3027  800,
3028  -400.,
3029  400.,
3030  "TimeMonitoring");
3031  histos->fill2DHist(fabs(globalPosition.z()),
3032  average_two,
3033  "seg_time_vs_absglobZ",
3034  "Segment time vs. abs(z position)",
3035  120,
3036  0.,
3037  1200.,
3038  800,
3039  -400.,
3040  400.,
3041  "TimeMonitoring");
3042 
3043  } //end segment loop
3044 
3045  //Now that the information for each segment we're interest in is stored, it is time to go through the pairs and make plots
3046  map<CSCDetId, float>::const_iterator it_outer; //for the outer loop
3047  map<CSCDetId, float>::const_iterator it_inner; //for the nested inner loop
3048  for (it_outer = segment_median_map.begin(); it_outer != segment_median_map.end(); it_outer++) {
3049  CSCDetId id_outer = it_outer->first;
3050  float t_outer = it_outer->second;
3051 
3052  //begin the inner loop
3053  for (it_inner = segment_median_map.begin(); it_inner != segment_median_map.end(); it_inner++) {
3054  CSCDetId id_inner = it_inner->first;
3055  float t_inner = it_inner->second;
3056 
3057  // we're looking at ordered pairs, so combinations will be double counted
3058  // (chamber a, chamber b) will be counted as well as (chamber b, chamber a)
3059  // We will avoid (chamber a, chamber a) with the following line
3060  if (chamberSerial(id_outer) == chamberSerial(id_inner))
3061  continue;
3062 
3063  // Calculate expected TOF (in ns units)
3064  // GlobalPoint gp_outer = segment_position_map.find(id_outer)->second;
3065  // GlobalPoint gp_inner = segment_position_map.find(id_inner)->second;
3066  // GlobalVector flight = gp_outer - gp_inner; //in cm
3067  // float TOF = flight.mag()/30.0; //to ns
3068 
3069  //Plot t(ME+) - t(ME-) for chamber pairs in the same stations and rings but opposite endcaps
3070  if (id_outer.endcap() == 1 && id_inner.endcap() == 2 && id_outer.station() == id_inner.station() &&
3071  id_outer.ring() == id_inner.ring()) {
3072  histos->fill1DHist(t_outer - t_inner,
3073  "diff_opposite_endcaps",
3074  "#Delta t [ME+]-[ME-] for chambers in same station and ring",
3075  800,
3076  -400.,
3077  400.,
3078  "TimeMonitoring");
3079  histos->fill1DHistByType(t_outer - t_inner,
3080  "diff_opposite_endcaps_byType",
3081  "#Delta t [ME+]-[ME-] for chambers in same station and ring",
3082  id_outer,
3083  800,
3084  -400.,
3085  400.,
3086  "TimeMonitoring");
3087  }
3088 
3089  } //end inner loop of segment pairs
3090  } //end outer loop of segment pairs
3091 
3092  //if the digis, return here
3093  if (!useDigis)
3094  return;
3095 
3096  //looking for the global trigger number
3097  vector<L1MuGMTReadoutRecord> L1Mrec = pCollection->getRecords();
3098  vector<L1MuGMTReadoutRecord>::const_iterator igmtrr;
3099  int L1GMT_BXN = -100;
3100  bool has_CSCTrigger = false;
3101  bool has_beamHaloTrigger = false;
3102  for (igmtrr = L1Mrec.begin(); igmtrr != L1Mrec.end(); igmtrr++) {
3103  std::vector<L1MuRegionalCand>::const_iterator iter1;
3104  std::vector<L1MuRegionalCand> rmc;
3105  // CSC
3106  int icsc = 0;
3107  rmc = igmtrr->getCSCCands();
3108  for (iter1 = rmc.begin(); iter1 != rmc.end(); iter1++) {
3109  if (!(*iter1).empty()) {
3110  icsc++;
3111  int kQuality = (*iter1).quality(); // kQuality = 1 means beam halo
3112  if (kQuality == 1)
3113  has_beamHaloTrigger = true;
3114  }
3115  }
3116  if (igmtrr->getBxInEvent() == 0 && icsc > 0) {
3117  //printf("L1 CSCCands exist. L1MuGMTReadoutRecord BXN = %d \n", igmtrr->getBxNr());
3118  L1GMT_BXN = igmtrr->getBxNr();
3119  has_CSCTrigger = true;
3120  } else if (igmtrr->getBxInEvent() == 0) {
3121  //printf("L1 CSCCands do not exist. L1MuGMTReadoutRecord BXN = %d \n", igmtrr->getBxNr());
3122  L1GMT_BXN = igmtrr->getBxNr();
3123  }
3124  }
3125 
3126  // *************************************************
3127  // *** ALCT Digis **********************************
3128  // *************************************************
3129 
3130  int n_alcts = 0;
3131  map<CSCDetId, int> ALCT_KeyWG_map; //structure for storing the key wire group for the first ALCT for each chamber
3132  for (CSCALCTDigiCollection::DigiRangeIterator j = alcts->begin(); j != alcts->end(); j++) {
3133  const CSCALCTDigiCollection::Range& range = (*j).second;
3134  const CSCDetId& idALCT = (*j).first;
3135  for (CSCALCTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
3136  // Valid digi in the chamber (or in neighbouring chamber)
3137  if ((*digiIt).isValid()) {
3138  n_alcts++;
3139  histos->fill1DHist((*digiIt).getBX(), "ALCT_getBX", "ALCT.getBX()", 11, -0.5, 10.5, "TimeMonitoring");
3140  histos->fill1DHist(
3141  (*digiIt).getFullBX(), "ALCT_getFullBX", "ALCT.getFullBX()", 3601, -0.5, 3600.5, "TimeMonitoring");
3142  //if we don't already have digi information stored for this chamber, then we fill it
3143  if (ALCT_KeyWG_map.find(idALCT.chamberId()) == ALCT_KeyWG_map.end()) {
3144  ALCT_KeyWG_map[idALCT.chamberId()] = (*digiIt).getKeyWG();
3145  //printf("I did fill ALCT info for Chamber %d %d %d %d \n",idALCT.chamberId().endcap(), idALCT.chamberId().station(), idALCT.chamberId().ring(), idALCT.chamberId().chamber());
3146  }
3147  }
3148  }
3149  }
3150 
3151  // *************************************************
3152  // *** CLCT Digis **********************************
3153  // *************************************************
3154  int n_clcts = 0;
3155  map<CSCDetId, int> CLCT_getFullBx_map; //structure for storing the pretrigger bxn for the first CLCT for each chamber
3156  for (CSCCLCTDigiCollection::DigiRangeIterator j = clcts->begin(); j != clcts->end(); j++) {
3157  const CSCCLCTDigiCollection::Range& range = (*j).second;
3158  const CSCDetId& idCLCT = (*j).first;
3159  for (CSCCLCTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
3160  // Valid digi in the chamber (or in neighbouring chamber)
3161  if ((*digiIt).isValid()) {
3162  n_clcts++;
3163  histos->fill1DHist((*digiIt).getBX(), "CLCT_getBX", "CLCT.getBX()", 11, -0.5, 10.5, "TimeMonitoring");
3164  histos->fill1DHist(
3165  (*digiIt).getFullBX(), "CLCT_getFullBX", "CLCT.getFullBX()", 3601, -0.5, 3600.5, "TimeMonitoring");
3166  //if we don't already have digi information stored for this chamber, then we fill it
3167  if (CLCT_getFullBx_map.find(idCLCT.chamberId()) == CLCT_getFullBx_map.end()) {
3168  CLCT_getFullBx_map[idCLCT.chamberId()] = (*digiIt).getFullBX();
3169  //printf("I did fill CLCT info for Chamber %d %d %d %d \n",idCLCT.chamberId().endcap(), idCLCT.chamberId().station(), idCLCT.chamberId().ring(), idCLCT.chamberId().chamber());
3170  }
3171  }
3172  }
3173  }
3174 
3175  // *************************************************
3176  // *** CorrelatedLCT Digis *************************
3177  // *************************************************
3178  int n_correlatedlcts = 0;
3179  for (CSCCorrelatedLCTDigiCollection::DigiRangeIterator j = correlatedlcts->begin(); j != correlatedlcts->end(); j++) {
3180  const CSCCorrelatedLCTDigiCollection::Range& range = (*j).second;
3181  for (CSCCorrelatedLCTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
3182  if ((*digiIt).isValid()) {
3183  n_correlatedlcts++;
3184  histos->fill1DHist(
3185  (*digiIt).getBX(), "CorrelatedLCTS_getBX", "CorrelatedLCT.getBX()", 11, -0.5, 10.5, "TimeMonitoring");
3186  }
3187  }
3188  }
3189 
3190  int nRecHits = recHits->size();
3191  int nSegments = cscSegments->size();
3192  if (has_CSCTrigger) {
3193  histos->fill1DHist(L1GMT_BXN, "BX_L1CSCCand", "BX of L1 CSC Cand", 4001, -0.5, 4000.5, "TimeMonitoring");
3194  histos->fill2DHist(L1GMT_BXN,
3195  n_alcts,
3196  "n_ALCTs_v_BX_L1CSCCand",
3197  "Number of ALCTs vs. BX of L1 CSC Cand",
3198  4001,
3199  -0.5,
3200  4000.5,
3201  51,
3202  -0.5,
3203  50.5,
3204  "TimeMonitoring");
3205  histos->fill2DHist(L1GMT_BXN,
3206  n_clcts,
3207  "n_CLCTs_v_BX_L1CSCCand",
3208  "Number of CLCTs vs. BX of L1 CSC Cand",
3209  4001,
3210  -0.5,
3211  4000.5,
3212  51,
3213  -0.5,
3214  50.5,
3215  "TimeMonitoring");
3216  histos->fill2DHist(L1GMT_BXN,
3217  n_correlatedlcts,
3218  "n_CorrelatedLCTs_v_BX_L1CSCCand",
3219  "Number of CorrelatedLCTs vs. BX of L1 CSC Cand",
3220  4001,
3221  -0.5,
3222  4000.5,
3223  51,
3224  -0.5,
3225  50.5,
3226  "TimeMonitoring");
3227  histos->fill2DHist(L1GMT_BXN,
3228  nRecHits,
3229  "n_RecHits_v_BX_L1CSCCand",
3230  "Number of RecHits vs. BX of L1 CSC Cand",
3231  4001,
3232  -0.5,
3233  4000.5,
3234  101,
3235  -0.5,
3236  100.5,
3237  "TimeMonitoring");
3238  histos->fill2DHist(L1GMT_BXN,
3239  nSegments,
3240  "n_Segments_v_BX_L1CSCCand",
3241  "Number of Segments vs. BX of L1 CSC Cand",
3242  4001,
3243  -0.5,
3244  4000.5,
3245  51,
3246  -0.5,
3247  50.5,
3248  "TimeMonitoring");
3249  }
3250  if (has_CSCTrigger && has_beamHaloTrigger) {
3251  histos->fill1DHist(
3252  L1GMT_BXN, "BX_L1CSCCand_w_beamHalo", "BX of L1 CSC (w beamHalo bit)", 4001, -0.5, 4000.5, "TimeMonitoring");
3253  histos->fill2DHist(L1GMT_BXN,
3254  n_alcts,
3255  "n_ALCTs_v_BX_L1CSCCand_w_beamHalo",
3256  "Number of ALCTs vs. BX of L1 CSC Cand (w beamHalo bit)",
3257  4001,
3258  -0.5,
3259  4000.5,
3260  51,
3261  -0.5,
3262  50.5,
3263  "TimeMonitoring");
3264  histos->fill2DHist(L1GMT_BXN,
3265  n_clcts,
3266  "n_CLCTs_v_BX_L1CSCCand_w_beamHalo",
3267  "Number of CLCTs vs. BX of L1 CSC Cand (w beamHalo bit)",
3268  4001,
3269  -0.5,
3270  4000.5,
3271  51,
3272  -0.5,
3273  50.5,
3274  "TimeMonitoring");
3275  histos->fill2DHist(L1GMT_BXN,
3276  n_correlatedlcts,
3277  "n_CorrelatedLCTs_v_BX_L1CSCCand_w_beamHalo",
3278  "Number of CorrelatedLCTs vs. BX of L1 CSC Cand (w beamHalo bit)",
3279  4001,
3280  -0.5,
3281  4000.5,
3282  51,
3283  -0.5,
3284  50.5,
3285  "TimeMonitoring");
3286  histos->fill2DHist(L1GMT_BXN,
3287  nRecHits,
3288  "n_RecHits_v_BX_L1CSCCand_w_beamHalo",
3289  "Number of RecHits vs. BX of L1 CSC Cand (w beamHalo bit)",
3290  4001,
3291  -0.5,
3292  4000.5,
3293  101,
3294  -0.5,
3295  100.5,
3296  "TimeMonitoring");
3297  histos->fill2DHist(L1GMT_BXN,
3298  nSegments,
3299  "n_Segments_v_BX_L1CSCCand_w_beamHalo",
3300  "Number of Segments vs. BX of L1 CSC Cand (w beamHalo bit)",
3301  4001,
3302  -0.5,
3303  4000.5,
3304  51,
3305  -0.5,
3306  50.5,
3307  "TimeMonitoring");
3308  }
3309 
3310  // *******************************************************************
3311  // Get information from the TMB header.
3312  // Can this eventually come out of the digis?
3313  // Taking code from EventFilter/CSCRawToDigis/CSCDCCUnpacker.cc
3314  // *******************************************************************
3315 
3316  edm::ESHandle<CSCCrateMap> hcrate;
3317  eventSetup.get<CSCCrateMapRcd>().get(hcrate);
3318  const CSCCrateMap* pcrate = hcrate.product();
3319 
3321  edm::Handle<FEDRawDataCollection> rawdata;
3322  event.getByToken(rd_token, rawdata);
3323  // If set selective unpacking mode
3324  // hardcoded examiner mask below to check for DCC and DDU level errors will be used first
3325  // then examinerMask for CSC level errors will be used during unpacking of each CSC block
3326  unsigned long dccBinCheckMask = 0x06080016;
3327  unsigned int examinerMask = 0x1FEBF3F6;
3328  unsigned int errorMask = 0x0;
3329 
3330  for (int id = FEDNumbering::MINCSCFEDID; id <= FEDNumbering::MAXCSCFEDID; ++id) {
3331  // loop over DCCs
3334 
3336  const FEDRawData& fedData = rawdata->FEDData(id);
3337  unsigned long length = fedData.size();
3338 
3339  if (length >= 32) {
3340  std::stringstream examiner_out, examiner_err;
3342  CSCDCCExaminer* examiner = new CSCDCCExaminer();
3343  if (examinerMask & 0x40000)
3344  examiner->crcCFEB(true);
3345  if (examinerMask & 0x8000)
3346  examiner->crcTMB(true);
3347  if (examinerMask & 0x0400)
3348  examiner->crcALCT(true);
3349  examiner->setMask(examinerMask);
3350  const short unsigned int* data = (short unsigned int*)fedData.data();
3351 
3352  bool goodEvent;
3353  if (examiner->check(data, long(fedData.size() / 2)) < 0) {
3354  goodEvent = false;
3355  } else {
3356  goodEvent = !(examiner->errors() & dccBinCheckMask);
3357  }
3358 
3359  if (goodEvent) {
3361  CSCDCCExaminer* ptrExaminer = examiner;
3362  CSCDCCEventData dccData((short unsigned int*)fedData.data(), ptrExaminer);
3363 
3365  const std::vector<CSCDDUEventData>& dduData = dccData.dduData();
3366 
3368  CSCDetId layer(1, 1, 1, 1, 1);
3369 
3370  for (unsigned int iDDU = 0; iDDU < dduData.size(); ++iDDU) { // loop over DDUs
3373  if (dduData[iDDU].trailer().errorstat() & errorMask) {
3374  LogTrace("CSCDCCUnpacker|CSCRawToDigi") << "DDU# " << iDDU << " has serious error - no digis unpacked! "
3375  << std::hex << dduData[iDDU].trailer().errorstat();
3376  continue; // to next iteration of DDU loop
3377  }
3378 
3380  const std::vector<CSCEventData>& cscData = dduData[iDDU].cscData();
3381  for (unsigned int iCSC = 0; iCSC < cscData.size(); ++iCSC) { // loop over CSCs
3382 
3383  int vmecrate = cscData[iCSC].dmbHeader()->crateID();
3384  int dmb = cscData[iCSC].dmbHeader()->dmbID();
3385 
3387  // SKIPPING MTCC redefinition of vmecrate
3388 
3389  int icfeb = 0;
3390  int ilayer = 0;
3391 
3392  if ((vmecrate >= 1) && (vmecrate <= 60) && (dmb >= 1) && (dmb <= 10) && (dmb != 6)) {
3393  layer = pcrate->detId(vmecrate, dmb, icfeb, ilayer);
3394  } else {
3395  LogTrace("CSCTimingAlignment|CSCDCCUnpacker|CSCRawToDigi") << " detID input out of range!!! ";
3396  LogTrace("CSCTimingAlignment|CSCDCCUnpacker|CSCRawToDigi")
3397  << " skipping chamber vme= " << vmecrate << " dmb= " << dmb;
3398  continue; // to next iteration of iCSC loop
3399  }
3400 
3402  int nalct = cscData[iCSC].dmbHeader()->nalct();
3403  bool goodALCT = false;
3404  //if (nalct&&(cscData[iCSC].dataPresent>>6&0x1)==1) {
3405  if (nalct && cscData[iCSC].alctHeader()) {
3406  if (cscData[iCSC].alctHeader()->check()) {
3407  goodALCT = true;
3408  }
3409  }
3410 
3412  int nclct = cscData[iCSC].dmbHeader()->nclct();
3413  bool goodTMB = false;
3414  if (nclct && cscData[iCSC].tmbData()) {
3415  if (cscData[iCSC].tmbHeader()->check()) {
3416  if (cscData[iCSC].comparatorData()->check())
3417  goodTMB = true;
3418  }
3419  }
3420 
3421  if (goodTMB && goodALCT) {
3422  if (ALCT_KeyWG_map.find(layer) == ALCT_KeyWG_map.end()) {
3423  printf("no ALCT info for Chamber %d %d %d %d \n",
3424  layer.endcap(),
3425  layer.station(),
3426  layer.ring(),
3427  layer.chamber());
3428  continue;
3429  }
3430  if (CLCT_getFullBx_map.find(layer) == CLCT_getFullBx_map.end()) {
3431  printf("no CLCT info for Chamber %d %d %d %d \n",
3432  layer.endcap(),
3433  layer.station(),
3434  layer.ring(),
3435  layer.chamber());
3436  continue;
3437  }
3438  int ALCT0Key = ALCT_KeyWG_map.find(layer)->second;
3439  int CLCTPretrigger = CLCT_getFullBx_map.find(layer)->second;
3440 
3441  const CSCTMBHeader* tmbHead = cscData[iCSC].tmbHeader();
3442 
3443  histos->fill1DHistByStation(tmbHead->BXNCount(),
3444  "TMB_BXNCount",
3445  "TMB_BXNCount",
3446  layer.chamberId(),
3447  3601,
3448  -0.5,
3449  3600.5,
3450  "TimeMonitoring");
3451  histos->fill1DHistByStation(tmbHead->ALCTMatchTime(),
3452  "TMB_ALCTMatchTime",
3453  "TMB_ALCTMatchTime",
3454  layer.chamberId(),
3455  7,
3456  -0.5,
3457  6.5,
3458  "TimeMonitoring");
3459 
3460  histos->fill1DHist(
3461  tmbHead->BXNCount(), "TMB_BXNCount", "TMB_BXNCount", 3601, -0.5, 3600.5, "TimeMonitoring");
3462  histos->fill1DHist(
3463  tmbHead->ALCTMatchTime(), "TMB_ALCTMatchTime", "TMB_ALCTMatchTime", 7, -0.5, 6.5, "TimeMonitoring");
3464 
3465  histos->fill1DHistByType(tmbHead->ALCTMatchTime(),
3466  "TMB_ALCTMatchTime",
3467  "TMB_ALCTMatchTime",
3468  layer.chamberId(),
3469  7,
3470  -0.5,
3471  6.5,
3472  "TimeMonitoring");
3473 
3474  histos->fillProfile(chamberSerial(layer.chamberId()),
3475  tmbHead->ALCTMatchTime(),
3476  "prof_TMB_ALCTMatchTime",
3477  "prof_TMB_ALCTMatchTime",
3478  601,
3479  -0.5,
3480  600.5,
3481  -0.5,
3482  7.5,
3483  "TimeMonitoring");
3484  histos->fillProfile(ALCT0Key,
3485  tmbHead->ALCTMatchTime(),
3486  "prof_TMB_ALCTMatchTime_v_ALCT0KeyWG",
3487  "prof_TMB_ALCTMatchTime_v_ALCT0KeyWG",
3488  128,
3489  -0.5,
3490  127.5,
3491  0,
3492  7,
3493  "TimeMonitoring");
3494  histos->fillProfileByType(ALCT0Key,
3495  tmbHead->ALCTMatchTime(),
3496  "prf_TMB_ALCTMatchTime_v_ALCT0KeyWG",
3497  "prf_TMB_ALCTMatchTime_v_ALCT0KeyWG",
3498  layer.chamberId(),
3499  128,
3500  -0.5,
3501  127.5,
3502  0,
3503  7,
3504  "TimeMonitoring");
3505 
3506  //Attempt to make a few sum plots
3507 
3508  int TMB_ALCT_rel_L1A = tmbHead->BXNCount() - (CLCTPretrigger + 2 + tmbHead->ALCTMatchTime());
3509  if (TMB_ALCT_rel_L1A > 3563)
3510  TMB_ALCT_rel_L1A = TMB_ALCT_rel_L1A - 3564;
3511  if (TMB_ALCT_rel_L1A < 0)
3512  TMB_ALCT_rel_L1A = TMB_ALCT_rel_L1A + 3564;
3513 
3514  //Plot TMB_ALCT_rel_L1A
3515  histos->fill1DHist(
3516  TMB_ALCT_rel_L1A, "h1D_TMB_ALCT_rel_L1A", "h1D_TMB_ALCT_rel_L1A", 11, 144.5, 155.5, "TimeMonitoring");
3517  histos->fill2DHist(chamberSerial(layer.chamberId()),
3518  TMB_ALCT_rel_L1A,
3519  "h2D_TMB_ALCT_rel_L1A",
3520  "h2D_TMB_ALCT_rel_L1A",
3521  601,
3522  -0.5,
3523  600.5,
3524  11,
3525  144.5,
3526  155.5,
3527  "TimeMonitoring");
3528  histos->fill2DHist(ringSerial(layer.chamberId()),
3529  TMB_ALCT_rel_L1A,
3530  "h2D_TMB_ALCT_rel_L1A_by_ring",
3531  "h2D_TMB_ALCT_rel_L1A_by_ring",
3532  19,
3533  -9.5,
3534  9.5,
3535  11,
3536  144.5,
3537  155.5,
3538  "TimeMonitoring");
3539  histos->fillProfile(chamberSerial(layer.chamberId()),
3540  TMB_ALCT_rel_L1A,
3541  "prof_TMB_ALCT_rel_L1A",
3542  "prof_TMB_ALCT_rel_L1A",
3543  601,
3544  -0.5,
3545  600.5,
3546  145,
3547  155,
3548  "TimeMonitoring");
3549  histos->fillProfile(ringSerial(layer.chamberId()),
3550  TMB_ALCT_rel_L1A,
3551  "prof_TMB_ALCT_rel_L1A_by_ring",
3552  "prof_TMB_ALCT_rel_L1A_by_ring",
3553  19,
3554  -9.5,
3555  9.5,
3556  145,
3557  155,
3558  "TimeMonitoring");
3559 
3560  histos->fill2DHist(ALCT0Key,
3561  TMB_ALCT_rel_L1A,
3562  "h2D_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3563  "h2D_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3564  128,
3565  -0.5,
3566  127.5,
3567  11,
3568  144.5,
3569  155.5,
3570  "TimeMonitoring");
3571  histos->fillProfile(ALCT0Key,
3572  TMB_ALCT_rel_L1A,
3573  "prof_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3574  "prof_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3575  128,
3576  -0.5,
3577  127.5,
3578  145,
3579  155,
3580  "TimeMonitoring");
3581  histos->fillProfileByType(ALCT0Key,
3582  TMB_ALCT_rel_L1A,
3583  "prf_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3584  "prf_TMB_ALCT_rel_L1A_v_ALCT0KeyWG",
3585  layer.chamberId(),
3586  128,
3587  -0.5,
3588  127.5,
3589  145,
3590  155,
3591  "TimeMonitoring");
3592  }
3593 
3594  } // end CSCData loop
3595  } // end ddu data loop
3596  } // end if goodEvent
3597  if (examiner != nullptr)
3598  delete examiner;
3599  } // end if non-zero fed data
3600  } // end DCC loop for NON-REFERENCE
3601 }
3602 
3603 void CSCValidation::endJob() { std::cout << "Events in " << nEventsAnalyzed << std::endl; }
3604 
3605 DEFINE_FWK_MODULE(CSCValidation);
CSCTMBHeader::BXNCount
uint16_t BXNCount() const
Definition: CSCTMBHeader.h:43
CSCDCCEventData
01/20/05 A.Tumanov
Definition: CSCDCCEventData.h:13
CSCValidation::getthisSignal
float getthisSignal(const CSCStripDigiCollection &stripdigis, CSCDetId idRH, int centerStrip)
Definition: CSCValidation.cc:2325
CSCDetId::chamber
int chamber() const
Definition: CSCDetId.h:62
CSCDCCExaminer::crcCFEB
void crcCFEB(bool enable)
Definition: CSCDCCExaminer.cc:45
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
LocalError::xx
float xx() const
Definition: LocalError.h:22
c
const edm::EventSetup & c
Definition: SiStripLAProfileBooker.cc:66
CSCValidation::fitX
float fitX(const CLHEP::HepMatrix &sp, const CLHEP::HepMatrix &ep)
Definition: CSCValidation.cc:1422
mps_fire.i
i
Definition: mps_fire.py:428
CSCSegment::localPosition
LocalPoint localPosition() const override
Definition: CSCSegment.h:39
CSCDetId
CSCDetId
Definition: RecoLocalMuon_RECO.doi:36
HLT_FULL_cff.pMin
tuple pMin
Definition: HLT_FULL_cff.py:13051
gpuClustering::id
uint16_t *__restrict__ id
Definition: gpuClusterChargeCut.h:20
TrapezoidalPlaneBounds::parameters
virtual const std::array< const float, 4 > parameters() const
Definition: TrapezoidalPlaneBounds.cc:49
CSCChamber::id
CSCDetId id() const
Get the (concrete) DetId.
Definition: CSCChamber.h:34
Reference_intrackfit_cff.endcap
endcap
Definition: Reference_intrackfit_cff.py:82
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
slope
static const double slope[3]
Definition: CastorTimeSlew.cc:6
CSCCrateMap::detId
CSCDetId detId(int vme, int dmb, int cfeb, int layer=0) const
Definition: CSCCrateMap.cc:9
GeomDet::toGlobal
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:49
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
TrackValidation_cff.pset
tuple pset
Definition: TrackValidation_cff.py:591
GlobalPoint
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
PixelTestBeamValidation_cfi.Efficiency
tuple Efficiency
Definition: PixelTestBeamValidation_cfi.py:87
CSCDCCEventData::dduData
const std::vector< CSCDDUEventData > & dduData() const
accessor to dduData
Definition: CSCDCCEventData.h:25
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:150
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
GeomDet::toLocal
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:58
CSCValidation::doNoiseHits
void doNoiseHits(edm::Handle< CSCRecHit2DCollection > recHits, edm::Handle< CSCSegmentCollection > cscSegments, edm::ESHandle< CSCGeometry > cscGeom, edm::Handle< CSCStripDigiCollection > strips)
Definition: CSCValidation.cc:1965
CSCValidation::endJob
void endJob() override
Definition: CSCValidation.cc:3603
CSCDBGains::gains
GainContainer gains
Definition: CSCDBGains.h:24
CSCDCCExaminer::errors
ExaminerStatusType errors(void) const
Definition: CSCDCCExaminer.h:170
CSCDCCExaminer::crcALCT
void crcALCT(bool enable)
Definition: CSCDCCExaminer.cc:29
CSCValidation::doRecHits
void doRecHits(edm::Handle< CSCRecHit2DCollection > recHits, edm::ESHandle< CSCGeometry > cscGeom)
Definition: CSCValidation.cc:896
CSCDetId::layer
int layer() const
Definition: CSCDetId.h:56
srCondWrite_cfg.deltaPhi
tuple deltaPhi
Definition: srCondWrite_cfg.py:54
CSCValidation::doEfficiencies
void doEfficiencies(edm::Handle< CSCWireDigiCollection > wires, edm::Handle< CSCStripDigiCollection > strips, edm::Handle< CSCRecHit2DCollection > recHits, edm::Handle< CSCSegmentCollection > cscSegments, edm::ESHandle< CSCGeometry > cscGeom)
Definition: CSCValidation.cc:1462
me0TriggerPseudoDigis_cff.nStrips
tuple nStrips
1.2 is to make the matching window safely the two nearest strips 0.35 is the size of an ME0 chamber i...
Definition: me0TriggerPseudoDigis_cff.py:26
FEDRawData::size
size_t size() const
Lenght of the data buffer in bytes.
Definition: FEDRawData.h:45
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
CSCValidation::withinSensitiveRegion
bool withinSensitiveRegion(LocalPoint localPos, const std::array< const float, 4 > &layerBounds, int station, int ring, float shiftFromEdge, float shiftFromDeadZone)
Definition: CSCValidation.cc:1797
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:234
CSCValidation::doTrigger
bool doTrigger(edm::Handle< L1MuGMTReadoutCollection > pCollection)
Definition: CSCValidation.cc:555
phase1PixelTopology::layer
constexpr std::array< uint8_t, layerIndexSize > layer
Definition: phase1PixelTopology.h:110
sistrip::SpyUtilities::range
const uint16_t range(const Frame &aFrame)
Definition: SiStripSpyUtilities.cc:76
PV3DBase::theta
Geom::Theta< T > theta() const
Definition: PV3DBase.h:72
CSCChamberSpecs::chamberTypeName
std::string chamberTypeName() const
Definition: CSCChamberSpecs.cc:156
CSCDetId::endcap
int endcap() const
Definition: CSCDetId.h:85
ztail.d
tuple d
Definition: ztail.py:151
CSCSegment::localDirection
LocalVector localDirection() const override
Local direction.
Definition: CSCSegment.h:42
edm::RangeMap::const_iterator
C::const_iterator const_iterator
constant access iterator type
Definition: RangeMap.h:43
CSCValidation::doCompTiming
void doCompTiming(const CSCComparatorDigiCollection &)
Definition: CSCValidation.cc:2782
dtT0Analyzer_cfg.rootFileName
tuple rootFileName
Definition: dtT0Analyzer_cfg.py:26
LocalError::xy
float xy() const
Definition: LocalError.h:23
CSCValidation::doPedestalNoise
void doPedestalNoise(edm::Handle< CSCStripDigiCollection > strips)
Definition: CSCValidation.cc:853
HLT_FULL_cff.chi2
tuple chi2
Definition: HLT_FULL_cff.py:10024
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:64
np
int np
Definition: AMPTWrapper.h:43
CSCIndexer::dbIndex
int dbIndex(const CSCDetId &id, int &channel)
Definition: CSCIndexer.cc:231
CSCChamber::specs
const CSCChamberSpecs * specs() const
Definition: CSCChamber.h:39
LocalError::yy
float yy() const
Definition: LocalError.h:24
CSCValidation::doStripDigis
void doStripDigis(edm::Handle< CSCStripDigiCollection > strips)
Definition: CSCValidation.cc:811
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
CSCValidation::getWidth
int getWidth(const CSCStripDigiCollection &stripdigis, CSCDetId idRH, int centerStrip)
Definition: CSCValidation.cc:2377
CSCValidation::histoEfficiency
void histoEfficiency(TH1F *readHisto, TH1F *writeHisto)
Definition: CSCValidation.cc:1778
DigiDM_cff.strips
tuple strips
#turn off noise in all subdetectors simHcalUnsuppressedDigis.doNoise = False mix.digitizers.hcal.doNoise = False simEcalUnsuppressedDigis.doNoise = False mix.digitizers.ecal.doNoise = False simEcalUnsuppressedDigis.doESNoise = False simSiPixelDigis.AddNoise = False mix.digitizers.pixel.AddNoise = False simSiStripDigis.Noise = False mix.digitizers.strip.AddNoise = False
Definition: DigiDM_cff.py:32
CSCValidation::ringSerial
int ringSerial(CSCDetId id)
Definition: CSCValidation.cc:1393
gpuVertexFinder::printf
printf("params %d %f %f %f\n", minT, eps, errmax, chi2max)
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
create_complex_test_file0_cfg.firstEvent
tuple firstEvent
Definition: create_complex_test_file0_cfg.py:18
if
if(conf_.getParameter< bool >("UseStripCablingDB"))
Definition: SiStripLAProfileBooker.cc:72
change_name.diff
list diff
Definition: change_name.py:13
CSCDCCExaminer::crcTMB
void crcTMB(bool enable)
Definition: CSCDCCExaminer.cc:37
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
CSCValidation::~CSCValidation
~CSCValidation() override
Destructor.
Definition: CSCValidation.cc:132
CSCDetId::chamberId
CSCDetId chamberId() const
Definition: CSCDetId.h:47
CSCChamber::layer
const CSCLayer * layer(CSCDetId id) const
Return the layer corresponding to the given id.
Definition: CSCChamber.cc:30
CSCValidation::filterEvents
bool filterEvents(edm::Handle< CSCRecHit2DCollection > recHits, edm::Handle< CSCSegmentCollection > cscSegments, edm::Handle< reco::TrackCollection > saMuons)
Definition: CSCValidation.cc:337
CSCValidation::doHLT
bool doHLT(edm::Handle< edm::TriggerResults > hltResults)
Definition: CSCValidation.cc:676
dtDQMClient_cfg.threshold
tuple threshold
Definition: dtDQMClient_cfg.py:15
ChiSquaredProbability
float ChiSquaredProbability(double chiSquared, double nrDOF)
Definition: ChiSquaredProbability.cc:13
kLayer
static const std::string kLayer("layer")
CSCValidation::analyze
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup) override
Perform the analysis.
Definition: CSCValidation.cc:176
CSCValidation::doSimHits
void doSimHits(edm::Handle< CSCRecHit2DCollection > recHits, edm::Handle< edm::PSimHitContainer > simHits)
Definition: CSCValidation.cc:1029
HLT_FULL_cff.muon
tuple muon
Definition: HLT_FULL_cff.py:11954
CSCTMBHeader::ALCTMatchTime
uint16_t ALCTMatchTime() const
Definition: CSCTMBHeader.h:44
CSCValidation::doGasGain
void doGasGain(const CSCWireDigiCollection &, const CSCStripDigiCollection &, const CSCRecHit2DCollection &)
Definition: CSCValidation.cc:2442
for
for(Iditer=Id.begin();Iditer!=Id.end();Iditer++)
Definition: SiStripLAProfileBooker.cc:79
CSCDetId::ring
int ring() const
Definition: CSCDetId.h:68
DetId
Definition: DetId.h:17
CSCValidation::getEfficiency
void getEfficiency(float bin, float Norm, std::vector< float > &eff)
Definition: CSCValidation.cc:1764
CSCValidation::doResolution
void doResolution(edm::Handle< CSCSegmentCollection > cscSegments, edm::ESHandle< CSCGeometry > cscGeom)
Definition: CSCValidation.cc:1175
HLT_FULL_cff.recHits
tuple recHits
Definition: HLT_FULL_cff.py:85944
CSCValidation::findNonAssociatedRecHits
void findNonAssociatedRecHits(edm::ESHandle< CSCGeometry > cscGeom, edm::Handle< CSCStripDigiCollection > strips)
Definition: CSCValidation.cc:2026
MuonDigiCollection::Range
std::pair< const_iterator, const_iterator > Range
Definition: MuonDigiCollection.h:95
MuonSubdetId::RPC
static constexpr int RPC
Definition: MuonSubdetId.h:13
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
digitizers_cfi.strip
strip
Definition: digitizers_cfi.py:45
MuonDigiCollection::const_iterator
std::vector< DigiType >::const_iterator const_iterator
Definition: MuonDigiCollection.h:94
trackerHits.simHits
tuple simHits
Definition: trackerHits.py:16
CSCDCCExaminer::check
int32_t check(const uint16_t *&buffer, int32_t length)
Definition: CSCDCCExaminer.cc:264
CSCValidation::doAFEBTiming
void doAFEBTiming(const CSCWireDigiCollection &)
Definition: CSCValidation.cc:2715
CSCValidation::doTimeMonitoring
void doTimeMonitoring(edm::Handle< CSCRecHit2DCollection > recHits, edm::Handle< CSCSegmentCollection > cscSegments, edm::Handle< CSCALCTDigiCollection > alcts, edm::Handle< CSCCLCTDigiCollection > clcts, edm::Handle< CSCCorrelatedLCTDigiCollection > correlatedlcts, edm::Handle< L1MuGMTReadoutCollection > pCollection, edm::ESHandle< CSCGeometry > cscGeom, const edm::EventSetup &eventSetup, const edm::Event &event)
Definition: CSCValidation.cc:2933
CSCGeometry::ChamberContainer
std::vector< const CSCChamber * > ChamberContainer
Definition: CSCGeometry.h:30
CSCValidation::getSignal
float getSignal(const CSCStripDigiCollection &stripdigis, CSCDetId idRH, int centerStrip)
Definition: CSCValidation.cc:1902
CSCDBPedestals::pedestals
PedestalContainer pedestals
Definition: CSCDBPedestals.h:26
data
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:79
CSCValidation::doStandalone
void doStandalone(edm::Handle< reco::TrackCollection > saMuons)
Definition: CSCValidation.cc:1256
CSCValidation::chamberSerial
int chamberSerial(CSCDetId id)
Definition: CSCValidation.cc:1358
edm::EventSetup::get
T get() const
Definition: EventSetup.h:88
edm::OwnVector::const_iterator
Definition: OwnVector.h:41
FEDRawData::data
const unsigned char * data() const
Return a const pointer to the beginning of the data buffer.
Definition: FEDRawData.cc:24
CSCValidation::doWireDigis
void doWireDigis(edm::Handle< CSCWireDigiCollection > wires)
Definition: CSCValidation.cc:776
CSCDetId::station
int station() const
Definition: CSCDetId.h:79
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144
CSCValidation::doADCTiming
void doADCTiming(const CSCRecHit2DCollection &)
Definition: CSCValidation.cc:2843
DigiDM_cff.wires
tuple wires
Definition: DigiDM_cff.py:33
CSCDBCrosstalk::crosstalk
CrosstalkContainer crosstalk
Definition: CSCDBCrosstalk.h:28
CSCDBNoiseMatrix::matrix
NoiseMatrixContainer matrix
Definition: CSCDBNoiseMatrix.h:24
edm::EventSetup::getHandle
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:157
dqmdumpme.last
tuple last
Definition: dqmdumpme.py:56
MuonSubdetId::DT
static constexpr int DT
Definition: MuonSubdetId.h:11
histoStyle.weight
int weight
Definition: histoStyle.py:51
CSCValidation::CSCValidation
CSCValidation(const edm::ParameterSet &pset)
Constructor.
Definition: CSCValidation.cc:16
cscSegments_cfi.cscSegments
tuple cscSegments
Definition: cscSegments_cfi.py:9
CSCValidation::doSegments
void doSegments(edm::Handle< CSCSegmentCollection > cscSegments, edm::ESHandle< CSCGeometry > cscGeom)
Definition: CSCValidation.cc:1081
MuonSubdetId::CSC
static constexpr int CSC
Definition: MuonSubdetId.h:12
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
CSCValidation::doOccupancies
void doOccupancies(edm::Handle< CSCStripDigiCollection > strips, edm::Handle< CSCWireDigiCollection > wires, edm::Handle< CSCRecHit2DCollection > recHits, edm::Handle< CSCSegmentCollection > cscSegments)
Definition: CSCValidation.cc:418
CSCLayer::geometry
const CSCLayerGeometry * geometry() const
Definition: CSCLayer.h:44
CSCDCCExaminer::setMask
void setMask(ExaminerMaskType mask)
Definition: CSCDCCExaminer.h:167
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
CSCValidation::doCalibrations
void doCalibrations(const edm::EventSetup &eventSetup)
Definition: CSCValidation.cc:693
Generated for CMSSW Reference Manual by doxygen 1.8.5