Newer
Older
TB_Chris / TbAnalysis / src / .svn / text-base / TbEffPur.cpp.svn-base
  1. #include <fstream>
  2.  
  3. // ROOT
  4. #include "TMath.h"
  5.  
  6. // Gaudi
  7. #include "GaudiKernel/PhysicalConstants.h"
  8. #include "GaudiUtils/HistoLabels.h"
  9.  
  10. // Tb/TbKernel
  11. #include "TbKernel/TbConstants.h"
  12. #include "TbKernel/TbModule.h"
  13.  
  14. // Local
  15. #include "TbEffPur.h"
  16.  
  17. using namespace Gaudi::Utils::Histos;
  18.  
  19. DECLARE_ALGORITHM_FACTORY(TbEffPur)
  20.  
  21. //=============================================================================
  22. // Standard constructor
  23. //=============================================================================
  24. TbEffPur::TbEffPur(const std::string& name, ISvcLocator* pSvcLocator)
  25. : TbAlgorithm(name, pSvcLocator),
  26. m_pitch(0.055),
  27. m_nDUTpixels(256),
  28. m_nTracks(0),
  29. m_nClusters(0),
  30. m_nTrackedClusters(0),
  31. m_nClustersPassedCentral(0),
  32. m_nTracksCentral(0),
  33. m_nClustersPassedCorner(0),
  34. m_nTracksCorner(0),
  35. m_eff(0.),
  36. m_pur(0.),
  37. m_deadAreaRadius(2.),
  38. m_trackAssociated(NULL)
  39. {
  40. declareProperty("TrackLocation",
  41. m_trackLocation = LHCb::TbTrackLocation::Default);
  42. declareProperty("VertexLocation",
  43. m_vertexLocation = LHCb::TbVertexLocation::Default);
  44. declareProperty("ClusterLocation",
  45. m_clusterLocation = LHCb::TbClusterLocation::Default);
  46.  
  47. // Parameters.
  48. declareProperty("DUTindex", m_DUTindex = 4);
  49. declareProperty("GlobalCutXLow", m_xLow = 3);
  50. declareProperty("GlobalCutXUp", m_xUp = 9);
  51. declareProperty("GlobalCutYLow", m_yLow = 3);
  52. declareProperty("GlobalCutYUp", m_yUp = 9);
  53. declareProperty("LocalProbabilityCut", m_probCut = 0.5);
  54.  
  55. declareProperty("RResidualCut", m_rResidualCut = 0.1);
  56. declareProperty("TResidualCut", m_tResidualCut = 100);
  57. declareProperty("ChargeCutLow", m_chargeCutLow = 0);
  58. // ... note this global cut enforced at different point to above.
  59. declareProperty("ChargeCutUp", m_chargeCutUp = 1000000);
  60. declareProperty("LitSquareSide", m_litSquareSide = 0);
  61.  
  62. declareProperty("ViewerOutput", m_viewerOutput = false);
  63. declareProperty("ViewerEventNum", m_viewerEvent = 100);
  64. declareProperty("TGap", m_tGap = 200);
  65. declareProperty("CorrelationTimeWindow", m_correlationTimeWindow = 100);
  66. declareProperty("ApplyVeto", m_applyVeto = true);
  67. declareProperty("TelescopeClusterVetoDelT", m_telescopeClusterVetoDelT = 30);
  68. declareProperty("EdgeVetoDistance", m_edgeVetoDistance = 0.025);
  69. m_nEvent = -1;
  70. }
  71.  
  72. //=============================================================================
  73. // Initialization
  74. //=============================================================================
  75. StatusCode TbEffPur::initialize() {
  76.  
  77. // Initialise the base class.
  78. StatusCode sc = TbAlgorithm::initialize();
  79. if (sc.isFailure()) return sc;
  80.  
  81. // Efficiency objects.
  82. m_effX = new TEfficiency("efficiencyX", "efficiencyX", 40, 0.0, 0.11);
  83. m_effY = new TEfficiency("efficiencyY", "efficiencyY", 40, 0.0, 0.11);
  84. m_effs = new TEfficiency("efficiency", "efficiency", 1, -0.5, 0.5);
  85. m_purs = new TEfficiency("purity", "Tb/TbEffPur/pur", 1, -0.5, 0.5);
  86. m_effHitmap = new TEfficiency("efficiencyHitmap", "efficiencyHitmap",
  87. 64, 0.0, 14.08, 64, 0.0, 14.08);
  88. m_purHitmap = new TEfficiency("puritHitmap", "purityHitmap",
  89. 64, 0.0, 14.08, 64, 0.0, 14.08);
  90.  
  91. int nBins = 50;
  92. m_effHitmapInterPixel = new TEfficiency("efficiencyHitmapInterPixel", "efficiencyHitmapInterPixel",
  93. nBins, 0.0, 0.055, nBins, 0.0, 0.055);
  94. m_purHitmapInterPixel = new TEfficiency("puritHitmapInterPixel", "purityHitmapInterPixel",
  95. nBins, 0.0, 0.055, nBins, 0.0, 0.055);
  96.  
  97. m_effHitmapInterPixelTriple = new TEfficiency("efficiencyHitmapInterPixelTriple", "efficiencyHitmapInterPixelTriple",
  98. 150, 0.0, 10*0.055, 20, 0.0, 2*0.055);
  99.  
  100. for (unsigned int i=1; i<5; i++) {
  101. std::string name = "efficiencyHitmapInterClusterSize" + std::to_string(i);
  102. TEfficiency * e = new TEfficiency(name.c_str(), name.c_str(),
  103. nBins, 0.0, 0.055, nBins, 0.0, 0.055);
  104. m_effHitmapInterPixelVsSizes.push_back(e);
  105. }
  106.  
  107. // Plots for post correlations.
  108. m_remainsCorrelationsX = book2D("remainsClusterCorrelationsX", "remainsClusterCorrelationsX", 0, 14, 200, 0, 14, 200);
  109. m_remainsCorrelationsY = book2D("remainsClusterCorrelationsY", "remainsClusterCorrelationsY", 0, 14, 200, 0, 14, 200);
  110. m_remainsDifferencesXY = book2D("remainsDifferencesXY", "remainsDifferencesXY", -1.5, 1.5, 150, -1.5, 1.5, 150);
  111. m_clusterRemainsPositionsLocal = book2D("clusterRemainsPositionsLocal", "clusterRemainsPositionsLocal", -5, 20, 600, -5, 20, 600);
  112. m_trackRemainsPositionsLocal = book2D("trackRemainsPositionsLocal", "trackRemainsPositionsLocal", 0, 256, 600, 0, 256, 600);
  113. m_clusterRemainsPositionsGlobal = book2D("clusterRemainsPositionsGlobal", "clusterRemainsPositionsGlobal", -5, 20, 600, -5, 20, 600);
  114. m_trackRemainsPositionsGlobal = book2D("trackRemainsPositionsGlobal", "trackRemainsPositionsGlobal", -5, 20, 600, -5, 20, 600);
  115. m_vetoTracksHitmap = book2D("vetoTrackHitmap", "vetoTrackHitmap", -5, 20, 600, -5, 20, 600);
  116. m_vetoClustersHitmap = book2D("vetoClusterHitmap", "vetoClusterHitmap", -5, 20, 600, -5, 20, 600);
  117. m_timeResidualVsColumn = book2D("timeResidualVsColumn", "timeResidualVsColumn", 0, 256, 256, -50, 50, 500);
  118.  
  119.  
  120. // Setup the tools.
  121. m_trackFit = tool<ITbTrackFit>("TbTrackFit", "Fitter", this);
  122. // h_effVsCharge = new TH2F("effVsCharge", "effVsCharge", 50, 0.5, 1, 20, 0.5, 20.5);
  123. // TFile * f = new TFile("/afs/cern.ch/user/d/dsaunder/cmtuser/KEPLER/KEPLER_HEAD/EDR_plots/EfficiencyResults/Eff-S22-Run6309-500V.root", "READ");
  124. // h_effInter = (TH2F*)f->Get("eff_histo");
  125. return StatusCode::SUCCESS;
  126.  
  127. }
  128.  
  129.  
  130. //=============================================================================
  131. // Finalizer
  132. //=============================================================================
  133. StatusCode TbEffPur::finalize() {
  134. m_effs->SetTotalEvents(0, m_nTracks);
  135. m_effs->SetPassedEvents(0, m_nTrackedClusters);
  136.  
  137. m_purs->SetTotalEvents(0, m_nClusters);
  138. m_purs->SetPassedEvents(0, m_nTrackedClusters);
  139.  
  140. m_eff = m_nTrackedClusters/(1.0*m_nTracks);
  141. m_pur = m_nTrackedClusters/(1.0*m_nClusters);
  142. std::cout<<"ith ROOT Efficiency: " << m_effs->GetEfficiency(0) << " + " << m_effs->GetEfficiencyErrorUp(0) <<" - "<< m_effs->GetEfficiencyErrorLow(0) <<std::endl;
  143. std::cout<<"ith ROOT Purity: " << m_purs->GetEfficiency(0) << " + " << m_purs->GetEfficiencyErrorUp(0) <<" - "<< m_purs->GetEfficiencyErrorLow(0) <<std::endl;
  144.  
  145. std::cout<<"Corner efficiency: "<<m_nClustersPassedCorner/(1.*m_nTracksCorner)<<std::endl;
  146. std::cout<<"Central efficiency: "<<m_nClustersPassedCentral/(1.*m_nTracksCentral)<<std::endl;
  147.  
  148. plot(m_eff, "Efficiency", "Efficiency", 0.0, 1, 1);
  149. plot(m_pur, "Purity", "Purity", 0.0, 1, 1);
  150. plot(m_effs->GetEfficiencyErrorUp(0), "EfficiencyError", "EfficiencyError", 0.0, 1, 1);
  151. plot(m_purs->GetEfficiencyErrorUp(0), "PurityError", "PurityError", 0.0, 1, 1);
  152.  
  153. TH2D * h = Gaudi::Utils::Aida2ROOT::aida2root(m_trackRemainsPositionsLocal);
  154. double maxBin = h->GetBinContent(h->GetMaximumBin());
  155. double minBin = h->GetBinContent(h->GetMinimumBin());
  156. for (int i=0; i<m_trackRemainsPositionsLocal->xAxis().bins(); i++) {
  157. for (int j=0; j<m_trackRemainsPositionsLocal->yAxis().bins(); j++) {
  158. plot(m_trackRemainsPositionsLocal->binHeight(i, j), "trackRemainsBins",
  159. "trackRemainsBins", minBin, maxBin, int(maxBin-minBin));
  160. }
  161. }
  162.  
  163. TFile * f = new TFile("EfficiencyResults.root", "RECREATE");
  164. //m_effs->Write();
  165.  
  166. m_effHitmapInterPixel->Write();
  167. m_effHitmapInterPixelTriple->Write();
  168. //h_effVsCharge->Write();
  169. //m_effs->CreateHistogram()->Write();
  170. m_effHitmapInterPixel->CreateHistogram()->Write();
  171. m_effHitmapInterPixelTriple->CreateHistogram()->Write();
  172.  
  173.  
  174.  
  175. // for (unsigned int i=0; i<m_effHitmapInterPixelVsSizes.size(); i++) {
  176. // m_effHitmapInterPixelVsSizes[i]->Write();
  177. // //m_effHitmapInterPixelVsSizes[i]->CreateHistogram()->Write();
  178. // }
  179. f->Close();
  180. return TbAlgorithm::finalize();
  181. }
  182.  
  183.  
  184. //=============================================================================
  185. // Main execution
  186. //=============================================================================
  187. StatusCode TbEffPur::execute()
  188. {
  189. m_nEvent++;
  190.  
  191. m_tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  192. //m_vertices = getIfExists<LHCb::TbVertices>(m_vertexLocation);
  193. m_clusters = getIfExists<LHCb::TbClusters>(m_clusterLocation + std::to_string(m_DUTindex));
  194.  
  195. if (m_clusters->size() < 100 || m_tracks->size() < 100) return StatusCode::SUCCESS;
  196.  
  197. // Matching.
  198. effPur();
  199. return StatusCode::SUCCESS;
  200. }
  201.  
  202.  
  203. //=============================================================================
  204. // Efficiency finding
  205. //=============================================================================
  206. void TbEffPur::effPur()
  207. {
  208. // Non-veto'd containers.
  209. std::vector<LHCb::TbCluster*> * cutClusters = new std::vector<LHCb::TbCluster*>();
  210. std::vector<LHCb::TbTrack*> cutTracks;
  211. int nClustersPreVeto = m_clusters->size();
  212. int nTracksPreVeto = m_tracks->size();
  213.  
  214.  
  215. // Apply veto.
  216. if (m_applyVeto) applyVeto(cutClusters, &cutTracks);
  217. else fillAllTrackClusters(cutClusters, &cutTracks);
  218. m_trackAssociated = new std::vector<bool>(cutTracks.size(), false);
  219. //std::cout<<"Fraction tracks not veto'd: "<<cutTracks.size()/(1.0*nTracksPreVeto)<<"\tEvent: "<<m_nEvent<<std::endl;
  220. //std::cout<<"Fraction clusters not veto'd: "<<cutClusters->size()/(1.0*nClustersPreVeto)<<"\tEvent: "<<m_nEvent<<std::endl;
  221. plot(cutClusters->size()/(1.0*nClustersPreVeto), "clusterFractionVetod", "clusterFractionVetod", 0.0, 1.0, 200);
  222. plot(cutTracks.size()/(1.0*nTracksPreVeto), "trackFractionVetod", "trackFractionVetod", 0.0, 1.0, 200);
  223.  
  224.  
  225. // Do matching.
  226. m_nClusters += cutClusters->size();
  227. m_nTracks += cutTracks.size();
  228. trackClusters(cutClusters, &cutTracks);
  229.  
  230.  
  231. // Viewer options.
  232. if (m_nEvent == m_viewerEvent) {
  233. if (m_viewerOutput) outputViewerData();
  234. for (LHCb::TbClusters::iterator iclust = cutClusters->begin();
  235. iclust != cutClusters->end(); iclust++) {
  236. // if (!(*iclust)->associated())
  237. // std::cout<<"Note: non-associated clusters at hTime: "<<(*iclust)->htime()<<std::endl;
  238. }
  239. }
  240.  
  241. // Correlate remains.
  242. //correlateRemains(&cutTracks, cutClusters);
  243. delete cutClusters;
  244. delete m_trackAssociated;
  245. }
  246.  
  247.  
  248. //=============================================================================
  249. // Global cuts
  250. //=============================================================================
  251. void TbEffPur::applyVeto(std::vector<LHCb::TbCluster*> * cutClusters,
  252. std::vector<LHCb::TbTrack*> * cutTracks)
  253. {
  254. // Gather all cluster times (over all planes) and sort it ___________________
  255. std::vector<double> allClusterTimes;
  256. for (unsigned int i=0; i<m_nPlanes; i++) {
  257. LHCb::TbClusters * clusters = getIfExists<LHCb::TbClusters>(m_clusterLocation + std::to_string(i));
  258. for (LHCb::TbClusters::const_iterator it = clusters->begin(); it != clusters->end(); ++it)
  259. allClusterTimes.push_back((*it)->htime());
  260. }
  261. std::sort(allClusterTimes.begin(), allClusterTimes.end());
  262.  
  263.  
  264.  
  265. // Get big enough gaps in pairs _____________________________________________
  266. std::vector<double> tGapCenters;
  267. tGapCenters.push_back(0);
  268. for (std::vector<double>::iterator itime = allClusterTimes.begin()+1;
  269. itime != allClusterTimes.end(); itime++) {
  270. plot((*itime) - (*(itime-1)), "allGapWidths", "allGapWidths", 0.0, 3000.0, 200);
  271. if ((*itime) - (*(itime-1)) > m_tGap) {
  272. tGapCenters.push_back(0.5*((*itime) + (*(itime-1))));
  273. plot((*itime) - (*(itime-1)), "bigGapWidths", "bigGapWidths", 0.0, 3000.0, 200);
  274. }
  275. }
  276. tGapCenters.push_back(allClusterTimes.back());
  277. counter("nGaps") += tGapCenters.size();
  278. plot(tGapCenters.size(), "nGaps", "nGaps", 0.0, 1000.0, 200);
  279.  
  280.  
  281.  
  282. // Loop over these sub events and veto ______________________________________
  283. std::vector<double>::iterator igap;
  284. for (igap = tGapCenters.begin(); igap != tGapCenters.end() - 1; igap++) {
  285. bool veto = false;
  286.  
  287. if (igap == tGapCenters.begin() || igap == tGapCenters.end() -1) veto = true;
  288.  
  289. // Loop tracks.
  290. if (!veto) {
  291. for (LHCb::TbTrack* track : *m_tracks) {
  292. if (track->htime() > (*igap) && track->htime() < (*(igap+1))) {
  293. // Veto on tracks outside cuts.
  294. const Gaudi::XYZPoint trackInterceptGlobal = geomSvc()->intercept(track, m_DUTindex);
  295. if (!globalCutPosition(trackInterceptGlobal)) veto = true;
  296. else if (outsideDUT(trackInterceptGlobal)) veto = true;
  297. else if (interceptDeadPixel(trackInterceptGlobal)) veto = true;
  298. }
  299. }
  300. }
  301.  
  302.  
  303. // Loop vertices.
  304. // if (!veto) {
  305. // for (LHCb::TbVertex* vertex : *m_vertices) {
  306. // if (vertex->htime() > (*igap) && vertex->htime() < (*(igap+1))) {
  307. // veto = true;
  308. // }
  309. // }
  310. // }
  311.  
  312.  
  313. // Loop DUT clusters.
  314. if (!veto) {
  315. for (const LHCb::TbCluster* cluster : *m_clusters) {
  316. if (cluster->htime() > (*igap) && cluster->htime() < (*(igap+1))) {
  317. // Veto on clusters outside cuts.
  318. Gaudi::XYZPoint clusterInterceptGlobal(cluster->x(), cluster->y(), cluster->z());
  319. if (!globalCutPosition(clusterInterceptGlobal)) veto = true;
  320. if (cluster->charge() > m_chargeCutUp || cluster->charge() < m_chargeCutLow) veto = true;
  321. }
  322. }
  323. }
  324.  
  325.  
  326. // Loop telescope clusters.
  327. if (!veto) {
  328. for (unsigned int i=0; i<m_nPlanes; i++) {
  329. if (i == m_DUTindex) continue;
  330. LHCb::TbClusters * clusters = getIfExists<LHCb::TbClusters>(m_clusterLocation + std::to_string(i));
  331. for (LHCb::TbClusters::const_iterator it = clusters->begin(); it != clusters->end()-1; ++it) {
  332. if ((*it)->htime() > (*igap) && (*it)->htime() < (*(igap+1))) {
  333. double delt = (*it)->htime() - (*(it+1))->htime();
  334. if (fabs(delt) < m_telescopeClusterVetoDelT) {
  335. veto = true;
  336. break;
  337. }
  338. }
  339. }
  340. }
  341. }
  342.  
  343.  
  344. if (!veto) fillTrackClusters(cutClusters, cutTracks, (*igap), (*(igap+1)));
  345. else counter("nGapsFiltered") += 1;
  346. }
  347. }
  348.  
  349.  
  350.  
  351. //=============================================================================
  352. // Correlate remains
  353. //=============================================================================
  354. void TbEffPur::correlateRemains(std::vector<LHCb::TbTrack*> * cutTracks,
  355. std::vector<LHCb::TbCluster*> * cutClusters)
  356. {
  357. for (std::vector<LHCb::TbCluster*>::iterator ic = cutClusters->begin();
  358. ic != cutClusters->end(); ++ic) {
  359. if (!(*ic)->associated()) {
  360. m_clusterRemainsPositionsGlobal->fill((*ic)->x(), (*ic)->y());
  361. m_clusterRemainsPositionsLocal->fill((*ic)->xloc(), (*ic)->yloc());
  362. plot((*ic)->charge(), "chargeClusterRemains", "chargeClusterRemains", 0.0, 1000.0, 125);
  363. }
  364. }
  365.  
  366. unsigned int i = 0;
  367. for (std::vector<LHCb::TbTrack*>::iterator itrack = cutTracks->begin();
  368. itrack != cutTracks->end(); itrack++) {
  369. if (!m_trackAssociated->at(i)) {
  370. Gaudi::XYZPoint trackIntercept = geomSvc()->intercept((*itrack), m_DUTindex);
  371. m_trackRemainsPositionsGlobal->fill(trackIntercept.x(), trackIntercept.y());
  372. auto interceptUL = geomSvc()->globalToLocal(trackIntercept, m_DUTindex);
  373. m_trackRemainsPositionsLocal->fill(interceptUL.x()/m_pitch - 0.5, interceptUL.y()/m_pitch - 0.5);
  374. }
  375. i++;
  376. }
  377.  
  378.  
  379. // Draws correlations plots between non-associated clusters and tracks.
  380. for (std::vector<LHCb::TbTrack*>::iterator itrack = cutTracks->begin();
  381. itrack != cutTracks->end(); itrack++) {
  382. for (std::vector<LHCb::TbCluster*>::iterator ic = cutClusters->begin();
  383. ic != cutClusters->end(); ++ic) {
  384. if ((*ic)->htime() < (*itrack)->htime() - m_correlationTimeWindow) continue;
  385. if (!(*ic)->associated()) {
  386. //if (m_nEvent == m_viewerEvent) std::cout<<"Impurity at time: "<<(*ic)->htime()<<std::endl;
  387. Gaudi::XYZPoint trackIntercept = geomSvc()->intercept((*itrack), m_DUTindex);
  388. m_remainsCorrelationsX->fill(trackIntercept.x(), (*ic)->x());
  389. m_remainsCorrelationsY->fill(trackIntercept.y(), (*ic)->y());
  390.  
  391. plot(trackIntercept.x() - (*ic)->x(), "remainClustDifferencesX", "remainsClustDifferencesX", -1.5, 1.5, 150);
  392. plot(trackIntercept.y() - (*ic)->y(), "remainClustDifferencesY", "remainClustDifferencesY", -1.5, 1.5, 150);
  393. plot((*itrack)->htime() - (*ic)->htime(), "remainClustDifferencesT", "remainClustDifferencesT", -50, 50, 150);
  394.  
  395. m_remainsDifferencesXY->fill(trackIntercept.y() - (*ic)->y(), trackIntercept.x() - (*ic)->x());
  396. }
  397. if ((*ic)->htime() > (*itrack)->htime() + m_correlationTimeWindow) break;
  398. }
  399. }
  400. }
  401.  
  402.  
  403. //=============================================================================
  404. //
  405. //=============================================================================
  406. void TbEffPur::fillAllTrackClusters(std::vector<LHCb::TbCluster*> * cutClusters,
  407. std::vector<LHCb::TbTrack*> * cutTracks)
  408. {
  409. for (LHCb::TbTrack* track : *m_tracks) cutTracks->push_back(track);
  410. for (std::vector<LHCb::TbCluster*>::iterator iClust = m_clusters->begin();
  411. iClust != m_clusters->end(); iClust++) cutClusters->push_back(*iClust);
  412. }
  413.  
  414.  
  415. //=============================================================================
  416. //
  417. //=============================================================================
  418. void TbEffPur::fillTrackClusters(std::vector<LHCb::TbCluster*> * cutClusters,
  419. std::vector<LHCb::TbTrack*> * cutTracks, double tlow, double tup)
  420. {
  421. // Push tracks and clusters inside this time window, which passed the veto.
  422. for (LHCb::TbTrack* track : *m_tracks) {
  423. if (track->htime() > tlow && track->htime() < tup) {
  424. Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
  425. cutTracks->push_back(track);
  426. m_vetoTracksHitmap->fill(trackIntercept.x(), trackIntercept.y());
  427.  
  428. Gaudi::XYZPoint trackInterceptLocal = geomSvc()->globalToLocal(trackIntercept, m_DUTindex);
  429. int row = trackInterceptLocal.y()/m_pitch - 0.5;
  430. int col = trackInterceptLocal.x()/m_pitch - 0.5;
  431.  
  432. if (pixelSvc()->isMasked(pixelSvc()->address(col, row), m_DUTindex)) {
  433. std::cout<<"Shouldn't be here!"<<std::endl;
  434. }
  435. }
  436. }
  437.  
  438. for (LHCb::TbCluster* cluster : *m_clusters) {
  439. if (cluster->htime() > tlow && cluster->htime() < tup) {
  440. cutClusters->push_back(cluster);
  441. m_vetoClustersHitmap->fill(cluster->x(), cluster->y());
  442. }
  443. }
  444. }
  445.  
  446.  
  447. //=============================================================================
  448. // Global cuts
  449. //=============================================================================
  450. void TbEffPur::trackClusters(std::vector<LHCb::TbCluster*> * cutClusters,
  451. std::vector<LHCb::TbTrack*> * cutTracks)
  452. {
  453. for (std::vector<LHCb::TbTrack*>::iterator itrack = cutTracks->begin();
  454. itrack != cutTracks->end(); itrack++) {
  455.  
  456. // Get local position of track.
  457. const auto interceptUG = geomSvc()->intercept((*itrack), m_DUTindex);
  458. const auto interceptUL = geomSvc()->globalToLocal(interceptUG, m_DUTindex);
  459.  
  460.  
  461. // Loop over clusters to find the closest.
  462. LHCb::TbCluster * closestCluster = NULL;
  463. double bestRadialDistance;
  464. for (std::vector<LHCb::TbCluster*>::iterator iclust = cutClusters->begin();
  465. iclust != cutClusters->end(); iclust++) {
  466. bool match = matchTrackToCluster((*iclust), (*itrack));
  467. double radialDistance = getRadialSeparation((*iclust), (*itrack));
  468. if (match) {
  469. if (!closestCluster) {
  470. closestCluster = (*iclust);
  471. bestRadialDistance = radialDistance;
  472. }
  473. else if (radialDistance<bestRadialDistance) {
  474. closestCluster = (*iclust);
  475. bestRadialDistance = radialDistance;
  476. }
  477. }
  478.  
  479. if ((*iclust)->htime() - (*itrack)->htime() > m_tResidualCut) break;
  480. }
  481.  
  482. double trackX = interceptUL.x();
  483. double trackY = interceptUL.y();
  484.  
  485. if (closestCluster != NULL && closestCluster->charge() > m_chargeCutLow) {
  486. closestCluster->setAssociated(true);
  487. m_nTrackedClusters++;
  488. m_effHitmap->Fill(true, trackX, trackY);
  489. m_effX->Fill(true, trackX);
  490. m_effY->Fill(true, trackY);
  491.  
  492. if (fmod(trackX, m_pitch)<0.0183 && fmod(trackY, m_pitch) < 0.0183) m_nClustersPassedCorner++;
  493. if (fmod(trackX, m_pitch) > 0.0183 &&
  494. fmod(trackX, m_pitch) < (2*0.0183) &&
  495. fmod(trackY, m_pitch) > 0.0183 &&
  496. fmod(trackY, m_pitch) < (2*0.0183)) m_nClustersPassedCentral++;
  497.  
  498.  
  499. //h_effVsCharge->Fill(h_effInter->Interpolate(fmod(trackX, m_pitch), fmod(trackY, m_pitch)), closestCluster->charge());
  500. // if (trackY > 9.25 && trackY < 10.25) m_effHitmapInterPixel->Fill(true, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
  501. // if (trackX > 28.05 &&
  502. // trackX < 28.6) m_effHitmapInterPixelTriple->Fill(true, trackX-28.05, fmod(trackY, m_pitch));
  503.  
  504.  
  505. m_effHitmapInterPixel->Fill(true, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
  506. m_effHitmapInterPixelTriple->Fill(true, trackX-28.05, fmod(trackY, m_pitch));
  507.  
  508.  
  509. if (closestCluster->size() <= m_effHitmapInterPixelVsSizes.size())
  510. m_effHitmapInterPixelVsSizes[closestCluster->size()-1]->Fill(true, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
  511. m_trackAssociated->at(itrack - cutTracks->begin()) = true;
  512. plot(closestCluster->x() - interceptUG.x(), "xResidualsClustMinusTrack", "xResidualsClustMinusTrack", -0.2, 0.2, 400);
  513. plot(closestCluster->y() - interceptUG.y(), "yResidualsClustMinusTrack", "yResidualsClustMinusTrack", -0.2, 0.2, 400);
  514. plot(closestCluster->htime() - (*itrack)->htime(), "hTimeResidualClustMinusTrack", "hTimeResidualClustMinusTrack", -50, 50, 400);
  515. if (closestCluster->size() == 1) {
  516. auto hits = closestCluster->hits();
  517. int col = (*hits.begin())->col();
  518. m_timeResidualVsColumn->fill(col, closestCluster->htime() - (*itrack)->htime());
  519. }
  520. if ((*itrack)->vertexed()) counter("nVerticesCorrelated") += 1;
  521. }
  522. else {
  523. m_effHitmap->Fill(false, trackX, trackY);
  524. // if (trackY > 9.25 && trackY < 10.25) m_effHitmapInterPixel->Fill(false, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
  525. // if (trackX > 28.05 &&
  526. // trackX < 28.6) m_effHitmapInterPixelTriple->Fill(false, trackX-28.05, fmod(trackY, m_pitch));
  527.  
  528. m_effHitmapInterPixel->Fill(false, fmod(trackX, m_pitch), fmod(trackY, m_pitch));
  529. m_effHitmapInterPixelTriple->Fill(false, trackX-28.05, fmod(trackY, m_pitch));
  530. m_effX->Fill(false, trackX);
  531. m_effY->Fill(false, trackY);
  532. //std::cout<<"Inefficiency at time: "<<(*itrack)->htime()<<std::endl;
  533. }
  534. if (fmod(trackX, m_pitch)<0.0183 && fmod(trackY, m_pitch) < 0.0183) m_nTracksCorner++;
  535. if (fmod(trackX, m_pitch) > 0.0183 &&
  536. fmod(trackX, m_pitch) < (2*0.0183) &&
  537. fmod(trackY, m_pitch) > 0.0183 &&
  538. fmod(trackY, m_pitch) < (2*0.0183)) m_nTracksCentral++;
  539. }
  540. }
  541.  
  542.  
  543.  
  544. //=============================================================================
  545. //
  546. //=============================================================================
  547. double TbEffPur::getRadialSeparation(LHCb::TbCluster * cluster,
  548. LHCb::TbTrack * track)
  549. {
  550. Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
  551. double xResidual = cluster->x() - trackIntercept.x();
  552. double yResidual = cluster->y() - trackIntercept.y();
  553.  
  554. double r2 = xResidual*xResidual + yResidual*yResidual;
  555. return pow(r2, 0.5);
  556. }
  557.  
  558.  
  559. //=============================================================================
  560. // Local cuts
  561. //=============================================================================
  562. bool TbEffPur::matchTrackToCluster(LHCb::TbCluster * cluster,
  563. LHCb::TbTrack * track)
  564. {
  565. Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
  566. double xResidual = cluster->x() - trackIntercept.x();
  567. double yResidual = cluster->y() - trackIntercept.y();
  568. double tResidual = cluster->htime() - track->htime();
  569. double delr = pow(pow(xResidual, 2) + pow(yResidual, 2), 0.5);
  570.  
  571. // if (track->vertexed()) plot(xResidual, "vertexResiduals/X", "vertexResiduals/X", -0.5, 0.5, 200);
  572. // if (track->vertexed()) plot(yResidual, "vertexResiduals/Y", "vertexResiduals/Y", -0.5, 0.5, 200);
  573. // if (track->vertexed()) plot(tResidual, "vertexResiduals/T", "vertexResiduals/T", -0.5, 0.5, 200);
  574. //
  575. // if (track->vertexed() && litPixel(cluster, track)) {
  576. // plot(xResidual, "vertexResiduals/Xlit", "vertexResiduals/Xlit", -0.5, 0.5, 200);
  577. // plot(yResidual, "vertexResiduals/Ylit", "vertexResiduals/Ylit", -0.5, 0.5, 200);
  578. // }
  579.  
  580.  
  581. if (fabs(tResidual) > m_tResidualCut) {
  582. counter("nTimeRejected") += 1;
  583. return false;
  584. }
  585.  
  586. if (delr > m_rResidualCut) {
  587. counter("nSpatialRejected") += 1;
  588. if (litPixel(cluster, track)) return true;
  589. else {
  590. counter("nSpatialAndLitRejected") += 1;
  591. return false;
  592. }
  593. }
  594.  
  595. // if (!litPixel(cluster, track)) return false;
  596.  
  597. return true;
  598. }
  599.  
  600.  
  601. //=============================================================================
  602. //
  603. //=============================================================================
  604. bool TbEffPur::interceptDeadPixel(Gaudi::XYZPoint trackInterceptGlobal)
  605. {
  606. // Find the row and column corresponding to the track intercept.
  607. Gaudi::XYZPoint trackInterceptLocal = geomSvc()->globalToLocal(trackInterceptGlobal, m_DUTindex);
  608. int row = trackInterceptLocal.y()/m_pitch - 0.5;
  609. int col = trackInterceptLocal.x()/m_pitch - 0.5;
  610.  
  611. for (int icol = col - 1; icol != col+2; icol++) {
  612. for (int irow = row - 1; irow != row+2; irow++) {
  613. if (icol >= 0 && icol <256 && irow>=0 && irow<256) {
  614. if (pixelSvc()->isMasked(pixelSvc()->address(icol, irow), m_DUTindex)) {
  615. return true;
  616. }
  617.  
  618. // if (icol == 17 && irow == 36) {
  619. // return true;
  620. // }
  621. // if (icol == 18 && irow == 36) {
  622. // return true;
  623. // }
  624. // if (icol == 53 && irow == 3) {
  625. // return true;
  626. // }
  627. // if (icol == 54 && irow == 3) {
  628. // return true;
  629. // }
  630. // if (icol == 73 && irow == 26) {
  631. // return true;
  632. // }
  633. // if (icol == 74 && irow == 26) {
  634. // return true;
  635. // }
  636. // if (icol == 19 && irow == 103) {
  637. // return true;
  638. // }
  639. // if (icol == 31 && irow == 106) {
  640. // return true;
  641. // }
  642. // if (icol == 32 && irow == 106) {
  643. // return true;
  644. // }
  645. // if (icol == 38 && irow == 108) {
  646. // return true;
  647. // }
  648. // if (icol == 63 && irow == 95) {
  649. // return true;
  650. // }
  651. // if (icol == 64 && irow == 95) {
  652. // return true;
  653. // }
  654. // if (icol == 103 && irow == 23) {
  655. // return true;
  656. // }
  657. // if (icol == 104 && irow == 23) {
  658. // return true;
  659. // }
  660. // if (icol == 115 && irow == 20) {
  661. // return true;
  662. // }
  663. // if (icol == 116 && irow == 94) {
  664. // return true;
  665. // }
  666. }
  667. }
  668. }
  669.  
  670.  
  671.  
  672. return false;
  673. }
  674.  
  675. //=============================================================================
  676. //
  677. //=============================================================================
  678. bool TbEffPur::litPixel(LHCb::TbCluster * cluster, LHCb::TbTrack * track)
  679. {
  680. // Find the row and column corresponding to the track intercept.
  681. Gaudi::XYZPoint trackIntercept = geomSvc()->intercept(track, m_DUTindex);
  682. Gaudi::XYZPoint pLocal = geomSvc()->globalToLocal(trackIntercept, m_DUTindex);
  683. int row = pLocal.y()/m_pitch - 0.5;
  684. int col = pLocal.x()/m_pitch - 0.5;
  685.  
  686. // See if within a pre-defined square of pixels.
  687. for (unsigned int i=0; i<cluster->size(); i++) {
  688. unsigned int delRow = abs(row - int(cluster->hits()[i]->row()));
  689. unsigned int delCol = abs(col - int(cluster->hits()[i]->col()));
  690.  
  691. if (delRow <= m_litSquareSide && delCol <= m_litSquareSide) {
  692. counter("nLitPixels") += 1;
  693. return true;
  694. }
  695. }
  696.  
  697. return false;
  698. }
  699.  
  700.  
  701. //=============================================================================
  702. // Excluding dead regions.
  703. //=============================================================================
  704. bool TbEffPur::globalCutPosition(Gaudi::XYZPoint pGlobal)
  705. {
  706. if (pGlobal.x() < m_xLow + m_edgeVetoDistance) return false;
  707. else if (pGlobal.x() > m_xUp - m_edgeVetoDistance) return false;
  708. else if (pGlobal.y() < m_yLow + m_edgeVetoDistance) return false;
  709. else if (pGlobal.y() > m_yUp - m_edgeVetoDistance) return false;
  710. return true; // Inside.
  711. }
  712.  
  713.  
  714. //=============================================================================
  715. // Excluding dead regions.
  716. //=============================================================================
  717. bool TbEffPur::outsideDUT(Gaudi::XYZPoint interceptUG)
  718. {
  719. const auto interceptUL = geomSvc()->globalToLocal(interceptUG, m_DUTindex);
  720. if (interceptUL.x() < 0 + m_edgeVetoDistance ||
  721. interceptUL.x() > 14.08 - m_edgeVetoDistance ||
  722. interceptUL.y() < 0 + m_edgeVetoDistance||
  723. interceptUL.y() > 14.08 - m_edgeVetoDistance) return true;
  724. return false;
  725. }
  726.  
  727.  
  728.  
  729. //=============================================================================
  730. // Viewer output.
  731. //=============================================================================
  732. //=============================================================================
  733. // Viewer output.
  734. //=============================================================================
  735. //=============================================================================
  736. // Viewer output.
  737. //=============================================================================
  738.  
  739. void TbEffPur::outputDeadRegion(unsigned int col, unsigned int row) {
  740. std::ofstream myfile;
  741. myfile.open("/afs/cern.ch/user/d/dsaunder/KeplerViewerData.dat", std::ofstream::app);
  742. std::string outputLine;
  743. for (int irow = std::max(int(row-m_deadAreaRadius), 0);
  744. irow < std::min(int(row+m_deadAreaRadius+1), m_nDUTpixels); irow++) {
  745. for (int icol = std::max(int(col-m_deadAreaRadius), 0);
  746. icol < std::min(int(col+m_deadAreaRadius+1), m_nDUTpixels); icol++) {
  747.  
  748.  
  749. outputLine = "DeadPixel ";
  750. const double xLocal = (col-icol) * m_pitch; // Dont want the middle!
  751. const double yLocal = (row-irow) * m_pitch; // Dont want the middle!
  752. Gaudi::XYZPoint pLocal(xLocal, yLocal, 0.);
  753.  
  754. Gaudi::XYZPoint posn = geomSvc()->localToGlobal(pLocal, m_DUTindex);
  755. outputLine += std::to_string(posn.x()) + " ";
  756. outputLine += std::to_string(posn.y()) + " ";
  757. outputLine += std::to_string(posn.z()) + " ";
  758.  
  759. Gaudi::XYZPoint posn2(pLocal.x() + 0.055, pLocal.y(), 0.);
  760. posn = geomSvc()->localToGlobal(posn2, m_DUTindex);
  761. outputLine += std::to_string(posn.x()) + " ";
  762. outputLine += std::to_string(posn.y()) + " ";
  763. outputLine += std::to_string(posn.z()) + " ";
  764.  
  765. Gaudi::XYZPoint posn3(pLocal.x() + 0.055, pLocal.y()+0.055, 0.);
  766. posn = geomSvc()->localToGlobal(posn3, m_DUTindex);
  767. outputLine += std::to_string(posn.x()) + " ";
  768. outputLine += std::to_string(posn.y()) + " ";
  769. outputLine += std::to_string(posn.z()) + " ";
  770.  
  771. Gaudi::XYZPoint posn4(pLocal.x(), pLocal.y()+0.055, 0.);
  772. posn = geomSvc()->localToGlobal(posn4, m_DUTindex);
  773. outputLine += std::to_string(posn.x()) + " ";
  774. outputLine += std::to_string(posn.y()) + " ";
  775. outputLine += std::to_string(posn.z()) + " ";
  776.  
  777. outputLine += "\n";
  778. myfile << outputLine;
  779. }
  780. }
  781. myfile.close();
  782. }
  783.  
  784.  
  785. //=============================================================================
  786. // Viewer output.
  787. //=============================================================================
  788.  
  789. void TbEffPur::outputViewerData()
  790. {
  791. std::ofstream myfile;
  792. myfile.open("/afs/cern.ch/user/d/dsaunder/KeplerViewerData.dat", std::ofstream::app);
  793. std::string outputLine;
  794. //
  795. // // First output the chips.
  796. // for (unsigned int i=0; i<m_nPlanes; i++) {
  797. // outputLine = "Chip ";
  798. // Gaudi::XYZPoint posn1(0., 14.08, 0.);
  799. // Gaudi::XYZPoint posn = geomSvc()->localToGlobal(posn1, i);
  800. // outputLine += std::to_string(posn.x()) + " ";
  801. // outputLine += std::to_string(posn.y()) + " ";
  802. // outputLine += std::to_string(posn.z()) + " ";
  803. //
  804. // Gaudi::XYZPoint posn2(14.08, 14.08, 0.);
  805. // posn = geomSvc()->localToGlobal(posn2, i);
  806. // outputLine += std::to_string(posn.x()) + " ";
  807. // outputLine += std::to_string(posn.y()) + " ";
  808. // outputLine += std::to_string(posn.z()) + " ";
  809. //
  810. // Gaudi::XYZPoint posn3(14.08, 0., 0.);
  811. // posn = geomSvc()->localToGlobal(posn3, i);
  812. // outputLine += std::to_string(posn.x()) + " ";
  813. // outputLine += std::to_string(posn.y()) + " ";
  814. // outputLine += std::to_string(posn.z()) + " ";
  815. //
  816. // Gaudi::XYZPoint posn4(0., 0., 0.);
  817. // posn = geomSvc()->localToGlobal(posn4, i);
  818. // outputLine += std::to_string(posn.x()) + " ";
  819. // outputLine += std::to_string(posn.y()) + " ";
  820. // outputLine += std::to_string(posn.z()) + " ";
  821. //
  822. // outputLine += "\n";
  823. // myfile << outputLine;
  824. // }
  825.  
  826.  
  827. // Clusters.
  828. auto ic = m_clusters->begin();
  829. const auto end = m_clusters->end();
  830. for (; ic != end; ++ic) {
  831. outputLine = "Cluster ";
  832. outputLine += std::to_string((*ic)->x()) + " ";
  833. outputLine += std::to_string((*ic)->y()) + " ";
  834. outputLine += std::to_string((*ic)->z()) + " ";
  835. outputLine += std::to_string((*ic)->htime()) + " ";
  836.  
  837. // if ((*ic)->endCluster() && (*ic)->vertexed()) outputLine += "5 \n";
  838. //if ((*ic)->vertexed()) outputLine += "4 \n";
  839. // else if ((*ic)->endCluster()) outputLine += "3 \n";
  840. if ((*ic)->associated()) outputLine += "2 \n";
  841. // else if ((*ic)->volumed()) outputLine += "1 \n";
  842.  
  843. else {
  844. outputLine += "0 \n";
  845. }
  846. myfile << outputLine;
  847. }
  848.  
  849. outputLine = "CentralRegion ";
  850. outputLine += std::to_string(m_xLow) + " ";
  851. outputLine += std::to_string(m_xUp) + " ";
  852. outputLine += std::to_string(m_yLow) + " ";
  853. outputLine += std::to_string(m_yUp) + " ";
  854. outputLine += std::to_string(geomSvc()->module(m_DUTindex)->z()) + " ";
  855. myfile << outputLine;
  856.  
  857. // // Its hits.
  858. // for (auto hit : (*ic)->hits()) {
  859. // outputLine = "Pixel ";
  860. // const double xLocal = (hit->col()) * m_pitch; // Dont want the middle!
  861. // const double yLocal = (hit->row()) * m_pitch; // Dont want the middle!
  862. // Gaudi::XYZPoint pLocal(xLocal, yLocal, 0.);
  863. //
  864. // Gaudi::XYZPoint posn = geomSvc()->localToGlobal(pLocal, i);
  865. // outputLine += std::to_string(posn.x()) + " ";
  866. // outputLine += std::to_string(posn.y()) + " ";
  867. // outputLine += std::to_string(posn.z()) + " ";
  868. //
  869. // Gaudi::XYZPoint posn2(pLocal.x() + 0.055, pLocal.y(), 0.);
  870. // posn = geomSvc()->localToGlobal(posn2, i);
  871. // outputLine += std::to_string(posn.x()) + " ";
  872. // outputLine += std::to_string(posn.y()) + " ";
  873. // outputLine += std::to_string(posn.z()) + " ";
  874. //
  875. // Gaudi::XYZPoint posn3(pLocal.x() + 0.055, pLocal.y()+0.055, 0.);
  876. // posn = geomSvc()->localToGlobal(posn3, i);
  877. // outputLine += std::to_string(posn.x()) + " ";
  878. // outputLine += std::to_string(posn.y()) + " ";
  879. // outputLine += std::to_string(posn.z()) + " ";
  880. //
  881. // Gaudi::XYZPoint posn4(pLocal.x(), pLocal.y()+0.055, 0.);
  882. // posn = geomSvc()->localToGlobal(posn4, i);
  883. // outputLine += std::to_string(posn.x()) + " ";
  884. // outputLine += std::to_string(posn.y()) + " ";
  885. // outputLine += std::to_string(posn.z()) + " ";
  886. //
  887. // outputLine += std::to_string(hit->htime()) + " ";
  888. // outputLine += std::to_string(hit->ToT()) + " ";
  889. //
  890. // outputLine += "\n";
  891. // myfile << outputLine;
  892. // }
  893.  
  894. myfile.close();
  895. }
  896.  
  897.  
  898. //=============================================================================
  899. // END
  900. //=============================================================================