- #include <fstream>
-
- // ROOT
- #include "TMath.h"
-
- // Gaudi
- #include "GaudiKernel/PhysicalConstants.h"
- #include "GaudiUtils/HistoLabels.h"
-
- // Tb/TbKernel
- #include "TbKernel/TbConstants.h"
- #include "TbKernel/TbModule.h"
-
- // Local
- #include "TbEffPur.h"
-
- using namespace Gaudi::Utils::Histos;
-
- DECLARE_ALGORITHM_FACTORY(TbEffPur)
-
- //=============================================================================
- // Standard constructor
- //=============================================================================
- TbEffPur::TbEffPur(const std::string& name, ISvcLocator* pSvcLocator)
- : TbAlgorithm(name, pSvcLocator),
- m_pitch(0.055),
- m_nDUTpixels(256),
- m_nTracks(0),
- m_nClusters(0),
- m_nTrackedClusters(0),
- m_nClustersPassedCentral(0),
- m_nTracksCentral(0),
- m_nClustersPassedCorner(0),
- m_nTracksCorner(0),
- m_eff(0.),
- m_pur(0.),
- m_deadAreaRadius(2.),
- m_trackAssociated(NULL)
- {
- declareProperty("TrackLocation",
- m_trackLocation = LHCb::TbTrackLocation::Default);
- declareProperty("VertexLocation",
- m_vertexLocation = LHCb::TbVertexLocation::Default);
- declareProperty("ClusterLocation",
- m_clusterLocation = LHCb::TbClusterLocation::Default);
-
- // Parameters.
- declareProperty("DUTindex", m_DUTindex = 4);
- declareProperty("GlobalCutXLow", m_xLow = 3);
- declareProperty("GlobalCutXUp", m_xUp = 9);
- declareProperty("GlobalCutYLow", m_yLow = 3);
- declareProperty("GlobalCutYUp", m_yUp = 9);
- declareProperty("LocalProbabilityCut", m_probCut = 0.5);
-
- declareProperty("RResidualCut", m_rResidualCut = 0.1);
- declareProperty("TResidualCut", m_tResidualCut = 100);
- declareProperty("ChargeCutLow", m_chargeCutLow = 0);
- // ... note this global cut enforced at different point to above.
- declareProperty("ChargeCutUp", m_chargeCutUp = 1000000);
- declareProperty("LitSquareSide", m_litSquareSide = 0);
-
- declareProperty("ViewerOutput", m_viewerOutput = false);
- declareProperty("ViewerEventNum", m_viewerEvent = 100);
- declareProperty("TGap", m_tGap = 200);
- declareProperty("CorrelationTimeWindow", m_correlationTimeWindow = 100);
- declareProperty("ApplyVeto", m_applyVeto = true);
- declareProperty("TelescopeClusterVetoDelT", m_telescopeClusterVetoDelT = 30);
- declareProperty("EdgeVetoDistance", m_edgeVetoDistance = 0.025);
- m_nEvent = -1;
- }
-
- //=============================================================================
- // Initialization
- //=============================================================================
- StatusCode TbEffPur::initialize() {
-
- // Initialise the base class.
- StatusCode sc = TbAlgorithm::initialize();
- if (sc.isFailure()) return sc;
-
- // Efficiency objects.
- m_effX = new TEfficiency("efficiencyX", "efficiencyX", 40, 0.0, 0.11);
- m_effY = new TEfficiency("efficiencyY", "efficiencyY", 40, 0.0, 0.11);
- m_effs = new TEfficiency("efficiency", "efficiency", 1, -0.5, 0.5);
- m_purs = new TEfficiency("purity", "Tb/TbEffPur/pur", 1, -0.5, 0.5);
- m_effHitmap = new TEfficiency("efficiencyHitmap", "efficiencyHitmap",
- 64, 0.0, 14.08, 64, 0.0, 14.08);
- m_purHitmap = new TEfficiency("puritHitmap", "purityHitmap",
- 64, 0.0, 14.08, 64, 0.0, 14.08);
-
- int nBins = 50;
- m_effHitmapInterPixel = new TEfficiency("efficiencyHitmapInterPixel", "efficiencyHitmapInterPixel",
- nBins, 0.0, 0.055, nBins, 0.0, 0.055);
- m_purHitmapInterPixel = new TEfficiency("puritHitmapInterPixel", "purityHitmapInterPixel",
- nBins, 0.0, 0.055, nBins, 0.0, 0.055);
-
- m_effHitmapInterPixelTriple = new TEfficiency("efficiencyHitmapInterPixelTriple", "efficiencyHitmapInterPixelTriple",
- 150, 0.0, 10*0.055, 20, 0.0, 2*0.055);
-
- for (unsigned int i=1; i<5; i++) {
- std::string name = "efficiencyHitmapInterClusterSize" + std::to_string(i);
- TEfficiency * e = new TEfficiency(name.c_str(), name.c_str(),
- nBins, 0.0, 0.055, nBins, 0.0, 0.055);
- m_effHitmapInterPixelVsSizes.push_back(e);
- }
-
- // Plots for post correlations.
- m_remainsCorrelationsX = book2D("remainsClusterCorrelationsX", "remainsClusterCorrelationsX", 0, 14, 200, 0, 14, 200);
- m_remainsCorrelationsY = book2D("remainsClusterCorrelationsY", "remainsClusterCorrelationsY", 0, 14, 200, 0, 14, 200);
- m_remainsDifferencesXY = book2D("remainsDifferencesXY", "remainsDifferencesXY", -1.5, 1.5, 150, -1.5, 1.5, 150);
- m_clusterRemainsPositionsLocal = book2D("clusterRemainsPositionsLocal", "clusterRemainsPositionsLocal", -5, 20, 600, -5, 20, 600);
- m_trackRemainsPositionsLocal = book2D("trackRemainsPositionsLocal", "trackRemainsPositionsLocal", 0, 256, 600, 0, 256, 600);
- m_clusterRemainsPositionsGlobal = book2D("clusterRemainsPositionsGlobal", "clusterRemainsPositionsGlobal", -5, 20, 600, -5, 20, 600);
- m_trackRemainsPositionsGlobal = book2D("trackRemainsPositionsGlobal", "trackRemainsPositionsGlobal", -5, 20, 600, -5, 20, 600);
- m_vetoTracksHitmap = book2D("vetoTrackHitmap", "vetoTrackHitmap", -5, 20, 600, -5, 20, 600);
- m_vetoClustersHitmap = book2D("vetoClusterHitmap", "vetoClusterHitmap", -5, 20, 600, -5, 20, 600);
- m_timeResidualVsColumn = book2D("timeResidualVsColumn", "timeResidualVsColumn", 0, 256, 256, -50, 50, 500);
-
-
- // Setup the tools.
- m_trackFit = tool<ITbTrackFit>("TbTrackFit", "Fitter", this);
- // h_effVsCharge = new TH2F("effVsCharge", "effVsCharge", 50, 0.5, 1, 20, 0.5, 20.5);
- // TFile * f = new TFile("/afs/cern.ch/user/d/dsaunder/cmtuser/KEPLER/KEPLER_HEAD/EDR_plots/EfficiencyResults/Eff-S22-Run6309-500V.root", "READ");
- // h_effInter = (TH2F*)f->Get("eff_histo");
- return StatusCode::SUCCESS;
-
- }
-
-
- //=============================================================================
- // Finalizer
- //=============================================================================
- StatusCode TbEffPur::finalize() {
- m_effs->SetTotalEvents(0, m_nTracks);
- m_effs->SetPassedEvents(0, m_nTrackedClusters);
-
- m_purs->SetTotalEvents(0, m_nClusters);
- m_purs->SetPassedEvents(0, m_nTrackedClusters);
-
- m_eff = m_nTrackedClusters/(1.0*m_nTracks);
- m_pur = m_nTrackedClusters/(1.0*m_nClusters);
- std::cout<<"ith ROOT Efficiency: " << m_effs->GetEfficiency(0) << " + " << m_effs->GetEfficiencyErrorUp(0) <<" - "<< m_effs->GetEfficiencyErrorLow(0) <<std::endl;
- std::cout<<"ith ROOT Purity: " << m_purs->GetEfficiency(0) << " + " << m_purs->GetEfficiencyErrorUp(0) <<" - "<< m_purs->GetEfficiencyErrorLow(0) <<std::endl;
-
- std::cout<<"Corner efficiency: "<<m_nClustersPassedCorner/(1.*m_nTracksCorner)<<std::endl;
- std::cout<<"Central efficiency: "<<m_nClustersPassedCentral/(1.*m_nTracksCentral)<<std::endl;
-
- plot(m_eff, "Efficiency", "Efficiency", 0.0, 1, 1);
- plot(m_pur, "Purity", "Purity", 0.0, 1, 1);
- plot(m_effs->GetEfficiencyErrorUp(0), "EfficiencyError", "EfficiencyError", 0.0, 1, 1);
- plot(m_purs->GetEfficiencyErrorUp(0), "PurityError", "PurityError", 0.0, 1, 1);
-
- TH2D * h = Gaudi::Utils::Aida2ROOT::aida2root(m_trackRemainsPositionsLocal);
- double maxBin = h->GetBinContent(h->GetMaximumBin());
- double minBin = h->GetBinContent(h->GetMinimumBin());
- for (int i=0; i<m_trackRemainsPositionsLocal->xAxis().bins(); i++) {
- for (int j=0; j<m_trackRemainsPositionsLocal->yAxis().bins(); j++) {
- plot(m_trackRemainsPositionsLocal->binHeight(i, j), "trackRemainsBins",
- "trackRemainsBins", minBin, maxBin, int(maxBin-minBin));
- }
- }
-
- TFile * f = new TFile("EfficiencyResults.root", "RECREATE");
- //m_effs->Write();
-
- m_effHitmapInterPixel->Write();
- m_effHitmapInterPixelTriple->Write();
- //h_effVsCharge->Write();
- //m_effs->CreateHistogram()->Write();
- m_effHitmapInterPixel->CreateHistogram()->Write();
- m_effHitmapInterPixelTriple->CreateHistogram()->Write();
-
-
-
- // for (unsigned int i=0; i<m_effHitmapInterPixelVsSizes.size(); i++) {
- // m_effHitmapInterPixelVsSizes[i]->Write();
- // //m_effHitmapInterPixelVsSizes[i]->CreateHistogram()->Write();
- // }
- f->Close();
- return TbAlgorithm::finalize();
- }
-
-
- //=============================================================================
- // Main execution
- //=============================================================================
- StatusCode TbEffPur::execute()
- {
- m_nEvent++;
-
- m_tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
- //m_vertices = getIfExists<LHCb::TbVertices>(m_vertexLocation);
- m_clusters = getIfExists<LHCb::TbClusters>(m_clusterLocation + std::to_string(m_DUTindex));
-
- if (m_clusters->size() < 100 || m_tracks->size() < 100) return StatusCode::SUCCESS;
-
- // Matching.
- effPur();
- return StatusCode::SUCCESS;
- }
-
-
- //=============================================================================
- // Efficiency finding
- //=============================================================================
- void TbEffPur::effPur()
- {
- // Non-veto'd containers.
- std::vector<LHCb::TbCluster*> * cutClusters = new std::vector<LHCb::TbCluster*>();
- std::vector<LHCb::TbTrack*> cutTracks;
- int nClustersPreVeto = m_clusters->size();
- int nTracksPreVeto = m_tracks->size();
-
-
- // Apply veto.
- if (m_applyVeto) applyVeto(cutClusters, &cutTracks);
- else fillAllTrackClusters(cutClusters, &cutTracks);
- m_trackAssociated = new std::vector<bool>(cutTracks.size(), false);
- //std::cout<<"Fraction tracks not veto'd: "<<cutTracks.size()/(1.0*nTracksPreVeto)<<"\tEvent: "<<m_nEvent<<std::endl;
- //std::cout<<"Fraction clusters not veto'd: "<<cutClusters->size()/(1.0*nClustersPreVeto)<<"\tEvent: "<<m_nEvent<<std::endl;
- plot(cutClusters->size()/(1.0*nClustersPreVeto), "clusterFractionVetod", "clusterFractionVetod", 0.0, 1.0, 200);
- plot(cutTracks.size()/(1.0*nTracksPreVeto), "trackFractionVetod", "trackFractionVetod", 0.0, 1.0, 200);
-
-
- // Do matching.
- m_nClusters += cutClusters->size();
- m_nTracks += cutTracks.size();
- trackClusters(cutClusters, &cutTracks);
-
-
- // Viewer options.
- if (m_nEvent == m_viewerEvent) {
- if (m_viewerOutput) outputViewerData();
- for (LHCb::TbClusters::iterator iclust = cutClusters->begin();
- iclust != cutClusters->end(); iclust++) {
- // if (!(*iclust)->associated())
- // std::cout<<"Note: non-associated clusters at hTime: "<<(*iclust)->htime()<<std::endl;
- }
- }
-
- // Correlate remains.
- //correlateRemains(&cutTracks, cutClusters);
- delete cutClusters;
- delete m_trackAssociated;
- }
-
-
- //=============================================================================
- // Global cuts
- //=============================================================================
- void TbEffPur::applyVeto(std::vector<LHCb::TbCluster*> * cutClusters,
- std::vector<LHCb::TbTrack*> * cutTracks)
- {
- // Gather all cluster times (over all planes) and sort it ___________________
- std::vector<double> allClusterTimes;
- for (unsigned int i=0; i<m_nPlanes; i++) {
- LHCb::TbClusters * clusters = getIfExists<LHCb::TbClusters>(m_clusterLocation + std::to_string(i));
- for (LHCb::TbClusters::const_iterator it = clusters->begin(); it != clusters->end(); ++it)
- allClusterTimes.push_back((*it)->htime());
- }
- std::sort(allClusterTimes.begin(), allClusterTimes.end());
-
-
-
- // Get big enough gaps in pairs _____________________________________________
- std::vector<double> tGapCenters;
- tGapCenters.push_back(0);
- for (std::vector<double>::iterator itime = allClusterTimes.begin()+1;
- itime != allClusterTimes.end(); itime++) {
- plot((*itime) - (*(itime-1)), "allGapWidths", "allGapWidths", 0.0, 3000.0, 200);
- if ((*itime) - (*(itime-1)) > m_tGap) {
- tGapCenters.push_back(0.5*((*itime) + (*(itime-1))));
- plot((*itime) - (*(itime-1)), "bigGapWidths", "bigGapWidths", 0.0, 3000.0, 200);
- }
- }
- tGapCenters.push_back(allClusterTimes.back());
- counter("nGaps") += tGapCenters.size();
- plot(tGapCenters.size(), "nGaps", "nGaps", 0.0, 1000.0, 200);
-
-
-
- // Loop over these sub events and veto ______________________________________
- std::vector<double>::iterator igap;
- for (igap = tGapCenters.begin(); igap != tGapCenters.end() - 1; igap++) {
- bool veto = false;
-
- if (igap == tGapCenters.begin() || igap == tGapCenters.end() -1) veto = true;
-
- // Loop tracks.
- if (!veto) {
- for (LHCb::TbTrack* track : *m_tracks) {
- if (track->htime() > (*igap) && track->htime() < (*(igap+1))) {
- // Veto on tracks outside cuts.
- const Gaudi::XYZPoint trackInterceptGlobal = geomSvc()->intercept(track, m_DUTindex);
- if (!globalCutPosition(trackInterceptGlobal)) veto = true;
- else if (outsideDUT(trackInterceptGlobal)) veto = true;
- else if (interceptDeadPixel(trackInterceptGlobal)) veto = true;
- }
- }
- }
-
-
- // Loop vertices.
- // if (!veto) {
- // for (LHCb::TbVertex* vertex : *m_vertices) {
- // if (vertex->htime() > (*igap) && vertex->htime() < (*(igap+1))) {
- // veto = true;
- // }
- // }
- // }
-
-
- // Loop DUT clusters.
- if (!veto) {
- for (const LHCb::TbCluster* cluster : *m_clusters) {
- if (cluster->htime() > (*igap) && cluster->htime() < (*(igap+1))) {
- // Veto on clusters outside cuts.
- Gaudi::XYZPoint clusterInterceptGlobal(cluster->x(), cluster->y(), cluster->z());
- if (!globalCutPosition(clusterInterceptGlobal)) veto = true;
- if (cluster->charge() > m_chargeCutUp || cluster->charge() < m_chargeCutLow) veto = true;
- }
- }
- }
-
-
- // Loop telescope clusters.
- if (!veto) {
- for (unsigned int i=0; i<m_nPlanes; i++) {
- if (i == m_DUTindex) continue;
- LHCb::TbClusters * clusters = getIfExists<LHCb::TbClusters>(m_clusterLocation + std::to_string(i));
- for (LHCb::TbClusters::const_iterator it = clusters->begin(); it != clusters->end()-1; ++it) {
- if ((*it)->htime() > (*igap) && (*it)->htime() < (*(igap+1))) {
- double delt = (*it)->htime() - (*(it+1))->htime();
- if (fabs(delt) < m_telescopeClusterVetoDelT) {
- veto = true;
- break;
- }
- }
- }
- }
- }
-
-
- if (!veto) fillTrackClusters(cutClusters, cutTracks, (*igap), (*(igap+1)));
- else counter("nGapsFiltered") += 1;
- }
- }
-
-
-
- //=============================================================================
- // Correlate remains
- //=============================================================================
- void TbEffPur::correlateRemains(std::vector<LHCb::TbTrack*> * cutTracks,
- std::vector<LHCb::TbCluster*> * cutClusters)
- {
- for (std::vector<LHCb::TbCluster*>::iterator ic = cutClusters->begin();
- ic != cutClusters->end(); ++ic) {
- if (!(*ic)->associated()) {
- m_clusterRemainsPositionsGlobal->fill((*ic)->x(), (*ic)->y());
- m_clusterRemainsPositionsLocal->fill((*ic)->xloc(), (*ic)->yloc());
- plot((*ic)->charge(), "chargeClusterRemains", "chargeClusterRemains", 0.0, 1000.0, 125);
- }
- }
-
- unsigned int i = 0;
- for (std::vector<LHCb::TbTrack*>::iterator itrack = cutTracks->begin();
- itrack != cutTracks->end(); itrack++) {
- if (!m_trackAssociated->at(i)) {
- Gaudi::XYZPoint trackIntercept = geomSvc()->intercept((*itrack), m_DUTindex);
- m_trackRemainsPositionsGlobal->fill(trackIntercept.x(), trackIntercept.y());
- auto interceptUL = geomSvc()->globalToLocal(trackIntercept, m_DUTindex);
- m_trackRemainsPositionsLocal->fill(interceptUL.x()/m_pitch - 0.5, interceptUL.y()/m_pitch - 0.5);
- }
- i++;
- }
-
-
- // Draws correlations plots between non-associated clusters and tracks.
- for (std::vector<LHCb::TbTrack*>::iterator itrack = cutTracks->begin();
- itrack != cutTracks->end(); itrack++) {
- for (std::vector<LHCb::TbCluster*>::iterator ic = cutClusters->begin();
- ic != cutClusters->end(); ++ic) {
- if ((*ic)->htime() < (*itrack)->htime() - m_correlationTimeWindow) continue;
- if (!(*ic)->associated()) {
- //if (m_nEvent == m_viewerEvent) std::cout<<"Impurity at time: "<<(*ic)->htime()<<std::endl;
- Gaudi::XYZPoint trackIntercept = geomSvc()->intercept((*itrack), m_DUTindex);
- m_remainsCorrelationsX->fill(trackIntercept.x(), (*ic)->x());
- m_remainsCorrelationsY->fill(trackIntercept.y(), (*ic)->y());
-
- plot(trackIntercept.x() - (*ic)->x(), "remainClustDifferencesX", "remainsClustDifferencesX", -1.5, 1.5, 150);
- plot(trackIntercept.y() - (*ic)->y(), "remainClustDifferencesY", "remainClustDifferencesY", -1.5, 1.5, 150);
- plot((*itrack)->htime() - (*ic)->htime(), "remainClustDifferencesT", "remainClustDifferencesT", -50, 50, 150);
-
- m_remainsDifferencesXY->fill(trackIntercept.y() - (*ic)->y(), trackIntercept.x() - (*ic)->x());
- }
- if ((*ic)->htime() > (*itrack)->htime() + m_correlationTimeWindow) break;
- }
- }
- }
-
-
- //=============================================================================
- //
- //=============================================================================
- void TbEffPur::fillAllTrackClusters(std::vector<LHCb::TbCluster*> * cutClusters,
- std::vector<LHCb::TbTrack*> * cutTracks)
- {
- for (LHCb::TbTrack* track : *m_tracks) cutTracks->push_back(track);
- for (std::vector<LHCb::TbCluster*>::iterator iClust = m_clusters->begin();
- iClust != m_clusters->end(); iClust++) cutClusters->push_back(*iClust);
- }
-
-
- //=============================================================================
- //
- //=============================================================================
- void TbEffPur::fillTrackClusters(std::vector<LHCb::TbCluster*> * cutClusters,
- std::vector<LHCb::TbTrack*> * cutTracks, double tlow, double tup)
- {
- // Push tracks and clusters inside this time window, which passed the veto.
- for (LHCb::TbTrack* track : *m_tracks) {
- if (track->htime() > tlow && track->htime() < tup) {
- Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
- cutTracks->push_back(track);
- m_vetoTracksHitmap->fill(trackIntercept.x(), trackIntercept.y());
-
- Gaudi::XYZPoint trackInterceptLocal = geomSvc()->globalToLocal(trackIntercept, m_DUTindex);
- int row = trackInterceptLocal.y()/m_pitch - 0.5;
- int col = trackInterceptLocal.x()/m_pitch - 0.5;
-
- if (pixelSvc()->isMasked(pixelSvc()->address(col, row), m_DUTindex)) {
- std::cout<<"Shouldn't be here!"<<std::endl;
- }
- }
- }
-
- for (LHCb::TbCluster* cluster : *m_clusters) {
- if (cluster->htime() > tlow && cluster->htime() < tup) {
- cutClusters->push_back(cluster);
- m_vetoClustersHitmap->fill(cluster->x(), cluster->y());
- }
- }
- }
-
-
- //=============================================================================
- // Global cuts
- //=============================================================================
- void TbEffPur::trackClusters(std::vector<LHCb::TbCluster*> * cutClusters,
- std::vector<LHCb::TbTrack*> * cutTracks)
- {
- for (std::vector<LHCb::TbTrack*>::iterator itrack = cutTracks->begin();
- itrack != cutTracks->end(); itrack++) {
-
- // Get local position of track.
- const auto interceptUG = geomSvc()->intercept((*itrack), m_DUTindex);
- const auto interceptUL = geomSvc()->globalToLocal(interceptUG, m_DUTindex);
-
-
- // Loop over clusters to find the closest.
- LHCb::TbCluster * closestCluster = NULL;
- double bestRadialDistance;
- for (std::vector<LHCb::TbCluster*>::iterator iclust = cutClusters->begin();
- iclust != cutClusters->end(); iclust++) {
- bool match = matchTrackToCluster((*iclust), (*itrack));
- double radialDistance = getRadialSeparation((*iclust), (*itrack));
- if (match) {
- if (!closestCluster) {
- closestCluster = (*iclust);
- bestRadialDistance = radialDistance;
- }
- else if (radialDistance<bestRadialDistance) {
- closestCluster = (*iclust);
- bestRadialDistance = radialDistance;
- }
- }
-
- if ((*iclust)->htime() - (*itrack)->htime() > m_tResidualCut) break;
- }
-
- double trackX = interceptUL.x();
- double trackY = interceptUL.y();
-
- if (closestCluster != NULL && closestCluster->charge() > m_chargeCutLow) {
- closestCluster->setAssociated(true);
- m_nTrackedClusters++;
- m_effHitmap->Fill(true, trackX, trackY);
- m_effX->Fill(true, trackX);
- m_effY->Fill(true, trackY);
-
- if (fmod(trackX, m_pitch)<0.0183 && fmod(trackY, m_pitch) < 0.0183) m_nClustersPassedCorner++;
- if (fmod(trackX, m_pitch) > 0.0183 &&
- fmod(trackX, m_pitch) < (2*0.0183) &&
- fmod(trackY, m_pitch) > 0.0183 &&
- fmod(trackY, m_pitch) < (2*0.0183)) m_nClustersPassedCentral++;
-
-
- //h_effVsCharge->Fill(h_effInter->Interpolate(fmod(trackX, m_pitch), fmod(trackY, m_pitch)), closestCluster->charge());
- // if (trackY > 9.25 && trackY < 10.25) m_effHitmapInterPixel->Fill(true, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
- // if (trackX > 28.05 &&
- // trackX < 28.6) m_effHitmapInterPixelTriple->Fill(true, trackX-28.05, fmod(trackY, m_pitch));
-
-
- m_effHitmapInterPixel->Fill(true, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
- m_effHitmapInterPixelTriple->Fill(true, trackX-28.05, fmod(trackY, m_pitch));
-
-
- if (closestCluster->size() <= m_effHitmapInterPixelVsSizes.size())
- m_effHitmapInterPixelVsSizes[closestCluster->size()-1]->Fill(true, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
- m_trackAssociated->at(itrack - cutTracks->begin()) = true;
- plot(closestCluster->x() - interceptUG.x(), "xResidualsClustMinusTrack", "xResidualsClustMinusTrack", -0.2, 0.2, 400);
- plot(closestCluster->y() - interceptUG.y(), "yResidualsClustMinusTrack", "yResidualsClustMinusTrack", -0.2, 0.2, 400);
- plot(closestCluster->htime() - (*itrack)->htime(), "hTimeResidualClustMinusTrack", "hTimeResidualClustMinusTrack", -50, 50, 400);
- if (closestCluster->size() == 1) {
- auto hits = closestCluster->hits();
- int col = (*hits.begin())->col();
- m_timeResidualVsColumn->fill(col, closestCluster->htime() - (*itrack)->htime());
- }
- if ((*itrack)->vertexed()) counter("nVerticesCorrelated") += 1;
- }
- else {
- m_effHitmap->Fill(false, trackX, trackY);
- // if (trackY > 9.25 && trackY < 10.25) m_effHitmapInterPixel->Fill(false, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
- // if (trackX > 28.05 &&
- // trackX < 28.6) m_effHitmapInterPixelTriple->Fill(false, trackX-28.05, fmod(trackY, m_pitch));
-
- m_effHitmapInterPixel->Fill(false, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
- m_effHitmapInterPixelTriple->Fill(false, trackX-28.05, fmod(trackY, m_pitch));
- m_effX->Fill(false, trackX);
- m_effY->Fill(false, trackY);
- //std::cout<<"Inefficiency at time: "<<(*itrack)->htime()<<std::endl;
- }
- if (fmod(trackX, m_pitch)<0.0183 && fmod(trackY, m_pitch) < 0.0183) m_nTracksCorner++;
- if (fmod(trackX, m_pitch) > 0.0183 &&
- fmod(trackX, m_pitch) < (2*0.0183) &&
- fmod(trackY, m_pitch) > 0.0183 &&
- fmod(trackY, m_pitch) < (2*0.0183)) m_nTracksCentral++;
- }
- }
-
-
-
- //=============================================================================
- //
- //=============================================================================
- double TbEffPur::getRadialSeparation(LHCb::TbCluster * cluster,
- LHCb::TbTrack * track)
- {
- Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
- double xResidual = cluster->x() - trackIntercept.x();
- double yResidual = cluster->y() - trackIntercept.y();
-
- double r2 = xResidual*xResidual + yResidual*yResidual;
- return pow(r2, 0.5);
- }
-
-
- //=============================================================================
- // Local cuts
- //=============================================================================
- bool TbEffPur::matchTrackToCluster(LHCb::TbCluster * cluster,
- LHCb::TbTrack * track)
- {
- Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
- double xResidual = cluster->x() - trackIntercept.x();
- double yResidual = cluster->y() - trackIntercept.y();
- double tResidual = cluster->htime() - track->htime();
- double delr = pow(pow(xResidual, 2) + pow(yResidual, 2), 0.5);
-
- // if (track->vertexed()) plot(xResidual, "vertexResiduals/X", "vertexResiduals/X", -0.5, 0.5, 200);
- // if (track->vertexed()) plot(yResidual, "vertexResiduals/Y", "vertexResiduals/Y", -0.5, 0.5, 200);
- // if (track->vertexed()) plot(tResidual, "vertexResiduals/T", "vertexResiduals/T", -0.5, 0.5, 200);
- //
- // if (track->vertexed() && litPixel(cluster, track)) {
- // plot(xResidual, "vertexResiduals/Xlit", "vertexResiduals/Xlit", -0.5, 0.5, 200);
- // plot(yResidual, "vertexResiduals/Ylit", "vertexResiduals/Ylit", -0.5, 0.5, 200);
- // }
-
-
- if (fabs(tResidual) > m_tResidualCut) {
- counter("nTimeRejected") += 1;
- return false;
- }
-
- if (delr > m_rResidualCut) {
- counter("nSpatialRejected") += 1;
- if (litPixel(cluster, track)) return true;
- else {
- counter("nSpatialAndLitRejected") += 1;
- return false;
- }
- }
-
- // if (!litPixel(cluster, track)) return false;
-
- return true;
- }
-
-
- //=============================================================================
- //
- //=============================================================================
- bool TbEffPur::interceptDeadPixel(Gaudi::XYZPoint trackInterceptGlobal)
- {
- // Find the row and column corresponding to the track intercept.
- Gaudi::XYZPoint trackInterceptLocal = geomSvc()->globalToLocal(trackInterceptGlobal, m_DUTindex);
- int row = trackInterceptLocal.y()/m_pitch - 0.5;
- int col = trackInterceptLocal.x()/m_pitch - 0.5;
-
- for (int icol = col - 1; icol != col+2; icol++) {
- for (int irow = row - 1; irow != row+2; irow++) {
- if (icol >= 0 && icol <256 && irow>=0 && irow<256) {
- if (pixelSvc()->isMasked(pixelSvc()->address(icol, irow), m_DUTindex)) {
- return true;
- }
-
- // if (icol == 17 && irow == 36) {
- // return true;
- // }
- // if (icol == 18 && irow == 36) {
- // return true;
- // }
- // if (icol == 53 && irow == 3) {
- // return true;
- // }
- // if (icol == 54 && irow == 3) {
- // return true;
- // }
- // if (icol == 73 && irow == 26) {
- // return true;
- // }
- // if (icol == 74 && irow == 26) {
- // return true;
- // }
- // if (icol == 19 && irow == 103) {
- // return true;
- // }
- // if (icol == 31 && irow == 106) {
- // return true;
- // }
- // if (icol == 32 && irow == 106) {
- // return true;
- // }
- // if (icol == 38 && irow == 108) {
- // return true;
- // }
- // if (icol == 63 && irow == 95) {
- // return true;
- // }
- // if (icol == 64 && irow == 95) {
- // return true;
- // }
- // if (icol == 103 && irow == 23) {
- // return true;
- // }
- // if (icol == 104 && irow == 23) {
- // return true;
- // }
- // if (icol == 115 && irow == 20) {
- // return true;
- // }
- // if (icol == 116 && irow == 94) {
- // return true;
- // }
- }
- }
- }
-
-
-
- return false;
- }
-
- //=============================================================================
- //
- //=============================================================================
- bool TbEffPur::litPixel(LHCb::TbCluster * cluster, LHCb::TbTrack * track)
- {
- // Find the row and column corresponding to the track intercept.
- Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
- Gaudi::XYZPoint pLocal = geomSvc()->globalToLocal(trackIntercept, m_DUTindex);
- int row = pLocal.y()/m_pitch - 0.5;
- int col = pLocal.x()/m_pitch - 0.5;
-
- // See if within a pre-defined square of pixels.
- for (unsigned int i=0; i<cluster->size(); i++) {
- unsigned int delRow = abs(row - int(cluster->hits()[i]->row()));
- unsigned int delCol = abs(col - int(cluster->hits()[i]->col()));
-
- if (delRow <= m_litSquareSide && delCol <= m_litSquareSide) {
- counter("nLitPixels") += 1;
- return true;
- }
- }
-
- return false;
- }
-
-
- //=============================================================================
- // Excluding dead regions.
- //=============================================================================
- bool TbEffPur::globalCutPosition(Gaudi::XYZPoint pGlobal)
- {
- if (pGlobal.x() < m_xLow + m_edgeVetoDistance) return false;
- else if (pGlobal.x() > m_xUp - m_edgeVetoDistance) return false;
- else if (pGlobal.y() < m_yLow + m_edgeVetoDistance) return false;
- else if (pGlobal.y() > m_yUp - m_edgeVetoDistance) return false;
- return true; // Inside.
- }
-
-
- //=============================================================================
- // Excluding dead regions.
- //=============================================================================
- bool TbEffPur::outsideDUT(Gaudi::XYZPoint interceptUG)
- {
- const auto interceptUL = geomSvc()->globalToLocal(interceptUG, m_DUTindex);
- if (interceptUL.x() < 0 + m_edgeVetoDistance ||
- interceptUL.x() > 14.08 - m_edgeVetoDistance ||
- interceptUL.y() < 0 + m_edgeVetoDistance||
- interceptUL.y() > 14.08 - m_edgeVetoDistance) return true;
- return false;
- }
-
-
-
- //=============================================================================
- // Viewer output.
- //=============================================================================
- //=============================================================================
- // Viewer output.
- //=============================================================================
- //=============================================================================
- // Viewer output.
- //=============================================================================
-
- void TbEffPur::outputDeadRegion(unsigned int col, unsigned int row) {
- std::ofstream myfile;
- myfile.open("/afs/cern.ch/user/d/dsaunder/KeplerViewerData.dat", std::ofstream::app);
- std::string outputLine;
- for (int irow = std::max(int(row-m_deadAreaRadius), 0);
- irow < std::min(int(row+m_deadAreaRadius+1), m_nDUTpixels); irow++) {
- for (int icol = std::max(int(col-m_deadAreaRadius), 0);
- icol < std::min(int(col+m_deadAreaRadius+1), m_nDUTpixels); icol++) {
-
-
- outputLine = "DeadPixel ";
- const double xLocal = (col-icol) * m_pitch; // Dont want the middle!
- const double yLocal = (row-irow) * m_pitch; // Dont want the middle!
- Gaudi::XYZPoint pLocal(xLocal, yLocal, 0.);
-
- Gaudi::XYZPoint posn = geomSvc()->localToGlobal(pLocal, m_DUTindex);
- outputLine += std::to_string(posn.x()) + " ";
- outputLine += std::to_string(posn.y()) + " ";
- outputLine += std::to_string(posn.z()) + " ";
-
- Gaudi::XYZPoint posn2(pLocal.x() + 0.055, pLocal.y(), 0.);
- posn = geomSvc()->localToGlobal(posn2, m_DUTindex);
- outputLine += std::to_string(posn.x()) + " ";
- outputLine += std::to_string(posn.y()) + " ";
- outputLine += std::to_string(posn.z()) + " ";
-
- Gaudi::XYZPoint posn3(pLocal.x() + 0.055, pLocal.y()+0.055, 0.);
- posn = geomSvc()->localToGlobal(posn3, m_DUTindex);
- outputLine += std::to_string(posn.x()) + " ";
- outputLine += std::to_string(posn.y()) + " ";
- outputLine += std::to_string(posn.z()) + " ";
-
- Gaudi::XYZPoint posn4(pLocal.x(), pLocal.y()+0.055, 0.);
- posn = geomSvc()->localToGlobal(posn4, m_DUTindex);
- outputLine += std::to_string(posn.x()) + " ";
- outputLine += std::to_string(posn.y()) + " ";
- outputLine += std::to_string(posn.z()) + " ";
-
- outputLine += "\n";
- myfile << outputLine;
- }
- }
- myfile.close();
- }
-
-
- //=============================================================================
- // Viewer output.
- //=============================================================================
-
- void TbEffPur::outputViewerData()
- {
- std::ofstream myfile;
- myfile.open("/afs/cern.ch/user/d/dsaunder/KeplerViewerData.dat", std::ofstream::app);
- std::string outputLine;
- //
- // // First output the chips.
- // for (unsigned int i=0; i<m_nPlanes; i++) {
- // outputLine = "Chip ";
- // Gaudi::XYZPoint posn1(0., 14.08, 0.);
- // Gaudi::XYZPoint posn = geomSvc()->localToGlobal(posn1, i);
- // outputLine += std::to_string(posn.x()) + " ";
- // outputLine += std::to_string(posn.y()) + " ";
- // outputLine += std::to_string(posn.z()) + " ";
- //
- // Gaudi::XYZPoint posn2(14.08, 14.08, 0.);
- // posn = geomSvc()->localToGlobal(posn2, i);
- // outputLine += std::to_string(posn.x()) + " ";
- // outputLine += std::to_string(posn.y()) + " ";
- // outputLine += std::to_string(posn.z()) + " ";
- //
- // Gaudi::XYZPoint posn3(14.08, 0., 0.);
- // posn = geomSvc()->localToGlobal(posn3, i);
- // outputLine += std::to_string(posn.x()) + " ";
- // outputLine += std::to_string(posn.y()) + " ";
- // outputLine += std::to_string(posn.z()) + " ";
- //
- // Gaudi::XYZPoint posn4(0., 0., 0.);
- // posn = geomSvc()->localToGlobal(posn4, i);
- // outputLine += std::to_string(posn.x()) + " ";
- // outputLine += std::to_string(posn.y()) + " ";
- // outputLine += std::to_string(posn.z()) + " ";
- //
- // outputLine += "\n";
- // myfile << outputLine;
- // }
-
-
- // Clusters.
- auto ic = m_clusters->begin();
- const auto end = m_clusters->end();
- for (; ic != end; ++ic) {
- outputLine = "Cluster ";
- outputLine += std::to_string((*ic)->x()) + " ";
- outputLine += std::to_string((*ic)->y()) + " ";
- outputLine += std::to_string((*ic)->z()) + " ";
- outputLine += std::to_string((*ic)->htime()) + " ";
-
- // if ((*ic)->endCluster() && (*ic)->vertexed()) outputLine += "5 \n";
- //if ((*ic)->vertexed()) outputLine += "4 \n";
- // else if ((*ic)->endCluster()) outputLine += "3 \n";
- if ((*ic)->associated()) outputLine += "2 \n";
- // else if ((*ic)->volumed()) outputLine += "1 \n";
-
- else {
- outputLine += "0 \n";
- }
- myfile << outputLine;
- }
-
- outputLine = "CentralRegion ";
- outputLine += std::to_string(m_xLow) + " ";
- outputLine += std::to_string(m_xUp) + " ";
- outputLine += std::to_string(m_yLow) + " ";
- outputLine += std::to_string(m_yUp) + " ";
- outputLine += std::to_string(geomSvc()->module(m_DUTindex)->z()) + " ";
- myfile << outputLine;
-
- // // Its hits.
- // for (auto hit : (*ic)->hits()) {
- // outputLine = "Pixel ";
- // const double xLocal = (hit->col()) * m_pitch; // Dont want the middle!
- // const double yLocal = (hit->row()) * m_pitch; // Dont want the middle!
- // Gaudi::XYZPoint pLocal(xLocal, yLocal, 0.);
- //
- // Gaudi::XYZPoint posn = geomSvc()->localToGlobal(pLocal, i);
- // outputLine += std::to_string(posn.x()) + " ";
- // outputLine += std::to_string(posn.y()) + " ";
- // outputLine += std::to_string(posn.z()) + " ";
- //
- // Gaudi::XYZPoint posn2(pLocal.x() + 0.055, pLocal.y(), 0.);
- // posn = geomSvc()->localToGlobal(posn2, i);
- // outputLine += std::to_string(posn.x()) + " ";
- // outputLine += std::to_string(posn.y()) + " ";
- // outputLine += std::to_string(posn.z()) + " ";
- //
- // Gaudi::XYZPoint posn3(pLocal.x() + 0.055, pLocal.y()+0.055, 0.);
- // posn = geomSvc()->localToGlobal(posn3, i);
- // outputLine += std::to_string(posn.x()) + " ";
- // outputLine += std::to_string(posn.y()) + " ";
- // outputLine += std::to_string(posn.z()) + " ";
- //
- // Gaudi::XYZPoint posn4(pLocal.x(), pLocal.y()+0.055, 0.);
- // posn = geomSvc()->localToGlobal(posn4, i);
- // outputLine += std::to_string(posn.x()) + " ";
- // outputLine += std::to_string(posn.y()) + " ";
- // outputLine += std::to_string(posn.z()) + " ";
- //
- // outputLine += std::to_string(hit->htime()) + " ";
- // outputLine += std::to_string(hit->ToT()) + " ";
- //
- // outputLine += "\n";
- // myfile << outputLine;
- // }
-
- myfile.close();
- }
-
-
- //=============================================================================
- // END
- //=============================================================================