File indexing completed on 2021-09-07 22:41:50
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043 #include "DPGAnalysis/Skims/interface/CSCSkim.h"
0044
0045 using namespace std;
0046 using namespace edm;
0047
0048
0049
0050
0051 CSCSkim::CSCSkim(const edm::ParameterSet& pset) : m_CSCGeomToken(esConsumes()) {
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062 wds_token = consumes<CSCWireDigiCollection>(edm::InputTag("simMuonCSCDigis", "MuonCSCWireDigi"));
0063 sds_token = consumes<CSCStripDigiCollection>(edm::InputTag("simMuonCSCDigis", "MuonCSCStripDigi"));
0064 wdr_token = consumes<CSCWireDigiCollection>(edm::InputTag("muonCSCDigis", "MuonCSCWireDigi"));
0065 sdr_token = consumes<CSCStripDigiCollection>(edm::InputTag("muonCSCDigis", "MuonCSCStripDigi"));
0066
0067 rh_token = consumes<CSCRecHit2DCollection>(pset.getParameter<InputTag>("cscRecHitTag"));
0068 seg_token = consumes<CSCSegmentCollection>(pset.getParameter<InputTag>("cscSegmentTag"));
0069
0070 sam_token = consumes<reco::TrackCollection>(pset.getParameter<InputTag>("SAMuonTag"));
0071 trk_token = consumes<reco::TrackCollection>(pset.getParameter<InputTag>("trackTag"));
0072 glm_token = consumes<reco::MuonCollection>(pset.getParameter<InputTag>("GLBMuonTag"));
0073
0074
0075 outputFileName = pset.getUntrackedParameter<std::string>("outputFileName", "outputSkim.root");
0076 histogramFileName = pset.getUntrackedParameter<std::string>("histogramFileName", "histos.root");
0077 typeOfSkim = pset.getUntrackedParameter<int>("typeOfSkim", 1);
0078 nLayersWithHitsMinimum = pset.getUntrackedParameter<int>("nLayersWithHitsMinimum", 3);
0079 minimumHitChambers = pset.getUntrackedParameter<int>("minimumHitChambers", 1);
0080 minimumSegments = pset.getUntrackedParameter<int>("minimumSegments", 3);
0081 demandChambersBothSides = pset.getUntrackedParameter<bool>("demandChambersBothSides", false);
0082 makeHistograms = pset.getUntrackedParameter<bool>("makeHistograms", false);
0083 makeHistogramsForMessyEvents = pset.getUntrackedParameter<bool>("makeHistogramsForMessyEvebts", false);
0084 whichEndcap = pset.getUntrackedParameter<int>("whichEndcap", 2);
0085 whichStation = pset.getUntrackedParameter<int>("whichStation", 3);
0086 whichRing = pset.getUntrackedParameter<int>("whichRing", 2);
0087 whichChamber = pset.getUntrackedParameter<int>("whichChamber", 24);
0088
0089
0090 pMin = pset.getUntrackedParameter<double>("pMin", 3.);
0091 zLengthMin = pset.getUntrackedParameter<double>("zLengthMin", 200.);
0092 nCSCHitsMin = pset.getUntrackedParameter<int>("nCSCHitsMin", 9);
0093 zInnerMax = pset.getUntrackedParameter<double>("zInnerMax", 9000.);
0094 nTrHitsMin = pset.getUntrackedParameter<int>("nTrHitsMin", 8);
0095 zLengthTrMin = pset.getUntrackedParameter<double>("zLengthTrMin", 180.);
0096 rExtMax = pset.getUntrackedParameter<double>("rExtMax", 3000.);
0097 redChiSqMax = pset.getUntrackedParameter<double>("redChiSqMax", 20.);
0098 nValidHitsMin = pset.getUntrackedParameter<int>("nValidHitsMin", 8);
0099
0100 LogInfo("[CSCSkim] Setup") << "\n\t===== CSCSkim =====\n"
0101 << "\t\ttype of skim ...............................\t" << typeOfSkim
0102 << "\t\tminimum number of layers with hits .........\t" << nLayersWithHitsMinimum
0103 << "\n\t\tminimum number of chambers w/ hit layers..\t" << minimumHitChambers
0104 << "\n\t\tminimum number of segments ...............\t" << minimumSegments
0105 << "\n\t\tdemand chambers on both sides.............\t" << demandChambersBothSides
0106 << "\n\t\tmake histograms...........................\t" << makeHistograms
0107 << "\n\t\t..for messy events........................\t" << makeHistogramsForMessyEvents
0108 << "\n\t===================\n\n";
0109 }
0110
0111
0112
0113
0114 CSCSkim::~CSCSkim() {}
0115
0116
0117
0118
0119 void CSCSkim::beginJob() {
0120
0121 nEventsAnalyzed = 0;
0122 nEventsSelected = 0;
0123 nEventsChambersBothSides = 0;
0124 nEventsOverlappingChambers = 0;
0125 nEventsMessy = 0;
0126 nEventsCertainChamber = 0;
0127 nEventsDTOverlap = 0;
0128 nEventsHaloLike = 0;
0129 nEventsLongSATrack = 0;
0130 nEventsForBFieldStudies = 0;
0131 iRun = 0;
0132 iEvent = 0;
0133
0134 if (makeHistograms || makeHistogramsForMessyEvents) {
0135
0136 theHistogramFile = new TFile(histogramFileName.c_str(), "RECREATE");
0137 theHistogramFile->cd();
0138
0139 if (makeHistograms) {
0140
0141 hxnRecHits = new TH1F("hxnRecHits", "n RecHits", 61, -0.5, 60.5);
0142 hxnSegments = new TH1F("hxnSegments", "n Segments", 11, -0.5, 10.5);
0143 hxnHitChambers = new TH1F("hxnHitsChambers", "n chambers with hits", 11, -0.5, 10.5);
0144 hxnRecHitsSel = new TH1F("hxnRecHitsSel", "n RecHits selected", 61, -0.5, 60.5);
0145
0146 xxP = new TH1F("xxP", "P global", 100, 0., 200.);
0147 xxnValidHits = new TH1F("xxnValidHits", "n valid hits global", 61, -0.5, 60.5);
0148 xxnTrackerHits = new TH1F("xxnTrackerHits", "n tracker hits global", 61, -0.5, 60.5);
0149 xxnCSCHits = new TH1F("xxnCSCHits", "n CSC hits global", 41, -0.5, 40.5);
0150 xxredChiSq = new TH1F("xxredChiSq", "red chisq global", 100, 0., 100.);
0151 }
0152 if (makeHistogramsForMessyEvents) {
0153
0154 mevnRecHits0 = new TH1F("mevnRecHits0", "n RecHits", 121, -0.5, 120.5);
0155 mevnChambers0 = new TH1F("mevnChambers0", "n chambers with hits", 21, -0.5, 20.5);
0156 mevnSegments0 = new TH1F("mevnSegments0", "n Segments", 21, -0.5, 20.5);
0157 mevnRecHits1 = new TH1F("mevnRecHits1", "n RecHits", 100, 0., 300.);
0158 mevnChambers1 = new TH1F("mevnChambers1", "n chambers with hits", 50, 0., 50.);
0159 mevnSegments1 = new TH1F("mevnSegments1", "n Segments", 30, 0., 30.);
0160 }
0161 }
0162 }
0163
0164
0165
0166
0167 void CSCSkim::endJob() {
0168
0169
0170 float fraction = 0.;
0171 if (nEventsAnalyzed > 0) {
0172 fraction = (float)nEventsSelected / (float)nEventsAnalyzed;
0173 }
0174
0175 LogInfo("[CSCSkim] Summary") << "\n\n\t====== CSCSkim ==========================================================\n"
0176 << "\t\ttype of skim ...............................\t" << typeOfSkim << "\n"
0177 << "\t\tevents analyzed ..............\t" << nEventsAnalyzed << "\n"
0178 << "\t\tevents selected ..............\t" << nEventsSelected
0179 << "\tfraction= " << fraction << std::endl
0180 << "\t\tevents chambers both sides ...\t" << nEventsChambersBothSides << "\n"
0181 << "\t\tevents w/ overlaps .......... \t" << nEventsOverlappingChambers << "\n"
0182 << "\t\tevents lots of hit chambers . \t" << nEventsMessy << "\n"
0183 << "\t\tevents from certain chamber . \t" << nEventsCertainChamber << "\n"
0184 << "\t\tevents in DT-CSC overlap .... \t" << nEventsDTOverlap << "\n"
0185 << "\t\tevents halo-like ............ \t" << nEventsHaloLike << "\n"
0186 << "\t\tevents w/ long SA track ..... \t" << nEventsLongSATrack << "\n"
0187 << "\t\tevents good for BField ..... \t" << nEventsForBFieldStudies << "\n"
0188 << "\t=========================================================================\n\n";
0189
0190 if (makeHistograms || makeHistogramsForMessyEvents) {
0191
0192 LogDebug("[CSCSkim]") << "======= write out my histograms ====\n";
0193 theHistogramFile->cd();
0194 if (makeHistograms) {
0195 hxnRecHits->Write();
0196 hxnSegments->Write();
0197 hxnHitChambers->Write();
0198 hxnRecHitsSel->Write();
0199 }
0200 if (makeHistogramsForMessyEvents) {
0201 mevnRecHits0->Write();
0202 mevnChambers0->Write();
0203 mevnSegments0->Write();
0204 mevnRecHits1->Write();
0205 mevnChambers1->Write();
0206 mevnSegments1->Write();
0207 }
0208 theHistogramFile->Close();
0209 }
0210 }
0211
0212
0213
0214
0215 bool CSCSkim::filter(edm::Event& event, const edm::EventSetup& eventSetup) {
0216
0217 nEventsAnalyzed++;
0218
0219 iRun = event.id().run();
0220 iEvent = event.id().event();
0221
0222 LogDebug("[CSCSkim] EventInfo") << "Run: " << iRun << "\tEvent: " << iEvent << "\tn Analyzed: " << nEventsAnalyzed;
0223
0224
0225 ESHandle<CSCGeometry> cscGeom = eventSetup.getHandle(m_CSCGeomToken);
0226
0227
0228 edm::Handle<CSCWireDigiCollection> wires;
0229 edm::Handle<CSCStripDigiCollection> strips;
0230
0231 if (event.eventAuxiliary().isRealData()) {
0232 event.getByToken(wdr_token, wires);
0233 event.getByToken(sdr_token, strips);
0234 } else {
0235 event.getByToken(wds_token, wires);
0236 event.getByToken(sds_token, strips);
0237 }
0238
0239
0240 Handle<CSCRecHit2DCollection> cscRecHits;
0241 event.getByToken(rh_token, cscRecHits);
0242
0243
0244 Handle<CSCSegmentCollection> cscSegments;
0245 event.getByToken(seg_token, cscSegments);
0246
0247
0248 Handle<reco::TrackCollection> saMuons;
0249 if (typeOfSkim == 8) {
0250 event.getByToken(sam_token, saMuons);
0251 }
0252
0253
0254 Handle<reco::TrackCollection> tracks;
0255 Handle<reco::MuonCollection> gMuons;
0256 if (typeOfSkim == 9) {
0257 event.getByToken(sam_token, saMuons);
0258 event.getByToken(trk_token, tracks);
0259 event.getByToken(glm_token, gMuons);
0260 }
0261
0262
0263
0264
0265
0266
0267 bool basicEvent = false;
0268 if (typeOfSkim == 1 || typeOfSkim == 2) {
0269 basicEvent = doCSCSkimming(cscRecHits, cscSegments);
0270 }
0271
0272
0273 bool goodOverlapEvent = false;
0274 if (typeOfSkim == 3) {
0275 goodOverlapEvent = doOverlapSkimming(cscSegments);
0276 if (goodOverlapEvent) {
0277 nEventsOverlappingChambers++;
0278 }
0279 }
0280
0281
0282 bool messyEvent = false;
0283 if (typeOfSkim == 4) {
0284 messyEvent = doMessyEventSkimming(cscRecHits, cscSegments);
0285 if (messyEvent) {
0286 nEventsMessy++;
0287 }
0288 }
0289
0290
0291 bool hasChamber = false;
0292 if (typeOfSkim == 5) {
0293 hasChamber = doCertainChamberSelection(wires, strips);
0294 if (hasChamber) {
0295 nEventsCertainChamber++;
0296 }
0297 }
0298
0299
0300 bool DTOverlapCandidate = false;
0301 if (typeOfSkim == 6) {
0302 DTOverlapCandidate = doDTOverlap(cscSegments);
0303 if (DTOverlapCandidate) {
0304 nEventsDTOverlap++;
0305 }
0306 }
0307
0308
0309 bool HaloLike = false;
0310 if (typeOfSkim == 7) {
0311 HaloLike = doHaloLike(cscSegments);
0312 if (HaloLike) {
0313 nEventsHaloLike++;
0314 }
0315 }
0316
0317
0318 bool LongSATrack = false;
0319 if (typeOfSkim == 8) {
0320 LongSATrack = doLongSATrack(saMuons);
0321 if (LongSATrack) {
0322 nEventsLongSATrack++;
0323 }
0324 }
0325
0326
0327 bool GoodForBFieldStudy = false;
0328 if (typeOfSkim == 9) {
0329 GoodForBFieldStudy = doBFieldStudySelection(saMuons, tracks, gMuons);
0330 if (GoodForBFieldStudy) {
0331 nEventsForBFieldStudies++;
0332 }
0333 }
0334
0335
0336 bool selectThisEvent = false;
0337 if (typeOfSkim == 1 || typeOfSkim == 2) {
0338 selectThisEvent = basicEvent;
0339 }
0340 if (typeOfSkim == 3) {
0341 selectThisEvent = goodOverlapEvent;
0342 }
0343 if (typeOfSkim == 4) {
0344 selectThisEvent = messyEvent;
0345 }
0346 if (typeOfSkim == 5) {
0347 selectThisEvent = hasChamber;
0348 }
0349 if (typeOfSkim == 6) {
0350 selectThisEvent = DTOverlapCandidate;
0351 }
0352 if (typeOfSkim == 7) {
0353 selectThisEvent = HaloLike;
0354 }
0355 if (typeOfSkim == 8) {
0356 selectThisEvent = LongSATrack;
0357 }
0358 if (typeOfSkim == 9) {
0359 selectThisEvent = GoodForBFieldStudy;
0360 }
0361
0362 if (selectThisEvent) {
0363 nEventsSelected++;
0364 }
0365
0366 return selectThisEvent;
0367 }
0368
0369
0370
0371
0372
0373
0374
0375 bool CSCSkim::doCSCSkimming(edm::Handle<CSCRecHit2DCollection> cscRecHits,
0376 edm::Handle<CSCSegmentCollection> cscSegments) {
0377
0378 int nRecHits = cscRecHits->size();
0379
0380
0381 int cntRecHit[600];
0382 for (int i = 0; i < 600; i++) {
0383 cntRecHit[i] = 0;
0384 }
0385
0386
0387
0388
0389
0390 CSCRecHit2DCollection::const_iterator recIt;
0391 for (recIt = cscRecHits->begin(); recIt != cscRecHits->end(); recIt++) {
0392
0393 CSCDetId idrec = (CSCDetId)(*recIt).cscDetId();
0394 int kEndcap = idrec.endcap();
0395 int kRing = idrec.ring();
0396 int kStation = idrec.station();
0397 int kChamber = idrec.chamber();
0398 int kLayer = idrec.layer();
0399
0400
0401 int kSerial = chamberSerial(kEndcap, kStation, kRing, kChamber);
0402
0403
0404
0405 int kDigit = (int)pow((float)10., (float)(kLayer - 1));
0406 cntRecHit[kSerial] += kDigit;
0407
0408 }
0409
0410
0411
0412
0413
0414 int nChambersWithMinimalHits = 0;
0415 int nChambersWithMinimalHitsPOS = 0;
0416 int nChambersWithMinimalHitsNEG = 0;
0417 if (nRecHits > 0) {
0418 for (int i = 0; i < 600; i++) {
0419 if (cntRecHit[i] > 0) {
0420 int nLayersWithHits = 0;
0421 float dummy = (float)cntRecHit[i];
0422 for (int j = 5; j > -1; j--) {
0423 float digit = dummy / pow((float)10., (float)j);
0424 int kCount = (int)digit;
0425 if (kCount > 0)
0426 nLayersWithHits++;
0427 dummy = dummy - ((float)kCount) * pow((float)10., (float)j);
0428 }
0429 if (nLayersWithHits > nLayersWithHitsMinimum) {
0430 if (i < 300) {
0431 nChambersWithMinimalHitsPOS++;
0432 } else {
0433 nChambersWithMinimalHitsNEG++;
0434 }
0435 }
0436 }
0437 }
0438 nChambersWithMinimalHits = nChambersWithMinimalHitsPOS + nChambersWithMinimalHitsNEG;
0439 }
0440
0441
0442 int nSegments = cscSegments->size();
0443
0444
0445
0446
0447
0448 if (makeHistograms) {
0449 hxnRecHits->Fill(nRecHits);
0450 if (nRecHits > 0) {
0451 hxnSegments->Fill(nSegments);
0452 hxnHitChambers->Fill(nChambersWithMinimalHits);
0453 }
0454 if (nChambersWithMinimalHits > 0) {
0455 hxnRecHitsSel->Fill(nRecHits);
0456 }
0457 }
0458
0459
0460
0461
0462 bool basicEvent = (nChambersWithMinimalHits >= minimumHitChambers) && (nSegments >= minimumSegments);
0463
0464 bool chambersOnBothSides =
0465 ((nChambersWithMinimalHitsPOS >= minimumHitChambers) && (nChambersWithMinimalHitsNEG >= minimumHitChambers));
0466
0467 if (chambersOnBothSides) {
0468 nEventsChambersBothSides++;
0469 }
0470
0471 bool selectEvent = false;
0472 if (typeOfSkim == 1) {
0473 selectEvent = basicEvent;
0474 }
0475 if (typeOfSkim == 2) {
0476 selectEvent = chambersOnBothSides;
0477 }
0478
0479
0480 LogDebug("[CSCSkim]") << "----- nRecHits = " << nRecHits
0481 << "\tnChambersWithMinimalHits = " << nChambersWithMinimalHits << "\tnSegments = " << nSegments
0482 << "\tselect? " << selectEvent << std::endl;
0483
0484
0485
0486
0487
0488
0489
0490
0491
0492
0493
0494
0495
0496
0497
0498
0499
0500
0501
0502
0503 return selectEvent;
0504 }
0505
0506
0507
0508
0509
0510 bool CSCSkim::doOverlapSkimming(edm::Handle<CSCSegmentCollection> cscSegments) {
0511 const int nhitsMinimum = 4;
0512 const float chisqMaximum = 100.;
0513 const int nAllMaximum = 3;
0514
0515
0516
0517
0518
0519 int nAll[600];
0520 int nGood[600];
0521 for (int i = 0; i < 600; i++) {
0522 nAll[i] = 0;
0523 nGood[i] = 0;
0524 }
0525
0526
0527
0528
0529 for (CSCSegmentCollection::const_iterator it = cscSegments->begin(); it != cscSegments->end(); it++) {
0530
0531 CSCDetId id = (CSCDetId)(*it).cscDetId();
0532 int kEndcap = id.endcap();
0533 int kStation = id.station();
0534 int kRing = id.ring();
0535 int kChamber = id.chamber();
0536 int kSerial = chamberSerial(kEndcap, kStation, kRing, kChamber);
0537
0538
0539 float chisq = (*it).chi2();
0540 int nhits = (*it).nRecHits();
0541
0542
0543 bool goodSegment = (nhits >= nhitsMinimum) && (chisq < chisqMaximum);
0544
0545
0546
0547
0548
0549
0550
0551
0552
0553
0554
0555
0556 nAll[kSerial - 1]++;
0557 if (goodSegment)
0558 nGood[kSerial]++;
0559
0560 }
0561
0562
0563
0564
0565
0566
0567 bool messyChamber = false;
0568 for (int i = 0; i < 600; i++) {
0569 if (nAll[i] > nAllMaximum)
0570 messyChamber = true;
0571 }
0572
0573
0574
0575 bool consecutiveChambers = false;
0576 for (int i = 0; i < 599; i++) {
0577 if ((nGood[i] > 0) && (nGood[i + 1] > 0))
0578 consecutiveChambers = true;
0579 }
0580
0581 bool selectThisEvent = !messyChamber && consecutiveChambers;
0582
0583 return selectThisEvent;
0584 }
0585
0586
0587
0588
0589
0590
0591
0592 bool CSCSkim::doMessyEventSkimming(edm::Handle<CSCRecHit2DCollection> cscRecHits,
0593 edm::Handle<CSCSegmentCollection> cscSegments) {
0594
0595 int nRecHits = cscRecHits->size();
0596
0597
0598 int cntRecHit[600];
0599 for (int i = 0; i < 600; i++) {
0600 cntRecHit[i] = 0;
0601 }
0602
0603
0604
0605
0606
0607 CSCRecHit2DCollection::const_iterator recIt;
0608 for (recIt = cscRecHits->begin(); recIt != cscRecHits->end(); recIt++) {
0609
0610 CSCDetId idrec = (CSCDetId)(*recIt).cscDetId();
0611 int kEndcap = idrec.endcap();
0612 int kRing = idrec.ring();
0613 int kStation = idrec.station();
0614 int kChamber = idrec.chamber();
0615 int kLayer = idrec.layer();
0616
0617
0618 int kSerial = chamberSerial(kEndcap, kStation, kRing, kChamber);
0619
0620
0621
0622 int kDigit = (int)pow((float)10., (float)(kLayer - 1));
0623 cntRecHit[kSerial] += kDigit;
0624
0625 }
0626
0627
0628
0629
0630
0631 int nChambersWithMinimalHits = 0;
0632 int nChambersWithMinimalHitsPOS = 0;
0633 int nChambersWithMinimalHitsNEG = 0;
0634 if (nRecHits > 0) {
0635 for (int i = 0; i < 600; i++) {
0636 if (cntRecHit[i] > 0) {
0637 int nLayersWithHits = 0;
0638 float dummy = (float)cntRecHit[i];
0639 for (int j = 5; j > -1; j--) {
0640 float digit = dummy / pow((float)10., (float)j);
0641 int kCount = (int)digit;
0642 if (kCount > 0)
0643 nLayersWithHits++;
0644 dummy = dummy - ((float)kCount) * pow((float)10., (float)j);
0645 }
0646 if (nLayersWithHits > nLayersWithHitsMinimum) {
0647 if (i < 300) {
0648 nChambersWithMinimalHitsPOS++;
0649 } else {
0650 nChambersWithMinimalHitsNEG++;
0651 }
0652 }
0653 }
0654 }
0655 nChambersWithMinimalHits = nChambersWithMinimalHitsPOS + nChambersWithMinimalHitsNEG;
0656 }
0657
0658
0659 int nSegments = cscSegments->size();
0660
0661
0662
0663
0664
0665 if (makeHistogramsForMessyEvents) {
0666 if (nRecHits > 8) {
0667 mevnRecHits0->Fill(nRecHits);
0668 mevnChambers0->Fill(nChambersWithMinimalHits);
0669 mevnSegments0->Fill(nSegments);
0670 }
0671 if (nRecHits > 54) {
0672 double dummy = (double)nRecHits;
0673 if (dummy > 299.9)
0674 dummy = 299.9;
0675 mevnRecHits1->Fill(dummy);
0676 dummy = (double)nChambersWithMinimalHits;
0677 if (dummy > 49.9)
0678 dummy = 49.9;
0679 mevnChambers1->Fill(dummy);
0680 dummy = (double)nSegments;
0681 if (dummy > 29.9)
0682 dummy = 29.9;
0683 mevnSegments1->Fill(dummy);
0684 }
0685 }
0686
0687
0688
0689
0690
0691 bool selectEvent = false;
0692 if ((nRecHits > 54) && (nChambersWithMinimalHits > 5)) {
0693 selectEvent = true;
0694 }
0695
0696
0697 LogDebug("[CSCSkim]") << "----- nRecHits = " << nRecHits
0698 << "\tnChambersWithMinimalHits = " << nChambersWithMinimalHits << "\tnSegments = " << nSegments
0699 << "\tselect? " << selectEvent << std::endl;
0700
0701
0702
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715
0716
0717
0718
0719
0720 return selectEvent;
0721 }
0722
0723
0724
0725
0726
0727
0728 bool CSCSkim::doCertainChamberSelection(edm::Handle<CSCWireDigiCollection> wires,
0729 edm::Handle<CSCStripDigiCollection> strips) {
0730
0731 bool certainChamberIsPresentInWires = false;
0732 for (CSCWireDigiCollection::DigiRangeIterator jw = wires->begin(); jw != wires->end(); jw++) {
0733 CSCDetId id = (CSCDetId)(*jw).first;
0734 int kEndcap = id.endcap();
0735 int kRing = id.ring();
0736 int kStation = id.station();
0737 int kChamber = id.chamber();
0738 if ((kEndcap == whichEndcap) && (kStation == whichStation) && (kRing == whichRing) && (kChamber == whichChamber)) {
0739 certainChamberIsPresentInWires = true;
0740 }
0741 }
0742
0743
0744 bool certainChamberIsPresentInStrips = false;
0745 for (CSCStripDigiCollection::DigiRangeIterator js = strips->begin(); js != strips->end(); js++) {
0746 CSCDetId id = (CSCDetId)(*js).first;
0747 int kEndcap = id.endcap();
0748 int kRing = id.ring();
0749 int kStation = id.station();
0750 int kChamber = id.chamber();
0751 if ((kEndcap == whichEndcap) && (kStation == whichStation) && (kRing == whichRing) && (kChamber == whichChamber)) {
0752 certainChamberIsPresentInStrips = true;
0753 }
0754 }
0755
0756 bool certainChamberIsPresent = certainChamberIsPresentInWires || certainChamberIsPresentInStrips;
0757
0758 return certainChamberIsPresent;
0759 }
0760
0761
0762
0763
0764
0765
0766 bool CSCSkim::doDTOverlap(Handle<CSCSegmentCollection> cscSegments) {
0767 const float chisqMax = 100.;
0768 const int nhitsMin = 5;
0769 const int maxNSegments = 3;
0770
0771
0772 bool DTOverlapCandidate = false;
0773 int cntMEP13[36];
0774 int cntMEN13[36];
0775 int cntMEP22[36];
0776 int cntMEN22[36];
0777 int cntMEP32[36];
0778 int cntMEN32[36];
0779 for (int i = 0; i < 36; ++i) {
0780 cntMEP13[i] = 0;
0781 cntMEN13[i] = 0;
0782 cntMEP22[i] = 0;
0783 cntMEN22[i] = 0;
0784 cntMEP32[i] = 0;
0785 cntMEN32[i] = 0;
0786 }
0787
0788
0789
0790
0791
0792 int nSegments = cscSegments->size();
0793 if (nSegments < 2)
0794 return DTOverlapCandidate;
0795
0796 for (CSCSegmentCollection::const_iterator it = cscSegments->begin(); it != cscSegments->end(); it++) {
0797
0798 CSCDetId id = (CSCDetId)(*it).cscDetId();
0799 int kEndcap = id.endcap();
0800 int kStation = id.station();
0801 int kRing = id.ring();
0802 int kChamber = id.chamber();
0803
0804 float chisq = (*it).chi2();
0805 int nhits = (*it).nRecHits();
0806 bool goodSegment = (chisq < chisqMax) && (nhits >= nhitsMin);
0807 if (goodSegment) {
0808 if ((kStation == 1) && (kRing == 3)) {
0809 if (kEndcap == 1) {
0810 cntMEP13[kChamber - 1]++;
0811 }
0812 if (kEndcap == 2) {
0813 cntMEN13[kChamber - 1]++;
0814 }
0815 }
0816 if ((kStation == 2) && (kRing == 2)) {
0817 if (kEndcap == 1) {
0818 cntMEP22[kChamber - 1]++;
0819 }
0820 if (kEndcap == 2) {
0821 cntMEN22[kChamber - 1]++;
0822 }
0823 }
0824 if ((kStation == 3) && (kRing == 2)) {
0825 if (kEndcap == 1) {
0826 cntMEP32[kChamber - 1]++;
0827 }
0828 if (kEndcap == 2) {
0829 cntMEN32[kChamber - 1]++;
0830 }
0831 }
0832 }
0833 }
0834
0835
0836
0837
0838 bool tooManySegments = false;
0839 for (int i = 0; i < 36; ++i) {
0840 if ((cntMEP13[i] > maxNSegments) || (cntMEN13[i] > maxNSegments) || (cntMEP22[i] > maxNSegments) ||
0841 (cntMEN22[i] > maxNSegments) || (cntMEP32[i] > maxNSegments) || (cntMEN32[i] > maxNSegments))
0842 tooManySegments = true;
0843 }
0844 if (tooManySegments) {
0845 return DTOverlapCandidate;
0846 }
0847
0848
0849
0850
0851 bool matchup = false;
0852 for (int i = 0; i < 36; ++i) {
0853 if ((cntMEP13[i] > 0) && (cntMEP22[i] + cntMEP32[i] > 0)) {
0854 matchup = true;
0855 }
0856 if ((cntMEN13[i] > 0) && (cntMEN22[i] + cntMEN32[i] > 0)) {
0857 matchup = true;
0858 }
0859 }
0860
0861
0862
0863
0864
0865
0866
0867
0868
0869
0870
0871
0872
0873
0874
0875 DTOverlapCandidate = matchup;
0876 return DTOverlapCandidate;
0877 }
0878
0879
0880
0881
0882
0883
0884
0885 bool CSCSkim::doHaloLike(Handle<CSCSegmentCollection> cscSegments) {
0886 const float chisqMax = 100.;
0887 const int nhitsMin = 5;
0888 const int maxNSegments = 3;
0889
0890
0891 bool HaloLike = false;
0892 int cntMEP11[36];
0893 int cntMEN11[36];
0894 int cntMEP12[36];
0895 int cntMEN12[36];
0896 int cntMEP21[36];
0897 int cntMEN21[36];
0898 int cntMEP31[36];
0899 int cntMEN31[36];
0900 int cntMEP41[36];
0901 int cntMEN41[36];
0902 for (int i = 0; i < 36; ++i) {
0903 cntMEP11[i] = 0;
0904 cntMEN11[i] = 0;
0905 cntMEP12[i] = 0;
0906 cntMEN12[i] = 0;
0907 cntMEP21[i] = 0;
0908 cntMEN21[i] = 0;
0909 cntMEP31[i] = 0;
0910 cntMEN31[i] = 0;
0911 cntMEP41[i] = 0;
0912 cntMEN41[i] = 0;
0913 }
0914
0915
0916
0917
0918 int nSegments = cscSegments->size();
0919 if (nSegments < 4)
0920 return HaloLike;
0921
0922 for (CSCSegmentCollection::const_iterator it = cscSegments->begin(); it != cscSegments->end(); it++) {
0923
0924 CSCDetId id = (CSCDetId)(*it).cscDetId();
0925 int kEndcap = id.endcap();
0926 int kStation = id.station();
0927 int kRing = id.ring();
0928 int kChamber = id.chamber();
0929
0930 float chisq = (*it).chi2();
0931 int nhits = (*it).nRecHits();
0932 bool goodSegment = (chisq < chisqMax) && (nhits >= nhitsMin);
0933 if (goodSegment) {
0934 if ((kStation == 1) && (kRing == 1)) {
0935 if (kEndcap == 1) {
0936 cntMEP11[kChamber - 1]++;
0937 }
0938 if (kEndcap == 2) {
0939 cntMEN11[kChamber - 1]++;
0940 }
0941 }
0942 if ((kStation == 1) && (kRing == 2)) {
0943 if (kEndcap == 1) {
0944 cntMEP12[kChamber - 1]++;
0945 }
0946 if (kEndcap == 2) {
0947 cntMEN12[kChamber - 1]++;
0948 }
0949 }
0950 if ((kStation == 2) && (kRing == 1)) {
0951 if (kEndcap == 1) {
0952 cntMEP21[kChamber - 1]++;
0953 }
0954 if (kEndcap == 2) {
0955 cntMEN21[kChamber - 1]++;
0956 }
0957 }
0958 if ((kStation == 3) && (kRing == 1)) {
0959 if (kEndcap == 1) {
0960 cntMEP31[kChamber - 1]++;
0961 }
0962 if (kEndcap == 2) {
0963 cntMEN31[kChamber - 1]++;
0964 }
0965 }
0966 if ((kStation == 4) && (kRing == 1)) {
0967 if (kEndcap == 1) {
0968 cntMEP41[kChamber - 1]++;
0969 }
0970 if (kEndcap == 2) {
0971 cntMEN41[kChamber - 1]++;
0972 }
0973 }
0974 }
0975 }
0976
0977
0978
0979
0980 bool tooManySegments = false;
0981 for (int i = 0; i < 36; ++i) {
0982 if ((cntMEP11[i] > 3 * maxNSegments) || (cntMEN11[i] > 3 * maxNSegments) || (cntMEP12[i] > maxNSegments) ||
0983 (cntMEN12[i] > maxNSegments) || (cntMEP21[i] > maxNSegments) || (cntMEN21[i] > maxNSegments) ||
0984 (cntMEP31[i] > maxNSegments) || (cntMEN31[i] > maxNSegments) || (cntMEP41[i] > maxNSegments) ||
0985 (cntMEN41[i] > maxNSegments))
0986 tooManySegments = true;
0987 }
0988 if (tooManySegments) {
0989 return HaloLike;
0990 }
0991
0992
0993
0994
0995 bool matchup = false;
0996 for (int i = 0; i < 36; ++i) {
0997 if ((cntMEP11[i] + cntMEP12[i] > 0) && (cntMEP21[i] > 0) && (cntMEP31[i] > 0) && (cntMEP41[i] > 0)) {
0998 matchup = true;
0999 }
1000 if ((cntMEN11[i] + cntMEN12[i] > 0) && (cntMEN21[i] > 0) && (cntMEN31[i] > 0) && (cntMEN41[i] > 0)) {
1001 matchup = true;
1002 }
1003 }
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036 HaloLike = matchup;
1037 return HaloLike;
1038 }
1039
1040
1041
1042
1043 bool CSCSkim::doLongSATrack(edm::Handle<reco::TrackCollection> saMuons) {
1044 const float zDistanceMax = 2500.;
1045 const float zDistanceMin = 700.;
1046 const int nCSCHitsMin = 25;
1047 const int nCSCHitsMax = 50;
1048 const float zInnerMax = 80000.;
1049
1050 const int nNiceMuonsMin = 1;
1051
1052
1053
1054
1055
1056 int nNiceMuons = 0;
1057
1058 for (reco::TrackCollection::const_iterator muon = saMuons->begin(); muon != saMuons->end(); ++muon) {
1059
1060 math::XYZVector innerMo = muon->innerMomentum();
1061 GlobalVector im(innerMo.x(), innerMo.y(), innerMo.z());
1062 math::XYZPoint innerPo = muon->innerPosition();
1063 GlobalPoint ip(innerPo.x(), innerPo.y(), innerPo.z());
1064 math::XYZPoint outerPo = muon->outerPosition();
1065 GlobalPoint op(outerPo.x(), outerPo.y(), outerPo.z());
1066 float zInner = ip.z();
1067 float zOuter = op.z();
1068 float zDistance = fabs(zOuter - zInner);
1069
1070
1071 int nDTHits = 0;
1072 int nCSCHits = 0;
1073 for (trackingRecHit_iterator hit = muon->recHitsBegin(); hit != muon->recHitsEnd(); ++hit) {
1074 const DetId detId((*hit)->geographicalId());
1075 if (detId.det() == DetId::Muon) {
1076 if (detId.subdetId() == MuonSubdetId::DT) {
1077
1078
1079 nDTHits++;
1080 } else if (detId.subdetId() == MuonSubdetId::CSC) {
1081
1082
1083 nCSCHits++;
1084 }
1085 }
1086 }
1087
1088
1089 if ((zDistance < zDistanceMax) && (zDistance > zDistanceMin) && (nCSCHits > nCSCHitsMin) &&
1090 (nCSCHits < nCSCHitsMax) && (min(fabs(zInner), fabs(zOuter)) < zInnerMax) &&
1091 (fabs(innerMo.z()) > 0.000000001)) {
1092 nNiceMuons++;
1093 }
1094 }
1095
1096 bool select = (nNiceMuons >= nNiceMuonsMin);
1097
1098 return select;
1099 }
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109 bool CSCSkim::doBFieldStudySelection(edm::Handle<reco::TrackCollection> saMuons,
1110 edm::Handle<reco::TrackCollection> tracks,
1111 edm::Handle<reco::MuonCollection> gMuons) {
1112 bool acceptThisEvent = false;
1113
1114
1115
1116
1117 int nGoodSAMuons = 0;
1118 for (reco::TrackCollection::const_iterator muon = saMuons->begin(); muon != saMuons->end(); ++muon) {
1119 float preco = muon->p();
1120
1121 math::XYZPoint innerPo = muon->innerPosition();
1122 GlobalPoint iPnt(innerPo.x(), innerPo.y(), innerPo.z());
1123 math::XYZPoint outerPo = muon->outerPosition();
1124 GlobalPoint oPnt(outerPo.x(), outerPo.y(), outerPo.z());
1125 float zLength = abs(iPnt.z() - oPnt.z());
1126
1127 math::XYZVector innerMom = muon->innerMomentum();
1128 GlobalVector iP(innerMom.x(), innerMom.y(), innerMom.z());
1129 math::XYZVector outerMom = muon->outerMomentum();
1130 GlobalVector oP(outerMom.x(), outerMom.y(), outerMom.z());
1131
1132 const float zRef = 300.;
1133 float xExt = 10000.;
1134 float yExt = 10000.;
1135 if (abs(oPnt.z()) < abs(iPnt.z())) {
1136 float deltaZ = 0.;
1137 if (oPnt.z() > 0) {
1138 deltaZ = zRef - oPnt.z();
1139 } else {
1140 deltaZ = -zRef - oPnt.z();
1141 }
1142 xExt = oPnt.x() + deltaZ * oP.x() / oP.z();
1143 yExt = oPnt.y() + deltaZ * oP.y() / oP.z();
1144 } else {
1145 float deltaZ = 0.;
1146 if (iPnt.z() > 0) {
1147 deltaZ = zRef - iPnt.z();
1148 } else {
1149 deltaZ = -zRef - iPnt.z();
1150 }
1151 xExt = iPnt.x() + deltaZ * iP.x() / iP.z();
1152 yExt = iPnt.y() + deltaZ * iP.y() / iP.z();
1153 }
1154 float rExt = sqrt(xExt * xExt + yExt * yExt);
1155
1156 int kHit = 0;
1157 int nDTHits = 0;
1158 int nCSCHits = 0;
1159 for (trackingRecHit_iterator hit = muon->recHitsBegin(); hit != muon->recHitsEnd(); ++hit) {
1160 ++kHit;
1161 const DetId detId((*hit)->geographicalId());
1162 if (detId.det() == DetId::Muon) {
1163 if (detId.subdetId() == MuonSubdetId::DT) {
1164 nDTHits++;
1165 } else if (detId.subdetId() == MuonSubdetId::CSC) {
1166 nCSCHits++;
1167 }
1168 }
1169 }
1170
1171 float zInner = -1.;
1172 if (nCSCHits >= nCSCHitsMin) {
1173 if (abs(iPnt.z()) < abs(iPnt.z())) {
1174 zInner = iPnt.z();
1175 } else {
1176 zInner = oPnt.z();
1177 }
1178 }
1179
1180 bool goodSAMuon = (preco > pMin) && (zLength > zLengthMin) && (nCSCHits >= nCSCHitsMin) && (zInner < zInnerMax) &&
1181 (rExt < rExtMax);
1182
1183 if (goodSAMuon) {
1184 nGoodSAMuons++;
1185 }
1186
1187 }
1188
1189
1190
1191
1192 int nGoodTracks = 0;
1193 for (reco::TrackCollection::const_iterator track = tracks->begin(); track != tracks->end(); ++track) {
1194 float preco = track->p();
1195 int n = track->recHitsSize();
1196
1197 math::XYZPoint innerPo = track->innerPosition();
1198 GlobalPoint iPnt(innerPo.x(), innerPo.y(), innerPo.z());
1199 math::XYZPoint outerPo = track->outerPosition();
1200 GlobalPoint oPnt(outerPo.x(), outerPo.y(), outerPo.z());
1201 float zLength = abs(iPnt.z() - oPnt.z());
1202
1203 math::XYZVector innerMom = track->innerMomentum();
1204 GlobalVector iP(innerMom.x(), innerMom.y(), innerMom.z());
1205 math::XYZVector outerMom = track->outerMomentum();
1206 GlobalVector oP(outerMom.x(), outerMom.y(), outerMom.z());
1207
1208 const float zRef = 300.;
1209 float xExt = 10000.;
1210 float yExt = 10000.;
1211 if (abs(oPnt.z()) > abs(iPnt.z())) {
1212 float deltaZ = 0.;
1213 if (oPnt.z() > 0) {
1214 deltaZ = zRef - oPnt.z();
1215 } else {
1216 deltaZ = -zRef - oPnt.z();
1217 }
1218 xExt = oPnt.x() + deltaZ * oP.x() / oP.z();
1219 yExt = oPnt.y() + deltaZ * oP.y() / oP.z();
1220 } else {
1221 float deltaZ = 0.;
1222 if (iPnt.z() > 0) {
1223 deltaZ = zRef - iPnt.z();
1224 } else {
1225 deltaZ = -zRef - iPnt.z();
1226 }
1227 xExt = iPnt.x() + deltaZ * iP.x() / iP.z();
1228 yExt = iPnt.y() + deltaZ * iP.y() / iP.z();
1229 }
1230 float rExt = sqrt(xExt * xExt + yExt * yExt);
1231
1232 bool goodTrack = (preco > pMin) && (n >= nTrHitsMin) && (zLength > zLengthTrMin) && (rExt < rExtMax);
1233
1234 if (goodTrack) {
1235 nGoodTracks++;
1236 }
1237
1238 }
1239
1240
1241
1242
1243 int nGoodGlobalMuons = 0;
1244 for (reco::MuonCollection::const_iterator global = gMuons->begin(); global != gMuons->end(); ++global) {
1245 if (global->isGlobalMuon()) {
1246 float pDef = global->p();
1247 float redChiSq = global->globalTrack()->normalizedChi2();
1248 const reco::HitPattern& hp = (global->globalTrack())->hitPattern();
1249
1250
1251 int nTrackerHits = hp.numberOfValidTrackerHits();
1252
1253
1254
1255 int nDTHits = 0;
1256 int nCSCHits = 0;
1257 for (trackingRecHit_iterator hit = (global->globalTrack())->recHitsBegin();
1258 hit != (global->globalTrack())->recHitsEnd();
1259 ++hit) {
1260 const DetId detId((*hit)->geographicalId());
1261 if (detId.det() == DetId::Muon) {
1262 if (detId.subdetId() == MuonSubdetId::DT) {
1263 nDTHits++;
1264 } else if (detId.subdetId() == MuonSubdetId::CSC) {
1265 nCSCHits++;
1266 }
1267 }
1268 }
1269
1270 bool goodGlobalMuon =
1271 (pDef > pMin) && (nTrackerHits >= nValidHitsMin) && (nCSCHits >= nCSCHitsMin) && (redChiSq < redChiSqMax);
1272
1273 if (goodGlobalMuon) {
1274 nGoodGlobalMuons++;
1275 }
1276
1277 }
1278 }
1279
1280
1281
1282
1283
1284 acceptThisEvent = ((nGoodSAMuons > 0) && (nGoodTracks > 0)) || (nGoodGlobalMuons > 0);
1285
1286 return acceptThisEvent;
1287 }
1288
1289
1290
1291
1292
1293 int CSCSkim::chamberSerial(int kEndcap, int kStation, int kRing, int kChamber) {
1294 int kSerial = kChamber;
1295 if (kStation == 1 && kRing == 1) {
1296 kSerial = kChamber;
1297 }
1298 if (kStation == 1 && kRing == 2) {
1299 kSerial = kChamber + 36;
1300 }
1301 if (kStation == 1 && kRing == 3) {
1302 kSerial = kChamber + 72;
1303 }
1304 if (kStation == 1 && kRing == 4) {
1305 kSerial = kChamber;
1306 }
1307 if (kStation == 2 && kRing == 1) {
1308 kSerial = kChamber + 108;
1309 }
1310 if (kStation == 2 && kRing == 2) {
1311 kSerial = kChamber + 126;
1312 }
1313 if (kStation == 3 && kRing == 1) {
1314 kSerial = kChamber + 162;
1315 }
1316 if (kStation == 3 && kRing == 2) {
1317 kSerial = kChamber + 180;
1318 }
1319 if (kStation == 4 && kRing == 1) {
1320 kSerial = kChamber + 216;
1321 }
1322 if (kStation == 4 && kRing == 2) {
1323 kSerial = kChamber + 234;
1324 }
1325 if (kEndcap == 2) {
1326 kSerial = kSerial + 300;
1327 }
1328 return kSerial;
1329 }
1330
1331
1332 DEFINE_FWK_MODULE(CSCSkim);