30 int iRun =
event.id().run();
31 int iEvent =
event.id().event();
34 printf(
"\n==enter==CSCEfficiency===== run %i\tevent %i\tn Analyzed %i\n",iRun,iEvent,
nEventsAnalyzed);
40 if (
printalot) printf(
"\tget handles for digi collections\n");
45 if (
printalot) printf(
"\tpass handles\n");
71 event.getByLabel(
tracksTag,trackCollectionH);
75 if (
printalot) printf(
"\tget the CSC geometry.\n");
83 bool triggerPassed =
true;
98 if(
theService->magneticField()->inTesla(gpZero).mag2()<0.1){
106 fillDigiInfo(alcts, clcts, correlatedlcts, wires, strips, simhits, rechits, segments, cscGeom);
110 event.getByLabel(muonTag_,muons);
113 event.getByLabel(
"offlineBeamSpot", beamSpotHandle);
116 std::vector <reco::MuonCollection::const_iterator> goodMuons_it;
117 unsigned int nPositiveZ = 0;
118 unsigned int nNegativeZ = 0;
119 float muonOuterZPosition = -99999.;
122 for ( reco::MuonCollection::const_iterator
muon = muons->begin();
muon != muons->end(); ++
muon ) {
126 " eta = "<<
muon->eta()<<
" phi = "<<
muon->phi()<<
129 muon->isGlobalMuon()<<
"/"<<
muon->isTrackerMuon()<<
"/"<<
muon->isStandAloneMuon()<<std::endl;
131 if(!(
muon->isTrackerMuon() &&
muon->isGlobalMuon())){
135 double relISO = (
muon->isolationR03().sumPt +
136 muon->isolationR03().emEt +
137 muon->isolationR03().hadEt)/
muon->track()->pt();
139 std::cout<<
" relISO = "<<relISO<<
" emVetoEt = "<<
muon->isolationR03().emVetoEt<<
" caloComp = "<<
148 if(
muon->track()->hitPattern().numberOfValidPixelHits()<1 ||
149 muon->track()->hitPattern().numberOfValidTrackerHits()<11 ||
150 muon->combinedMuon()->hitPattern().numberOfValidMuonHits()<1 ||
151 muon->combinedMuon()->normalizedChi2()>10. ||
152 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.){
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<<
" dPt = "<<
216 track->pt()-goodMuons_it[iM]->track()->pt()<<std::endl;
218 if(deltaR>0.01 || fabs(
track->pt()-goodMuons_it[iM]->track()->pt())>0.1 ){
226 muonOuterZPosition = goodMuons_it[iM]->combinedMuon()->outerPosition().z();
233 std::cout<<
" failed: trackOK "<<std::endl;
240 if(trackCollection.
size()>2){
244 if(!
i && 2==trackCollection.
size()){
247 if(
track->outerPosition().z()*trackTwo->outerPosition().z()>0){
254 std::cout<<
"i track = "<<
i<<
" P = "<<
track->p()<<
" chi2/ndf = "<<
track->normalizedChi2()<<
" nSeg = "<<segments->size()<<std::endl;
255 std::cout<<
"quality undef/loose/tight/high/confirmed/goodIt/size "<<
266 std::cout<<
" track inner position = "<<
track->innerPosition()<<
" outer position = "<<
track->outerPosition()<<std::endl;
267 std::cout<<
"track eta (outer) = "<<
track->outerPosition().eta()<<
" phi (outer) = "<<
268 track->outerPosition().phi()<<std::endl;
269 if(fabs(
track->innerPosition().z())>500.){
273 if(fabs(
track->outerPosition().z())>500.){
291 float dpT_ov_pT = 0.;
292 if(fabs(
track->pt())>0.001){
304 if(!segments->size()){
319 CLHEP::Hep3Vector r3T_inner(
track->innerPosition().x(),
track->innerPosition().y(),
track->innerPosition().z());
320 CLHEP::Hep3Vector r3T(
track->outerPosition().x(),
track->outerPosition().y(),
track->outerPosition().z());
323 CLHEP::Hep3Vector p3T(
track->outerMomentum().x(),
track->outerMomentum().y(),
track->outerMomentum().z());
324 CLHEP::Hep3Vector p3_propagated, r3_propagated;
327 cov_propagated *= 1
e-20;
332 std::cout<<
" dump the very first FTS = "<<debug.
dumpFTS(ftsStart)<<std::endl;
337 if(
track->outerPosition().z()>0){
345 std::vector< CSCDetId > refME;
346 for(
int iS=1;iS<5;++iS){
347 for(
int iR=1;iR<4;++iR){
351 else if(4==iS && iR>1){
354 refME.push_back(
CSCDetId(endcap, iS, iR, chamber));
358 for(
size_t iSt = 0; iSt<refME.size();++iSt){
360 std::cout<<
"loop iStatation = "<<iSt<<std::endl;
361 std::cout<<
"refME[iSt]: st = "<<refME[iSt].station()<<
" rg = "<<refME[iSt].ring()<<std::endl;
363 std::map <std::string, bool> chamberTypes;
364 chamberTypes[
"ME11"] =
false;
365 chamberTypes[
"ME12"] =
false;
366 chamberTypes[
"ME13"] =
false;
367 chamberTypes[
"ME21"] =
false;
368 chamberTypes[
"ME22"] =
false;
369 chamberTypes[
"ME31"] =
false;
370 chamberTypes[
"ME32"] =
false;
371 chamberTypes[
"ME41"] =
false;
372 const CSCChamber* cscChamber_base = cscGeom->chamber(refME[iSt].chamberId());
375 std::cout<<
" base iStation : eta = "<<cscGeom->idToDet(detId)->surface().position().eta()<<
" phi = "<<
376 cscGeom->idToDet(detId)->surface().position().phi() <<
" y = " <<cscGeom->idToDet(detId)->surface().position().y()<<std::endl;
381 tSOSDest =
propagate(ftsStart, cscGeom->idToDet(detId)->surface());
385 getFromFTS(ftsStart, p3_propagated, r3_propagated, charge, cov_propagated);
386 float feta = fabs(r3_propagated.eta());
387 float phi = r3_propagated.phi();
391 map<std::string,bool>::iterator iter;
394 for( iter = chamberTypes.begin(); iter != chamberTypes.end(); iter++ ) {
397 if(iter->second && (iterations-1)==int(iSt)){
399 std::cout<<
" Chamber type "<< iter->first<<
" is a candidate..."<<std::endl;
400 std::cout<<
" station() = "<< refME[iSt].station()<<
" ring() = "<<refME[iSt].ring()<<
" iSt = "<<iSt<<std::endl;
402 std::vector <int> coupleOfChambers;
406 for(
size_t iCh =0;iCh<coupleOfChambers.size();++iCh){
408 if (
printalot)
std::cout<<
" Check chamber N = "<<coupleOfChambers.at(iCh)<<std::endl;;
412 [refME[iSt].
ring()-1]
413 [coupleOfChambers.at(iCh)-
FirstCh])){
417 const CSCChamber* cscChamber = cscGeom->chamber(theCSCId.chamberId());
418 const BoundPlane bpCh = cscGeom->idToDet(cscChamber->geographicalId())->surface();
420 float dz = fabs(bpCh.
position().
z() - zFTS);
421 float zDistInner =
track->innerPosition().z() - bpCh.
position().
z();
422 float zDistOuter =
track->outerPosition().z() - bpCh.
position().
z();
427 if(!
isIPdata && (zDistInner*zDistOuter>0. || fabs(zDistInner)<15. || 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());
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;
464 std::cout<<
"Surface to propagate to: pos = "<<cscChamber->layer(iLayer+1)->surface().position()<<
" eta = "
465 <<cscChamber->layer(iLayer+1)->surface().position().eta()<<
" phi = "
466 <<cscChamber->layer(iLayer+1)->surface().position().phi()<<std::endl;
469 tSOSDest =
propagate(ftsInit, cscChamber->layer(iLayer+1)->surface());
473 getFromFTS(ftsInit, p3_propagated, r3_propagated, charge, cov_propagated);
476 if (
printalot)
std::cout<<
"Propagation between layers not successful - notValid TSOS"<<std::endl;
477 extrapolationPassed =
false;
482 if(extrapolationPassed){
483 GlobalPoint theExtrapolationPoint(r3_propagated.x(),r3_propagated.y(),r3_propagated.z());
484 LocalPoint theLocalPoint = cscChamber->layer(iLayer+1)->toLocal(theExtrapolationPoint);
486 inDeadZone = ( inDeadZone ||
488 refME[iSt].station(), refME[iSt].ring()));
490 std::cout<<
" Candidate chamber: extrapolated LocalPoint = "<<theLocalPoint<<
"inDeadZone = "<<inDeadZone<<std::endl;
535 if (
printalot) printf(
"==exit===CSCEfficiency===== run %i\tevent %i\n\n",iRun,iEvent);
edm::InputTag corrlctDigiTag_
edm::InputTag segmentDigiTag_
bool applyTrigger(edm::Handle< edm::TriggerResults > &hltR, const edm::TriggerNames &triggerNames)
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
bool recSimHitEfficiency(CSCDetId &id, FreeTrajectoryState &ftsChamber)
Strings const & triggerNames() const
bool getAbsoluteEfficiency
float caloCompatibility(const reco::Muon &muon)
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
MuonServiceProxy * theService
FreeTrajectoryState * freeState(bool withErrors=true) const
edm::InputTag rechitDigiTag_
DetId geographicalId() const
The label of this GeomDet.
edm::InputTag hlTriggerResults_
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)
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
edm::InputTag clctDigiTag_
bool inSensitiveLocalRegion(double xLocal, double yLocal, int station, int ring)
bool emptyChambers[2][4][4][NumCh]
const Point & position() const
position
edm::InputTag stripDigiTag_
edm::InputTag alctDigiTag_
edm::InputTag wireDigiTag_
const PositionType & position() const
unsigned int minTrackHits
Power< A, B >::type pow(const A &a, const B &b)
TrajectoryStateOnSurface propagate(FreeTrajectoryState &ftsStart, const BoundPlane &bp)