Newer
Older
TB_Chris / TbUT / src / .svn / text-base / TbUTCMSIterativelyPerBeetle.cpp.svn-base
  1. #include "TbUTCMSIterativelyPerBeetle.h"
  2. #include <iostream>
  3. #include <cmath>
  4.  
  5. using namespace TbUT;
  6. using namespace std;
  7.  
  8. CMSIterativelyPerBeetle::CMSIterativelyPerBeetle(IChannelMaskProvider& p_masksProvider):
  9. m_masksProvider(p_masksProvider),
  10. m_channelNumber(RawData<>::getnChannelNumber()),
  11. m_channelPerBeetle(32)
  12. {
  13. cout<<"create iteratively correlation per beetle"<<endl;
  14. initializeCorrectionAndHitMaps();
  15. }
  16.  
  17.  
  18. void CMSIterativelyPerBeetle::initializeCorrectionAndHitMaps()
  19. {
  20. for (int channel=0;channel<m_channelNumber;channel+=m_channelPerBeetle){
  21. m_correctionPerBeetle.insert(std::make_pair(channel,0));
  22. m_hitThresholdPerBeetle.insert(std::make_pair(channel,200));
  23. }
  24. }
  25.  
  26. void CMSIterativelyPerBeetle::processEvent(RawData<>* p_data, RawData<double> **p_output)
  27. {
  28. int numberOfIteration=2;
  29. RawData<double>* tmpData;
  30. // first iteration
  31. calculateCorrection<int>(p_data);
  32. removeCM<int, double>(p_data,p_output);
  33. // rest of the iterations
  34. for(int iteration=1;iteration<numberOfIteration; iteration++){
  35. tmpData=new RawData<double>(**p_output);
  36. calculateCorrection<double>(tmpData);
  37. removeCM<double, double>(tmpData,p_output);
  38. }
  39. resetHitThresholds();
  40. }
  41.  
  42. template<typename DATA_TYPE>
  43. void CMSIterativelyPerBeetle::calculateCorrection(RawData<DATA_TYPE>* p_inputData)
  44. {
  45. for(auto& mapIt : m_correctionPerBeetle)
  46. {
  47. int usedChannels=0;
  48. double rmsPerBeetle=0;
  49. for(int channel=mapIt.first;channel<mapIt.first+m_channelPerBeetle;channel++)
  50. {
  51. double signal=static_cast<double>(p_inputData->getSignal(channel));
  52. if(!m_masksProvider.isMasked(channel) && abs(signal) <m_hitThresholdPerBeetle[mapIt.first]){
  53. mapIt.second+=signal;
  54. rmsPerBeetle+=signal*signal;
  55. usedChannels++;
  56. }
  57. }
  58. if(usedChannels) mapIt.second/=static_cast<double>(usedChannels);
  59. if(usedChannels) rmsPerBeetle/=static_cast<double>(usedChannels);
  60. rmsPerBeetle-=mapIt.second*mapIt.second;
  61. double rmsMultiplicity=4;
  62. m_hitThresholdPerBeetle[mapIt.first]=rmsMultiplicity*sqrt(rmsPerBeetle);
  63. }
  64. }
  65.  
  66. template<typename INPUT_DATA_TYPE, typename OUTPUT_TYPE_NAME>
  67. void CMSIterativelyPerBeetle::removeCM(RawData<INPUT_DATA_TYPE>* p_data, RawData<OUTPUT_TYPE_NAME> **p_output)
  68. {
  69. for(auto& mapIt : m_correctionPerBeetle)
  70. {
  71. for(int channel=mapIt.first;channel<mapIt.first+m_channelPerBeetle;channel++)
  72. {
  73. if(m_masksProvider.isMasked(channel)){
  74. OUTPUT_TYPE_NAME signalMaskedChannel=0;
  75. (*p_output)->setSignal(signalMaskedChannel);
  76. }else{
  77. OUTPUT_TYPE_NAME l_channelSignal=p_data->getSignal(channel)-mapIt.second;
  78. (*p_output)->setSignal(l_channelSignal);
  79. }
  80. }
  81.  
  82. }
  83. (*p_output)->setTDC(p_data->getTDC());
  84. (*p_output)->setTime(p_data->getTime());
  85.  
  86. }
  87.  
  88. void CMSIterativelyPerBeetle::resetHitThresholds()
  89. {
  90. for(auto& mapIt : m_hitThresholdPerBeetle) mapIt.second=200;
  91. }