Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:38 

0001 /*
0002  *  Routine to calculate CSC efficiencies
0003  *  Comments about the program logic are denoted by //----
0004  *
0005  *  Stoyan Stoynev, Northwestern University.
0006  */
0007 
0008 #include "RecoLocalMuon/CSCEfficiency/src/CSCEfficiency.h"
0009 
0010 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0011 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0012 #include "FWCore/Common/interface/TriggerNames.h"
0013 #include "DataFormats/MuonReco/interface/Muon.h"
0014 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0015 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0016 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0017 
0018 #include "DataFormats/Provenance/interface/RunLumiEventNumber.h"
0019 
0020 using namespace std;
0021 
0022 bool CSCEfficiency::filter(edm::Event &event, const edm::EventSetup &eventSetup) {
0023   passTheEvent = false;
0024   DataFlow->Fill(0.);
0025   MuonPatternRecoDumper debug;
0026 
0027   //---- increment counter
0028   nEventsAnalyzed++;
0029   // printalot debug output
0030   printalot = (nEventsAnalyzed < int(printout_NEvents));  //
0031   edm::RunNumber_t const iRun = event.id().run();
0032   edm::EventNumber_t const iEvent = event.id().event();
0033   if (0 == fmod(double(nEventsAnalyzed), double(1000))) {
0034     if (printalot) {
0035       printf("\n==enter==CSCEfficiency===== run %u\tevent %llu\tn Analyzed %i\n", iRun, iEvent, nEventsAnalyzed);
0036     }
0037   }
0038   theService->update(eventSetup);
0039   //---- These declarations create handles to the types of records that you want
0040   //---- to retrieve from event "e".
0041   if (printalot)
0042     printf("\tget handles for digi collections\n");
0043 
0044   //---- Pass the handle to the method "getByType", which is used to retrieve
0045   //---- one and only one instance of the type in question out of event "e". If
0046   //---- zero or more than one instance exists in the event an exception is thrown.
0047   if (printalot)
0048     printf("\tpass handles\n");
0049   edm::Handle<CSCALCTDigiCollection> alcts;
0050   edm::Handle<CSCCLCTDigiCollection> clcts;
0051   edm::Handle<CSCCorrelatedLCTDigiCollection> correlatedlcts;
0052   edm::Handle<CSCWireDigiCollection> wires;
0053   edm::Handle<CSCStripDigiCollection> strips;
0054   edm::Handle<CSCRecHit2DCollection> rechits;
0055   edm::Handle<CSCSegmentCollection> segments;
0056   edm::Handle<edm::View<reco::Track> > trackCollectionH;
0057   edm::Handle<edm::PSimHitContainer> simhits;
0058 
0059   if (useDigis) {
0060     event.getByToken(wd_token, wires);
0061     event.getByToken(sd_token, strips);
0062     event.getByToken(al_token, alcts);
0063     event.getByToken(cl_token, clcts);
0064     event.getByToken(co_token, correlatedlcts);
0065   }
0066   if (!isData) {
0067     event.getByToken(sh_token, simhits);
0068   }
0069   event.getByToken(rh_token, rechits);
0070   event.getByToken(se_token, segments);
0071   event.getByToken(tk_token, trackCollectionH);
0072   const edm::View<reco::Track> trackCollection = *(trackCollectionH.product());
0073 
0074   //---- Get the CSC Geometry :
0075   if (printalot)
0076     printf("\tget the CSC geometry.\n");
0077   edm::ESHandle<CSCGeometry> cscGeom = eventSetup.getHandle(geomToken_);
0078 
0079   // use theTrackingGeometry instead of cscGeom?
0080   //edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry = eventSetup.getHandle(trackingGeomToken_);
0081 
0082   bool triggerPassed = true;
0083   if (useTrigger) {
0084     // access the trigger information
0085     // trigger names can be find in HLTrigger/Configuration/python/HLT_2E30_cff.py (or?)
0086     // get hold of TriggerResults
0087     edm::Handle<edm::TriggerResults> hltR;
0088     event.getByToken(ht_token, hltR);
0089     const edm::TriggerNames &triggerNames = event.triggerNames(*hltR);
0090     triggerPassed = applyTrigger(hltR, triggerNames);
0091   }
0092   if (!triggerPassed) {
0093     return triggerPassed;
0094   }
0095   DataFlow->Fill(1.);
0096   GlobalPoint gpZero(0., 0., 0.);
0097   if (theService->magneticField()->inTesla(gpZero).mag2() < 0.1) {
0098     magField = false;
0099   } else {
0100     magField = true;
0101   }
0102 
0103   //---- store info from digis
0104   fillDigiInfo(alcts, clcts, correlatedlcts, wires, strips, simhits, rechits, segments, cscGeom);
0105   //
0106   edm::Handle<reco::MuonCollection> muons;
0107   edm::InputTag muonTag_("muons");
0108   event.getByLabel(muonTag_, muons);
0109 
0110   edm::Handle<reco::BeamSpot> beamSpotHandle;
0111   event.getByLabel("offlineBeamSpot", beamSpotHandle);
0112   reco::BeamSpot vertexBeamSpot = *beamSpotHandle;
0113   //
0114   std::vector<reco::MuonCollection::const_iterator> goodMuons_it;
0115   unsigned int nPositiveZ = 0;
0116   unsigned int nNegativeZ = 0;
0117   float muonOuterZPosition = -99999.;
0118   if (isIPdata) {
0119     if (printalot)
0120       std::cout << " muons.size() = " << muons->size() << std::endl;
0121     for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0122       DataFlow->Fill(31.);
0123       if (printalot) {
0124         std::cout << "  iMuon = " << muon - muons->begin() << " charge = " << muon->charge() << " p = " << muon->p()
0125                   << " pt = " << muon->pt() << " eta = " << muon->eta() << " phi = " << muon->phi()
0126                   << " matches = " << muon->matches().size()
0127                   << " matched Seg = " << muon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration)
0128                   << " GLB/TR/STA = " << muon->isGlobalMuon() << "/" << muon->isTrackerMuon() << "/"
0129                   << muon->isStandAloneMuon() << std::endl;
0130       }
0131       if (!(muon->isTrackerMuon() && muon->isGlobalMuon())) {
0132         continue;
0133       }
0134       DataFlow->Fill(32.);
0135       double relISO =
0136           (muon->isolationR03().sumPt + muon->isolationR03().emEt + muon->isolationR03().hadEt) / muon->track()->pt();
0137       if (printalot) {
0138         std::cout << " relISO = " << relISO << " emVetoEt = " << muon->isolationR03().emVetoEt
0139                   << " caloComp = " << muon::caloCompatibility(*(muon))
0140                   << " dxy = " << fabs(muon->track()->dxy(vertexBeamSpot.position())) << std::endl;
0141       }
0142       if (
0143           //relISO>0.1 || muon::caloCompatibility(*(muon))<.90 ||
0144           fabs(muon->track()->dxy(vertexBeamSpot.position())) > 0.2 || muon->pt() < 6.) {
0145         continue;
0146       }
0147       DataFlow->Fill(33.);
0148       if (muon->track()->hitPattern().numberOfValidPixelHits() < 1 ||
0149           muon->track()->hitPattern().numberOfValidTrackerHits() < 11 ||
0150           muon->combinedMuon()->hitPattern().numberOfValidMuonHits() < 1 ||
0151           muon->combinedMuon()->normalizedChi2() > 10. || muon->numberOfMatches() < 2) {
0152         continue;
0153       }
0154       DataFlow->Fill(34.);
0155       float zOuter = muon->combinedMuon()->outerPosition().z();
0156       float rhoOuter = muon->combinedMuon()->outerPosition().rho();
0157       bool passDepth = true;
0158       // barrel region
0159       //if ( fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 480.){
0160       if (fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 540.) {
0161         passDepth = false;
0162       }
0163       // endcap region
0164       //else if( fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.){
0165       else if (fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.) {
0166         passDepth = false;
0167       }
0168       // overlap region
0169       //else if ( fabs(zOuter) > 680. && fabs(zOuter) < 730. && rhoOuter < 480.){
0170       else if (fabs(zOuter) > 680. && fabs(zOuter) < 880. && rhoOuter < 540.) {
0171         passDepth = false;
0172       }
0173       if (!passDepth) {
0174         continue;
0175       }
0176       DataFlow->Fill(35.);
0177       goodMuons_it.push_back(muon);
0178       if (muon->track()->momentum().z() > 0.) {
0179         ++nPositiveZ;
0180       }
0181       if (muon->track()->momentum().z() < 0.) {
0182         ++nNegativeZ;
0183       }
0184     }
0185   }
0186 
0187   //
0188 
0189   if (printalot)
0190     std::cout << "Start track loop over " << trackCollection.size() << " tracks" << std::endl;
0191   for (edm::View<reco::Track>::size_type i = 0; i < trackCollection.size(); ++i) {
0192     DataFlow->Fill(2.);
0193     edm::RefToBase<reco::Track> track(trackCollectionH, i);
0194     //std::cout<<" iTR = "<<i<<" eta = "<<track->eta()<<" phi = "<<track->phi()<<std::cout<<" pt = "<<track->pt()<<std::endl;
0195     if (isIPdata) {
0196       if (printalot) {
0197         std::cout << " nNegativeZ = " << nNegativeZ << " nPositiveZ = " << nPositiveZ << std::endl;
0198       }
0199       if (nNegativeZ > 1 || nPositiveZ > 1) {
0200         break;
0201       }
0202       bool trackOK = false;
0203       if (printalot) {
0204         std::cout << " goodMuons_it.size() = " << goodMuons_it.size() << std::endl;
0205       }
0206       for (size_t iM = 0; iM < goodMuons_it.size(); ++iM) {
0207         //std::cout<<" iM = "<<iM<<" eta = "<<goodMuons_it[iM]->track()->eta()<<
0208         //" phi = "<<goodMuons_it[iM]->track()->phi()<<
0209         //" pt = "<<goodMuons_it[iM]->track()->pt()<<std::endl;
0210         float deltaR = pow(track->phi() - goodMuons_it[iM]->track()->phi(), 2) +
0211                        pow(track->eta() - goodMuons_it[iM]->track()->eta(), 2);
0212         deltaR = sqrt(deltaR);
0213         if (printalot) {
0214           std::cout << " TR mu match to a tr: deltaR = " << deltaR
0215                     << " dPt = " << track->pt() - goodMuons_it[iM]->track()->pt() << std::endl;
0216         }
0217         if (deltaR > 0.01 || fabs(track->pt() - goodMuons_it[iM]->track()->pt()) > 0.1) {
0218           continue;
0219         } else {
0220           trackOK = true;
0221           if (printalot) {
0222             std::cout << " trackOK " << std::endl;
0223           }
0224           muonOuterZPosition = goodMuons_it[iM]->combinedMuon()->outerPosition().z();
0225           break;
0226           //++nChosenTracks;
0227         }
0228       }
0229       if (!trackOK) {
0230         if (printalot) {
0231           std::cout << " failed: trackOK " << std::endl;
0232         }
0233         continue;
0234       }
0235     } else {
0236       //---- Do we need a better "clean track" definition?
0237       if (trackCollection.size() > 2) {
0238         break;
0239       }
0240       DataFlow->Fill(3.);
0241       if (!i && 2 == trackCollection.size()) {
0242         edm::View<reco::Track>::size_type tType = 1;
0243         edm::RefToBase<reco::Track> trackTwo(trackCollectionH, tType);
0244         if (track->outerPosition().z() * trackTwo->outerPosition().z() > 0) {  // in one and the same "endcap"
0245           break;
0246         }
0247       }
0248     }
0249     DataFlow->Fill(4.);
0250     if (printalot) {
0251       std::cout << "i track = " << i << " P = " << track->p() << " chi2/ndf = " << track->normalizedChi2()
0252                 << " nSeg = " << segments->size() << std::endl;
0253       std::cout << "quality undef/loose/tight/high/confirmed/goodIt/size " << track->quality(reco::Track::undefQuality)
0254                 << "/" << track->quality(reco::Track::loose) << "/" << track->quality(reco::Track::tight) << "/"
0255                 << track->quality(reco::Track::highPurity) << "/" << track->quality(reco::Track::confirmed) << "/"
0256                 << track->quality(reco::Track::goodIterative) << "/" << track->quality(reco::Track::qualitySize)
0257                 << std::endl;
0258       std::cout << " pt = " << track->pt() << " +-" << track->ptError() << " q/pt = " << track->qoverp() << " +- "
0259                 << track->qoverpError() << std::endl;
0260       //std::cout<<" const Pmin = "<<minTrackMomentum<<" pMax = "<<maxTrackMomentum<<" maxNormChi2 = "<<maxNormChi2<<std::endl;
0261       std::cout << " track inner position = " << track->innerPosition()
0262                 << " outer position = " << track->outerPosition() << std::endl;
0263       std::cout << "track eta (outer) = " << track->outerPosition().eta()
0264                 << " phi (outer) = " << track->outerPosition().phi() << std::endl;
0265       if (fabs(track->innerPosition().z()) > 500.) {
0266         DetId innerDetId(track->innerDetId());
0267         std::cout << " dump inner state MUON detid  = " << debug.dumpMuonId(innerDetId) << std::endl;
0268       }
0269       if (fabs(track->outerPosition().z()) > 500.) {
0270         DetId outerDetId(track->outerDetId());
0271         std::cout << " dump outer state MUON detid  = " << debug.dumpMuonId(outerDetId) << std::endl;
0272       }
0273 
0274       std::cout << " nHits = " << track->found() << std::endl;
0275       /*
0276       trackingRecHit_iterator rhbegin = track->recHitsBegin();
0277       trackingRecHit_iterator rhend = track->recHitsEnd();
0278       int iRH = 0;
0279       for(trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit){
0280         const GeomDet* geomDet = theTrackingGeometry->idToDet((*recHit)->geographicalId());
0281     std::cout<<"hit "<<iRH<<" loc pos = " <<(*recHit)->localPosition()<<
0282       " glob pos = " <<geomDet->toGlobal((*recHit)->localPosition())<<std::endl;
0283         ++iRH;
0284       }
0285       */
0286     }
0287     float dpT_ov_pT = 0.;
0288     if (fabs(track->pt()) > 0.001) {
0289       dpT_ov_pT = track->ptError() / track->pt();
0290     }
0291     //---- These define a "good" track
0292     if (track->normalizedChi2() > maxNormChi2) {  // quality
0293       break;
0294     }
0295     DataFlow->Fill(5.);
0296     if (track->found() < minTrackHits) {  // enough data points
0297       break;
0298     }
0299     DataFlow->Fill(6.);
0300     if (!segments->size()) {  // better have something in the CSC
0301       break;
0302     }
0303     DataFlow->Fill(7.);
0304     if (magField && (track->p() < minP || track->p() > maxP)) {  // proper energy range
0305       break;
0306     }
0307     DataFlow->Fill(8.);
0308     if (magField && (dpT_ov_pT > 0.5)) {  // not too crazy uncertainty
0309       break;
0310     }
0311     DataFlow->Fill(9.);
0312 
0313     passTheEvent = true;
0314     if (printalot)
0315       std::cout << "good Track" << std::endl;
0316     CLHEP::Hep3Vector r3T_inner(track->innerPosition().x(), track->innerPosition().y(), track->innerPosition().z());
0317     CLHEP::Hep3Vector r3T(track->outerPosition().x(), track->outerPosition().y(), track->outerPosition().z());
0318     chooseDirection(r3T_inner, r3T);  // for non-IP
0319 
0320     CLHEP::Hep3Vector p3T(track->outerMomentum().x(), track->outerMomentum().y(), track->outerMomentum().z());
0321     CLHEP::Hep3Vector p3_propagated, r3_propagated;
0322     AlgebraicSymMatrix66 cov_propagated, covT;
0323     covT *= 1e-20;
0324     cov_propagated *= 1e-20;
0325     int charge = track->charge();
0326     FreeTrajectoryState ftsStart = getFromCLHEP(p3T, r3T, charge, covT, &*(theService->magneticField()));
0327     if (printalot) {
0328       std::cout << " p = " << track->p() << " norm chi2 = " << track->normalizedChi2() << std::endl;
0329       std::cout << " dump the very first FTS  = " << debug.dumpFTS(ftsStart) << std::endl;
0330     }
0331     TrajectoryStateOnSurface tSOSDest;
0332     int endcap = 0;
0333     //---- which endcap to look at
0334     if (track->outerPosition().z() > 0) {
0335       endcap = 1;
0336     } else {
0337       endcap = 2;
0338     }
0339     int chamber = 1;
0340     //---- a "refference" CSCDetId for each ring
0341     std::vector<CSCDetId> refME;
0342     for (int iS = 1; iS < 5; ++iS) {
0343       for (int iR = 1; iR < 4; ++iR) {
0344         if (1 != iS && iR > 2) {
0345           continue;
0346         } else if (4 == iS && iR > 1) {
0347           continue;
0348         }
0349         refME.push_back(CSCDetId(endcap, iS, iR, chamber));
0350       }
0351     }
0352     //---- loop over the "refference" CSCDetIds
0353     for (size_t iSt = 0; iSt < refME.size(); ++iSt) {
0354       if (printalot) {
0355         std::cout << "loop iStatation = " << iSt << std::endl;
0356         std::cout << "refME[iSt]: st = " << refME[iSt].station() << " rg = " << refME[iSt].ring() << std::endl;
0357       }
0358       std::map<std::string, bool> chamberTypes;
0359       chamberTypes["ME11"] = false;
0360       chamberTypes["ME12"] = false;
0361       chamberTypes["ME13"] = false;
0362       chamberTypes["ME21"] = false;
0363       chamberTypes["ME22"] = false;
0364       chamberTypes["ME31"] = false;
0365       chamberTypes["ME32"] = false;
0366       chamberTypes["ME41"] = false;
0367       const CSCChamber *cscChamber_base = cscGeom->chamber(refME[iSt].chamberId());
0368       DetId detId = cscChamber_base->geographicalId();
0369       if (printalot) {
0370         std::cout << " base iStation : eta = " << cscGeom->idToDet(detId)->surface().position().eta()
0371                   << " phi = " << cscGeom->idToDet(detId)->surface().position().phi()
0372                   << " y = " << cscGeom->idToDet(detId)->surface().position().y() << std::endl;
0373         std::cout << " dump base iStation detid  = " << debug.dumpMuonId(detId) << std::endl;
0374         std::cout << " dump FTS start  = " << debug.dumpFTS(ftsStart) << std::endl;
0375       }
0376       //---- propagate to this ME
0377       tSOSDest = propagate(ftsStart, cscGeom->idToDet(detId)->surface());
0378       if (tSOSDest.isValid()) {
0379         ftsStart = *tSOSDest.freeState();
0380         if (printalot)
0381           std::cout << "  dump FTS end   = " << debug.dumpFTS(ftsStart) << std::endl;
0382         getFromFTS(ftsStart, p3_propagated, r3_propagated, charge, cov_propagated);
0383         float feta = fabs(r3_propagated.eta());
0384         float phi = r3_propagated.phi();
0385         //---- which rings are (possibly) penetrated
0386         ringCandidates(refME[iSt].station(), feta, chamberTypes);
0387 
0388         map<std::string, bool>::iterator iter;
0389         int iterations = 0;
0390         //---- loop over ring candidates
0391         for (iter = chamberTypes.begin(); iter != chamberTypes.end(); iter++) {
0392           ++iterations;
0393           //---- is this ME a machinig candidate station
0394           if (iter->second && (iterations - 1) == int(iSt)) {
0395             if (printalot) {
0396               std::cout << " Chamber type " << iter->first << " is a candidate..." << std::endl;
0397               std::cout << " station() = " << refME[iSt].station() << " ring() = " << refME[iSt].ring()
0398                         << " iSt = " << iSt << std::endl;
0399             }
0400             std::vector<int> coupleOfChambers;
0401             //---- which chamber (and its closes neighbor) is penetrated by the track - candidates
0402             chamberCandidates(refME[iSt].station(), refME[iSt].ring(), phi, coupleOfChambers);
0403             //---- loop over the two chamber candidates
0404             for (size_t iCh = 0; iCh < coupleOfChambers.size(); ++iCh) {
0405               DataFlow->Fill(11.);
0406               if (printalot)
0407                 std::cout << " Check chamber N = " << coupleOfChambers.at(iCh) << std::endl;
0408               ;
0409               if ((!getAbsoluteEfficiency) &&
0410                   (true == emptyChambers[refME[iSt].endcap() - 1][refME[iSt].station() - 1][refME[iSt].ring() - 1]
0411                                         [coupleOfChambers.at(iCh) - FirstCh])) {
0412                 continue;
0413               }
0414               CSCDetId theCSCId(refME[iSt].endcap(), refME[iSt].station(), refME[iSt].ring(), coupleOfChambers.at(iCh));
0415               const CSCChamber *cscChamber = cscGeom->chamber(theCSCId.chamberId());
0416               const BoundPlane bpCh = cscGeom->idToDet(cscChamber->geographicalId())->surface();
0417               float zFTS = ftsStart.position().z();
0418               float dz = fabs(bpCh.position().z() - zFTS);
0419               float zDistInner = track->innerPosition().z() - bpCh.position().z();
0420               float zDistOuter = track->outerPosition().z() - bpCh.position().z();
0421               //---- only detectors between the inner and outer points of the track are considered for non IP-data
0422               if (printalot) {
0423                 std::cout << " zIn = " << track->innerPosition().z() << " zOut = " << track->outerPosition().z()
0424                           << " zSurf = " << bpCh.position().z() << std::endl;
0425               }
0426               if (!isIPdata && (zDistInner * zDistOuter > 0. || fabs(zDistInner) < 15. ||
0427                                 fabs(zDistOuter) < 15.)) {  // for non IP-data
0428                 if (printalot) {
0429                   std::cout << " Not an intermediate (as defined) point... Skip." << std::endl;
0430                 }
0431                 continue;
0432               }
0433               if (isIPdata && fabs(track->eta()) < 1.8) {
0434                 if (fabs(muonOuterZPosition) - fabs(bpCh.position().z()) < 0 ||
0435                     fabs(muonOuterZPosition - bpCh.position().z()) < 15.) {
0436                   continue;
0437                 }
0438               }
0439               DataFlow->Fill(13.);
0440               //---- propagate to the chamber (from this ME) if it is a different surface (odd/even chambers)
0441               if (dz > 0.1) {  // i.e. non-zero (float 0 check is bad)
0442                 //if(fabs(zChanmber - zFTS ) > 0.1){
0443                 tSOSDest = propagate(ftsStart, cscGeom->idToDet(cscChamber->geographicalId())->surface());
0444                 if (tSOSDest.isValid()) {
0445                   ftsStart = *tSOSDest.freeState();
0446                 } else {
0447                   if (printalot)
0448                     std::cout << "TSOS not valid! Break." << std::endl;
0449                   break;
0450                 }
0451               } else {
0452                 if (printalot)
0453                   std::cout << " info: dz<0.1" << std::endl;
0454               }
0455               DataFlow->Fill(15.);
0456               FreeTrajectoryState ftsInit = ftsStart;
0457               bool inDeadZone = false;
0458               //---- loop over the 6 layers
0459               for (int iLayer = 0; iLayer < 6; ++iLayer) {
0460                 bool extrapolationPassed = true;
0461                 if (printalot) {
0462                   std::cout << " iLayer = " << iLayer << "   dump FTS init  = " << debug.dumpFTS(ftsInit) << std::endl;
0463                   std::cout << " dump detid  = " << debug.dumpMuonId(cscChamber->geographicalId()) << std::endl;
0464                   std::cout << "Surface to propagate to:  pos = " << cscChamber->layer(iLayer + 1)->surface().position()
0465                             << " eta = " << cscChamber->layer(iLayer + 1)->surface().position().eta()
0466                             << " phi = " << cscChamber->layer(iLayer + 1)->surface().position().phi() << std::endl;
0467                 }
0468                 //---- propagate to this layer
0469                 tSOSDest = propagate(ftsInit, cscChamber->layer(iLayer + 1)->surface());
0470                 if (tSOSDest.isValid()) {
0471                   ftsInit = *tSOSDest.freeState();
0472                   if (printalot)
0473                     std::cout << " Propagation between layers successful:  dump FTS end  = " << debug.dumpFTS(ftsInit)
0474                               << std::endl;
0475                   getFromFTS(ftsInit, p3_propagated, r3_propagated, charge, cov_propagated);
0476                 } else {
0477                   if (printalot)
0478                     std::cout << "Propagation between layers not successful - notValid TSOS" << std::endl;
0479                   extrapolationPassed = false;
0480                   inDeadZone = true;
0481                 }
0482                 //}
0483                 //---- Extrapolation passed? For each layer?
0484                 if (extrapolationPassed) {
0485                   GlobalPoint theExtrapolationPoint(r3_propagated.x(), r3_propagated.y(), r3_propagated.z());
0486                   LocalPoint theLocalPoint = cscChamber->layer(iLayer + 1)->toLocal(theExtrapolationPoint);
0487                   //std::cout<<" Candidate chamber: extrapolated LocalPoint = "<<theLocalPoint<<std::endl;
0488                   inDeadZone = (inDeadZone ||
0489                                 !inSensitiveLocalRegion(
0490                                     theLocalPoint.x(), theLocalPoint.y(), refME[iSt].station(), refME[iSt].ring()));
0491                   if (printalot) {
0492                     std::cout << " Candidate chamber: extrapolated LocalPoint = " << theLocalPoint
0493                               << "inDeadZone = " << inDeadZone << std::endl;
0494                   }
0495                   //---- break if in dead zone for any layer ("clean" tracks)
0496                   if (inDeadZone) {
0497                     break;
0498                   }
0499                 } else {
0500                   break;
0501                 }
0502               }
0503               DataFlow->Fill(17.);
0504               //---- Is a track in a sensitive area for each layer?
0505               if (!inDeadZone) {  //---- for any layer
0506                 DataFlow->Fill(19.);
0507                 if (printalot)
0508                   std::cout << "Do efficiencies..." << std::endl;
0509                 //---- Do efficiencies
0510                 // angle cuts applied (if configured)
0511                 bool angle_flag = true;
0512                 angle_flag = efficienciesPerChamber(theCSCId, cscChamber, ftsStart);
0513                 if (useDigis && angle_flag) {
0514                   stripWire_Efficiencies(theCSCId, ftsStart);
0515                 }
0516                 if (angle_flag) {
0517                   recHitSegment_Efficiencies(theCSCId, cscChamber, ftsStart);
0518                   if (!isData) {
0519                     recSimHitEfficiency(theCSCId, ftsStart);
0520                   }
0521                 }
0522               } else {
0523                 if (printalot)
0524                   std::cout << " Not in active area for all layers" << std::endl;
0525               }
0526             }
0527             if (tSOSDest.isValid()) {
0528               ftsStart = *tSOSDest.freeState();
0529             }
0530           }
0531         }
0532       } else {
0533         if (printalot)
0534           std::cout << " TSOS not valid..." << std::endl;
0535       }
0536     }
0537   }
0538   //---- End
0539   if (printalot)
0540     printf("==exit===CSCEfficiency===== run %u\tevent %llu\n\n", iRun, iEvent);
0541   return passTheEvent;
0542 }
0543 
0544 //
0545 bool CSCEfficiency::inSensitiveLocalRegion(double xLocal, double yLocal, int station, int ring) {
0546   //---- Good region means sensitive area of a chamber. "Local" stands for the local system
0547   bool pass = false;
0548   std::vector<double> chamberBounds(3);  // the sensitive area
0549   float y_center = 99999.;
0550   //---- hardcoded... not good
0551   if (station > 1 && station < 5) {
0552     if (2 == ring) {
0553       chamberBounds[0] = 66.46 / 2;   // (+-)x1 shorter
0554       chamberBounds[1] = 127.15 / 2;  // (+-)x2 longer
0555       chamberBounds[2] = 323.06 / 2;
0556       y_center = -0.95;
0557     } else {
0558       if (2 == station) {
0559         chamberBounds[0] = 54.00 / 2;   // (+-)x1 shorter
0560         chamberBounds[1] = 125.71 / 2;  // (+-)x2 longer
0561         chamberBounds[2] = 189.66 / 2;
0562         y_center = -0.955;
0563       } else if (3 == station) {
0564         chamberBounds[0] = 61.40 / 2;   // (+-)x1 shorter
0565         chamberBounds[1] = 125.71 / 2;  // (+-)x2 longer
0566         chamberBounds[2] = 169.70 / 2;
0567         y_center = -0.97;
0568       } else if (4 == station) {
0569         chamberBounds[0] = 69.01 / 2;   // (+-)x1 shorter
0570         chamberBounds[1] = 125.65 / 2;  // (+-)x2 longer
0571         chamberBounds[2] = 149.42 / 2;
0572         y_center = -0.94;
0573       }
0574     }
0575   } else if (1 == station) {
0576     if (3 == ring) {
0577       chamberBounds[0] = 63.40 / 2;  // (+-)x1 shorter
0578       chamberBounds[1] = 92.10 / 2;  // (+-)x2 longer
0579       chamberBounds[2] = 164.16 / 2;
0580       y_center = -1.075;
0581     } else if (2 == ring) {
0582       chamberBounds[0] = 51.00 / 2;  // (+-)x1 shorter
0583       chamberBounds[1] = 83.74 / 2;  // (+-)x2 longer
0584       chamberBounds[2] = 174.49 / 2;
0585       y_center = -0.96;
0586     } else {                        // to be investigated
0587       chamberBounds[0] = 30. / 2;   //40./2; // (+-)x1 shorter
0588       chamberBounds[1] = 60. / 2;   //100./2; // (+-)x2 longer
0589       chamberBounds[2] = 160. / 2;  //142./2;
0590       y_center = 0.;
0591     }
0592   }
0593   double yUp = chamberBounds[2] + y_center;
0594   double yDown = -chamberBounds[2] + y_center;
0595   double xBound1Shifted = chamberBounds[0] - distanceFromDeadZone;  //
0596   double xBound2Shifted = chamberBounds[1] - distanceFromDeadZone;  //
0597   double lineSlope = (yUp - yDown) / (xBound2Shifted - xBound1Shifted);
0598   double lineConst = yUp - lineSlope * xBound2Shifted;
0599   double yBoundary = lineSlope * abs(xLocal) + lineConst;
0600   pass = checkLocal(yLocal, yBoundary, station, ring);
0601   return pass;
0602 }
0603 
0604 bool CSCEfficiency::checkLocal(double yLocal, double yBoundary, int station, int ring) {
0605   //---- check if it is in a good local region (sensitive area - geometrical and HV boundaries excluded)
0606   bool pass = false;
0607   std::vector<float> deadZoneCenter(6);
0608   const float deadZoneHalf = 0.32 * 7 / 2;              // wire spacing * (wires missing + 1)/2
0609   float cutZone = deadZoneHalf + distanceFromDeadZone;  //cm
0610   //---- hardcoded... not good
0611   if (station > 1 && station < 5) {
0612     if (2 == ring) {
0613       deadZoneCenter[0] = -162.48;
0614       deadZoneCenter[1] = -81.8744;
0615       deadZoneCenter[2] = -21.18165;
0616       deadZoneCenter[3] = 39.51105;
0617       deadZoneCenter[4] = 100.2939;
0618       deadZoneCenter[5] = 160.58;
0619 
0620       if (yLocal > yBoundary && ((yLocal > deadZoneCenter[0] + cutZone && yLocal < deadZoneCenter[1] - cutZone) ||
0621                                  (yLocal > deadZoneCenter[1] + cutZone && yLocal < deadZoneCenter[2] - cutZone) ||
0622                                  (yLocal > deadZoneCenter[2] + cutZone && yLocal < deadZoneCenter[3] - cutZone) ||
0623                                  (yLocal > deadZoneCenter[3] + cutZone && yLocal < deadZoneCenter[4] - cutZone) ||
0624                                  (yLocal > deadZoneCenter[4] + cutZone && yLocal < deadZoneCenter[5] - cutZone))) {
0625         pass = true;
0626       }
0627     } else if (1 == ring) {
0628       if (2 == station) {
0629         deadZoneCenter[0] = -95.94;
0630         deadZoneCenter[1] = -27.47;
0631         deadZoneCenter[2] = 33.67;
0632         deadZoneCenter[3] = 93.72;
0633       } else if (3 == station) {
0634         deadZoneCenter[0] = -85.97;
0635         deadZoneCenter[1] = -36.21;
0636         deadZoneCenter[2] = 23.68;
0637         deadZoneCenter[3] = 84.04;
0638       } else if (4 == station) {
0639         deadZoneCenter[0] = -75.82;
0640         deadZoneCenter[1] = -26.14;
0641         deadZoneCenter[2] = 23.85;
0642         deadZoneCenter[3] = 73.91;
0643       }
0644       if (yLocal > yBoundary && ((yLocal > deadZoneCenter[0] + cutZone && yLocal < deadZoneCenter[1] - cutZone) ||
0645                                  (yLocal > deadZoneCenter[1] + cutZone && yLocal < deadZoneCenter[2] - cutZone) ||
0646                                  (yLocal > deadZoneCenter[2] + cutZone && yLocal < deadZoneCenter[3] - cutZone))) {
0647         pass = true;
0648       }
0649     }
0650   } else if (1 == station) {
0651     if (3 == ring) {
0652       deadZoneCenter[0] = -83.155;
0653       deadZoneCenter[1] = -22.7401;
0654       deadZoneCenter[2] = 27.86665;
0655       deadZoneCenter[3] = 81.005;
0656       if (yLocal > yBoundary && ((yLocal > deadZoneCenter[0] + cutZone && yLocal < deadZoneCenter[1] - cutZone) ||
0657                                  (yLocal > deadZoneCenter[1] + cutZone && yLocal < deadZoneCenter[2] - cutZone) ||
0658                                  (yLocal > deadZoneCenter[2] + cutZone && yLocal < deadZoneCenter[3] - cutZone))) {
0659         pass = true;
0660       }
0661     } else if (2 == ring) {
0662       deadZoneCenter[0] = -86.285;
0663       deadZoneCenter[1] = -32.88305;
0664       deadZoneCenter[2] = 32.867423;
0665       deadZoneCenter[3] = 88.205;
0666       if (yLocal > (yBoundary) && ((yLocal > deadZoneCenter[0] + cutZone && yLocal < deadZoneCenter[1] - cutZone) ||
0667                                    (yLocal > deadZoneCenter[1] + cutZone && yLocal < deadZoneCenter[2] - cutZone) ||
0668                                    (yLocal > deadZoneCenter[2] + cutZone && yLocal < deadZoneCenter[3] - cutZone))) {
0669         pass = true;
0670       }
0671     } else {
0672       deadZoneCenter[0] = -81.0;
0673       deadZoneCenter[1] = 81.0;
0674       if (yLocal > (yBoundary) && ((yLocal > deadZoneCenter[0] + cutZone && yLocal < deadZoneCenter[1] - cutZone))) {
0675         pass = true;
0676       }
0677     }
0678   }
0679   return pass;
0680 }
0681 
0682 void CSCEfficiency::fillDigiInfo(edm::Handle<CSCALCTDigiCollection> &alcts,
0683                                  edm::Handle<CSCCLCTDigiCollection> &clcts,
0684                                  edm::Handle<CSCCorrelatedLCTDigiCollection> &correlatedlcts,
0685                                  edm::Handle<CSCWireDigiCollection> &wires,
0686                                  edm::Handle<CSCStripDigiCollection> &strips,
0687                                  edm::Handle<edm::PSimHitContainer> &simhits,
0688                                  edm::Handle<CSCRecHit2DCollection> &rechits,
0689                                  edm::Handle<CSCSegmentCollection> &segments,
0690                                  edm::ESHandle<CSCGeometry> &cscGeom) {
0691   for (int iE = 0; iE < 2; iE++) {
0692     for (int iS = 0; iS < 4; iS++) {
0693       for (int iR = 0; iR < 4; iR++) {
0694         for (int iC = 0; iC < NumCh; iC++) {
0695           allSegments[iE][iS][iR][iC].clear();
0696           allCLCT[iE][iS][iR][iC] = allALCT[iE][iS][iR][iC] = allCorrLCT[iE][iS][iR][iC] = false;
0697           for (int iL = 0; iL < 6; iL++) {
0698             allStrips[iE][iS][iR][iC][iL].clear();
0699             allWG[iE][iS][iR][iC][iL].clear();
0700             allRechits[iE][iS][iR][iC][iL].clear();
0701             allSimhits[iE][iS][iR][iC][iL].clear();
0702           }
0703         }
0704       }
0705     }
0706   }
0707   //
0708   if (useDigis) {
0709     fillLCT_info(alcts, clcts, correlatedlcts);
0710     fillWG_info(wires, cscGeom);
0711     fillStrips_info(strips);
0712   }
0713   fillRechitsSegments_info(rechits, segments, cscGeom);
0714   if (!isData) {
0715     fillSimhit_info(simhits);
0716   }
0717 }
0718 
0719 void CSCEfficiency::fillLCT_info(edm::Handle<CSCALCTDigiCollection> &alcts,
0720                                  edm::Handle<CSCCLCTDigiCollection> &clcts,
0721                                  edm::Handle<CSCCorrelatedLCTDigiCollection> &correlatedlcts) {
0722   //---- ALCTDigis
0723   int nSize = 0;
0724   for (CSCALCTDigiCollection::DigiRangeIterator j = alcts->begin(); j != alcts->end(); j++) {
0725     ++nSize;
0726     const CSCDetId &id = (*j).first;
0727     const CSCALCTDigiCollection::Range &range = (*j).second;
0728     for (CSCALCTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
0729       // Valid digi in the chamber (or in neighbouring chamber)
0730       if ((*digiIt).isValid()) {
0731         allALCT[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh] = true;
0732       }
0733     }  // for digis in layer
0734   }  // end of for (j=...
0735   ALCTPerEvent->Fill(nSize);
0736   //---- CLCTDigis
0737   nSize = 0;
0738   for (CSCCLCTDigiCollection::DigiRangeIterator j = clcts->begin(); j != clcts->end(); j++) {
0739     ++nSize;
0740     const CSCDetId &id = (*j).first;
0741     std::vector<CSCCLCTDigi>::const_iterator digiIt = (*j).second.first;
0742     std::vector<CSCCLCTDigi>::const_iterator last = (*j).second.second;
0743     for (; digiIt != last; ++digiIt) {
0744       // Valid digi in the chamber (or in neighbouring chamber)
0745       if ((*digiIt).isValid()) {
0746         allCLCT[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh] = true;
0747       }
0748     }
0749   }
0750   CLCTPerEvent->Fill(nSize);
0751   //---- CorrLCTDigis
0752   for (CSCCorrelatedLCTDigiCollection::DigiRangeIterator j = correlatedlcts->begin(); j != correlatedlcts->end(); j++) {
0753     const CSCDetId &id = (*j).first;
0754     std::vector<CSCCorrelatedLCTDigi>::const_iterator digiIt = (*j).second.first;
0755     std::vector<CSCCorrelatedLCTDigi>::const_iterator last = (*j).second.second;
0756     for (; digiIt != last; ++digiIt) {
0757       // Valid digi in the chamber (or in neighbouring chamber)
0758       if ((*digiIt).isValid()) {
0759         allCorrLCT[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh] = true;
0760       }
0761     }
0762   }
0763 }
0764 //
0765 void CSCEfficiency::fillWG_info(edm::Handle<CSCWireDigiCollection> &wires, edm::ESHandle<CSCGeometry> &cscGeom) {
0766   //---- WIRE GROUPS
0767   for (CSCWireDigiCollection::DigiRangeIterator j = wires->begin(); j != wires->end(); j++) {
0768     CSCDetId id = (CSCDetId)(*j).first;
0769     const CSCLayer *layer_p = cscGeom->layer(id);
0770     const CSCLayerGeometry *layerGeom = layer_p->geometry();
0771     //
0772     std::vector<CSCWireDigi>::const_iterator digiItr = (*j).second.first;
0773     std::vector<CSCWireDigi>::const_iterator last = (*j).second.second;
0774     //
0775     for (; digiItr != last; ++digiItr) {
0776       std::pair<int, float> WG_pos(digiItr->getWireGroup(), layerGeom->yOfWireGroup(digiItr->getWireGroup()));
0777       std::pair<std::pair<int, float>, int> LayerSignal(WG_pos, digiItr->getTimeBin());
0778 
0779       //---- AllWG contains basic information about WG (WG number and Y-position, time bin)
0780       allWG[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh][id.layer() - 1].push_back(
0781           LayerSignal);
0782       if (printalot) {
0783         //std::cout<<" WG check : "<<std::endl;
0784         //printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),id.layer());
0785         //std::cout<<" WG size = "<<allWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh]
0786         //[id.layer()-1].size()<<std::endl;
0787       }
0788     }
0789   }
0790 }
0791 void CSCEfficiency::fillStrips_info(edm::Handle<CSCStripDigiCollection> &strips) {
0792   //---- STRIPS
0793   const float threshold = 13.3;
0794   for (CSCStripDigiCollection::DigiRangeIterator j = strips->begin(); j != strips->end(); j++) {
0795     CSCDetId id = (CSCDetId)(*j).first;
0796     int largestADCValue = -1;
0797     std::vector<CSCStripDigi>::const_iterator digiItr = (*j).second.first;
0798     std::vector<CSCStripDigi>::const_iterator last = (*j).second.second;
0799     for (; digiItr != last; ++digiItr) {
0800       int maxADC = largestADCValue;
0801       int myStrip = digiItr->getStrip();
0802       std::vector<int> myADCVals = digiItr->getADCCounts();
0803       float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]);
0804       float peakADC = -1000.;
0805       for (int myADCVal : myADCVals) {
0806         float diff = (float)myADCVal - thisPedestal;
0807         if (diff > threshold) {
0808           if (myADCVal > largestADCValue)
0809             largestADCValue = myADCVal;
0810           if (diff > peakADC)
0811             peakADC = diff;
0812         }
0813       }
0814       if (largestADCValue > maxADC) {  // FIX IT!!!
0815         std::pair<int, float> LayerSignal(myStrip, peakADC);
0816         //---- AllStrips contains basic information about strips
0817         //---- (strip number and peak signal for most significant strip in the layer)
0818         allStrips[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - 1][id.layer() - 1].clear();
0819         allStrips[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - 1][id.layer() - 1].push_back(
0820             LayerSignal);
0821       }
0822     }
0823   }
0824 }
0825 void CSCEfficiency::fillSimhit_info(edm::Handle<edm::PSimHitContainer> &simhits) {
0826   //---- SIMHITS
0827   edm::PSimHitContainer::const_iterator dSHsimIter;
0828   for (dSHsimIter = simhits->begin(); dSHsimIter != simhits->end(); dSHsimIter++) {
0829     // Get DetID for this simHit:
0830     CSCDetId sId = (CSCDetId)(*dSHsimIter).detUnitId();
0831     std::pair<LocalPoint, int> simHitPos((*dSHsimIter).localPosition(), (*dSHsimIter).particleType());
0832     allSimhits[sId.endcap() - 1][sId.station() - 1][sId.ring() - 1][sId.chamber() - FirstCh][sId.layer() - 1].push_back(
0833         simHitPos);
0834   }
0835 }
0836 //
0837 void CSCEfficiency::fillRechitsSegments_info(edm::Handle<CSCRecHit2DCollection> &rechits,
0838                                              edm::Handle<CSCSegmentCollection> &segments,
0839                                              edm::ESHandle<CSCGeometry> &cscGeom) {
0840   //---- RECHITS AND SEGMENTS
0841   //---- Loop over rechits
0842   if (printalot) {
0843     //printf("\tGet the recHits collection.\t ");
0844     printf("  The size of the rechit collection is %i\n", int(rechits->size()));
0845     //printf("\t...start loop over rechits...\n");
0846   }
0847   recHitsPerEvent->Fill(rechits->size());
0848   //---- Build iterator for rechits and loop :
0849   CSCRecHit2DCollection::const_iterator recIt;
0850   for (recIt = rechits->begin(); recIt != rechits->end(); recIt++) {
0851     //---- Find chamber with rechits in CSC
0852     CSCDetId id = (CSCDetId)(*recIt).cscDetId();
0853     if (printalot) {
0854       const CSCLayer *csclayer = cscGeom->layer(id);
0855       LocalPoint rhitlocal = (*recIt).localPosition();
0856       LocalError rerrlocal = (*recIt).localPositionError();
0857       GlobalPoint rhitglobal = csclayer->toGlobal(rhitlocal);
0858       printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
0859              id.endcap(),
0860              id.station(),
0861              id.ring(),
0862              id.chamber(),
0863              id.layer());
0864       printf("\t\tx,y,z: %f, %f, %f\texx,eey,exy: %f, %f, %f\tglobal x,y,z: %f, %f, %f \n",
0865              rhitlocal.x(),
0866              rhitlocal.y(),
0867              rhitlocal.z(),
0868              rerrlocal.xx(),
0869              rerrlocal.yy(),
0870              rerrlocal.xy(),
0871              rhitglobal.x(),
0872              rhitglobal.y(),
0873              rhitglobal.z());
0874     }
0875     std::pair<LocalPoint, bool> recHitPos((*recIt).localPosition(), false);
0876     allRechits[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh][id.layer() - 1].push_back(
0877         recHitPos);
0878   }
0879   //---- "Empty" chambers
0880   for (int iE = 0; iE < 2; iE++) {
0881     for (int iS = 0; iS < 4; iS++) {
0882       for (int iR = 0; iR < 4; iR++) {
0883         for (int iC = 0; iC < NumCh; iC++) {
0884           int numLayers = 0;
0885           for (int iL = 0; iL < 6; iL++) {
0886             if (!allRechits[iE][iS][iR][iC][iL].empty()) {
0887               ++numLayers;
0888             }
0889           }
0890           if (numLayers > 1) {
0891             emptyChambers[iE][iS][iR][iC] = false;
0892           } else {
0893             emptyChambers[iE][iS][iR][iC] = true;
0894           }
0895         }
0896       }
0897     }
0898   }
0899 
0900   //
0901   if (printalot) {
0902     printf("  The size of the segment collection is %i\n", int(segments->size()));
0903     //printf("\t...start loop over segments...\n");
0904   }
0905   segmentsPerEvent->Fill(segments->size());
0906   for (CSCSegmentCollection::const_iterator it = segments->begin(); it != segments->end(); it++) {
0907     CSCDetId id = (CSCDetId)(*it).cscDetId();
0908     StHist[id.endcap() - 1][id.station() - 1].segmentChi2_ndf->Fill((*it).chi2() / (*it).degreesOfFreedom());
0909     StHist[id.endcap() - 1][id.station() - 1].hitsInSegment->Fill((*it).nRecHits());
0910     if (printalot) {
0911       printf("\tendcap/station/ring/chamber: %i %i %i %i\n", id.endcap(), id.station(), id.ring(), id.chamber());
0912       std::cout << "\tposition(loc) = " << (*it).localPosition() << " error(loc) = " << (*it).localPositionError()
0913                 << std::endl;
0914       std::cout << "\t chi2/ndf = " << (*it).chi2() / (*it).degreesOfFreedom() << " nhits = " << (*it).nRecHits()
0915                 << std::endl;
0916     }
0917     allSegments[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh].push_back(
0918         make_pair((*it).localPosition(), (*it).localDirection()));
0919 
0920     //---- try to get the CSC recHits that contribute to this segment.
0921     //if (printalot) printf("\tGet the recHits for this segment.\t");
0922     std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits();
0923     int nRH = (*it).nRecHits();
0924     if (printalot) {
0925       printf("\tGet the recHits for this segment.\t");
0926       printf("    nRH = %i\n", nRH);
0927     }
0928     //---- Find which of the rechits in the chamber is in the segment
0929     int layerRH = 0;
0930     for (vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) {
0931       ++layerRH;
0932       CSCDetId idRH = (CSCDetId)(*iRH).cscDetId();
0933       if (printalot) {
0934         printf("\t%i RH\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
0935                layerRH,
0936                idRH.endcap(),
0937                idRH.station(),
0938                idRH.ring(),
0939                idRH.chamber(),
0940                idRH.layer());
0941       }
0942       for (size_t jRH = 0;
0943            jRH <
0944            allRechits[idRH.endcap() - 1][idRH.station() - 1][idRH.ring() - 1][idRH.chamber() - FirstCh][idRH.layer() - 1]
0945                .size();
0946            ++jRH) {
0947         float xDiff = iRH->localPosition().x() - allRechits[idRH.endcap() - 1][idRH.station() - 1][idRH.ring() - 1]
0948                                                            [idRH.chamber() - FirstCh][idRH.layer() - 1][jRH]
0949                                                                .first.x();
0950         float yDiff = iRH->localPosition().y() - allRechits[idRH.endcap() - 1][idRH.station() - 1][idRH.ring() - 1]
0951                                                            [idRH.chamber() - FirstCh][idRH.layer() - 1][jRH]
0952                                                                .first.y();
0953         if (fabs(xDiff) < 0.0001 && fabs(yDiff) < 0.0001) {
0954           std::pair<LocalPoint, bool> recHitPos(allRechits[idRH.endcap() - 1][idRH.station() - 1][idRH.ring() - 1]
0955                                                           [idRH.chamber() - FirstCh][idRH.layer() - 1][jRH]
0956                                                               .first,
0957                                                 true);
0958           allRechits[idRH.endcap() - 1][idRH.station() - 1][idRH.ring() - 1][idRH.chamber() - FirstCh][idRH.layer() - 1]
0959                     [jRH] = recHitPos;
0960           if (printalot) {
0961             std::cout << " number of the rechit (from zero) in the segment = " << jRH << std::endl;
0962           }
0963         }
0964       }
0965     }
0966   }
0967 }
0968 //
0969 
0970 void CSCEfficiency::ringCandidates(int station, float feta, std::map<std::string, bool> &chamberTypes) {
0971   // yeah, hardcoded again...
0972   switch (station) {
0973     case 1:
0974       if (feta > 0.85 && feta < 1.18) {  //ME13
0975         chamberTypes["ME13"] = true;
0976       }
0977       if (feta > 1.18 && feta < 1.7) {  //ME12
0978         chamberTypes["ME12"] = true;
0979       }
0980       if (feta > 1.5 && feta < 2.45) {  //ME11
0981         chamberTypes["ME11"] = true;
0982       }
0983       break;
0984     case 2:
0985       if (feta > 0.95 && feta < 1.6) {  //ME22
0986         chamberTypes["ME22"] = true;
0987       }
0988       if (feta > 1.55 && feta < 2.45) {  //ME21
0989         chamberTypes["ME21"] = true;
0990       }
0991       break;
0992     case 3:
0993       if (feta > 1.08 && feta < 1.72) {  //ME32
0994         chamberTypes["ME32"] = true;
0995       }
0996       if (feta > 1.69 && feta < 2.45) {  //ME31
0997         chamberTypes["ME31"] = true;
0998       }
0999       break;
1000     case 4:
1001       if (feta > 1.78 && feta < 2.45) {  //ME41
1002         chamberTypes["ME41"] = true;
1003       }
1004       break;
1005     default:
1006       break;
1007   }
1008 }
1009 //
1010 void CSCEfficiency::chamberCandidates(int station, int ring, float phi, std::vector<int> &coupleOfChambers) {
1011   coupleOfChambers.clear();
1012   // -pi< phi<+pi
1013   float phi_zero = 0.;  // check! the phi at the "edge" of Ch 1
1014   float phi_const = 2. * M_PI / 36.;
1015   int last_chamber = 36;
1016   int first_chamber = 1;
1017   if (1 != station && 1 == ring) {  // 18 chambers in the ring
1018     phi_const *= 2;
1019     last_chamber /= 2;
1020   }
1021   if (phi < 0.) {
1022     if (printalot)
1023       std::cout << " info: negative phi = " << phi << std::endl;
1024     phi += 2 * M_PI;
1025   }
1026   float chamber_float = (phi - phi_zero) / phi_const;
1027   int chamber_int = int(chamber_float);
1028   if (chamber_float - float(chamber_int) - 0.5 < 0.) {
1029     if (0 != chamber_int) {
1030       coupleOfChambers.push_back(chamber_int);
1031     } else {
1032       coupleOfChambers.push_back(last_chamber);
1033     }
1034     coupleOfChambers.push_back(chamber_int + 1);
1035 
1036   } else {
1037     coupleOfChambers.push_back(chamber_int + 1);
1038     if (last_chamber != chamber_int + 1) {
1039       coupleOfChambers.push_back(chamber_int + 2);
1040     } else {
1041       coupleOfChambers.push_back(first_chamber);
1042     }
1043   }
1044   if (printalot)
1045     std::cout << " phi = " << phi << " phi_zero = " << phi_zero << " phi_const = " << phi_const
1046               << " candidate chambers: first ch = " << coupleOfChambers[0] << " second ch = " << coupleOfChambers[1]
1047               << std::endl;
1048 }
1049 
1050 //
1051 bool CSCEfficiency::efficienciesPerChamber(CSCDetId &id,
1052                                            const CSCChamber *cscChamber,
1053                                            FreeTrajectoryState &ftsChamber) {
1054   int ec, st, rg, ch, secondRing;
1055   returnTypes(id, ec, st, rg, ch, secondRing);
1056 
1057   LocalVector localDir = cscChamber->toLocal(ftsChamber.momentum());
1058   if (printalot) {
1059     std::cout << " global dir = " << ftsChamber.momentum() << std::endl;
1060     std::cout << " local dir = " << localDir << std::endl;
1061     std::cout << " local theta = " << localDir.theta() << std::endl;
1062   }
1063   float dxdz = localDir.x() / localDir.z();
1064   float dydz = localDir.y() / localDir.z();
1065   if (2 == st || 3 == st) {
1066     if (printalot) {
1067       std::cout << "st 3 or 4 ... flip dy/dz" << std::endl;
1068     }
1069     dydz = -dydz;
1070   }
1071   if (printalot) {
1072     std::cout << "dy/dz = " << dydz << std::endl;
1073   }
1074   // Apply angle cut
1075   bool out = true;
1076   if (applyIPangleCuts) {
1077     if (dydz > local_DY_DZ_Max || dydz < local_DY_DZ_Min || fabs(dxdz) > local_DX_DZ_Max) {
1078       out = false;
1079     }
1080   }
1081 
1082   // Segments
1083   bool firstCondition = !allSegments[ec][st][rg][ch].empty() ? true : false;
1084   bool secondCondition = false;
1085   //---- ME1 is special as usual - ME1a and ME1b are actually one chamber
1086   if (secondRing > -1) {
1087     secondCondition = !allSegments[ec][st][secondRing][ch].empty() ? true : false;
1088   }
1089   if (firstCondition || secondCondition) {
1090     if (out) {
1091       ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(1);
1092     }
1093   } else {
1094     if (out) {
1095       ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(0);
1096     }
1097   }
1098 
1099   if (useDigis) {
1100     // ALCTs
1101     firstCondition = allALCT[ec][st][rg][ch];
1102     secondCondition = false;
1103     if (secondRing > -1) {
1104       secondCondition = allALCT[ec][st][secondRing][ch];
1105     }
1106     if (firstCondition || secondCondition) {
1107       if (out) {
1108         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(3);
1109       }
1110       // always apply partial angle cuts for this kind of histos
1111       if (fabs(dxdz) < local_DX_DZ_Max) {
1112         StHist[ec][st].EfficientALCT_momTheta->Fill(ftsChamber.momentum().theta());
1113         ChHist[ec][st][rg][ch].EfficientALCT_dydz->Fill(dydz);
1114       }
1115     } else {
1116       if (out) {
1117         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(2);
1118       }
1119       if (fabs(dxdz) < local_DX_DZ_Max) {
1120         StHist[ec][st].InefficientALCT_momTheta->Fill(ftsChamber.momentum().theta());
1121         ChHist[ec][st][rg][ch].InefficientALCT_dydz->Fill(dydz);
1122       }
1123       if (printalot) {
1124         std::cout << " missing ALCT (dy/dz = " << dydz << ")";
1125         printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n", ec + 1, st + 1, rg + 1, ch + 1);
1126       }
1127     }
1128 
1129     // CLCTs
1130     firstCondition = allCLCT[ec][st][rg][ch];
1131     secondCondition = false;
1132     if (secondRing > -1) {
1133       secondCondition = allCLCT[ec][st][secondRing][ch];
1134     }
1135     if (firstCondition || secondCondition) {
1136       if (out) {
1137         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(5);
1138       }
1139       if (dydz < local_DY_DZ_Max && dydz > local_DY_DZ_Min) {
1140         StHist[ec][st].EfficientCLCT_momPhi->Fill(ftsChamber.momentum().phi());  // - phi chamber...
1141         ChHist[ec][st][rg][ch].EfficientCLCT_dxdz->Fill(dxdz);
1142       }
1143     } else {
1144       if (out) {
1145         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(4);
1146       }
1147       if (dydz < local_DY_DZ_Max && dydz > local_DY_DZ_Min) {
1148         StHist[ec][st].InefficientCLCT_momPhi->Fill(ftsChamber.momentum().phi());  // - phi chamber...
1149         ChHist[ec][st][rg][ch].InefficientCLCT_dxdz->Fill(dxdz);
1150       }
1151       if (printalot) {
1152         std::cout << " missing CLCT  (dx/dz = " << dxdz << ")";
1153         printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n", ec + 1, st + 1, rg + 1, ch + 1);
1154       }
1155     }
1156     if (out) {
1157       // CorrLCTs
1158       firstCondition = allCorrLCT[ec][st][rg][ch];
1159       secondCondition = false;
1160       if (secondRing > -1) {
1161         secondCondition = allCorrLCT[ec][st][secondRing][ch];
1162       }
1163       if (firstCondition || secondCondition) {
1164         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(7);
1165       } else {
1166         ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(6);
1167       }
1168     }
1169   }
1170   return out;
1171 }
1172 
1173 //
1174 bool CSCEfficiency::stripWire_Efficiencies(CSCDetId &id, FreeTrajectoryState &ftsChamber) {
1175   int ec, st, rg, ch, secondRing;
1176   returnTypes(id, ec, st, rg, ch, secondRing);
1177 
1178   bool firstCondition, secondCondition;
1179   int missingLayers_s = 0;
1180   int missingLayers_wg = 0;
1181   for (int iLayer = 0; iLayer < 6; iLayer++) {
1182     //----Strips
1183     if (printalot) {
1184       printf("\t%i swEff: \tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
1185              iLayer + 1,
1186              id.endcap(),
1187              id.station(),
1188              id.ring(),
1189              id.chamber(),
1190              iLayer + 1);
1191       std::cout << " size S = "
1192                 << allStrips[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh][iLayer].size()
1193                 << "size W = "
1194                 << allWG[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh][iLayer].size()
1195                 << std::endl;
1196     }
1197     firstCondition = !allStrips[ec][st][rg][ch][iLayer].empty() ? true : false;
1198     //allSegments[ec][st][rg][ch].size() ? true : false;
1199     secondCondition = false;
1200     if (secondRing > -1) {
1201       secondCondition = !allStrips[ec][st][secondRing][ch][iLayer].empty() ? true : false;
1202     }
1203     if (firstCondition || secondCondition) {
1204       ChHist[ec][st][rg][ch].EfficientStrips->Fill(iLayer + 1);
1205     } else {
1206       if (printalot) {
1207         std::cout << "missing strips ";
1208         printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
1209                id.endcap(),
1210                id.station(),
1211                id.ring(),
1212                id.chamber(),
1213                iLayer + 1);
1214       }
1215     }
1216     // Wires
1217     firstCondition = !allWG[ec][st][rg][ch][iLayer].empty() ? true : false;
1218     secondCondition = false;
1219     if (secondRing > -1) {
1220       secondCondition = !allWG[ec][st][secondRing][ch][iLayer].empty() ? true : false;
1221     }
1222     if (firstCondition || secondCondition) {
1223       ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(iLayer + 1);
1224     } else {
1225       if (printalot) {
1226         std::cout << "missing wires ";
1227         printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
1228                id.endcap(),
1229                id.station(),
1230                id.ring(),
1231                id.chamber(),
1232                iLayer + 1);
1233       }
1234     }
1235   }
1236   // Normalization
1237   if (6 != missingLayers_s) {
1238     ChHist[ec][st][rg][ch].EfficientStrips->Fill(8);
1239   }
1240   if (6 != missingLayers_wg) {
1241     ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(8);
1242   }
1243   ChHist[ec][st][rg][ch].EfficientStrips->Fill(9);
1244   ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(9);
1245   //
1246   ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(1);
1247   if (missingLayers_s != missingLayers_wg) {
1248     ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(2);
1249     if (6 == missingLayers_wg) {
1250       ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(3);
1251       ChHist[ec][st][rg][ch].NoWires_momTheta->Fill(ftsChamber.momentum().theta());
1252     }
1253     if (6 == missingLayers_s) {
1254       ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(4);
1255       ChHist[ec][st][rg][ch].NoStrips_momPhi->Fill(ftsChamber.momentum().theta());
1256     }
1257   } else if (6 == missingLayers_s) {
1258     ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(5);
1259   }
1260 
1261   return true;
1262 }
1263 //
1264 bool CSCEfficiency::recSimHitEfficiency(CSCDetId &id, FreeTrajectoryState &ftsChamber) {
1265   int ec, st, rg, ch, secondRing;
1266   returnTypes(id, ec, st, rg, ch, secondRing);
1267   bool firstCondition, secondCondition;
1268   for (int iLayer = 0; iLayer < 6; iLayer++) {
1269     firstCondition = !allSimhits[ec][st][rg][ch][iLayer].empty() ? true : false;
1270     secondCondition = false;
1271     int thisRing = rg;
1272     if (secondRing > -1) {
1273       secondCondition = !allSimhits[ec][st][secondRing][ch][iLayer].empty() ? true : false;
1274       if (secondCondition) {
1275         thisRing = secondRing;
1276       }
1277     }
1278     if (firstCondition || secondCondition) {
1279       for (size_t iSH = 0; iSH < allSimhits[ec][st][thisRing][ch][iLayer].size(); ++iSH) {
1280         if (13 == fabs(allSimhits[ec][st][thisRing][ch][iLayer][iSH].second)) {
1281           ChHist[ec][st][rg][ch].SimSimhits->Fill(iLayer + 1);
1282           if (!allRechits[ec][st][thisRing][ch][iLayer].empty()) {
1283             ChHist[ec][st][rg][ch].SimRechits->Fill(iLayer + 1);
1284           }
1285           break;
1286         }
1287       }
1288       //---- Next is not too usefull...
1289       /*
1290       for(unsigned int iSimHits=0;
1291       iSimHits<allSimhits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size();
1292       iSimHits++){
1293     ChHist[ec][st][rg][id.chamber()-FirstCh].SimSimhits_each->Fill(iLayer+1);
1294       }
1295       for(unsigned int iRecHits=0;
1296       iRecHits<allRechits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size();
1297       iRecHits++){
1298     ChHist[ec][st][rg][id.chamber()-FirstCh].SimRechits_each->Fill(iLayer+1);
1299       }
1300       */
1301       //
1302     }
1303   }
1304   return true;
1305 }
1306 
1307 //
1308 bool CSCEfficiency::recHitSegment_Efficiencies(CSCDetId &id,
1309                                                const CSCChamber *cscChamber,
1310                                                FreeTrajectoryState &ftsChamber) {
1311   int ec, st, rg, ch, secondRing;
1312   returnTypes(id, ec, st, rg, ch, secondRing);
1313   bool firstCondition, secondCondition;
1314 
1315   std::vector<bool> missingLayers_rh(6);
1316   std::vector<int> usedInSegment(6);
1317   // Rechits
1318   if (printalot)
1319     std::cout << "RecHits eff" << std::endl;
1320   for (int iLayer = 0; iLayer < 6; ++iLayer) {
1321     firstCondition = !allRechits[ec][st][rg][ch][iLayer].empty() ? true : false;
1322     secondCondition = false;
1323     int thisRing = rg;
1324     if (secondRing > -1) {
1325       secondCondition = !allRechits[ec][st][secondRing][ch][iLayer].empty() ? true : false;
1326       if (secondCondition) {
1327         thisRing = secondRing;
1328       }
1329     }
1330     if (firstCondition || secondCondition) {
1331       ChHist[ec][st][rg][ch].EfficientRechits_good->Fill(iLayer + 1);
1332       for (size_t iR = 0; iR < allRechits[ec][st][thisRing][ch][iLayer].size(); ++iR) {
1333         if (allRechits[ec][st][thisRing][ch][iLayer][iR].second) {
1334           usedInSegment[iLayer] = 1;
1335           break;
1336         } else {
1337           usedInSegment[iLayer] = -1;
1338         }
1339       }
1340     } else {
1341       missingLayers_rh[iLayer] = true;
1342       if (printalot) {
1343         std::cout << "missing rechits ";
1344         printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",
1345                id.endcap(),
1346                id.station(),
1347                id.ring(),
1348                id.chamber(),
1349                iLayer + 1);
1350       }
1351     }
1352   }
1353   GlobalVector globalDir;
1354   GlobalPoint globalPos;
1355   // Segments
1356   firstCondition = !allSegments[ec][st][rg][ch].empty() ? true : false;
1357   secondCondition = false;
1358   int secondSize = 0;
1359   int thisRing = rg;
1360   if (secondRing > -1) {
1361     secondCondition = !allSegments[ec][st][secondRing][ch].empty() ? true : false;
1362     secondSize = allSegments[ec][st][secondRing][ch].size();
1363     if (secondCondition) {
1364       thisRing = secondRing;
1365     }
1366   }
1367   if (firstCondition || secondCondition) {
1368     if (printalot)
1369       std::cout << "segments - start ec = " << ec << " st = " << st << " rg = " << rg << " ch = " << ch << std::endl;
1370     StHist[ec][st].EfficientSegments_XY->Fill(ftsChamber.position().x(), ftsChamber.position().y());
1371     if (1 == allSegments[ec][st][rg][ch].size() + secondSize) {
1372       globalDir = cscChamber->toGlobal(allSegments[ec][st][thisRing][ch][0].second);
1373       globalPos = cscChamber->toGlobal(allSegments[ec][st][thisRing][ch][0].first);
1374       StHist[ec][st].EfficientSegments_eta->Fill(fabs(ftsChamber.position().eta()));
1375       double residual =
1376           sqrt(pow(ftsChamber.position().x() - globalPos.x(), 2) + pow(ftsChamber.position().y() - globalPos.y(), 2) +
1377                pow(ftsChamber.position().z() - globalPos.z(), 2));
1378       if (printalot)
1379         std::cout << " fts.position() = " << ftsChamber.position() << " segPos = " << globalPos << " res = " << residual
1380                   << std::endl;
1381       StHist[ec][st].ResidualSegments->Fill(residual);
1382     }
1383     for (int iLayer = 0; iLayer < 6; ++iLayer) {
1384       if (printalot)
1385         std::cout << " iLayer = " << iLayer << " usedInSegment = " << usedInSegment[iLayer] << std::endl;
1386       if (0 != usedInSegment[iLayer]) {
1387         if (-1 == usedInSegment[iLayer]) {
1388           ChHist[ec][st][rg][ch].InefficientSingleHits->Fill(iLayer + 1);
1389         }
1390         ChHist[ec][st][rg][ch].AllSingleHits->Fill(iLayer + 1);
1391       }
1392       firstCondition = !allRechits[ec][st][rg][ch][iLayer].empty() ? true : false;
1393       secondCondition = false;
1394       if (secondRing > -1) {
1395         secondCondition = !allRechits[ec][st][secondRing][ch][iLayer].empty() ? true : false;
1396       }
1397       float stripAngle = 99999.;
1398       std::vector<float> posXY(2);
1399       bool oneSegment = false;
1400       if (1 == allSegments[ec][st][rg][ch].size() + secondSize) {
1401         oneSegment = true;
1402         const BoundPlane bp = cscChamber->layer(iLayer + 1)->surface();
1403         linearExtrapolation(globalPos, globalDir, bp.position().z(), posXY);
1404         GlobalPoint gp_extrapol(posXY.at(0), posXY.at(1), bp.position().z());
1405         const LocalPoint lp_extrapol = cscChamber->layer(iLayer + 1)->toLocal(gp_extrapol);
1406         posXY.at(0) = lp_extrapol.x();
1407         posXY.at(1) = lp_extrapol.y();
1408         int nearestStrip = cscChamber->layer(iLayer + 1)->geometry()->nearestStrip(lp_extrapol);
1409         stripAngle = cscChamber->layer(iLayer + 1)->geometry()->stripAngle(nearestStrip) - M_PI / 2.;
1410       }
1411       if (firstCondition || secondCondition) {
1412         ChHist[ec][st][rg][ch].EfficientRechits_inSegment->Fill(iLayer + 1);
1413         if (oneSegment) {
1414           ChHist[ec][st][rg][ch].Y_EfficientRecHits_inSegment[iLayer]->Fill(posXY.at(1));
1415           ChHist[ec][st][rg][ch].Phi_EfficientRecHits_inSegment[iLayer]->Fill(stripAngle);
1416         }
1417       } else {
1418         if (oneSegment) {
1419           ChHist[ec][st][rg][ch].Y_InefficientRecHits_inSegment[iLayer]->Fill(posXY.at(1));
1420           ChHist[ec][st][rg][ch].Phi_InefficientRecHits_inSegment[iLayer]->Fill(stripAngle);
1421         }
1422       }
1423     }
1424   } else {
1425     StHist[ec][st].InefficientSegments_XY->Fill(ftsChamber.position().x(), ftsChamber.position().y());
1426     if (printalot) {
1427       std::cout << "missing segment " << std::endl;
1428       printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n", id.endcap(), id.station(), id.ring(), id.chamber());
1429       std::cout << " fts.position() = " << ftsChamber.position() << std::endl;
1430     }
1431   }
1432   // Normalization
1433   ChHist[ec][st][rg][ch].EfficientRechits_good->Fill(8);
1434   if (allSegments[ec][st][rg][ch].size() + secondSize < 2) {
1435     StHist[ec][st].AllSegments_eta->Fill(fabs(ftsChamber.position().eta()));
1436   }
1437   ChHist[ec][st][rg][id.chamber() - FirstCh].EfficientRechits_inSegment->Fill(9);
1438 
1439   return true;
1440 }
1441 //
1442 void CSCEfficiency::returnTypes(CSCDetId &id, int &ec, int &st, int &rg, int &ch, int &secondRing) {
1443   ec = id.endcap() - 1;
1444   st = id.station() - 1;
1445   rg = id.ring() - 1;
1446   secondRing = -1;
1447   if (1 == id.station() && (4 == id.ring() || 1 == id.ring())) {
1448     rg = 0;
1449     secondRing = 3;
1450   }
1451   ch = id.chamber() - FirstCh;
1452 }
1453 
1454 //
1455 void CSCEfficiency::getFromFTS(const FreeTrajectoryState &fts,
1456                                CLHEP::Hep3Vector &p3,
1457                                CLHEP::Hep3Vector &r3,
1458                                int &charge,
1459                                AlgebraicSymMatrix66 &cov) {
1460   GlobalVector p3GV = fts.momentum();
1461   GlobalPoint r3GP = fts.position();
1462 
1463   p3.set(p3GV.x(), p3GV.y(), p3GV.z());
1464   r3.set(r3GP.x(), r3GP.y(), r3GP.z());
1465 
1466   charge = fts.charge();
1467   cov = fts.hasError() ? fts.cartesianError().matrix() : AlgebraicSymMatrix66();
1468 }
1469 
1470 FreeTrajectoryState CSCEfficiency::getFromCLHEP(const CLHEP::Hep3Vector &p3,
1471                                                 const CLHEP::Hep3Vector &r3,
1472                                                 int charge,
1473                                                 const AlgebraicSymMatrix66 &cov,
1474                                                 const MagneticField *field) {
1475   GlobalVector p3GV(p3.x(), p3.y(), p3.z());
1476   GlobalPoint r3GP(r3.x(), r3.y(), r3.z());
1477   GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, field);
1478 
1479   CartesianTrajectoryError tCov(cov);
1480 
1481   return cov.kRows == 6 ? FreeTrajectoryState(tPars, tCov) : FreeTrajectoryState(tPars);
1482 }
1483 
1484 void CSCEfficiency::linearExtrapolation(GlobalPoint initialPosition,
1485                                         GlobalVector initialDirection,
1486                                         float zSurface,
1487                                         std::vector<float> &posZY) {
1488   double paramLine = lineParameter(initialPosition.z(), zSurface, initialDirection.z());
1489   double xPosition = extrapolate1D(initialPosition.x(), initialDirection.x(), paramLine);
1490   double yPosition = extrapolate1D(initialPosition.y(), initialDirection.y(), paramLine);
1491   posZY.clear();
1492   posZY.push_back(xPosition);
1493   posZY.push_back(yPosition);
1494 }
1495 //
1496 double CSCEfficiency::extrapolate1D(double initPosition, double initDirection, double parameterOfTheLine) {
1497   double extrapolatedPosition = initPosition + initDirection * parameterOfTheLine;
1498   return extrapolatedPosition;
1499 }
1500 //
1501 double CSCEfficiency::lineParameter(double initZPosition, double destZPosition, double initZDirection) {
1502   double paramLine = (destZPosition - initZPosition) / initZDirection;
1503   return paramLine;
1504 }
1505 //
1506 void CSCEfficiency::chooseDirection(CLHEP::Hep3Vector &innerPosition, CLHEP::Hep3Vector &outerPosition) {
1507   //---- Be careful with trigger conditions too
1508   if (!isIPdata) {
1509     float dy = outerPosition.y() - innerPosition.y();
1510     float dz = outerPosition.z() - innerPosition.z();
1511     if (isBeamdata) {
1512       if (dz > 0) {
1513         alongZ = true;
1514       } else {
1515         alongZ = false;
1516       }
1517     } else {  //cosmics
1518       if (dy / dz > 0) {
1519         alongZ = false;
1520       } else {
1521         alongZ = true;
1522       }
1523     }
1524   }
1525 }
1526 //
1527 const Propagator *CSCEfficiency::propagator(std::string propagatorName) const {
1528   return &*theService->propagator(propagatorName);
1529 }
1530 
1531 //
1532 TrajectoryStateOnSurface CSCEfficiency::propagate(FreeTrajectoryState &ftsStart, const BoundPlane &bpDest) {
1533   TrajectoryStateOnSurface tSOSDest;
1534   std::string propagatorName;
1535   /*
1536 // it would work if cosmic muons had properly assigned direction...
1537   bool dzPositive = bpDest.position().z() - ftsStart.position().z() > 0 ? true : false;
1538  //---- Be careful with trigger conditions too
1539   if(!isIPdata){
1540     bool rightDirection = !(alongZ^dzPositive);
1541     if(rightDirection){
1542       if(printalot) std::cout<<" propagate along momentum"<<std::endl;
1543       propagatorName = "SteppingHelixPropagatorAlong";
1544     }
1545     else{
1546       if(printalot) std::cout<<" propagate opposite momentum"<<std::endl;
1547       propagatorName = "SteppingHelixPropagatorOpposite";
1548     }
1549   }
1550   else{
1551     if(printalot) std::cout<<" propagate any (momentum)"<<std::endl;
1552     propagatorName = "SteppingHelixPropagatorAny";
1553   }
1554 */
1555   propagatorName = "SteppingHelixPropagatorAny";
1556   tSOSDest = propagator(propagatorName)->propagate(ftsStart, bpDest);
1557   return tSOSDest;
1558 }
1559 //
1560 bool CSCEfficiency::applyTrigger(edm::Handle<edm::TriggerResults> &hltR, const edm::TriggerNames &triggerNames) {
1561   bool triggerPassed = true;
1562   std::vector<std::string> hlNames = triggerNames.triggerNames();
1563   pointToTriggers.clear();
1564   for (size_t imyT = 0; imyT < myTriggers.size(); ++imyT) {
1565     for (size_t iT = 0; iT < hlNames.size(); ++iT) {
1566       //std::cout<<" iT = "<<iT<<" hlNames[iT] = "<<hlNames[iT]<<
1567       //" : wasrun = "<<hltR->wasrun(iT)<<" accept = "<<
1568       //     hltR->accept(iT)<<" !error = "<<
1569       //    !hltR->error(iT)<<std::endl;
1570       if (!imyT) {
1571         if (hltR->wasrun(iT) && hltR->accept(iT) && !hltR->error(iT)) {
1572           TriggersFired->Fill(iT);
1573         }
1574       }
1575       if (hlNames[iT] == myTriggers[imyT]) {
1576         pointToTriggers.push_back(iT);
1577         if (imyT) {
1578           break;
1579         }
1580       }
1581     }
1582   }
1583   if (pointToTriggers.size() != myTriggers.size()) {
1584     pointToTriggers.clear();
1585     if (printalot) {
1586       std::cout << " Not all trigger names found - all trigger specifications will be ignored. Check your cfg file!"
1587                 << std::endl;
1588     }
1589   } else {
1590     if (!pointToTriggers.empty()) {
1591       if (printalot) {
1592         std::cout << "The following triggers will be required in the event: " << std::endl;
1593         for (size_t imyT = 0; imyT < pointToTriggers.size(); ++imyT) {
1594           std::cout << "  " << hlNames[pointToTriggers[imyT]];
1595         }
1596         std::cout << std::endl;
1597         std::cout << " in condition (AND/OR) : " << !andOr << "/" << andOr << std::endl;
1598       }
1599     }
1600   }
1601 
1602   if (hltR.isValid()) {
1603     if (pointToTriggers.empty()) {
1604       if (printalot) {
1605         std::cout
1606             << " No triggers specified in the configuration or all ignored - no trigger information will be considered"
1607             << std::endl;
1608       }
1609     }
1610     for (size_t imyT = 0; imyT < pointToTriggers.size(); ++imyT) {
1611       if (hltR->wasrun(pointToTriggers[imyT]) && hltR->accept(pointToTriggers[imyT]) &&
1612           !hltR->error(pointToTriggers[imyT])) {
1613         triggerPassed = true;
1614         if (andOr) {
1615           break;
1616         }
1617       } else {
1618         triggerPassed = false;
1619         if (!andOr) {
1620           triggerPassed = false;
1621           break;
1622         }
1623       }
1624     }
1625   } else {
1626     if (printalot) {
1627       std::cout << " TriggerResults handle returns invalid state?! No trigger information will be considered"
1628                 << std::endl;
1629     }
1630   }
1631   if (printalot) {
1632     std::cout << " Trigger passed: " << triggerPassed << std::endl;
1633   }
1634   return triggerPassed;
1635 }
1636 //
1637 
1638 // Constructor
1639 CSCEfficiency::CSCEfficiency(const edm::ParameterSet &pset) {
1640   // const float Xmin = -70;
1641   //const float Xmax = 70;
1642   //const int nXbins = int(4.*(Xmax - Xmin));
1643   const float Ymin = -165;
1644   const float Ymax = 165;
1645   const int nYbins = int((Ymax - Ymin) / 2);
1646   const float Layer_min = -0.5;
1647   const float Layer_max = 9.5;
1648   const int nLayer_bins = int(Layer_max - Layer_min);
1649   //
1650 
1651   //---- Get the input parameters
1652   printout_NEvents = pset.getUntrackedParameter<unsigned int>("printout_NEvents", 0);
1653   rootFileName = pset.getUntrackedParameter<string>("rootFileName", "cscHists.root");
1654 
1655   isData = pset.getUntrackedParameter<bool>("runOnData", true);                             //
1656   isIPdata = pset.getUntrackedParameter<bool>("IPdata", false);                             //
1657   isBeamdata = pset.getUntrackedParameter<bool>("Beamdata", false);                         //
1658   getAbsoluteEfficiency = pset.getUntrackedParameter<bool>("getAbsoluteEfficiency", true);  //
1659   useDigis = pset.getUntrackedParameter<bool>("useDigis", true);                            //
1660   distanceFromDeadZone = pset.getUntrackedParameter<double>("distanceFromDeadZone", 10.);   //
1661   minP = pset.getUntrackedParameter<double>("minP", 20.);                                   //
1662   maxP = pset.getUntrackedParameter<double>("maxP", 100.);                                  //
1663   maxNormChi2 = pset.getUntrackedParameter<double>("maxNormChi2", 3.);                      //
1664   minTrackHits = pset.getUntrackedParameter<unsigned int>("minTrackHits", 10);              //
1665 
1666   applyIPangleCuts = pset.getUntrackedParameter<bool>("applyIPangleCuts", false);  //
1667   local_DY_DZ_Max = pset.getUntrackedParameter<double>("local_DY_DZ_Max", -0.1);   //
1668   local_DY_DZ_Min = pset.getUntrackedParameter<double>("local_DY_DZ_Min", -0.8);   //
1669   local_DX_DZ_Max = pset.getUntrackedParameter<double>("local_DX_DZ_Max", 0.2);    //
1670 
1671   sd_token = consumes<CSCStripDigiCollection>(pset.getParameter<edm::InputTag>("stripDigiTag"));
1672   wd_token = consumes<CSCWireDigiCollection>(pset.getParameter<edm::InputTag>("wireDigiTag"));
1673   al_token = consumes<CSCALCTDigiCollection>(pset.getParameter<edm::InputTag>("alctDigiTag"));
1674   cl_token = consumes<CSCCLCTDigiCollection>(pset.getParameter<edm::InputTag>("clctDigiTag"));
1675   co_token = consumes<CSCCorrelatedLCTDigiCollection>(pset.getParameter<edm::InputTag>("corrlctDigiTag"));
1676   rh_token = consumes<CSCRecHit2DCollection>(pset.getParameter<edm::InputTag>("rechitTag"));
1677   se_token = consumes<CSCSegmentCollection>(pset.getParameter<edm::InputTag>("segmentTag"));
1678   tk_token = consumes<edm::View<reco::Track> >(pset.getParameter<edm::InputTag>("tracksTag"));
1679   sh_token = consumes<edm::PSimHitContainer>(pset.getParameter<edm::InputTag>("simHitTag"));
1680 
1681   geomToken_ = esConsumes<CSCGeometry, MuonGeometryRecord>();
1682 
1683   edm::ParameterSet serviceParameters = pset.getParameter<edm::ParameterSet>("ServiceParameters");
1684   // maybe use the service for getting magnetic field, propagators, etc. ...
1685   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
1686 
1687   // Trigger
1688   useTrigger = pset.getUntrackedParameter<bool>("useTrigger", false);
1689 
1690   ht_token = consumes<edm::TriggerResults>(pset.getParameter<edm::InputTag>("HLTriggerResults"));
1691 
1692   myTriggers = pset.getParameter<std::vector<std::string> >("myTriggers");
1693   andOr = pset.getUntrackedParameter<bool>("andOr");
1694   pointToTriggers.clear();
1695 
1696   //---- set counter to zero
1697   nEventsAnalyzed = 0;
1698   //---- set presence of magnetic field
1699   magField = true;
1700   //
1701   std::string Path = "AllChambers/";
1702   std::string FullName;
1703   //---- File with output histograms
1704   theFile = new TFile(rootFileName.c_str(), "RECREATE");
1705   theFile->cd();
1706   //---- Book histograms for the analysis
1707   char SpecName[60];
1708 
1709   sprintf(SpecName, "DataFlow");
1710   DataFlow = new TH1F(SpecName, "Data flow;condition number;entries", 40, -0.5, 39.5);
1711   //
1712   sprintf(SpecName, "TriggersFired");
1713   TriggersFired = new TH1F(SpecName, "Triggers fired;trigger number;entries", 140, -0.5, 139.5);
1714   //
1715   int Chan = 50;
1716   float minChan = -0.5;
1717   float maxChan = 49.5;
1718   //
1719   sprintf(SpecName, "ALCTPerEvent");
1720   ALCTPerEvent = new TH1F(SpecName, "ALCTs per event;N digis;entries", Chan, minChan, maxChan);
1721   //
1722   sprintf(SpecName, "CLCTPerEvent");
1723   CLCTPerEvent = new TH1F(SpecName, "CLCTs per event;N digis;entries", Chan, minChan, maxChan);
1724   //
1725   sprintf(SpecName, "recHitsPerEvent");
1726   recHitsPerEvent = new TH1F(SpecName, "RecHits per event;N digis;entries", 150, -0.5, 149.5);
1727   //
1728   sprintf(SpecName, "segmentsPerEvent");
1729   segmentsPerEvent = new TH1F(SpecName, "segments per event;N digis;entries", Chan, minChan, maxChan);
1730   //
1731   //---- Book groups of histograms (for any chamber)
1732 
1733   map<std::string, bool>::iterator iter;
1734   for (int ec = 0; ec < 2; ++ec) {
1735     for (int st = 0; st < 4; ++st) {
1736       theFile->cd();
1737       sprintf(SpecName, "Stations__E%d_S%d", ec + 1, st + 1);
1738       theFile->mkdir(SpecName);
1739       theFile->cd(SpecName);
1740 
1741       //
1742       sprintf(SpecName, "segmentChi2_ndf_St%d", st + 1);
1743       StHist[ec][st].segmentChi2_ndf = new TH1F(SpecName, "Chi2/ndf of a segment;chi2/ndf;entries", 100, 0., 20.);
1744       //
1745       sprintf(SpecName, "hitsInSegment_St%d", st + 1);
1746       StHist[ec][st].hitsInSegment = new TH1F(SpecName, "Number of hits in a segment;nHits;entries", 7, -0.5, 6.5);
1747       //
1748       Chan = 170;
1749       minChan = 0.85;
1750       maxChan = 2.55;
1751       //
1752       sprintf(SpecName, "AllSegments_eta_St%d", st + 1);
1753       StHist[ec][st].AllSegments_eta = new TH1F(SpecName, "All segments in eta;eta;entries", Chan, minChan, maxChan);
1754       //
1755       sprintf(SpecName, "EfficientSegments_eta_St%d", st + 1);
1756       StHist[ec][st].EfficientSegments_eta =
1757           new TH1F(SpecName, "Efficient segments in eta;eta;entries", Chan, minChan, maxChan);
1758       //
1759       sprintf(SpecName, "ResidualSegments_St%d", st + 1);
1760       StHist[ec][st].ResidualSegments = new TH1F(SpecName, "Residual (segments);residual,cm;entries", 75, 0., 15.);
1761       //
1762       Chan = 200;
1763       minChan = -800.;
1764       maxChan = 800.;
1765       int Chan2 = 200;
1766       float minChan2 = -800.;
1767       float maxChan2 = 800.;
1768 
1769       sprintf(SpecName, "EfficientSegments_XY_St%d", st + 1);
1770       StHist[ec][st].EfficientSegments_XY =
1771           new TH2F(SpecName, "Efficient segments in XY;X;Y", Chan, minChan, maxChan, Chan2, minChan2, maxChan2);
1772       sprintf(SpecName, "InefficientSegments_XY_St%d", st + 1);
1773       StHist[ec][st].InefficientSegments_XY =
1774           new TH2F(SpecName, "Inefficient segments in XY;X;Y", Chan, minChan, maxChan, Chan2, minChan2, maxChan2);
1775       //
1776       Chan = 80;
1777       minChan = 0;
1778       maxChan = 3.2;
1779       sprintf(SpecName, "EfficientALCT_momTheta_St%d", st + 1);
1780       StHist[ec][st].EfficientALCT_momTheta =
1781           new TH1F(SpecName, "Efficient ALCT in theta (momentum);theta, rad;entries", Chan, minChan, maxChan);
1782       //
1783       sprintf(SpecName, "InefficientALCT_momTheta_St%d", st + 1);
1784       StHist[ec][st].InefficientALCT_momTheta =
1785           new TH1F(SpecName, "Inefficient ALCT in theta (momentum);theta, rad;entries", Chan, minChan, maxChan);
1786       //
1787       Chan = 160;
1788       minChan = -3.2;
1789       maxChan = 3.2;
1790       sprintf(SpecName, "EfficientCLCT_momPhi_St%d", st + 1);
1791       StHist[ec][st].EfficientCLCT_momPhi =
1792           new TH1F(SpecName, "Efficient CLCT in phi (momentum);phi, rad;entries", Chan, minChan, maxChan);
1793       //
1794       sprintf(SpecName, "InefficientCLCT_momPhi_St%d", st + 1);
1795       StHist[ec][st].InefficientCLCT_momPhi =
1796           new TH1F(SpecName, "Inefficient CLCT in phi (momentum);phi, rad;entries", Chan, minChan, maxChan);
1797       //
1798       theFile->cd();
1799       for (int rg = 0; rg < 3; ++rg) {
1800         if (0 != st && rg > 1) {
1801           continue;
1802         } else if (1 == rg && 3 == st) {
1803           continue;
1804         }
1805         for (int iChamber = FirstCh; iChamber < FirstCh + NumCh; iChamber++) {
1806           if (0 != st && 0 == rg && iChamber > 18) {
1807             continue;
1808           }
1809           theFile->cd();
1810           sprintf(SpecName, "Chambers__E%d_S%d_R%d_Chamber_%d", ec + 1, st + 1, rg + 1, iChamber);
1811           theFile->mkdir(SpecName);
1812           theFile->cd(SpecName);
1813           //
1814 
1815           sprintf(SpecName, "EfficientRechits_inSegment_Ch%d", iChamber);
1816           ChHist[ec][st][rg][iChamber - FirstCh].EfficientRechits_inSegment = new TH1F(
1817               SpecName, "Existing RecHit given a segment;layers (1-6);entries", nLayer_bins, Layer_min, Layer_max);
1818           //
1819           sprintf(SpecName, "InefficientSingleHits_Ch%d", iChamber);
1820           ChHist[ec][st][rg][iChamber - FirstCh].InefficientSingleHits = new TH1F(
1821               SpecName, "Single RecHits not in the segment;layers (1-6);entries ", nLayer_bins, Layer_min, Layer_max);
1822           //
1823           sprintf(SpecName, "AllSingleHits_Ch%d", iChamber);
1824           ChHist[ec][st][rg][iChamber - FirstCh].AllSingleHits = new TH1F(
1825               SpecName, "Single RecHits given a segment; layers (1-6);entries", nLayer_bins, Layer_min, Layer_max);
1826           //
1827           sprintf(SpecName, "digiAppearanceCount_Ch%d", iChamber);
1828           ChHist[ec][st][rg][iChamber - FirstCh].digiAppearanceCount =
1829               new TH1F(SpecName,
1830                        "Digi appearance (no-yes): segment(0,1), ALCT(2,3), CLCT(4,5), CorrLCT(6,7); digi type;entries",
1831                        8,
1832                        -0.5,
1833                        7.5);
1834           //
1835           Chan = 100;
1836           minChan = -1.1;
1837           maxChan = 0.9;
1838           sprintf(SpecName, "EfficientALCT_dydz_Ch%d", iChamber);
1839           ChHist[ec][st][rg][iChamber - FirstCh].EfficientALCT_dydz =
1840               new TH1F(SpecName, "Efficient ALCT; local dy/dz (ME 3 and 4 flipped);entries", Chan, minChan, maxChan);
1841           //
1842           sprintf(SpecName, "InefficientALCT_dydz_Ch%d", iChamber);
1843           ChHist[ec][st][rg][iChamber - FirstCh].InefficientALCT_dydz =
1844               new TH1F(SpecName, "Inefficient ALCT; local dy/dz (ME 3 and 4 flipped);entries", Chan, minChan, maxChan);
1845           //
1846           Chan = 100;
1847           minChan = -1.;
1848           maxChan = 1.0;
1849           sprintf(SpecName, "EfficientCLCT_dxdz_Ch%d", iChamber);
1850           ChHist[ec][st][rg][iChamber - FirstCh].EfficientCLCT_dxdz =
1851               new TH1F(SpecName, "Efficient CLCT; local dxdz;entries", Chan, minChan, maxChan);
1852           //
1853           sprintf(SpecName, "InefficientCLCT_dxdz_Ch%d", iChamber);
1854           ChHist[ec][st][rg][iChamber - FirstCh].InefficientCLCT_dxdz =
1855               new TH1F(SpecName, "Inefficient CLCT; local dxdz;entries", Chan, minChan, maxChan);
1856           //
1857           sprintf(SpecName, "EfficientRechits_good_Ch%d", iChamber);
1858           ChHist[ec][st][rg][iChamber - FirstCh].EfficientRechits_good = new TH1F(
1859               SpecName, "Existing RecHit - sensitive area only;layers (1-6);entries", nLayer_bins, Layer_min, Layer_max);
1860           //
1861           sprintf(SpecName, "EfficientStrips_Ch%d", iChamber);
1862           ChHist[ec][st][rg][iChamber - FirstCh].EfficientStrips =
1863               new TH1F(SpecName, "Existing strip;layer (1-6); entries", nLayer_bins, Layer_min, Layer_max);
1864           //
1865           sprintf(SpecName, "EfficientWireGroups_Ch%d", iChamber);
1866           ChHist[ec][st][rg][iChamber - FirstCh].EfficientWireGroups =
1867               new TH1F(SpecName, "Existing WireGroups;layer (1-6); entries ", nLayer_bins, Layer_min, Layer_max);
1868           //
1869           sprintf(SpecName, "StripWiresCorrelations_Ch%d", iChamber);
1870           ChHist[ec][st][rg][iChamber - FirstCh].StripWiresCorrelations =
1871               new TH1F(SpecName, "StripWire correlations;; entries ", 5, 0.5, 5.5);
1872           //
1873           Chan = 80;
1874           minChan = 0;
1875           maxChan = 3.2;
1876           sprintf(SpecName, "NoWires_momTheta_Ch%d", iChamber);
1877           ChHist[ec][st][rg][iChamber - FirstCh].NoWires_momTheta =
1878               new TH1F(SpecName,
1879                        "No wires (all strips present) - in theta (momentum);theta, rad;entries",
1880                        Chan,
1881                        minChan,
1882                        maxChan);
1883           //
1884           Chan = 160;
1885           minChan = -3.2;
1886           maxChan = 3.2;
1887           sprintf(SpecName, "NoStrips_momPhi_Ch%d", iChamber);
1888           ChHist[ec][st][rg][iChamber - FirstCh].NoStrips_momPhi = new TH1F(
1889               SpecName, "No strips (all wires present) - in phi (momentum);phi, rad;entries", Chan, minChan, maxChan);
1890           //
1891           for (int iLayer = 0; iLayer < 6; iLayer++) {
1892             sprintf(SpecName, "Y_InefficientRecHits_inSegment_Ch%d_L%d", iChamber, iLayer);
1893             ChHist[ec][st][rg][iChamber - FirstCh].Y_InefficientRecHits_inSegment.push_back(
1894                 new TH1F(SpecName,
1895                          "Missing RecHit/layer in a segment (local system, whole chamber);Y, cm; entries",
1896                          nYbins,
1897                          Ymin,
1898                          Ymax));
1899             //
1900             sprintf(SpecName, "Y_EfficientRecHits_inSegment_Ch%d_L%d", iChamber, iLayer);
1901             ChHist[ec][st][rg][iChamber - FirstCh].Y_EfficientRecHits_inSegment.push_back(
1902                 new TH1F(SpecName,
1903                          "Efficient (extrapolated from the segment) RecHit/layer in a segment (local system, whole "
1904                          "chamber);Y, cm; entries",
1905                          nYbins,
1906                          Ymin,
1907                          Ymax));
1908             //
1909             Chan = 200;
1910             minChan = -0.2;
1911             maxChan = 0.2;
1912             sprintf(SpecName, "Phi_InefficientRecHits_inSegment_Ch%d_L%d", iChamber, iLayer);
1913             ChHist[ec][st][rg][iChamber - FirstCh].Phi_InefficientRecHits_inSegment.push_back(
1914                 new TH1F(SpecName,
1915                          "Missing RecHit/layer in a segment (local system, whole chamber);Phi, rad; entries",
1916                          Chan,
1917                          minChan,
1918                          maxChan));
1919             //
1920             sprintf(SpecName, "Phi_EfficientRecHits_inSegment_Ch%d_L%d", iChamber, iLayer);
1921             ChHist[ec][st][rg][iChamber - FirstCh].Phi_EfficientRecHits_inSegment.push_back(
1922                 new TH1F(SpecName,
1923                          "Efficient (extrapolated from the segment) in a segment (local system, whole chamber);Phi, "
1924                          "rad; entries",
1925                          Chan,
1926                          minChan,
1927                          maxChan));
1928           }
1929           //
1930           sprintf(SpecName, "Sim_Rechits_Ch%d", iChamber);
1931           ChHist[ec][st][rg][iChamber - FirstCh].SimRechits =
1932               new TH1F(SpecName, "Existing RecHit (Sim);layers (1-6);entries", nLayer_bins, Layer_min, Layer_max);
1933           //
1934           sprintf(SpecName, "Sim_Simhits_Ch%d", iChamber);
1935           ChHist[ec][st][rg][iChamber - FirstCh].SimSimhits =
1936               new TH1F(SpecName, "Existing SimHit (Sim);layers (1-6);entries", nLayer_bins, Layer_min, Layer_max);
1937           //
1938           /*
1939       sprintf(SpecName,"Sim_Rechits_each_Ch%d",iChamber);
1940       ChHist[ec][st][rg][iChamber-FirstCh].SimRechits_each =
1941         new TH1F(SpecName,"Existing RecHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
1942       //
1943       sprintf(SpecName,"Sim_Simhits_each_Ch%d",iChamber);
1944       ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits_each =
1945         new TH1F(SpecName,"Existing SimHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max);
1946       */
1947           theFile->cd();
1948         }
1949       }
1950     }
1951   }
1952 }
1953 
1954 // Destructor
1955 CSCEfficiency::~CSCEfficiency() {
1956   if (theService)
1957     delete theService;
1958   // Write the histos to a file
1959   theFile->cd();
1960   //
1961   char SpecName[60];
1962   std::vector<float> bins, Efficiency, EffError;
1963   std::vector<float> eff(2);
1964 
1965   //---- loop over chambers
1966   std::map<std::string, bool> chamberTypes;
1967   chamberTypes["ME11"] = false;
1968   chamberTypes["ME12"] = false;
1969   chamberTypes["ME13"] = false;
1970   chamberTypes["ME21"] = false;
1971   chamberTypes["ME22"] = false;
1972   chamberTypes["ME31"] = false;
1973   chamberTypes["ME32"] = false;
1974   chamberTypes["ME41"] = false;
1975 
1976   map<std::string, bool>::iterator iter;
1977   std::cout << " Writing proper histogram structure (patience)..." << std::endl;
1978   for (int ec = 0; ec < 2; ++ec) {
1979     for (int st = 0; st < 4; ++st) {
1980       snprintf(SpecName, sizeof(SpecName), "Stations__E%d_S%d", ec + 1, st + 1);
1981       theFile->cd(SpecName);
1982       StHist[ec][st].segmentChi2_ndf->Write();
1983       StHist[ec][st].hitsInSegment->Write();
1984       StHist[ec][st].AllSegments_eta->Write();
1985       StHist[ec][st].EfficientSegments_eta->Write();
1986       StHist[ec][st].ResidualSegments->Write();
1987       StHist[ec][st].EfficientSegments_XY->Write();
1988       StHist[ec][st].InefficientSegments_XY->Write();
1989       StHist[ec][st].EfficientALCT_momTheta->Write();
1990       StHist[ec][st].InefficientALCT_momTheta->Write();
1991       StHist[ec][st].EfficientCLCT_momPhi->Write();
1992       StHist[ec][st].InefficientCLCT_momPhi->Write();
1993       for (int rg = 0; rg < 3; ++rg) {
1994         if (0 != st && rg > 1) {
1995           continue;
1996         } else if (1 == rg && 3 == st) {
1997           continue;
1998         }
1999         for (int iChamber = FirstCh; iChamber < FirstCh + NumCh; iChamber++) {
2000           if (0 != st && 0 == rg && iChamber > 18) {
2001             continue;
2002           }
2003           snprintf(SpecName, sizeof(SpecName), "Chambers__E%d_S%d_R%d_Chamber_%d", ec + 1, st + 1, rg + 1, iChamber);
2004           theFile->cd(SpecName);
2005 
2006           ChHist[ec][st][rg][iChamber - FirstCh].EfficientRechits_inSegment->Write();
2007           ChHist[ec][st][rg][iChamber - FirstCh].AllSingleHits->Write();
2008           ChHist[ec][st][rg][iChamber - FirstCh].digiAppearanceCount->Write();
2009           ChHist[ec][st][rg][iChamber - FirstCh].EfficientALCT_dydz->Write();
2010           ChHist[ec][st][rg][iChamber - FirstCh].InefficientALCT_dydz->Write();
2011           ChHist[ec][st][rg][iChamber - FirstCh].EfficientCLCT_dxdz->Write();
2012           ChHist[ec][st][rg][iChamber - FirstCh].InefficientCLCT_dxdz->Write();
2013           ChHist[ec][st][rg][iChamber - FirstCh].InefficientSingleHits->Write();
2014           ChHist[ec][st][rg][iChamber - FirstCh].EfficientRechits_good->Write();
2015           ChHist[ec][st][rg][iChamber - FirstCh].EfficientStrips->Write();
2016           ChHist[ec][st][rg][iChamber - FirstCh].StripWiresCorrelations->Write();
2017           ChHist[ec][st][rg][iChamber - FirstCh].NoWires_momTheta->Write();
2018           ChHist[ec][st][rg][iChamber - FirstCh].NoStrips_momPhi->Write();
2019           ChHist[ec][st][rg][iChamber - FirstCh].EfficientWireGroups->Write();
2020           for (unsigned int iLayer = 0; iLayer < 6; iLayer++) {
2021             ChHist[ec][st][rg][iChamber - FirstCh].Y_InefficientRecHits_inSegment[iLayer]->Write();
2022             ChHist[ec][st][rg][iChamber - FirstCh].Y_EfficientRecHits_inSegment[iLayer]->Write();
2023             ChHist[ec][st][rg][iChamber - FirstCh].Phi_InefficientRecHits_inSegment[iLayer]->Write();
2024             ChHist[ec][st][rg][iChamber - FirstCh].Phi_EfficientRecHits_inSegment[iLayer]->Write();
2025           }
2026           ChHist[ec][st][rg][iChamber - FirstCh].SimRechits->Write();
2027           ChHist[ec][st][rg][iChamber - FirstCh].SimSimhits->Write();
2028           /*
2029       ChHist[ec][st][rg][iChamber-FirstCh].SimRechits_each->Write();
2030       ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits_each->Write();
2031       */
2032           //
2033           theFile->cd(SpecName);
2034           theFile->cd();
2035         }
2036       }
2037     }
2038   }
2039   //
2040   snprintf(SpecName, sizeof(SpecName), "AllChambers");
2041   theFile->mkdir(SpecName);
2042   theFile->cd(SpecName);
2043   DataFlow->Write();
2044   TriggersFired->Write();
2045   ALCTPerEvent->Write();
2046   CLCTPerEvent->Write();
2047   recHitsPerEvent->Write();
2048   segmentsPerEvent->Write();
2049   //
2050   theFile->cd(SpecName);
2051   //---- Close the file
2052   theFile->Close();
2053 }
2054 
2055 DEFINE_FWK_MODULE(CSCEfficiency);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!