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 tracklet part is not yet touched.
  78 The default raw version is 2.
  79
  80                                                                  Ken Oyama
  81 */
  82
  83 #include <fstream>
  84
  85 #include <TMath.h>
  86
  87 #include "AliLog.h"
  88
  89 #include "AliTRDmcmSim.h"
  90 #include "AliTRDfeeParam.h"
  91 #include "AliTRDSimParam.h"
  92 #include "AliTRDgeometry.h"
  93 #include "AliTRDcalibDB.h"
  94
  95 ClassImp(AliTRDmcmSim)
  96
  97 //_____________________________________________________________________________
  98 AliTRDmcmSim::AliTRDmcmSim() :TObject()
  99   ,fInitialized(kFALSE)
 100   ,fChaId(-1)
 101   ,fSector(-1)
 102   ,fStack(-1)
 103   ,fLayer(-1)
 104   ,fRobPos(-1)
 105   ,fMcmPos(-1)
 106   ,fNADC(-1)
 107   ,fNTimeBin(-1)
 108   ,fRow(-1)
 109   ,fADCR(NULL)
 110   ,fADCF(NULL)
 111   ,fZSM(NULL)
 112   ,fZSM1Dim(NULL)
 113   ,fFeeParam(NULL)
 114   ,fSimParam(NULL)
 115   ,fCal(NULL)
 116   ,fGeo(NULL)
 117 {
 118   //
 119   // AliTRDmcmSim default constructor
 120   //
 121
 122   // By default, nothing is initialized.
 123   // It is necessary to issue Init before use.
 124 }
 125
 126 //_____________________________________________________________________________
 127 AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m)
 128   :TObject(m)
 129   ,fInitialized(kFALSE)
 130   ,fChaId(-1)
 131   ,fSector(-1)
 132   ,fStack(-1)
 133   ,fLayer(-1)
 134   ,fRobPos(-1)
 135   ,fMcmPos(-1)
 136   ,fNADC(-1)
 137   ,fNTimeBin(-1)
 138   ,fRow(-1)
 139   ,fADCR(NULL)
 140   ,fADCF(NULL)
 141   ,fZSM(NULL)
 142   ,fZSM1Dim(NULL)
 143   ,fFeeParam(NULL)
 144   ,fSimParam(NULL)
 145   ,fCal(NULL)
 146   ,fGeo(NULL)
 147 {
 148   //
 149   // AliTRDmcmSim copy constructor
 150   //
 151
 152   // By default, nothing is initialized.
 153   // It is necessary to issue Init before use.
 154 }
 155
 156 //_____________________________________________________________________________
 157 AliTRDmcmSim::~AliTRDmcmSim()
 158 {
 159   //
 160   // AliTRDmcmSim destructor
 161   //
 162
 163   if( fADCR != NULL ) {
 164     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 165       delete [] fADCR[iadc];
 166       delete [] fADCF[iadc];
 167       delete [] fZSM [iadc];
 168     }
 169     delete [] fADCR;
 170     delete [] fADCF;
 171     delete [] fZSM;
 172     delete [] fZSM1Dim;
 173   }
 174   delete fGeo;
 175
 176 }
 177
 178 //_____________________________________________________________________________
 179 AliTRDmcmSim &AliTRDmcmSim::operator=(const AliTRDmcmSim &m)
 180 {
 181   //
 182   // Assignment operator
 183   //
 184
 185   if (this != &m) {
 186     ((AliTRDmcmSim &) m).Copy(*this);
 187   }
 188   return *this;
 189
 190 }
 191
 192 //_____________________________________________________________________________
 193 void AliTRDmcmSim::Copy(TObject &m) const
 194 {
 195   //
 196   // Copy function
 197   //
 198
 199   ((AliTRDmcmSim &) m).fInitialized   = 0;
 200   ((AliTRDmcmSim &) m).fChaId         = 0;
 201   ((AliTRDmcmSim &) m).fSector        = 0;
 202   ((AliTRDmcmSim &) m).fStack         = 0;
 203   ((AliTRDmcmSim &) m).fLayer         = 0;
 204   ((AliTRDmcmSim &) m).fRobPos        = 0;
 205   ((AliTRDmcmSim &) m).fMcmPos        = 0;
 206   ((AliTRDmcmSim &) m).fNADC          = 0;
 207   ((AliTRDmcmSim &) m).fNTimeBin      = 0;
 208   ((AliTRDmcmSim &) m).fRow           = 0;
 209   ((AliTRDmcmSim &) m).fADCR          = 0;
 210   ((AliTRDmcmSim &) m).fADCF          = 0;
 211   ((AliTRDmcmSim &) m).fZSM           = 0;
 212   ((AliTRDmcmSim &) m).fZSM1Dim       = 0;
 213   ((AliTRDmcmSim &) m).fFeeParam      = 0;
 214   ((AliTRDmcmSim &) m).fSimParam      = 0;
 215   ((AliTRDmcmSim &) m).fCal           = 0;
 216   ((AliTRDmcmSim &) m).fGeo           = 0;
 217
 218 }
 219
 220 //_____________________________________________________________________________
 221 void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos )
 222 {
 223   //
 224   // Initialize the class with new geometry information
 225   // fADC array will be reused with filled by zero
 226   //
 227
 228   fFeeParam      = AliTRDfeeParam::Instance();
 229   fSimParam      = AliTRDSimParam::Instance();
 230   fCal           = AliTRDcalibDB::Instance();
 231   fGeo           = new AliTRDgeometry();
 232   fChaId         = chaId;
 233   fSector        = fGeo->GetSector( fChaId );
 234   fStack         = fGeo->GetChamber( fChaId );
 235   fLayer         = fGeo->GetPlane( fChaId );
 236   fRobPos        = robPos;
 237   fMcmPos        = mcmPos;
 238   fNADC          = fFeeParam->GetNadcMcm();
 239   fNTimeBin      = fCal->GetNumberOfTimeBins();
 240   fRow           = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
 241
 242   // Allocate ADC data memory if not yet done
 243   if( fADCR == NULL ) {
 244     fADCR    = new Int_t *[fNADC];
 245     fADCF    = new Int_t *[fNADC];
 246     fZSM     = new Int_t *[fNADC];
 247     fZSM1Dim = new Int_t  [fNADC];
 248     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 249       fADCR[iadc] = new Int_t[fNTimeBin];
 250       fADCF[iadc] = new Int_t[fNTimeBin];
 251       fZSM [iadc] = new Int_t[fNTimeBin];
 252     }
 253   }
 254
 255   // Initialize ADC data
 256   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 257     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 258       fADCR[iadc][it] = 0;
 259       fADCF[iadc][it] = 0;
 260       fZSM [iadc][it] = 1;   // Default unread = 1
 261     }
 262     fZSM1Dim[iadc] = 1;      // Default unread = 1
 263   }
 264
 265   fInitialized = kTRUE;
 266 }
 267
 268 //_____________________________________________________________________________
 269 Bool_t AliTRDmcmSim::CheckInitialized()
 270 {
 271   //
 272   // Check whether object is initialized
 273   //
 274
 275   if( ! fInitialized ) {
 276     AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
 277   }
 278   return fInitialized;
 279 }
 280
 281 //_____________________________________________________________________________
 282 void AliTRDmcmSim::SetData( Int_t iadc, Int_t *adc )
 283 {
 284   //
 285   // Store ADC data into array of raw data
 286   //
 287
 288   if( !CheckInitialized() ) return;
 289
 290   if( iadc < 0 || iadc >= fNADC ) {
 291     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
 292     return;
 293   }
 294
 295   for( int it = 0 ;  it < fNTimeBin ; it++ ) {
 296     fADCR[iadc][it] = (Int_t)(adc[it]);
 297   }
 298 }
 299
 300 //_____________________________________________________________________________
 301 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
 302 {
 303   //
 304   // Store ADC data into array of raw data
 305   //
 306
 307   if( !CheckInitialized() ) return;
 308
 309   if( iadc < 0 || iadc >= fNADC ) {
 310     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
 311     return;
 312   }
 313
 314   fADCR[iadc][it] = adc;
 315 }
 316
 317 //_____________________________________________________________________________
 318 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
 319 {
 320   //
 321   // Store ADC data into array of raw data
 322   //
 323
 324   if( !CheckInitialized() ) return;
 325
 326   if( iadc < 0 || iadc >= fNADC ) {
 327     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
 328     return;
 329   }
 330
 331   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 332     fADCR[iadc][it] = fSimParam->GetADCbaseline();
 333   }
 334 }
 335
 336 //_____________________________________________________________________________
 337 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
 338 {
 339   //
 340   // Return column id of the pad for the given ADC channel
 341   //
 342
 343   if( !CheckInitialized() ) return -1;
 344
 345   return fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
 346 }
 347
 348 //_____________________________________________________________________________
 349 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize )
 350 {
 351   //
 352   // Produce raw data stream from this MCM and put in buf
 353   // Returns number of words filled, or negative value
 354   // with -1 * number of overflowed words
 355   //
 356
 357   UInt_t  x;
 358   UInt_t  iEv = 0;
 359   Int_t   nw  = 0;  // Number of written words
 360   Int_t   of  = 0;  // Number of overflowed words
 361   Int_t   rawVer   = fFeeParam->GetRAWversion();
 362   Int_t **adc;
 363
 364   if( !CheckInitialized() ) return 0;
 365
 366   if( fFeeParam->GetRAWstoreRaw() ) {
 367     adc = fADCR;
 368   } else {
 369     adc = fADCF;
 370   }
 371
 372   // Produce MCM header
 373   x = ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
 374   if (nw < maxSize) {
 375     buf[nw++] = x;
 376   }
 377   else {
 378     of++;
 379   }
 380
 381   // Produce ADC mask
 382   if( rawVer >= 3 ) {
 383     x = 0;
 384     for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
 385       if( fZSM1Dim[iAdc] == 0 ) { //  0 means not suppressed
 386         x = x | (1 << iAdc);
 387       }
 388     }
 389     if (nw < maxSize) {
 390       buf[nw++] = x;
 391     }
 392     else {
 393       of++;
 394     }
 395   }
 396
 397   // Produce ADC data. 3 timebins are packed into one 32 bits word
 398   // In this version, different ADC channel will NOT share the same word
 399
 400   UInt_t aa=0, a1=0, a2=0, a3=0;
 401
 402   for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
 403     if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // suppressed
 404     aa = !(iAdc & 1) + 2;
 405     for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
 406       a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] : 0;
 407       a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
 408       a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
 409       x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
 410       if (nw < maxSize) {
 411         buf[nw++] = x;
 412       }
 413       else {
 414         of++;
 415       }
 416     }
 417   }
 418
 419   if( of != 0 ) return -of; else return nw;
 420 }
 421
 422 //_____________________________________________________________________________
 423 void AliTRDmcmSim::Filter()
 424 {
 425   //
 426   // Apply digital filter
 427   //
 428
 429   if( !CheckInitialized() ) return;
 430
 431   // Initialize filtered data array with raw data
 432   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 433     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 434       fADCF[iadc][it] = fADCR[iadc][it];
 435     }
 436   }
 437
 438   // Then apply fileters one by one to filtered data array
 439   if( fFeeParam->IsPFon() ) FilterPedestal();
 440   if( fFeeParam->IsGFon() ) FilterGain();
 441   if( fFeeParam->IsTFon() ) FilterTail();
 442 }
 443
 444 //_____________________________________________________________________________
 445 void AliTRDmcmSim::FilterPedestal()
 446 {
 447   //
 448   // Apply pedestal filter
 449   //
 450
 451   Int_t ap = fSimParam->GetADCbaseline();      // ADC instrinsic pedestal
 452   Int_t ep = fFeeParam->GetPFeffectPedestal(); // effective pedestal
 453   //Int_t tc = fFeeParam->GetPFtimeConstant();   // this makes no sense yet
 454
 455   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 456     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 457       fADCF[iadc][it] = fADCF[iadc][it] - ap + ep;
 458     }
 459   }
 460 }
 461
 462 //_____________________________________________________________________________
 463 void AliTRDmcmSim::FilterGain()
 464 {
 465   //
 466   // Apply gain filter (not implemented)
 467   // Later it will be implemented because gain digital filiter will
 468   // increase noise level.
 469   //
 470
 471 }
 472
 473 //_____________________________________________________________________________
 474 void AliTRDmcmSim::FilterTail()
 475 {
 476   //
 477   // Apply exponential tail filter (Bogdan's version)
 478   //
 479
 480   Double_t *dtarg  = new Double_t[fNTimeBin];
 481   Int_t    *itarg  = new Int_t[fNTimeBin];
 482   Int_t     nexp   = fFeeParam->GetTFnExp();
 483   Int_t     tftype = fFeeParam->GetTFtype();
 484
 485   switch( tftype ) {
 486
 487   case 0: // Exponential Filter Analog Bogdan
 488     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 489       FilterSimDeConvExpA( fADCF[iCol], dtarg, fNTimeBin, nexp);
 490       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 491         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
 492       }
 493     }
 494     break;
 495
 496   case 1: // Exponential filter digital Bogdan
 497     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 498       FilterSimDeConvExpD( fADCF[iCol], itarg, fNTimeBin, nexp);
 499       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 500         fADCF[iCol][iTime] = itarg[iTime];
 501       }
 502     }
 503     break;
 504
 505   case 2: // Exponential filter Marian special
 506     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
 507       FilterSimDeConvExpMI( fADCF[iCol], dtarg, fNTimeBin);
 508       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
 509         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
 510       }
 511     }
 512     break;
 513
 514   default:
 515     AliError(Form("Invalid filter type %d ! \n", tftype ));
 516     break;
 517   }
 518
 519   delete dtarg;
 520   delete itarg;
 521 }
 522
 523 //_____________________________________________________________________________
 524 void AliTRDmcmSim::ZSMapping()
 525 {
 526   //
 527   // Zero Suppression Mapping implemented in TRAP chip
 528   //
 529   // See detail TRAP manual "Data Indication" section:
 530   // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
 531   //
 532
 533   Int_t eBIS = fFeeParam->GetEBsglIndThr();       // TRAP default = 0x4  (Tis=4)
 534   Int_t eBIT = fFeeParam->GetEBsumIndThr();       // TRAP default = 0x28 (Tit=40)
 535   Int_t eBIL = fFeeParam->GetEBindLUT();          // TRAP default = 0xf0
 536                                                   // (lookup table accept (I2,I1,I0)=(111)
 537                                                   // or (110) or (101) or (100))
 538   Int_t eBIN = fFeeParam->GetEBignoreNeighbour(); // TRAP default = 1 (no neighbor sensitivity)
 539   Int_t ep   = AliTRDfeeParam::GetPFeffectPedestal();
 540
 541   if( !CheckInitialized() ) return;
 542
 543   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
 544     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 545
 546       // Get ADC data currently in filter buffer
 547       Int_t ap = fADCF[iadc-1][it] - ep; // previous
 548       Int_t ac = fADCF[iadc  ][it] - ep; // current
 549       Int_t an = fADCF[iadc+1][it] - ep; // next
 550
 551       // evaluate three conditions
 552       Int_t i0 = ( ac >=  ap && ac >=  an ) ? 0 : 1; // peak center detection
 553       Int_t i1 = ( ap + ac + an > eBIT )    ? 0 : 1; // cluster
 554       Int_t i2 = ( ac > eBIS )              ? 0 : 1; // absolute large peak
 555
 556       Int_t i = i2 * 4 + i1 * 2 + i0;    // Bit position in lookup table
 557       Int_t d = (eBIL >> i) & 1;         // Looking up  (here d=0 means true
 558                                          // and d=1 means false according to TRAP manual)
 559
 560       fZSM[iadc][it] &= d;
 561       if( eBIN == 0 ) {  // turn on neighboring ADCs
 562         fZSM[iadc-1][it] &= d;
 563         fZSM[iadc+1][it] &= d;
 564       }
 565
 566     }
 567   }
 568
 569   // do 1 dim projection
 570   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 571     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 572       fZSM1Dim[iadc] &= fZSM[iadc][it];
 573     }
 574   }
 575
 576 }
 577
 578 //_____________________________________________________________________________
 579 void AliTRDmcmSim::DumpData( char *f, char *target )
 580 {
 581   //
 582   // Dump data stored (for debugging).
 583   // target should contain one or multiple of the following characters
 584   //   R   for raw data
 585   //   F   for filtered data
 586   //   Z   for zero suppression map
 587   //   S   Raw dat astream
 588   // other characters are simply ignored
 589   //
 590
 591   UInt_t tempbuf[1024];
 592
 593   if( !CheckInitialized() ) return;
 594
 595   std::ofstream of( f, std::ios::out | std::ios::app );
 596   of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
 597              fChaId, fSector, fStack, fLayer, fRobPos, fMcmPos );
 598
 599   for( int t=0 ; target[t] != 0 ; t++ ) {
 600     switch( target[t] ) {
 601     case 'R' :
 602     case 'r' :
 603       of << Form("fADCR (raw ADC data)\n");
 604       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 605         of << Form("  ADC %02d: ", iadc);
 606         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 607           of << Form("% 4d",  fADCR[iadc][it]);
 608         }
 609         of << Form("\n");
 610       }
 611       break;
 612     case 'F' :
 613     case 'f' :
 614       of << Form("fADCF (filtered ADC data)\n");
 615       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 616         of << Form("  ADC %02d: ", iadc);
 617         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 618           of << Form("% 4d",  fADCF[iadc][it]);
 619         }
 620         of << Form("\n");
 621       }
 622       break;
 623     case 'Z' :
 624     case 'z' :
 625       of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
 626       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
 627         of << Form("  ADC %02d: ", iadc);
 628         if( fZSM1Dim[iadc] == 0 ) { of << " R   " ; } else { of << " .   "; } // R:read .:suppressed
 629         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
 630           if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
 631         }
 632         of << Form("\n");
 633       }
 634       break;
 635     case 'S' :
 636     case 's' :
 637       Int_t s = ProduceRawStream( tempbuf, 1024 );
 638       of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
 639       of << Form("  address  data\n");
 640       for( int i = 0 ; i < s ; i++ ) {
 641         of << Form("  %04x     %08x\n", i, tempbuf[i]);
 642       }
 643     }
 644   }
 645 }
 646
 647 //_____________________________________________________________________________
 648 void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target
 649                                      , Int_t n, Int_t nexp)
 650 {
 651   //
 652   // Exponential filter "analog"
 653   // source will not be changed
 654   //
 655
 656   Int_t    i = 0;
 657   Int_t    k = 0;
 658   Double_t reminder[2];
 659   Double_t correction;
 660   Double_t result;
 661   Double_t rates[2];
 662   Double_t coefficients[2];
 663
 664   // Initialize (coefficient = alpha, rates = lambda)
 665   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
 666
 667   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
 668   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
 669   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
 670   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
 671
 672   coefficients[0] = c1;
 673   coefficients[1] = c2;
 674
 675   Double_t dt = 0.1;
 676   rates[0] = TMath::Exp(-dt/(r1));
 677   rates[1] = TMath::Exp(-dt/(r2));
 678
 679   // Attention: computation order is important
 680   correction = 0.0;
 681   for (k = 0; k < nexp; k++) {
 682     reminder[k] = 0.0;
 683   }
 684
 685   for (i = 0; i < n; i++) {
 686
 687     result    = ((Double_t)source[i] - correction);    // no rescaling
 688     target[i] = result;
 689
 690     for (k = 0; k < nexp; k++) {
 691       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
 692     }
 693
 694     correction = 0.0;
 695     for (k = 0; k < nexp; k++) {
 696       correction += reminder[k];
 697     }
 698   }
 699 }
 700
 701 //_____________________________________________________________________________
 702 void AliTRDmcmSim::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n
 703                                      , Int_t nexp)
 704 {
 705   //
 706   // Exponential filter "digital"
 707   // source will not be changed
 708   //
 709
 710   Int_t i        = 0;
 711   Int_t fAlphaL  = 0;
 712   Int_t fAlphaS  = 0;
 713   Int_t fTailPed = 0;
 714   Int_t iAlphaL  = 0;
 715   Int_t iAlphaS  = 0;
 716
 717   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
 718   // initialize (coefficient = alpha, rates = lambda)
 719
 720   Double_t dt = 0.1;
 721   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
 722   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
 723   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
 724   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
 725
 726   Int_t fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
 727   Int_t fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
 728   Int_t iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; //  9 bit paramter + fixed bits
 729   Int_t iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; //  9 bit paramter + fixed bits
 730
 731   if (nexp == 1) {
 732     fAlphaL = (Int_t) (c1 * 2048.0);
 733     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
 734   }
 735   if (nexp == 2) {
 736     fAlphaL = (Int_t) (c1 * 2048.0);
 737     fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
 738     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
 739     iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400;              // 10 bit paramter + fixed bits
 740   }
 741
 742   Double_t iAl = iAlphaL  / 2048.0;            // alpha  L: correspondence to floating point numbers
 743   Double_t iAs = iAlphaS  / 2048.0;            // alpha  S: correspondence to floating point numbers
 744   Double_t iLl = iLambdaL / 2048.0;            // lambda L: correspondence to floating point numbers
 745   Double_t iLs = iLambdaS / 2048.0;            // lambda S: correspondence to floating point numbers
 746
 747   Int_t h1;
 748   Int_t h2;
 749   Int_t rem1;
 750   Int_t rem2;
 751   Int_t correction;
 752   Int_t result;
 753   Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2;
 754
 755   Double_t xi = 1 - (iLl*iAs + iLs*iAl);             // Calculation of equilibrium values of the
 756   rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
 757   rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
 758
 759   // further initialization
 760   if ((rem1 + rem2) > 0x0FFF) {
 761     correction = 0x0FFF;
 762   }
 763   else {
 764     correction = (rem1 + rem2) & 0x0FFF;
 765   }
 766
 767   fTailPed = iFactor - correction;
 768
 769   for (i = 0; i < n; i++) {
 770
 771     result = (source[i]  - correction);
 772     if (result < 0) { // Too much undershoot
 773       result = 0;
 774     }
 775
 776     target[i] = result;
 777
 778     h1 = (rem1 + ((iAlphaL * result) >> 11));
 779     if (h1 > 0x0FFF) {
 780       h1 = 0x0FFF;
 781     }
 782     else {
 783       h1 &= 0x0FFF;
 784     }
 785
 786     h2 = (rem2 + ((iAlphaS * result) >> 11));
 787     if (h2 > 0x0FFF) {
 788       h2 = 0x0FFF;
 789     }
 790     else {
 791       h2 &= 0x0FFF;
 792     }
 793
 794     rem1 = (iLambdaL * h1 ) >> 11;
 795     rem2 = (iLambdaS * h2 ) >> 11;
 796
 797     if ((rem1 + rem2) > 0x0FFF) {
 798       correction = 0x0FFF;
 799     }
 800     else {
 801       correction = (rem1 + rem2) & 0x0FFF;
 802     }
 803
 804   }
 805
 806 }
 807
 808 //_____________________________________________________________________________
 809 void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target
 810                                       , Int_t n)
 811 {
 812   //
 813   // Exponential filter (M. Ivanov)
 814   // source will not be changed
 815   //
 816
 817   Int_t i = 0;
 818   Double_t sig1[100];
 819   Double_t sig2[100];
 820   Double_t sig3[100];
 821
 822   for (i = 0; i < n; i++) {
 823     sig1[i] = (Double_t)source[i];
 824   }
 825
 826   Float_t dt      = 0.1;
 827   Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt;
 828   Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt;
 829
 830   FilterSimTailMakerSpline( sig1, sig2, lambda0, n);
 831   FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n);
 832
 833   for (i = 0; i < n; i++) {
 834     target[i] = sig3[i];
 835   }
 836
 837 }
 838
 839 //______________________________________________________________________________
 840 void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout
 841                                           , Double_t lambda, Int_t n)
 842 {
 843   //
 844   // Special filter (M. Ivanov)
 845   //
 846
 847   Int_t    i = 0;
 848   Double_t l = TMath::Exp(-lambda*0.5);
 849   Double_t in[1000];
 850   Double_t out[1000];
 851
 852   // Initialize in[] and out[] goes 0 ... 2*n+19
 853   for (i = 0; i < n*2+20; i++) {
 854     in[i] = out[i] = 0;
 855   }
 856
 857   // in[] goes 0, 1
 858   in[0] = ampin[0];
 859   in[1] = (ampin[0] + ampin[1]) * 0.5;
 860
 861   // Add charge to the end
 862   for (i = 0; i < 22; i++) {
 863     in[2*(n-1)+i] = ampin[n-1]; // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
 864   }
 865
 866   // Use arithmetic mean
 867   for (i = 1; i < n-1; i++) {
 868     in[2*i]   = ampin[i];    // in[] goes 2, 3, ... , 2*n-4, 2*n-3
 869     in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
 870   }
 871
 872   Double_t temp;
 873   out[2*n]    = in[2*n];
 874   temp        = 0;
 875   for (i = 2*n; i >= 0; i--) {
 876     out[i]    = in[i] + temp;
 877     temp      = l*(temp+in[i]);
 878   }
 879
 880   for (i = 0; i < n; i++){
 881     //ampout[i] = out[2*i+1];  // org
 882     ampout[i] = out[2*i];
 883   }
 884
 885 }
 886
 887 //______________________________________________________________________________
 888 void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout
 889                                             , Double_t norm, Double_t lambda
 890                                             , Int_t n)
 891 {
 892   //
 893   // Special filter (M. Ivanov)
 894   //
 895
 896   Int_t    i = 0;
 897
 898   Double_t l = TMath::Exp(-lambda*0.5);
 899   Double_t k = l*(1.0 - norm*lambda*0.5);
 900   Double_t in[1000];
 901   Double_t out[1000];
 902
 903   // Initialize in[] and out[] goes 0 ... 2*n+19
 904   for (i = 0; i < n*2+20; i++) {
 905     in[i] = out[i] = 0;
 906   }
 907
 908   // in[] goes 0, 1
 909   in[0] = ampin[0];
 910   in[1] = (ampin[0]+ampin[1])*0.5;
 911
 912   // Add charge to the end
 913   for (i =-2; i < 22; i++) {
 914     // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
 915     in[2*(n-1)+i] = ampin[n-1];
 916   }
 917
 918   for (i = 1; i < n-2; i++) {
 919     // in[] goes 2, 3, ... , 2*n-6, 2*n-5
 920     in[2*i]    = ampin[i];
 921     in[2*i+1]  = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
 922     //in[2*i+1]  = ((ampin[i]+ampin[i+1]))/2.0;
 923   }
 924
 925   Double_t temp;
 926   out[0] = in[0];
 927   temp   = in[0];
 928   for (i = 1; i <= 2*n; i++) {
 929     out[i] = in[i] + (k-l)*temp;
 930     temp   = in[i] +  k   *temp;
 931   }
 932
 933   for (i = 0; i < n; i++) {
 934     //ampout[i] = out[2*i+1];  // org
 935     //ampout[i] = TMath::Max(out[2*i+1],0.0);  // org
 936     ampout[i] = TMath::Max(out[2*i],0.0);
 937   }
 938 }
 939