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 //                                                                           //
  18 //  TRD MCM (Multi Chip Module) simulator                                    //
  19 //                                                                           //
  20 ///////////////////////////////////////////////////////////////////////////////
  21
  22 /* $Id$ */
  23
  24 /*
  25
  26   New release on 2007/08/17
  27
  28 AliTRDmcmSim is now stably working and zero suppression function seems ok.
  29 From now, the default version of raw data is set to 3 in AliTRDfeeParam.
  30
  31 The following internal parameters were abolished because it is useless and
  32 made trouble:
  33
  34    fColOfADCbeg
  35    fColOfADCend
  36
  37 GetCol member was modified accordingly.
  38
  39 New member function DumpData was prepared for diagnostics.
  40
  41 ZSMapping member function was debugged. It was causing crash due to
  42 wrong indexing in 1 dimensional numbering. Also code was shaped up better.
  43
  44 */
  45
  46 /*Semi-final version of TRD raw data simulation code with zero suppression (ZS)
  47 similar to TRD FEE. ZS is realized by the class group:
  48
  49   AliTRDfeeParam
  50   AliTRDmcmSim
  51   AliTRDrawData
  52
  53 AliTRDfeeParam has been modified to have more parameters like raw data
  54 production version and so on. AliTRDmcmSim is new class and this is the core
  55 of MCM (PASA+TRAP) simulator. It has still very simple function and it will be
  56 another project to improve this to make it closer to the reall FEE.
  57 AliTRDrawData has been modified to use new class AliTRDmcmSim.
  58
  59 These modifications were tested on Aug. 02 HEAD version that code itself
  60 compiles. I'm sure there must be still bugs and we need testing by as many as
  61 possible persons now. Especially it seems HLT part is impacted by problems
  62 because some parameters were moved from AliTRDrawData to AliTRDfeeParam (like
  63 fRawVersion disappeared from AliTRDrawData).
  64
  65 In TRD definition, we have now 4 raw data versions.
  66
  67   0 very old offline version (by Bogdan)
  68   1 test version (no zero suppression)
  69   2 final version (no zero suppression)
  70   3 test version (with zero suppression)
  71
  72 The default is still set to 2 in AliTRDfeeParam::fgkRAWversion and it uses
  73 previously existing codes. If you set this to 3, AliTRDrawData changes behavior
  74 to use AliTRDmcmSim with ZS.
  75
  76 Plan is after we make sure it works stably, we delete AliTRDmcm which is obsolete.
  77 However it still take time because tarcklet part is not yet touched.
  78 The default raw version is 2.
  79
  80                                                                  Ken Oyama
  81 */
  82
  83 #include <fstream>
  84 #include <TMath.h>
  85 #include "AliLog.h"
  86 #include "AliTRDmcmSim.h"
  87 #include "AliTRDfeeParam.h"
  88 #include "AliTRDSimParam.h"
  89 #include "AliTRDgeometry.h"
  90 #include "AliTRDcalibDB.h"
  91
  92 ClassImp(AliTRDmcmSim)
  93
  94 //_____________________________________________________________________________
  95 AliTRDmcmSim::AliTRDmcmSim() :TObject()
  96   ,fInitialized(kFALSE)
  97   ,fFeeParam(NULL)
  98   ,fSimParam(NULL)
  99   ,fCal(NULL)
 100   ,fGeo(NULL)
 101   ,fChaId(-1)
 102   ,fSector(-1)
 103   ,fStack(-1)
 104   ,fLayer(-1)
 105   ,fRobPos(-1)
 106   ,fMcmPos(-1)
 107   ,fNADC(-1)
 108   ,fNTimeBin(-1)
 109   ,fRow(-1)
 110   ,fADCR(NULL)
 111   ,fADCF(NULL)
 112   ,fZSM(NULL)
 113   ,fZSM1Dim(NULL)
 114 {
 115   //
 116   // AliTRDmcmSim default constructor
 117   //
 118
 119   // By default, nothing is initialized.
 120   // It is necessary to issue Init before use.
 121 }
 122
 123 //_____________________________________________________________________________
 124 AliTRDmcmSim::~AliTRDmcmSim()
 125 {
 126   //
 127   // AliTRDmcmSim destructor
 128   //
 129   if( fADCR != NULL ) {
 130     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 131       delete fADCR[iadc];
 132       delete fADCF[iadc];
 133       delete fZSM [iadc];
 134     }
 135     delete fADCR;
 136     delete fADCF;
 137     delete fZSM;
 138     delete fZSM1Dim;
 139   }
 140   delete fGeo;
 141 }
 142
 143 //_____________________________________________________________________________
 144 void AliTRDmcmSim::Init( Int_t cha_id, Int_t rob_pos, Int_t mcm_pos )
 145 {
 146   // Initialize the class with new geometry information
 147   // fADC array will be reused with filled by zero
 148
 149   fFeeParam      = AliTRDfeeParam::Instance();
 150   fSimParam      = AliTRDSimParam::Instance();
 151   fCal           = AliTRDcalibDB::Instance();
 152   fGeo           = new AliTRDgeometry();
 153   fChaId         = cha_id;
 154   fSector        = fGeo->GetSector( fChaId );
 155   fStack         = fGeo->GetChamber( fChaId );
 156   fLayer         = fGeo->GetPlane( fChaId );
 157   fRobPos        = rob_pos;
 158   fMcmPos        = mcm_pos;
 159   fNADC          = fFeeParam->GetNadcMcm();
 160   fNTimeBin      = fCal->GetNumberOfTimeBins();
 161   fRow           = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
 162
 163   // Allocate ADC data memory if not yet done
 164   if( fADCR == NULL ) {
 165     fADCR    = new Int_t *[fNADC];
 166     fADCF    = new Int_t *[fNADC];
 167     fZSM     = new Int_t *[fNADC];
 168     fZSM1Dim = new Int_t  [fNADC];
 169     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 170       fADCR[iadc] = new Int_t[fNTimeBin];
 171       fADCF[iadc] = new Int_t[fNTimeBin];
 172       fZSM [iadc] = new Int_t[fNTimeBin];
 173     }
 174   }
 175
 176   // Initialize ADC data
 177   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 178     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 179       fADCR[iadc][it] = 0;
 180       fADCF[iadc][it] = 0;
 181       fZSM [iadc][it] = 1;   // Default unread = 1
 182     }
 183     fZSM1Dim[iadc] = 1;      // Default unread = 1
 184   }
 185
 186   fInitialized = kTRUE;
 187 }
 188
 189 //_____________________________________________________________________________
 190 Bool_t AliTRDmcmSim::CheckInitialized()
 191 {
 192   if( ! fInitialized ) {
 193     AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
 194   }
 195   return fInitialized;
 196 }
 197
 198 //_____________________________________________________________________________
 199 void AliTRDmcmSim::SetData( Int_t iadc, Int_t *adc )
 200 {
 201   // Store ADC data into array of raw data
 202
 203   if( !CheckInitialized() ) return;
 204
 205   if( iadc < 0 || iadc >= fNADC ) {
 206     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
 207     return;
 208   }
 209
 210   for( int it = 0 ;  it < fNTimeBin ; it++ ) {
 211     fADCR[iadc][it] = (Int_t)(adc[it]);
 212   }
 213 }
 214
 215 //_____________________________________________________________________________
 216 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
 217 {
 218   // Store ADC data into array of raw data
 219
 220   if( !CheckInitialized() ) return;
 221
 222   if( iadc < 0 || iadc >= fNADC ) {
 223     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
 224     return;
 225   }
 226
 227   fADCR[iadc][it] = adc;
 228 }
 229
 230 //_____________________________________________________________________________
 231 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
 232 {
 233   // Store ADC data into array of raw data
 234
 235   if( !CheckInitialized() ) return;
 236
 237   if( iadc < 0 || iadc >= fNADC ) {
 238     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
 239     return;
 240   }
 241
 242   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 243     fADCR[iadc][it] = fSimParam->GetADCbaseline();
 244   }
 245 }
 246
 247 //_____________________________________________________________________________
 248 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
 249 {
 250   // Return column id of the pad for the given ADC channel
 251   if( !CheckInitialized() ) return -1;
 252
 253   return fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
 254 }
 255
 256
 257
 258
 259 //_____________________________________________________________________________
 260 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize )
 261 {
 262   // Produce raw data stream from this MCM and put in buf
 263   // Returns number of words filled, or negative value with -1 * number of overflowed words
 264
 265   UInt_t  x;
 266   UInt_t  iEv = 0;
 267   Int_t   nw  = 0;  // Number of written words
 268   Int_t   of  = 0;  // Number of overflowed words
 269   Int_t   rawVer   = fFeeParam->GetRAWversion();
 270   Int_t **adc;
 271
 272   if( !CheckInitialized() ) return 0;
 273
 274   if( fFeeParam->GetRAWstoreRaw() ) {
 275     adc = fADCR;
 276   } else {
 277     adc = fADCF;
 278   }
 279
 280   // Produce MCM header
 281   x = ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
 282   if (nw < maxSize) {
 283     buf[nw++] = x;
 284   }
 285   else {
 286     of++;
 287   }
 288
 289   // Produce ADC mask
 290   if( rawVer >= 3 ) {
 291     x = 0;
 292     for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
 293       if( fZSM1Dim[iAdc] == 0 ) { //  0 means not suppressed
 294         x = x | (1 << iAdc);
 295       }
 296     }
 297     if (nw < maxSize) {
 298       buf[nw++] = x;
 299     }
 300     else {
 301       of++;
 302     }
 303   }
 304
 305   // Produce ADC data. 3 timebins are packed into one 32 bits word
 306   // In this version, different ADC channel will NOT share the same word
 307
 308   UInt_t aa=0, a1=0, a2=0, a3=0;
 309
 310   for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
 311     if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // suppressed
 312     aa = !(iAdc & 1) + 2;
 313     for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
 314       a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] : 0;
 315       a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
 316       a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
 317       x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
 318       if (nw < maxSize) {
 319         buf[nw++] = x;
 320       }
 321       else {
 322         of++;
 323       }
 324     }
 325   }
 326
 327   if( of != 0 ) return -of; else return nw;
 328 }
 329
 330
 331 //_____________________________________________________________________________
 332 void AliTRDmcmSim::Filter()
 333 {
 334   // Apply digital filter
 335
 336   if( !CheckInitialized() ) return;
 337
 338   // Initialize filtered data array with raw data
 339   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 340     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 341       fADCF[iadc][it] = fADCR[iadc][it];
 342     }
 343   }
 344
 345   // Then apply fileters one by one to filtered data array
 346   if( fFeeParam->isPFon() ) FilterPedestal();
 347   if( fFeeParam->isGFon() ) FilterGain();
 348   if( fFeeParam->isTFon() ) FilterTail();
 349 }
 350
 351 //_____________________________________________________________________________
 352 void AliTRDmcmSim::FilterPedestal()
 353 {
 354   // Apply pedestal filter
 355
 356   Int_t ap = fSimParam->GetADCbaseline();      // ADC instrinsic pedestal
 357   Int_t ep = fFeeParam->GetPFeffectPedestal(); // effective pedestal
 358   Int_t tc = fFeeParam->GetPFtimeConstant();   // this makes no sense yet
 359
 360   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 361     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 362       fADCF[iadc][it] = fADCF[iadc][it] - ap + ep;
 363     }
 364   }
 365 }
 366
 367 //_____________________________________________________________________________
 368 void AliTRDmcmSim::FilterGain()
 369 {
 370   // Apply gain filter (not implemented)
 371   // Later it will be implemented because gain digital filiter will
 372   // increase noise level.
 373 }
 374
 375 //_____________________________________________________________________________
 376 void AliTRDmcmSim::FilterTail()
 377 {
 378   // Apply exponential tail filter (Bogdan's version)
 379
 380   Double_t *dtarg  = new Double_t[fNTimeBin];
 381   Int_t    *itarg  = new Int_t[fNTimeBin];
 382   Int_t     nexp   = fFeeParam->GetTFnExp();
 383   Int_t     tftype = fFeeParam->GetTFtype();
 384
 385   switch( tftype ) {
 386
 387   case 0: // Exponential Filter Analog Bogdan
 388     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 389       FilterSimDeConvExpA( fADCF[iCol], dtarg, fNTimeBin, nexp);
 390       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 391         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
 392       }
 393     }
 394     break;
 395
 396   case 1: // Exponential filter digital Bogdan
 397     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 398       FilterSimDeConvExpD( fADCF[iCol], itarg, fNTimeBin, nexp);
 399       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 400         fADCF[iCol][iTime] = itarg[iTime];
 401       }
 402     }
 403     break;
 404
 405   case 2: // Exponential filter Marian special
 406     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 407       FilterSimDeConvExpMI( fADCF[iCol], dtarg, fNTimeBin);
 408       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 409         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
 410       }
 411     }
 412     break;
 413
 414   default:
 415     AliError(Form("Invalid filter type %d ! \n", tftype ));
 416     break;
 417   }
 418
 419   delete dtarg;
 420   delete itarg;
 421 }
 422
 423 //_____________________________________________________________________________
 424 void AliTRDmcmSim::ZSMapping()
 425 {
 426   // Zero Suppression Mapping implemented in TRAP chip
 427   //
 428   // See detail TRAP manual "Data Indication" section:
 429   // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
 430
 431   Int_t EBIS = fFeeParam->GetEBsglIndThr();       // TRAP default = 0x4  (Tis=4)
 432   Int_t EBIT = fFeeParam->GetEBsumIndThr();       // TRAP default = 0x28 (Tit=40)
 433   Int_t EBIL = fFeeParam->GetEBindLUT();          // TRAP default = 0xf0 (lookup table accept (I2,I1,I0)=(111) or (110) or (101) or (100))
 434   Int_t EBIN = fFeeParam->GetEBignoreNeighbour(); // TRAP default = 1 (no neighbor sensitivity)
 435   Int_t ep   = AliTRDfeeParam::GetPFeffectPedestal();
 436
 437   if( !CheckInitialized() ) return;
 438
 439   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
 440     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 441
 442       // Get ADC data currently in filter buffer
 443       Int_t Ap = fADCF[iadc-1][it] - ep; // previous
 444       Int_t Ac = fADCF[iadc  ][it] - ep; // current
 445       Int_t An = fADCF[iadc+1][it] - ep; // next
 446
 447       // evaluate three conditions
 448       Int_t I0 = ( Ac >=  Ap && Ac >=  An ) ? 0 : 1; // peak center detection
 449       Int_t I1 = ( Ap + Ac + An > EBIT )    ? 0 : 1; // cluster
 450       Int_t I2 = ( Ac > EBIS )              ? 0 : 1; // absolute large peak
 451
 452       Int_t I = I2 * 4 + I1 * 2 + I0;    // Bit position in lookup table
 453       Int_t D = (EBIL >> I) & 1;         // Looking up  (here D=0 means true and D=1 means false according to TRAP manual)
 454
 455       fZSM[iadc][it] &= D;
 456       if( EBIN == 0 ) {  // turn on neighboring ADCs
 457         fZSM[iadc-1][it] &= D;
 458         fZSM[iadc+1][it] &= D;
 459       }
 460     }
 461   }
 462
 463   // do 1 dim projection
 464   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 465     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 466       fZSM1Dim[iadc] &= fZSM[iadc][it];
 467     }
 468   }
 469 }
 470
 471 //_____________________________________________________________________________
 472 void AliTRDmcmSim::DumpData( char *f, char *target )
 473 {
 474   // Dump data stored (for debugging).
 475   // target should contain one or multiple of the following characters
 476   //   R   for raw data
 477   //   F   for filtered data
 478   //   Z   for zero suppression map
 479   //   S   Raw dat astream
 480   // other characters are simply ignored
 481   UInt_t tempbuf[1024];
 482
 483   if( !CheckInitialized() ) return;
 484
 485   std::ofstream of( f, std::ios::out | std::ios::app );
 486   of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
 487              fChaId, fSector, fStack, fLayer, fRobPos, fMcmPos );
 488
 489   for( int t=0 ; target[t] != 0 ; t++ ) {
 490     switch( target[t] ) {
 491     case 'R' :
 492     case 'r' :
 493       of << Form("fADCR (raw ADC data)\n");
 494       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 495         of << Form("  ADC %02d: ", iadc);
 496         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 497           of << Form("% 4d",  fADCR[iadc][it]);
 498         }
 499         of << Form("\n");
 500       }
 501       break;
 502     case 'F' :
 503     case 'f' :
 504       of << Form("fADCF (filtered ADC data)\n");
 505       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 506         of << Form("  ADC %02d: ", iadc);
 507         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 508           of << Form("% 4d",  fADCF[iadc][it]);
 509         }
 510         of << Form("\n");
 511       }
 512       break;
 513     case 'Z' :
 514     case 'z' :
 515       of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
 516       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 517         of << Form("  ADC %02d: ", iadc);
 518         if( fZSM1Dim[iadc] == 0 ) { of << " R   " ; } else { of << " .   "; } // R:read .:suppressed
 519         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 520           if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
 521         }
 522         of << Form("\n");
 523       }
 524       break;
 525     case 'S' :
 526     case 's' :
 527       Int_t s = ProduceRawStream( tempbuf, 1024 );
 528       of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
 529       of << Form("  address  data\n");
 530       for( int i = 0 ; i < s ; i++ ) {
 531         of << Form("  %04x     %08x\n", i, tempbuf[i]);
 532       }
 533     }
 534   }
 535 }
 536
 537 //_____________________________________________________________________________
 538 void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target, Int_t n, Int_t nexp)
 539 {
 540   //
 541   // Exponential filter "analog"
 542   // source will not be changed
 543   //
 544
 545   Int_t    i = 0;
 546   Int_t    k = 0;
 547   Double_t reminder[2];
 548   Double_t correction;
 549   Double_t result;
 550   Double_t rates[2];
 551   Double_t coefficients[2];
 552
 553   // Initialize (coefficient = alpha, rates = lambda)
 554   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
 555
 556   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
 557   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
 558   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
 559   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
 560
 561   coefficients[0] = c1;
 562   coefficients[1] = c2;
 563
 564   Double_t dt = 0.1;
 565   rates[0] = TMath::Exp(-dt/(r1));
 566   rates[1] = TMath::Exp(-dt/(r2));
 567
 568   // Attention: computation order is important
 569   correction = 0.0;
 570   for (k = 0; k < nexp; k++) {
 571     reminder[k] = 0.0;
 572   }
 573
 574   for (i = 0; i < n; i++) {
 575
 576     result    = ((Double_t)source[i] - correction);    // no rescaling
 577     target[i] = result;
 578
 579     for (k = 0; k < nexp; k++) {
 580       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
 581     }
 582
 583     correction = 0.0;
 584     for (k = 0; k < nexp; k++) {
 585       correction += reminder[k];
 586     }
 587   }
 588 }
 589
 590 //_____________________________________________________________________________
 591 void AliTRDmcmSim::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n, Int_t nexp)
 592 {
 593   //
 594   // Exponential filter "digital"
 595   // source will not be changed
 596   //
 597
 598   Int_t i        = 0;
 599   Int_t fAlphaL  = 0;
 600   Int_t fAlphaS  = 0;
 601   Int_t fTailPed = 0;
 602   Int_t iAlphaL  = 0;
 603   Int_t iAlphaS  = 0;
 604
 605   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
 606   // initialize (coefficient = alpha, rates = lambda)
 607
 608   Double_t dt = 0.1;
 609   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
 610   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
 611   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
 612   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
 613
 614   Int_t fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
 615   Int_t fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
 616   Int_t iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; //  9 bit paramter + fixed bits
 617   Int_t iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; //  9 bit paramter + fixed bits
 618
 619   if (nexp == 1) {
 620     fAlphaL = (Int_t) (c1 * 2048.0);
 621     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
 622   }
 623   if (nexp == 2) {
 624     fAlphaL = (Int_t) (c1 * 2048.0);
 625     fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
 626     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
 627     iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400;              // 10 bit paramter + fixed bits
 628   }
 629
 630   Double_t iAl = iAlphaL  / 2048.0;            // alpha  L: correspondence to floating point numbers
 631   Double_t iAs = iAlphaS  / 2048.0;            // alpha  S: correspondence to floating point numbers
 632   Double_t iLl = iLambdaL / 2048.0;            // lambda L: correspondence to floating point numbers
 633   Double_t iLs = iLambdaS / 2048.0;            // lambda S: correspondence to floating point numbers
 634
 635   Int_t h1;
 636   Int_t h2;
 637   Int_t rem1;
 638   Int_t rem2;
 639   Int_t correction;
 640   Int_t result;
 641   Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2;
 642
 643   Double_t xi = 1 - (iLl*iAs + iLs*iAl);             // Calculation of equilibrium values of the
 644   rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
 645   rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
 646
 647   // further initialization
 648   if ((rem1 + rem2) > 0x0FFF) {
 649     correction = 0x0FFF;
 650   }
 651   else {
 652     correction = (rem1 + rem2) & 0x0FFF;
 653   }
 654
 655   fTailPed = iFactor - correction;
 656
 657   for (i = 0; i < n; i++) {
 658
 659     result = (source[i]  - correction);
 660     if (result < 0) { // Too much undershoot
 661       result = 0;
 662     }
 663
 664     target[i] = result;
 665
 666     h1 = (rem1 + ((iAlphaL * result) >> 11));
 667     if (h1 > 0x0FFF) {
 668       h1 = 0x0FFF;
 669     }
 670     else {
 671       h1 &= 0x0FFF;
 672     }
 673
 674     h2 = (rem2 + ((iAlphaS * result) >> 11));
 675     if (h2 > 0x0FFF) {
 676       h2 = 0x0FFF;
 677     }
 678     else {
 679       h2 &= 0x0FFF;
 680     }
 681
 682     rem1 = (iLambdaL * h1 ) >> 11;
 683     rem2 = (iLambdaS * h2 ) >> 11;
 684
 685     if ((rem1 + rem2) > 0x0FFF) {
 686       correction = 0x0FFF;
 687     }
 688     else {
 689       correction = (rem1 + rem2) & 0x0FFF;
 690     }
 691
 692   }
 693 }
 694
 695 //_____________________________________________________________________________
 696 void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target, Int_t n)
 697 {
 698   //
 699   // Exponential filter (M. Ivanov)
 700   // source will not be changed
 701   //
 702
 703   Int_t i = 0;
 704   Double_t sig1[100];
 705   Double_t sig2[100];
 706   Double_t sig3[100];
 707
 708   for (i = 0; i < n; i++) {
 709     sig1[i] = (Double_t)source[i];
 710   }
 711
 712   Float_t dt      = 0.1;
 713   Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt;
 714   Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt;
 715
 716   FilterSimTailMakerSpline( sig1, sig2, lambda0, n);
 717   FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n);
 718
 719   for (i = 0; i < n; i++) {
 720     target[i] = sig3[i];
 721   }
 722
 723 }
 724
 725 //______________________________________________________________________________
 726 void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n)
 727 {
 728   //
 729   // Special filter (M. Ivanov)
 730   //
 731
 732   Int_t    i = 0;
 733   Double_t l = TMath::Exp(-lambda*0.5);
 734   Double_t in[1000];
 735   Double_t out[1000];
 736
 737   // Initialize in[] and out[] goes 0 ... 2*n+19
 738   for (i = 0; i < n*2+20; i++) {
 739     in[i] = out[i] = 0;
 740   }
 741
 742   // in[] goes 0, 1
 743   in[0] = ampin[0];
 744   in[1] = (ampin[0] + ampin[1]) * 0.5;
 745
 746   // Add charge to the end
 747   for (i = 0; i < 22; i++) {
 748     in[2*(n-1)+i] = ampin[n-1]; // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
 749   }
 750
 751   // Use arithmetic mean
 752   for (i = 1; i < n-1; i++) {
 753     in[2*i]   = ampin[i];    // in[] goes 2, 3, ... , 2*n-4, 2*n-3
 754     in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
 755   }
 756
 757   Double_t temp;
 758   out[2*n]    = in[2*n];
 759   temp        = 0;
 760   for (i = 2*n; i >= 0; i--) {
 761     out[i]    = in[i] + temp;
 762     temp      = l*(temp+in[i]);
 763   }
 764
 765   for (i = 0; i < n; i++){
 766     //ampout[i] = out[2*i+1];  // org
 767     ampout[i] = out[2*i];
 768   }
 769
 770 }
 771
 772 //______________________________________________________________________________
 773 void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout, Double_t norm, Double_t lambda, Int_t n)
 774 {
 775   //
 776   // Special filter (M. Ivanov)
 777   //
 778
 779   Int_t    i = 0;
 780
 781   Double_t l = TMath::Exp(-lambda*0.5);
 782   Double_t k = l*(1.0 - norm*lambda*0.5);
 783   Double_t in[1000];
 784   Double_t out[1000];
 785
 786   // Initialize in[] and out[] goes 0 ... 2*n+19
 787   for (i = 0; i < n*2+20; i++) {
 788     in[i] = out[i] = 0;
 789   }
 790
 791   // in[] goes 0, 1
 792   in[0] = ampin[0];
 793   in[1] = (ampin[0]+ampin[1])*0.5;
 794
 795   // Add charge to the end
 796   for (i =-2; i < 22; i++) {
 797     // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
 798     in[2*(n-1)+i] = ampin[n-1];
 799   }
 800
 801   for (i = 1; i < n-2; i++) {
 802     // in[] goes 2, 3, ... , 2*n-6, 2*n-5
 803     in[2*i]    = ampin[i];
 804     in[2*i+1]  = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
 805     //in[2*i+1]  = ((ampin[i]+ampin[i+1]))/2.0;
 806   }
 807
 808   Double_t temp;
 809   out[0] = in[0];
 810   temp   = in[0];
 811   for (i = 1; i <= 2*n; i++) {
 812     out[i] = in[i] + (k-l)*temp;
 813     temp   = in[i] +  k   *temp;
 814   }
 815
 816   for (i = 0; i < n; i++) {
 817     //ampout[i] = out[2*i+1];  // org
 818     //ampout[i] = TMath::Max(out[2*i+1],0.0);  // org
 819     ampout[i] = TMath::Max(out[2*i],0.0);
 820   }
 821 }
 822
 823 // EOF