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