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