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