- // Gaudi
- #include "GaudiKernel/PhysicalConstants.h"
- // TODO (hschindl)
- #include "GaudiKernel/IEventProcessor.h"
-
- // Tb/TbEvent
- #include "Event/TbTrack.h"
- #include "Event/TbCluster.h"
-
- // Tb/TbKernel
- #include "TbKernel/TbConstants.h"
- #include "TbKernel/TbModule.h"
-
- // Local
- #include "TbEfficiency.h"
-
- DECLARE_ALGORITHM_FACTORY(TbEfficiency)
-
- //=============================================================================
- // Standard constructor
- //=============================================================================
- TbEfficiency::TbEfficiency(const std::string& name,
- ISvcLocator* pSvcLocator)
- : TbAlgorithm(name, pSvcLocator),
- m_nTracksConsidered(0),
- m_nTracksAssociated(0),
- m_event(0),
- m_pitch(0.055) {
-
- declareProperty("TrackLocation",
- m_trackLocation = LHCb::TbTrackLocation::Default);
- declareProperty("ClusterLocation",
- m_clusterLocation = LHCb::TbClusterLocation::Default);
- declareProperty("DUT", m_dut = 4);
- declareProperty("CheckHitDUT", m_checkHitDUT = true);
- declareProperty("CheckHitAlivePixel", m_checkHitAlivePixel = true);
- declareProperty("PointingResAllowance", m_pointingResAllowance = 0.01);
- declareProperty("PointingResAllowanceDeadPixel", m_pointingResAllowance_deadPixels = 1);
- declareProperty("TakeDeadPixelsFromFile", m_takeDeadPixelsFromFile = true);
- declareProperty("nTotalTracks", m_nTotalTracks = 0);
- declareProperty("MaxChi", m_maxChi = 0);
- }
-
- //=============================================================================
- // Finalization
- //=============================================================================
- StatusCode TbEfficiency::finalize() {
-
- info() <<"DUT Efficiency: "<<m_nTracksAssociated/(1.*m_nTracksConsidered)<<endmsg;
- info()<<"nTracksAssociated: "<<m_nTracksAssociated<<endmsg;
- info()<<"nTracksConsidered: "<<m_nTracksConsidered<<endmsg;
- info()<<"cf: "<<m_eff->GetEfficiency(1)<<"\t ±"<<m_eff->GetEfficiencyErrorLow(1)<<endmsg;
-
- TFile * f = new TFile("EfficiencyResults.root", "RECREATE");
- h_row->Write();
- h_col->Write();
- m_eff->Write();
- h_hitmap->Write();
- h_pixel->Write();
- h_pixel2x2->Write();
-
- m_eff->CreateGraph()->Write();
- h_hitmap->CreateHistogram()->Write();
- h_col->CreateGraph()->Write();
- h_row->CreateGraph()->Write();
- h_pixel->CreateHistogram()->Write();
- h_pixel2x2->CreateHistogram()->Write();
- f->Close();
-
- return TbAlgorithm::finalize();
- }
-
-
- //=============================================================================
- // Initialization
- //=============================================================================
- StatusCode TbEfficiency::initialize() {
- // Initialise the base class.
- StatusCode sc = TbAlgorithm::initialize();
- if (sc.isFailure()) return sc;
-
- uint nBins = 80;
- h_row = new TEfficiency("row", "row", 256, 0, 256);
- h_col = new TEfficiency("col", "col", 3*256+2, 0, 3*256+2);
- m_eff = new TEfficiency("eff", "eff", 1, 0, 1);
- h_hitmap = new TEfficiency("hitmap", "hitmap", 3*256+2, 0, 3*256+2, 256, 0, 256);
- h_pixel = new TEfficiency("pixel", "pixel", nBins, 0, m_pitch, nBins, 0, m_pitch);
- h_pixel2x2 = new TEfficiency("pixel2x2", "pixel2x2", 2*nBins, 0, 2*m_pitch, 2*nBins, 0, 2*m_pitch);
-
- if (m_takeDeadPixelsFromFile) {
- TFile * fIn = new TFile("telPlaneRef.root", "READ");
- const std::string name = "Tb/TbHitMonitor/HitMap/Plane" + std::to_string(m_dut);
- m_deadPixelMap = (TH2F*)fIn->Get(name.c_str());
- info() <<"DeadPixelMap name: "<<m_deadPixelMap->GetName()<<endmsg;
- }
-
- return StatusCode::SUCCESS;
- }
-
- //=============================================================================
- // Main execution
- //=============================================================================
- StatusCode TbEfficiency::execute() {
- // Grab the tracks.
- const LHCb::TbTracks* tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
- if (!tracks) {
- return Error("No tracks in " + m_trackLocation);
- }
-
- // Loop over the tracks.
- for (LHCb::TbTrack* track : *tracks) {
- if (m_nTotalTracks != 0 && m_nTracksConsidered >= m_nTotalTracks) {
- SmartIF<IEventProcessor> app(serviceLocator()->service("ApplicationMgr"));
- if (app) app->stopRun();
- return StatusCode::SUCCESS;
- }
-
- // Find the local position of the track intercept with the DUT.
- const Gaudi::XYZPoint interceptUG = geomSvc()->intercept(track, m_dut);
- const auto interceptUL = geomSvc()->globalToLocal(interceptUG, m_dut);
-
- if (!passSelectionCriteria(track, interceptUL)) continue;
- m_nTracksConsidered++;
- bool pass = false;
- if (track->associatedClusters().size() > 0) pass = true;
-
- int row = interceptUL.y()/m_pitch - 0.5;
- int col = interceptUL.x()/m_pitch - 0.5;
-
- h_pixel->Fill(pass, fmod(interceptUL.x(), m_pitch), fmod(interceptUL.y(), m_pitch));
- h_pixel2x2->Fill(pass, fmod(interceptUL.x(), 2*m_pitch), fmod(interceptUL.y(), 2*m_pitch));
- h_row->Fill(pass, row);
- h_col->Fill(pass, col);
- h_hitmap->Fill(pass, col, row);
- m_eff->Fill(pass, 0);
- if (pass) m_nTracksAssociated++;
- if (!pass && m_event == 1050) info()<<track->htime()<<endmsg;
- }
-
- m_event++;
- return StatusCode::SUCCESS;
- }
-
-
- //=============================================================================
-
- bool TbEfficiency::passSelectionCriteria(LHCb::TbTrack * track, Gaudi::XYZPoint interceptUL)
- {
- if (m_maxChi > 0 && track->chi2PerNdof() > m_maxChi) return false;
- if (m_checkHitDUT && !passedThroughDUT(interceptUL)) {
- return false;
- }
- if (m_checkHitAlivePixel && !passedThroughAlivePixel(interceptUL)) return false;
-
- return true;
- }
-
-
- //=============================================================================
-
- bool TbEfficiency::passedThroughDUT(Gaudi::XYZPoint interceptUL)
- {
- double xLow = 0.0 + m_pointingResAllowance;
- double yLow = 0.0 + m_pointingResAllowance;
- double yUp = 14.08 - m_pointingResAllowance;
- double xUp = 14.08;
-
- if (geomSvc()->modules().at(m_dut)->nChips() == 1)
- xUp = 14.08 - m_pointingResAllowance;
- else if (geomSvc()->modules().at(m_dut)->nChips() == 3)
- xUp = 42.35 - m_pointingResAllowance;
- else error()<<"Unknown device!"<<endmsg;
-
- if (interceptUL.x() < xLow ||
- interceptUL.x() > xUp ||
- interceptUL.y() < yLow ||
- interceptUL.y() > yUp) {
- counter("nTracksNotStrikingDUT") += 1;
- return false;
- }
- return true;
- }
-
-
- //=============================================================================
-
- bool TbEfficiency::passedThroughAlivePixel(Gaudi::XYZPoint interceptUL)
- {
- int row = interceptUL.y()/0.055 - 0.5;
- int col = interceptUL.x()/0.055 - 0.5;
- int one = 1;
- for (int iRow = row - m_pointingResAllowance_deadPixels;
- iRow != row+m_pointingResAllowance_deadPixels+one; iRow++) {
- for (int iCol = col - m_pointingResAllowance_deadPixels;
- iCol != col + m_pointingResAllowance_deadPixels+one; iCol++) {
-
- int colLim;
- if (geomSvc()->modules().at(m_dut)->nChips() == 3) colLim = 3*256+2;
- else colLim = 256;
- if (iRow < 0 || iCol < 0 || iRow >= 256 || iCol >= colLim) continue;
- if (pixelSvc()->isMasked(pixelSvc()->address(iCol, iRow), m_dut)) {
- counter("nTracksStrikingMaskedPixels") += 1;
- return false;
- }
- else if (m_takeDeadPixelsFromFile && m_deadPixelMap->GetBinContent(iCol, iRow) == 0) {
- counter("nTracksStrikingDeadPixels") += 1;
- return false;
- }
-
- }
- }
- return true;
- }
-
-
- //=============================================================================