[u/mrichter/AliRoot.git] / STEER / AliSignalProcesor.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

#include <TF1.h>
#include <TMath.h>

#include "AliSignalProcesor.h"
#include "AliMathBase.h"


ClassImp(AliSignalProcesor)


Double_t asymgauss(Double_t* x, Double_t* par)
{
  // par[0] = normalization
  // par[1] = mean
  // par[2] = sigma
  // norm0  = 1
  // par[3] = lambda0
  // par[4] = norm1
  // par[5] = lambda1
  //
  
  Double_t par1save = par[1];    
  Double_t par2save = par[2];
  Double_t par3save = par[3];
  Double_t par5save = par[5];
  Double_t dx   = x[0]-par1save;
  //
  //
  Double_t  sigma2  = par2save*par2save;
  Double_t  sqrt2   = TMath::Sqrt(2.);
  if (-par3save*(dx-0.5*par3save*sigma2)>100) return 0;   // avoid overflow
  if (-par5save*(dx-0.5*par5save*sigma2)>100) return 0;   // avoid overflow 
  if (TMath::Abs(par[4])>1) return 0;
  Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
    *(1-AliMathBase::ErfFast((par3save*sigma2-dx)/(sqrt2*par2save)));

  Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
    *(1-AliMathBase::ErfFast((par5save*sigma2-dx)/(sqrt2*par2save)));


  return par[0]*(exp1+par[4]*exp2);
}

Double_t asymgaussN(Double_t* x, Double_t* par)
{
  // par[0] = normalization
  // par[1] = mean
  // par[2] = sigma
  // norm0  = 1
  // par[3] = lambda0
  // par[4] = norm1
  // par[5] = lambda1
  //
  
  Double_t par1save = par[1];    
  Double_t par2save = par[2];
  Double_t par3save = par[3];
  Double_t par5save = par[5];
  Double_t dx   = x[0]-par1save;
  //
  //
  Double_t  sigma2  = par2save*par2save;
  Double_t  sqrt2   = TMath::Sqrt(2.);
  if (-par3save*(dx-0.5*par3save*sigma2)>100) return 0;   // avoid overflow
  if (-par5save*(dx-0.5*par5save*sigma2)>100) return 0;   // avoid overflow 
  if (TMath::Abs(par[4])>=1) return 0;
  
  Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
    *0.5*(1-AliMathBase::ErfFast((par3save*sigma2-dx)/(sqrt2*par2save)));

  Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
    *0.5*(1-AliMathBase::ErfFast((par5save*sigma2-dx)/(sqrt2*par2save)));


  return par[0]*(1.*exp1+par[4]*exp2)/(1.+par[4]);
}


TF1 * AliSignalProcesor::GetAsymGauss()
{
  TF1 * f1 = new TF1("asymg",asymgaussN,-10,40,6);
  return f1;
}


void AliSignalProcesor::SplineSmoother(const Double_t *ampin, Double_t *ampout, Int_t n) const
{  
  //
  //
  Float_t in[10000];
  Float_t out[10000];
  in[0] = ampin[0];
  in[1] = (ampin[0]+ampin[1])*0.5;
  in[2*(n-1)]    = ampin[n-1];
  in[2*(n-1)+1]  = ampin[n-1];
  //
  // add charge to the end
  for (Int_t i=0;i<10;i++){
    in[2*(n-1)+i]=ampin[n-1];
  }

  //
  for (Int_t i=1;i<n-1;i++){
    in[2*i]    = ampin[i];
    in[2*i+1]  = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16.;    
  }
  //
  out[0] = in[0];
  for (Int_t i=1;i<=2*n;i++){
    out[i]  = (9.*(in[i]+in[i+1])-(in[i-1]+in[i+2]))/16.;    
  }
  //
  //
  for (int i=0;i<n;i++){
    ampout[i]      = out[2*i+1]; 
  }
}


void AliSignalProcesor::TailCancelationALTRO(const Double_t *ampin, Double_t *ampout, Float_t k, Float_t l, Int_t n) const
{
  //
  // ALTRO
  Float_t temp;
  ampout[0]  = ampin[0];
  temp = ampin[0];
  for (int i=1;i<n;i++){
    ampout[i]   = ampin[i]   + (k-l)*temp;
    temp        = ampin[i]   +  k*temp;
  }
}

//
//
void AliSignalProcesor::TailCancelationTRD(const Double_t *ampin, Double_t *ampout, Float_t r, Float_t c, 
		      Int_t n) const
{
  //TRD
  //
  Double_t reminder=0;
  //
  for (Int_t i=0; i<n; i++){
    ampout[i] = ampin[i]-reminder;
    //
    reminder  = r*(reminder+c*ampout[i]);
  }
  
}

