/
githubmirror
/
cmssw
Обзор
Документация
Войти
/
githubmirror
/
cmssw
Код
Запросы
0
Пакеты
0
Релизы
0
Аналитика
Безопасность
master
RecoLocalMuon/CSCEfficiency/src/CSCEfficiency.cc
2 055 строк
85 KB
Shahzad Malik Muzaffar
[RECONSTRUCTION] [LLVM21] Apply clang-tidy changes
26 ноя 2025, 17:10
26 ноя 2025, 17:10
7863ba0
Код
Авторство
О чём код?
/* * Routine to calculate CSC efficiencies * Comments about the program logic are denoted by //---- * * Stoyan Stoynev, Northwestern University. */ #include "RecoLocalMuon/CSCEfficiency/src/CSCEfficiency.h" #include "TrackingTools/Records/interface/TrackingComponentsRecord.h" #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h" #include "FWCore/Common/interface/TriggerNames.h" #include "DataFormats/MuonReco/interface/Muon.h" #include "DataFormats/MuonReco/interface/MuonFwd.h" #include "DataFormats/MuonReco/interface/MuonSelectors.h" #include "DataFormats/BeamSpot/interface/BeamSpot.h" #include "DataFormats/Provenance/interface/RunLumiEventNumber.h" using namespace std; bool CSCEfficiency::filter(edm::Event &event, const edm::EventSetup &eventSetup) { passTheEvent = false; DataFlow->Fill(0.); MuonPatternRecoDumper debug; //---- increment counter nEventsAnalyzed++; // printalot debug output printalot = (nEventsAnalyzed < int(printout_NEvents)); // edm::RunNumber_t const iRun = event.id().run(); edm::EventNumber_t const iEvent = event.id().event(); if (0 == fmod(double(nEventsAnalyzed), double(1000))) { if (printalot) { printf("\n==enter==CSCEfficiency===== run %u\tevent %llu\tn Analyzed %i\n", iRun, iEvent, nEventsAnalyzed); } } theService->update(eventSetup); //---- These declarations create handles to the types of records that you want //---- to retrieve from event "e". if (printalot) printf("\tget handles for digi collections\n"); //---- Pass the handle to the method "getByType", which is used to retrieve //---- one and only one instance of the type in question out of event "e". If //---- zero or more than one instance exists in the event an exception is thrown. if (printalot) printf("\tpass handles\n"); edm::Handle<CSCALCTDigiCollection> alcts; edm::Handle<CSCCLCTDigiCollection> clcts; edm::Handle<CSCCorrelatedLCTDigiCollection> correlatedlcts; edm::Handle<CSCWireDigiCollection> wires; edm::Handle<CSCStripDigiCollection> strips; edm::Handle<CSCRecHit2DCollection> rechits; edm::Handle<CSCSegmentCollection> segments; edm::Handle<edm::View<reco::Track> > trackCollectionH; edm::Handle<edm::PSimHitContainer> simhits; if (useDigis) { event.getByToken(wd_token, wires); event.getByToken(sd_token, strips); event.getByToken(al_token, alcts); event.getByToken(cl_token, clcts); event.getByToken(co_token, correlatedlcts); } if (!isData) { event.getByToken(sh_token, simhits); } event.getByToken(rh_token, rechits); event.getByToken(se_token, segments); event.getByToken(tk_token, trackCollectionH); const edm::View<reco::Track> trackCollection = *(trackCollectionH.product()); //---- Get the CSC Geometry : if (printalot) printf("\tget the CSC geometry.\n"); edm::ESHandle<CSCGeometry> cscGeom = eventSetup.getHandle(geomToken_); // use theTrackingGeometry instead of cscGeom? //edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry = eventSetup.getHandle(trackingGeomToken_); bool triggerPassed = true; if (useTrigger) { // access the trigger information // trigger names can be find in HLTrigger/Configuration/python/HLT_2E30_cff.py (or?) // get hold of TriggerResults edm::Handle<edm::TriggerResults> hltR; event.getByToken(ht_token, hltR); const edm::TriggerNames &triggerNames = event.triggerNames(*hltR); triggerPassed = applyTrigger(hltR, triggerNames); } if (!triggerPassed) { return triggerPassed; } DataFlow->Fill(1.); GlobalPoint gpZero(0., 0., 0.); if (theService->magneticField()->inTesla(gpZero).mag2() < 0.1) { magField = false; } else { magField = true; } //---- store info from digis fillDigiInfo(alcts, clcts, correlatedlcts, wires, strips, simhits, rechits, segments, cscGeom); // edm::Handle<reco::MuonCollection> muons; edm::InputTag muonTag_("muons"); event.getByLabel(muonTag_, muons); edm::Handle<reco::BeamSpot> beamSpotHandle; event.getByLabel("offlineBeamSpot", beamSpotHandle); reco::BeamSpot vertexBeamSpot = *beamSpotHandle; // std::vector<reco::MuonCollection::const_iterator> goodMuons_it; unsigned int nPositiveZ = 0; unsigned int nNegativeZ = 0; float muonOuterZPosition = -99999.; if (isIPdata) { if (printalot) std::cout << " muons.size() = " << muons->size() << std::endl; for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) { DataFlow->Fill(31.); if (printalot) { std::cout << " iMuon = " << muon - muons->begin() << " charge = " << muon->charge() << " p = " << muon->p() << " pt = " << muon->pt() << " eta = " << muon->eta() << " phi = " << muon->phi() << " matches = " << muon->matches().size() << " matched Seg = " << muon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration) << " GLB/TR/STA = " << muon->isGlobalMuon() << "/" << muon->isTrackerMuon() << "/" << muon->isStandAloneMuon() << std::endl; } if (!(muon->isTrackerMuon() && muon->isGlobalMuon())) { continue; } DataFlow->Fill(32.); double relISO = (muon->isolationR03().sumPt + muon->isolationR03().emEt + muon->isolationR03().hadEt) / muon->track()->pt(); if (printalot) { std::cout << " relISO = " << relISO << " emVetoEt = " << muon->isolationR03().emVetoEt << " caloComp = " << muon::caloCompatibility(*(muon)) << " dxy = " << fabs(muon->track()->dxy(vertexBeamSpot.position())) << std::endl; } if ( //relISO>0.1 || muon::caloCompatibility(*(muon))<.90 || fabs(muon->track()->dxy(vertexBeamSpot.position())) > 0.2 || muon->pt() < 6.) { continue; } DataFlow->Fill(33.); if (muon->track()->hitPattern().numberOfValidPixelHits() < 1 || muon->track()->hitPattern().numberOfValidTrackerHits() < 11 || muon->combinedMuon()->hitPattern().numberOfValidMuonHits() < 1 || muon->combinedMuon()->normalizedChi2() > 10. || muon->numberOfMatches() < 2) { continue; } DataFlow->Fill(34.); float zOuter = muon->combinedMuon()->outerPosition().z(); float rhoOuter = muon->combinedMuon()->outerPosition().rho(); bool passDepth = true; // barrel region //if ( fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 480.){ if (fabs(zOuter) < 660. && rhoOuter > 400. && rhoOuter < 540.) { passDepth = false; } // endcap region //else if( fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.){ else if (fabs(zOuter) > 550. && fabs(zOuter) < 650. && rhoOuter < 300.) { passDepth = false; } // overlap region //else if ( fabs(zOuter) > 680. && fabs(zOuter) < 730. && rhoOuter < 480.){ else if (fabs(zOuter) > 680. && fabs(zOuter) < 880. && rhoOuter < 540.) { passDepth = false; } if (!passDepth) { continue; } DataFlow->Fill(35.); goodMuons_it.push_back(muon); if (muon->track()->momentum().z() > 0.) { ++nPositiveZ; } if (muon->track()->momentum().z() < 0.) { ++nNegativeZ; } } } // if (printalot) std::cout << "Start track loop over " << trackCollection.size() << " tracks" << std::endl; for (edm::View<reco::Track>::size_type i = 0; i < trackCollection.size(); ++i) { DataFlow->Fill(2.); edm::RefToBase<reco::Track> track(trackCollectionH, i); //std::cout<<" iTR = "<<i<<" eta = "<<track->eta()<<" phi = "<<track->phi()<<std::cout<<" pt = "<<track->pt()<<std::endl; if (isIPdata) { if (printalot) { std::cout << " nNegativeZ = " << nNegativeZ << " nPositiveZ = " << nPositiveZ << std::endl; } if (nNegativeZ > 1 || nPositiveZ > 1) { break; } bool trackOK = false; if (printalot) { std::cout << " goodMuons_it.size() = " << goodMuons_it.size() << std::endl; } for (size_t iM = 0; iM < goodMuons_it.size(); ++iM) { //std::cout<<" iM = "<<iM<<" eta = "<<goodMuons_it[iM]->track()->eta()<< //" phi = "<<goodMuons_it[iM]->track()->phi()<< //" pt = "<<goodMuons_it[iM]->track()->pt()<<std::endl; float deltaR = pow(track->phi() - goodMuons_it[iM]->track()->phi(), 2) + pow(track->eta() - goodMuons_it[iM]->track()->eta(), 2); deltaR = sqrt(deltaR); if (printalot) { std::cout << " TR mu match to a tr: deltaR = " << deltaR << " dPt = " << track->pt() - goodMuons_it[iM]->track()->pt() << std::endl; } if (deltaR > 0.01 || fabs(track->pt() - goodMuons_it[iM]->track()->pt()) > 0.1) { continue; } else { trackOK = true; if (printalot) { std::cout << " trackOK " << std::endl; } muonOuterZPosition = goodMuons_it[iM]->combinedMuon()->outerPosition().z(); break; //++nChosenTracks; } } if (!trackOK) { if (printalot) { std::cout << " failed: trackOK " << std::endl; } continue; } } else { //---- Do we need a better "clean track" definition? if (trackCollection.size() > 2) { break; } DataFlow->Fill(3.); if (!i && 2 == trackCollection.size()) { edm::View<reco::Track>::size_type tType = 1; edm::RefToBase<reco::Track> trackTwo(trackCollectionH, tType); if (track->outerPosition().z() * trackTwo->outerPosition().z() > 0) { // in one and the same "endcap" break; } } } DataFlow->Fill(4.); if (printalot) { std::cout << "i track = " << i << " P = " << track->p() << " chi2/ndf = " << track->normalizedChi2() << " nSeg = " << segments->size() << std::endl; std::cout << "quality undef/loose/tight/high/confirmed/goodIt/size " << track->quality(reco::Track::undefQuality) << "/" << track->quality(reco::Track::loose) << "/" << track->quality(reco::Track::tight) << "/" << track->quality(reco::Track::highPurity) << "/" << track->quality(reco::Track::confirmed) << "/" << track->quality(reco::Track::goodIterative) << "/" << track->quality(reco::Track::qualitySize) << std::endl; std::cout << " pt = " << track->pt() << " +-" << track->ptError() << " q/pt = " << track->qoverp() << " +- " << track->qoverpError() << std::endl; //std::cout<<" const Pmin = "<<minTrackMomentum<<" pMax = "<<maxTrackMomentum<<" maxNormChi2 = "<<maxNormChi2<<std::endl; std::cout << " track inner position = " << track->innerPosition() << " outer position = " << track->outerPosition() << std::endl; std::cout << "track eta (outer) = " << track->outerPosition().eta() << " phi (outer) = " << track->outerPosition().phi() << std::endl; if (fabs(track->innerPosition().z()) > 500.) { DetId innerDetId(track->innerDetId()); std::cout << " dump inner state MUON detid = " << debug.dumpMuonId(innerDetId) << std::endl; } if (fabs(track->outerPosition().z()) > 500.) { DetId outerDetId(track->outerDetId()); std::cout << " dump outer state MUON detid = " << debug.dumpMuonId(outerDetId) << std::endl; } std::cout << " nHits = " << track->found() << std::endl; /* trackingRecHit_iterator rhbegin = track->recHitsBegin(); trackingRecHit_iterator rhend = track->recHitsEnd(); int iRH = 0; for(trackingRecHit_iterator recHit = rhbegin; recHit != rhend; ++recHit){ const GeomDet* geomDet = theTrackingGeometry->idToDet((*recHit)->geographicalId()); std::cout<<"hit "<<iRH<<" loc pos = " <<(*recHit)->localPosition()<< " glob pos = " <<geomDet->toGlobal((*recHit)->localPosition())<<std::endl; ++iRH; } */ } float dpT_ov_pT = 0.; if (fabs(track->pt()) > 0.001) { dpT_ov_pT = track->ptError() / track->pt(); } //---- These define a "good" track if (track->normalizedChi2() > maxNormChi2) { // quality break; } DataFlow->Fill(5.); if (track->found() < minTrackHits) { // enough data points break; } DataFlow->Fill(6.); if (!segments->size()) { // better have something in the CSC break; } DataFlow->Fill(7.); if (magField && (track->p() < minP || track->p() > maxP)) { // proper energy range break; } DataFlow->Fill(8.); if (magField && (dpT_ov_pT > 0.5)) { // not too crazy uncertainty break; } DataFlow->Fill(9.); passTheEvent = true; if (printalot) std::cout << "good Track" << std::endl; CLHEP::Hep3Vector r3T_inner(track->innerPosition().x(), track->innerPosition().y(), track->innerPosition().z()); CLHEP::Hep3Vector r3T(track->outerPosition().x(), track->outerPosition().y(), track->outerPosition().z()); chooseDirection(r3T_inner, r3T); // for non-IP CLHEP::Hep3Vector p3T(track->outerMomentum().x(), track->outerMomentum().y(), track->outerMomentum().z()); CLHEP::Hep3Vector p3_propagated, r3_propagated; AlgebraicSymMatrix66 cov_propagated, covT; covT *= 1e-20; cov_propagated *= 1e-20; int charge = track->charge(); FreeTrajectoryState ftsStart = getFromCLHEP(p3T, r3T, charge, covT, &*(theService->magneticField())); if (printalot) { std::cout << " p = " << track->p() << " norm chi2 = " << track->normalizedChi2() << std::endl; std::cout << " dump the very first FTS = " << debug.dumpFTS(ftsStart) << std::endl; } TrajectoryStateOnSurface tSOSDest; int endcap = 0; //---- which endcap to look at if (track->outerPosition().z() > 0) { endcap = 1; } else { endcap = 2; } int chamber = 1; //---- a "refference" CSCDetId for each ring std::vector<CSCDetId> refME; for (int iS = 1; iS < 5; ++iS) { for (int iR = 1; iR < 4; ++iR) { if (1 != iS && iR > 2) { continue; } else if (4 == iS && iR > 1) { continue; } refME.push_back(CSCDetId(endcap, iS, iR, chamber)); } } //---- loop over the "refference" CSCDetIds for (size_t iSt = 0; iSt < refME.size(); ++iSt) { if (printalot) { std::cout << "loop iStatation = " << iSt << std::endl; std::cout << "refME[iSt]: st = " << refME[iSt].station() << " rg = " << refME[iSt].ring() << std::endl; } std::map<std::string, bool> chamberTypes; chamberTypes["ME11"] = false; chamberTypes["ME12"] = false; chamberTypes["ME13"] = false; chamberTypes["ME21"] = false; chamberTypes["ME22"] = false; chamberTypes["ME31"] = false; chamberTypes["ME32"] = false; chamberTypes["ME41"] = false; const CSCChamber *cscChamber_base = cscGeom->chamber(refME[iSt].chamberId()); DetId detId = cscChamber_base->geographicalId(); if (printalot) { std::cout << " base iStation : eta = " << cscGeom->idToDet(detId)->surface().position().eta() << " phi = " << cscGeom->idToDet(detId)->surface().position().phi() << " y = " << cscGeom->idToDet(detId)->surface().position().y() << std::endl; std::cout << " dump base iStation detid = " << debug.dumpMuonId(detId) << std::endl; std::cout << " dump FTS start = " << debug.dumpFTS(ftsStart) << std::endl; } //---- propagate to this ME tSOSDest = propagate(ftsStart, cscGeom->idToDet(detId)->surface()); if (tSOSDest.isValid()) { ftsStart = *tSOSDest.freeState(); if (printalot) std::cout << " dump FTS end = " << debug.dumpFTS(ftsStart) << std::endl; getFromFTS(ftsStart, p3_propagated, r3_propagated, charge, cov_propagated); float feta = fabs(r3_propagated.eta()); float phi = r3_propagated.phi(); //---- which rings are (possibly) penetrated ringCandidates(refME[iSt].station(), feta, chamberTypes); map<std::string, bool>::iterator iter; int iterations = 0; //---- loop over ring candidates for (iter = chamberTypes.begin(); iter != chamberTypes.end(); iter++) { ++iterations; //---- is this ME a machinig candidate station if (iter->second && (iterations - 1) == int(iSt)) { if (printalot) { std::cout << " Chamber type " << iter->first << " is a candidate..." << std::endl; std::cout << " station() = " << refME[iSt].station() << " ring() = " << refME[iSt].ring() << " iSt = " << iSt << std::endl; } std::vector<int> coupleOfChambers; //---- which chamber (and its closes neighbor) is penetrated by the track - candidates chamberCandidates(refME[iSt].station(), refME[iSt].ring(), phi, coupleOfChambers); //---- loop over the two chamber candidates for (size_t iCh = 0; iCh < coupleOfChambers.size(); ++iCh) { DataFlow->Fill(11.); if (printalot) std::cout << " Check chamber N = " << coupleOfChambers.at(iCh) << std::endl; ; if ((!getAbsoluteEfficiency) && (true == emptyChambers[refME[iSt].endcap() - 1][refME[iSt].station() - 1][refME[iSt].ring() - 1] [coupleOfChambers.at(iCh) - FirstCh])) { continue; } CSCDetId theCSCId(refME[iSt].endcap(), refME[iSt].station(), refME[iSt].ring(), coupleOfChambers.at(iCh)); const CSCChamber *cscChamber = cscGeom->chamber(theCSCId.chamberId()); const BoundPlane bpCh = cscGeom->idToDet(cscChamber->geographicalId())->surface(); float zFTS = ftsStart.position().z(); float dz = fabs(bpCh.position().z() - zFTS); float zDistInner = track->innerPosition().z() - bpCh.position().z(); float zDistOuter = track->outerPosition().z() - bpCh.position().z(); //---- only detectors between the inner and outer points of the track are considered for non IP-data if (printalot) { std::cout << " zIn = " << track->innerPosition().z() << " zOut = " << track->outerPosition().z() << " zSurf = " << bpCh.position().z() << std::endl; } if (!isIPdata && (zDistInner * zDistOuter > 0. || fabs(zDistInner) < 15. || fabs(zDistOuter) < 15.)) { // for non IP-data if (printalot) { std::cout << " Not an intermediate (as defined) point... Skip." << std::endl; } continue; } if (isIPdata && fabs(track->eta()) < 1.8) { if (fabs(muonOuterZPosition) - fabs(bpCh.position().z()) < 0 || fabs(muonOuterZPosition - bpCh.position().z()) < 15.) { continue; } } DataFlow->Fill(13.); //---- propagate to the chamber (from this ME) if it is a different surface (odd/even chambers) if (dz > 0.1) { // i.e. non-zero (float 0 check is bad) //if(fabs(zChanmber - zFTS ) > 0.1){ tSOSDest = propagate(ftsStart, cscGeom->idToDet(cscChamber->geographicalId())->surface()); if (tSOSDest.isValid()) { ftsStart = *tSOSDest.freeState(); } else { if (printalot) std::cout << "TSOS not valid! Break." << std::endl; break; } } else { if (printalot) std::cout << " info: dz<0.1" << std::endl; } DataFlow->Fill(15.); FreeTrajectoryState ftsInit = ftsStart; bool inDeadZone = false; //---- loop over the 6 layers for (int iLayer = 0; iLayer < 6; ++iLayer) { bool extrapolationPassed = true; if (printalot) { std::cout << " iLayer = " << iLayer << " dump FTS init = " << debug.dumpFTS(ftsInit) << std::endl; std::cout << " dump detid = " << debug.dumpMuonId(cscChamber->geographicalId()) << std::endl; std::cout << "Surface to propagate to: pos = " << cscChamber->layer(iLayer + 1)->surface().position() << " eta = " << cscChamber->layer(iLayer + 1)->surface().position().eta() << " phi = " << cscChamber->layer(iLayer + 1)->surface().position().phi() << std::endl; } //---- propagate to this layer tSOSDest = propagate(ftsInit, cscChamber->layer(iLayer + 1)->surface()); if (tSOSDest.isValid()) { ftsInit = *tSOSDest.freeState(); if (printalot) std::cout << " Propagation between layers successful: dump FTS end = " << debug.dumpFTS(ftsInit) << std::endl; getFromFTS(ftsInit, p3_propagated, r3_propagated, charge, cov_propagated); } else { if (printalot) std::cout << "Propagation between layers not successful - notValid TSOS" << std::endl; extrapolationPassed = false; inDeadZone = true; } //} //---- Extrapolation passed? For each layer? if (extrapolationPassed) { GlobalPoint theExtrapolationPoint(r3_propagated.x(), r3_propagated.y(), r3_propagated.z()); LocalPoint theLocalPoint = cscChamber->layer(iLayer + 1)->toLocal(theExtrapolationPoint); //std::cout<<" Candidate chamber: extrapolated LocalPoint = "<<theLocalPoint<<std::endl; inDeadZone = (inDeadZone || !inSensitiveLocalRegion( theLocalPoint.x(), theLocalPoint.y(), refME[iSt].station(), refME[iSt].ring())); if (printalot) { std::cout << " Candidate chamber: extrapolated LocalPoint = " << theLocalPoint << "inDeadZone = " << inDeadZone << std::endl; } //---- break if in dead zone for any layer ("clean" tracks) if (inDeadZone) { break; } } else { break; } } DataFlow->Fill(17.); //---- Is a track in a sensitive area for each layer? if (!inDeadZone) { //---- for any layer DataFlow->Fill(19.); if (printalot) std::cout << "Do efficiencies..." << std::endl; //---- Do efficiencies // angle cuts applied (if configured) bool angle_flag = true; angle_flag = efficienciesPerChamber(theCSCId, cscChamber, ftsStart); if (useDigis && angle_flag) { stripWire_Efficiencies(theCSCId, ftsStart); } if (angle_flag) { recHitSegment_Efficiencies(theCSCId, cscChamber, ftsStart); if (!isData) { recSimHitEfficiency(theCSCId, ftsStart); } } } else { if (printalot) std::cout << " Not in active area for all layers" << std::endl; } } if (tSOSDest.isValid()) { ftsStart = *tSOSDest.freeState(); } } } } else { if (printalot) std::cout << " TSOS not valid..." << std::endl; } } } //---- End if (printalot) printf("==exit===CSCEfficiency===== run %u\tevent %llu\n\n", iRun, iEvent); return passTheEvent; } // bool CSCEfficiency::inSensitiveLocalRegion(double xLocal, double yLocal, int station, int ring) { //---- Good region means sensitive area of a chamber. "Local" stands for the local system bool pass = false; std::vector<double> chamberBounds(3); // the sensitive area float y_center = 99999.; //---- hardcoded... not good if (station > 1 && station < 5) { if (2 == ring) { chamberBounds[0] = 66.46 / 2; // (+-)x1 shorter chamberBounds[1] = 127.15 / 2; // (+-)x2 longer chamberBounds[2] = 323.06 / 2; y_center = -0.95; } else { if (2 == station) { chamberBounds[0] = 54.00 / 2; // (+-)x1 shorter chamberBounds[1] = 125.71 / 2; // (+-)x2 longer chamberBounds[2] = 189.66 / 2; y_center = -0.955; } else if (3 == station) { chamberBounds[0] = 61.40 / 2; // (+-)x1 shorter chamberBounds[1] = 125.71 / 2; // (+-)x2 longer chamberBounds[2] = 169.70 / 2; y_center = -0.97; } else if (4 == station) { chamberBounds[0] = 69.01 / 2; // (+-)x1 shorter chamberBounds[1] = 125.65 / 2; // (+-)x2 longer chamberBounds[2] = 149.42 / 2; y_center = -0.94; } } } else if (1 == station) { if (3 == ring) { chamberBounds[0] = 63.40 / 2; // (+-)x1 shorter chamberBounds[1] = 92.10 / 2; // (+-)x2 longer chamberBounds[2] = 164.16 / 2; y_center = -1.075; } else if (2 == ring) { chamberBounds[0] = 51.00 / 2; // (+-)x1 shorter chamberBounds[1] = 83.74 / 2; // (+-)x2 longer chamberBounds[2] = 174.49 / 2; y_center = -0.96; } else { // to be investigated chamberBounds[0] = 30. / 2; //40./2; // (+-)x1 shorter chamberBounds[1] = 60. / 2; //100./2; // (+-)x2 longer chamberBounds[2] = 160. / 2; //142./2; y_center = 0.; } } double yUp = chamberBounds[2] + y_center; double yDown = -chamberBounds[2] + y_center; double xBound1Shifted = chamberBounds[0] - distanceFromDeadZone; // double xBound2Shifted = chamberBounds[1] - distanceFromDeadZone; // double lineSlope = (yUp - yDown) / (xBound2Shifted - xBound1Shifted); double lineConst = yUp - lineSlope * xBound2Shifted; double yBoundary = lineSlope * abs(xLocal) + lineConst; pass = checkLocal(yLocal, yBoundary, station, ring); return pass; } bool CSCEfficiency::checkLocal(double yLocal, double yBoundary, int station, int ring) { //---- check if it is in a good local region (sensitive area - geometrical and HV boundaries excluded) bool pass = false; std::vector<float> deadZoneCenter(6); const float deadZoneHalf = 0.32 * 7 / 2; // wire spacing * (wires missing + 1)/2 float cutZone = deadZoneHalf + distanceFromDeadZone; //cm //---- hardcoded... not good if (station > 1 && station < 5) { if (2 == ring) { deadZoneCenter[0] = -162.48; deadZoneCenter[1] = -81.8744; deadZoneCenter[2] = -21.18165; deadZoneCenter[3] = 39.51105; deadZoneCenter[4] = 100.2939; deadZoneCenter[5] = 160.58; if (yLocal > yBoundary && ((yLocal > deadZoneCenter[0] + cutZone && yLocal < deadZoneCenter[1] - cutZone) || (yLocal > deadZoneCenter[1] + cutZone && yLocal < deadZoneCenter[2] - cutZone) || (yLocal > deadZoneCenter[2] + cutZone && yLocal < deadZoneCenter[3] - cutZone) || (yLocal > deadZoneCenter[3] + cutZone && yLocal < deadZoneCenter[4] - cutZone) || (yLocal > deadZoneCenter[4] + cutZone && yLocal < deadZoneCenter[5] - cutZone))) { pass = true; } } else if (1 == ring) { if (2 == station) { deadZoneCenter[0] = -95.94; deadZoneCenter[1] = -27.47; deadZoneCenter[2] = 33.67; deadZoneCenter[3] = 93.72; } else if (3 == station) { deadZoneCenter[0] = -85.97; deadZoneCenter[1] = -36.21; deadZoneCenter[2] = 23.68; deadZoneCenter[3] = 84.04; } else if (4 == station) { deadZoneCenter[0] = -75.82; deadZoneCenter[1] = -26.14; deadZoneCenter[2] = 23.85; deadZoneCenter[3] = 73.91; } if (yLocal > yBoundary && ((yLocal > deadZoneCenter[0] + cutZone && yLocal < deadZoneCenter[1] - cutZone) || (yLocal > deadZoneCenter[1] + cutZone && yLocal < deadZoneCenter[2] - cutZone) || (yLocal > deadZoneCenter[2] + cutZone && yLocal < deadZoneCenter[3] - cutZone))) { pass = true; } } } else if (1 == station) { if (3 == ring) { deadZoneCenter[0] = -83.155; deadZoneCenter[1] = -22.7401; deadZoneCenter[2] = 27.86665; deadZoneCenter[3] = 81.005; if (yLocal > yBoundary && ((yLocal > deadZoneCenter[0] + cutZone && yLocal < deadZoneCenter[1] - cutZone) || (yLocal > deadZoneCenter[1] + cutZone && yLocal < deadZoneCenter[2] - cutZone) || (yLocal > deadZoneCenter[2] + cutZone && yLocal < deadZoneCenter[3] - cutZone))) { pass = true; } } else if (2 == ring) { deadZoneCenter[0] = -86.285; deadZoneCenter[1] = -32.88305; deadZoneCenter[2] = 32.867423; deadZoneCenter[3] = 88.205; if (yLocal > (yBoundary) && ((yLocal > deadZoneCenter[0] + cutZone && yLocal < deadZoneCenter[1] - cutZone) || (yLocal > deadZoneCenter[1] + cutZone && yLocal < deadZoneCenter[2] - cutZone) || (yLocal > deadZoneCenter[2] + cutZone && yLocal < deadZoneCenter[3] - cutZone))) { pass = true; } } else { deadZoneCenter[0] = -81.0; deadZoneCenter[1] = 81.0; if (yLocal > (yBoundary) && ((yLocal > deadZoneCenter[0] + cutZone && yLocal < deadZoneCenter[1] - cutZone))) { pass = true; } } } return pass; } void CSCEfficiency::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) { for (int iE = 0; iE < 2; iE++) { for (int iS = 0; iS < 4; iS++) { for (int iR = 0; iR < 4; iR++) { for (int iC = 0; iC < NumCh; iC++) { allSegments[iE][iS][iR][iC].clear(); allCLCT[iE][iS][iR][iC] = allALCT[iE][iS][iR][iC] = allCorrLCT[iE][iS][iR][iC] = false; for (int iL = 0; iL < 6; iL++) { allStrips[iE][iS][iR][iC][iL].clear(); allWG[iE][iS][iR][iC][iL].clear(); allRechits[iE][iS][iR][iC][iL].clear(); allSimhits[iE][iS][iR][iC][iL].clear(); } } } } } // if (useDigis) { fillLCT_info(alcts, clcts, correlatedlcts); fillWG_info(wires, cscGeom); fillStrips_info(strips); } fillRechitsSegments_info(rechits, segments, cscGeom); if (!isData) { fillSimhit_info(simhits); } } void CSCEfficiency::fillLCT_info(edm::Handle<CSCALCTDigiCollection> &alcts, edm::Handle<CSCCLCTDigiCollection> &clcts, edm::Handle<CSCCorrelatedLCTDigiCollection> &correlatedlcts) { //---- ALCTDigis int nSize = 0; for (CSCALCTDigiCollection::DigiRangeIterator j = alcts->begin(); j != alcts->end(); j++) { ++nSize; const CSCDetId &id = (*j).first; const CSCALCTDigiCollection::Range &range = (*j).second; for (CSCALCTDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) { // Valid digi in the chamber (or in neighbouring chamber) if ((*digiIt).isValid()) { allALCT[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh] = true; } } // for digis in layer } // end of for (j=... ALCTPerEvent->Fill(nSize); //---- CLCTDigis nSize = 0; for (CSCCLCTDigiCollection::DigiRangeIterator j = clcts->begin(); j != clcts->end(); j++) { ++nSize; const CSCDetId &id = (*j).first; std::vector<CSCCLCTDigi>::const_iterator digiIt = (*j).second.first; std::vector<CSCCLCTDigi>::const_iterator last = (*j).second.second; for (; digiIt != last; ++digiIt) { // Valid digi in the chamber (or in neighbouring chamber) if ((*digiIt).isValid()) { allCLCT[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh] = true; } } } CLCTPerEvent->Fill(nSize); //---- CorrLCTDigis for (CSCCorrelatedLCTDigiCollection::DigiRangeIterator j = correlatedlcts->begin(); j != correlatedlcts->end(); j++) { const CSCDetId &id = (*j).first; std::vector<CSCCorrelatedLCTDigi>::const_iterator digiIt = (*j).second.first; std::vector<CSCCorrelatedLCTDigi>::const_iterator last = (*j).second.second; for (; digiIt != last; ++digiIt) { // Valid digi in the chamber (or in neighbouring chamber) if ((*digiIt).isValid()) { allCorrLCT[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh] = true; } } } } // void CSCEfficiency::fillWG_info(edm::Handle<CSCWireDigiCollection> &wires, edm::ESHandle<CSCGeometry> &cscGeom) { //---- WIRE GROUPS for (CSCWireDigiCollection::DigiRangeIterator j = wires->begin(); j != wires->end(); j++) { CSCDetId id = (CSCDetId)(*j).first; const CSCLayer *layer_p = cscGeom->layer(id); const CSCLayerGeometry *layerGeom = layer_p->geometry(); // std::vector<CSCWireDigi>::const_iterator digiItr = (*j).second.first; std::vector<CSCWireDigi>::const_iterator last = (*j).second.second; // for (; digiItr != last; ++digiItr) { std::pair<int, float> WG_pos(digiItr->getWireGroup(), layerGeom->yOfWireGroup(digiItr->getWireGroup())); std::pair<std::pair<int, float>, int> LayerSignal(WG_pos, digiItr->getTimeBin()); //---- AllWG contains basic information about WG (WG number and Y-position, time bin) allWG[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh][id.layer() - 1].push_back( LayerSignal); if (printalot) { //std::cout<<" WG check : "<<std::endl; //printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n",id.endcap(),id.station(),id.ring(),id.chamber(),id.layer()); //std::cout<<" WG size = "<<allWG[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh] //[id.layer()-1].size()<<std::endl; } } } } void CSCEfficiency::fillStrips_info(edm::Handle<CSCStripDigiCollection> &strips) { //---- STRIPS const float threshold = 13.3; for (CSCStripDigiCollection::DigiRangeIterator j = strips->begin(); j != strips->end(); j++) { CSCDetId id = (CSCDetId)(*j).first; int largestADCValue = -1; std::vector<CSCStripDigi>::const_iterator digiItr = (*j).second.first; std::vector<CSCStripDigi>::const_iterator last = (*j).second.second; for (; digiItr != last; ++digiItr) { int maxADC = largestADCValue; int myStrip = digiItr->getStrip(); std::vector<int> myADCVals = digiItr->getADCCounts(); float thisPedestal = 0.5 * (float)(myADCVals[0] + myADCVals[1]); float peakADC = -1000.; for (int myADCVal : myADCVals) { float diff = (float)myADCVal - thisPedestal; if (diff > threshold) { if (myADCVal > largestADCValue) largestADCValue = myADCVal; if (diff > peakADC) peakADC = diff; } } if (largestADCValue > maxADC) { // FIX IT!!! std::pair<int, float> LayerSignal(myStrip, peakADC); //---- AllStrips contains basic information about strips //---- (strip number and peak signal for most significant strip in the layer) allStrips[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - 1][id.layer() - 1].clear(); allStrips[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - 1][id.layer() - 1].push_back( LayerSignal); } } } } void CSCEfficiency::fillSimhit_info(edm::Handle<edm::PSimHitContainer> &simhits) { //---- SIMHITS edm::PSimHitContainer::const_iterator dSHsimIter; for (dSHsimIter = simhits->begin(); dSHsimIter != simhits->end(); dSHsimIter++) { // Get DetID for this simHit: CSCDetId sId = (CSCDetId)(*dSHsimIter).detUnitId(); std::pair<LocalPoint, int> simHitPos((*dSHsimIter).localPosition(), (*dSHsimIter).particleType()); allSimhits[sId.endcap() - 1][sId.station() - 1][sId.ring() - 1][sId.chamber() - FirstCh][sId.layer() - 1].push_back( simHitPos); } } // void CSCEfficiency::fillRechitsSegments_info(edm::Handle<CSCRecHit2DCollection> &rechits, edm::Handle<CSCSegmentCollection> &segments, edm::ESHandle<CSCGeometry> &cscGeom) { //---- RECHITS AND SEGMENTS //---- Loop over rechits if (printalot) { //printf("\tGet the recHits collection.\t "); printf(" The size of the rechit collection is %i\n", int(rechits->size())); //printf("\t...start loop over rechits...\n"); } recHitsPerEvent->Fill(rechits->size()); //---- Build iterator for rechits and loop : CSCRecHit2DCollection::const_iterator recIt; for (recIt = rechits->begin(); recIt != rechits->end(); recIt++) { //---- Find chamber with rechits in CSC CSCDetId id = (CSCDetId)(*recIt).cscDetId(); if (printalot) { const CSCLayer *csclayer = cscGeom->layer(id); LocalPoint rhitlocal = (*recIt).localPosition(); LocalError rerrlocal = (*recIt).localPositionError(); GlobalPoint rhitglobal = csclayer->toGlobal(rhitlocal); printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n", id.endcap(), id.station(), id.ring(), id.chamber(), id.layer()); printf("\t\tx,y,z: %f, %f, %f\texx,eey,exy: %f, %f, %f\tglobal x,y,z: %f, %f, %f \n", rhitlocal.x(), rhitlocal.y(), rhitlocal.z(), rerrlocal.xx(), rerrlocal.yy(), rerrlocal.xy(), rhitglobal.x(), rhitglobal.y(), rhitglobal.z()); } std::pair<LocalPoint, bool> recHitPos((*recIt).localPosition(), false); allRechits[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh][id.layer() - 1].push_back( recHitPos); } //---- "Empty" chambers for (int iE = 0; iE < 2; iE++) { for (int iS = 0; iS < 4; iS++) { for (int iR = 0; iR < 4; iR++) { for (int iC = 0; iC < NumCh; iC++) { int numLayers = 0; for (int iL = 0; iL < 6; iL++) { if (!allRechits[iE][iS][iR][iC][iL].empty()) { ++numLayers; } } if (numLayers > 1) { emptyChambers[iE][iS][iR][iC] = false; } else { emptyChambers[iE][iS][iR][iC] = true; } } } } } // if (printalot) { printf(" The size of the segment collection is %i\n", int(segments->size())); //printf("\t...start loop over segments...\n"); } segmentsPerEvent->Fill(segments->size()); for (CSCSegmentCollection::const_iterator it = segments->begin(); it != segments->end(); it++) { CSCDetId id = (CSCDetId)(*it).cscDetId(); StHist[id.endcap() - 1][id.station() - 1].segmentChi2_ndf->Fill((*it).chi2() / (*it).degreesOfFreedom()); StHist[id.endcap() - 1][id.station() - 1].hitsInSegment->Fill((*it).nRecHits()); if (printalot) { printf("\tendcap/station/ring/chamber: %i %i %i %i\n", id.endcap(), id.station(), id.ring(), id.chamber()); std::cout << "\tposition(loc) = " << (*it).localPosition() << " error(loc) = " << (*it).localPositionError() << std::endl; std::cout << "\t chi2/ndf = " << (*it).chi2() / (*it).degreesOfFreedom() << " nhits = " << (*it).nRecHits() << std::endl; } allSegments[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh].push_back( make_pair((*it).localPosition(), (*it).localDirection())); //---- try to get the CSC recHits that contribute to this segment. //if (printalot) printf("\tGet the recHits for this segment.\t"); std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits(); int nRH = (*it).nRecHits(); if (printalot) { printf("\tGet the recHits for this segment.\t"); printf(" nRH = %i\n", nRH); } //---- Find which of the rechits in the chamber is in the segment int layerRH = 0; for (vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); iRH++) { ++layerRH; CSCDetId idRH = (CSCDetId)(*iRH).cscDetId(); if (printalot) { printf("\t%i RH\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n", layerRH, idRH.endcap(), idRH.station(), idRH.ring(), idRH.chamber(), idRH.layer()); } for (size_t jRH = 0; jRH < allRechits[idRH.endcap() - 1][idRH.station() - 1][idRH.ring() - 1][idRH.chamber() - FirstCh][idRH.layer() - 1] .size(); ++jRH) { float xDiff = iRH->localPosition().x() - allRechits[idRH.endcap() - 1][idRH.station() - 1][idRH.ring() - 1] [idRH.chamber() - FirstCh][idRH.layer() - 1][jRH] .first.x(); float yDiff = iRH->localPosition().y() - allRechits[idRH.endcap() - 1][idRH.station() - 1][idRH.ring() - 1] [idRH.chamber() - FirstCh][idRH.layer() - 1][jRH] .first.y(); if (fabs(xDiff) < 0.0001 && fabs(yDiff) < 0.0001) { std::pair<LocalPoint, bool> recHitPos(allRechits[idRH.endcap() - 1][idRH.station() - 1][idRH.ring() - 1] [idRH.chamber() - FirstCh][idRH.layer() - 1][jRH] .first, true); allRechits[idRH.endcap() - 1][idRH.station() - 1][idRH.ring() - 1][idRH.chamber() - FirstCh][idRH.layer() - 1] [jRH] = recHitPos; if (printalot) { std::cout << " number of the rechit (from zero) in the segment = " << jRH << std::endl; } } } } } } // void CSCEfficiency::ringCandidates(int station, float feta, std::map<std::string, bool> &chamberTypes) { // yeah, hardcoded again... switch (station) { case 1: if (feta > 0.85 && feta < 1.18) { //ME13 chamberTypes["ME13"] = true; } if (feta > 1.18 && feta < 1.7) { //ME12 chamberTypes["ME12"] = true; } if (feta > 1.5 && feta < 2.45) { //ME11 chamberTypes["ME11"] = true; } break; case 2: if (feta > 0.95 && feta < 1.6) { //ME22 chamberTypes["ME22"] = true; } if (feta > 1.55 && feta < 2.45) { //ME21 chamberTypes["ME21"] = true; } break; case 3: if (feta > 1.08 && feta < 1.72) { //ME32 chamberTypes["ME32"] = true; } if (feta > 1.69 && feta < 2.45) { //ME31 chamberTypes["ME31"] = true; } break; case 4: if (feta > 1.78 && feta < 2.45) { //ME41 chamberTypes["ME41"] = true; } break; default: break; } } // void CSCEfficiency::chamberCandidates(int station, int ring, float phi, std::vector<int> &coupleOfChambers) { coupleOfChambers.clear(); // -pi< phi<+pi float phi_zero = 0.; // check! the phi at the "edge" of Ch 1 float phi_const = 2. * M_PI / 36.; int last_chamber = 36; int first_chamber = 1; if (1 != station && 1 == ring) { // 18 chambers in the ring phi_const *= 2; last_chamber /= 2; } if (phi < 0.) { if (printalot) std::cout << " info: negative phi = " << phi << std::endl; phi += 2 * M_PI; } float chamber_float = (phi - phi_zero) / phi_const; int chamber_int = int(chamber_float); if (chamber_float - float(chamber_int) - 0.5 < 0.) { if (0 != chamber_int) { coupleOfChambers.push_back(chamber_int); } else { coupleOfChambers.push_back(last_chamber); } coupleOfChambers.push_back(chamber_int + 1); } else { coupleOfChambers.push_back(chamber_int + 1); if (last_chamber != chamber_int + 1) { coupleOfChambers.push_back(chamber_int + 2); } else { coupleOfChambers.push_back(first_chamber); } } if (printalot) std::cout << " phi = " << phi << " phi_zero = " << phi_zero << " phi_const = " << phi_const << " candidate chambers: first ch = " << coupleOfChambers[0] << " second ch = " << coupleOfChambers[1] << std::endl; } // bool CSCEfficiency::efficienciesPerChamber(CSCDetId &id, const CSCChamber *cscChamber, FreeTrajectoryState &ftsChamber) { int ec, st, rg, ch, secondRing; returnTypes(id, ec, st, rg, ch, secondRing); LocalVector localDir = cscChamber->toLocal(ftsChamber.momentum()); if (printalot) { std::cout << " global dir = " << ftsChamber.momentum() << std::endl; std::cout << " local dir = " << localDir << std::endl; std::cout << " local theta = " << localDir.theta() << std::endl; } float dxdz = localDir.x() / localDir.z(); float dydz = localDir.y() / localDir.z(); if (2 == st || 3 == st) { if (printalot) { std::cout << "st 3 or 4 ... flip dy/dz" << std::endl; } dydz = -dydz; } if (printalot) { std::cout << "dy/dz = " << dydz << std::endl; } // Apply angle cut bool out = true; if (applyIPangleCuts) { if (dydz > local_DY_DZ_Max || dydz < local_DY_DZ_Min || fabs(dxdz) > local_DX_DZ_Max) { out = false; } } // Segments bool firstCondition = !allSegments[ec][st][rg][ch].empty() ? true : false; bool secondCondition = false; //---- ME1 is special as usual - ME1a and ME1b are actually one chamber if (secondRing > -1) { secondCondition = !allSegments[ec][st][secondRing][ch].empty() ? true : false; } if (firstCondition || secondCondition) { if (out) { ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(1); } } else { if (out) { ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(0); } } if (useDigis) { // ALCTs firstCondition = allALCT[ec][st][rg][ch]; secondCondition = false; if (secondRing > -1) { secondCondition = allALCT[ec][st][secondRing][ch]; } if (firstCondition || secondCondition) { if (out) { ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(3); } // always apply partial angle cuts for this kind of histos if (fabs(dxdz) < local_DX_DZ_Max) { StHist[ec][st].EfficientALCT_momTheta->Fill(ftsChamber.momentum().theta()); ChHist[ec][st][rg][ch].EfficientALCT_dydz->Fill(dydz); } } else { if (out) { ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(2); } if (fabs(dxdz) < local_DX_DZ_Max) { StHist[ec][st].InefficientALCT_momTheta->Fill(ftsChamber.momentum().theta()); ChHist[ec][st][rg][ch].InefficientALCT_dydz->Fill(dydz); } if (printalot) { std::cout << " missing ALCT (dy/dz = " << dydz << ")"; printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n", ec + 1, st + 1, rg + 1, ch + 1); } } // CLCTs firstCondition = allCLCT[ec][st][rg][ch]; secondCondition = false; if (secondRing > -1) { secondCondition = allCLCT[ec][st][secondRing][ch]; } if (firstCondition || secondCondition) { if (out) { ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(5); } if (dydz < local_DY_DZ_Max && dydz > local_DY_DZ_Min) { StHist[ec][st].EfficientCLCT_momPhi->Fill(ftsChamber.momentum().phi()); // - phi chamber... ChHist[ec][st][rg][ch].EfficientCLCT_dxdz->Fill(dxdz); } } else { if (out) { ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(4); } if (dydz < local_DY_DZ_Max && dydz > local_DY_DZ_Min) { StHist[ec][st].InefficientCLCT_momPhi->Fill(ftsChamber.momentum().phi()); // - phi chamber... ChHist[ec][st][rg][ch].InefficientCLCT_dxdz->Fill(dxdz); } if (printalot) { std::cout << " missing CLCT (dx/dz = " << dxdz << ")"; printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n", ec + 1, st + 1, rg + 1, ch + 1); } } if (out) { // CorrLCTs firstCondition = allCorrLCT[ec][st][rg][ch]; secondCondition = false; if (secondRing > -1) { secondCondition = allCorrLCT[ec][st][secondRing][ch]; } if (firstCondition || secondCondition) { ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(7); } else { ChHist[ec][st][rg][ch].digiAppearanceCount->Fill(6); } } } return out; } // bool CSCEfficiency::stripWire_Efficiencies(CSCDetId &id, FreeTrajectoryState &ftsChamber) { int ec, st, rg, ch, secondRing; returnTypes(id, ec, st, rg, ch, secondRing); bool firstCondition, secondCondition; int missingLayers_s = 0; int missingLayers_wg = 0; for (int iLayer = 0; iLayer < 6; iLayer++) { //----Strips if (printalot) { printf("\t%i swEff: \tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n", iLayer + 1, id.endcap(), id.station(), id.ring(), id.chamber(), iLayer + 1); std::cout << " size S = " << allStrips[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh][iLayer].size() << "size W = " << allWG[id.endcap() - 1][id.station() - 1][id.ring() - 1][id.chamber() - FirstCh][iLayer].size() << std::endl; } firstCondition = !allStrips[ec][st][rg][ch][iLayer].empty() ? true : false; //allSegments[ec][st][rg][ch].size() ? true : false; secondCondition = false; if (secondRing > -1) { secondCondition = !allStrips[ec][st][secondRing][ch][iLayer].empty() ? true : false; } if (firstCondition || secondCondition) { ChHist[ec][st][rg][ch].EfficientStrips->Fill(iLayer + 1); } else { if (printalot) { std::cout << "missing strips "; printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n", id.endcap(), id.station(), id.ring(), id.chamber(), iLayer + 1); } } // Wires firstCondition = !allWG[ec][st][rg][ch][iLayer].empty() ? true : false; secondCondition = false; if (secondRing > -1) { secondCondition = !allWG[ec][st][secondRing][ch][iLayer].empty() ? true : false; } if (firstCondition || secondCondition) { ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(iLayer + 1); } else { if (printalot) { std::cout << "missing wires "; printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n", id.endcap(), id.station(), id.ring(), id.chamber(), iLayer + 1); } } } // Normalization if (6 != missingLayers_s) { ChHist[ec][st][rg][ch].EfficientStrips->Fill(8); } if (6 != missingLayers_wg) { ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(8); } ChHist[ec][st][rg][ch].EfficientStrips->Fill(9); ChHist[ec][st][rg][ch].EfficientWireGroups->Fill(9); // ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(1); if (missingLayers_s != missingLayers_wg) { ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(2); if (6 == missingLayers_wg) { ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(3); ChHist[ec][st][rg][ch].NoWires_momTheta->Fill(ftsChamber.momentum().theta()); } if (6 == missingLayers_s) { ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(4); ChHist[ec][st][rg][ch].NoStrips_momPhi->Fill(ftsChamber.momentum().theta()); } } else if (6 == missingLayers_s) { ChHist[ec][st][rg][ch].StripWiresCorrelations->Fill(5); } return true; } // bool CSCEfficiency::recSimHitEfficiency(CSCDetId &id, FreeTrajectoryState &ftsChamber) { int ec, st, rg, ch, secondRing; returnTypes(id, ec, st, rg, ch, secondRing); bool firstCondition, secondCondition; for (int iLayer = 0; iLayer < 6; iLayer++) { firstCondition = !allSimhits[ec][st][rg][ch][iLayer].empty() ? true : false; secondCondition = false; int thisRing = rg; if (secondRing > -1) { secondCondition = !allSimhits[ec][st][secondRing][ch][iLayer].empty() ? true : false; if (secondCondition) { thisRing = secondRing; } } if (firstCondition || secondCondition) { for (size_t iSH = 0; iSH < allSimhits[ec][st][thisRing][ch][iLayer].size(); ++iSH) { if (13 == fabs(allSimhits[ec][st][thisRing][ch][iLayer][iSH].second)) { ChHist[ec][st][rg][ch].SimSimhits->Fill(iLayer + 1); if (!allRechits[ec][st][thisRing][ch][iLayer].empty()) { ChHist[ec][st][rg][ch].SimRechits->Fill(iLayer + 1); } break; } } //---- Next is not too usefull... /* for(unsigned int iSimHits=0; iSimHits<allSimhits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size(); iSimHits++){ ChHist[ec][st][rg][id.chamber()-FirstCh].SimSimhits_each->Fill(iLayer+1); } for(unsigned int iRecHits=0; iRecHits<allRechits[id.endcap()-1][id.station()-1][id.ring()-1][id.chamber()-FirstCh][iLayer].size(); iRecHits++){ ChHist[ec][st][rg][id.chamber()-FirstCh].SimRechits_each->Fill(iLayer+1); } */ // } } return true; } // bool CSCEfficiency::recHitSegment_Efficiencies(CSCDetId &id, const CSCChamber *cscChamber, FreeTrajectoryState &ftsChamber) { int ec, st, rg, ch, secondRing; returnTypes(id, ec, st, rg, ch, secondRing); bool firstCondition, secondCondition; std::vector<bool> missingLayers_rh(6); std::vector<int> usedInSegment(6); // Rechits if (printalot) std::cout << "RecHits eff" << std::endl; for (int iLayer = 0; iLayer < 6; ++iLayer) { firstCondition = !allRechits[ec][st][rg][ch][iLayer].empty() ? true : false; secondCondition = false; int thisRing = rg; if (secondRing > -1) { secondCondition = !allRechits[ec][st][secondRing][ch][iLayer].empty() ? true : false; if (secondCondition) { thisRing = secondRing; } } if (firstCondition || secondCondition) { ChHist[ec][st][rg][ch].EfficientRechits_good->Fill(iLayer + 1); for (size_t iR = 0; iR < allRechits[ec][st][thisRing][ch][iLayer].size(); ++iR) { if (allRechits[ec][st][thisRing][ch][iLayer][iR].second) { usedInSegment[iLayer] = 1; break; } else { usedInSegment[iLayer] = -1; } } } else { missingLayers_rh[iLayer] = true; if (printalot) { std::cout << "missing rechits "; printf("\t\tendcap/station/ring/chamber/layer: %i/%i/%i/%i/%i\n", id.endcap(), id.station(), id.ring(), id.chamber(), iLayer + 1); } } } GlobalVector globalDir; GlobalPoint globalPos; // Segments firstCondition = !allSegments[ec][st][rg][ch].empty() ? true : false; secondCondition = false; int secondSize = 0; int thisRing = rg; if (secondRing > -1) { secondCondition = !allSegments[ec][st][secondRing][ch].empty() ? true : false; secondSize = allSegments[ec][st][secondRing][ch].size(); if (secondCondition) { thisRing = secondRing; } } if (firstCondition || secondCondition) { if (printalot) std::cout << "segments - start ec = " << ec << " st = " << st << " rg = " << rg << " ch = " << ch << std::endl; StHist[ec][st].EfficientSegments_XY->Fill(ftsChamber.position().x(), ftsChamber.position().y()); if (1 == allSegments[ec][st][rg][ch].size() + secondSize) { globalDir = cscChamber->toGlobal(allSegments[ec][st][thisRing][ch][0].second); globalPos = cscChamber->toGlobal(allSegments[ec][st][thisRing][ch][0].first); StHist[ec][st].EfficientSegments_eta->Fill(fabs(ftsChamber.position().eta())); double residual = sqrt(pow(ftsChamber.position().x() - globalPos.x(), 2) + pow(ftsChamber.position().y() - globalPos.y(), 2) + pow(ftsChamber.position().z() - globalPos.z(), 2)); if (printalot) std::cout << " fts.position() = " << ftsChamber.position() << " segPos = " << globalPos << " res = " << residual << std::endl; StHist[ec][st].ResidualSegments->Fill(residual); } for (int iLayer = 0; iLayer < 6; ++iLayer) { if (printalot) std::cout << " iLayer = " << iLayer << " usedInSegment = " << usedInSegment[iLayer] << std::endl; if (0 != usedInSegment[iLayer]) { if (-1 == usedInSegment[iLayer]) { ChHist[ec][st][rg][ch].InefficientSingleHits->Fill(iLayer + 1); } ChHist[ec][st][rg][ch].AllSingleHits->Fill(iLayer + 1); } firstCondition = !allRechits[ec][st][rg][ch][iLayer].empty() ? true : false; secondCondition = false; if (secondRing > -1) { secondCondition = !allRechits[ec][st][secondRing][ch][iLayer].empty() ? true : false; } float stripAngle = 99999.; std::vector<float> posXY(2); bool oneSegment = false; if (1 == allSegments[ec][st][rg][ch].size() + secondSize) { oneSegment = true; const BoundPlane bp = cscChamber->layer(iLayer + 1)->surface(); linearExtrapolation(globalPos, globalDir, bp.position().z(), posXY); GlobalPoint gp_extrapol(posXY.at(0), posXY.at(1), bp.position().z()); const LocalPoint lp_extrapol = cscChamber->layer(iLayer + 1)->toLocal(gp_extrapol); posXY.at(0) = lp_extrapol.x(); posXY.at(1) = lp_extrapol.y(); int nearestStrip = cscChamber->layer(iLayer + 1)->geometry()->nearestStrip(lp_extrapol); stripAngle = cscChamber->layer(iLayer + 1)->geometry()->stripAngle(nearestStrip) - M_PI / 2.; } if (firstCondition || secondCondition) { ChHist[ec][st][rg][ch].EfficientRechits_inSegment->Fill(iLayer + 1); if (oneSegment) { ChHist[ec][st][rg][ch].Y_EfficientRecHits_inSegment[iLayer]->Fill(posXY.at(1)); ChHist[ec][st][rg][ch].Phi_EfficientRecHits_inSegment[iLayer]->Fill(stripAngle); } } else { if (oneSegment) { ChHist[ec][st][rg][ch].Y_InefficientRecHits_inSegment[iLayer]->Fill(posXY.at(1)); ChHist[ec][st][rg][ch].Phi_InefficientRecHits_inSegment[iLayer]->Fill(stripAngle); } } } } else { StHist[ec][st].InefficientSegments_XY->Fill(ftsChamber.position().x(), ftsChamber.position().y()); if (printalot) { std::cout << "missing segment " << std::endl; printf("\t\tendcap/station/ring/chamber: %i/%i/%i/%i\n", id.endcap(), id.station(), id.ring(), id.chamber()); std::cout << " fts.position() = " << ftsChamber.position() << std::endl; } } // Normalization ChHist[ec][st][rg][ch].EfficientRechits_good->Fill(8); if (allSegments[ec][st][rg][ch].size() + secondSize < 2) { StHist[ec][st].AllSegments_eta->Fill(fabs(ftsChamber.position().eta())); } ChHist[ec][st][rg][id.chamber() - FirstCh].EfficientRechits_inSegment->Fill(9); return true; } // void CSCEfficiency::returnTypes(CSCDetId &id, int &ec, int &st, int &rg, int &ch, int &secondRing) { ec = id.endcap() - 1; st = id.station() - 1; rg = id.ring() - 1; secondRing = -1; if (1 == id.station() && (4 == id.ring() || 1 == id.ring())) { rg = 0; secondRing = 3; } ch = id.chamber() - FirstCh; } // void CSCEfficiency::getFromFTS(const FreeTrajectoryState &fts, CLHEP::Hep3Vector &p3, CLHEP::Hep3Vector &r3, int &charge, AlgebraicSymMatrix66 &cov) { GlobalVector p3GV = fts.momentum(); GlobalPoint r3GP = fts.position(); p3.set(p3GV.x(), p3GV.y(), p3GV.z()); r3.set(r3GP.x(), r3GP.y(), r3GP.z()); charge = fts.charge(); cov = fts.hasError() ? fts.cartesianError().matrix() : AlgebraicSymMatrix66(); } FreeTrajectoryState CSCEfficiency::getFromCLHEP(const CLHEP::Hep3Vector &p3, const CLHEP::Hep3Vector &r3, int charge, const AlgebraicSymMatrix66 &cov, const MagneticField *field) { GlobalVector p3GV(p3.x(), p3.y(), p3.z()); GlobalPoint r3GP(r3.x(), r3.y(), r3.z()); GlobalTrajectoryParameters tPars(r3GP, p3GV, charge, field); CartesianTrajectoryError tCov(cov); return cov.kRows == 6 ? FreeTrajectoryState(tPars, tCov) : FreeTrajectoryState(tPars); } void CSCEfficiency::linearExtrapolation(GlobalPoint initialPosition, GlobalVector initialDirection, float zSurface, std::vector<float> &posZY) { double paramLine = lineParameter(initialPosition.z(), zSurface, initialDirection.z()); double xPosition = extrapolate1D(initialPosition.x(), initialDirection.x(), paramLine); double yPosition = extrapolate1D(initialPosition.y(), initialDirection.y(), paramLine); posZY.clear(); posZY.push_back(xPosition); posZY.push_back(yPosition); } // double CSCEfficiency::extrapolate1D(double initPosition, double initDirection, double parameterOfTheLine) { double extrapolatedPosition = initPosition + initDirection * parameterOfTheLine; return extrapolatedPosition; } // double CSCEfficiency::lineParameter(double initZPosition, double destZPosition, double initZDirection) { double paramLine = (destZPosition - initZPosition) / initZDirection; return paramLine; } // void CSCEfficiency::chooseDirection(CLHEP::Hep3Vector &innerPosition, CLHEP::Hep3Vector &outerPosition) { //---- Be careful with trigger conditions too if (!isIPdata) { float dy = outerPosition.y() - innerPosition.y(); float dz = outerPosition.z() - innerPosition.z(); if (isBeamdata) { if (dz > 0) { alongZ = true; } else { alongZ = false; } } else { //cosmics if (dy / dz > 0) { alongZ = false; } else { alongZ = true; } } } } // const Propagator *CSCEfficiency::propagator(std::string propagatorName) const { return &*theService->propagator(propagatorName); } // TrajectoryStateOnSurface CSCEfficiency::propagate(FreeTrajectoryState &ftsStart, const BoundPlane &bpDest) { TrajectoryStateOnSurface tSOSDest; std::string propagatorName; /* // it would work if cosmic muons had properly assigned direction... bool dzPositive = bpDest.position().z() - ftsStart.position().z() > 0 ? true : false; //---- Be careful with trigger conditions too if(!isIPdata){ bool rightDirection = !(alongZ^dzPositive); if(rightDirection){ if(printalot) std::cout<<" propagate along momentum"<<std::endl; propagatorName = "SteppingHelixPropagatorAlong"; } else{ if(printalot) std::cout<<" propagate opposite momentum"<<std::endl; propagatorName = "SteppingHelixPropagatorOpposite"; } } else{ if(printalot) std::cout<<" propagate any (momentum)"<<std::endl; propagatorName = "SteppingHelixPropagatorAny"; } */ propagatorName = "SteppingHelixPropagatorAny"; tSOSDest = propagator(propagatorName)->propagate(ftsStart, bpDest); return tSOSDest; } // bool CSCEfficiency::applyTrigger(edm::Handle<edm::TriggerResults> &hltR, const edm::TriggerNames &triggerNames) { bool triggerPassed = true; const std::vector<std::string> &hlNames = triggerNames.triggerNames(); pointToTriggers.clear(); for (size_t imyT = 0; imyT < myTriggers.size(); ++imyT) { for (size_t iT = 0; iT < hlNames.size(); ++iT) { //std::cout<<" iT = "<<iT<<" hlNames[iT] = "<<hlNames[iT]<< //" : wasrun = "<<hltR->wasrun(iT)<<" accept = "<< // hltR->accept(iT)<<" !error = "<< // !hltR->error(iT)<<std::endl; if (!imyT) { if (hltR->wasrun(iT) && hltR->accept(iT) && !hltR->error(iT)) { TriggersFired->Fill(iT); } } if (hlNames[iT] == myTriggers[imyT]) { pointToTriggers.push_back(iT); if (imyT) { break; } } } } if (pointToTriggers.size() != myTriggers.size()) { pointToTriggers.clear(); if (printalot) { std::cout << " Not all trigger names found - all trigger specifications will be ignored. Check your cfg file!" << std::endl; } } else { if (!pointToTriggers.empty()) { if (printalot) { std::cout << "The following triggers will be required in the event: " << std::endl; for (size_t imyT = 0; imyT < pointToTriggers.size(); ++imyT) { std::cout << " " << hlNames[pointToTriggers[imyT]]; } std::cout << std::endl; std::cout << " in condition (AND/OR) : " << !andOr << "/" << andOr << std::endl; } } } if (hltR.isValid()) { if (pointToTriggers.empty()) { if (printalot) { std::cout << " No triggers specified in the configuration or all ignored - no trigger information will be considered" << std::endl; } } for (size_t imyT = 0; imyT < pointToTriggers.size(); ++imyT) { if (hltR->wasrun(pointToTriggers[imyT]) && hltR->accept(pointToTriggers[imyT]) && !hltR->error(pointToTriggers[imyT])) { triggerPassed = true; if (andOr) { break; } } else { triggerPassed = false; if (!andOr) { triggerPassed = false; break; } } } } else { if (printalot) { std::cout << " TriggerResults handle returns invalid state?! No trigger information will be considered" << std::endl; } } if (printalot) { std::cout << " Trigger passed: " << triggerPassed << std::endl; } return triggerPassed; } // // Constructor CSCEfficiency::CSCEfficiency(const edm::ParameterSet &pset) { // const float Xmin = -70; //const float Xmax = 70; //const int nXbins = int(4.*(Xmax - Xmin)); const float Ymin = -165; const float Ymax = 165; const int nYbins = int((Ymax - Ymin) / 2); const float Layer_min = -0.5; const float Layer_max = 9.5; const int nLayer_bins = int(Layer_max - Layer_min); // //---- Get the input parameters printout_NEvents = pset.getUntrackedParameter<unsigned int>("printout_NEvents", 0); rootFileName = pset.getUntrackedParameter<string>("rootFileName", "cscHists.root"); isData = pset.getUntrackedParameter<bool>("runOnData", true); // isIPdata = pset.getUntrackedParameter<bool>("IPdata", false); // isBeamdata = pset.getUntrackedParameter<bool>("Beamdata", false); // getAbsoluteEfficiency = pset.getUntrackedParameter<bool>("getAbsoluteEfficiency", true); // useDigis = pset.getUntrackedParameter<bool>("useDigis", true); // distanceFromDeadZone = pset.getUntrackedParameter<double>("distanceFromDeadZone", 10.); // minP = pset.getUntrackedParameter<double>("minP", 20.); // maxP = pset.getUntrackedParameter<double>("maxP", 100.); // maxNormChi2 = pset.getUntrackedParameter<double>("maxNormChi2", 3.); // minTrackHits = pset.getUntrackedParameter<unsigned int>("minTrackHits", 10); // applyIPangleCuts = pset.getUntrackedParameter<bool>("applyIPangleCuts", false); // local_DY_DZ_Max = pset.getUntrackedParameter<double>("local_DY_DZ_Max", -0.1); // local_DY_DZ_Min = pset.getUntrackedParameter<double>("local_DY_DZ_Min", -0.8); // local_DX_DZ_Max = pset.getUntrackedParameter<double>("local_DX_DZ_Max", 0.2); // sd_token = consumes<CSCStripDigiCollection>(pset.getParameter<edm::InputTag>("stripDigiTag")); wd_token = consumes<CSCWireDigiCollection>(pset.getParameter<edm::InputTag>("wireDigiTag")); al_token = consumes<CSCALCTDigiCollection>(pset.getParameter<edm::InputTag>("alctDigiTag")); cl_token = consumes<CSCCLCTDigiCollection>(pset.getParameter<edm::InputTag>("clctDigiTag")); co_token = consumes<CSCCorrelatedLCTDigiCollection>(pset.getParameter<edm::InputTag>("corrlctDigiTag")); rh_token = consumes<CSCRecHit2DCollection>(pset.getParameter<edm::InputTag>("rechitTag")); se_token = consumes<CSCSegmentCollection>(pset.getParameter<edm::InputTag>("segmentTag")); tk_token = consumes<edm::View<reco::Track> >(pset.getParameter<edm::InputTag>("tracksTag")); sh_token = consumes<edm::PSimHitContainer>(pset.getParameter<edm::InputTag>("simHitTag")); geomToken_ = esConsumes<CSCGeometry, MuonGeometryRecord>(); edm::ParameterSet serviceParameters = pset.getParameter<edm::ParameterSet>("ServiceParameters"); // maybe use the service for getting magnetic field, propagators, etc. ... theService = new MuonServiceProxy(serviceParameters, consumesCollector()); // Trigger useTrigger = pset.getUntrackedParameter<bool>("useTrigger", false); ht_token = consumes<edm::TriggerResults>(pset.getParameter<edm::InputTag>("HLTriggerResults")); myTriggers = pset.getParameter<std::vector<std::string> >("myTriggers"); andOr = pset.getUntrackedParameter<bool>("andOr"); pointToTriggers.clear(); //---- set counter to zero nEventsAnalyzed = 0; //---- set presence of magnetic field magField = true; // std::string Path = "AllChambers/"; std::string FullName; //---- File with output histograms theFile = new TFile(rootFileName.c_str(), "RECREATE"); theFile->cd(); //---- Book histograms for the analysis char SpecName[60]; sprintf(SpecName, "DataFlow"); DataFlow = new TH1F(SpecName, "Data flow;condition number;entries", 40, -0.5, 39.5); // sprintf(SpecName, "TriggersFired"); TriggersFired = new TH1F(SpecName, "Triggers fired;trigger number;entries", 140, -0.5, 139.5); // int Chan = 50; float minChan = -0.5; float maxChan = 49.5; // sprintf(SpecName, "ALCTPerEvent"); ALCTPerEvent = new TH1F(SpecName, "ALCTs per event;N digis;entries", Chan, minChan, maxChan); // sprintf(SpecName, "CLCTPerEvent"); CLCTPerEvent = new TH1F(SpecName, "CLCTs per event;N digis;entries", Chan, minChan, maxChan); // sprintf(SpecName, "recHitsPerEvent"); recHitsPerEvent = new TH1F(SpecName, "RecHits per event;N digis;entries", 150, -0.5, 149.5); // sprintf(SpecName, "segmentsPerEvent"); segmentsPerEvent = new TH1F(SpecName, "segments per event;N digis;entries", Chan, minChan, maxChan); // //---- Book groups of histograms (for any chamber) map<std::string, bool>::iterator iter; for (int ec = 0; ec < 2; ++ec) { for (int st = 0; st < 4; ++st) { theFile->cd(); sprintf(SpecName, "Stations__E%d_S%d", ec + 1, st + 1); theFile->mkdir(SpecName); theFile->cd(SpecName); // sprintf(SpecName, "segmentChi2_ndf_St%d", st + 1); StHist[ec][st].segmentChi2_ndf = new TH1F(SpecName, "Chi2/ndf of a segment;chi2/ndf;entries", 100, 0., 20.); // sprintf(SpecName, "hitsInSegment_St%d", st + 1); StHist[ec][st].hitsInSegment = new TH1F(SpecName, "Number of hits in a segment;nHits;entries", 7, -0.5, 6.5); // Chan = 170; minChan = 0.85; maxChan = 2.55; // sprintf(SpecName, "AllSegments_eta_St%d", st + 1); StHist[ec][st].AllSegments_eta = new TH1F(SpecName, "All segments in eta;eta;entries", Chan, minChan, maxChan); // sprintf(SpecName, "EfficientSegments_eta_St%d", st + 1); StHist[ec][st].EfficientSegments_eta = new TH1F(SpecName, "Efficient segments in eta;eta;entries", Chan, minChan, maxChan); // sprintf(SpecName, "ResidualSegments_St%d", st + 1); StHist[ec][st].ResidualSegments = new TH1F(SpecName, "Residual (segments);residual,cm;entries", 75, 0., 15.); // Chan = 200; minChan = -800.; maxChan = 800.; int Chan2 = 200; float minChan2 = -800.; float maxChan2 = 800.; sprintf(SpecName, "EfficientSegments_XY_St%d", st + 1); StHist[ec][st].EfficientSegments_XY = new TH2F(SpecName, "Efficient segments in XY;X;Y", Chan, minChan, maxChan, Chan2, minChan2, maxChan2); sprintf(SpecName, "InefficientSegments_XY_St%d", st + 1); StHist[ec][st].InefficientSegments_XY = new TH2F(SpecName, "Inefficient segments in XY;X;Y", Chan, minChan, maxChan, Chan2, minChan2, maxChan2); // Chan = 80; minChan = 0; maxChan = 3.2; sprintf(SpecName, "EfficientALCT_momTheta_St%d", st + 1); StHist[ec][st].EfficientALCT_momTheta = new TH1F(SpecName, "Efficient ALCT in theta (momentum);theta, rad;entries", Chan, minChan, maxChan); // sprintf(SpecName, "InefficientALCT_momTheta_St%d", st + 1); StHist[ec][st].InefficientALCT_momTheta = new TH1F(SpecName, "Inefficient ALCT in theta (momentum);theta, rad;entries", Chan, minChan, maxChan); // Chan = 160; minChan = -3.2; maxChan = 3.2; sprintf(SpecName, "EfficientCLCT_momPhi_St%d", st + 1); StHist[ec][st].EfficientCLCT_momPhi = new TH1F(SpecName, "Efficient CLCT in phi (momentum);phi, rad;entries", Chan, minChan, maxChan); // sprintf(SpecName, "InefficientCLCT_momPhi_St%d", st + 1); StHist[ec][st].InefficientCLCT_momPhi = new TH1F(SpecName, "Inefficient CLCT in phi (momentum);phi, rad;entries", Chan, minChan, maxChan); // theFile->cd(); for (int rg = 0; rg < 3; ++rg) { if (0 != st && rg > 1) { continue; } else if (1 == rg && 3 == st) { continue; } for (int iChamber = FirstCh; iChamber < FirstCh + NumCh; iChamber++) { if (0 != st && 0 == rg && iChamber > 18) { continue; } theFile->cd(); sprintf(SpecName, "Chambers__E%d_S%d_R%d_Chamber_%d", ec + 1, st + 1, rg + 1, iChamber); theFile->mkdir(SpecName); theFile->cd(SpecName); // sprintf(SpecName, "EfficientRechits_inSegment_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].EfficientRechits_inSegment = new TH1F( SpecName, "Existing RecHit given a segment;layers (1-6);entries", nLayer_bins, Layer_min, Layer_max); // sprintf(SpecName, "InefficientSingleHits_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].InefficientSingleHits = new TH1F( SpecName, "Single RecHits not in the segment;layers (1-6);entries ", nLayer_bins, Layer_min, Layer_max); // sprintf(SpecName, "AllSingleHits_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].AllSingleHits = new TH1F( SpecName, "Single RecHits given a segment; layers (1-6);entries", nLayer_bins, Layer_min, Layer_max); // sprintf(SpecName, "digiAppearanceCount_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].digiAppearanceCount = new TH1F(SpecName, "Digi appearance (no-yes): segment(0,1), ALCT(2,3), CLCT(4,5), CorrLCT(6,7); digi type;entries", 8, -0.5, 7.5); // Chan = 100; minChan = -1.1; maxChan = 0.9; sprintf(SpecName, "EfficientALCT_dydz_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].EfficientALCT_dydz = new TH1F(SpecName, "Efficient ALCT; local dy/dz (ME 3 and 4 flipped);entries", Chan, minChan, maxChan); // sprintf(SpecName, "InefficientALCT_dydz_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].InefficientALCT_dydz = new TH1F(SpecName, "Inefficient ALCT; local dy/dz (ME 3 and 4 flipped);entries", Chan, minChan, maxChan); // Chan = 100; minChan = -1.; maxChan = 1.0; sprintf(SpecName, "EfficientCLCT_dxdz_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].EfficientCLCT_dxdz = new TH1F(SpecName, "Efficient CLCT; local dxdz;entries", Chan, minChan, maxChan); // sprintf(SpecName, "InefficientCLCT_dxdz_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].InefficientCLCT_dxdz = new TH1F(SpecName, "Inefficient CLCT; local dxdz;entries", Chan, minChan, maxChan); // sprintf(SpecName, "EfficientRechits_good_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].EfficientRechits_good = new TH1F( SpecName, "Existing RecHit - sensitive area only;layers (1-6);entries", nLayer_bins, Layer_min, Layer_max); // sprintf(SpecName, "EfficientStrips_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].EfficientStrips = new TH1F(SpecName, "Existing strip;layer (1-6); entries", nLayer_bins, Layer_min, Layer_max); // sprintf(SpecName, "EfficientWireGroups_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].EfficientWireGroups = new TH1F(SpecName, "Existing WireGroups;layer (1-6); entries ", nLayer_bins, Layer_min, Layer_max); // sprintf(SpecName, "StripWiresCorrelations_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].StripWiresCorrelations = new TH1F(SpecName, "StripWire correlations;; entries ", 5, 0.5, 5.5); // Chan = 80; minChan = 0; maxChan = 3.2; sprintf(SpecName, "NoWires_momTheta_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].NoWires_momTheta = new TH1F(SpecName, "No wires (all strips present) - in theta (momentum);theta, rad;entries", Chan, minChan, maxChan); // Chan = 160; minChan = -3.2; maxChan = 3.2; sprintf(SpecName, "NoStrips_momPhi_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].NoStrips_momPhi = new TH1F( SpecName, "No strips (all wires present) - in phi (momentum);phi, rad;entries", Chan, minChan, maxChan); // for (int iLayer = 0; iLayer < 6; iLayer++) { sprintf(SpecName, "Y_InefficientRecHits_inSegment_Ch%d_L%d", iChamber, iLayer); ChHist[ec][st][rg][iChamber - FirstCh].Y_InefficientRecHits_inSegment.push_back( new TH1F(SpecName, "Missing RecHit/layer in a segment (local system, whole chamber);Y, cm; entries", nYbins, Ymin, Ymax)); // sprintf(SpecName, "Y_EfficientRecHits_inSegment_Ch%d_L%d", iChamber, iLayer); ChHist[ec][st][rg][iChamber - FirstCh].Y_EfficientRecHits_inSegment.push_back( new TH1F(SpecName, "Efficient (extrapolated from the segment) RecHit/layer in a segment (local system, whole " "chamber);Y, cm; entries", nYbins, Ymin, Ymax)); // Chan = 200; minChan = -0.2; maxChan = 0.2; sprintf(SpecName, "Phi_InefficientRecHits_inSegment_Ch%d_L%d", iChamber, iLayer); ChHist[ec][st][rg][iChamber - FirstCh].Phi_InefficientRecHits_inSegment.push_back( new TH1F(SpecName, "Missing RecHit/layer in a segment (local system, whole chamber);Phi, rad; entries", Chan, minChan, maxChan)); // sprintf(SpecName, "Phi_EfficientRecHits_inSegment_Ch%d_L%d", iChamber, iLayer); ChHist[ec][st][rg][iChamber - FirstCh].Phi_EfficientRecHits_inSegment.push_back( new TH1F(SpecName, "Efficient (extrapolated from the segment) in a segment (local system, whole chamber);Phi, " "rad; entries", Chan, minChan, maxChan)); } // sprintf(SpecName, "Sim_Rechits_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].SimRechits = new TH1F(SpecName, "Existing RecHit (Sim);layers (1-6);entries", nLayer_bins, Layer_min, Layer_max); // sprintf(SpecName, "Sim_Simhits_Ch%d", iChamber); ChHist[ec][st][rg][iChamber - FirstCh].SimSimhits = new TH1F(SpecName, "Existing SimHit (Sim);layers (1-6);entries", nLayer_bins, Layer_min, Layer_max); // /* sprintf(SpecName,"Sim_Rechits_each_Ch%d",iChamber); ChHist[ec][st][rg][iChamber-FirstCh].SimRechits_each = new TH1F(SpecName,"Existing RecHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max); // sprintf(SpecName,"Sim_Simhits_each_Ch%d",iChamber); ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits_each = new TH1F(SpecName,"Existing SimHit (Sim), each;layers (1-6);entries",nLayer_bins,Layer_min,Layer_max); */ theFile->cd(); } } } } } // Destructor CSCEfficiency::~CSCEfficiency() { if (theService) delete theService; // Write the histos to a file theFile->cd(); // char SpecName[60]; std::vector<float> bins, Efficiency, EffError; std::vector<float> eff(2); //---- loop over chambers std::map<std::string, bool> chamberTypes; chamberTypes["ME11"] = false; chamberTypes["ME12"] = false; chamberTypes["ME13"] = false; chamberTypes["ME21"] = false; chamberTypes["ME22"] = false; chamberTypes["ME31"] = false; chamberTypes["ME32"] = false; chamberTypes["ME41"] = false; map<std::string, bool>::iterator iter; std::cout << " Writing proper histogram structure (patience)..." << std::endl; for (int ec = 0; ec < 2; ++ec) { for (int st = 0; st < 4; ++st) { snprintf(SpecName, sizeof(SpecName), "Stations__E%d_S%d", ec + 1, st + 1); theFile->cd(SpecName); StHist[ec][st].segmentChi2_ndf->Write(); StHist[ec][st].hitsInSegment->Write(); StHist[ec][st].AllSegments_eta->Write(); StHist[ec][st].EfficientSegments_eta->Write(); StHist[ec][st].ResidualSegments->Write(); StHist[ec][st].EfficientSegments_XY->Write(); StHist[ec][st].InefficientSegments_XY->Write(); StHist[ec][st].EfficientALCT_momTheta->Write(); StHist[ec][st].InefficientALCT_momTheta->Write(); StHist[ec][st].EfficientCLCT_momPhi->Write(); StHist[ec][st].InefficientCLCT_momPhi->Write(); for (int rg = 0; rg < 3; ++rg) { if (0 != st && rg > 1) { continue; } else if (1 == rg && 3 == st) { continue; } for (int iChamber = FirstCh; iChamber < FirstCh + NumCh; iChamber++) { if (0 != st && 0 == rg && iChamber > 18) { continue; } snprintf(SpecName, sizeof(SpecName), "Chambers__E%d_S%d_R%d_Chamber_%d", ec + 1, st + 1, rg + 1, iChamber); theFile->cd(SpecName); ChHist[ec][st][rg][iChamber - FirstCh].EfficientRechits_inSegment->Write(); ChHist[ec][st][rg][iChamber - FirstCh].AllSingleHits->Write(); ChHist[ec][st][rg][iChamber - FirstCh].digiAppearanceCount->Write(); ChHist[ec][st][rg][iChamber - FirstCh].EfficientALCT_dydz->Write(); ChHist[ec][st][rg][iChamber - FirstCh].InefficientALCT_dydz->Write(); ChHist[ec][st][rg][iChamber - FirstCh].EfficientCLCT_dxdz->Write(); ChHist[ec][st][rg][iChamber - FirstCh].InefficientCLCT_dxdz->Write(); ChHist[ec][st][rg][iChamber - FirstCh].InefficientSingleHits->Write(); ChHist[ec][st][rg][iChamber - FirstCh].EfficientRechits_good->Write(); ChHist[ec][st][rg][iChamber - FirstCh].EfficientStrips->Write(); ChHist[ec][st][rg][iChamber - FirstCh].StripWiresCorrelations->Write(); ChHist[ec][st][rg][iChamber - FirstCh].NoWires_momTheta->Write(); ChHist[ec][st][rg][iChamber - FirstCh].NoStrips_momPhi->Write(); ChHist[ec][st][rg][iChamber - FirstCh].EfficientWireGroups->Write(); for (unsigned int iLayer = 0; iLayer < 6; iLayer++) { ChHist[ec][st][rg][iChamber - FirstCh].Y_InefficientRecHits_inSegment[iLayer]->Write(); ChHist[ec][st][rg][iChamber - FirstCh].Y_EfficientRecHits_inSegment[iLayer]->Write(); ChHist[ec][st][rg][iChamber - FirstCh].Phi_InefficientRecHits_inSegment[iLayer]->Write(); ChHist[ec][st][rg][iChamber - FirstCh].Phi_EfficientRecHits_inSegment[iLayer]->Write(); } ChHist[ec][st][rg][iChamber - FirstCh].SimRechits->Write(); ChHist[ec][st][rg][iChamber - FirstCh].SimSimhits->Write(); /* ChHist[ec][st][rg][iChamber-FirstCh].SimRechits_each->Write(); ChHist[ec][st][rg][iChamber-FirstCh].SimSimhits_each->Write(); */ // theFile->cd(SpecName); theFile->cd(); } } } } // snprintf(SpecName, sizeof(SpecName), "AllChambers"); theFile->mkdir(SpecName); theFile->cd(SpecName); DataFlow->Write(); TriggersFired->Write(); ALCTPerEvent->Write(); CLCTPerEvent->Write(); recHitsPerEvent->Write(); segmentsPerEvent->Write(); // theFile->cd(SpecName); //---- Close the file theFile->Close(); } DEFINE_FWK_MODULE(CSCEfficiency);