AliTRDmcmSim:
[u/mrichter/AliRoot.git] / TRD / AliTRDmcmSim.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //  TRD MCM (Multi Chip Module) simulator                                    //
21 //  which simulates the TRAP processing after the AD-conversion.             //
22 //  The relevant parameters (i.e. configuration settings of the TRAP)        //
23 //  are taken from AliTRDtrapConfig.                                         //
24 //                                                                           //
25 ///////////////////////////////////////////////////////////////////////////////
26
27 #include <iostream>
28 #include <iomanip>
29
30 #include "TCanvas.h"
31 #include "TH1F.h"
32 #include "TH2F.h"
33 #include "TGraph.h"
34 #include "TLine.h"
35 #include "TRandom.h"
36 #include "TClonesArray.h"
37
38 #include "AliLog.h"
39 #include "AliRunLoader.h"
40 #include "AliLoader.h"
41
42 #include "AliTRDfeeParam.h"
43 #include "AliTRDtrapConfig.h"
44 #include "AliTRDdigitsManager.h"
45 #include "AliTRDarrayADC.h"
46 #include "AliTRDarrayDictionary.h"
47 #include "AliTRDtrackletMCM.h"
48 #include "AliTRDmcmSim.h"
49
50 ClassImp(AliTRDmcmSim)
51
52 Bool_t AliTRDmcmSim::fgApplyCut = kTRUE;
53 Int_t  AliTRDmcmSim::fgAddBaseline = 0;
54
55 const Int_t AliTRDmcmSim::fgkFormatIndex = std::ios_base::xalloc(); 
56
57 const Int_t AliTRDmcmSim::fgkNADC = AliTRDfeeParam::GetNadcMcm();
58 const UShort_t AliTRDmcmSim::fgkFPshifts[4] = {11, 14, 17, 21}; 
59
60
61 AliTRDmcmSim::AliTRDmcmSim() : 
62   TObject(),
63   fInitialized(kFALSE),
64   fDetector(-1),
65   fRobPos(-1),
66   fMcmPos(-1),
67   fRow (-1),
68   fNTimeBin(-1),
69   fADCR(NULL),
70   fADCF(NULL),
71   fMCMT(NULL),
72   fTrackletArray(NULL),
73   fZSMap(NULL),
74   fFeeParam(NULL),
75   fTrapConfig(NULL),
76   fDigitsManager(NULL),
77   fPedAcc(NULL),
78   fGainCounterA(NULL),
79   fGainCounterB(NULL),
80   fTailAmplLong(NULL),
81   fTailAmplShort(NULL),
82   fNHits(0),
83   fFitReg(NULL)
84 {
85   //
86   // AliTRDmcmSim default constructor
87   // By default, nothing is initialized.
88   // It is necessary to issue Init before use.
89 }
90
91 AliTRDmcmSim::~AliTRDmcmSim() 
92 {
93   //
94   // AliTRDmcmSim destructor
95   //
96
97   if(fInitialized) {
98     for( Int_t iAdc = 0 ; iAdc < fgkNADC; iAdc++ ) {
99       delete [] fADCR[iAdc];
100       delete [] fADCF[iAdc];
101     }
102     delete [] fADCR;
103     delete [] fADCF;
104     delete [] fZSMap;
105     delete [] fMCMT;
106  
107     delete [] fPedAcc;
108     delete [] fGainCounterA;
109     delete [] fGainCounterB;
110     delete [] fTailAmplLong;
111     delete [] fTailAmplShort;
112     delete [] fFitReg;
113     
114     fTrackletArray->Delete();
115     delete fTrackletArray;
116   }
117 }
118
119 void AliTRDmcmSim::Init( Int_t det, Int_t robPos, Int_t mcmPos, Bool_t /* newEvent */ ) 
120 {
121   //
122   // Initialize the class with new MCM position information
123   // memory is allocated in the first initialization
124   //
125    
126   if (!fInitialized) {
127     fFeeParam      = AliTRDfeeParam::Instance();
128     fTrapConfig    = AliTRDtrapConfig::Instance();
129   }
130
131   fDetector      = det;
132   fRobPos        = robPos;
133   fMcmPos        = mcmPos;
134   fNTimeBin      = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kC13CPUA);
135   fRow           = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
136   
137   if (!fInitialized) {
138     fADCR    = new Int_t *[fgkNADC];
139     fADCF    = new Int_t *[fgkNADC];
140     fZSMap   = new Int_t  [fgkNADC];
141     fGainCounterA = new UInt_t[fgkNADC];
142     fGainCounterB = new UInt_t[fgkNADC];
143     for( Int_t iAdc = 0 ; iAdc < fgkNADC; iAdc++ ) {
144       fADCR[iAdc] = new Int_t[fNTimeBin];
145       fADCF[iAdc] = new Int_t[fNTimeBin];
146     }
147     
148     // filter registers
149     fPedAcc = new UInt_t[fgkNADC]; // accumulator for pedestal filter
150     fTailAmplLong = new UShort_t[fgkNADC];
151     fTailAmplShort = new UShort_t[fgkNADC];
152     
153     // tracklet calculation
154     fFitReg = new FitReg_t[fgkNADC]; 
155     fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fgkMaxTracklets);
156     
157     fMCMT = new UInt_t[fgkMaxTracklets];
158   }
159
160   fInitialized = kTRUE;
161
162   Reset();
163 }
164
165 void AliTRDmcmSim::Reset()
166 {
167   // Resets the data values and internal filter registers
168   // by re-initialising them
169
170   if( !CheckInitialized() ) 
171     return;
172
173   for( Int_t iAdc = 0 ; iAdc < fgkNADC; iAdc++ ) {
174     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
175       fADCR[iAdc][it] = 0;
176       fADCF[iAdc][it] = 0;
177     }
178     fZSMap[iAdc] = -1;      // Default unread, low active bit mask
179     fGainCounterA[iAdc] = 0;
180     fGainCounterB[iAdc] = 0;
181   }
182   
183   for(Int_t i = 0; i < fgkMaxTracklets; i++) {
184     fMCMT[i] = 0;
185   }
186
187   for (Int_t iDict = 0; iDict < 3; iDict++)
188     fDict[iDict] = 0x0;
189   
190   FilterPedestalInit();
191   FilterGainInit();
192   FilterTailInit();
193 }
194
195 void AliTRDmcmSim::SetNTimebins(Int_t ntimebins) 
196 {
197   // Reallocate memory if a change in the number of timebins 
198   // is needed (should not be the case for real data)
199
200   if( !CheckInitialized() ) 
201     return;
202
203   fNTimeBin = ntimebins;
204   for( Int_t iAdc = 0 ; iAdc < fgkNADC; iAdc++ ) {
205     delete fADCR[iAdc];
206     delete fADCF[iAdc];
207     fADCR[iAdc] = new Int_t[fNTimeBin];
208     fADCF[iAdc] = new Int_t[fNTimeBin];
209   }
210 }
211
212 Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob, Int_t mcm) 
213 {
214   // loads the ADC data as obtained from the digitsManager for the specified MCM.
215   // This method is meant for rare execution, e.g. in the visualization. When called 
216   // frequently use SetData(...) instead. 
217
218   Init(det, rob, mcm);
219
220   if (!runloader) {
221     AliError("No Runloader given");
222     return kFALSE;
223   }
224
225   AliLoader *trdLoader = runloader->GetLoader("TRDLoader");
226   if (!trdLoader) {
227     AliError("Could not get TRDLoader");
228     return kFALSE;
229   }
230
231   Bool_t retval = kTRUE;
232   trdLoader->LoadDigits();
233   fDigitsManager = 0x0;
234   AliTRDdigitsManager *digMgr = new AliTRDdigitsManager();
235   digMgr->SetSDigits(0);
236   digMgr->CreateArrays();
237   digMgr->ReadDigits(trdLoader->TreeD());
238   AliTRDarrayADC *digits = (AliTRDarrayADC*) digMgr->GetDigits(det);
239   if (digits->HasData()) {
240     digits->Expand();
241
242     if (fNTimeBin != digits->GetNtime()) {
243       AliWarning(Form("Changing no. of timebins from %i to %i", fNTimeBin, digits->GetNtime()));
244       SetNTimebins(digits->GetNtime());
245     }
246
247     SetData(digits);
248   }
249   else 
250     retval = kFALSE;
251   
252   delete digMgr;
253   
254   return retval;
255 }
256
257 void AliTRDmcmSim::NoiseTest(Int_t nsamples, Int_t mean, Int_t sigma, Int_t inputGain, Int_t inputTail)
258 {
259   // This function can be used to test the filters. 
260   // It feeds nsamples of ADC values with a gaussian distribution specified by mean and sigma.
261   // The filter chain implemented here consists of:
262   // Pedestal -> Gain -> Tail
263   // With inputGain and inputTail the input to the gain and tail filter, respectively, 
264   // can be chosen where 
265   // 0: noise input
266   // 1: pedestal output
267   // 2: gain output
268   // The input has to be chosen from a stage before. 
269   // The filter behaviour is controlled by the TRAP parameters from AliTRDtrapConfig in the 
270   // same way as in normal simulation.
271   // The functions produces four histograms with the values at the different stages.
272
273   if( !CheckInitialized() ) 
274     return;
275
276   TString nameInputGain;
277   TString nameInputTail; 
278
279   switch (inputGain) {
280       case 0:
281         nameInputGain = "Noise";
282         break;
283
284       case 1:
285         nameInputGain = "Pedestal";
286         break;
287
288       default:
289         AliError("Undefined input to tail cancellation filter");
290         return;
291   }
292
293   switch (inputTail) {
294       case 0:
295         nameInputTail = "Noise";
296         break;
297
298       case 1:
299         nameInputTail = "Pedestal";
300         break;
301
302       case 2:
303         nameInputTail = "Gain";
304         break;
305
306       default:
307         AliError("Undefined input to tail cancellation filter");
308         return;
309   }
310
311   TH1F *h   = new TH1F("noise", "Gaussian Noise;sample;ADC count",
312                        nsamples, 0, nsamples);
313   TH1F *hfp = new TH1F("ped", "Noise #rightarrow Pedestal filter;sample;ADC count", nsamples, 0, nsamples);
314   TH1F *hfg = new TH1F("gain", 
315                        (nameInputGain + "#rightarrow Gain;sample;ADC count").Data(), 
316                        nsamples, 0, nsamples);
317   TH1F *hft = new TH1F("tail", 
318                        (nameInputTail + "#rightarrow Tail;sample;ADC count").Data(), 
319                        nsamples, 0, nsamples);
320   h->SetStats(kFALSE);
321   hfp->SetStats(kFALSE);
322   hfg->SetStats(kFALSE);
323   hft->SetStats(kFALSE);
324   
325   Int_t value;  // ADC count with noise (10 bit)
326   Int_t valuep; // pedestal filter output (12 bit)
327   Int_t valueg; // gain filter output (12 bit)
328   Int_t valuet; // tail filter value (12 bit)
329   
330   for (Int_t i = 0; i < nsamples; i++) {
331     value = (Int_t) gRandom->Gaus(mean, sigma);  // generate noise with gaussian distribution 
332     h->SetBinContent(i, value);
333
334     valuep = FilterPedestalNextSample(1, 0, ((Int_t) value) << 2);
335     
336     if (inputGain == 0)
337       valueg = FilterGainNextSample(1, ((Int_t) value) << 2);
338     else 
339       valueg = FilterGainNextSample(1, valuep); 
340     
341     if (inputTail == 0)
342       valuet = FilterTailNextSample(1, ((Int_t) value) << 2);
343     else if (inputTail == 1)
344       valuet = FilterTailNextSample(1, valuep); 
345     else
346       valuet = FilterTailNextSample(1, valueg); 
347
348     hfp->SetBinContent(i, valuep >> 2);
349     hfg->SetBinContent(i, valueg >> 2);
350     hft->SetBinContent(i, valuet >> 2);
351   }
352
353   TCanvas *c = new TCanvas; 
354   c->Divide(2,2);
355   c->cd(1);
356   h->Draw();
357   c->cd(2);
358   hfp->Draw();
359   c->cd(3);
360   hfg->Draw();
361   c->cd(4);
362   hft->Draw();
363 }
364
365 Bool_t AliTRDmcmSim::CheckInitialized() const
366 {
367   //
368   // Check whether object is initialized
369   //
370
371   if( ! fInitialized ) 
372     AliError(Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
373
374   return fInitialized;
375 }
376
377 void AliTRDmcmSim::Print(Option_t* const option) const
378 {
379   // Prints the data stored and/or calculated for this MCM.
380   // The output is controlled by option which can be a sequence of any of 
381   // the following characters:
382   // R - prints raw ADC data
383   // F - prints filtered data 
384   // H - prints detected hits
385   // T - prints found tracklets
386   // The later stages are only meaningful after the corresponding calculations 
387   // have been performed.
388
389   if ( !CheckInitialized() ) 
390     return;
391
392   printf("MCM %i on ROB %i in detector %i\n", fMcmPos, fRobPos, fDetector);
393
394   TString opt = option;
395   if (opt.Contains("R") || opt.Contains("F")) {
396     std::cout << *this;
397   }
398
399   if (opt.Contains("H")) {
400     printf("Found %i hits:\n", fNHits);
401     for (Int_t iHit = 0; iHit < fNHits; iHit++) {
402       printf("Hit %3i in timebin %2i, ADC %2i has charge %3i and position %3i\n",
403              iHit,  fHits[iHit].fTimebin, fHits[iHit].fChannel, fHits[iHit].fQtot, fHits[iHit].fYpos);
404     }
405   }
406
407   if (opt.Contains("T")) {
408     printf("Tracklets:\n");
409     for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntriesFast(); iTrkl++) {
410       printf("tracklet %i: 0x%08x\n", iTrkl, ((AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl])->GetTrackletWord());
411     }
412   }
413 }
414
415 void AliTRDmcmSim::Draw(Option_t* const option) 
416 {
417   // Plots the data stored in a 2-dim. timebin vs. ADC channel plot.
418   // The option selects what data is plotted and can be a sequence of 
419   // the following characters:
420   // R - plot raw data (default)
421   // F - plot filtered data (meaningless if R is specified)
422   // In addition to the ADC values:
423   // H - plot hits 
424   // T - plot tracklets
425
426   if( !CheckInitialized() ) 
427     return;
428
429   TString opt = option;
430
431   TH2F *hist = new TH2F("mcmdata", Form("Data of MCM %i on ROB %i in detector %i", \
432                                         fMcmPos, fRobPos, fDetector), \
433                         fgkNADC, -0.5, fgkNADC-.5, fNTimeBin, -.5, fNTimeBin-.5);
434   hist->GetXaxis()->SetTitle("ADC Channel");
435   hist->GetYaxis()->SetTitle("Timebin");
436   hist->SetStats(kFALSE);
437
438   if (opt.Contains("R")) {
439     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
440       for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
441         hist->SetBinContent(iAdc+1, iTimeBin+1, fADCR[iAdc][iTimeBin] >> fgkAddDigits);
442       }
443     }
444   }
445   else {
446     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
447       for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
448         hist->SetBinContent(iAdc+1, iTimeBin+1, fADCF[iAdc][iTimeBin] >> fgkAddDigits);
449       }
450     }
451   }
452   hist->Draw("colz");
453
454   if (opt.Contains("H")) {
455     TGraph *grHits = new TGraph();
456     for (Int_t iHit = 0; iHit < fNHits; iHit++) {
457       grHits->SetPoint(iHit, 
458                        fHits[iHit].fChannel + 1 + fHits[iHit].fYpos/256., 
459                        fHits[iHit].fTimebin);
460     }
461     grHits->Draw("*");
462   }
463
464   if (opt.Contains("T")) {
465     TLine *trklLines = new TLine[4];
466     for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntries(); iTrkl++) {
467       AliTRDtrackletMCM *trkl = (AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl];
468       Float_t padWidth = 0.635 + 0.03 * (fDetector % 6);
469       Float_t offset   = padWidth/256. * ((((((fRobPos & 0x1) << 2) + (fMcmPos & 0x3)) * 18) << 8) - ((18*4*2 - 18*2 - 3) << 7)); // revert adding offset in FitTracklet
470       Int_t   ndrift   = fTrapConfig->GetDmem(0xc025, fDetector, fRobPos, fMcmPos) >> 5; 
471       Float_t slope    = trkl->GetdY() * 140e-4 / ndrift; 
472
473       Int_t t0 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS);
474       Int_t t1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE);
475
476       trklLines[iTrkl].SetX1((offset - (trkl->GetY() - slope * t0)) / padWidth); // ??? sign?
477       trklLines[iTrkl].SetY1(t0);
478       trklLines[iTrkl].SetX2((offset - (trkl->GetY() - slope * t1)) / padWidth); // ??? sign?
479       trklLines[iTrkl].SetY2(t1);
480       trklLines[iTrkl].SetLineColor(2);
481       trklLines[iTrkl].SetLineWidth(2);
482       printf("Tracklet %i: y = %f, dy = %f, offset = %f\n", iTrkl, trkl->GetY(), (trkl->GetdY() * 140e-4), offset);
483       trklLines[iTrkl].Draw();
484     }
485   }
486 }
487
488 void AliTRDmcmSim::SetData( Int_t adc, Int_t* const data )
489 {
490   //
491   // Store ADC data into array of raw data
492   //
493
494   if( !CheckInitialized() ) return;
495
496   if( adc < 0 || adc >= fgkNADC ) {
497     AliError(Form ("Error: ADC %i is out of range (0 .. %d).", adc, fgkNADC-1));
498     return;
499   }
500
501   for( Int_t it = 0 ;  it < fNTimeBin ; it++ ) {
502     fADCR[adc][it] = (Int_t) (data[it]) << fgkAddDigits;
503     fADCF[adc][it] = (Int_t) (data[it]) << fgkAddDigits;
504   }
505 }
506
507 void AliTRDmcmSim::SetData( Int_t adc, Int_t it, Int_t data )
508 {
509   //
510   // Store ADC data into array of raw data
511   //
512
513   if( !CheckInitialized() ) return;
514
515   if( adc < 0 || adc >= fgkNADC ) {
516     AliError(Form ("Error: ADC %i is out of range (0 .. %d).", adc, fgkNADC-1));
517     return;
518   }
519
520   fADCR[adc][it] = data << fgkAddDigits;
521   fADCF[adc][it] = data << fgkAddDigits;
522 }
523
524 void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager *digitsManager)
525 {
526   // Set the ADC data from an AliTRDarrayADC
527
528   if( !CheckInitialized() ) 
529     return;
530
531   fDigitsManager = digitsManager;
532   if (fDigitsManager) {
533     for (Int_t iDict = 0; iDict < 3; iDict++) {
534       AliTRDarrayDictionary *newDict = (AliTRDarrayDictionary*) fDigitsManager->GetDictionary(fDetector, iDict);
535       if (fDict[iDict] != 0x0 && newDict != 0x0) {
536         
537         if (fDict[iDict] == newDict)
538           continue;
539
540         fDict[iDict] = newDict;
541         
542         if (fDict[iDict]->GetDim() == 0) {
543           AliError(Form("Dictionary %i of det. %i has dim. 0", fDetector, iDict));
544           continue;
545         }
546         fDict[iDict]->Expand(); 
547       }
548       else {
549         fDict[iDict] = newDict;
550         if (fDict[iDict])
551           fDict[iDict]->Expand();
552       }
553     }
554   }
555
556   if (fNTimeBin != adcArray->GetNtime())
557     SetNTimebins(adcArray->GetNtime());
558   
559   Int_t offset = (fMcmPos % 4 + 1) * 21 + (fRobPos % 2) * 84 - 1;
560
561   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
562     for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
563       Int_t value = adcArray->GetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin);
564       if (value < 0 || (offset - iAdc < 1) || (offset - iAdc > 165)) {
565         fADCR[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP) + (fgAddBaseline << fgkAddDigits);
566         fADCF[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
567       }
568       else {
569         fZSMap[iAdc] = 0;
570         fADCR[iAdc][iTimeBin] = (value << fgkAddDigits) + (fgAddBaseline << fgkAddDigits);
571         fADCF[iAdc][iTimeBin] = (value << fgkAddDigits) + (fgAddBaseline << fgkAddDigits);
572       }
573     }
574   }
575 }
576
577 void AliTRDmcmSim::SetDataPedestal( Int_t adc )
578 {
579   //
580   // Store ADC data into array of raw data
581   //
582
583   if( !CheckInitialized() ) 
584     return;
585
586   if( adc < 0 || adc >= fgkNADC ) {
587     return;
588   }
589
590   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
591     fADCR[adc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP) + (fgAddBaseline << fgkAddDigits);
592     fADCF[adc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
593   }
594 }
595
596 Bool_t AliTRDmcmSim::GetHit(Int_t index, Int_t &channel, Int_t &timebin, Int_t &qtot, Int_t &ypos, Float_t &y, Int_t &label) const
597 {
598   // retrieve the MC hit information (not available in TRAP hardware)
599
600   if (index < 0 || index >= fNHits)
601     return kFALSE;
602   
603   channel = fHits[index].fChannel;
604   timebin = fHits[index].fTimebin;
605   qtot    = fHits[index].fQtot;
606   ypos    = fHits[index].fYpos;
607   y       = (Float_t) ((((((fRobPos & 0x1) << 2) + (fMcmPos & 0x3)) * 18) << 8) - ((18*4*2 - 18*2 - 1) << 7) -
608                         (channel << 8) - ypos) 
609     * (0.635 + 0.03 * (fDetector % 6))
610     / 256.0;
611   label   = fHits[index].fLabel;
612
613   return kTRUE;
614 }
615
616 Int_t AliTRDmcmSim::GetCol( Int_t adc )
617 {
618   //
619   // Return column id of the pad for the given ADC channel
620   //
621
622   if( !CheckInitialized() ) 
623     return -1;
624
625   Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adc);
626   if (col < 0 || col >= fFeeParam->GetNcol()) 
627     return -1;
628   else 
629     return col;
630 }
631
632 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t bufSize, UInt_t iEv) const
633 {
634   //
635   // Produce raw data stream from this MCM and put in buf
636   // Returns number of words filled, or negative value 
637   // with -1 * number of overflowed words
638   //
639
640   if( !CheckInitialized() ) 
641     return 0;
642
643   UInt_t  x;
644   Int_t   nw  = 0;  // Number of written words
645   Int_t   of  = 0;  // Number of overflowed words
646   Int_t   rawVer   = fFeeParam->GetRAWversion();
647   Int_t **adc;
648   Int_t   nActiveADC = 0;       // number of activated ADC bits in a word
649
650   if( !CheckInitialized() ) 
651     return 0;
652
653   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
654     adc = fADCR;
655   else 
656     adc = fADCF;
657   
658   // Produce MCM header
659   x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
660
661   if (nw < bufSize) {
662     buf[nw++] = x;
663   }
664   else {
665     of++;
666   }
667
668   // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
669   //                            n : unused , c : ADC count, m : selected ADCs
670   if( rawVer >= 3 ) {
671     x = 0;
672     for( Int_t iAdc = 0 ; iAdc < fgkNADC ; iAdc++ ) {
673       if( ~fZSMap[iAdc] != 0 ) { //  0 means not suppressed
674                 x = x | (1 << (iAdc+4) );       // last 4 digit reserved for 1100=0xc
675                 nActiveADC++;           // number of 1 in mmm....m
676       }
677     }
678         x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC;   // nn = 01, ccccc are inverted, 0xc=1100
679
680     if (nw < bufSize) {
681       buf[nw++] = x;
682     }
683     else {
684       of++;
685     }
686   }
687
688   // Produce ADC data. 3 timebins are packed into one 32 bits word
689   // In this version, different ADC channel will NOT share the same word
690
691   UInt_t aa=0, a1=0, a2=0, a3=0;
692
693   for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
694     if( rawVer>= 3 && ~fZSMap[iAdc] == 0 ) continue; // Zero Suppression, 0 means not suppressed
695     aa = !(iAdc & 1) + 2;
696     for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
697       a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] >> fgkAddDigits : 0;
698       a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] >> fgkAddDigits : 0;
699       a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] >> fgkAddDigits : 0;
700       x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
701       if (nw < bufSize) {
702         buf[nw++] = x;
703       }
704       else {
705         of++;
706       }
707     }
708   }
709
710   if( of != 0 ) return -of; else return nw;
711 }
712
713 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t bufSize )
714 {
715   //
716   // Produce tracklet data stream from this MCM and put in buf
717   // Returns number of words filled, or negative value 
718   // with -1 * number of overflowed words
719   //
720
721   if( !CheckInitialized() ) 
722     return 0;
723
724   Int_t   nw  = 0;  // Number of written words
725   Int_t   of  = 0;  // Number of overflowed words
726     
727   // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM 
728   // fMCMT is filled continuously until no more tracklet words available
729
730   for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
731     if (nw < bufSize) 
732       buf[nw++] = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet])->GetTrackletWord();
733     else 
734       of++;
735   }
736   
737   if( of != 0 ) return -of; else return nw;
738 }
739
740 void AliTRDmcmSim::Filter()
741 {
742   //
743   // Filter the raw ADC values. The active filter stages and their
744   // parameters are taken from AliTRDtrapConfig.
745   // The raw data is stored separate from the filtered data. Thus, 
746   // it is possible to run the filters on a set of raw values 
747   // sequentially for parameter tuning.
748   //
749
750   if( !CheckInitialized() ) 
751     return;
752
753   // Apply filters sequentially. Bypass is handled by filters
754   // since counters and internal registers may be updated even 
755   // if the filter is bypassed.
756   // The first filter takes the data from fADCR and 
757   // outputs to fADCF. 
758   
759   // Non-linearity filter not implemented.
760   FilterPedestal();
761   FilterGain();
762   FilterTail();
763   // Crosstalk filter not implemented.
764 }
765
766 void AliTRDmcmSim::FilterPedestalInit(Int_t baseline) 
767 {
768   // Initializes the pedestal filter assuming that the input has 
769   // been constant for a long time (compared to the time constant).
770
771   UShort_t    fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
772
773   for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++)
774     fPedAcc[iAdc] = (baseline << 2) * (1 << fgkFPshifts[fptc]); 
775 }
776
777 UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
778 {
779   // Returns the output of the pedestal filter given the input value.
780   // The output depends on the internal registers and, thus, the 
781   // history of the filter.
782
783   UShort_t    fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
784   UShort_t    fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
785   UShort_t    fpby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPBY); // 0..1 bypass, active low
786
787   UShort_t accumulatorShifted;
788   Int_t correction;
789   UShort_t inpAdd;
790   
791   inpAdd = value + fpnp;
792
793   accumulatorShifted = (fPedAcc[adc] >> fgkFPshifts[fptc]) & 0x3FF;   // 10 bits
794   if (timebin == 0) // the accumulator is disabled in the drift time
795   {
796     correction = (value & 0x3FF) - accumulatorShifted;
797     fPedAcc[adc] = (fPedAcc[adc] + correction) & 0x7FFFFFFF;             // 31 bits
798   }
799   
800   if (fpby == 0)
801     return value;
802
803   if (inpAdd <= accumulatorShifted)
804     return 0;
805   else
806   {
807     inpAdd = inpAdd - accumulatorShifted;
808     if (inpAdd > 0xFFF) 
809       return 0xFFF;
810     else 
811       return inpAdd;
812   }
813 }
814
815 void AliTRDmcmSim::FilterPedestal()
816 {
817   //
818   // Apply pedestal filter
819   //
820   // As the first filter in the chain it reads data from fADCR 
821   // and outputs to fADCF. 
822   // It has only an effect if previous samples have been fed to 
823   // find the pedestal. Currently, the simulation assumes that 
824   // the input has been stable for a sufficiently long time.
825
826   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
827     for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
828       fADCF[iAdc][iTimeBin] = FilterPedestalNextSample(iAdc, iTimeBin, fADCR[iAdc][iTimeBin]);
829     }
830   }
831 }
832
833 void AliTRDmcmSim::FilterGainInit()
834 {
835   // Initializes the gain filter. In this case, only threshold 
836   // counters are reset.
837
838   for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
839     // these are counters which in hardware continue 
840     // until maximum or reset
841     fGainCounterA[iAdc] = 0;
842     fGainCounterB[iAdc] = 0;
843   }
844 }
845
846 UShort_t AliTRDmcmSim::FilterGainNextSample(Int_t adc, UShort_t value)
847 {
848   // Apply the gain filter to the given value.
849   // BEGIN_LATEX O_{i}(t) = #gamma_{i} * I_{i}(t) + a_{i} END_LATEX
850   // The output depends on the internal registers and, thus, the 
851   // history of the filter.
852
853   UShort_t    fgby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGBY); // bypass, active low
854   UShort_t    fgf  = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGF0 + adc)); // 0x700 + (0 & 0x1ff);
855   UShort_t    fga  = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGA0 + adc)); // 40;
856   UShort_t    fgta = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTA); // 20;
857   UShort_t    fgtb = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTB); // 2060;
858
859   UInt_t corr; // corrected value
860
861   value &= 0xFFF;
862   corr = (value * fgf) >> 11;
863   corr = corr > 0xfff ? 0xfff : corr;
864   corr = AddUintClipping(corr, fga, 12);
865
866   // Update threshold counters 
867   // not really useful as they are cleared with every new event
868   if (!((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF)))
869   // stop when full
870   {
871     if (corr >= fgtb) 
872       fGainCounterB[adc]++;
873     else if (corr >= fgta) 
874       fGainCounterA[adc]++;
875   }
876
877   if (fgby == 1)
878     return corr; 
879   else
880     return value;
881 }
882
883 void AliTRDmcmSim::FilterGain()
884 {
885   // Read data from fADCF and apply gain filter.
886
887   for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
888     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
889         fADCF[iAdc][iTimeBin] = FilterGainNextSample(iAdc, fADCF[iAdc][iTimeBin]);
890     }
891   }
892 }
893
894 void AliTRDmcmSim::FilterTailInit(Int_t baseline)
895 {
896   // Initializes the tail filter assuming that the input has 
897   // been at the baseline value (configured by FTFP) for a 
898   // sufficiently long time.
899
900   // exponents and weight calculated from configuration
901   UShort_t    alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
902   UShort_t    lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
903   UShort_t    lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
904
905   Float_t lambdaL = lambdaLong  * 1.0 / (1 << 11);
906   Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
907   Float_t alphaL  = alphaLong   * 1.0 / (1 << 11);
908   Float_t qup, qdn;
909   qup = (1 - lambdaL) * (1 - lambdaS);
910   qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
911   Float_t kdc = qup/qdn;
912
913   Float_t kt, ql, qs;
914   UShort_t aout;
915
916   if (baseline < 0)
917     baseline = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP);
918   
919   ql = lambdaL * (1 - lambdaS) *      alphaL;
920   qs = lambdaS * (1 - lambdaL) * (1 - alphaL);
921
922   for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
923     Int_t value = baseline & 0xFFF;
924     Int_t corr = (value * fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGF0 + iAdc))) >> 11;
925     corr = corr > 0xfff ? 0xfff : corr;
926     corr = AddUintClipping(corr, fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGA0 + iAdc)), 12);
927
928     kt = kdc * baseline;
929     aout = baseline - (UShort_t) kt;
930
931     fTailAmplLong[iAdc]  = (UShort_t) (aout * ql / (ql + qs));
932     fTailAmplShort[iAdc] = (UShort_t) (aout * qs / (ql + qs));
933   }
934 }
935
936 UShort_t AliTRDmcmSim::FilterTailNextSample(Int_t adc, UShort_t value)
937 {
938   // Returns the output of the tail filter for the given input value. 
939   // The output depends on the internal registers and, thus, the 
940   // history of the filter.
941
942   // exponents and weight calculated from configuration
943   UShort_t    alphaLong   = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL);                          // the weight of the long component
944   UShort_t    lambdaLong  = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier of the long component
945   UShort_t    lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier of the short component
946
947   // intermediate signals
948   UInt_t   aDiff;
949   UInt_t   alInpv;
950   UShort_t aQ;
951   UInt_t   tmp;
952   
953   UShort_t inpVolt = value & 0xFFF;    // 12 bits
954       
955   // add the present generator outputs
956   aQ = AddUintClipping(fTailAmplLong[adc], fTailAmplShort[adc], 12);
957
958   // calculate the difference between the input and the generated signal
959   if (inpVolt > aQ) 
960     aDiff = inpVolt - aQ;
961   else                
962     aDiff = 0;
963   
964   // the inputs to the two generators, weighted
965   alInpv = (aDiff * alphaLong) >> 11;
966   
967   // the new values of the registers, used next time
968   // long component
969   tmp = AddUintClipping(fTailAmplLong[adc], alInpv, 12);
970   tmp =  (tmp * lambdaLong) >> 11;
971   fTailAmplLong[adc] = tmp & 0xFFF;
972   // short component
973   tmp = AddUintClipping(fTailAmplShort[adc], aDiff - alInpv, 12);
974   tmp =  (tmp * lambdaShort) >> 11;
975   fTailAmplShort[adc] = tmp & 0xFFF;
976   
977   // the output of the filter
978   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY) == 0) // bypass mode, active low
979     return value;
980   else
981     return aDiff;
982 }
983
984 void AliTRDmcmSim::FilterTail()
985 {
986   // Apply tail cancellation filter to all data. 
987
988   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
989     for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
990       fADCF[iAdc][iTimeBin] = FilterTailNextSample(iAdc, fADCF[iAdc][iTimeBin]);
991     }
992   }
993 }
994
995 void AliTRDmcmSim::ZSMapping()
996 {
997   //
998   // Zero Suppression Mapping implemented in TRAP chip
999   // only implemented for up to 30 timebins
1000   //
1001   // See detail TRAP manual "Data Indication" section:
1002   // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
1003   //
1004
1005   if( !CheckInitialized() ) 
1006     return;
1007
1008   Int_t eBIS = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIS); 
1009   Int_t eBIT = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIT); 
1010   Int_t eBIL = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIL); 
1011   Int_t eBIN = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIN); 
1012
1013   Int_t **adc = fADCF;
1014
1015   for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) 
1016     fZSMap[iAdc] = -1;
1017
1018   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1019     Int_t iAdc; // current ADC channel
1020     Int_t ap;
1021     Int_t ac;
1022     Int_t an;
1023     Int_t mask;
1024     Int_t supp; // suppression of the current channel (low active)
1025     
1026     // ----- first channel -----
1027     iAdc = 0;
1028     
1029     ap = 0;               // previous
1030     ac = adc[iAdc  ][it]; // current
1031     an = adc[iAdc+1][it]; // next
1032     
1033     mask  = ( ac >=  ap && ac >=  an ) ? 0 : 0x1; // peak center detection
1034     mask += ( ap + ac + an > eBIT )    ? 0 : 0x2; // cluster
1035     mask += ( ac > eBIS )              ? 0 : 0x4; // absolute large peak
1036     
1037     supp = (eBIL >> mask) & 1;
1038     
1039     fZSMap[iAdc] &= ~((1-supp) << it);
1040     if( eBIN == 0 ) {  // neighbour sensitivity
1041       fZSMap[iAdc+1] &= ~((1-supp) << it);
1042     }
1043     
1044     // ----- last channel -----
1045     iAdc = fgkNADC - 1;
1046     
1047     ap = adc[iAdc-1][it]; // previous
1048     ac = adc[iAdc  ][it]; // current
1049     an = 0;               // next
1050     
1051     mask  = ( ac >=  ap && ac >=  an ) ? 0 : 0x1; // peak center detection
1052     mask += ( ap + ac + an > eBIT )    ? 0 : 0x2; // cluster
1053     mask += ( ac > eBIS )              ? 0 : 0x4; // absolute large peak
1054     
1055     supp = (eBIL >> mask) & 1;
1056     
1057     fZSMap[iAdc] &= ~((1-supp) << it);
1058     if( eBIN == 0 ) {  // neighbour sensitivity
1059       fZSMap[iAdc-1] &= ~((1-supp) << it);
1060     }
1061     
1062     // ----- middle channels -----
1063     for( iAdc = 1 ; iAdc < fgkNADC-1; iAdc++ ) {
1064       ap = adc[iAdc-1][it]; // previous
1065       ac = adc[iAdc  ][it]; // current
1066       an = adc[iAdc+1][it]; // next
1067       
1068       mask  = ( ac >=  ap && ac >=  an ) ? 0 : 0x1; // peak center detection
1069       mask += ( ap + ac + an > eBIT )    ? 0 : 0x2; // cluster
1070       mask += ( ac > eBIS )              ? 0 : 0x4; // absolute large peak
1071       
1072       supp = (eBIL >> mask) & 1;
1073       
1074       fZSMap[iAdc] &= ~((1-supp) << it);
1075       if( eBIN == 0 ) {  // neighbour sensitivity
1076         fZSMap[iAdc-1] &= ~((1-supp) << it);
1077         fZSMap[iAdc+1] &= ~((1-supp) << it);
1078       }
1079     }
1080
1081   }
1082 }
1083
1084 void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label) 
1085 {
1086   // Add the given hit to the fit register which is lateron used for 
1087   // the tracklet calculation. 
1088   // In addition to the fit sums in the fit register MC information 
1089   // is stored.
1090
1091   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)) && 
1092       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0)))
1093     fFitReg[adc].fQ0 += qtot;
1094   
1095   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) && 
1096       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)))
1097     fFitReg[adc].fQ1 += qtot;
1098   
1099   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) && 
1100       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE)))
1101   {
1102     fFitReg[adc].fSumX  += timebin;
1103     fFitReg[adc].fSumX2 += timebin*timebin;
1104     fFitReg[adc].fNhits++;
1105     fFitReg[adc].fSumY  += ypos;
1106     fFitReg[adc].fSumY2 += ypos*ypos;
1107     fFitReg[adc].fSumXY += timebin*ypos;
1108   }
1109
1110   // register hits (MC info)
1111   fHits[fNHits].fChannel = adc;
1112   fHits[fNHits].fQtot = qtot;
1113   fHits[fNHits].fYpos = ypos;
1114   fHits[fNHits].fTimebin = timebin;
1115   fHits[fNHits].fLabel = label;
1116   fNHits++;
1117 }
1118
1119 void AliTRDmcmSim::CalcFitreg() 
1120 {
1121   // Preprocessing.
1122   // Detect the hits and fill the fit registers.
1123   // Requires 12-bit data from fADCF which means Filter() 
1124   // has to be called before even if all filters are bypassed.
1125
1126   //??? to be clarified:
1127   UInt_t adcMask = 0xffffffff;
1128   
1129   UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp;
1130   Short_t ypos, fromLeft, fromRight, found;
1131   UShort_t qTotal[19]; // the last is dummy
1132   UShort_t marked[6], qMarked[6], worse1, worse2;
1133   
1134   timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS); 
1135   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0) 
1136       < timebin1)
1137     timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0);
1138   timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE); 
1139   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1) 
1140       > timebin2)
1141     timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1);
1142
1143   // reset the fit registers
1144   fNHits = 0; 
1145   for (adcch = 0; adcch < fgkNADC-2; adcch++) // due to border channels
1146   {
1147     fFitReg[adcch].fNhits = 0;
1148     fFitReg[adcch].fQ0    = 0;
1149     fFitReg[adcch].fQ1    = 0;
1150     fFitReg[adcch].fSumX  = 0;
1151     fFitReg[adcch].fSumY  = 0;
1152     fFitReg[adcch].fSumX2 = 0;
1153     fFitReg[adcch].fSumY2 = 0;
1154     fFitReg[adcch].fSumXY = 0;
1155   }
1156   
1157   for (timebin = timebin1; timebin < timebin2; timebin++)
1158   {
1159     // first find the hit candidates and store the total cluster charge in qTotal array
1160     // in case of not hit store 0 there.
1161     for (adcch = 0; adcch < fgkNADC-2; adcch++) {
1162       if ( ( (adcMask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS
1163       {
1164         adcLeft  = fADCF[adcch  ][timebin];
1165         adcCentral  = fADCF[adcch+1][timebin];
1166         adcRight = fADCF[adcch+2][timebin];
1167         if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY) == 1) 
1168           hitQual = ( (adcLeft * adcRight) < 
1169                        (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT) * adcCentral) );
1170         else            
1171           hitQual = 1;
1172         // The accumulated charge is with the pedestal!!!
1173         qtotTemp = adcLeft + adcCentral + adcRight;
1174         if ( (hitQual) &&
1175              (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) &&
1176              (adcLeft <= adcCentral) &&
1177              (adcCentral > adcRight) )
1178           qTotal[adcch] = qtotTemp;
1179         else
1180           qTotal[adcch] = 0;
1181       }
1182       else
1183         qTotal[adcch] = 0; //jkl
1184       if (qTotal[adcch] != 0) 
1185         AliDebug(10,Form("ch %2d   qTotal %5d",adcch, qTotal[adcch]));
1186     }
1187
1188     fromLeft = -1;
1189     adcch = 0;
1190     found = 0;
1191     marked[4] = 19; // invalid channel
1192     marked[5] = 19; // invalid channel
1193     qTotal[19] = 0;
1194     while ((adcch < 16) && (found < 3))
1195     {
1196       if (qTotal[adcch] > 0)
1197       {
1198         fromLeft = adcch;
1199         marked[2*found+1]=adcch;
1200         found++;
1201       }
1202       adcch++;
1203     }
1204     
1205     fromRight = -1;
1206     adcch = 18;
1207     found = 0;
1208     while ((adcch > 2) && (found < 3))
1209     {
1210       if (qTotal[adcch] > 0)
1211       {
1212         marked[2*found]=adcch;
1213         found++;
1214         fromRight = adcch;
1215       }
1216       adcch--;
1217     }
1218
1219     AliDebug(10,Form("Fromleft=%d, Fromright=%d",fromLeft, fromRight));
1220     // here mask the hit candidates in the middle, if any
1221     if ((fromLeft >= 0) && (fromRight >= 0) && (fromLeft < fromRight))
1222       for (adcch = fromLeft+1; adcch < fromRight; adcch++)
1223         qTotal[adcch] = 0;
1224     
1225     found = 0;
1226     for (adcch = 0; adcch < 19; adcch++)
1227       if (qTotal[adcch] > 0) found++;
1228     // NOT READY
1229
1230     if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates!
1231     {
1232       if (marked[4] == marked[5]) marked[5] = 19;
1233       for (found=0; found<6; found++)
1234       {
1235         qMarked[found] = qTotal[marked[found]] >> 4;
1236         AliDebug(10,Form("ch_%d qTotal %d qTotals %d",marked[found],qTotal[marked[found]],qMarked[found]));
1237       }
1238       
1239       Sort6To2Worst(marked[0], marked[3], marked[4], marked[1], marked[2], marked[5],
1240                     qMarked[0],
1241                     qMarked[3],
1242                     qMarked[4],
1243                     qMarked[1],
1244                     qMarked[2],
1245                     qMarked[5],
1246                     &worse1, &worse2);
1247       // Now mask the two channels with the smallest charge
1248       if (worse1 < 19)
1249       {
1250         qTotal[worse1] = 0;
1251         AliDebug(10,Form("Kill ch %d\n",worse1));
1252       }
1253       if (worse2 < 19)
1254       {
1255         qTotal[worse2] = 0;
1256         AliDebug(10,Form("Kill ch %d\n",worse2));
1257       }
1258     }
1259     
1260     for (adcch = 0; adcch < 19; adcch++) {
1261       if (qTotal[adcch] > 0) // the channel is marked for processing
1262       {
1263         adcLeft  = fADCF[adcch  ][timebin];
1264         adcCentral  = fADCF[adcch+1][timebin];
1265         adcRight = fADCF[adcch+2][timebin];
1266         // hit detected, in TRAP we have 4 units and a hit-selection, here we proceed all channels!
1267         // subtract the pedestal TPFP, clipping instead of wrapping
1268         
1269         Int_t regTPFP = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP);
1270         AliDebug(10, Form("Hit found, time=%d, adcch=%d/%d/%d, adc values=%d/%d/%d, regTPFP=%d, TPHT=%d\n",
1271                timebin, adcch, adcch+1, adcch+2, adcLeft, adcCentral, adcRight, regTPFP, 
1272                fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)));
1273
1274         if (adcLeft  < regTPFP) adcLeft  = 0; else adcLeft  -= regTPFP;
1275         if (adcCentral  < regTPFP) adcCentral  = 0; else adcCentral  -= regTPFP;
1276         if (adcRight < regTPFP) adcRight = 0; else adcRight -= regTPFP;
1277
1278         // Calculate the center of gravity
1279         // checking for adcCentral != 0 (in case of "bad" configuration)
1280         if (adcCentral == 0)
1281           continue;
1282         ypos = 128*(adcLeft - adcRight) / adcCentral;
1283         if (ypos < 0) ypos = -ypos;
1284         // make the correction using the position LUT
1285         ypos = ypos + fTrapConfig->GetTrapReg((AliTRDtrapConfig::TrapReg_t) (AliTRDtrapConfig::kTPL00 + (ypos & 0x7F))); 
1286         if (adcLeft > adcRight) ypos = -ypos;
1287
1288         // label calculation
1289         Int_t mcLabel = -1;
1290         if (fDigitsManager) {
1291           Int_t label[9] = { 0 }; // up to 9 different labels possible
1292           Int_t count[9] = { 0 };
1293           Int_t maxIdx = -1;
1294           Int_t maxCount = 0;
1295           Int_t nLabels = 0;
1296           Int_t padcol[3]; 
1297           padcol[0] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch);
1298           padcol[1] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch+1);
1299           padcol[2] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch+2);
1300           Int_t padrow = fFeeParam->GetPadRowFromMCM(fRobPos, fMcmPos);
1301           for (Int_t iDict = 0; iDict < 3; iDict++) {
1302             if (!fDict[iDict])
1303               continue;
1304             for (Int_t iPad = 0; iPad < 3; iPad++) {
1305               if (padcol[iPad] < 0) 
1306                 continue;
1307               Int_t currLabel = fDict[iDict]->GetData(padrow, padcol[iPad], timebin); //fDigitsManager->GetTrack(iDict, padrow, padcol, timebin, fDetector);
1308               AliDebug(10, Form("Read label: %4i for det: %3i, row: %i, col: %i, tb: %i\n", currLabel, fDetector, padrow, padcol[iPad], timebin));
1309               for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
1310                 if (currLabel == label[iLabel]) {
1311                   count[iLabel]++;
1312                   if (count[iLabel] > maxCount) {
1313                     maxCount = count[iLabel];
1314                     maxIdx = iLabel;
1315                   }
1316                   currLabel = 0;
1317                   break;
1318                 }
1319               } 
1320               if (currLabel > 0) {
1321                 label[nLabels++] = currLabel;
1322               }
1323             }
1324           }
1325           if (maxIdx >= 0)
1326             mcLabel = label[maxIdx];
1327         }
1328
1329         // add the hit to the fitregister
1330         AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, mcLabel);
1331       }
1332     }
1333   }
1334
1335   for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
1336     if (fFitReg[iAdc].fNhits != 0) {
1337       AliDebug(2, Form("fitreg[%i]: nHits = %i, sumX = %i, sumY = %i, sumX2 = %i, sumY2 = %i, sumXY = %i", iAdc,
1338                        fFitReg[iAdc].fNhits,
1339                        fFitReg[iAdc].fSumX,
1340                        fFitReg[iAdc].fSumY,
1341                        fFitReg[iAdc].fSumX2,
1342                        fFitReg[iAdc].fSumY2,
1343                        fFitReg[iAdc].fSumXY
1344                  ));
1345     }
1346   }
1347 }
1348
1349 void AliTRDmcmSim::TrackletSelection() 
1350 {
1351   // Select up to 4 tracklet candidates from the fit registers  
1352   // and assign them to the CPUs.
1353
1354   UShort_t adcIdx, i, j, ntracks, tmp;
1355   UShort_t trackletCand[18][2]; // store the adcch[0] and number of hits[1] for all tracklet candidates
1356
1357   ntracks = 0;
1358   for (adcIdx = 0; adcIdx < 18; adcIdx++) // ADCs
1359     if ( (fFitReg[adcIdx].fNhits 
1360           >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL)) &&
1361          (fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits
1362           >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT)))
1363     {
1364       trackletCand[ntracks][0] = adcIdx;
1365       trackletCand[ntracks][1] = fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits;
1366       AliDebug(10,Form("%d  %2d %4d\n", ntracks, trackletCand[ntracks][0], trackletCand[ntracks][1]));
1367       ntracks++;
1368     };
1369
1370   for (i=0; i<ntracks;i++) 
1371     AliDebug(10,Form("%d %d %d\n",i,trackletCand[i][0], trackletCand[i][1]));
1372
1373   if (ntracks > 4)
1374   {
1375     // primitive sorting according to the number of hits
1376     for (j = 0; j < (ntracks-1); j++)
1377     {
1378       for (i = j+1; i < ntracks; i++)
1379       {
1380         if ( (trackletCand[j][1]  < trackletCand[i][1]) ||
1381              ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) )
1382         {
1383           // swap j & i
1384           tmp = trackletCand[j][1];
1385           trackletCand[j][1] = trackletCand[i][1];
1386           trackletCand[i][1] = tmp;
1387           tmp = trackletCand[j][0];
1388           trackletCand[j][0] = trackletCand[i][0];
1389           trackletCand[i][0] = tmp;
1390         }
1391       }
1392     }
1393     ntracks = 4; // cut the rest, 4 is the max
1394   }
1395   // else is not necessary to sort
1396   
1397   // now sort, so that the first tracklet going to CPU0 corresponds to the highest adc channel - as in the TRAP
1398   for (j = 0; j < (ntracks-1); j++)
1399   {
1400     for (i = j+1; i < ntracks; i++)
1401     {
1402       if (trackletCand[j][0] < trackletCand[i][0])
1403       {
1404         // swap j & i
1405         tmp = trackletCand[j][1];
1406         trackletCand[j][1] = trackletCand[i][1];
1407         trackletCand[i][1] = tmp;
1408         tmp = trackletCand[j][0];
1409         trackletCand[j][0] = trackletCand[i][0];
1410         trackletCand[i][0] = tmp;
1411       }
1412     }
1413   }
1414   for (i = 0; i < ntracks; i++)  // CPUs with tracklets.
1415     fFitPtr[i] = trackletCand[i][0]; // pointer to the left channel with tracklet for CPU[i]
1416   for (i = ntracks; i < 4; i++)  // CPUs without tracklets
1417     fFitPtr[i] = 31;            // pointer to the left channel with tracklet for CPU[i] = 31 (invalid)
1418   AliDebug(10,Form("found %i tracklet candidates\n", ntracks));
1419   for (i = 0; i < 4; i++)
1420     AliDebug(10,Form("fitPtr[%i]: %i\n", i, fFitPtr[i]));
1421 }
1422
1423 void AliTRDmcmSim::FitTracklet()
1424 {
1425   // Perform the actual tracklet fit based on the fit sums 
1426   // which have been filled in the fit registers. 
1427
1428   // parameters in fitred.asm (fit program)
1429   Int_t decPlaces = 5;
1430   Int_t rndAdd = 0;
1431   if (decPlaces >  1) 
1432     rndAdd = (1 << (decPlaces-1)) + 1;
1433   else if (decPlaces == 1)
1434     rndAdd = 1;
1435   Int_t ndriftDp = 5;  // decimal places for drift time
1436   Long64_t shift = ((Long64_t) 1 << 32);
1437
1438   // calculated in fitred.asm
1439   Int_t padrow = ((fRobPos >> 1) << 2) | (fMcmPos >> 2);
1440   Int_t yoffs = (((((fRobPos & 0x1) << 2) + (fMcmPos & 0x3)) * 18) << 8) - 
1441     ((18*4*2 - 18*2 - 1) << 7);
1442   yoffs = yoffs << decPlaces; // holds position of ADC channel 1
1443   Int_t layer = fDetector % 6;
1444   UInt_t scaleY = (UInt_t) ((0.635 + 0.03 * layer)/(256.0 * 160.0e-4) * shift);
1445   UInt_t scaleD = (UInt_t) ((0.635 + 0.03 * layer)/(256.0 * 140.0e-4) * shift);
1446
1447   Int_t deflCorr = fTrapConfig->GetDmem(0xc022, fDetector, fRobPos, fMcmPos);
1448   Int_t ndrift   = fTrapConfig->GetDmem(0xc025, fDetector, fRobPos, fMcmPos); 
1449
1450   // local variables for calculation
1451   Long64_t mult, temp, denom; //???
1452   UInt_t q0, q1, qTotal;          // charges in the two windows and total charge
1453   UShort_t nHits;                 // number of hits
1454   Int_t slope, offset;            // slope and offset of the tracklet
1455   Int_t sumX, sumY, sumXY, sumX2; // fit sums from fit registers
1456   //int32_t SumY2;                // not used in the current TRAP program
1457   FitReg_t *fit0, *fit1;          // pointers to relevant fit registers
1458   
1459 //  const uint32_t OneDivN[32] = {  // 2**31/N : exactly like in the TRAP, the simple division here gives the same result!
1460 //      0x00000000, 0x80000000, 0x40000000, 0x2AAAAAA0, 0x20000000, 0x19999990, 0x15555550, 0x12492490,
1461 //      0x10000000, 0x0E38E380, 0x0CCCCCC0, 0x0BA2E8B0, 0x0AAAAAA0, 0x09D89D80, 0x09249240, 0x08888880,
1462 //      0x08000000, 0x07878780, 0x071C71C0, 0x06BCA1A0, 0x06666660, 0x06186180, 0x05D17450, 0x0590B210,
1463 //      0x05555550, 0x051EB850, 0x04EC4EC0, 0x04BDA120, 0x04924920, 0x0469EE50, 0x04444440, 0x04210840};
1464
1465   for (Int_t cpu = 0; cpu < 4; cpu++) {
1466     if (fFitPtr[cpu] == 31)
1467     {
1468       fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker(); 
1469     }
1470     else
1471     {
1472       fit0 = &fFitReg[fFitPtr[cpu]  ];
1473       fit1 = &fFitReg[fFitPtr[cpu]+1]; // next channel
1474
1475       mult = 1;
1476       mult = mult << (32 + decPlaces);
1477       mult = -mult;
1478
1479       // Merging
1480       nHits   = fit0->fNhits + fit1->fNhits; // number of hits
1481       sumX    = fit0->fSumX  + fit1->fSumX;
1482       sumX2   = fit0->fSumX2 + fit1->fSumX2;
1483       denom   = nHits*sumX2 - sumX*sumX;
1484
1485       mult    = mult / denom; // exactly like in the TRAP program
1486       q0      = fit0->fQ0    + fit1->fQ0;
1487       q1      = fit0->fQ1    + fit1->fQ1;
1488       sumY    = fit0->fSumY  + fit1->fSumY  + 256*fit1->fNhits;
1489       sumXY   = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX;
1490
1491       slope   = nHits*sumXY - sumX * sumY;
1492       offset  = sumX2*sumY  - sumX * sumXY;
1493       temp    = mult * slope;
1494       slope   = temp >> 32; // take the upper 32 bits
1495       slope   = -slope;
1496       temp    = mult * offset;
1497       offset  = temp >> 32; // take the upper 32 bits
1498
1499       offset = offset + yoffs;
1500       AliDebug(10, Form("slope = %i, slope * ndrift = %i, deflCorr: %i", 
1501                        slope, slope * ndrift, deflCorr));
1502       slope  = ((slope * ndrift) >> ndriftDp) + deflCorr;
1503       offset = offset - (fFitPtr[cpu] << (8 + decPlaces));
1504       
1505       temp    = slope;
1506       temp    = temp * scaleD;
1507       slope   = (temp >> 32);
1508       temp    = offset;
1509       temp    = temp * scaleY;
1510       offset  = (temp >> 32);
1511         
1512       // rounding, like in the TRAP
1513       slope   = (slope  + rndAdd) >> decPlaces;
1514       offset  = (offset + rndAdd) >> decPlaces;
1515
1516       AliDebug(5, Form("Det: %3i, ROB: %i, MCM: %2i: deflection: %i, min: %i, max: %i", 
1517                        fDetector, fRobPos, fMcmPos, slope, 
1518                        fTrapConfig->GetDmem(0xc030 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos), 
1519                        fTrapConfig->GetDmem(0xc031 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos)));
1520
1521       Bool_t rejected = kFALSE;
1522       // deflection range table from DMEM
1523       if ((slope < fTrapConfig->GetDmem(0xc030 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos)) || 
1524           (slope > fTrapConfig->GetDmem(0xc031 + 2*fFitPtr[cpu], fDetector, fRobPos, fMcmPos)))
1525         rejected = kTRUE;
1526
1527       if (rejected && GetApplyCut())
1528       {
1529         fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1530       }
1531       else
1532       {
1533         if (slope > 63 || slope < -64) { // wrapping in TRAP!
1534           AliError(Form("Overflow in slope: %i, tracklet discarded!", slope));
1535           fMCMT[cpu] = 0x10001000;
1536           continue;
1537         }
1538
1539         slope   = slope  &   0x7F; // 7 bit
1540         
1541         if (offset > 0xfff || offset < -0xfff) 
1542           AliWarning("Overflow in offset");
1543         offset  = offset & 0x1FFF; // 13 bit
1544
1545         qTotal = 0; // set to zero as long as no reasonable PID calculation is available 
1546                     // before: GetPID(q0/length/fgChargeNorm, q1/length/fgChargeNorm);
1547
1548         if (qTotal > 0xff)
1549           AliWarning("Overflow in charge");
1550         qTotal  = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program
1551         
1552         // assemble and store the tracklet word
1553         fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
1554
1555         // calculate MC label
1556         Int_t mcLabel = -1;
1557         Int_t nHits0 = 0;
1558         Int_t nHits1 = 0;
1559         if (fDigitsManager) {
1560           Int_t label[30] = {0}; // up to 30 different labels possible
1561           Int_t count[30] = {0};
1562           Int_t maxIdx = -1;
1563           Int_t maxCount = 0;
1564           Int_t nLabels = 0;
1565           for (Int_t iHit = 0; iHit < fNHits; iHit++) {
1566             if ((fHits[iHit].fChannel - fFitPtr[cpu] < 0) ||
1567                 (fHits[iHit].fChannel - fFitPtr[cpu] > 1))
1568               continue;
1569
1570             // counting contributing hits
1571             if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0) &&
1572                 fHits[iHit].fTimebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0))
1573               nHits0++;
1574             if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1) &&
1575                 fHits[iHit].fTimebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1))
1576               nHits1++;
1577
1578             Int_t currLabel = fHits[iHit].fLabel;
1579             for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
1580               if (currLabel == label[iLabel]) {
1581                 count[iLabel]++;
1582                 if (count[iLabel] > maxCount) {
1583                   maxCount = count[iLabel];
1584                   maxIdx = iLabel;
1585                 }
1586                 currLabel = 0;
1587                 break;
1588               }
1589             }
1590             if (currLabel > 0) {
1591               label[nLabels++] = currLabel;
1592             }
1593           }
1594           if (maxIdx >= 0)
1595             mcLabel = label[maxIdx];
1596         }
1597         new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
1598         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetLabel(mcLabel);
1599
1600        
1601         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits(fit0->fNhits + fit1->fNhits);
1602         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits0(nHits0);
1603         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits1(nHits1);
1604         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ0(q0);
1605         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ1(q1);
1606       }
1607     }
1608   }
1609 }
1610
1611 void AliTRDmcmSim::Tracklet()
1612 {
1613   // Run the tracklet calculation by calling sequentially:
1614   // CalcFitreg(); TrackletSelection(); FitTracklet()
1615   // and store the tracklets 
1616
1617   if (!fInitialized) {
1618     AliError("Called uninitialized! Nothing done!");
1619     return;
1620   }
1621
1622   fTrackletArray->Delete();
1623
1624   CalcFitreg();
1625   if (fNHits == 0)
1626     return;
1627   TrackletSelection();
1628   FitTracklet();
1629 }
1630
1631 Bool_t AliTRDmcmSim::StoreTracklets() 
1632 {
1633   // store the found tracklets via the loader
1634
1635   if (fTrackletArray->GetEntriesFast() == 0) 
1636     return kTRUE;
1637
1638   AliRunLoader *rl = AliRunLoader::Instance();
1639   AliDataLoader *dl = 0x0;
1640   if (rl)
1641     dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets");
1642   if (!dl) {
1643     AliError("Could not get the tracklets data loader!");
1644     return kFALSE;
1645   }
1646
1647   TTree *trackletTree = dl->Tree();
1648   if (!trackletTree) {
1649     dl->MakeTree();
1650     trackletTree = dl->Tree();
1651   }
1652   
1653   AliTRDtrackletMCM *trkl = 0x0;
1654   TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
1655   if (!trkbranch)
1656     trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
1657   
1658   for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
1659     trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]);
1660     trkbranch->SetAddress(&trkl);
1661     trkbranch->Fill();
1662   }
1663
1664   return kTRUE;
1665 }
1666
1667 void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
1668 {
1669   // write back the processed data configured by EBSF
1670   // EBSF = 1: unfiltered data; EBSF = 0: filtered data
1671   // zero-suppressed valued are written as -1 to digits
1672
1673   if( !CheckInitialized() ) 
1674     return;
1675
1676   Int_t offset = (fMcmPos % 4 + 1) * 21 + (fRobPos % 2) * 84 - 1;
1677
1678   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
1679   {
1680     for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
1681       if (~fZSMap[iAdc] == 0) {
1682         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1683           digits->SetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin, -1);
1684         }
1685       }
1686     }
1687   }
1688   else {
1689     for (Int_t iAdc = 0; iAdc < fgkNADC; iAdc++) {
1690       if (~fZSMap[iAdc] != 0) {
1691         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1692           digits->SetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin, (fADCF[iAdc][iTimeBin] >> fgkAddDigits) - fgAddBaseline);
1693         }
1694       }
1695       else {
1696         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1697           digits->SetDataByAdcCol(GetRow(), offset - iAdc, iTimeBin, -1);
1698         }
1699       }
1700     }
1701   }
1702 }
1703
1704 // help functions, to be cleaned up
1705
1706 UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
1707 {
1708   // 
1709   // This function adds a and b (unsigned) and clips to 
1710   // the specified number of bits. 
1711   //  
1712
1713   UInt_t sum = a + b;
1714   if (nbits < 32)
1715   {
1716     UInt_t maxv = (1 << nbits) - 1;;
1717     if (sum > maxv) 
1718       sum = maxv;
1719   }
1720   else
1721   {
1722     if ((sum < a) || (sum < b)) 
1723       sum = 0xFFFFFFFF;
1724   }
1725   return sum;
1726 }
1727
1728 void AliTRDmcmSim::Sort2(UShort_t  idx1i, UShort_t  idx2i, \
1729                             UShort_t  val1i, UShort_t  val2i, \
1730                             UShort_t *idx1o, UShort_t *idx2o, \
1731                             UShort_t *val1o, UShort_t *val2o) const
1732 {
1733   // sorting for tracklet selection
1734
1735     if (val1i > val2i)
1736     {
1737         *idx1o = idx1i;
1738         *idx2o = idx2i;
1739         *val1o = val1i;
1740         *val2o = val2i;
1741     }
1742     else
1743     {
1744         *idx1o = idx2i;
1745         *idx2o = idx1i;
1746         *val1o = val2i;
1747         *val2o = val1i;
1748     }
1749 }
1750
1751 void AliTRDmcmSim::Sort3(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, \
1752                             UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, \
1753                             UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, \
1754                             UShort_t *val1o, UShort_t *val2o, UShort_t *val3o)
1755 {
1756   // sorting for tracklet selection
1757
1758     Int_t sel;
1759
1760
1761     if (val1i > val2i) sel=4; else sel=0;
1762     if (val2i > val3i) sel=sel + 2;
1763     if (val3i > val1i) sel=sel + 1;
1764     switch(sel)
1765     {
1766         case 6 : // 1 >  2  >  3            => 1 2 3
1767         case 0 : // 1 =  2  =  3            => 1 2 3 : in this case doesn't matter, but so is in hardware!
1768             *idx1o = idx1i;
1769             *idx2o = idx2i;
1770             *idx3o = idx3i;
1771             *val1o = val1i;
1772             *val2o = val2i;
1773             *val3o = val3i;
1774             break;
1775
1776         case 4 : // 1 >  2, 2 <= 3, 3 <= 1  => 1 3 2
1777             *idx1o = idx1i;
1778             *idx2o = idx3i;
1779             *idx3o = idx2i;
1780             *val1o = val1i;
1781             *val2o = val3i;
1782             *val3o = val2i;
1783             break;
1784
1785         case 2 : // 1 <= 2, 2 > 3, 3 <= 1   => 2 1 3
1786             *idx1o = idx2i;
1787             *idx2o = idx1i;
1788             *idx3o = idx3i;
1789             *val1o = val2i;
1790             *val2o = val1i;
1791             *val3o = val3i;
1792             break;
1793
1794         case 3 : // 1 <= 2, 2 > 3, 3  > 1   => 2 3 1
1795             *idx1o = idx2i;
1796             *idx2o = idx3i;
1797             *idx3o = idx1i;
1798             *val1o = val2i;
1799             *val2o = val3i;
1800             *val3o = val1i;
1801             break;
1802
1803         case 1 : // 1 <= 2, 2 <= 3, 3 > 1   => 3 2 1
1804             *idx1o = idx3i;
1805             *idx2o = idx2i;
1806             *idx3o = idx1i;
1807             *val1o = val3i;
1808             *val2o = val2i;
1809             *val3o = val1i;
1810         break;
1811
1812         case 5 : // 1 > 2, 2 <= 3, 3 >  1   => 3 1 2
1813             *idx1o = idx3i;
1814             *idx2o = idx1i;
1815             *idx3o = idx2i;
1816             *val1o = val3i;
1817             *val2o = val1i;
1818             *val3o = val2i;
1819         break;
1820
1821         default: // the rest should NEVER happen!
1822             AliError("ERROR in Sort3!!!\n");
1823         break;
1824     }
1825 }
1826
1827 void AliTRDmcmSim::Sort6To4(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, UShort_t  idx4i, UShort_t  idx5i, UShort_t  idx6i, \
1828                                UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, UShort_t  val4i, UShort_t  val5i, UShort_t  val6i, \
1829                                UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, UShort_t *idx4o, \
1830                                UShort_t *val1o, UShort_t *val2o, UShort_t *val3o, UShort_t *val4o)
1831 {
1832   // sorting for tracklet selection
1833
1834     UShort_t idx21s, idx22s, idx23s, dummy;
1835     UShort_t val21s, val22s, val23s;
1836     UShort_t idx23as, idx23bs;
1837     UShort_t val23as, val23bs;
1838
1839     Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
1840                  idx1o, &idx21s, &idx23as,
1841                  val1o, &val21s, &val23as);
1842
1843     Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1844                  idx2o, &idx22s, &idx23bs,
1845                  val2o, &val22s, &val23bs);
1846
1847     Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, &dummy, &val23s, &dummy);
1848
1849     Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
1850                  idx3o, idx4o, &dummy,
1851                  val3o, val4o, &dummy);
1852
1853 }
1854
1855 void AliTRDmcmSim::Sort6To2Worst(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, UShort_t  idx4i, UShort_t  idx5i, UShort_t  idx6i, \
1856                                     UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, UShort_t  val4i, UShort_t  val5i, UShort_t  val6i, \
1857                                     UShort_t *idx5o, UShort_t *idx6o)
1858 {
1859   // sorting for tracklet selection
1860
1861     UShort_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5;
1862     UShort_t val21s, val22s, val23s;
1863     UShort_t idx23as, idx23bs;
1864     UShort_t val23as, val23bs;
1865
1866     Sort3(idx1i, idx2i,   idx3i, val1i, val2i, val3i,
1867                  &dummy1, &idx21s, &idx23as,
1868                  &dummy2, &val21s, &val23as);
1869
1870     Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1871                  &dummy1, &idx22s, &idx23bs,
1872                  &dummy2, &val22s, &val23bs);
1873
1874     Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, idx5o, &val23s, &dummy1);
1875
1876     Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
1877                  &dummy1, &dummy2, idx6o,
1878                  &dummy3, &dummy4, &dummy5);
1879 }
1880
1881
1882 // ----- I/O implementation -----
1883
1884 ostream& AliTRDmcmSim::text(ostream& os)
1885 {
1886   // manipulator to activate output in text format (default)
1887
1888   os.iword(fgkFormatIndex) = 0;
1889   return os;
1890 }
1891
1892 ostream& AliTRDmcmSim::cfdat(ostream& os)
1893 {
1894   // manipulator to activate output in CFDAT format 
1895   // to send to the FEE via SCSN
1896
1897   os.iword(fgkFormatIndex) = 1; 
1898   return os;
1899 }
1900
1901 ostream& AliTRDmcmSim::raw(ostream& os)
1902 {
1903   // manipulator to activate output as raw data dump
1904
1905   os.iword(fgkFormatIndex) = 2;
1906   return os;
1907 }
1908
1909 ostream& operator<<(ostream& os, const AliTRDmcmSim& mcm)
1910 {
1911   // output implementation
1912   
1913   // no output for non-initialized MCM
1914   if (!mcm.CheckInitialized())
1915     return os;
1916
1917   // ----- human-readable output -----
1918   if (os.iword(AliTRDmcmSim::fgkFormatIndex) == 0) {
1919     
1920     os << "MCM " << mcm.fMcmPos << " on ROB " << mcm.fRobPos << 
1921       " in detector " << mcm.fDetector << std::endl;
1922     
1923     os << "----- Unfiltered ADC data (10 bit) -----" << std::endl;
1924     os << "ch    ";
1925     for (Int_t iChannel = 0; iChannel < mcm.fgkNADC; iChannel++) 
1926       os << std::setw(5) << iChannel;
1927     os << std::endl;
1928     for (Int_t iTimeBin = 0; iTimeBin < mcm.fNTimeBin; iTimeBin++) {
1929       os << "tb " << std::setw(2) << iTimeBin << ":";
1930       for (Int_t iChannel = 0; iChannel < mcm.fgkNADC; iChannel++) {
1931         os << std::setw(5) << (mcm.fADCR[iChannel][iTimeBin] >> mcm.fgkAddDigits);
1932       }
1933       os << std::endl;
1934     }
1935     
1936     os << "----- Filtered ADC data (10+2 bit) -----" << std::endl;
1937     os << "ch    ";
1938     for (Int_t iChannel = 0; iChannel < mcm.fgkNADC; iChannel++) 
1939       os << std::setw(4) << iChannel
1940          << ((~mcm.fZSMap[iChannel] != 0) ? "!" : " ");
1941     os << std::endl;
1942     for (Int_t iTimeBin = 0; iTimeBin < mcm.fNTimeBin; iTimeBin++) {
1943       os << "tb " << std::setw(2) << iTimeBin << ":";
1944       for (Int_t iChannel = 0; iChannel < mcm.fgkNADC; iChannel++) {
1945         os << std::setw(4) << (mcm.fADCF[iChannel][iTimeBin])
1946            << (((mcm.fZSMap[iChannel] & (1 << iTimeBin)) == 0) ? "!" : " ");
1947       }
1948       os << std::endl;
1949     }
1950   }
1951
1952   // ----- CFDAT output -----
1953   else if(os.iword(AliTRDmcmSim::fgkFormatIndex) == 1) {
1954     Int_t dest       = 127;
1955     Int_t addrOffset = 0x2000;
1956     Int_t addrStep   = 0x80;
1957     
1958     for (Int_t iTimeBin = 0; iTimeBin < mcm.fNTimeBin; iTimeBin++) {
1959       for (Int_t iChannel = 0; iChannel < mcm.fgkNADC; iChannel++) {
1960         os << std::setw(5) << 10 
1961            << std::setw(5) << addrOffset + iChannel * addrStep + iTimeBin 
1962            << std::setw(5) << (mcm.fADCF[iChannel][iTimeBin])
1963            << std::setw(5) << dest << std::endl;
1964       }
1965       os << std::endl;
1966     }
1967   }
1968
1969   // ----- raw data ouptut -----
1970   else if (os.iword(AliTRDmcmSim::fgkFormatIndex) == 2) {
1971     Int_t   bufSize   = 300;
1972     UInt_t *buf       = new UInt_t[bufSize];
1973     
1974     Int_t bufLength   = mcm.ProduceRawStream(&buf[0], bufSize);
1975     
1976     for (Int_t i = 0; i < bufLength; i++) 
1977       std::cout << "0x" << std::hex << buf[i] << std::endl;
1978     
1979     delete [] buf;
1980   }
1981
1982   else {
1983     os << "unknown format set" << std::endl;
1984   }
1985
1986   return os;
1987 }