Newer
Older
TB_Chris / TbAnalysis / src / .svn / text-base / TbTimewalkMonitor.cpp.svn-base
  1. // Gaudi
  2. #include "GaudiKernel/PhysicalConstants.h"
  3. #include "GaudiUtils/HistoLabels.h"
  4. #include "GaudiUtils/Aida2ROOT.h"
  5.  
  6. // Tb/TbEvent
  7. #include "Event/TbTrack.h"
  8. #include "Event/TbCluster.h"
  9.  
  10. // Tb/TbKernel
  11. #include "TbKernel/TbConstants.h"
  12. #include "TbKernel/TbModule.h"
  13.  
  14. // Local
  15. #include "TbTimewalkMonitor.h"
  16.  
  17. // ROOT
  18. #include "TH1D.h"
  19.  
  20. using namespace Gaudi::Utils::Histos;
  21.  
  22. DECLARE_ALGORITHM_FACTORY(TbTimewalkMonitor)
  23.  
  24. //=============================================================================
  25. // Standard constructor
  26. //=============================================================================
  27. TbTimewalkMonitor::TbTimewalkMonitor(const std::string& name,
  28. ISvcLocator* pSvcLocator)
  29. : TbAlgorithm(name, pSvcLocator) {
  30.  
  31. declareProperty("TrackLocation",
  32. m_trackLocation = LHCb::TbTrackLocation::Default);
  33. declareProperty("WidthMin", m_widthMin = 37 );
  34. declareProperty("WidthMax", m_widthMax = 42 );
  35. }
  36.  
  37. //=============================================================================
  38. // Initialization
  39. //=============================================================================
  40. StatusCode TbTimewalkMonitor::initialize() {
  41.  
  42. // Initialise the base class.
  43. StatusCode sc = TbAlgorithm::initialize();
  44. if (sc.isFailure()) return sc;
  45.  
  46. // Book the histograms.
  47. for (unsigned int i = 0; i < m_nPlanes; ++i) {
  48. const std::string plane = std::to_string(i);
  49. const std::string title = geomSvc()->modules().at(i)->id();
  50. m_space.push_back(bookProfile1D("SpatialResidual/Plane"+plane,title,-0.5,255.5,256));
  51.  
  52. m_spd2D.push_back(book2D("spd2D/Plane"+plane,title,-0.5,255.5,256,-100,100,20));
  53. m_twd.push_back(book2D("twd/Plane"+plane,title,1.,101.,100,0.,62.5,40));
  54. m_twdQ.push_back(book2D("twdQ/Plane"+plane,title,1000,10000,100,0.,62.5,40));
  55. m_twdQL.push_back(book2D("twdQL/Plane"+plane,title,1000,10000,100,0.,62.5,40));
  56. m_twdQR.push_back(book2D("twdQR/Plane"+plane,title,1000,10000,100,0.,62.5,40));
  57.  
  58. setAxisLabels( m_twd[i], "ToT", "#Delta t [ns]");
  59.  
  60. m_twdL.push_back(book2D("twdL/Plane"+plane,title,0.,100.,100,0.,62.5,40));
  61. m_twdR.push_back(book2D("twdR/Plane"+plane,title,0.,100.,100,0.,62.5,40));
  62. m_LRSYNC.emplace_back(bookProfile1D("LR_SYNC/Plane"+plane,title,-0.5,255.5,256));
  63. m_UDSYNC.emplace_back(bookProfile1D("UD_SYNC/Plane"+plane,title,-0.5,255.5,256));
  64.  
  65. m_quad.emplace_back(bookProfile1D("QUAD/Plane"+plane,title,-0.5,255.5,256));
  66.  
  67. m_inscol.emplace_back(bookProfile1D("Scol/Plane"+plane,title,-0.5,127.5,128));
  68. m_interscol.emplace_back(bookProfile1D("InterScol/Plane"+plane,title,-0.5,127.5,128));
  69.  
  70. m_timewalk.push_back(bookProfile1D("Timewalk/Plane"+plane,title, 1., 401., 400));
  71. m_dtDist.push_back( book1D( "dtDist/Plane"+plane,title, 0., 90., 58 ) );
  72. m_cDist.push_back( book1D( "cDist/Plane"+plane,title,0.,90.,58 ) );
  73. m_timewalkQ.push_back(bookProfile1D("TimewalkQ/Plane"+plane,title,0,20000,200));
  74. setAxisLabels(m_timewalkQ[i], "charge [e^{-}]","#Delta t [ns]");
  75. setAxisLabels(m_timewalk[i], "ToT", "#Delta t [ns]");
  76. }
  77.  
  78. return StatusCode::SUCCESS;
  79. }
  80.  
  81. StatusCode TbTimewalkMonitor::finalize() {
  82.  
  83. for( unsigned int i = 0 ; i < m_nPlanes; ++i ){
  84. TH1D* tdist = Gaudi::Utils::Aida2ROOT::aida2root( m_dtDist[i] );
  85. AIDA::IHistogram1D* cdist = m_cDist[i];
  86. int total = 0;
  87. for (int i=0;i<tdist->GetNbinsX();i++)
  88. {
  89. total += tdist->GetBinContent(i);
  90. cdist->fill( tdist->GetBinCenter(i) , total /tdist->GetEntries() ) ;
  91. }
  92. }
  93. return StatusCode::SUCCESS;
  94. }
  95.  
  96. //=============================================================================
  97. // Main execution
  98. //=============================================================================
  99. StatusCode TbTimewalkMonitor::execute() {
  100.  
  101. // Grab the tracks.
  102.  
  103. const LHCb::TbTracks* tracks = getIfExists<LHCb::TbTracks>(m_trackLocation);
  104. if (!tracks) return StatusCode::SUCCESS;
  105. for (const LHCb::TbTrack* track : *tracks) {
  106. // info() << track << endmsg;
  107. const SmartRefVector<LHCb::TbCluster>& clusters = track->clusters();
  108. const auto ref = std::min_element(
  109. clusters.begin(), clusters.end(),
  110. [](const LHCb::TbCluster* h1,
  111. const LHCb::TbCluster* h2) { return h1->htime() < h2->htime(); });
  112. const double tref = (*ref)->htime();
  113. const double syncRef = (*clusters.rbegin())->htime();
  114. unsigned int p = (*clusters.rbegin())->plane();
  115.  
  116. for (auto it = clusters.cbegin(), end = clusters.cend(); it != end; ++it){
  117.  
  118. auto hits = (*it)->hits();
  119. for (const auto hit : hits) {
  120. if( p==5 ){
  121. m_space[(*it)->plane()]->fill( hit->col(), hit->htime() - syncRef );
  122. m_spd2D[(*it)->plane()]->fill( hit->col(), hit->htime() - syncRef );
  123. }
  124. m_timewalk[(*it)->plane()]->fill(hit->ToT(), hit->htime() - tref);
  125. m_twd[(*it)->plane()]->fill(hit->ToT(), hit->htime() - tref);
  126. m_timewalkQ[(*it)->plane()]->fill(hit->charge(), hit->htime() - tref );
  127. m_twdQ[(*it)->plane()]->fill(hit->charge(), hit->htime() - tref );
  128. m_dtDist[(*it)->plane()]->fill( hit->htime() - tref );
  129. }
  130. if( (*it)->size() == 2 ){
  131. auto minmax = std::minmax_element(hits.begin(), hits.end(),[](LHCb::TbHit* h1, LHCb::TbHit* h2) { return h1->col() < h2->col(); });
  132. if( (*minmax.first)->col() != (*minmax.second)->col() ){
  133. m_twdL[(*it)->plane()]->fill( (*minmax.first)->ToT(), (*minmax.first)->htime() - tref );
  134. m_twdR[(*it)->plane()]->fill( (*minmax.second)->ToT(), (*minmax.second)->htime() - tref );
  135. m_twdQL[(*it)->plane()]->fill( (*minmax.first)->charge(), (*minmax.first)->htime() - tref );
  136. m_twdQR[(*it)->plane()]->fill( (*minmax.second)->charge(), (*minmax.second)->htime() - tref );
  137. }
  138. }
  139. }
  140. const SmartRefVector<LHCb::TbCluster>& associatedClusters = track->associatedClusters();
  141.  
  142. for (auto it = associatedClusters.cbegin(), end = associatedClusters.cend(); it != end; ++it){
  143. auto hits = (*it)->hits();
  144. //info() << "Filling associated clusters for " << track << endmsg;
  145. for (const auto hit : (*it)->hits()) {
  146. m_timewalk[(*it)->plane()]->fill(hit->ToT(), hit->htime() - tref);
  147. m_twd[(*it)->plane()]->fill(hit->ToT(), hit->htime() - tref);
  148. //m_twdQ[(*it)->plane()]->fill(hit->charge(), hit->htime() - tref );
  149. //m_dtDist[(*it)->plane()]->fill( hit->htime() - tref );
  150. m_timewalkQ[(*it)->plane()]->fill(hit->charge(), hit->htime() - tref );
  151. //if( (*it)->size() == 1 ) m_timewalkOneHit[(*it)->plane()]->fill(hit->ToT(), ( hit->htime() - tref ));
  152. //if( (*it)->size() == 2 ) m_timewalkTwoHit[(*it)->plane()]->fill(hit->ToT(), ( hit->htime() - tref ));
  153. }
  154. /*
  155. if( (*it)->size() == 2 ){
  156. auto minmax = std::minmax_element(hits.begin(), hits.end(),[](LHCb::TbHit* h1, LHCb::TbHit* h2) { return h1->col() < h2->col(); });
  157. if( (*minmax.first)->col() != (*minmax.second)->col() ){
  158. m_twdL[(*it)->plane()]->fill( (*minmax.first)->ToT(), (*minmax.first)->htime() - tref );
  159. m_twdR[(*it)->plane()]->fill( (*minmax.second)->ToT(), (*minmax.second)->htime() - tref );
  160. m_twdQL[(*it)->plane()]->fill( (*minmax.first)->charge(), (*minmax.first)->htime() - tref );
  161. m_twdQR[(*it)->plane()]->fill( (*minmax.second)->charge(), (*minmax.second)->htime() - tref );
  162. }
  163. }
  164. */
  165. //info() << "Filled associated clusters " << endmsg;
  166. }
  167. }
  168.  
  169. return StatusCode::SUCCESS;
  170. }
  171. /*
  172. for(unsigned int plane=0;plane<m_nPlanes; ++plane ){
  173. const LHCb::TbClusters* clusters =
  174. getIfExists<LHCb::TbClusters>(LHCb::TbClusterLocation::Default+std::to_string(plane));
  175. for( auto& cluster : *clusters ){
  176. if( cluster->size() == 2 ){
  177. auto hits = cluster->hits();
  178. auto minmax = std::minmax_element(hits.begin(), hits.end(),
  179. [](LHCb::TbHit* h1, LHCb::TbHit* h2) { return h1->col() < h2->col(); });
  180. const LHCb::TbHit* left = (*minmax.first);
  181. const LHCb::TbHit* right = (*minmax.second);
  182. if( left->col() != right->col() ){
  183. m_LRSYNC[plane]->fill( left->col() , left->htime() - right->htime() );
  184. if( int(left->col() /2 ) == int(right->col() /2 ) )
  185. m_inscol[plane]->fill( int(left->col()/2) , left->htime() - right->htime() );
  186. else
  187. m_interscol[plane]->fill( int(left->col()/2) , left->htime() - right->htime() );
  188. }
  189. else {
  190. minmax = std::minmax_element(hits.begin(), hits.end(),
  191. [](LHCb::TbHit* h1, LHCb::TbHit* h2) { return h1->row() < h2->row(); });
  192. m_UDSYNC[plane]->fill( (*minmax.first)->row(),
  193. (*minmax.first)->htime() - (*minmax.second)->htime());
  194. }
  195. }
  196. if( cluster->size() == 4 ){
  197. auto hits = cluster->hits();
  198. auto minmax = std::minmax_element(hits.begin(), hits.end(),
  199. [](LHCb::TbHit* h1, LHCb::TbHit* h2) { return h1->col() < h2->col(); });
  200. const LHCb::TbHit* left = (*minmax.first);
  201. const LHCb::TbHit* right = (*minmax.second);
  202. if( left->col() == right->col() + 4 ){
  203. m_quad[plane]->fill( left->col() , left->htime() - right->htime() );
  204. }
  205. }
  206.  
  207. }
  208. }
  209. return StatusCode::SUCCESS;
  210. }
  211. */