35 printf(
"\n==enter==CSCEfficiency===== run %u\tevent %llu\tn Analyzed %i\n", iRun,
iEvent,
nEventsAnalyzed);
42 printf(
"\tget handles for digi collections\n");
48 printf(
"\tpass handles\n");
64 event.getByToken(
co_token, correlatedlcts);
70 event.getByToken(
se_token, segments);
71 event.getByToken(
tk_token, trackCollectionH);
76 printf(
"\tget the CSC geometry.\n");
82 bool triggerPassed =
true;
97 if (
theService->magneticField()->inTesla(gpZero).mag2() < 0.1) {
108 event.getByLabel(muonTag_,
muons);
111 event.getByLabel(
"offlineBeamSpot", beamSpotHandle);
114 std::vector<reco::MuonCollection::const_iterator> goodMuons_it;
115 unsigned int nPositiveZ = 0;
116 unsigned int nNegativeZ = 0;
117 float muonOuterZPosition = -99999.;
125 <<
" pt = " <<
muon->pt() <<
" eta = " <<
muon->eta() <<
" phi = " <<
muon->phi()
126 <<
" matches = " <<
muon->matches().size()
128 <<
" GLB/TR/STA = " <<
muon->isGlobalMuon() <<
"/" <<
muon->isTrackerMuon() <<
"/" 129 <<
muon->isStandAloneMuon() << std::endl;
131 if (!(
muon->isTrackerMuon() &&
muon->isGlobalMuon())) {
136 (
muon->isolationR03().sumPt +
muon->isolationR03().emEt +
muon->isolationR03().hadEt) /
muon->track()->pt();
138 std::cout <<
" relISO = " << relISO <<
" emVetoEt = " <<
muon->isolationR03().emVetoEt
140 <<
" dxy = " << fabs(
muon->track()->dxy(vertexBeamSpot.
position())) << std::endl;
144 fabs(
muon->track()->dxy(vertexBeamSpot.
position())) > 0.2 ||
muon->pt() < 6.) {
148 if (
muon->track()->hitPattern().numberOfValidPixelHits() < 1 ||
149 muon->track()->hitPattern().numberOfValidTrackerHits() < 11 ||
150 muon->combinedMuon()->hitPattern().numberOfValidMuonHits() < 1 ||
151 muon->combinedMuon()->normalizedChi2() > 10. ||
muon->numberOfMatches() < 2) {
155 float zOuter =
muon->combinedMuon()->outerPosition().z();
156 float rhoOuter =
muon->combinedMuon()->outerPosition().rho();
157 bool passDepth =
true;
160 if (fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 540.) {
165 else if (fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.) {
170 else if (fabs(zOuter) > 680. && fabs(zOuter) < 880. && rhoOuter < 540.) {
177 goodMuons_it.push_back(
muon);
178 if (
muon->track()->momentum().z() > 0.) {
181 if (
muon->track()->momentum().z() < 0.) {
197 std::cout <<
" nNegativeZ = " << nNegativeZ <<
" nPositiveZ = " << nPositiveZ << std::endl;
199 if (nNegativeZ > 1 || nPositiveZ > 1) {
202 bool trackOK =
false;
204 std::cout <<
" goodMuons_it.size() = " << goodMuons_it.size() << std::endl;
206 for (
size_t iM = 0; iM < goodMuons_it.size(); ++iM) {
210 float deltaR =
pow(
track->phi() - goodMuons_it[iM]->track()->phi(), 2) +
211 pow(
track->eta() - goodMuons_it[iM]->track()->eta(), 2);
215 <<
" dPt = " <<
track->pt() - goodMuons_it[iM]->track()->pt() << std::endl;
217 if (
deltaR > 0.01 || fabs(
track->pt() - goodMuons_it[iM]->track()->pt()) > 0.1) {
224 muonOuterZPosition = goodMuons_it[iM]->combinedMuon()->outerPosition().z();
231 std::cout <<
" failed: trackOK " << std::endl;
244 if (
track->outerPosition().z() * trackTwo->outerPosition().z() > 0) {
251 std::cout <<
"i track = " <<
i <<
" P = " <<
track->p() <<
" chi2/ndf = " <<
track->normalizedChi2()
252 <<
" nSeg = " << segments->size() << std::endl;
259 <<
track->qoverpError() << std::endl;
262 <<
" outer position = " <<
track->outerPosition() << std::endl;
263 std::cout <<
"track eta (outer) = " <<
track->outerPosition().eta()
264 <<
" phi (outer) = " <<
track->outerPosition().phi() << std::endl;
265 if (fabs(
track->innerPosition().z()) > 500.) {
267 std::cout <<
" dump inner state MUON detid = " <<
debug.dumpMuonId(innerDetId) << std::endl;
269 if (fabs(
track->outerPosition().z()) > 500.) {
271 std::cout <<
" dump outer state MUON detid = " <<
debug.dumpMuonId(outerDetId) << std::endl;
287 float dpT_ov_pT = 0.;
288 if (fabs(
track->pt()) > 0.001) {
300 if (!segments->size()) {
308 if (
magField && (dpT_ov_pT > 0.5)) {
316 CLHEP::Hep3Vector r3T_inner(
track->innerPosition().x(),
track->innerPosition().y(),
track->innerPosition().z());
317 CLHEP::Hep3Vector r3T(
track->outerPosition().x(),
track->outerPosition().y(),
track->outerPosition().z());
320 CLHEP::Hep3Vector p3T(
track->outerMomentum().x(),
track->outerMomentum().y(),
track->outerMomentum().z());
321 CLHEP::Hep3Vector p3_propagated, r3_propagated;
324 cov_propagated *= 1
e-20;
328 std::cout <<
" p = " <<
track->p() <<
" norm chi2 = " <<
track->normalizedChi2() << std::endl;
329 std::cout <<
" dump the very first FTS = " <<
debug.dumpFTS(ftsStart) << std::endl;
334 if (
track->outerPosition().z() > 0) {
341 std::vector<CSCDetId> refME;
342 for (
int iS = 1; iS < 5; ++iS) {
343 for (
int iR = 1; iR < 4; ++iR) {
344 if (1 != iS && iR > 2) {
346 }
else if (4 == iS && iR > 1) {
353 for (
size_t iSt = 0; iSt < refME.size(); ++iSt) {
355 std::cout <<
"loop iStatation = " << iSt << std::endl;
356 std::cout <<
"refME[iSt]: st = " << refME[iSt].station() <<
" rg = " << refME[iSt].ring() << std::endl;
358 std::map<std::string, bool> chamberTypes;
359 chamberTypes[
"ME11"] =
false;
360 chamberTypes[
"ME12"] =
false;
361 chamberTypes[
"ME13"] =
false;
362 chamberTypes[
"ME21"] =
false;
363 chamberTypes[
"ME22"] =
false;
364 chamberTypes[
"ME31"] =
false;
365 chamberTypes[
"ME32"] =
false;
366 chamberTypes[
"ME41"] =
false;
373 std::cout <<
" dump base iStation detid = " <<
debug.dumpMuonId(detId) << std::endl;
374 std::cout <<
" dump FTS start = " <<
debug.dumpFTS(ftsStart) << std::endl;
381 std::cout <<
" dump FTS end = " <<
debug.dumpFTS(ftsStart) << std::endl;
383 float feta = fabs(r3_propagated.eta());
384 float phi = r3_propagated.phi();
388 map<std::string, bool>::iterator iter;
391 for (iter = chamberTypes.begin(); iter != chamberTypes.end(); iter++) {
394 if (iter->second && (iterations - 1) ==
int(iSt)) {
396 std::cout <<
" Chamber type " << iter->first <<
" is a candidate..." << std::endl;
397 std::cout <<
" station() = " << refME[iSt].station() <<
" ring() = " << refME[iSt].ring()
398 <<
" iSt = " << iSt << std::endl;
400 std::vector<int> coupleOfChambers;
404 for (
size_t iCh = 0; iCh < coupleOfChambers.size(); ++iCh) {
407 std::cout <<
" Check chamber N = " << coupleOfChambers.at(iCh) << std::endl;
411 [coupleOfChambers.at(iCh) -
FirstCh])) {
416 const BoundPlane bpCh = cscGeom->
idToDet(cscChamber->geographicalId())->surface();
418 float dz = fabs(bpCh.position().z() - zFTS);
419 float zDistInner =
track->innerPosition().z() - bpCh.position().z();
420 float zDistOuter =
track->outerPosition().z() - bpCh.position().z();
423 std::cout <<
" zIn = " <<
track->innerPosition().z() <<
" zOut = " <<
track->outerPosition().z()
424 <<
" zSurf = " << bpCh.position().z() << std::endl;
426 if (!
isIPdata && (zDistInner * zDistOuter > 0. || fabs(zDistInner) < 15. ||
427 fabs(zDistOuter) < 15.)) {
429 std::cout <<
" Not an intermediate (as defined) point... Skip." << std::endl;
434 if (fabs(muonOuterZPosition) - fabs(bpCh.position().z()) < 0 ||
435 fabs(muonOuterZPosition - bpCh.position().z()) < 15.) {
443 tSOSDest =
propagate(ftsStart, cscGeom->
idToDet(cscChamber->geographicalId())->surface());
448 std::cout <<
"TSOS not valid! Break." << std::endl;
453 std::cout <<
" info: dz<0.1" << std::endl;
457 bool inDeadZone =
false;
459 for (
int iLayer = 0; iLayer < 6; ++iLayer) {
460 bool extrapolationPassed =
true;
462 std::cout <<
" iLayer = " << iLayer <<
" dump FTS init = " <<
debug.dumpFTS(ftsInit) << std::endl;
463 std::cout <<
" dump detid = " <<
debug.dumpMuonId(cscChamber->geographicalId()) << std::endl;
464 std::cout <<
"Surface to propagate to: pos = " << cscChamber->layer(iLayer + 1)->surface().position()
465 <<
" eta = " << cscChamber->layer(iLayer + 1)->surface().position().eta()
466 <<
" phi = " << cscChamber->layer(iLayer + 1)->surface().position().phi() << std::endl;
469 tSOSDest =
propagate(ftsInit, cscChamber->layer(iLayer + 1)->surface());
473 std::cout <<
" Propagation between layers successful: dump FTS end = " <<
debug.dumpFTS(ftsInit)
478 std::cout <<
"Propagation between layers not successful - notValid TSOS" << std::endl;
479 extrapolationPassed =
false;
484 if (extrapolationPassed) {
485 GlobalPoint theExtrapolationPoint(r3_propagated.x(), r3_propagated.y(), r3_propagated.z());
486 LocalPoint theLocalPoint = cscChamber->layer(iLayer + 1)->toLocal(theExtrapolationPoint);
488 inDeadZone = (inDeadZone ||
490 theLocalPoint.x(), theLocalPoint.y(), refME[iSt].station(), refME[iSt].ring()));
492 std::cout <<
" Candidate chamber: extrapolated LocalPoint = " << theLocalPoint
493 <<
"inDeadZone = " << inDeadZone << std::endl;
508 std::cout <<
"Do efficiencies..." << std::endl;
511 bool angle_flag =
true;
524 std::cout <<
" Not in active area for all layers" << std::endl;
534 std::cout <<
" TSOS not valid..." << std::endl;
540 printf(
"==exit===CSCEfficiency===== run %u\tevent %llu\n\n", iRun,
iEvent);
bool applyTrigger(edm::Handle< edm::TriggerResults > &hltR, const edm::TriggerNames &triggerNames)
edm::EDGetTokenT< CSCStripDigiCollection > sd_token
const CSCChamber * chamber(CSCDetId id) const
Return the chamber corresponding to given DetId.
const Point & position() const
position
void chooseDirection(CLHEP::Hep3Vector &innerPosition, CLHEP::Hep3Vector &outerPosition)
void ringCandidates(int station, float absEta, std::map< std::string, bool > &chamberTypes)
Geom::Phi< T > phi() const
T const * product() const
unsigned long long EventNumber_t
edm::EDGetTokenT< CSCCorrelatedLCTDigiCollection > co_token
bool recSimHitEfficiency(CSCDetId &id, FreeTrajectoryState &ftsChamber)
edm::EDGetTokenT< CSCCLCTDigiCollection > cl_token
bool emptyChambers[2][4][4][(36 - 1+1)]
bool getAbsoluteEfficiency
float caloCompatibility(const reco::Muon &muon)
edm::EDGetTokenT< CSCSegmentCollection > se_token
FreeTrajectoryState getFromCLHEP(const CLHEP::Hep3Vector &p3, const CLHEP::Hep3Vector &r3, int charge, const AlgebraicSymMatrix66 &cov, const MagneticField *field)
GlobalPoint position() const
edm::ESGetToken< CSCGeometry, MuonGeometryRecord > geomToken_
edm::EDGetTokenT< edm::PSimHitContainer > sh_token
MuonServiceProxy * theService
edm::EDGetTokenT< CSCWireDigiCollection > wd_token
void getFromFTS(const FreeTrajectoryState &fts, CLHEP::Hep3Vector &p3, CLHEP::Hep3Vector &r3, int &charge, AlgebraicSymMatrix66 &cov)
DetId geographicalId() const
The label of this GeomDet.
bool efficienciesPerChamber(CSCDetId &id, const CSCChamber *cscChamber, FreeTrajectoryState &ftsChamber)
void chamberCandidates(int station, int ring, float phi, std::vector< int > &coupleOfChambers)
bool stripWire_Efficiencies(CSCDetId &cscDetId, FreeTrajectoryState &ftsChamber)
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
edm::EDGetTokenT< CSCRecHit2DCollection > rh_token
const Plane & surface() const
The nominal surface of the GeomDet.
void fillDigiInfo(edm::Handle< CSCALCTDigiCollection > &alcts, edm::Handle< CSCCLCTDigiCollection > &clcts, edm::Handle< CSCCorrelatedLCTDigiCollection > &correlatedlcts, edm::Handle< CSCWireDigiCollection > &wires, edm::Handle< CSCStripDigiCollection > &strips, edm::Handle< edm::PSimHitContainer > &simhits, edm::Handle< CSCRecHit2DCollection > &rechits, edm::Handle< CSCSegmentCollection > &segments, edm::ESHandle< CSCGeometry > &cscGeom)
const PositionType & position() const
bool recHitSegment_Efficiencies(CSCDetId &cscDetId, const CSCChamber *cscChamber, FreeTrajectoryState &ftsChamber)
unsigned int printout_NEvents
bool inSensitiveLocalRegion(double xLocal, double yLocal, int station, int ring)
FreeTrajectoryState const * freeState(bool withErrors=true) const
strips
#turn off noise in all subdetectors simHcalUnsuppressedDigis.doNoise = False mix.digitizers.hcal.doNoise = False simEcalUnsuppressedDigis.doNoise = False mix.digitizers.ecal.doNoise = False simEcalUnsuppressedDigis.doESNoise = False simSiPixelDigis.AddNoise = False mix.digitizers.pixel.AddNoise = False simSiStripDigis.Noise = False mix.digitizers.strip.AddNoise = False
unsigned int minTrackHits
edm::EDGetTokenT< CSCALCTDigiCollection > al_token
edm::EDGetTokenT< edm::TriggerResults > ht_token
Power< A, B >::type pow(const A &a, const B &b)
TrajectoryStateOnSurface propagate(FreeTrajectoryState &ftsStart, const BoundPlane &bp)
edm::EDGetTokenT< edm::View< reco::Track > > tk_token
const GeomDet * idToDet(DetId) const override