running local for SP and LYZEP and new histograms for QC
[u/mrichter/AliRoot.git] / TRD / AliTRDmcm.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 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 //                                                                        //
21 //  Multi Chip Module exponential filter and tracklet finder              //
22 //                                                                        //
23 //                                                                        //
24 ////////////////////////////////////////////////////////////////////////////
25
26 #include <TMath.h>
27
28 #include "AliLog.h"
29
30 #include "AliTRDmcm.h"
31 #include "AliTRDtrigParam.h"
32
33 ClassImp(AliTRDmcm)
34
35 //_____________________________________________________________________________
36 AliTRDmcm::AliTRDmcm() 
37   :TObject()
38   ,fNtrk(0)
39   ,fRobId(0)
40   ,fChaId(0)
41   ,fRow(0)
42   ,fColFirst(0)
43   ,fColLast(0)
44   ,fTime1(0)
45   ,fTime2(0)
46   ,fClusThr(0)
47   ,fPadThr(0)
48   ,fNtrkSeeds(0)
49   ,fR1(0)
50   ,fR2(0)
51   ,fC1(0)
52   ,fC2(0)
53   ,fPedestal(0)
54   ,fId(0)
55 {
56   //
57   // AliTRDmcm default constructor
58   //
59
60   Int_t i = 0;
61   Int_t j = 0;
62
63   for (i = 0; i < kMaxTrackletsPerMCM; i++) {
64     fTrkIndex[i] = 0;
65     fSeedCol[i]  = -1;
66   }
67   for (i = 0; i < kMcmCol; i++) {
68     fPadHits[i]  = 0;
69     for (j = 0; j < kMcmTBmax; j++) {
70       fADC[i][j]    = 0.0;
71       fIsClus[i][j] = kFALSE;
72     }
73   }
74
75 }
76
77 //_____________________________________________________________________________
78 AliTRDmcm::AliTRDmcm(Int_t id) 
79   :TObject()
80   ,fNtrk(0)
81   ,fRobId(0)
82   ,fChaId(0)
83   ,fRow(0)
84   ,fColFirst(0)
85   ,fColLast(0)
86   ,fTime1(0)
87   ,fTime2(0)
88   ,fClusThr(0)
89   ,fPadThr(0)
90   ,fNtrkSeeds(0)
91   ,fR1(0)
92   ,fR2(0)
93   ,fC1(0)
94   ,fC2(0)
95   ,fPedestal(0)
96   ,fId(id)
97 {
98   //
99   // AliTRDmcm constructor
100   //
101
102   Int_t i = 0;
103   Int_t j = 0;
104
105   for (i = 0; i < kMaxTrackletsPerMCM; i++) {
106     fTrkIndex[i] = 0;
107     fSeedCol[i]  = -1;
108   }
109   for (i = 0; i < kMcmCol; i++) {
110     fPadHits[i]  = 0;
111     for (j = 0; j < kMcmTBmax; j++) {
112       fADC[i][j]    = 0.0;
113       fIsClus[i][j] = kFALSE;
114     }
115   }
116
117   fTime1   = AliTRDtrigParam::Instance()->GetTime1();
118   fTime2   = AliTRDtrigParam::Instance()->GetTime2();
119   fClusThr = AliTRDtrigParam::Instance()->GetClusThr();
120   fPadThr  = AliTRDtrigParam::Instance()->GetPadThr();
121   
122   AliTRDtrigParam::Instance()->GetFilterParam(fR1,fR2,fC1,fC2,fPedestal);
123
124 }
125
126 //_____________________________________________________________________________
127 AliTRDmcm::AliTRDmcm(const AliTRDmcm &m) 
128   :TObject(m)
129   ,fNtrk(m.fNtrk)
130   ,fRobId(m.fRobId)
131   ,fChaId(m.fChaId)
132   ,fRow(m.fRow)
133   ,fColFirst(m.fColFirst)
134   ,fColLast(m.fColLast)
135   ,fTime1(m.fTime1)
136   ,fTime2(m.fTime2)
137   ,fClusThr(m.fClusThr)
138   ,fPadThr(m.fPadThr)
139   ,fNtrkSeeds(m.fNtrkSeeds)
140   ,fR1(m.fR1)
141   ,fR2(m.fR2)
142   ,fC1(m.fC1)
143   ,fC2(m.fC2)
144   ,fPedestal(m.fPedestal)
145   ,fId(m.fId)
146 {
147   //
148   // AliTRDmcm copy constructor
149   //
150
151   Int_t i = 0;
152   Int_t j = 0;
153
154   for (i = 0; i < kMaxTrackletsPerMCM; i++) {
155     ((AliTRDmcm &) m).fTrkIndex[i] = 0;
156     ((AliTRDmcm &) m).fSeedCol[i]  = -1;
157   }
158   for (i = 0; i < kMcmCol; i++) {
159     ((AliTRDmcm &) m).fPadHits[i]  = 0;
160     for (j = 0; j < kMcmTBmax; j++) {
161       ((AliTRDmcm &) m).fADC[i][j]    = 0.0;
162       ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
163     }
164   }
165
166 }
167
168 //_____________________________________________________________________________
169 AliTRDmcm::~AliTRDmcm() 
170 {
171   //
172   // AliTRDmcm destructor
173   //
174
175 }
176
177 //_____________________________________________________________________________
178 AliTRDmcm &AliTRDmcm::operator=(const AliTRDmcm &m)
179 {
180   //
181   // Assignment operator
182   //
183
184   if (this != &m) ((AliTRDmcm &) m).Copy(*this); 
185
186   return *this;
187
188 }
189
190 //_____________________________________________________________________________
191 void AliTRDmcm::Copy(TObject &m) const
192 {
193   //
194   // Copy function
195   //
196
197   Int_t i = 0;
198   Int_t j = 0;
199
200   ((AliTRDmcm &) m).fNtrk       = fNtrk;
201   ((AliTRDmcm &) m).fRobId      = fRobId;
202   ((AliTRDmcm &) m).fChaId      = fChaId;
203   ((AliTRDmcm &) m).fRow        = fRow;
204   ((AliTRDmcm &) m).fColFirst   = fColFirst;
205   ((AliTRDmcm &) m).fColLast    = fColLast;
206   ((AliTRDmcm &) m).fTime1      = fTime1;
207   ((AliTRDmcm &) m).fTime2      = fTime2;
208   ((AliTRDmcm &) m).fClusThr    = fClusThr;
209   ((AliTRDmcm &) m).fPadThr     = fPadThr;
210   ((AliTRDmcm &) m).fNtrkSeeds  = fNtrkSeeds;
211   ((AliTRDmcm &) m).fR1         = fR1;
212   ((AliTRDmcm &) m).fR2         = fR2;
213   ((AliTRDmcm &) m).fC1         = fC1;
214   ((AliTRDmcm &) m).fC2         = fC2;
215   ((AliTRDmcm &) m).fPedestal   = fPedestal;
216   ((AliTRDmcm &) m).fId         = fId;
217
218   for (i = 0; i < kMaxTrackletsPerMCM; i++) {
219     ((AliTRDmcm &) m).fTrkIndex[i] = 0;
220     ((AliTRDmcm &) m).fSeedCol[i]  = -1;
221   }
222   for (i = 0; i < kMcmCol; i++) {
223     ((AliTRDmcm &) m).fPadHits[i]  = 0;
224     for (j = 0; j < kMcmTBmax; j++) {
225       ((AliTRDmcm &) m).fADC[i][j]    = 0.0;
226       ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
227     }
228   }
229
230 }
231
232 //_____________________________________________________________________________
233 void AliTRDmcm::AddTrk(Int_t id) 
234 {
235   //
236   // Add a tracklet index
237   //
238
239   fTrkIndex[fNtrk] = id;
240   fNtrk++;
241
242   return;
243
244 }
245
246 //_____________________________________________________________________________
247 void AliTRDmcm::Reset()
248 {
249   //
250   // Reset MCM data
251   //
252
253   Int_t i = 0;
254   Int_t j = 0;
255
256   for (i = 0; i < kMcmCol; i++) {
257     fPadHits[i] = 0;
258     for (j = 0; j < kMcmTBmax; j++) {
259       fADC[i][j]    = 0.0;
260       fIsClus[i][j] = kFALSE;
261     }
262   }
263   for (i = 0; i < kMaxTrackletsPerMCM; i++) {
264     fSeedCol[i] = -1;
265   }
266   
267 }
268
269 //_____________________________________________________________________________
270 Bool_t AliTRDmcm::Run()
271 {
272   //
273   // Run MCM
274   //
275
276   AliDebug(2,(Form("Run MCM %d\n",Id())));
277
278   Int_t   iTime = 0;
279   Int_t   iCol  = 0;
280   Int_t   iPlus = 0;
281   Int_t   i     = 0;
282   Int_t   j     = 0;
283
284   Float_t amp[3] = { 0.0, 0.0, 0.0 };
285   Int_t   nClus;
286   Int_t   clusCol[kMcmCol/2];
287   Float_t clusAmp[kMcmCol/2];
288   Float_t veryLarge;
289   Int_t   clusMin = -1;
290   
291   // Main TB loop
292   for (iTime = fTime1; iTime <= fTime2; iTime++) {  
293
294     // Find clusters...
295     nClus = 0;
296     for (iCol = 1; iCol < (kMcmCol-1); iCol++) {
297       amp[0] = fADC[iCol-1][iTime];
298       amp[1] = fADC[iCol  ][iTime];
299       amp[2] = fADC[iCol+1][iTime];
300       if (IsCluster(amp)) {
301         fIsClus[iCol][iTime] = kTRUE;
302         clusCol[nClus]       = iCol;
303         clusAmp[nClus]       = amp[0]+amp[1]+amp[2];
304         nClus++;
305         if (nClus == kMcmCol/2) {
306           AliWarning(Form("Too many clusters in time bin %2d MCM %d...\n",iTime,Id()));
307           break;
308         }
309       }
310     }
311
312     // ...but no more than six...
313     if (nClus > (Int_t) kSelClus) {
314       for (j = kSelClus/2; j < nClus-kSelClus/2; j++) {
315         fIsClus[clusCol[j]][iTime] = kFALSE;
316       } 
317     }
318
319     // ...and take the largest four.
320     Int_t nClusPlus = nClus - kMaxClus;
321     for (iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
322       veryLarge = 1.E+10;
323       for (i = 0; i < nClus; i++) {
324         if (fIsClus[clusCol[i]][iTime]) {
325           if (clusAmp[i] <= veryLarge) {
326             veryLarge = clusAmp[i];
327             clusMin = i;
328           }
329         }
330       }
331       fIsClus[clusCol[clusMin]][iTime] = kFALSE;
332     }
333
334     AddTimeBin(iTime);
335
336   }  // end main TB loop
337     
338   if ((fNtrkSeeds = CreateSeeds())) {
339     return kTRUE;
340   }
341
342   return kFALSE;
343
344 }
345
346 //_____________________________________________________________________________
347 Int_t AliTRDmcm::CreateSeeds()
348 {
349   //
350   // Make column seeds (from Falk Lesser, ex KIP)
351   //
352
353   Int_t i      = 0;
354   Int_t j      = 0;
355   Int_t nSeeds = 0;
356
357   AliDebug(2,Form("AliTRDmcm::CreateSeeds MCM %d \n",Id()));
358
359   // Working array for hit sums
360   Int_t fHit2padSum[2][kMcmCol];            
361
362   // Initialize the array
363   for (i = 0; i < 2; i++) {
364     for (j = 0; j < kMcmCol; j++ ) {
365       if (i == 0) {
366         fHit2padSum[i][j] = j; 
367       } else {
368         fHit2padSum[i][j] = -1; 
369       }
370     }
371   }
372
373   Int_t sum10 = AliTRDtrigParam::Instance()->GetSum10();
374   Int_t sum12 = AliTRDtrigParam::Instance()->GetSum12();
375
376   // Build the 2padSum
377   Int_t nsum2seed = 0;
378   for (i = 0; i < kMcmCol; i++) {
379     if (i < (kMcmCol-1)) {
380       if ((fPadHits[i] >= sum10) && ((fPadHits[i] + fPadHits[i+1]) >= sum12)) {
381         fHit2padSum[1][i] = fPadHits[i] + fPadHits[i+1]; 
382       } 
383       else {
384         fHit2padSum[1][i] = -1;
385       }
386     } 
387     else {
388       if (fPadHits[i] >= sum12) {
389         fHit2padSum[1][i] = fPadHits[i]; 
390       } 
391       else {
392         fHit2padSum[1][i] = -1;
393       }
394     }
395     if (fHit2padSum[1][i] > 0) {
396       nsum2seed++;
397     }
398   }
399
400   // sort the sums in decreasing order of the amplitude 
401   Sort(kMcmCol,&fHit2padSum[0][0],&fHit2padSum[1][0],1);
402
403   // arrange (maximum number of) candidates in increasing order of the column number
404   nSeeds = TMath::Min(nsum2seed,kMaxTrackletsPerMCM);
405   Sort(nSeeds,&fHit2padSum[1][0],&fHit2padSum[0][0],0);
406
407   for (i = 0; i < nSeeds; i++) {
408     fSeedCol[i] = fHit2padSum[0][i];
409   }
410
411   // reject multiple found tracklets
412   Int_t imax = nSeeds - 1;
413   for (i = 0; i < imax; i++) {
414
415     if ((fHit2padSum[0][i]+1) == fHit2padSum[0][i+1]) {
416       nSeeds--;
417       if (fHit2padSum[1][i] >= fHit2padSum[1][i+1]) {
418         AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]));
419         fSeedCol[i+1] = -1;
420       } 
421       else {
422         AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]));
423         fSeedCol[i]   = -1;
424       }
425     }
426
427   }
428
429   return nSeeds;
430
431 }
432
433 //_____________________________________________________________________________
434 void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir) const
435 {
436   //
437   // Sort two parallel vectors (x1[nel], x2[nel]) after the second one (x2)
438   // in the direction: dir = 0 ascending order
439   //                   dir = 1 descending order
440   //
441
442   Int_t  i = 0;
443   Bool_t sort;
444   Int_t  tmp1;
445   Int_t  tmp2;
446
447   if (dir == 0) {
448
449     do { 
450       sort = kTRUE;
451       for (i = 0; i < (nel-1); i++) {
452         if (x2[i+1] < x2[i]) {
453           tmp2    = x2[i]; 
454           x2[i]   = x2[i+1]; 
455           x2[i+1] = tmp2;
456           tmp1    = x1[i]; 
457           x1[i]   = x1[i+1]; 
458           x1[i+1] = tmp1;
459           sort    = kFALSE;
460         }
461       }
462     } while (sort == kFALSE);
463
464   }
465
466   if (dir == 1) {
467
468     do { 
469       sort = kTRUE;
470       for (i = 0; i < (nel-1); i++) {
471         if (x2[i+1] > x2[i]) {
472           tmp2    = x2[i]; 
473           x2[i]   = x2[i+1]; 
474           x2[i+1] = tmp2;
475           tmp1    = x1[i]; 
476           x1[i]   = x1[i+1]; 
477           x1[i+1] = tmp1;
478           sort    = kFALSE;
479         }
480       }
481     } while (sort == kFALSE);
482
483   }
484
485 }
486
487 //_____________________________________________________________________________
488 void AliTRDmcm::AddTimeBin(const Int_t iTime)
489 {
490   //
491   // Build column seeds
492   //
493
494   for (Int_t iPad = 1; iPad < (kMcmCol-1); iPad++) {
495     if (fIsClus[iPad][iTime]) {
496       fPadHits[iPad]++;
497     }
498   }
499
500 }
501
502 //_____________________________________________________________________________
503 Bool_t AliTRDmcm::IsCluster(Float_t amp[3]) const
504 {
505   //
506   // Find if the amplitudes amp[0], amp[1], amp[2] are a cluster
507   //
508
509   // -> shape
510   if (amp[0] > amp[1] || amp[2] > amp[1]) {
511     return kFALSE;
512   }
513
514   // -> cluster amplitude
515   if ((amp[0]+amp[1]+amp[2]) < fClusThr) {
516     return kFALSE;
517   }
518
519   // -> pad amplitude
520   if (amp[0] < fPadThr && amp[2] < fPadThr) {
521     return kFALSE;
522   }
523
524   return kTRUE;
525
526 }
527
528 //_____________________________________________________________________________
529 void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
530 {
531   //
532   // Exponential filter
533   //
534
535   Int_t iCol  = 0;
536   Int_t iTime = 0;
537
538   Double_t sour[kMcmTBmax];
539   Double_t dtarg[kMcmTBmax];
540   Int_t    itarg[kMcmTBmax];
541
542   switch(ftype) {
543
544    case 0:
545
546     for (iCol = 0; iCol < kMcmCol; iCol++) {
547       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
548         sour[iTime] = fADC[iCol][iTime];
549       }
550       DeConvExpA(sour,dtarg,kMcmTBmax,nexp);
551       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
552         fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
553       }
554     }
555     break;
556
557   case 1:
558
559     for (iCol = 0; iCol < kMcmCol; iCol++) {
560       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
561         sour[iTime] = fADC[iCol][iTime];
562       }
563       DeConvExpD(sour,itarg,kMcmTBmax,nexp);
564       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
565         fADC[iCol][iTime] = itarg[iTime];
566       }
567     }
568     break;
569
570   case 2:
571
572     for (iCol = 0; iCol < kMcmCol; iCol++) {
573       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
574         sour[iTime] = fADC[iCol][iTime];
575       }
576       DeConvExpMI(sour,dtarg,kMcmTBmax);
577       for (iTime = 0; iTime < kMcmTBmax; iTime++) {
578         fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
579       }
580     }
581     break;
582
583   default:
584
585     AliError(Form("Invalid filter type %d ! \n",ftype));
586     return;
587
588   }
589
590 }
591
592 //_____________________________________________________________________________
593 void AliTRDmcm::DeConvExpA(Double_t *source, Double_t *target, Int_t n, Int_t nexp) 
594 {
595   //
596   // Exponential filter "analog"
597   //
598
599   Int_t    i = 0;
600   Int_t    k = 0;
601   Double_t reminder[2];
602   Double_t correction;
603   Double_t result;
604   Double_t rates[2];
605   Double_t coefficients[2];
606
607   // Initialize (coefficient = alpha, rates = lambda)
608   
609   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
610   Double_t r1, r2, c1, c2;
611   r1 = (Double_t)fR1;
612   r2 = (Double_t)fR2;
613   c1 = (Double_t)fC1;
614   c2 = (Double_t)fC2;
615   
616   coefficients[0] = c1;
617   coefficients[1] = c2;
618
619   Double_t dt = 0.1;
620   rates[0] = TMath::Exp(-dt/(r1));
621   rates[1] = TMath::Exp(-dt/(r2));
622
623   // Attention: computation order is important
624   correction = 0.0;
625   for (k = 0; k < nexp; k++) {
626     reminder[k] = 0.0;
627   }
628     
629   for (i = 0; i < n; i++) {
630
631     result    = (source[i] - correction);    // no rescaling
632     target[i] = result;
633     
634     for (k = 0; k < nexp; k++) {
635       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
636     }
637       
638     correction = 0.0;
639     for (k = 0; k < nexp; k++) {
640       correction += reminder[k];
641     }
642
643   }
644   
645 }
646
647 //_____________________________________________________________________________
648 void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp) 
649 {
650   //
651   // Exponential filter "digital"
652   //
653
654   Int_t i = 0;
655
656   Int_t fAlphaL;
657   Int_t fAlphaS;
658   Int_t fLambdaL;
659   Int_t fLambdaS;
660   Int_t fTailPed;
661
662   Int_t iAlphaL;
663   Int_t iAlphaS;
664   Int_t iLambdaL;
665   Int_t iLambdaS;
666
667   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
668   // initialize (coefficient = alpha, rates = lambda)
669
670   fLambdaL = 0; 
671   fAlphaL  = 0; 
672   fLambdaS = 0; 
673   fAlphaS  = 0;
674   iLambdaL = 0; 
675   iAlphaL  = 0; 
676   iLambdaS = 0; 
677   iAlphaS  = 0;
678
679   Double_t dt = 0.1;
680
681   Double_t r1;
682   Double_t r2;
683   Double_t c1;
684   Double_t c2;
685   r1 = (Double_t) fR1;
686   r2 = (Double_t) fR2;
687   c1 = (Double_t) fC1;
688   c2 = (Double_t) fC2;
689
690   fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
691   fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
692   iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600;     //  9 bit paramter + fixed bits
693   iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200;     //  9 bit paramter + fixed bits
694
695   if (nexp == 1) {
696     fAlphaL = (Int_t) (c1 * 2048.0);
697     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
698   }
699   if (nexp == 2) {
700     fAlphaL = (Int_t) (c1 * 2048.0);
701     fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
702     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
703     iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400;              // 10 bit paramter + fixed bits
704   }
705   
706   Double_t iAl = iAlphaL  / 2048.0;            // alpha L: correspondence to floating point numbers
707   Double_t iAs = iAlphaS  / 2048.0;            // alpha S: correspondence to floating point numbers
708   Double_t iLl = iLambdaL / 2048.0;            // lambda L: correspondence to floating point numbers
709   Double_t iLs = iLambdaS / 2048.0;            // lambda S: correspondence to floating point numbers
710
711   Int_t h1;
712   Int_t h2;
713   Int_t rem1;
714   Int_t rem2;
715   Int_t correction;
716   Int_t result;
717   Int_t iFactor = ((Int_t) fPedestal) << 2;
718
719   Double_t xi = 1 - (iLl*iAs + iLs*iAl);             // Calculation of equilibrium values of the
720   rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
721   rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
722   
723   // further initialization
724   if ((rem1 + rem2) > 0x0FFF) {
725     correction = 0x0FFF;
726   } 
727   else {
728     correction = (rem1 + rem2) & 0x0FFF;
729   }
730
731   fTailPed = iFactor - correction;
732
733   for (i = 0; i < n; i++) {
734
735     result = ((Int_t)source[i] - correction);
736     if (result < 0) {                           
737       result = 0;
738     }
739
740     target[i] = result;
741                                                         
742     h1 = (rem1 + ((iAlphaL * result) >> 11));
743     if (h1 > 0x0FFF) {
744       h1 = 0x0FFF;
745     } 
746     else {
747       h1 &= 0x0FFF;
748     }
749
750     h2 = (rem2 + ((iAlphaS * result) >> 11));
751     if (h2 > 0x0FFF) {
752       h2 = 0x0FFF;
753     } 
754     else {
755       h2 &= 0x0FFF;
756     }
757   
758     rem1 = (iLambdaL * h1 ) >> 11;
759     rem2 = (iLambdaS * h2 ) >> 11;
760     
761     if ((rem1 + rem2) > 0x0FFF) {
762       correction = 0x0FFF;
763     } 
764     else {
765       correction = (rem1 + rem2) & 0x0FFF;
766     }
767
768   }
769
770 }
771
772 //_____________________________________________________________________________
773 void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n) 
774 {
775   //
776   // Exponential filter (M. Ivanov)
777   //
778
779   Int_t i = 0;
780
781   Double_t sig1[100];
782   Double_t sig2[100];
783   Double_t sig3[100];
784
785   for (i = 0; i < n; i++) {
786     sig1[i] = source[i];
787   }
788
789   Float_t dt = 0.1;
790
791   Float_t lambda0 = (1.0 / fR2) * dt;
792   Float_t lambda1 = (1.0 / fR1) * dt;
793
794   TailMakerSpline(sig1,sig2,lambda0,n);
795   TailCancelationMI(sig2,sig3,0.7,lambda1,n);
796
797   for (i = 0; i < n; i++) {
798     target[i] = sig3[i];
799   }
800
801 }
802
803 //______________________________________________________________________________
804 void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n) 
805 {
806   //
807   // Special filter (M. Ivanov)
808   //
809
810   Int_t    i = 0;
811
812   Double_t l = TMath::Exp(-lambda*0.5);
813   Double_t in[1000];
814   Double_t out[1000];
815
816   // Initialize in[] and out[] goes 0 ... 2*n+19
817   for (i = 0; i < n*2+20; i++) {
818     in[i]  = 0;
819     out[i] = 0;
820   }
821
822   // in[] goes 0, 1
823   in[0] = ampin[0];
824   in[1] = (ampin[0] + ampin[1]) * 0.5;
825    
826   // Add charge to the end
827   for (i = 0; i < 22; i++) {
828     // in[] goes 2*n-2, 2*n-1, ... , 2*n+19 
829     in[2*(n-1)+i] = ampin[n-1];
830   }
831
832   // Use arithmetic mean
833   for (i = 1; i < n-1; i++) {
834     // in[] goes 2, 3, ... , 2*n-4, 2*n-3
835     in[2*i]   = ampin[i];
836     in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
837   }
838
839   Double_t temp;
840   out[2*n]    = in[2*n];
841   temp        = 0;
842   for (i = 2*n; i >= 0; i--) {
843     out[i]    = in[i] + temp;
844     temp      = l*(temp+in[i]);
845   }
846
847   for (i = 0; i < n; i++){
848     //ampout[i] = out[2*i+1];  // org
849     ampout[i] = out[2*i];
850   }
851
852 }
853
854 //______________________________________________________________________________
855 void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout
856                                 , Double_t norm, Double_t lambda, Int_t n) 
857 {
858   //
859   // Special filter (M. Ivanov)
860   //
861
862   Int_t    i = 0;
863
864   Double_t l = TMath::Exp(-lambda*0.5);
865   Double_t k = l*(1.0 - norm*lambda*0.5);
866   Double_t in[1000];
867   Double_t out[1000];
868
869   // Initialize in[] and out[] goes 0 ... 2*n+19
870   for (i = 0; i < n*2+20; i++) {
871     in[i]  = 0;
872     out[i] = 0;
873   }
874
875   // in[] goes 0, 1
876   in[0] = ampin[0];
877   in[1] = (ampin[0]+ampin[1])*0.5;
878
879   // Add charge to the end
880   for (i =-2; i < 22; i++) {
881     // in[] goes 2*n-4, 2*n-3, ... , 2*n+19 
882     in[2*(n-1)+i] = ampin[n-1];
883   }
884
885   for (i = 1; i < n-2; i++) {
886     // in[] goes 2, 3, ... , 2*n-6, 2*n-5
887     in[2*i]    = ampin[i];
888     in[2*i+1]  = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
889     //in[2*i+1]  = ((ampin[i]+ampin[i+1]))/2.0;
890   }
891
892   Double_t temp;
893   out[0] = in[0];
894   temp   = in[0];
895   for (i = 1; i <= 2*n; i++) {
896     out[i] = in[i] + (k-l)*temp;
897     temp   = in[i] +  k   *temp;
898   }
899
900   for (i = 0; i < n; i++) {
901     //ampout[i] = out[2*i+1];  // org
902     //ampout[i] = TMath::Max(out[2*i+1],0.0);  // org
903     ampout[i] = TMath::Max(out[2*i],0.0);
904   }
905
906 }
907