void AliSignalProcesor::TailMaker(const Double_t *ampin, Double_t *ampout, Float_t lambda, 
		      Int_t n) const
{
  
  Double_t l = TMath::Exp(-lambda);
  //
  Float_t temp=0;
  for (Int_t i=n-1; i>0; i--){
    ampout[i] = ampin[i]+temp;
    //
    temp  = l*(temp+ampin[i]);
  }
}

void AliSignalProcesor::TailCancelationALTRO1(Double_t *ampin, Double_t *ampout, Float_t norm, 
					    Float_t lambda, Int_t n)
{
  
  Double_t l = TMath::Exp(-lambda);
  Double_t k = l*(1.-norm*lambda);

  return TailCancelationALTRO(ampin,ampout,k,l,n);
}


void AliSignalProcesor::TailCancelationTRD1(Double_t *ampin, Double_t *ampout, Float_t norm, 
					  Float_t lambda, Int_t n)
{
  //
  //
  Double_t r = TMath::Exp(-lambda);
  Double_t c = norm*lambda;
  return TailCancelationTRD(ampin,ampout,r,c,n);
}


void AliSignalProcesor::TailCancelationMI(const Double_t *ampin, Double_t *ampout, Float_t norm, 
					Float_t lambda, Int_t n) const
{
  
  Double_t l = TMath::Exp(-lambda*0.5);
  Double_t k = l*(1.-norm*lambda*0.5);
  //
  //
  Float_t in[10000];
  Float_t out[10000];
  for (Int_t i=0;i<n*2+20;i++){
    in[i] = 0;
    out[i]= 0;
  }
  in[0] = ampin[0];
  in[1] = (ampin[0]+ampin[1])*0.5;
  in[2*(n-1)]    = ampin[n-1];
  in[2*(n-1)+1]  = ampin[n-1];
  //
  for (Int_t i=1;i<n-2;i++){
    in[2*i]    = ampin[i];
    in[2*i+1]  = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16;    
  }
  //
  Float_t temp;
  out[0]     = in[0];
  temp       = in[0];
  for (int i=1;i<=2*n;i++){
    out[i]      = in[i]   + (k-l)*temp;
    temp        = in[i]   +  k*temp;
  }
  //
  //
  for (int i=0;i<n;i++){
    ampout[i]      = out[2*i+1];
  }
}


void AliSignalProcesor::TailMakerSpline(const Double_t *ampin, Double_t *ampout, Float_t lambda, 
		      Int_t n) const
{
  
  Double_t l = TMath::Exp(-lambda*0.5);
  //
  //
  Float_t in[10000];
  Float_t out[10000];
  for (Int_t i=0;i<n*2+20;i++){
    in[i] = 0;
    out[i]= 0;
  }
  in[0] = ampin[0];
  in[1] = (ampin[0]+ampin[1])*0.5;
  in[2*(n-1)]    = ampin[n-1];
  in[2*(n-1)+1]  = ampin[n-1];
  //
  // add charge to the end
  for (Int_t i=0;i<10;i++){
    in[2*(n-1)+i]=ampin[n-1];
  }

  //
  for (Int_t i=1;i<n-2;i++){
    in[2*i]    = ampin[i];
    in[2*i+1]  = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16;    
  }
  //
  //
  Float_t temp;
  out[2*n]     = in[2*n];
  temp         = 0;
  for (int i=2*n;i>=0;i--){
    out[i]      = in[i]   + temp;
    temp        = l*(temp+in[i]);
  }
  //
  //
  for (int i=0;i<n;i++){
    ampout[i]      = out[2*i+1]; 
  }
}
Commit	Line	Data
090026bf	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	#include <TF1.h>
	19	#include <TMath.h>
	20
