1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////////
18 // TRD MCM (Multi Chip Module) simulator //
20 ///////////////////////////////////////////////////////////////////////////////
26 New release on 2007/08/17
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.
31 The following internal parameters were abolished because it is useless and
37 GetCol member was modified accordingly.
39 New member function DumpData was prepared for diagnostics.
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.
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:
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.
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).
65 In TRD definition, we have now 4 raw data versions.
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)
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.
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.
89 #include "AliTRDmcmSim.h"
90 #include "AliTRDfeeParam.h"
91 #include "AliTRDSimParam.h"
92 #include "AliTRDgeometry.h"
93 #include "AliTRDcalibDB.h"
95 ClassImp(AliTRDmcmSim)
97 //_____________________________________________________________________________
98 AliTRDmcmSim::AliTRDmcmSim() :TObject()
119 // AliTRDmcmSim default constructor
122 // By default, nothing is initialized.
123 // It is necessary to issue Init before use.
126 //_____________________________________________________________________________
127 AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m)
129 ,fInitialized(kFALSE)
149 // AliTRDmcmSim copy constructor
152 // By default, nothing is initialized.
153 // It is necessary to issue Init before use.
156 //_____________________________________________________________________________
157 AliTRDmcmSim::~AliTRDmcmSim()
160 // AliTRDmcmSim destructor
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];
178 //_____________________________________________________________________________
179 AliTRDmcmSim &AliTRDmcmSim::operator=(const AliTRDmcmSim &m)
182 // Assignment operator
186 ((AliTRDmcmSim &) m).Copy(*this);
192 //_____________________________________________________________________________
193 void AliTRDmcmSim::Copy(TObject &m) const
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;
220 //_____________________________________________________________________________
221 void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos )
224 // Initialize the class with new geometry information
225 // fADC array will be reused with filled by zero
228 fFeeParam = AliTRDfeeParam::Instance();
229 fSimParam = AliTRDSimParam::Instance();
230 fCal = AliTRDcalibDB::Instance();
231 fGeo = new AliTRDgeometry();
233 fSector = fGeo->GetSector( fChaId );
234 fStack = fGeo->GetChamber( fChaId );
235 fLayer = fGeo->GetPlane( fChaId );
238 fNADC = fFeeParam->GetNadcMcm();
239 fNTimeBin = fCal->GetNumberOfTimeBins();
240 fRow = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
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];
255 // Initialize ADC data
256 for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
257 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
260 fZSM [iadc][it] = 1; // Default unread = 1
262 fZSM1Dim[iadc] = 1; // Default unread = 1
265 fInitialized = kTRUE;
268 //_____________________________________________________________________________
269 Bool_t AliTRDmcmSim::CheckInitialized()
272 // Check whether object is initialized
275 if( ! fInitialized ) {
276 AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
281 //_____________________________________________________________________________
282 void AliTRDmcmSim::SetData( Int_t iadc, Int_t *adc )
285 // Store ADC data into array of raw data
288 if( !CheckInitialized() ) return;
290 if( iadc < 0 || iadc >= fNADC ) {
291 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
295 for( int it = 0 ; it < fNTimeBin ; it++ ) {
296 fADCR[iadc][it] = (Int_t)(adc[it]);
300 //_____________________________________________________________________________
301 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
304 // Store ADC data into array of raw data
307 if( !CheckInitialized() ) return;
309 if( iadc < 0 || iadc >= fNADC ) {
310 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
314 fADCR[iadc][it] = adc;
317 //_____________________________________________________________________________
318 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
321 // Store ADC data into array of raw data
324 if( !CheckInitialized() ) return;
326 if( iadc < 0 || iadc >= fNADC ) {
327 //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
331 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
332 fADCR[iadc][it] = fSimParam->GetADCbaseline();
336 //_____________________________________________________________________________
337 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
340 // Return column id of the pad for the given ADC channel
343 if( !CheckInitialized() ) return -1;
345 return fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
348 //_____________________________________________________________________________
349 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize )
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
359 Int_t nw = 0; // Number of written words
360 Int_t of = 0; // Number of overflowed words
361 Int_t rawVer = fFeeParam->GetRAWversion();
364 if( !CheckInitialized() ) return 0;
366 if( fFeeParam->GetRAWstoreRaw() ) {
372 // Produce MCM header
373 x = ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
384 for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
385 if( fZSM1Dim[iAdc] == 0 ) { // 0 means not suppressed
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
400 UInt_t aa=0, a1=0, a2=0, a3=0;
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;
419 if( of != 0 ) return -of; else return nw;
422 //_____________________________________________________________________________
423 void AliTRDmcmSim::Filter()
426 // Apply digital filter
429 if( !CheckInitialized() ) return;
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];
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();
444 //_____________________________________________________________________________
445 void AliTRDmcmSim::FilterPedestal()
448 // Apply pedestal filter
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
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;
462 //_____________________________________________________________________________
463 void AliTRDmcmSim::FilterGain()
466 // Apply gain filter (not implemented)
467 // Later it will be implemented because gain digital filiter will
468 // increase noise level.
473 //_____________________________________________________________________________
474 void AliTRDmcmSim::FilterTail()
477 // Apply exponential tail filter (Bogdan's version)
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();
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]);
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];
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]);
515 AliError(Form("Invalid filter type %d ! \n", tftype ));
523 //_____________________________________________________________________________
524 void AliTRDmcmSim::ZSMapping()
527 // Zero Suppression Mapping implemented in TRAP chip
529 // See detail TRAP manual "Data Indication" section:
530 // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
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();
541 if( !CheckInitialized() ) return;
543 for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
544 for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
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
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
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)
561 if( eBIN == 0 ) { // turn on neighboring ADCs
562 fZSM[iadc-1][it] &= d;
563 fZSM[iadc+1][it] &= d;
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];
578 //_____________________________________________________________________________
579 void AliTRDmcmSim::DumpData( char *f, char *target )
582 // Dump data stored (for debugging).
583 // target should contain one or multiple of the following characters
585 // F for filtered data
586 // Z for zero suppression map
588 // other characters are simply ignored
591 UInt_t tempbuf[1024];
593 if( !CheckInitialized() ) return;
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 );
599 for( int t=0 ; target[t] != 0 ; t++ ) {
600 switch( target[t] ) {
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]);
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]);
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
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]);
647 //_____________________________________________________________________________
648 void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target
649 , Int_t n, Int_t nexp)
652 // Exponential filter "analog"
653 // source will not be changed
658 Double_t reminder[2];
662 Double_t coefficients[2];
664 // Initialize (coefficient = alpha, rates = lambda)
665 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
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();
672 coefficients[0] = c1;
673 coefficients[1] = c2;
676 rates[0] = TMath::Exp(-dt/(r1));
677 rates[1] = TMath::Exp(-dt/(r2));
679 // Attention: computation order is important
681 for (k = 0; k < nexp; k++) {
685 for (i = 0; i < n; i++) {
687 result = ((Double_t)source[i] - correction); // no rescaling
690 for (k = 0; k < nexp; k++) {
691 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
695 for (k = 0; k < nexp; k++) {
696 correction += reminder[k];
701 //_____________________________________________________________________________
702 void AliTRDmcmSim::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n
706 // Exponential filter "digital"
707 // source will not be changed
717 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
718 // initialize (coefficient = alpha, rates = lambda)
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();
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
732 fAlphaL = (Int_t) (c1 * 2048.0);
733 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
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
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
753 Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2;
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));
759 // further initialization
760 if ((rem1 + rem2) > 0x0FFF) {
764 correction = (rem1 + rem2) & 0x0FFF;
767 fTailPed = iFactor - correction;
769 for (i = 0; i < n; i++) {
771 result = (source[i] - correction);
772 if (result < 0) { // Too much undershoot
778 h1 = (rem1 + ((iAlphaL * result) >> 11));
786 h2 = (rem2 + ((iAlphaS * result) >> 11));
794 rem1 = (iLambdaL * h1 ) >> 11;
795 rem2 = (iLambdaS * h2 ) >> 11;
797 if ((rem1 + rem2) > 0x0FFF) {
801 correction = (rem1 + rem2) & 0x0FFF;
808 //_____________________________________________________________________________
809 void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target
813 // Exponential filter (M. Ivanov)
814 // source will not be changed
822 for (i = 0; i < n; i++) {
823 sig1[i] = (Double_t)source[i];
827 Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt;
828 Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt;
830 FilterSimTailMakerSpline( sig1, sig2, lambda0, n);
831 FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n);
833 for (i = 0; i < n; i++) {
839 //______________________________________________________________________________
840 void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout
841 , Double_t lambda, Int_t n)
844 // Special filter (M. Ivanov)
848 Double_t l = TMath::Exp(-lambda*0.5);
852 // Initialize in[] and out[] goes 0 ... 2*n+19
853 for (i = 0; i < n*2+20; i++) {
859 in[1] = (ampin[0] + ampin[1]) * 0.5;
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
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.;
875 for (i = 2*n; i >= 0; i--) {
876 out[i] = in[i] + temp;
877 temp = l*(temp+in[i]);
880 for (i = 0; i < n; i++){
881 //ampout[i] = out[2*i+1]; // org
882 ampout[i] = out[2*i];
887 //______________________________________________________________________________
888 void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout
889 , Double_t norm, Double_t lambda
893 // Special filter (M. Ivanov)
898 Double_t l = TMath::Exp(-lambda*0.5);
899 Double_t k = l*(1.0 - norm*lambda*0.5);
903 // Initialize in[] and out[] goes 0 ... 2*n+19
904 for (i = 0; i < n*2+20; i++) {
910 in[1] = (ampin[0]+ampin[1])*0.5;
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];
918 for (i = 1; i < n-2; i++) {
919 // in[] goes 2, 3, ... , 2*n-6, 2*n-5
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;
928 for (i = 1; i <= 2*n; i++) {
929 out[i] = in[i] + (k-l)*temp;
930 temp = in[i] + k *temp;
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);