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