ec143d80	21	#include "AliSignalProcesor.h"
bb7e41dd	22	#include "AliMathBase.h"
ec143d80	23
	24
	25	ClassImp(AliSignalProcesor)
	26
	27
	28	Double_t asymgauss(Double_t* x, Double_t* par)
	29	{
	30	// par[0] = normalization
	31	// par[1] = mean
	32	// par[2] = sigma
	33	// norm0 = 1
	34	// par[3] = lambda0
	35	// par[4] = norm1
	36	// par[5] = lambda1
	37	//
	38
	39	Double_t par1save = par[1];
	40	Double_t par2save = par[2];
	41	Double_t par3save = par[3];
	42	Double_t par5save = par[5];
	43	Double_t dx = x[0]-par1save;
	44	//
	45	//
	46	Double_t sigma2 = par2save*par2save;
	47	Double_t sqrt2 = TMath::Sqrt(2.);
	48	if (-par3save(dx-0.5par3save*sigma2)>100) return 0; // avoid overflow
	49	if (-par5save(dx-0.5par5save*sigma2)>100) return 0; // avoid overflow
	50	if (TMath::Abs(par[4])>1) return 0;
	51	Double_t exp1 = par3saveTMath::Exp(-par3save(dx-0.5par3savesigma2))
bb7e41dd	52	(1-AliMathBase::ErfFast((par3savesigma2-dx)/(sqrt2*par2save)));
ec143d80	53
ec143d80	54	Double_t exp2 = par5saveTMath::Exp(-par5save(dx-0.5par5savesigma2))
bb7e41dd	55	(1-AliMathBase::ErfFast((par5savesigma2-dx)/(sqrt2*par2save)));
ec143d80	56
	57
	58	return par[0](exp1+par[4]exp2);
	59	}
	60
	61	Double_t asymgaussN(Double_t* x, Double_t* par)
	62	{
	63	// par[0] = normalization
	64	// par[1] = mean
	65	// par[2] = sigma
	66	// norm0 = 1
	67	// par[3] = lambda0
	68	// par[4] = norm1
	69	// par[5] = lambda1
	70	//
	71
	72	Double_t par1save = par[1];
	73	Double_t par2save = par[2];
	74	Double_t par3save = par[3];
	75	Double_t par5save = par[5];
	76	Double_t dx = x[0]-par1save;
	77	//
	78	//
	79	Double_t sigma2 = par2save*par2save;
	80	Double_t sqrt2 = TMath::Sqrt(2.);
	81	if (-par3save(dx-0.5par3save*sigma2)>100) return 0; // avoid overflow
	82	if (-par5save(dx-0.5par5save*sigma2)>100) return 0; // avoid overflow
	83	if (TMath::Abs(par[4])>=1) return 0;
	84
	85	Double_t exp1 = par3saveTMath::Exp(-par3save(dx-0.5par3savesigma2))
bb7e41dd	86	0.5(1-AliMathBase::ErfFast((par3savesigma2-dx)/(sqrt2par2save)));
ec143d80	87
ec143d80	88	Double_t exp2 = par5saveTMath::Exp(-par5save(dx-0.5par5savesigma2))
bb7e41dd	89	0.5(1-AliMathBase::ErfFast((par5savesigma2-dx)/(sqrt2par2save)));
ec143d80	90
	91
	92	return par[0](1.exp1+par[4]*exp2)/(1.+par[4]);
	93	}
	94
	95
	96	TF1 * AliSignalProcesor::GetAsymGauss()
	97	{
	98	TF1 * f1 = new TF1("asymg",asymgaussN,-10,40,6);
	99	return f1;
	100	}
	101
	102
	103
