Adding libSTEERBase to the list of libs needed to build the static DA executable
[u/mrichter/AliRoot.git] / TRD / AliTRDmcmSim.cxx
CommitLineData
dfd03fc3 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
33ClassImp(AliTRDmcmSim)
34
35//_____________________________________________________________________________
36AliTRDmcmSim::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//_____________________________________________________________________________
65AliTRDmcmSim::~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//_____________________________________________________________________________
85void 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//_____________________________________________________________________________
131void 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//_____________________________________________________________________________
151void 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//_____________________________________________________________________________
169void 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//_____________________________________________________________________________
189Int_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//_____________________________________________________________________________
200Int_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//_____________________________________________________________________________
270void 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//_____________________________________________________________________________
293void 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//_____________________________________________________________________________
309void AliTRDmcmSim::FilterGain()
310{
311 // Apply gain filter (not implemented)
312}
313
314//_____________________________________________________________________________
315void 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//_____________________________________________________________________________
363void 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//_____________________________________________________________________________
403void 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//_____________________________________________________________________________
456void 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//_____________________________________________________________________________
567void 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//______________________________________________________________________________
597void 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//______________________________________________________________________________
644void 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