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