Skims/src/CSCSkim.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    CSCSkim
0004 // Class:      CSCSkim
0005 //
0006 /**\class CSCSkim CSCSkim.cc RecoLocalMuon/CSCSkim/src/CSCSkim.cc
0007
0008  Description: Offline skim module for CSC cosmic ray data
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Michael Schmitt
0015 //         Created:  Sat Jul 12 17:43:33 CEST 2008
0016 //
0017 //
0018 //======================================================================
0019 //
0020 // CSCSkim:
0021 //
0022 // A simple skim module for extracting generally useful events from
0023 // the cosmic-ray runs (CRUZET-n and CRAFT).  The selected events
0024 // should be useful for most CSC-DPG and muon-POG studies.  However,
0025 // the selection requirements may bias the sample with respect to
0026 // trigger requirements and certain noise and efficiency-related
0027 // studies.
0028 //
0029 // Types of skims:   (typeOfSkim control word)
0030 //    1  = loose skim demanding hit chambers and/or segments
0031 //    2  = ask for hit chambers in both endcaps
0032 //    3  = segments in neighboring chambers - good for alignment
0033 //    4  = messy events
0034 //    5  = select events with DIGIs from one particular chamber
0035 //    6  = overlap with DT
0036 //    7  = nearly horizontal track going through ME1/1,2/1,3/1,4/1
0037 //    8  = ask for one long cosmic stand-alone muon track
0038 //    9  = selection for magnetic field studies
0039 //
0040 //
0041 //======================================================================
0042
0043 #include "DPGAnalysis/Skims/interface/CSCSkim.h"
0044
0045 using namespace std;
0046 using namespace edm;
0047
0048 //===================
0049 //  CONSTRUCTOR
0050 //===================
0051 CSCSkim::CSCSkim(const edm::ParameterSet& pset) : m_CSCGeomToken(esConsumes()) {
0052   // tokens from tags
0053
0054   // Really should define the wire and digi tags in config, but for now, to avoid having to modify
0055   // multiple config files, just hard-code those tags, to be equivalent to the pre-consumes code
0056
0057   //  wds_token = consumes<CSCWireDigiCollection>(pset.getParameter<InputTag>("simWireDigiTag"));
0058   //  sds_token = consumes<CSCStripDigiCollection>(pset.getParameter<InputTag>("simStripDigiTag"));
0059   //  wdr_token = consumes<CSCWireDigiCollection>(pset.getParameter<InputTag>("wireDigiTag"));
0060   //  sdr_token = consumes<CSCStripDigiCollection>(pset.getParameter<InputTag>("stripDigiTag"));
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   // Get the various input parameters
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   // for BStudy selection (skim type 9)
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 //  DESTRUCTOR
0113 //===================
0114 CSCSkim::~CSCSkim() {}
0115
0116 //================
0117 //  BEGIN JOB
0118 //================
0119 void CSCSkim::beginJob() {
0120   // set counters to zero
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     // Create the root file for the histograms
0136     theHistogramFile = new TFile(histogramFileName.c_str(), "RECREATE");
0137     theHistogramFile->cd();
0138
0139     if (makeHistograms) {
0140       // book histograms for the skimming module
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       // book histograms for the messy event skimming module
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 //  END JOB
0166 //================
0167 void CSCSkim::endJob() {
0168   // Write out results
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     // Write the histos to file
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 //  FILTER MAIN
0214 //================
0215 bool CSCSkim::filter(edm::Event& event, const edm::EventSetup& eventSetup) {
0216   // increment counter
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   // Get the CSC Geometry :
0225   ESHandle<CSCGeometry> cscGeom = eventSetup.getHandle(m_CSCGeomToken);
0226
0227   // Get the DIGI collections
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   // Get the RecHits collection :
0240   Handle<CSCRecHit2DCollection> cscRecHits;
0241   event.getByToken(rh_token, cscRecHits);
0242
0243   // get CSC segment collection
0244   Handle<CSCSegmentCollection> cscSegments;
0245   event.getByToken(seg_token, cscSegments);
0246
0247   // get the cosmic muons collection
0248   Handle<reco::TrackCollection> saMuons;
0249   if (typeOfSkim == 8) {
0250     event.getByToken(sam_token, saMuons);
0251   }
0252
0253   // get the stand-alone muons collection
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   // evaluate the skimming routines
0264   //======================================
0265
0266   // basic skimming
0267   bool basicEvent = false;
0268   if (typeOfSkim == 1 || typeOfSkim == 2) {
0269     basicEvent = doCSCSkimming(cscRecHits, cscSegments);
0270   }
0271
0272   // overlapping chamber skim
0273   bool goodOverlapEvent = false;
0274   if (typeOfSkim == 3) {
0275     goodOverlapEvent = doOverlapSkimming(cscSegments);
0276     if (goodOverlapEvent) {
0277       nEventsOverlappingChambers++;
0278     }
0279   }
0280
0281   // messy events skim
0282   bool messyEvent = false;
0283   if (typeOfSkim == 4) {
0284     messyEvent = doMessyEventSkimming(cscRecHits, cscSegments);
0285     if (messyEvent) {
0286       nEventsMessy++;
0287     }
0288   }
0289
0290   // select events with DIGIs in a certain chamber
0291   bool hasChamber = false;
0292   if (typeOfSkim == 5) {
0293     hasChamber = doCertainChamberSelection(wires, strips);
0294     if (hasChamber) {
0295       nEventsCertainChamber++;
0296     }
0297   }
0298
0299   // select events in the DT-CSC overlap region
0300   bool DTOverlapCandidate = false;
0301   if (typeOfSkim == 6) {
0302     DTOverlapCandidate = doDTOverlap(cscSegments);
0303     if (DTOverlapCandidate) {
0304       nEventsDTOverlap++;
0305     }
0306   }
0307
0308   // select halo-like events
0309   bool HaloLike = false;
0310   if (typeOfSkim == 7) {
0311     HaloLike = doHaloLike(cscSegments);
0312     if (HaloLike) {
0313       nEventsHaloLike++;
0314     }
0315   }
0316
0317   // select long cosmic tracks
0318   bool LongSATrack = false;
0319   if (typeOfSkim == 8) {
0320     LongSATrack = doLongSATrack(saMuons);
0321     if (LongSATrack) {
0322       nEventsLongSATrack++;
0323     }
0324   }
0325
0326   // select events suitable for a B-field study.  They have tracks in the tracker.
0327   bool GoodForBFieldStudy = false;
0328   if (typeOfSkim == 9) {
0329     GoodForBFieldStudy = doBFieldStudySelection(saMuons, tracks, gMuons);
0330     if (GoodForBFieldStudy) {
0331       nEventsForBFieldStudies++;
0332     }
0333   }
0334
0335   // set filter flag
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 // CSC Skimming
0372 //
0373 // ==============================================
0374
0375 bool CSCSkim::doCSCSkimming(edm::Handle<CSCRecHit2DCollection> cscRecHits,
0376                             edm::Handle<CSCSegmentCollection> cscSegments) {
0377   // how many RecHits in the collection?
0378   int nRecHits = cscRecHits->size();
0379
0380   // zero the recHit counter
0381   int cntRecHit[600];
0382   for (int i = 0; i < 600; i++) {
0383     cntRecHit[i] = 0;
0384   }
0385
0386   // ---------------------
0387   // Loop over rechits
0388   // ---------------------
0389
0390   CSCRecHit2DCollection::const_iterator recIt;
0391   for (recIt = cscRecHits->begin(); recIt != cscRecHits->end(); recIt++) {
0392     // which chamber is it?
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     // compute chamber serial number
0401     int kSerial = chamberSerial(kEndcap, kStation, kRing, kChamber);
0402
0403     // increment recHit counter
0404     //     (each layer is represented by a different power of 10)
0405     int kDigit = (int)pow((float)10., (float)(kLayer - 1));
0406     cntRecHit[kSerial] += kDigit;
0407
0408   }  //end rechit loop
0409
0410   // ------------------------------------------------------
0411   // Are there chambers with the minimum number of hits?
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   // how many Segments?
0442   int nSegments = cscSegments->size();
0443
0444   // ----------------------
0445   // fill histograms
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   // set the filter flag
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   // debug
0480   LogDebug("[CSCSkim]") << "----- nRecHits = " << nRecHits
0481                         << "\tnChambersWithMinimalHits = " << nChambersWithMinimalHits << "\tnSegments = " << nSegments
0482                         << "\tselect? " << selectEvent << std::endl;
0483
0484   /*
0485   if ((nChambersWithMinimalHitsPOS >= minimumHitChambers) && (nChambersWithMinimalHitsNEG >= minimumHitChambers)) {
0486     std::cout << "\n==========================================================================\n"
0487      << "\tinteresting event - chambers hit on both sides\n"
0488      << "\t  " <<  nEventsAnalyzed
0489      << "\trun " << iRun << "\tevent " << iEvent << std::endl;
0490     std::cout << "----- nRecHits = " << nRecHits
0491      << "\tnChambersWithMinimalHits = " << nChambersWithMinimalHits
0492      << "\tnSegments = " << nSegments
0493      << "\tselect? " << selectEvent << std::endl;
0494     for (int i = 0; i < 600; i++) {
0495       if (cntRecHit[i] > 0) {
0496     cout << "\t\t" << i << "\tcntRecHit= " << cntRecHit[i] << std::endl;
0497       }
0498     }
0499     std::cout << "==========================================================================\n\n" ;
0500   }
0501   */
0502
0503   return selectEvent;
0504 }
0505
0506 //-------------------------------------------------------------------------------
0507 // A module to select events useful in aligning chambers relative to each other
0508 // using the overlap regions at the edges (in Xlocal) of the chamber.
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   // how many Segments?
0516   //  int nSegments = cscSegments->size();
0517
0518   // zero arrays
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   // loop over segments
0528   // -----------------------
0529   for (CSCSegmentCollection::const_iterator it = cscSegments->begin(); it != cscSegments->end(); it++) {
0530     // which chamber?
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     // segment information
0539     float chisq = (*it).chi2();
0540     int nhits = (*it).nRecHits();
0541
0542     // is this a good segment?
0543     bool goodSegment = (nhits >= nhitsMinimum) && (chisq < chisqMaximum);
0544
0545     /*
0546     LocalPoint localPos = (*it).localPosition();
0547     float segX     = localPos.x();
0548     float segY     = localPos.y();
0549     std::cout << "E/S/R/Ch: " << kEndcap << "/" << kStation << "/" << kRing << "/" << kChamber
0550      << "\tnhits/chisq: " << nhits << "/" << chisq
0551      << "\tX/Y: " << segX << "/" << segY
0552      << "\tgood? " << goodSegment << std::endl;
0553     */
0554
0555     // count
0556     nAll[kSerial - 1]++;
0557     if (goodSegment)
0558       nGood[kSerial]++;
0559
0560   }  // end loop over segments
0561
0562   //----------------------
0563   // select the event
0564   //----------------------
0565
0566   // does any chamber have too many segments?
0567   bool messyChamber = false;
0568   for (int i = 0; i < 600; i++) {
0569     if (nAll[i] > nAllMaximum)
0570       messyChamber = true;
0571   }
0572
0573   // are there consecutive chambers with good segments
0574   // (This is a little sloppy but is probably fine for skimming...)
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 // This module selects events with a large numbere
0589 // of recHits and larger number of chambers with hits.
0590 //
0591 //============================================================
0592 bool CSCSkim::doMessyEventSkimming(edm::Handle<CSCRecHit2DCollection> cscRecHits,
0593                                    edm::Handle<CSCSegmentCollection> cscSegments) {
0594   // how many RecHits in the collection?
0595   int nRecHits = cscRecHits->size();
0596
0597   // zero the recHit counter
0598   int cntRecHit[600];
0599   for (int i = 0; i < 600; i++) {
0600     cntRecHit[i] = 0;
0601   }
0602
0603   // ---------------------
0604   // Loop over rechits
0605   // ---------------------
0606
0607   CSCRecHit2DCollection::const_iterator recIt;
0608   for (recIt = cscRecHits->begin(); recIt != cscRecHits->end(); recIt++) {
0609     // which chamber is it?
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     // compute chamber serial number
0618     int kSerial = chamberSerial(kEndcap, kStation, kRing, kChamber);
0619
0620     // increment recHit counter
0621     //     (each layer is represented by a different power of 10)
0622     int kDigit = (int)pow((float)10., (float)(kLayer - 1));
0623     cntRecHit[kSerial] += kDigit;
0624
0625   }  //end rechit loop
0626
0627   // ------------------------------------------------------
0628   // Are there chambers with the minimum number of hits?
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   // how many Segments?
0659   int nSegments = cscSegments->size();
0660
0661   // ----------------------
0662   // fill histograms
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   // set the filter flag
0689   // ----------------------
0690
0691   bool selectEvent = false;
0692   if ((nRecHits > 54) && (nChambersWithMinimalHits > 5)) {
0693     selectEvent = true;
0694   }
0695
0696   // debug
0697   LogDebug("[CSCSkim]") << "----- nRecHits = " << nRecHits
0698                         << "\tnChambersWithMinimalHits = " << nChambersWithMinimalHits << "\tnSegments = " << nSegments
0699                         << "\tselect? " << selectEvent << std::endl;
0700
0701   /*
0702   if (selectEvent) {
0703     std::cout << "\n==========================================================================\n"
0704      << "\tmessy event!\n"
0705      << "\t  " <<  nEventsAnalyzed
0706      << "\trun " << iRun << "\tevent " << iEvent << std::endl;
0707     std::cout << "----- nRecHits = " << nRecHits
0708      << "\tnChambersWithMinimalHits = " << nChambersWithMinimalHits
0709      << "\tnSegments = " << nSegments
0710      << "\tselect? " << selectEvent << std::endl;
0711     for (int i = 0; i < 600; i++) {
0712       if (cntRecHit[i] > 0) {
0713     cout << "\t\t" << i << "\tcntRecHit= " << cntRecHit[i] << std::endl;
0714       }
0715     }
0716     std::cout << "==========================================================================\n\n" ;
0717   }
0718   */
0719
0720   return selectEvent;
0721 }
0722
0723 //============================================================
0724 //
0725 // Select events with DIGIs are a particular chamber.
0726 //
0727 //============================================================
0728 bool CSCSkim::doCertainChamberSelection(edm::Handle<CSCWireDigiCollection> wires,
0729                                         edm::Handle<CSCStripDigiCollection> strips) {
0730   // Loop through the wire DIGIs, looking for a match
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   }  // end wire loop
0742
0743   // Loop through the strip DIGIs, looking for a match
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 // Select events which *might* probe the DT-CSC overlap region.
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   // initialize
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   // loop over segments
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     // which chamber?
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     // segment information
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     }  // this is a good segment
0833   }    // end loop over segments
0834
0835   // ---------------------------------------------
0836   // veto messy events
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   // check for relevant matchup of segments
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   if (matchup) {
0862     std::cout << "\tYYY looks like a good event.  Select!\n";
0863     std::cout << "-- pos endcap --\n"
0864      << "ME1/3: ";
0865     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEP13[k];}
0866     std::cout << "\nME2/2: ";
0867     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEP22[k];}
0868     std::cout << "\nME3/2: ";
0869     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEP32[k];}
0870     std::cout << std::endl;
0871   }
0872   */
0873
0874   // set the selection flag
0875   DTOverlapCandidate = matchup;
0876   return DTOverlapCandidate;
0877 }
0878
0879 //============================================================
0880 //
0881 // Select events which register in the inner parts of
0882 // stations 1, 2, 3 and 4.
0883 //
0884 //============================================================
0885 bool CSCSkim::doHaloLike(Handle<CSCSegmentCollection> cscSegments) {
0886   const float chisqMax = 100.;
0887   const int nhitsMin = 5;      // on a segment
0888   const int maxNSegments = 3;  // in a chamber
0889
0890   // initialize
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   // loop over segments
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     // which chamber?
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     // segment information
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     }  // this is a good segment
0975   }    // end loop over segments
0976
0977   // ---------------------------------------------
0978   // veto messy events
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   // check for relevant matchup of segments
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   if (matchup) {
1006     std::cout << "\tYYY looks like a good event.  Select!\n";
1007     std::cout << "-- pos endcap --\n"
1008      << "ME1/1: ";
1009     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEP11[k];}
1010     std::cout << "\nME1/2: ";
1011     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEP12[k];}
1012     std::cout << "\nME2/1: ";
1013     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEP21[k];}
1014     std::cout << "\nME3/1: ";
1015     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEP31[k];}
1016     std::cout << "\nME4/1: ";
1017     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEP41[k];}
1018     std::cout << std::endl;
1019     std::cout << "-- neg endcap --\n"
1020      << "ME1/1: ";
1021     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEN11[k];}
1022     std::cout << "\nME1/2: ";
1023     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEN12[k];}
1024     std::cout << "\nME2/1: ";
1025     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEN21[k];}
1026     std::cout << "\nME3/1: ";
1027     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEN31[k];}
1028     std::cout << "\nME4/1: ";
1029     for (int k=0; k<36; ++k) {std::cout << " " << setw(3) << cntMEN41[k];}
1030     std::cout << std::endl;
1031     std::cout << "\tn Analyzed = " << nEventsAnalyzed << "\tn Halo-like = " << nEventsHaloLike << std::endl;
1032   }
1033   */
1034
1035   // set the selection flag
1036   HaloLike = matchup;
1037   return HaloLike;
1038 }
1039
1040 //--------------------------------------------------------------
1041 // select events with at least one "long" stand-alone muon
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   // Loop through the track collection and test each one
1054   //
1055
1056   int nNiceMuons = 0;
1057
1058   for (reco::TrackCollection::const_iterator muon = saMuons->begin(); muon != saMuons->end(); ++muon) {
1059     // basic information
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     // loop over hits
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           //DTChamberId dtId(detId.rawId());
1078           //int chamberId = dtId.sector();
1079           nDTHits++;
1080         } else if (detId.subdetId() == MuonSubdetId::CSC) {
1081           //CSCDetId cscId(detId.rawId());
1082           //int chamberId = cscId.chamber();
1083           nCSCHits++;
1084         }
1085       }
1086     }
1087
1088     // is this a nice muon?
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 // Select events which are good for B-field studies.
1104 //
1105 // These events have a good track in the tracker.
1106 //
1107 //  D.Dibur and M.Schmitt
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   // examine the stand-alone tracks
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     }  // end loop over hits
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   }  // end loop over stand-alone muon collection
1188
1189   //-----------------------------------
1190   // examine the tracker tracks
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   }  // end loop over tracker tracks
1239
1240   //-----------------------------------
1241   // examine the global muons
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       // int nTotalHits = hp.numberOfHits();
1250       //    int nValidHits = hp.numberOfValidHits();
1251       int nTrackerHits = hp.numberOfValidTrackerHits();
1252       // int nPixelHits   = hp.numberOfValidPixelHits();
1253       // int nStripHits   = hp.numberOfValidStripHits();
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       }  // end loop over hits
1269
1270       bool goodGlobalMuon =
1271           (pDef > pMin) && (nTrackerHits >= nValidHitsMin) && (nCSCHits >= nCSCHitsMin) && (redChiSq < redChiSqMax);
1272
1273       if (goodGlobalMuon) {
1274         nGoodGlobalMuons++;
1275       }
1276
1277     }  // this is a global muon
1278   }    // end loop over stand-alone muon collection
1279
1280   //-----------------------------------
1281   // do we accept this event?
1282   //-----------------------------------
1283
1284   acceptThisEvent = ((nGoodSAMuons > 0) && (nGoodTracks > 0)) || (nGoodGlobalMuons > 0);
1285
1286   return acceptThisEvent;
1287 }
1288
1289 //--------------------------------------------------------------
1290 // Compute a serial number for the chamber.
1291 // This is useful when filling histograms and working with arrays.
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   }  // one day...
1325   if (kEndcap == 2) {
1326     kSerial = kSerial + 300;
1327   }
1328   return kSerial;
1329 }
1330
1331 //define this as a plug-in
1332 DEFINE_FWK_MODULE(CSCSkim);