0ee094cc	104	void AliSignalProcesor::SplineSmoother(const Double_t ampin, Double_t ampout, Int_t n) const
ec143d80	105	{
	106	//
	107	//
	108	Float_t in[10000];
	109	Float_t out[10000];
	110	in[0] = ampin[0];
	111	in[1] = (ampin[0]+ampin[1])*0.5;
	112	in[2*(n-1)] = ampin[n-1];
	113	in[2*(n-1)+1] = ampin[n-1];
	114	//
	115	// add charge to the end
	116	for (Int_t i=0;i<10;i++){
	117	in[2*(n-1)+i]=ampin[n-1];
	118	}
	119
	120	//
	121	for (Int_t i=1;i<n-1;i++){
	122	in[2*i] = ampin[i];
	123	in[2i+1] = (9.(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16.;
	124	}
	125	//
	126	out[0] = in[0];
	127	for (Int_t i=1;i<=2*n;i++){
	128	out[i] = (9.*(in[i]+in[i+1])-(in[i-1]+in[i+2]))/16.;
	129	}
	130	//
	131	//
	132	for (int i=0;i<n;i++){
	133	ampout[i] = out[2*i+1];
	134	}
	135	}
	136
	137
	138
	139
0ee094cc	140	void AliSignalProcesor::TailCancelationALTRO(const Double_t ampin, Double_t ampout, Float_t k, Float_t l, Int_t n) const
ec143d80	141	{
	142	//
	143	// ALTRO
	144	Float_t temp;
	145	ampout[0] = ampin[0];
	146	temp = ampin[0];
	147	for (int i=1;i<n;i++){
0ee094cc	148	ampout[i] = ampin[i] + (k-l)*temp;
0ee094cc	149	temp = ampin[i] + k*temp;
ec143d80	150	}
	151	}
	152
	153	//
	154	//
0ee094cc	155	void AliSignalProcesor::TailCancelationTRD(const Double_t ampin, Double_t ampout, Float_t r, Float_t c,
0ee094cc	156	Int_t n) const
ec143d80	157	{
	158	//TRD
	159	//
	160	Double_t reminder=0;
	161	//
	162	for (Int_t i=0; i<n; i++){
	163	ampout[i] = ampin[i]-reminder;
	164	//
	165	reminder = r(reminder+campout[i]);
	166	}
	167
	168	}
	169
0ee094cc	170	void AliSignalProcesor::TailMaker(const Double_t ampin, Double_t ampout, Float_t lambda,
0ee094cc	171	Int_t n) const
ec143d80	172	{
	173
	174	Double_t l = TMath::Exp(-lambda);
	175	//
	176	Float_t temp=0;
	177	for (Int_t i=n-1; i>0; i--){
	178	ampout[i] = ampin[i]+temp;
	179	//
	180	temp = l*(temp+ampin[i]);
	181	}
	182	}
	183
	184	void AliSignalProcesor::TailCancelationALTRO1(Double_t ampin, Double_t ampout, Float_t norm,
	185	Float_t lambda, Int_t n)
	186	{
	187
	188	Double_t l = TMath::Exp(-lambda);
	189	Double_t k = l(1.-normlambda);
	190
	191	return TailCancelationALTRO(ampin,ampout,k,l,n);
	192	}
	193
	194
	195	void AliSignalProcesor::TailCancelationTRD1(Double_t ampin, Double_t ampout, Float_t norm,
	196	Float_t lambda, Int_t n)
	197	{
	198	//
	199	//
	200	Double_t r = TMath::Exp(-lambda);
	201	Double_t c = norm*lambda;
	202	return TailCancelationTRD(ampin,ampout,r,c,n);
	203	}
	204
	205
	206
	207
0ee094cc	208	void AliSignalProcesor::TailCancelationMI(const Double_t ampin, Double_t ampout, Float_t norm,
0ee094cc	209	Float_t lambda, Int_t n) const
ec143d80	210	{
ec143d80	211
0ee094cc	212	Double_t l = TMath::Exp(-lambda*0.5);
0ee094cc	213	Double_t k = l(1.-normlambda*0.5);
ec143d80	214	//
	215	//
	216	Float_t in[10000];
	217	Float_t out[10000];
	218	for (Int_t i=0;i<n*2+20;i++){
	219	in[i] = 0;
	220	out[i]= 0;
	221	}
	222	in[0] = ampin[0];
	223	in[1] = (ampin[0]+ampin[1])*0.5;
	224	in[2*(n-1)] = ampin[n-1];
	225	in[2*(n-1)+1] = ampin[n-1];
	226	//
	227	for (Int_t i=1;i<n-2;i++){
	228	in[2*i] = ampin[i];
	229	in[2i+1] = (9.(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16;
	230	}
	231	//
	232	Float_t temp;
	233	out[0] = in[0];
	234	temp = in[0];
	235	for (int i=1;i<=2*n;i++){
0ee094cc	236	out[i] = in[i] + (k-l)*temp;
0ee094cc	237	temp = in[i] + k*temp;
ec143d80	238	}
	239	//
	240	//
	241	for (int i=0;i<n;i++){
	242	ampout[i] = out[2*i+1];
	243	}
	244	}
	245
	246
	247
	248
	249
0ee094cc	250	void AliSignalProcesor::TailMakerSpline(const Double_t ampin, Double_t ampout, Float_t lambda,
0ee094cc	251	Int_t n) const
ec143d80	252	{
	253
	254	Double_t l = TMath::Exp(-lambda*0.5);
	255	//
	256	//
	257	Float_t in[10000];
	258	Float_t out[10000];
	259	for (Int_t i=0;i<n*2+20;i++){
	260	in[i] = 0;
	261	out[i]= 0;
	262	}
	263	in[0] = ampin[0];
	264	in[1] = (ampin[0]+ampin[1])*0.5;
	265	in[2*(n-1)] = ampin[n-1];
	266	in[2*(n-1)+1] = ampin[n-1];
	267	//
	268	// add charge to the end
	269	for (Int_t i=0;i<10;i++){
	270	in[2*(n-1)+i]=ampin[n-1];
	271	}
	272
	273	//
	274	for (Int_t i=1;i<n-2;i++){
	275	in[2*i] = ampin[i];
	276	in[2i+1] = (9.(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16;
	277	}
	278	//
	279	//
	280	Float_t temp;
	281	out[2n] = in[2n];
	282	temp = 0;
	283	for (int i=2*n;i>=0;i--){
	284	out[i] = in[i] + temp;
	285	temp = l*(temp+in[i]);
	286	}
	287	//
	288	//
	289	for (int i=0;i<n;i++){
	290	ampout[i] = out[2*i+1];
	291	}
	292	}