35 printf(
"\n==enter==CSCEfficiency===== run %u\tevent %llu\tn Analyzed %i\n",iRun,iEvent,
nEventsAnalyzed);
41 if (
printalot) printf(
"\tget handles for digi collections\n");
46 if (
printalot) printf(
"\tpass handles\n");
59 event.getByToken(
sd_token, strips );
62 event.getByToken(
co_token, correlatedlcts );
65 event.getByToken(
sh_token, simhits );
67 event.getByToken(
rh_token, rechits );
68 event.getByToken(
se_token, segments );
69 event.getByToken(
tk_token, trackCollectionH );
73 if (
printalot) printf(
"\tget the CSC geometry.\n");
81 bool triggerPassed =
true;
96 if(
theService->magneticField()->inTesla(gpZero).mag2()<0.1){
104 fillDigiInfo(alcts, clcts, correlatedlcts, wires, strips, simhits, rechits, segments, cscGeom);
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.;
120 for ( reco::MuonCollection::const_iterator
muon = muons->begin();
muon != muons->end(); ++
muon ) {
124 " eta = "<<
muon->eta()<<
" phi = "<<
muon->phi()<<
127 muon->isGlobalMuon()<<
"/"<<
muon->isTrackerMuon()<<
"/"<<
muon->isStandAloneMuon()<<std::endl;
129 if(!(
muon->isTrackerMuon() &&
muon->isGlobalMuon())){
133 double relISO = (
muon->isolationR03().sumPt +
134 muon->isolationR03().emEt +
135 muon->isolationR03().hadEt)/
muon->track()->pt();
137 std::cout<<
" relISO = "<<relISO<<
" emVetoEt = "<<
muon->isolationR03().emVetoEt<<
" caloComp = "<<
146 if(
muon->track()->hitPattern().numberOfValidPixelHits()<1 ||
147 muon->track()->hitPattern().numberOfValidTrackerHits()<11 ||
148 muon->combinedMuon()->hitPattern().numberOfValidMuonHits()<1 ||
149 muon->combinedMuon()->normalizedChi2()>10. ||
150 muon->numberOfMatches()<2){
154 float zOuter =
muon->combinedMuon()->outerPosition().z();
155 float rhoOuter =
muon->combinedMuon()->outerPosition().rho();
156 bool passDepth =
true;
159 if ( fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 540.){
164 else if( fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.){
169 else if ( fabs(zOuter) > 680. && fabs(zOuter) < 880. && rhoOuter < 540.){
176 goodMuons_it.push_back(
muon);
177 if(
muon->track()->momentum().z()>0.){
180 if(
muon->track()->momentum().z()<0.){
196 std::cout<<
" nNegativeZ = "<<nNegativeZ<<
" nPositiveZ = "<<nPositiveZ<<std::endl;
198 if(nNegativeZ>1 || nPositiveZ>1){
201 bool trackOK =
false;
203 std::cout<<
" goodMuons_it.size() = "<<goodMuons_it.size()<<std::endl;
205 for(
size_t iM=0;iM<goodMuons_it.size();++iM){
209 float deltaR =
pow(
track->phi()-goodMuons_it[iM]->track()->phi(),2) +
210 pow(
track->eta()-goodMuons_it[iM]->track()->eta(),2);
211 deltaR =
sqrt(deltaR);
213 std::cout<<
" TR mu match to a tr: deltaR = "<<deltaR<<
" dPt = "<<
214 track->pt()-goodMuons_it[iM]->track()->pt()<<std::endl;
216 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;
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()<<
" nSeg = "<<segments->size()<<std::endl;
253 std::cout<<
"quality undef/loose/tight/high/confirmed/goodIt/size "<<
264 std::cout<<
" track inner position = "<<
track->innerPosition()<<
" outer position = "<<
track->outerPosition()<<std::endl;
265 std::cout<<
"track eta (outer) = "<<
track->outerPosition().eta()<<
" phi (outer) = "<<
266 track->outerPosition().phi()<<std::endl;
267 if(fabs(
track->innerPosition().z())>500.){
271 if(fabs(
track->outerPosition().z())>500.){
289 float dpT_ov_pT = 0.;
290 if(fabs(
track->pt())>0.001){
302 if(!segments->size()){
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;
330 std::cout<<
" dump the very first FTS = "<<debug.
dumpFTS(ftsStart)<<std::endl;
335 if(
track->outerPosition().z()>0){
343 std::vector< CSCDetId > refME;
344 for(
int iS=1;iS<5;++iS){
345 for(
int iR=1;iR<4;++iR){
349 else if(4==iS && iR>1){
352 refME.push_back(
CSCDetId(endcap, iS, iR, chamber));
356 for(
size_t iSt = 0; iSt<refME.size();++iSt){
358 std::cout<<
"loop iStatation = "<<iSt<<std::endl;
359 std::cout<<
"refME[iSt]: st = "<<refME[iSt].station()<<
" rg = "<<refME[iSt].ring()<<std::endl;
361 std::map <std::string, bool> chamberTypes;
362 chamberTypes[
"ME11"] =
false;
363 chamberTypes[
"ME12"] =
false;
364 chamberTypes[
"ME13"] =
false;
365 chamberTypes[
"ME21"] =
false;
366 chamberTypes[
"ME22"] =
false;
367 chamberTypes[
"ME31"] =
false;
368 chamberTypes[
"ME32"] =
false;
369 chamberTypes[
"ME41"] =
false;
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()<<
" iSt = "<<iSt<<std::endl;
400 std::vector <int> coupleOfChambers;
404 for(
size_t iCh =0;iCh<coupleOfChambers.size();++iCh){
406 if (
printalot)
std::cout<<
" Check chamber N = "<<coupleOfChambers.at(iCh)<<std::endl;;
410 [refME[iSt].
ring()-1]
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()<<
" zSurf = "<<bpCh.position().z()<<std::endl;
425 if(!
isIPdata && (zDistInner*zDistOuter>0. || fabs(zDistInner)<15. || fabs(zDistOuter)<15.)){
427 std::cout<<
" Not an intermediate (as defined) point... Skip."<<std::endl;
432 if(fabs(muonOuterZPosition) - fabs(bpCh.position().z())<0 ||
433 fabs(muonOuterZPosition-bpCh.position().z())<15.){
441 tSOSDest =
propagate(ftsStart, cscGeom->
idToDet(cscChamber->geographicalId())->surface());
455 bool inDeadZone =
false;
457 for(
int iLayer = 0;iLayer<6;++iLayer){
458 bool extrapolationPassed =
true;
460 std::cout<<
" iLayer = "<<iLayer<<
" dump FTS init = "<<debug.
dumpFTS(ftsInit)<<std::endl;
462 std::cout<<
"Surface to propagate to: pos = "<<cscChamber->layer(iLayer+1)->surface().position()<<
" eta = " 463 <<cscChamber->layer(iLayer+1)->surface().position().eta()<<
" phi = " 464 <<cscChamber->layer(iLayer+1)->surface().position().phi()<<std::endl;
467 tSOSDest =
propagate(ftsInit, cscChamber->layer(iLayer+1)->surface());
471 getFromFTS(ftsInit, p3_propagated, r3_propagated, charge, cov_propagated);
474 if (
printalot)
std::cout<<
"Propagation between layers not successful - notValid TSOS"<<std::endl;
475 extrapolationPassed =
false;
480 if(extrapolationPassed){
481 GlobalPoint theExtrapolationPoint(r3_propagated.x(),r3_propagated.y(),r3_propagated.z());
482 LocalPoint theLocalPoint = cscChamber->layer(iLayer+1)->toLocal(theExtrapolationPoint);
484 inDeadZone = ( inDeadZone ||
486 refME[iSt].station(), refME[iSt].ring()));
488 std::cout<<
" Candidate chamber: extrapolated LocalPoint = "<<theLocalPoint<<
"inDeadZone = "<<inDeadZone<<std::endl;
533 if (
printalot) 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 GeomDet * idToDet(DetId) const override
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
Geom::Phi< T > phi() const
unsigned long long EventNumber_t
edm::EDGetTokenT< CSCCorrelatedLCTDigiCollection > co_token
bool recSimHitEfficiency(CSCDetId &id, FreeTrajectoryState &ftsChamber)
edm::EDGetTokenT< CSCCLCTDigiCollection > cl_token
const Plane & surface() const
The nominal surface of the GeomDet.
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)
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
const CSCChamber * chamber(CSCDetId id) const
Return the chamber corresponding to given DetId.
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
const PositionType & position() 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