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