Ssemi-final version of TRD raw data simulation
[u/mrichter/AliRoot.git] / TRD / AliTRDmcmSim.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 */
19
20 ///////////////////////////////////////////////////////////////////////////////
21 //                                                                           //
22 //  TRD MCM (Multi Chip Module) simulator                                    //
23 //                                                                           //
24 ///////////////////////////////////////////////////////////////////////////////
25
26 #include <TMath.h>
27
28 #include "AliLog.h"
29 #include "AliTRDmcmSim.h"
30 #include "AliTRDfeeParam.h"
31 #include "AliTRDgeometry.h"
32
33 ClassImp(AliTRDmcmSim)
34
35 //_____________________________________________________________________________
36 AliTRDmcmSim::AliTRDmcmSim() :TObject()
37   ,fInitialized(kFALSE)
38   ,fFeeParam(NULL)
39   ,fGeo(NULL)
40   ,fChaId(-1)
41   ,fSector(-1)
42   ,fStack(-1)
43   ,fLayer(-1)
44   ,fRobPos(-1)
45   ,fMcmPos(-1)
46   ,fNADC(-1)
47   ,fNTimeBin(-1)
48   ,fRow(-1)
49   ,fColOfADCbeg(-1)
50   ,fColOfADCend(-1)
51   ,fADCR(NULL)
52   ,fADCF(NULL)
53   ,fZSM(NULL)
54   ,fZSM1Dim(NULL)
55 {
56   //
57   // AliTRDmcmSim default constructor
58   //
59
60   // By default, nothing is initialized.
61   // It is necessary to issue Init before use.
62 }
63
64 //_____________________________________________________________________________
65 AliTRDmcmSim::~AliTRDmcmSim() 
66 {
67   //
68   // AliTRDmcmSim destructor
69   //
70   if( fADCR != NULL ) {
71     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
72       delete fADCR[iadc];
73       delete fADCF[iadc];
74       delete fZSM [iadc];
75     }
76     delete fADCR;
77     delete fADCF;
78     delete fZSM;
79     delete fZSM1Dim;
80   }
81   delete fGeo;
82 }
83
84 //_____________________________________________________________________________
85 void AliTRDmcmSim::Init( Int_t cha_id, Int_t rob_pos, Int_t mcm_pos )
86 {
87   // Initialize the class with new geometry information
88   // fADC array will be reused with filled by zero
89
90   fFeeParam      = AliTRDfeeParam::Instance();
91   fGeo           = new AliTRDgeometry();
92   fChaId         = cha_id;
93   fSector        = fGeo->GetSector( fChaId );
94   fStack         = fGeo->GetChamber( fChaId );
95   fLayer         = fGeo->GetPlane( fChaId );
96   fRobPos        = rob_pos;
97   fMcmPos        = mcm_pos;
98   fNADC          = fFeeParam->GetNadcMcm();
99   fNTimeBin      = fFeeParam->GetNtimebin();
100   fRow           = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
101   fColOfADCbeg   = fFeeParam->GetPadColFromADC( fRobPos, fMcmPos, 0 );
102   fColOfADCend   = fFeeParam->GetPadColFromADC( fRobPos, fMcmPos, fNADC-1 );
103
104   // Allocate ADC data memory if not yet done
105   if( fADCR == NULL ) {
106     fADCR    = new Int_t *[fNADC];
107     fADCF    = new Int_t *[fNADC];
108     fZSM     = new Int_t *[fNADC];
109     fZSM1Dim = new Int_t  [fNADC];
110     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
111       fADCR[iadc] = new Int_t[fNTimeBin];
112       fADCF[iadc] = new Int_t[fNTimeBin];
113       fZSM [iadc] = new Int_t[fNTimeBin];
114     }
115   }
116
117   // Initialize ADC data
118   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
119     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
120       fADCR[iadc][it] = 0;
121       fADCF[iadc][it] = 0;
122       fZSM [iadc][it] = 1;   // Default unread = 1
123     }
124     fZSM1Dim[iadc] = 1;      // Default unread = 1
125   }
126
127   fInitialized = kTRUE;
128 }
129
130 //_____________________________________________________________________________
131 void AliTRDmcmSim::SetData( Int_t iadc, Int_t *adc )
132 {
133   // Store ADC data into array of raw data
134
135   if( ! fInitialized ) {
136     //Log (Form ("Error: AliTRDmcmSim is not initialized but setData is called."));
137     return;
138   }
139
140   if( iadc < 0 || iadc >= fNADC ) {
141     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
142     return;
143   }
144
145   for( int it = 0 ;  it < fNTimeBin ; it++ ) {
146     fADCR[iadc][it] = (Int_t)(adc[it]);
147   }
148 }
149
150 //_____________________________________________________________________________
151 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
152 {
153   // Store ADC data into array of raw data
154
155   if( ! fInitialized ) {
156     //Log (Form ("Error: AliTRDmcmSim is not initialized but setData is called."));
157     return;
158   }
159
160   if( iadc < 0 || iadc >= fNADC ) {
161     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
162     return;
163   }
164
165   fADCR[iadc][it] = adc;
166 }
167
168 //_____________________________________________________________________________
169 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
170 {
171   // Store ADC data into array of raw data
172
173   if( ! fInitialized ) {
174     //Log (Form ("Error: AliTRDmcmSim is not initialized but setData is called."));
175     return;
176   }
177
178   if( iadc < 0 || iadc >= fNADC ) {
179     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
180     return;
181   }
182
183   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
184     fADCR[iadc][it] = fFeeParam->GetADCpedestal();
185   }
186 }
187
188 //_____________________________________________________________________________
189 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
190 {
191   // Return column id of the pad for the given ADC channel
192
193   return (fColOfADCbeg - iadc);
194 }
195
196
197
198
199 //_____________________________________________________________________________
200 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize )
201 {
202   // Produce raw data stream from this MCM and put in buf
203   // Returns number of words filled, or negative value with -1 * number of overflowed words
204
205   UInt_t  x;
206   UInt_t  iEv = 0;
207   Int_t   nw  = 0;  // Number of written words
208   Int_t   of  = 0;  // Number of overflowed words
209   Int_t   rawVer   = fFeeParam->GetRAWversion();
210   Int_t **adc;
211
212   if( fFeeParam->GetRAWstoreRaw() ) {
213     adc = fADCR;
214   } else {
215     adc = fADCF;
216   }
217
218   // Produce MCM header
219   x = ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
220   if (nw < maxSize) {
221     buf[nw++] = x;
222   }
223   else {
224     of++;
225   }
226
227   // Produce ADC mask
228   if( rawVer >= 3 ) {
229     x = 0;
230     for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
231       if( fZSM1Dim[iAdc] == 0 ) { //  0 means not suppressed
232         x = x | (1 << iAdc);
233       }
234     }
235     if (nw < maxSize) {
236       buf[nw++] = x;
237     }
238     else {
239       of++;
240     }
241   }
242
243   // Produce ADC data. 3 timebins are packed into one 32 bits word
244   // In this version, different ADC channel will NOT share the same word
245
246   UInt_t aa=0, a1=0, a2=0, a3=0;
247
248   for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
249     if( rawVer>= 3 && fZSM1Dim[iAdc] == 1 ) continue; // suppressed
250     aa = !(iAdc & 1) + 2;
251     for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
252       a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] : 0;
253       a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
254       a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
255     }
256     x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
257     if (nw < maxSize) {
258       buf[nw++] = x;
259     }
260     else {
261       of++;
262     }
263   }
264
265   if( of != 0 ) return -of; else return nw;
266 }
267
268
269 //_____________________________________________________________________________
270 void AliTRDmcmSim::Filter()
271 {
272   // Apply digital filter
273
274   if( ! fInitialized ) {
275     // Log (Form ("Error: AliTRDmcmSim is not initialized but setData is called."));
276     return;
277   }
278
279   // Initialize filtered data array with raw data
280   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
281     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
282       fADCF[iadc][it] = fADCR[iadc][it]; 
283     }
284   }
285
286   // Then apply fileters one by one to filtered data array
287   if( fFeeParam->isPFon() ) FilterPedestal();
288   if( fFeeParam->isGFon() ) FilterGain();
289   if( fFeeParam->isTFon() ) FilterTail();
290 }
291
292 //_____________________________________________________________________________
293 void AliTRDmcmSim::FilterPedestal()
294 {
295   // Apply pedestal filter
296
297   Int_t ap = fFeeParam->GetADCpedestal();      // ADC instrinsic pedestal
298   Int_t ep = fFeeParam->GetPFeffectPedestal(); // effective pedestal
299   Int_t tc = fFeeParam->GetPFtimeConstant();   // this makes no sense yet
300
301   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
302     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
303       fADCF[iadc][it] = fADCF[iadc][it] - ap + ep;
304     }
305   }
306 }
307
308 //_____________________________________________________________________________
309 void AliTRDmcmSim::FilterGain()
310 {
311   // Apply gain filter (not implemented)
312 }
313
314 //_____________________________________________________________________________
315 void AliTRDmcmSim::FilterTail()
316 {
317   // Apply exponential tail filter (Bogdan's version)
318
319   Double_t *dtarg  = new Double_t[fNTimeBin];
320   Int_t    *itarg  = new Int_t[fNTimeBin];
321   Int_t     nexp   = fFeeParam->GetTFnExp();
322   Int_t     tftype = fFeeParam->GetTFtype();
323
324   switch( tftype ) {
325     
326   case 0: // Exponential Filter Analog Bogdan
327     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
328       FilterSimDeConvExpA( fADCF[iCol], dtarg, fNTimeBin, nexp);
329       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
330         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
331       }
332     }
333     break;
334
335   case 1: // Exponential filter digital Bogdan
336     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
337       FilterSimDeConvExpD( fADCF[iCol], itarg, fNTimeBin, nexp);
338       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
339         fADCF[iCol][iTime] = itarg[iTime];
340       }
341     }
342     break;
343     
344   case 2: // Exponential filter Marian special
345     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
346       FilterSimDeConvExpMI( fADCF[iCol], dtarg, fNTimeBin);
347       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
348         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
349       }
350     }
351     break;
352     
353   default:
354     AliError(Form("Invalid filter type %d ! \n", tftype ));
355     break;
356   }
357
358   delete dtarg;
359   delete itarg;
360 }
361
362 //_____________________________________________________________________________
363 void AliTRDmcmSim::ZSMapping()
364 {
365   // Zero Suppression Mapping implemented in TRAP chip
366   //
367   // See detail TRAP manual "Data Indication" section:
368   // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
369
370   Int_t EBIS = fFeeParam->GetEBsglIndThr();       // TRAP default = 0x4  (Tis=4)
371   Int_t EBIT = fFeeParam->GetEBsumIndThr();       // TRAP default = 0x28 (Tit=40)
372   Int_t EBIL = fFeeParam->GetEBindLUT();          // TRAP default = 0xf0 (lookup table accept (I2,I1,I0)=(111) or (110) or (101) or (100))
373   Int_t EBIN = fFeeParam->GetEBignoreNeighbour(); // TRAP default = 1 (no neighbor sensitivity)
374
375   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
376     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
377
378       // evaluate three conditions
379       Int_t I0 = ( fADCF[iadc][it] >=  fADCF[iadc-1][it] && fADCF[iadc][it] >=  fADCF[iadc+1][it] ) ? 0 : 1; // peak center detection
380       Int_t I1 = ( (fADCF[iadc-1][it] + fADCF[iadc][it] + fADCF[iadc+1][it]) > EBIT )               ? 0 : 1; // cluster
381       Int_t I2 = ( fADCF[iadc][it] > EBIS )                                                         ? 0 : 1; // absolute large peak
382
383       Int_t I = I2 * 4 + I1 * 2 + I0; // Bit position in lookup table
384       Int_t D = (EBIL >> I) & 1;      // Looking up  (here D=0 means true and D=1 means false according to TRAP manual)
385
386       fZSM[iadc][it] &= D;
387       if( EBIN == 0 ) {  // turn on neighboring ADCs
388         fZSM[iadc-1][it] &= D;
389         fZSM[iadc+1][it] &= D;
390       }
391     }
392   }
393
394   // do 1 dim projection
395   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
396     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
397       fZSM1Dim[iadc] &= fZSM[iadc+1][it];
398     }
399   }
400 }
401
402 //_____________________________________________________________________________
403 void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target, Int_t n, Int_t nexp) 
404 {
405   //
406   // Exponential filter "analog"
407   // source will not be changed
408   //
409
410   Int_t    i = 0;
411   Int_t    k = 0;
412   Double_t reminder[2];
413   Double_t correction;
414   Double_t result;
415   Double_t rates[2];
416   Double_t coefficients[2];
417
418   // Initialize (coefficient = alpha, rates = lambda)
419   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
420
421   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
422   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
423   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
424   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
425   
426   coefficients[0] = c1;
427   coefficients[1] = c2;
428
429   Double_t dt = 0.1;
430   rates[0] = TMath::Exp(-dt/(r1));
431   rates[1] = TMath::Exp(-dt/(r2));
432
433   // Attention: computation order is important
434   correction = 0.0;
435   for (k = 0; k < nexp; k++) {
436     reminder[k] = 0.0;
437   }
438     
439   for (i = 0; i < n; i++) {
440
441     result    = ((Double_t)source[i] - correction);    // no rescaling
442     target[i] = result;
443     
444     for (k = 0; k < nexp; k++) {
445       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
446     }
447       
448     correction = 0.0;
449     for (k = 0; k < nexp; k++) {
450       correction += reminder[k];
451     }
452   }
453 }
454
455 //_____________________________________________________________________________
456 void AliTRDmcmSim::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n, Int_t nexp) 
457 {
458   //
459   // Exponential filter "digital"
460   // source will not be changed
461   //
462
463   Int_t i = 0;
464
465   Int_t fAlphaL  = 0;
466   Int_t fAlphaS  = 0;
467   Int_t fLambdaL = 0;
468   Int_t fLambdaS = 0;
469   Int_t fTailPed = 0;
470
471   Int_t iAlphaL  = 0;
472   Int_t iAlphaS  = 0;
473   Int_t iLambdaL = 0;
474   Int_t iLambdaS = 0;
475
476   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
477   // initialize (coefficient = alpha, rates = lambda)
478
479   Double_t dt = 0.1;
480   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
481   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
482   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
483   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
484
485   fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
486   fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
487   iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600;     //  9 bit paramter + fixed bits
488   iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200;     //  9 bit paramter + fixed bits
489
490   if (nexp == 1) {
491     fAlphaL = (Int_t) (c1 * 2048.0);
492     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
493   }
494   if (nexp == 2) {
495     fAlphaL = (Int_t) (c1 * 2048.0);
496     fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
497     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
498     iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400;              // 10 bit paramter + fixed bits
499   }
500   
501   Double_t iAl = iAlphaL  / 2048.0;            // alpha  L: correspondence to floating point numbers
502   Double_t iAs = iAlphaS  / 2048.0;            // alpha  S: correspondence to floating point numbers
503   Double_t iLl = iLambdaL / 2048.0;            // lambda L: correspondence to floating point numbers
504   Double_t iLs = iLambdaS / 2048.0;            // lambda S: correspondence to floating point numbers
505
506   Int_t h1;
507   Int_t h2;
508   Int_t rem1;
509   Int_t rem2;
510   Int_t correction;
511   Int_t result;
512   Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2;
513
514   Double_t xi = 1 - (iLl*iAs + iLs*iAl);             // Calculation of equilibrium values of the
515   rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
516   rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
517   
518   // further initialization
519   if ((rem1 + rem2) > 0x0FFF) {
520     correction = 0x0FFF;
521   } 
522   else {
523     correction = (rem1 + rem2) & 0x0FFF;
524   }
525
526   fTailPed = iFactor - correction;
527
528   for (i = 0; i < n; i++) {
529
530     result = (source[i] - correction);
531     if (result < 0) {                           
532       result = 0;
533     }
534
535     target[i] = result;
536                                                         
537     h1 = (rem1 + ((iAlphaL * result) >> 11));
538     if (h1 > 0x0FFF) {
539       h1 = 0x0FFF;
540     } 
541     else {
542       h1 &= 0x0FFF;
543     }
544
545     h2 = (rem2 + ((iAlphaS * result) >> 11));
546     if (h2 > 0x0FFF) {
547       h2 = 0x0FFF;
548     } 
549     else {
550       h2 &= 0x0FFF;
551     }
552   
553     rem1 = (iLambdaL * h1 ) >> 11;
554     rem2 = (iLambdaS * h2 ) >> 11;
555     
556     if ((rem1 + rem2) > 0x0FFF) {
557       correction = 0x0FFF;
558     } 
559     else {
560       correction = (rem1 + rem2) & 0x0FFF;
561     }
562
563   }
564 }
565
566 //_____________________________________________________________________________
567 void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target, Int_t n) 
568 {
569   //
570   // Exponential filter (M. Ivanov)
571   // source will not be changed
572   //
573
574   Int_t i = 0;
575   Double_t sig1[100];
576   Double_t sig2[100];
577   Double_t sig3[100];
578
579   for (i = 0; i < n; i++) {
580     sig1[i] = (Double_t)source[i];
581   }
582
583   Float_t dt      = 0.1;
584   Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt;
585   Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt;
586
587   FilterSimTailMakerSpline( sig1, sig2, lambda0, n);
588   FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n);
589
590   for (i = 0; i < n; i++) {
591     target[i] = sig3[i];
592   }
593
594 }
595
596 //______________________________________________________________________________
597 void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n) 
598 {
599   //
600   // Special filter (M. Ivanov)
601   //
602
603   Int_t    i = 0;
604   Double_t l = TMath::Exp(-lambda*0.5);
605   Double_t in[1000];
606   Double_t out[1000];
607
608   // Initialize in[] and out[] goes 0 ... 2*n+19
609   for (i = 0; i < n*2+20; i++) {
610     in[i] = out[i] = 0;
611   }
612
613   // in[] goes 0, 1
614   in[0] = ampin[0];
615   in[1] = (ampin[0] + ampin[1]) * 0.5;
616    
617   // Add charge to the end
618   for (i = 0; i < 22; i++) {
619     in[2*(n-1)+i] = ampin[n-1]; // in[] goes 2*n-2, 2*n-1, ... , 2*n+19 
620   }
621
622   // Use arithmetic mean
623   for (i = 1; i < n-1; i++) {
624     in[2*i]   = ampin[i];    // in[] goes 2, 3, ... , 2*n-4, 2*n-3
625     in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
626   }
627
628   Double_t temp;
629   out[2*n]    = in[2*n];
630   temp        = 0;
631   for (i = 2*n; i >= 0; i--) {
632     out[i]    = in[i] + temp;
633     temp      = l*(temp+in[i]);
634   }
635
636   for (i = 0; i < n; i++){
637     //ampout[i] = out[2*i+1];  // org
638     ampout[i] = out[2*i];
639   }
640
641 }
642
643 //______________________________________________________________________________
644 void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout, Double_t norm, Double_t lambda, Int_t n) 
645 {
646   //
647   // Special filter (M. Ivanov)
648   //
649
650   Int_t    i = 0;
651
652   Double_t l = TMath::Exp(-lambda*0.5);
653   Double_t k = l*(1.0 - norm*lambda*0.5);
654   Double_t in[1000];
655   Double_t out[1000];
656
657   // Initialize in[] and out[] goes 0 ... 2*n+19
658   for (i = 0; i < n*2+20; i++) {
659     in[i] = out[i] = 0;
660   }
661
662   // in[] goes 0, 1
663   in[0] = ampin[0];
664   in[1] = (ampin[0]+ampin[1])*0.5;
665
666   // Add charge to the end
667   for (i =-2; i < 22; i++) {
668     // in[] goes 2*n-4, 2*n-3, ... , 2*n+19 
669     in[2*(n-1)+i] = ampin[n-1];
670   }
671
672   for (i = 1; i < n-2; i++) {
673     // in[] goes 2, 3, ... , 2*n-6, 2*n-5
674     in[2*i]    = ampin[i];
675     in[2*i+1]  = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
676     //in[2*i+1]  = ((ampin[i]+ampin[i+1]))/2.0;
677   }
678
679   Double_t temp;
680   out[0] = in[0];
681   temp   = in[0];
682   for (i = 1; i <= 2*n; i++) {
683     out[i] = in[i] + (k-l)*temp;
684     temp   = in[i] +  k   *temp;
685   }
686
687   for (i = 0; i < n; i++) {
688     //ampout[i] = out[2*i+1];  // org
689     //ampout[i] = TMath::Max(out[2*i+1],0.0);  // org
690     ampout[i] = TMath::Max(out[2*i],0.0);
691   }
692 }
693
694 // EOF