CPulseFitter.C
//-----------------------------------------------------------------------------
// $Header$
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class ...
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include <strstream.h>
#include <fstream.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "CCoGPulse.h"
#include "CPulseFitter.h"
CPulseFitter::CPulseFitter(double s1, double s2, double s3)
{
for(int i=0; i<=255; i++) derivative[i]=0;
minRegionLength = 4;
epsilon = 0.01;
deltaMax = 10.0 ;
maxIterations = 30;
calibratedSigma0 = s1; //96 data 1.889; //95 data 1.5725
calibratedSigma1 = s2; //96 data 0.8489e-3; //95 data 0.90795e-3
calibratedSigmaDeviation = s3; //96 data 2.058; //95 data 1.8419
//96 data 0.603e-3; //95 data 0.42846e-3
}
CPulseFitter::~CPulseFitter()
{
}
CPulseFitter & CPulseFitter::operator = (const CPulseFitter &)
{
return *this;
}
CBoolean CPulseFitter::fitAndSplit(CCyclicCollection<CCell>& cells,
CList<CCoGPulse>& pulses, int ipulse)
{
float startPosition[2], startAmplitude[2];
ianode = pulses(ipulse)->getAnode();
tMin = pulses(ipulse)->getTStart();
tMax = pulses(ipulse)->getTStop();
CCoGPulse *newPulse = 0;
int regionLength = tMax-tMin+1;
int thisMax;
int lowBin, higBin;
getMaxima(cells);
for( int i=1; i<=numberOfMaxima; i++) testSaturation(cells, i);
CBoolean tooShort = (regionLength<minRegionLength);
if(numberOfMaxima==1)
{
thisMax=1;
CBoolean chooppedPulse = (maxPos[thisMax]==tMin || maxPos[thisMax]==tMax); //scanner's fault
if(isSaturated[thisMax]==1 || tooShort)
{
getMomenta(cells);
fillPulse(pulses(ipulse), pulseCoG, pulseRMS, pulseArea);
}
else if(param3Regression(cells, tMin, tMax))
{
if(singlePulseMaxSigma(gaussTime0, gaussSigma) || chooppedPulse)
{
fillPulse(pulses(ipulse), gaussTime0, gaussSigma, gaussArea);
}
else
{
getMomenta(cells);
finalSigma[1] =calibratedPulseSigma(maxPos[thisMax]);
finalSigma[2] =finalSigma[1];
// startPosition[1] =gaussTime0;
// startAmplitude[1] =gaussAmpl;
// startPosition[2] =sqrt(abs(pulseArea/(pulseArea-gaussArea)*pulseRMS
// -gaussArea/(pulseArea-gaussArea)*(pow(finalSigma[1],2)
// +pow(gaussTime0-pulseCoG,2))-pow(finalSigma[2],2)))+pulseCoG;
// startAmplitude[2] =(pulseArea-gaussArea)/2.507/finalSigma[2];
int signumOfPulseCMS = pulseCMS>0 ? 1 : -1;
startPosition[1] =gaussTime0-5*(pulseCMS<0)-0.3*signumOfPulseCMS;
startAmplitude[1] =gaussAmpl*(0.9-0.6*(pulseCMS<0));
startPosition[2] =gaussTime0+5.0*(pulseCMS>0)-signumOfPulseCMS;
startAmplitude[2] =gaussAmpl*(0.9-0.6*(pulseCMS>0));
/*********************************************/
doDoubleFit(cells, startPosition, startAmplitude );
/*********************************************/
fillPulse(pulses(ipulse), finalPos[1], finalSigma[1], 2.5*finalSigma[1]*finalAmpl[1]);
newPulse = new CCoGPulse(ianode, tMin);
newPulse->setTStop(tMax);
fillPulse(newPulse, finalPos[2], finalSigma[2], 2.5*finalSigma[2]*finalAmpl[2]);
pulses.insert(newPulse);
}
}
else
{
getMomenta(cells);
fillPulse(pulses(ipulse), pulseCoG, pulseRMS, pulseArea);
}
}
else
if(numberOfMaxima==2)
{
for(int m=1; m<=2; m++){
finalSigma[m] = calibratedPulseSigma(maxPos[m]);
if(maxPos[m]-2 > tMin) lowBin = maxPos[m]-2; else lowBin = tMin;
if(maxPos[m]+2 < tMax) higBin = maxPos[m]+2; else higBin = tMax;
if(param3Regression(cells, lowBin, higBin))
{
startPosition[m] = gaussTime0;
startAmplitude[m] = gaussAmpl;
}
else
{
startPosition[m] = maxPos[m];
startAmplitude[m] = maxAmp[m];
}
}
/*********************************************/
doDoubleFit(cells, startPosition, startAmplitude );
/*********************************************/
fillPulse(pulses(ipulse), finalPos[1], finalSigma[1], 2.5*finalSigma[1]*finalAmpl[1]);
newPulse = new CCoGPulse(ianode, tMin);
newPulse->setTStop(tMax);
fillPulse(newPulse, finalPos[2], finalSigma[2], 2.5*finalSigma[2]*finalAmpl[2]);
pulses.insert(newPulse);
}
else
{
for(int m=1; m<=numberOfMaxima; m++){
finalSigma[1] = calibratedPulseSigma(maxPos[m]);
if(maxPos[m]-2 > tMin) lowBin = maxPos[m]-2; else lowBin = tMin;
if(maxPos[m]+2 < tMax) higBin = maxPos[m]+2; else higBin = tMax;
param3Regression(cells, lowBin, higBin);
if(m==1)
{
fillPulse(pulses(ipulse), gaussTime0, gaussSigma, gaussArea);
}
else
{
newPulse = new CCoGPulse(ianode, tMin);
newPulse->setTStop(tMax);
fillPulse(newPulse, gaussTime0, gaussSigma, gaussArea);
pulses.insert(newPulse);
}
}
}
return True;
}
//---------------------------------------------------------------------------------
// scan shape fo maxima, saturation etc.
//---------------------------------------------------------------------------------
void CPulseFitter::getMaxima(CCyclicCollection<CCell>& cell)
{
for(int j=tMin; j<tMax;j++)
derivative[j] = cell(ianode,j+1)->getAmp() - cell(ianode,j)->getAmp();
derivative[tMax]= -cell(ianode,tMax)->getAmp();
numberOfMaxima=0;
for(j=0; j<=20; j++) { maxPos[j] = 0; maxAmp[j] = 0; isSaturated[j] = 0; }
lowBorder = tMin-1;
upBorder = tMax;
for(j=lowBorder;j<=upBorder;j++){
pdown=j+1;
scanDownVicinity(cell);
while(!(downScan[0]>=0 && downScan[1]<0 && downScan[1]<0) && pdown<=upBorder+1){
pdown++;
scanDownVicinity(cell);
}
pup =pdown-1;
scanUpVicinity(cell);
while(!(upScan[0]>0 && upScan[1]>0 && upScan[2]<=0) && pup>lowBorder && pup>maxPos[numberOfMaxima]){
pup--;
scanUpVicinity(cell);
}
if(pup>maxPos[numberOfMaxima] && pdown<upBorder+2 ){
numberOfMaxima++;
maxPos[numberOfMaxima]=(pup+pdown)/2+1;
maxAmp[numberOfMaxima]=cell(ianode,maxPos[numberOfMaxima])->getAmp();
j=pdown+1;
}
}
}
void CPulseFitter::getMomenta(CCyclicCollection<CCell>& cell)
{
pulseCoG=0.0; pulseRMS=0.0; pulseCMS=0.0; pulseArea=0.0;
double pulseCoG_x=0.0, pulseCoG_xh2=0.0, pulseCoG_xh3=0.0;
float amplitude;
for(int j=tMin;j<=tMax;j++){
amplitude=cell(ianode,j)->getAmp();
pulseArea += amplitude;
pulseCoG_x += j*amplitude;
pulseCoG_xh2 += pow(j,2)*amplitude;
pulseCoG_xh3 += pow(j,3)*amplitude;
}
if(pulseArea>0 && tMax-tMin>=2){
pulseCoG = pulseCoG_x/pulseArea;
pulseRMS = (pulseCoG_xh2 - pow(pulseCoG_x,2)/pulseArea)/pulseArea;
pulseCMS = pulseCoG_xh3/pulseArea - 3*pulseCoG_x*pulseCoG_xh2/pow(pulseArea,2) +
2*pow(pulseCoG_x/pulseArea,3);
} else {
pulseCoG = tMin + (tMax-tMin)/2.0;
pulseRMS = 0.0;
pulseCMS = 0.0;
}
}
void CPulseFitter::scanUpVicinity(CCyclicCollection<CCell>& cell)
{
for(int scan=0;scan<3;scan++) {
int pup_v = pup-1+scan;
if(pup_v<tMin) upScan[scan] = cell(ianode,tMin)->getAmp(); else
if(pup_v>tMax) upScan[scan] = -cell(ianode,tMax)->getAmp(); else
upScan[scan]=derivative[pup_v];
}
}
void CPulseFitter::scanDownVicinity(CCyclicCollection<CCell>& cell)
{
for(int scan=0;scan<3;scan++) {
int pdown_v = pdown-1+scan;
if(pdown_v<tMin) downScan[scan] = cell(ianode,tMin)->getAmp(); else
if(pdown_v>tMax) downScan[scan] = -cell(ianode,tMax)->getAmp(); else
downScan[scan]=derivative[pdown_v];
}
}
void CPulseFitter::testSaturation(CCyclicCollection<CCell>& cell, int posId)
{
int pointer=maxPos[posId]-1;
int noEqual=1;
while(pointer>=tMin && cell(ianode,pointer)->getAmp()==maxAmp[posId]) {
pointer--;
noEqual++;
}
pointer=maxPos[posId]+1;
while(pointer<=tMax && cell(ianode,pointer)->getAmp()==maxAmp[posId]) {
pointer++;
noEqual++;
}
if(noEqual>=3) isSaturated[posId]=1; else isSaturated[posId]=0;
}
//---------------------------------------------------------------------------------
// Gauss regression
//---------------------------------------------------------------------------------
CBoolean CPulseFitter::param3Regression(CCyclicCollection<CCell>& cell, int first, int last)
{
//calculates the pulse parameters using gauss regression method
double s0 = 0.0; double s1 = 0.0; double s2 = 0.0;
double s3 = 0.0; double s4 = 0.0;
double sy = 0.0; double sxy = 0.0;
double sx2y = 0.0; double sy2 = 0.0;
gaussSigma = -1.0;
gaussTime0 = -1.0;
gaussAmpl = -1.0;
for(int i=0;i<=255;i++){w[i]=0; z[i]=0.0; subtrac[i]=0.0;}
for(i=first;i<=last;i++){
w[i]=cell(ianode,i)->getAmp();
if (cell(ianode,i)->getAmp()>0.0)
z[i]=log(cell(ianode,i)->getAmp());
else
z[i]=0;
s1=s1+i*w[i];
s2=s2+pow(i,2)*w[i];
s3=s3+pow(i,3)*w[i];
s4=s4+pow(i,4)*w[i];
sy=sy+z[i]*w[i];
sy2=sy2+pow(z[i],2)*w[i];
sxy=sxy+i*z[i]*w[i];
sx2y=sx2y+pow(i,2)*z[i]*w[i];
s0=s0+w[i];
}
double s1kv=s1*s1;
double s12=s1*s2;
double s13=s1*s3;
double s14=s1*s4;
double s2kv=s2*s2;
double s23=s2*s3;
double s24=s2*s4;
double s3kv=s3*s3;
double del=s0*(s24-s3kv)+s1*(s23-s14)+s2*(s13-s2kv);
if(del>0.0) {
cc=((s24-s3kv)*sy+(s23-s14)*sxy+(s13-s2kv)*sx2y)/del;
bb=((s23-s14)*sy+(s0*s4-s2kv)*sxy+(s12-s0*s3)*sx2y)/del;
aa=((s13-s2kv)*sy+(s12-s0*s3)*sxy+(s0*s2-s1kv)*sx2y)/del;
if (aa<0.0){
gaussSigma=sqrt(-0.5/aa);
gaussTime0=-bb/2.0/aa;
double determ=(cc-pow(bb,2)/4./aa);
if(abs(determ)>1.0e-10 && abs(determ)<50.0) gaussAmpl=exp(determ);
}
}
//--------------------------------------------------------
// third momentum assymetry
// assym=(s3/s0-3*s2*s1/s0**2+2*(s1/s0)**3)/gaussSigma**3x1
//--------------------------------------------------------
if(gaussAmpl<0 || gaussTime0<0 || gaussSigma<0) return False;
else {
gaussArea=2.50663*gaussSigma*gaussAmpl; //sqrt(2*pi)
return True;
}
}
void CPulseFitter::fillPulse(CCoGPulse* pulse, double newPulsePos, double newPulseRMS, double newPulseArea)
{
pulse->setCenter(newPulsePos);
pulse->setRmsTime(newPulseRMS);
pulse->setAmp(newPulseArea); // integral amplitude
pulse->setIsUsedForHit(False); // not yet used for hit
}
//---------------------------------------------------------------------------------
// RobustFit the kvazi-linearisation method
//---------------------------------------------------------------------------------
void CPulseFitter::doDoubleFit(CCyclicCollection<CCell>& cell, float* startPos, float* startAmpl )
{
float oldPos[3],oldAmpl[3];
finalPos[1] = 0.0; finalPos[2] = 0.0;
finalAmpl[1] = 0.0; finalAmpl[2] = 0.0;
numberOfIteration = 0;
for(int i=tMin;i<=tMax;i++)
if(cell(ianode,i)->getAmp()>0.0) w[i]=1.0/cell(ianode,i)->getAmp(); else w[i]=0.0;
for(int ndx=1;ndx<=2;ndx++){
oldPos[ndx] =startPos[ndx];
oldAmpl[ndx]=startAmpl[ndx];
}
for(int iter=1; iter<=maxIterations; iter++){
getNewAmpl(cell, oldPos, finalSigma);
getNewPos( cell, oldPos, newAmpl, finalSigma);
double delta1 = abs(oldPos[1]-newPos[1])+abs(oldPos[2]-newPos[2]);
double delta2 = abs(startPos[1]-newPos[1])+abs(startPos[2]-newPos[2]);
if( delta1 > epsilon && delta2 < deltaMax ){
for(ndx=1;ndx<=2;ndx++){
oldPos[ndx]=newPos[ndx];
oldAmpl[ndx]=newAmpl[ndx];
}
} else break;
}
for(ndx=1;ndx<=2;ndx++){
if(newPos[ndx]>0 && newPos[ndx]<256)
finalPos[ndx] =newPos[ndx]; else finalPos[ndx] =startPos[ndx];
if(newAmpl[ndx]>0 && newAmpl[ndx]<500)
finalAmpl[ndx]=newAmpl[ndx]; else finalAmpl[ndx]=startAmpl[ndx];
}
numberOfIteration=iter;
}
void CPulseFitter::getNewAmpl(CCyclicCollection<CCell>& cell, float* pos, float* sigma)
{
double e1kv=0, e2kv=0, e12=0, ye1=0, ye2=0, arg=0;
float anode;
for(int i=tMin;i<=tMax;i++){
anode=cell(ianode,i)->getAmp();
arg =pow((i-pos[1])/sigma[1],2);
e1kv +=w[i]*exp(-arg);
arg /=2.0;
ye1 += w[i]*anode*exp(-arg);
arg =pow((i-pos[2])/sigma[2],2);
e2kv += w[i]*exp(-arg);
arg /=2.0;
ye2 += w[i]*anode*exp(-arg);
arg =(pow((i-pos[1])/sigma[1],2)+pow((i-pos[2])/sigma[2],2))*0.5;
e12 += w[i]*exp(-arg);
}
double deter = e1kv*e2kv-e12*e12;
if(deter != 0.0){
newAmpl[1]=(ye1*e2kv-e12*ye2)/deter;
newAmpl[2]=(ye2*e1kv-e12*ye1)/deter;
} else {
newAmpl[1]=0.0;
newAmpl[2]=0.0;
}
};
void CPulseFitter::getNewPos(CCyclicCollection<CCell>& cell, float* pos, float* ampl, float* sigma)
{
double d0=0,arg1=0,arg2=0;
double bc=0,bd=0,be=0,bf=0,bg=0,bh=0;
for(int i=tMin; i<=tMax; i++){
arg1 =pow((i-pos[1])/sigma[1],2);
arg2 =pow((i-pos[2])/sigma[2],2);
d0 =cell(ianode,i)->getAmp()-ampl[1]*exp(-0.5*arg1)-ampl[2]*exp(-0.5*arg2);
bc +=(d0*(1+arg1)+ampl[1]*arg1*exp(-0.5*arg1))*exp(-0.5*arg1);
bd +=(i-pos[1])*(i-pos[2])*exp(-0.5*arg1)*exp(-0.5*arg2);
bf +=(d0*(1+arg2)+ampl[2]*arg2*exp(-0.5*arg2))*exp(-0.5*arg2);
bg +=d0*(i-pos[1])*exp(-0.5*arg1);
bh +=d0*(i-pos[2])*exp(-0.5*arg2);
}
be = ampl[1]/pow(sigma[1],2)*bd;
bd *= ampl[2]/pow(sigma[2],2);
double deter = bc*bf-bd*be;
if(deter != 0.0){
newPos[1]=(bg*bf-bd*bh)/deter+pos[1];
newPos[2]=(bc*bh-bg*be)/deter+pos[2];
} else {
newPos[1]=0.0;
newPos[2]=0.0;
}
};
float CPulseFitter::calibratedPulseSigma(float time )
{
if(time<0 || time>255) {
cerr << "time bin out of scope" <<endl;
cerr << "In CSidcRobustFit::calibratedPulseSigma" <<endl;
return 0;
}
return calibratedSigma0+calibratedSigma1*time;
}
CBoolean CPulseFitter::singlePulseMaxSigma(float time, float sigma)
{
// cut on three sigma !!
return sigma < (calibratedPulseSigma(time) + 3*calibratedSigmaDeviation);
}