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