Update master to aliroot
[u/mrichter/AliRoot.git] / STEER / STEER / AliSignalProcesor.cxx
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
21 #include "AliSignalProcesor.h"
22 #include "AliMathBase.h"
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.5*par3save*sigma2)>100) return 0;   // avoid overflow
49   if (-par5save*(dx-0.5*par5save*sigma2)>100) return 0;   // avoid overflow 
50   if (TMath::Abs(par[4])>1) return 0;
51   Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
52     *(1-AliMathBase::ErfFast((par3save*sigma2-dx)/(sqrt2*par2save)));
53
54   Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
55     *(1-AliMathBase::ErfFast((par5save*sigma2-dx)/(sqrt2*par2save)));
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.5*par3save*sigma2)>100) return 0;   // avoid overflow
82   if (-par5save*(dx-0.5*par5save*sigma2)>100) return 0;   // avoid overflow 
83   if (TMath::Abs(par[4])>=1) return 0;
84   
85   Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
86     *0.5*(1-AliMathBase::ErfFast((par3save*sigma2-dx)/(sqrt2*par2save)));
87
88   Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
89     *0.5*(1-AliMathBase::ErfFast((par5save*sigma2-dx)/(sqrt2*par2save)));
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
104 void AliSignalProcesor::SplineSmoother(const Double_t *ampin, Double_t *ampout, Int_t n) const
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[2*i+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
140 void AliSignalProcesor::TailCancelationALTRO(const Double_t *ampin, Double_t *ampout, Float_t k, Float_t l, Int_t n) const
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++){
148     ampout[i]   = ampin[i]   + (k-l)*temp;
149     temp        = ampin[i]   +  k*temp;
150   }
151 }
152
153 //
154 //
155 void AliSignalProcesor::TailCancelationTRD(const Double_t *ampin, Double_t *ampout, Float_t r, Float_t c, 
156                       Int_t n) const
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+c*ampout[i]);
166   }
167   
168 }
169
170 void AliSignalProcesor::TailMaker(const Double_t *ampin, Double_t *ampout, Float_t lambda, 
171                       Int_t n) const
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.-norm*lambda);
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
208 void AliSignalProcesor::TailCancelationMI(const Double_t *ampin, Double_t *ampout, Float_t norm, 
209                                         Float_t lambda, Int_t n) const
210 {
211   
212   Double_t l = TMath::Exp(-lambda*0.5);
213   Double_t k = l*(1.-norm*lambda*0.5);
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[2*i+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++){
236     out[i]      = in[i]   + (k-l)*temp;
237     temp        = in[i]   +  k*temp;
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
250 void AliSignalProcesor::TailMakerSpline(const Double_t *ampin, Double_t *ampout, Float_t lambda, 
251                       Int_t n) const
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[2*i+1]  = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16;    
277   }
278   //
279   //
280   Float_t temp;
281   out[2*n]     = in[2*n];
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 }