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");
83 bool triggerPassed =
true;
98 if (
theService->magneticField()->inTesla(gpZero).mag2() < 0.1) {
105 fillDigiInfo(alcts, clcts, correlatedlcts, wires, strips, simhits, rechits, segments, cscGeom);
109 event.getByLabel(muonTag_, muons);
112 event.getByLabel(
"offlineBeamSpot", beamSpotHandle);
115 std::vector<reco::MuonCollection::const_iterator> goodMuons_it;
116 unsigned int nPositiveZ = 0;
117 unsigned int nNegativeZ = 0;
118 float muonOuterZPosition = -99999.;
121 std::cout <<
" muons.size() = " << muons->size() << std::endl;
122 for (reco::MuonCollection::const_iterator
muon = muons->begin();
muon != muons->end(); ++
muon) {
125 std::cout <<
" iMuon = " <<
muon - muons->begin() <<
" charge = " <<
muon->charge() <<
" p = " <<
muon->p()
126 <<
" pt = " <<
muon->pt() <<
" eta = " <<
muon->eta() <<
" phi = " <<
muon->phi()
127 <<
" matches = " <<
muon->matches().size()
129 <<
" GLB/TR/STA = " <<
muon->isGlobalMuon() <<
"/" <<
muon->isTrackerMuon() <<
"/"
130 <<
muon->isStandAloneMuon() << std::endl;
132 if (!(
muon->isTrackerMuon() &&
muon->isGlobalMuon())) {
137 (
muon->isolationR03().sumPt +
muon->isolationR03().emEt +
muon->isolationR03().hadEt) /
muon->track()->pt();
139 std::cout <<
" relISO = " << relISO <<
" emVetoEt = " <<
muon->isolationR03().emVetoEt
141 <<
" dxy = " << fabs(
muon->track()->dxy(vertexBeamSpot.
position())) << std::endl;
145 fabs(
muon->track()->dxy(vertexBeamSpot.
position())) > 0.2 ||
muon->pt() < 6.) {
149 if (
muon->track()->hitPattern().numberOfValidPixelHits() < 1 ||
150 muon->track()->hitPattern().numberOfValidTrackerHits() < 11 ||
151 muon->combinedMuon()->hitPattern().numberOfValidMuonHits() < 1 ||
152 muon->combinedMuon()->normalizedChi2() > 10. ||
muon->numberOfMatches() < 2) {
156 float zOuter =
muon->combinedMuon()->outerPosition().z();
157 float rhoOuter =
muon->combinedMuon()->outerPosition().rho();
158 bool passDepth =
true;
161 if (fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 540.) {
166 else if (fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.) {
171 else if (fabs(zOuter) > 680. && fabs(zOuter) < 880. && rhoOuter < 540.) {
178 goodMuons_it.push_back(
muon);
179 if (
muon->track()->momentum().z() > 0.) {
182 if (
muon->track()->momentum().z() < 0.) {
191 std::cout <<
"Start track loop over " << trackCollection.
size() <<
" tracks" << std::endl;
198 std::cout <<
" nNegativeZ = " << nNegativeZ <<
" nPositiveZ = " << nPositiveZ << std::endl;
200 if (nNegativeZ > 1 || nPositiveZ > 1) {
203 bool trackOK =
false;
205 std::cout <<
" goodMuons_it.size() = " << goodMuons_it.size() << std::endl;
207 for (
size_t iM = 0; iM < goodMuons_it.size(); ++iM) {
211 float deltaR =
pow(
track->phi() - goodMuons_it[iM]->track()->phi(), 2) +
212 pow(
track->eta() - goodMuons_it[iM]->track()->eta(), 2);
213 deltaR =
sqrt(deltaR);
215 std::cout <<
" TR mu match to a tr: deltaR = " << deltaR
216 <<
" dPt = " <<
track->pt() - goodMuons_it[iM]->track()->pt() << std::endl;
218 if (deltaR > 0.01 || fabs(
track->pt() - goodMuons_it[iM]->track()->pt()) > 0.1) {
225 muonOuterZPosition = goodMuons_it[iM]->combinedMuon()->outerPosition().z();
232 std::cout <<
" failed: trackOK " << std::endl;
238 if (trackCollection.
size() > 2) {
242 if (!
i && 2 == trackCollection.
size()) {
245 if (
track->outerPosition().z() * trackTwo->outerPosition().z() > 0) {
252 std::cout <<
"i track = " <<
i <<
" P = " <<
track->p() <<
" chi2/ndf = " <<
track->normalizedChi2()
253 <<
" nSeg = " << segments->size() << std::endl;
260 <<
track->qoverpError() << std::endl;
263 <<
" outer position = " <<
track->outerPosition() << std::endl;
264 std::cout <<
"track eta (outer) = " <<
track->outerPosition().eta()
265 <<
" phi (outer) = " <<
track->outerPosition().phi() << std::endl;
266 if (fabs(
track->innerPosition().z()) > 500.) {
268 std::cout <<
" dump inner state MUON detid = " << debug.
dumpMuonId(innerDetId) << std::endl;
270 if (fabs(
track->outerPosition().z()) > 500.) {
272 std::cout <<
" dump outer state MUON detid = " << debug.
dumpMuonId(outerDetId) << std::endl;
288 float dpT_ov_pT = 0.;
289 if (fabs(
track->pt()) > 0.001) {
301 if (!segments->size()) {
309 if (
magField && (dpT_ov_pT > 0.5)) {
317 CLHEP::Hep3Vector r3T_inner(
track->innerPosition().x(),
track->innerPosition().y(),
track->innerPosition().z());
318 CLHEP::Hep3Vector r3T(
track->outerPosition().x(),
track->outerPosition().y(),
track->outerPosition().z());
321 CLHEP::Hep3Vector p3T(
track->outerMomentum().x(),
track->outerMomentum().y(),
track->outerMomentum().z());
322 CLHEP::Hep3Vector p3_propagated, r3_propagated;
325 cov_propagated *= 1
e-20;
329 std::cout <<
" p = " <<
track->p() <<
" norm chi2 = " <<
track->normalizedChi2() << std::endl;
330 std::cout <<
" dump the very first FTS = " << debug.
dumpFTS(ftsStart) << std::endl;
335 if (
track->outerPosition().z() > 0) {
342 std::vector<CSCDetId> refME;
343 for (
int iS = 1; iS < 5; ++iS) {
344 for (
int iR = 1; iR < 4; ++iR) {
345 if (1 != iS && iR > 2) {
347 }
else if (4 == iS && iR > 1) {
350 refME.push_back(
CSCDetId(endcap, iS, iR, chamber));
354 for (
size_t iSt = 0; iSt < refME.size(); ++iSt) {
356 std::cout <<
"loop iStatation = " << iSt << std::endl;
357 std::cout <<
"refME[iSt]: st = " << refME[iSt].station() <<
" rg = " << refME[iSt].ring() << std::endl;
359 std::map<std::string, bool> chamberTypes;
360 chamberTypes[
"ME11"] =
false;
361 chamberTypes[
"ME12"] =
false;
362 chamberTypes[
"ME13"] =
false;
363 chamberTypes[
"ME21"] =
false;
364 chamberTypes[
"ME22"] =
false;
365 chamberTypes[
"ME31"] =
false;
366 chamberTypes[
"ME32"] =
false;
367 chamberTypes[
"ME41"] =
false;
368 const CSCChamber *cscChamber_base = cscGeom->chamber(refME[iSt].chamberId());
371 std::cout <<
" base iStation : eta = " << cscGeom->idToDet(detId)->surface().position().eta()
372 <<
" phi = " << cscGeom->idToDet(detId)->surface().position().phi()
373 <<
" y = " << cscGeom->idToDet(detId)->surface().position().y() << std::endl;
378 tSOSDest =
propagate(ftsStart, cscGeom->idToDet(detId)->surface());
383 getFromFTS(ftsStart, p3_propagated, r3_propagated, charge, cov_propagated);
384 float feta = fabs(r3_propagated.eta());
385 float phi = r3_propagated.phi();
389 map<std::string, bool>::iterator iter;
392 for (iter = chamberTypes.begin(); iter != chamberTypes.end(); iter++) {
395 if (iter->second && (iterations - 1) == int(iSt)) {
397 std::cout <<
" Chamber type " << iter->first <<
" is a candidate..." << std::endl;
398 std::cout <<
" station() = " << refME[iSt].station() <<
" ring() = " << refME[iSt].ring()
399 <<
" iSt = " << iSt << std::endl;
401 std::vector<int> coupleOfChambers;
405 for (
size_t iCh = 0; iCh < coupleOfChambers.size(); ++iCh) {
408 std::cout <<
" Check chamber N = " << coupleOfChambers.at(iCh) << std::endl;
412 [coupleOfChambers.at(iCh) -
FirstCh])) {
416 const CSCChamber *cscChamber = cscGeom->chamber(theCSCId.chamberId());
417 const BoundPlane bpCh = cscGeom->idToDet(cscChamber->geographicalId())->surface();
419 float dz = fabs(bpCh.position().z() - zFTS);
420 float zDistInner =
track->innerPosition().z() - bpCh.position().z();
421 float zDistOuter =
track->outerPosition().z() - bpCh.position().z();
424 std::cout <<
" zIn = " <<
track->innerPosition().z() <<
" zOut = " <<
track->outerPosition().z()
425 <<
" zSurf = " << bpCh.position().z() << std::endl;
427 if (!
isIPdata && (zDistInner * zDistOuter > 0. || fabs(zDistInner) < 15. ||
428 fabs(zDistOuter) < 15.)) {
430 std::cout <<
" Not an intermediate (as defined) point... Skip." << std::endl;
435 if (fabs(muonOuterZPosition) - fabs(bpCh.position().z()) < 0 ||
436 fabs(muonOuterZPosition - bpCh.position().z()) < 15.) {
444 tSOSDest =
propagate(ftsStart, cscGeom->idToDet(cscChamber->geographicalId())->surface());
449 std::cout <<
"TSOS not valid! Break." << std::endl;
454 std::cout <<
" info: dz<0.1" << std::endl;
458 bool inDeadZone =
false;
460 for (
int iLayer = 0; iLayer < 6; ++iLayer) {
461 bool extrapolationPassed =
true;
463 std::cout <<
" iLayer = " << iLayer <<
" dump FTS init = " << debug.
dumpFTS(ftsInit) << std::endl;
464 std::cout <<
" dump detid = " << debug.
dumpMuonId(cscChamber->geographicalId()) << std::endl;
465 std::cout <<
"Surface to propagate to: pos = " << cscChamber->layer(iLayer + 1)->surface().position()
466 <<
" eta = " << cscChamber->layer(iLayer + 1)->surface().position().eta()
467 <<
" phi = " << cscChamber->layer(iLayer + 1)->surface().position().phi() << std::endl;
470 tSOSDest =
propagate(ftsInit, cscChamber->layer(iLayer + 1)->surface());
474 std::cout <<
" Propagation between layers successful: dump FTS end = " << debug.
dumpFTS(ftsInit)
476 getFromFTS(ftsInit, p3_propagated, r3_propagated, charge, cov_propagated);
479 std::cout <<
"Propagation between layers not successful - notValid TSOS" << std::endl;
480 extrapolationPassed =
false;
485 if (extrapolationPassed) {
486 GlobalPoint theExtrapolationPoint(r3_propagated.x(), r3_propagated.y(), r3_propagated.z());
487 LocalPoint theLocalPoint = cscChamber->layer(iLayer + 1)->toLocal(theExtrapolationPoint);
489 inDeadZone = (inDeadZone ||
491 theLocalPoint.x(), theLocalPoint.y(), refME[iSt].station(), refME[iSt].ring()));
493 std::cout <<
" Candidate chamber: extrapolated LocalPoint = " << theLocalPoint
494 <<
"inDeadZone = " << inDeadZone << std::endl;
509 std::cout <<
"Do efficiencies..." << std::endl;
512 bool angle_flag =
true;
525 std::cout <<
" Not in active area for all layers" << std::endl;
535 std::cout <<
" TSOS not valid..." << std::endl;
541 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
void chooseDirection(CLHEP::Hep3Vector &innerPosition, CLHEP::Hep3Vector &outerPosition)
void ringCandidates(int station, float absEta, std::map< std::string, bool > &chamberTypes)
unsigned long long EventNumber_t
edm::EDGetTokenT< CSCCorrelatedLCTDigiCollection > co_token
bool recSimHitEfficiency(CSCDetId &id, FreeTrajectoryState &ftsChamber)
edm::EDGetTokenT< CSCCLCTDigiCollection > cl_token
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)
std::string dumpMuonId(const DetId &id) const
std::string dumpFTS(const FreeTrajectoryState &fts) const
edm::ESGetToken< CSCGeometry, MuonGeometryRecord > geomToken_
edm::EDGetTokenT< edm::PSimHitContainer > sh_token
MuonServiceProxy * theService
edm::EDGetTokenT< CSCWireDigiCollection > wd_token
tuple strips
#turn off noise in all subdetectors simHcalUnsuppressedDigis.doNoise = False mix.digitizers.hcal.doNoise = False simEcalUnsuppressedDigis.doNoise = False mix.digitizers.ecal.doNoise = False simEcalUnsuppressedDigis.doESNoise = False simSiPixelDigis.AddNoise = False mix.digitizers.pixel.AddNoise = False simSiStripDigis.Noise = False mix.digitizers.strip.AddNoise = False
FreeTrajectoryState const * freeState(bool withErrors=true) const
printf("params %d %f %f %f\n", minT, eps, errmax, chi2max)
DetId geographicalId() const
The label of this GeomDet.
void getFromFTS(const FreeTrajectoryState &fts, CLHEP::Hep3Vector &p3, CLHEP::Hep3Vector &r3, int &charge, AlgebraicSymMatrix66 &cov)
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
GlobalPoint position() const
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)
T const * product() 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)
bool emptyChambers[2][4][4][(36-1+1)]
const Point & position() const
position
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
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