28 int iRun =
event.id().run();
29 int iEvent =
event.id().event();
32 printf(
"\n==enter==CSCEfficiency===== run %i\tevent %i\tn Analyzed %i\n",iRun,iEvent,
nEventsAnalyzed);
38 if (
printalot) printf(
"\tget handles for digi collections\n");
43 if (
printalot) printf(
"\tpass handles\n");
56 event.getByToken(
sd_token, strips );
59 event.getByToken(
co_token, correlatedlcts );
62 event.getByToken(
sh_token, simhits );
64 event.getByToken(
rh_token, rechits );
65 event.getByToken(
se_token, segments );
66 event.getByToken(
tk_token, trackCollectionH );
70 if (
printalot) printf(
"\tget the CSC geometry.\n");
78 bool triggerPassed =
true;
93 if(
theService->magneticField()->inTesla(gpZero).mag2()<0.1){
101 fillDigiInfo(alcts, clcts, correlatedlcts, wires, strips, simhits, rechits, segments, cscGeom);
105 event.getByLabel(muonTag_,muons);
108 event.getByLabel(
"offlineBeamSpot", beamSpotHandle);
111 std::vector <reco::MuonCollection::const_iterator> goodMuons_it;
112 unsigned int nPositiveZ = 0;
113 unsigned int nNegativeZ = 0;
114 float muonOuterZPosition = -99999.;
117 for ( reco::MuonCollection::const_iterator
muon = muons->begin();
muon != muons->end(); ++
muon ) {
121 " eta = "<<
muon->eta()<<
" phi = "<<
muon->phi()<<
124 muon->isGlobalMuon()<<
"/"<<
muon->isTrackerMuon()<<
"/"<<
muon->isStandAloneMuon()<<std::endl;
126 if(!(
muon->isTrackerMuon() &&
muon->isGlobalMuon())){
130 double relISO = (
muon->isolationR03().sumPt +
131 muon->isolationR03().emEt +
132 muon->isolationR03().hadEt)/
muon->track()->pt();
134 std::cout<<
" relISO = "<<relISO<<
" emVetoEt = "<<
muon->isolationR03().emVetoEt<<
" caloComp = "<<
143 if(
muon->track()->hitPattern().numberOfValidPixelHits()<1 ||
144 muon->track()->hitPattern().numberOfValidTrackerHits()<11 ||
145 muon->combinedMuon()->hitPattern().numberOfValidMuonHits()<1 ||
146 muon->combinedMuon()->normalizedChi2()>10. ||
147 muon->numberOfMatches()<2){
151 float zOuter =
muon->combinedMuon()->outerPosition().z();
152 float rhoOuter =
muon->combinedMuon()->outerPosition().rho();
153 bool passDepth =
true;
156 if ( fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 540.){
161 else if( fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.){
166 else if ( fabs(zOuter) > 680. && fabs(zOuter) < 880. && rhoOuter < 540.){
173 goodMuons_it.push_back(
muon);
174 if(
muon->track()->momentum().z()>0.){
177 if(
muon->track()->momentum().z()<0.){
193 std::cout<<
" nNegativeZ = "<<nNegativeZ<<
" nPositiveZ = "<<nPositiveZ<<std::endl;
195 if(nNegativeZ>1 || nPositiveZ>1){
198 bool trackOK =
false;
200 std::cout<<
" goodMuons_it.size() = "<<goodMuons_it.size()<<std::endl;
202 for(
size_t iM=0;iM<goodMuons_it.size();++iM){
206 float deltaR =
pow(track->phi()-goodMuons_it[iM]->track()->phi(),2) +
207 pow(track->eta()-goodMuons_it[iM]->track()->eta(),2);
208 deltaR =
sqrt(deltaR);
210 std::cout<<
" TR mu match to a tr: deltaR = "<<deltaR<<
" dPt = "<<
211 track->pt()-goodMuons_it[iM]->track()->pt()<<std::endl;
213 if(deltaR>0.01 || fabs(track->pt()-goodMuons_it[iM]->track()->pt())>0.1 ){
221 muonOuterZPosition = goodMuons_it[iM]->combinedMuon()->outerPosition().z();
228 std::cout<<
" failed: trackOK "<<std::endl;
235 if(trackCollection.
size()>2){
239 if(!
i && 2==trackCollection.
size()){
242 if(track->outerPosition().z()*trackTwo->outerPosition().z()>0){
249 std::cout<<
"i track = "<<
i<<
" P = "<<track->p()<<
" chi2/ndf = "<<track->normalizedChi2()<<
" nSeg = "<<segments->size()<<std::endl;
250 std::cout<<
"quality undef/loose/tight/high/confirmed/goodIt/size "<<
259 std::cout<<
" pt = "<< track->pt()<<
" +-"<<track->ptError()<<
" q/pt = "<<track->qoverp()<<
" +- "<<track->qoverpError()<<std::endl;
261 std::cout<<
" track inner position = "<<track->innerPosition()<<
" outer position = "<<track->outerPosition()<<std::endl;
262 std::cout<<
"track eta (outer) = "<<track->outerPosition().eta()<<
" phi (outer) = "<<
263 track->outerPosition().phi()<<std::endl;
264 if(fabs(track->innerPosition().z())>500.){
265 DetId innerDetId(track->innerDetId());
268 if(fabs(track->outerPosition().z())>500.){
269 DetId outerDetId(track->outerDetId());
273 std::cout<<
" nHits = "<<track->found()<<std::endl;
286 float dpT_ov_pT = 0.;
287 if(fabs(track->pt())>0.001){
288 dpT_ov_pT = track->ptError()/ track->pt();
299 if(!segments->size()){
314 CLHEP::Hep3Vector r3T_inner(track->innerPosition().x(),track->innerPosition().y(),track->innerPosition().z());
315 CLHEP::Hep3Vector r3T(track->outerPosition().x(),track->outerPosition().y(),track->outerPosition().z());
318 CLHEP::Hep3Vector p3T(track->outerMomentum().x(),track->outerMomentum().y(),track->outerMomentum().z());
319 CLHEP::Hep3Vector p3_propagated, r3_propagated;
322 cov_propagated *= 1
e-20;
323 int charge = track->charge();
326 std::cout<<
" p = "<<track->p()<<
" norm chi2 = "<<track->normalizedChi2()<<std::endl;
327 std::cout<<
" dump the very first FTS = "<<debug.
dumpFTS(ftsStart)<<std::endl;
332 if(track->outerPosition().z()>0){
340 std::vector< CSCDetId > refME;
341 for(
int iS=1;iS<5;++iS){
342 for(
int iR=1;iR<4;++iR){
346 else if(4==iS && iR>1){
349 refME.push_back(
CSCDetId(endcap, iS, iR, chamber));
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;
367 const CSCChamber* cscChamber_base = cscGeom->chamber(refME[iSt].chamberId());
370 std::cout<<
" base iStation : eta = "<<cscGeom->idToDet(detId)->surface().position().eta()<<
" phi = "<<
371 cscGeom->idToDet(detId)->surface().position().phi() <<
" y = " <<cscGeom->idToDet(detId)->surface().position().y()<<std::endl;
376 tSOSDest =
propagate(ftsStart, cscGeom->idToDet(detId)->surface());
380 getFromFTS(ftsStart, p3_propagated, r3_propagated, charge, cov_propagated);
381 float feta = fabs(r3_propagated.eta());
382 float phi = r3_propagated.phi();
386 map<std::string,bool>::iterator
iter;
389 for( iter = chamberTypes.begin(); iter != chamberTypes.end(); iter++ ) {
392 if(iter->second && (iterations-1)==int(iSt)){
394 std::cout<<
" Chamber type "<< iter->first<<
" is a candidate..."<<std::endl;
395 std::cout<<
" station() = "<< refME[iSt].station()<<
" ring() = "<<refME[iSt].ring()<<
" iSt = "<<iSt<<std::endl;
397 std::vector <int> coupleOfChambers;
401 for(
size_t iCh =0;iCh<coupleOfChambers.size();++iCh){
403 if (
printalot)
std::cout<<
" Check chamber N = "<<coupleOfChambers.at(iCh)<<std::endl;;
407 [refME[iSt].
ring()-1]
408 [coupleOfChambers.at(iCh)-
FirstCh])){
412 const CSCChamber* cscChamber = cscGeom->chamber(theCSCId.chamberId());
413 const BoundPlane bpCh = cscGeom->idToDet(cscChamber->geographicalId())->surface();
415 float dz = fabs(bpCh.position().z() - zFTS);
416 float zDistInner = track->innerPosition().z() - bpCh.position().z();
417 float zDistOuter = track->outerPosition().z() - bpCh.position().z();
420 std::cout<<
" zIn = "<<track->innerPosition().z()<<
" zOut = "<<track->outerPosition().z()<<
" zSurf = "<<bpCh.position().z()<<std::endl;
422 if(!
isIPdata && (zDistInner*zDistOuter>0. || fabs(zDistInner)<15. || fabs(zDistOuter)<15.)){
424 std::cout<<
" Not an intermediate (as defined) point... Skip."<<std::endl;
428 if(
isIPdata && fabs(track->eta())<1.8){
429 if(fabs(muonOuterZPosition) - fabs(bpCh.position().z())<0 ||
430 fabs(muonOuterZPosition-bpCh.position().z())<15.){
438 tSOSDest =
propagate(ftsStart, cscGeom->idToDet(cscChamber->geographicalId())->surface());
452 bool inDeadZone =
false;
454 for(
int iLayer = 0;iLayer<6;++iLayer){
455 bool extrapolationPassed =
true;
457 std::cout<<
" iLayer = "<<iLayer<<
" dump FTS init = "<<debug.
dumpFTS(ftsInit)<<std::endl;
459 std::cout<<
"Surface to propagate to: pos = "<<cscChamber->layer(iLayer+1)->surface().position()<<
" eta = "
460 <<cscChamber->layer(iLayer+1)->surface().position().eta()<<
" phi = "
461 <<cscChamber->layer(iLayer+1)->surface().position().phi()<<std::endl;
464 tSOSDest =
propagate(ftsInit, cscChamber->layer(iLayer+1)->surface());
468 getFromFTS(ftsInit, p3_propagated, r3_propagated, charge, cov_propagated);
471 if (
printalot)
std::cout<<
"Propagation between layers not successful - notValid TSOS"<<std::endl;
472 extrapolationPassed =
false;
477 if(extrapolationPassed){
478 GlobalPoint theExtrapolationPoint(r3_propagated.x(),r3_propagated.y(),r3_propagated.z());
479 LocalPoint theLocalPoint = cscChamber->layer(iLayer+1)->toLocal(theExtrapolationPoint);
481 inDeadZone = ( inDeadZone ||
483 refME[iSt].station(), refME[iSt].ring()));
485 std::cout<<
" Candidate chamber: extrapolated LocalPoint = "<<theLocalPoint<<
"inDeadZone = "<<inDeadZone<<std::endl;
530 if (
printalot) printf(
"==exit===CSCEfficiency===== run %i\tevent %i\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)
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
edm::EDGetTokenT< CSCCorrelatedLCTDigiCollection > co_token
bool recSimHitEfficiency(CSCDetId &id, FreeTrajectoryState &ftsChamber)
edm::EDGetTokenT< CSCCLCTDigiCollection > cl_token
Strings const & triggerNames() const
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::EDGetTokenT< edm::PSimHitContainer > sh_token
MuonServiceProxy * theService
edm::EDGetTokenT< CSCWireDigiCollection > wd_token
FreeTrajectoryState const * freeState(bool withErrors=true) const
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)
double deltaR(double eta1, double eta2, double phi1, double phi2)
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)
bool recHitSegment_Efficiencies(CSCDetId &cscDetId, const CSCChamber *cscChamber, FreeTrajectoryState &ftsChamber)
unsigned int printout_NEvents
T const * product() const
bool inSensitiveLocalRegion(double xLocal, double yLocal, int station, int ring)
bool emptyChambers[2][4][4][(36-1+1)]
const Point & position() const
position
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
Public Member Functions inherited from edm::EDFilter
1.8.5
