- Three classes by MinJung Kweon AliHFEpriVtx, AliHFEsecVtx and AliHFEmcQA for primar...
[u/mrichter/AliRoot.git] / 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
23
24 ClassImp(AliSignalProcesor)
25
26
27 Double_t asymgauss(Double_t* x, Double_t* par)
28 {
29   // par[0] = normalization
30   // par[1] = mean
31   // par[2] = sigma
32   // norm0  = 1
33   // par[3] = lambda0
34   // par[4] = norm1
35   // par[5] = lambda1
36   //
37   
38   Double_t par1save = par[1];    
39   Double_t par2save = par[2];
40   Double_t par3save = par[3];
41   Double_t par5save = par[5];
42   Double_t dx   = x[0]-par1save;
43   //
44   //
45   Double_t  sigma2  = par2save*par2save;
46   Double_t  sqrt2   = TMath::Sqrt(2.);
47   if (-par3save*(dx-0.5*par3save*sigma2)>100) return 0;   // avoid overflow
48   if (-par5save*(dx-0.5*par5save*sigma2)>100) return 0;   // avoid overflow 
49   if (TMath::Abs(par[4])>1) return 0;
50   Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
51     *(1-TMath::Erf((par3save*sigma2-dx)/(sqrt2*par2save)));
52
53   Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
54     *(1-TMath::Erf((par5save*sigma2-dx)/(sqrt2*par2save)));
55
56
57   return par[0]*(exp1+par[4]*exp2);
58 }
59
60 Double_t asymgaussN(Double_t* x, Double_t* par)
61 {
62   // par[0] = normalization
63   // par[1] = mean
64   // par[2] = sigma
65   // norm0  = 1
66   // par[3] = lambda0
67   // par[4] = norm1
68   // par[5] = lambda1
69   //
70   
71   Double_t par1save = par[1];    
72   Double_t par2save = par[2];
73   Double_t par3save = par[3];
74   Double_t par5save = par[5];
75   Double_t dx   = x[0]-par1save;
76   //
77   //
78   Double_t  sigma2  = par2save*par2save;
79   Double_t  sqrt2   = TMath::Sqrt(2.);
80   if (-par3save*(dx-0.5*par3save*sigma2)>100) return 0;   // avoid overflow
81   if (-par5save*(dx-0.5*par5save*sigma2)>100) return 0;   // avoid overflow 
82   if (TMath::Abs(par[4])>=1) return 0;
83   
84   Double_t  exp1    = par3save*TMath::Exp(-par3save*(dx-0.5*par3save*sigma2))
85     *0.5*(1-TMath::Erf((par3save*sigma2-dx)/(sqrt2*par2save)));
86
87   Double_t  exp2    = par5save*TMath::Exp(-par5save*(dx-0.5*par5save*sigma2))
88     *0.5*(1-TMath::Erf((par5save*sigma2-dx)/(sqrt2*par2save)));
89
90
91   return par[0]*(1.*exp1+par[4]*exp2)/(1.+par[4]);
92 }
93
94
95 TF1 * AliSignalProcesor::GetAsymGauss()
96 {
97   TF1 * f1 = new TF1("asymg",asymgaussN,-10,40,6);
98   return f1;
99 }
100
101
102
103 void AliSignalProcesor::SplineSmoother(Double_t *ampin, Double_t *ampout, Int_t n)
104 {  
105   //
106   //
107   Float_t in[10000];
108   Float_t out[10000];
109   in[0] = ampin[0];
110   in[1] = (ampin[0]+ampin[1])*0.5;
111   in[2*(n-1)]    = ampin[n-1];
112   in[2*(n-1)+1]  = ampin[n-1];
113   //
114   // add charge to the end
115   for (Int_t i=0;i<10;i++){
116     in[2*(n-1)+i]=ampin[n-1];
117   }
118
119   //
120   for (Int_t i=1;i<n-1;i++){
121     in[2*i]    = ampin[i];
122     in[2*i+1]  = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16.;    
123   }
124   //
125   out[0] = in[0];
126   for (Int_t i=1;i<=2*n;i++){
127     out[i]  = (9.*(in[i]+in[i+1])-(in[i-1]+in[i+2]))/16.;    
128   }
129   //
130   //
131   for (int i=0;i<n;i++){
132     ampout[i]      = out[2*i+1]; 
133   }
134 }
135
136
137
138
139 void AliSignalProcesor::TailCancelationALTRO(Double_t *ampin, Double_t *ampout, Float_t K, Float_t L, 
140                       Int_t n)
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(Double_t *ampin, Double_t *ampout, Float_t r, Float_t c, 
156                       Int_t n)
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(Double_t *ampin, Double_t *ampout, Float_t lambda, 
171                       Int_t n)
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(Double_t *ampin, Double_t *ampout, Float_t norm, 
209                                         Float_t lambda, Int_t n)
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(Double_t *ampin, Double_t *ampout, Float_t lambda, 
251                       Int_t n)
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 }