Newer
Older
TB_Chris / TbUT / src / TbUTCMSIterativelyPerBeetle.cpp
#include "TbUTCMSIterativelyPerBeetle.h"
#include <iostream>
#include <cmath>

using namespace TbUT;
using namespace std;

CMSIterativelyPerBeetle::CMSIterativelyPerBeetle(IChannelMaskProvider& p_masksProvider):
        m_masksProvider(p_masksProvider),
        m_channelNumber(RawData<>::getnChannelNumber()),
		m_channelPerBeetle(32)
{
    cout<<"create iteratively correlation per beetle"<<endl;
    initializeCorrectionAndHitMaps();
}


void CMSIterativelyPerBeetle::initializeCorrectionAndHitMaps()
{
	for (int channel=0;channel<m_channelNumber;channel+=m_channelPerBeetle){
		m_correctionPerBeetle.insert(std::make_pair(channel,0));
		m_hitThresholdPerBeetle.insert(std::make_pair(channel,200));
	}
}

void CMSIterativelyPerBeetle::processEvent(RawData<>* p_data, RawData<double> **p_output)
{
	int numberOfIteration=2;
	RawData<double>* tmpData;
	// first iteration
	calculateCorrection<int>(p_data);
	removeCM<int, double>(p_data,p_output);
	// rest of the iterations
	for(int iteration=1;iteration<numberOfIteration; iteration++){
		tmpData=new RawData<double>(**p_output);
		calculateCorrection<double>(tmpData);
		removeCM<double, double>(tmpData,p_output);
	}
	resetHitThresholds();
}

template<typename DATA_TYPE>
void CMSIterativelyPerBeetle::calculateCorrection(RawData<DATA_TYPE>* p_inputData)
{
	for(auto& mapIt : m_correctionPerBeetle)
	{
		int usedChannels=0;
		double rmsPerBeetle=0;
		for(int channel=mapIt.first;channel<mapIt.first+m_channelPerBeetle;channel++)
		{
			double signal=static_cast<double>(p_inputData->getSignal(channel));
			if(!m_masksProvider.isMasked(channel) && abs(signal) <m_hitThresholdPerBeetle[mapIt.first]){
				mapIt.second+=signal;
				rmsPerBeetle+=signal*signal;
				usedChannels++;
			}
		}
		if(usedChannels) mapIt.second/=static_cast<double>(usedChannels);
		if(usedChannels) rmsPerBeetle/=static_cast<double>(usedChannels);
		rmsPerBeetle-=mapIt.second*mapIt.second;
		double rmsMultiplicity=4;
		m_hitThresholdPerBeetle[mapIt.first]=rmsMultiplicity*sqrt(rmsPerBeetle);
	}
}

template<typename INPUT_DATA_TYPE, typename OUTPUT_TYPE_NAME>
void CMSIterativelyPerBeetle::removeCM(RawData<INPUT_DATA_TYPE>* p_data, RawData<OUTPUT_TYPE_NAME> **p_output)
{
	for(auto& mapIt : m_correctionPerBeetle)
	{
		for(int channel=mapIt.first;channel<mapIt.first+m_channelPerBeetle;channel++)
		{
			if(m_masksProvider.isMasked(channel)){
				OUTPUT_TYPE_NAME signalMaskedChannel=0;
				(*p_output)->setSignal(signalMaskedChannel);
			}else{
				OUTPUT_TYPE_NAME l_channelSignal=p_data->getSignal(channel)-mapIt.second;
				(*p_output)->setSignal(l_channelSignal);
			}
		}

	}
	 (*p_output)->setTDC(p_data->getTDC());
	 (*p_output)->setTime(p_data->getTime());

}

void CMSIterativelyPerBeetle::resetHitThresholds()
{
	for(auto& mapIt : m_hitThresholdPerBeetle) mapIt.second=200;
}