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