Update by Jochen to use new AliTRDmcmSim implementation
[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() ) 
554     return -1;
555
556   Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
557   if (col < 0 || col >= fFeeParam->GetNcol()) 
558     return -1;
559   else 
560     return col;
561 }
562
563 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv)
564 {
565   //
566   // Produce raw data stream from this MCM and put in buf
567   // Returns number of words filled, or negative value 
568   // with -1 * number of overflowed words
569   //
570
571   UInt_t  x;
572   Int_t   nw  = 0;  // Number of written words
573   Int_t   of  = 0;  // Number of overflowed words
574   Int_t   rawVer   = fFeeParam->GetRAWversion();
575   Int_t **adc;
576   Int_t   nActiveADC = 0;       // number of activated ADC bits in a word
577
578   if( !CheckInitialized() ) return 0;
579
580   if( fFeeParam->GetRAWstoreRaw() ) {
581     adc = fADCR;
582   } else {
583     adc = fADCF;
584   }
585
586   // Produce MCM header
587   x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
588
589   if (nw < maxSize) {
590     buf[nw++] = x;
591     //printf("\nMCM header: %X ",x);
592   }
593   else {
594     of++;
595   }
596
597   // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
598   //                            n : unused , c : ADC count, m : selected ADCs
599   if( rawVer >= 3 ) {
600     x = 0;
601     for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
602       if( fZSM1Dim[iAdc] == 0 ) { //  0 means not suppressed
603                 x = x | (1 << (iAdc+4) );       // last 4 digit reserved for 1100=0xc
604                 nActiveADC++;           // number of 1 in mmm....m
605       }
606     }
607         x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC;   // nn = 01, ccccc are inverted, 0xc=1100
608         //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x );
609
610     if (nw < maxSize) {
611       buf[nw++] = x;
612       //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
613     }
614     else {
615       of++;
616     }
617   }
618
619   // Produce ADC data. 3 timebins are packed into one 32 bits word
620   // In this version, different ADC channel will NOT share the same word
621
622   UInt_t aa=0, a1=0, a2=0, a3=0;
623
624   for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
625     if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
626     aa = !(iAdc & 1) + 2;
627     for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
628       a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] >> fgkAddDigits : 0;
629       a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] >> fgkAddDigits : 0;
630       a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] >> fgkAddDigits : 0;
631       x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
632       if (nw < maxSize) {
633         buf[nw++] = x;
634         //printf("%08X ",x);
635       }
636       else {
637         of++;
638       }
639     }
640   }
641
642   if( of != 0 ) return -of; else return nw;
643 }
644
645 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
646 {
647   //
648   // Produce tracklet data stream from this MCM and put in buf
649   // Returns number of words filled, or negative value 
650   // with -1 * number of overflowed words
651   //
652
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   for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
662     if (nw < maxSize) 
663       buf[nw++] = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet])->GetTrackletWord();
664     else 
665       of++;
666   }
667   
668   if( of != 0 ) return -of; else return nw;
669 }
670
671 void AliTRDmcmSim::Filter()
672 {
673   //
674   // Filter the raw ADC values. The active filter stages and their
675   // parameters are taken from AliTRDtrapConfig.
676   // The raw data is stored separate from the filtered data. Thus, 
677   // it is possible to run the filters on a set of raw values 
678   // sequentially for parameter tuning.
679   //
680
681   if( !CheckInitialized() ) {
682     AliError("got called before initialization! Nothing done!");
683     return;
684   }
685
686   // Apply filters sequentially. Bypass is handled by filters
687   // since counters and internal registers may be updated even 
688   // if the filter is bypassed.
689   // The first filter takes the data from fADCR and 
690   // outputs to fADCF. 
691   
692   // Non-linearity filter not implemented.
693   FilterPedestal();
694   FilterGain();
695   FilterTail();
696   // Crosstalk filter not implemented.
697 }
698
699 void AliTRDmcmSim::FilterPedestalInit() 
700 {
701   // Initializes the pedestal filter assuming that the input has 
702   // been constant for a long time (compared to the time constant).
703
704 //  UShort_t    fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
705   UShort_t    fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
706   UShort_t    shifts[4] = {11, 14, 17, 21}; //??? where to take shifts from?
707
708   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++)
709     fPedAcc[iAdc] = (fSimParam->GetADCbaseline() << 2) * (1<<shifts[fptc]);
710 }
711
712 UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
713 {
714   // Returns the output of the pedestal filter given the input value.
715   // The output depends on the internal registers and, thus, the 
716   // history of the filter.
717
718   UShort_t    fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
719   UShort_t    fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
720   UShort_t    fpby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPBY); // 0..1 the bypass, active low
721   UShort_t    shifts[4] = {11, 14, 17, 21}; //??? where to come from
722
723   UShort_t accumulatorShifted;
724   Int_t correction;
725   UShort_t inpAdd;
726   
727   inpAdd = value + fpnp;
728
729   if (fpby == 0) //??? before or after update of accumulator
730     return value;
731
732   accumulatorShifted = (fPedAcc[adc] >> shifts[fptc]) & 0x3FF;   // 10 bits
733   if (timebin == 0) // the accumulator is disabled in the drift time
734   {
735     correction = (value & 0x3FF) - accumulatorShifted;
736     fPedAcc[adc] = (fPedAcc[adc] + correction) & 0x7FFFFFFF;             // 31 bits
737   }
738   
739   if (inpAdd <= accumulatorShifted)
740     return 0;
741   else
742   {
743     inpAdd = inpAdd - accumulatorShifted;
744     if (inpAdd > 0xFFF) 
745       return 0xFFF;
746     else 
747       return inpAdd;
748   }
749 }
750
751 void AliTRDmcmSim::FilterPedestal()
752 {
753   //
754   // Apply pedestal filter
755   //
756   // As the first filter in the chain it reads data from fADCR 
757   // and outputs to fADCF. 
758   // It has only an effect if previous samples have been fed to 
759   // find the pedestal. Currently, the simulation assumes that 
760   // the input has been stable for a sufficiently long time.
761
762   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
763     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
764       fADCF[iAdc][iTimeBin] = FilterPedestalNextSample(iAdc, iTimeBin, fADCR[iAdc][iTimeBin]);
765     }
766   }
767 }
768
769 void AliTRDmcmSim::FilterGainInit()
770 {
771   // Initializes the gain filter. In this case, only threshold 
772   // counters are reset.
773
774   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
775     // these are counters which in hardware continue 
776     // until maximum or reset
777     fGainCounterA[iAdc] = 0;
778     fGainCounterB[iAdc] = 0;
779   }
780 }
781
782 UShort_t AliTRDmcmSim::FilterGainNextSample(Int_t adc, UShort_t value)
783 {
784   // Apply the gain filter to the given value.
785   // BEGIN_LATEX O_{i}(t) = #gamma_{i} * I_{i}(t) + a_{i} END_LATEX
786   // The output depends on the internal registers and, thus, the 
787   // history of the filter.
788
789   UShort_t    fgby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGBY); // bypass, active low
790   UShort_t    fgf  = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGF0 + adc)); // 0x700 + (0 & 0x1ff);
791   UShort_t    fga  = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGA0 + adc)); // 40;
792   UShort_t    fgta = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTA); // 20;
793   UShort_t    fgtb = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTB); // 2060;
794
795   UInt_t tmp;
796
797   value &= 0xFFF;
798   tmp = (value * fgf) >> 11;
799   if (tmp > 0xFFF) tmp = 0xFFF;
800
801   if (fgby == 1)
802     value = AddUintClipping(tmp, fga, 12);
803
804   // Update threshold counters 
805   // not really useful as they are cleared with every new event
806   if ((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF))
807   {
808     if (value >= fgtb) 
809       fGainCounterB[adc]++;
810     else if (value >= fgta) 
811       fGainCounterA[adc]++;
812   }
813
814   return value;
815 }
816
817 void AliTRDmcmSim::FilterGain()
818 {
819   // Read data from fADCF and apply gain filter.
820
821   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
822     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
823         fADCF[iAdc][iTimeBin] = FilterGainNextSample(iAdc, fADCF[iAdc][iTimeBin]);
824     }
825   }
826 }
827
828 void AliTRDmcmSim::FilterTailInit(Int_t baseline)
829 {
830   // Initializes the tail filter assuming that the input has 
831   // been at the baseline value (configured by FTFP) for a 
832   // sufficiently long time.
833
834   // exponents and weight calculated from configuration
835   UShort_t    alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
836   UShort_t    lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
837   UShort_t    lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
838
839   Float_t lambdaL = lambdaLong  * 1.0 / (1 << 11);
840   Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
841   Float_t alphaL  = alphaLong   * 1.0 / (1 << 11);
842   Float_t qup, qdn;
843   qup = (1 - lambdaL) * (1 - lambdaS);
844   qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
845   Float_t kdc = qup/qdn;
846
847   Float_t kt, ql, qs;
848   UShort_t aout;
849   
850   kt = kdc * baseline;
851   aout = baseline - (UShort_t) kt;
852   ql = lambdaL * (1 - lambdaS) *      alphaL;
853   qs = lambdaS * (1 - lambdaL) * (1 - alphaL);
854
855   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
856     fTailAmplLong[iAdc]  = (UShort_t) (aout * ql / (ql + qs));
857     fTailAmplShort[iAdc] = (UShort_t) (aout * qs / (ql + qs));
858   }
859 }
860
861 UShort_t AliTRDmcmSim::FilterTailNextSample(Int_t adc, UShort_t value)
862 {
863   // Returns the output of the tail filter for the given input value. 
864   // The output depends on the internal registers and, thus, the 
865   // history of the filter.
866
867   // exponents and weight calculated from configuration
868   UShort_t    alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
869   UShort_t    lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
870   UShort_t    lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
871
872   Float_t lambdaL = lambdaLong  * 1.0 / (1 << 11);
873   Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
874   Float_t alphaL  = alphaLong   * 1.0 / (1 << 11);
875   Float_t qup, qdn;
876   qup = (1 - lambdaL) * (1 - lambdaS);
877   qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
878 //  Float_t kdc = qup/qdn;
879
880   UInt_t aDiff;
881   UInt_t alInpv;
882   UShort_t aQ;
883   UInt_t tmp;
884   
885   UShort_t inpVolt = value & 0xFFF;    // 12 bits
886       
887   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY) == 0) // bypass mode, active low
888     return value;
889   else
890   {   
891     // add the present generator outputs
892     aQ = AddUintClipping(fTailAmplLong[adc], fTailAmplShort[adc], 12);
893
894     // calculate the difference between the input the generated signal
895     if (inpVolt > aQ) 
896       aDiff = inpVolt - aQ;
897     else                
898       aDiff = 0;
899
900     // the inputs to the two generators, weighted
901     alInpv = (aDiff * alphaLong) >> 11;
902
903     // the new values of the registers, used next time
904     // long component
905     tmp = AddUintClipping(fTailAmplLong[adc], alInpv, 12);
906     tmp =  (tmp * lambdaLong) >> 11;
907     fTailAmplLong[adc] = tmp & 0xFFF;
908     // short component
909     tmp = AddUintClipping(fTailAmplShort[adc], aDiff - alInpv, 12);
910     tmp =  (tmp * lambdaShort) >> 11;
911     fTailAmplShort[adc] = tmp & 0xFFF;
912
913     // the output of the filter
914     return aDiff;
915   }
916 }
917
918 void AliTRDmcmSim::FilterTail()
919 {
920   // Apply tail filter
921
922   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
923     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
924       fADCF[iAdc][iTimeBin] = FilterTailNextSample(iAdc, fADCF[iAdc][iTimeBin]);
925     }
926   }
927 }
928
929 void AliTRDmcmSim::ZSMapping()
930 {
931   //
932   // Zero Suppression Mapping implemented in TRAP chip
933   //
934   // See detail TRAP manual "Data Indication" section:
935   // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
936   //
937
938   //??? values should come from TRAPconfig
939   Int_t eBIS = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIS); // TRAP default = 0x4  (Tis=4)
940   Int_t eBIT = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIT); // TRAP default = 0x28 (Tit=40)
941   Int_t eBIL = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIL); // TRAP default = 0xf0
942                                                                  // (lookup table accept (I2,I1,I0)=(111)
943                                                                  // or (110) or (101) or (100))
944   Int_t eBIN = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIN); // TRAP default = 1 (no neighbor sensitivity)
945   Int_t ep   = 0; // fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); //??? really subtracted here
946
947   Int_t **adc = fADCF;
948
949   if( !CheckInitialized() ) {
950     AliError("got called uninitialized! Nothing done!");    
951     return;
952   }
953
954   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
955     for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
956
957       // Get ADC data currently in filter buffer
958       Int_t ap = adc[iadc-1][it] - ep; // previous
959       Int_t ac = adc[iadc  ][it] - ep; // current
960       Int_t an = adc[iadc+1][it] - ep; // next
961
962       // evaluate three conditions
963       Int_t i0 = ( ac >=  ap && ac >=  an ) ? 0 : 1; // peak center detection
964       Int_t i1 = ( ap + ac + an > eBIT )    ? 0 : 1; // cluster
965       Int_t i2 = ( ac > eBIS )              ? 0 : 1; // absolute large peak
966
967       Int_t i = i2 * 4 + i1 * 2 + i0;    // Bit position in lookup table
968       Int_t d = (eBIL >> i) & 1;         // Looking up  (here d=0 means true
969                                          // and d=1 means false according to TRAP manual)
970
971       fZSM[iadc][it] &= d;
972       if( eBIN == 0 ) {  // turn on neighboring ADCs
973         fZSM[iadc-1][it] &= d;
974         fZSM[iadc+1][it] &= d;
975       }
976     }
977   }
978
979   // do 1 dim projection
980   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
981     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
982       fZSM1Dim[iadc] &= fZSM[iadc][it];
983     }
984   }
985 }
986
987 void AliTRDmcmSim::DumpData( char *f, char *target )
988 {
989   //
990   // Dump data stored (for debugging).
991   // target should contain one or multiple of the following characters
992   //   R   for raw data
993   //   F   for filtered data
994   //   Z   for zero suppression map
995   //   S   Raw dat astream
996   // other characters are simply ignored
997   //
998
999   UInt_t tempbuf[1024];
1000
1001   if( !CheckInitialized() ) return;
1002
1003   std::ofstream of( f, std::ios::out | std::ios::app );
1004   of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
1005              fDetector, fGeo->GetSector(fDetector), fGeo->GetStack(fDetector), 
1006              fGeo->GetSector(fDetector), fRobPos, fMcmPos );
1007
1008   for( int t=0 ; target[t] != 0 ; t++ ) {
1009     switch( target[t] ) {
1010     case 'R' :
1011     case 'r' :
1012       of << Form("fADCR (raw ADC data)\n");
1013       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1014         of << Form("  ADC %02d: ", iadc);
1015         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1016           of << Form("% 4d",  fADCR[iadc][it]);
1017         }
1018         of << Form("\n");
1019       }
1020       break;
1021     case 'F' :
1022     case 'f' :
1023       of << Form("fADCF (filtered ADC data)\n");
1024       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1025         of << Form("  ADC %02d: ", iadc);
1026         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1027           of << Form("% 4d",  fADCF[iadc][it]);
1028         }
1029         of << Form("\n");
1030       }
1031       break;
1032     case 'Z' :
1033     case 'z' :
1034       of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
1035       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1036         of << Form("  ADC %02d: ", iadc);
1037         if( fZSM1Dim[iadc] == 0 ) { of << " R   " ; } else { of << " .   "; } // R:read .:suppressed
1038         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1039           if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
1040         }
1041         of << Form("\n");
1042       }
1043       break;
1044     case 'S' :
1045     case 's' :
1046       Int_t s = ProduceRawStream( tempbuf, 1024 ); 
1047       of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
1048       of << Form("  address  data\n");
1049       for( int i = 0 ; i < s ; i++ ) {
1050         of << Form("  %04x     %08x\n", i, tempbuf[i]);
1051       }
1052     }
1053   }
1054 }
1055
1056 void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label) 
1057 {
1058   // Add the given hit to the fit register which is lateron used for 
1059   // the tracklet calculation. 
1060   // In addition to the fit sums in the fit register MC information 
1061   // is stored.
1062
1063   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)) && 
1064       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0)))
1065     fFitReg[adc].fQ0 += qtot;
1066   
1067   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) && 
1068       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)))
1069     fFitReg[adc].fQ1 += qtot;
1070   
1071   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) && 
1072       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE)))
1073   {
1074     fFitReg[adc].fSumX  += timebin;
1075     fFitReg[adc].fSumX2 += timebin*timebin;
1076     fFitReg[adc].fNhits++;
1077     fFitReg[adc].fSumY  += ypos;
1078     fFitReg[adc].fSumY2 += ypos*ypos;
1079     fFitReg[adc].fSumXY += timebin*ypos;
1080   }
1081
1082   // register hits (MC info)
1083   fHits[fNHits].fChannel = adc;
1084   fHits[fNHits].fQtot = qtot;
1085   fHits[fNHits].fYpos = ypos;
1086   fHits[fNHits].fTimebin = timebin;
1087   fHits[fNHits].fLabel = label;
1088   fNHits++;
1089 }
1090
1091 void AliTRDmcmSim::CalcFitreg() 
1092 {
1093   // Preprocessing.
1094   // Detect the hits and fill the fit registers.
1095   // Requires 12-bit data from fADCF which means Filter() 
1096   // has to be called before even if all filters are bypassed.
1097
1098   //???
1099   // TRAP parameters:
1100   const uint16_t lutPos[128] = {   // move later to some other file
1101     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,
1102     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,
1103     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,
1104     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};
1105   
1106   //??? to be clarified:
1107   UInt_t adcMask = 0xfffff;
1108   
1109   UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp;
1110   Short_t ypos, fromLeft, fromRight, found;
1111   UShort_t qTotal[19]; // the last is dummy
1112   UShort_t marked[6], qMarked[6], worse1, worse2;
1113   
1114   timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS); 
1115   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0) 
1116       < timebin1)
1117     timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0);
1118   timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE); 
1119   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1) 
1120       > timebin2)
1121     timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1);
1122
1123   // reset the fit registers
1124   fNHits = 0; 
1125   for (adcch = 0; adcch < fNADC-2; adcch++) // due to border channels
1126   {
1127     fFitReg[adcch].fNhits = 0;
1128     fFitReg[adcch].fQ0    = 0;
1129     fFitReg[adcch].fQ1    = 0;
1130     fFitReg[adcch].fSumX  = 0;
1131     fFitReg[adcch].fSumY  = 0;
1132     fFitReg[adcch].fSumX2 = 0;
1133     fFitReg[adcch].fSumY2 = 0;
1134     fFitReg[adcch].fSumXY = 0;
1135   }
1136   
1137   for (timebin = timebin1; timebin < timebin2; timebin++)
1138   {
1139     // first find the hit candidates and store the total cluster charge in qTotal array
1140     // in case of not hit store 0 there.
1141     for (adcch = 0; adcch < fNADC-2; adcch++) {
1142       if ( ( (adcMask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS
1143       {
1144         adcLeft  = fADCF[adcch  ][timebin];
1145         adcCentral  = fADCF[adcch+1][timebin];
1146         adcRight = fADCF[adcch+2][timebin];
1147         if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY) == 1) 
1148           hitQual = ( (adcLeft * adcRight) < 
1149                        (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT) * adcCentral) );
1150         else            
1151           hitQual = 1;
1152         // The accumulated charge is with the pedestal!!!
1153         qtotTemp = adcLeft + adcCentral + adcRight;
1154         if ( (hitQual) &&
1155              (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) &&
1156              (adcLeft <= adcCentral) &&
1157              (adcCentral > adcRight) )
1158           qTotal[adcch] = qtotTemp;
1159         else
1160           qTotal[adcch] = 0;
1161         //printf("ch %2d   qTotal %5d\n",adcch, qTotal[adcch]);
1162       }
1163       else
1164         qTotal[adcch] = 0; //jkl
1165     }
1166
1167     fromLeft = -1;
1168     adcch = 0;
1169     found = 0;
1170     marked[4] = 19; // invalid channel
1171     marked[5] = 19; // invalid channel
1172     qTotal[19] = 0;
1173     while ((adcch < 16) && (found < 3))
1174     {
1175       if (qTotal[adcch] > 0)
1176       {
1177         fromLeft = adcch;
1178         marked[2*found+1]=adcch;
1179         found++;
1180       }
1181       adcch++;
1182     }
1183     
1184     fromRight = -1;
1185     adcch = 18;
1186     found = 0;
1187     while ((adcch > 2) && (found < 3))
1188     {
1189       if (qTotal[adcch] > 0)
1190       {
1191         marked[2*found]=adcch;
1192         found++;
1193         fromRight = adcch;
1194       }
1195       adcch--;
1196     }
1197
1198     //printf("Fromleft=%d, Fromright=%d\n",fromLeft, fromRight);
1199     // here mask the hit candidates in the middle, if any
1200     if ((fromLeft >= 0) && (fromRight >= 0) && (fromLeft < fromRight))
1201       for (adcch = fromLeft+1; adcch < fromRight; adcch++)
1202         qTotal[adcch] = 0;
1203     
1204     found = 0;
1205     for (adcch = 0; adcch < 19; adcch++)
1206       if (qTotal[adcch] > 0) found++;
1207     // NOT READY
1208
1209     if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates!
1210     {
1211       if (marked[4] == marked[5]) marked[5] = 19;
1212       for (found=0; found<6; found++)
1213       {
1214         qMarked[found] = qTotal[marked[found]] >> 4;
1215         //printf("ch_%d qTotal %d qTotals %d |",marked[found],qTotal[marked[found]],qMarked[found]);
1216       }
1217       //printf("\n");
1218       
1219       Sort6To2Worst(marked[0], marked[3], marked[4], marked[1], marked[2], marked[5],
1220                     qMarked[0],
1221                     qMarked[3],
1222                     qMarked[4],
1223                     qMarked[1],
1224                     qMarked[2],
1225                     qMarked[5],
1226                     &worse1, &worse2);
1227       // Now mask the two channels with the smallest charge
1228       if (worse1 < 19)
1229       {
1230         qTotal[worse1] = 0;
1231         //printf("Kill ch %d\n",worse1);
1232       }
1233       if (worse2 < 19)
1234       {
1235         qTotal[worse2] = 0;
1236         //printf("Kill ch %d\n",worse2);
1237       }
1238     }
1239     
1240     for (adcch = 0; adcch < 19; adcch++) {
1241       if (qTotal[adcch] > 0) // the channel is marked for processing
1242       {
1243         adcLeft  = fADCF[adcch  ][timebin];
1244         adcCentral  = fADCF[adcch+1][timebin];
1245         adcRight = fADCF[adcch+2][timebin];
1246         // hit detected, in TRAP we have 4 units and a hit-selection, here we proceed all channels!
1247         // subtract the pedestal TPFP, clipping instead of wrapping
1248         
1249         Int_t regTPFP = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP);
1250 //        printf("Hit found, time=%d, adcch=%d/%d/%d, adc values=%d/%d/%d, regTPFP=%d, TPHT=%d\n",
1251 //               timebin, adcch, adcch+1, adcch+2, adcLeft, adcCentral, adcRight, regTPFP, 
1252 //               fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT));
1253
1254         if (adcLeft  < regTPFP) adcLeft  = 0; else adcLeft  -= regTPFP;
1255         if (adcCentral  < regTPFP) adcCentral  = 0; else adcCentral  -= regTPFP;
1256         if (adcRight < regTPFP) adcRight = 0; else adcRight -= regTPFP;
1257
1258         // Calculate the center of gravity
1259         // checking for adcCentral != 0 (in case of "bad" configuration)
1260         if (adcCentral == 0)
1261           continue;
1262         ypos = 128*(adcLeft - adcRight) / adcCentral;
1263         if (ypos < 0) ypos = -ypos;
1264         // make the correction using the LUT
1265         ypos = ypos + lutPos[ypos & 0x7F];
1266         if (adcLeft > adcRight) ypos = -ypos;
1267         AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, -1);
1268       }
1269     }
1270   }
1271 }
1272
1273 void AliTRDmcmSim::TrackletSelection() 
1274 {
1275   // Select up to 4 tracklet candidates from the fit registers  
1276   // and assign them to the CPUs.
1277
1278   UShort_t adcIdx, i, j, ntracks, tmp;
1279   UShort_t trackletCand[18][2]; // store the adcch[0] and number of hits[1] for all tracklet candidates
1280
1281   ntracks = 0;
1282   for (adcIdx = 0; adcIdx < 18; adcIdx++) // ADCs
1283     if ( (fFitReg[adcIdx].fNhits 
1284           >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL)) &&
1285          (fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits
1286           >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT)))
1287     {
1288       trackletCand[ntracks][0] = adcIdx;
1289       trackletCand[ntracks][1] = fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits;
1290       //printf("%d  %2d %4d\n", ntracks, trackletCand[ntracks][0], trackletCand[ntracks][1]);
1291       ntracks++;
1292     };
1293
1294   // for (i=0; i<ntracks;i++) printf("%d %d %d\n",i,trackletCand[i][0], trackletCand[i][1]);
1295
1296   if (ntracks > 4)
1297   {
1298     // primitive sorting according to the number of hits
1299     for (j = 0; j < (ntracks-1); j++)
1300     {
1301       for (i = j+1; i < ntracks; i++)
1302       {
1303         if ( (trackletCand[j][1]  < trackletCand[i][1]) ||
1304              ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) )
1305         {
1306           // swap j & i
1307           tmp = trackletCand[j][1];
1308           trackletCand[j][1] = trackletCand[i][1];
1309           trackletCand[i][1] = tmp;
1310           tmp = trackletCand[j][0];
1311           trackletCand[j][0] = trackletCand[i][0];
1312           trackletCand[i][0] = tmp;
1313         }
1314       }
1315     }
1316     ntracks = 4; // cut the rest, 4 is the max
1317   }
1318   // else is not necessary to sort
1319   
1320   // now sort, so that the first tracklet going to CPU0 corresponds to the highest adc channel - as in the TRAP
1321   for (j = 0; j < (ntracks-1); j++)
1322   {
1323     for (i = j+1; i < ntracks; i++)
1324     {
1325       if (trackletCand[j][0] < trackletCand[i][0])
1326       {
1327         // swap j & i
1328         tmp = trackletCand[j][1];
1329         trackletCand[j][1] = trackletCand[i][1];
1330         trackletCand[i][1] = tmp;
1331         tmp = trackletCand[j][0];
1332         trackletCand[j][0] = trackletCand[i][0];
1333         trackletCand[i][0] = tmp;
1334       }
1335     }
1336   }
1337   for (i = 0; i < ntracks; i++)  // CPUs with tracklets.
1338     fFitPtr[i] = trackletCand[i][0]; // pointer to the left channel with tracklet for CPU[i]
1339   for (i = ntracks; i < 4; i++)  // CPUs without tracklets
1340     fFitPtr[i] = 31;            // pointer to the left channel with tracklet for CPU[i] = 31 (invalid)
1341 //  printf("found %i tracklet candidates\n", ntracks);
1342 }
1343
1344 void AliTRDmcmSim::FitTracklet()
1345 {
1346   // Perform the actual tracklet fit based on the fit sums 
1347   // which have been filled in the fit registers. 
1348
1349   // parameters in fitred.asm (fit program)
1350   Int_t decPlaces = 5;
1351   Int_t rndAdd = 0;
1352   if (decPlaces >  1) 
1353     rndAdd = (1 << (decPlaces-1)) + 1;
1354   else if (decPlaces == 1)
1355     rndAdd = 1;
1356
1357   // should come from trapConfig (DMEM) 
1358   AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
1359   Long64_t shift = ((Long64_t) 1 << 32);
1360   UInt_t scaleY = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 160e-4)));
1361   UInt_t scaleD = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 140e-4)));
1362   int padrow = fFeeParam->GetPadRowFromMCM(fRobPos, fMcmPos);
1363   int yoffs  = (fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19) - fFeeParam->GetNcol()/2) << (8 + decPlaces); 
1364   int ndrift = 20; //??? value in simulation?
1365   int deflCorr = 0; // -370;
1366   int minslope = -10000; // no pt-cut so far
1367   int maxslope =  10000; // no pt-cut so far
1368
1369   // local variables for calculation
1370   Long64_t mult, temp, denom; //???
1371   UInt_t q0, q1, qTotal;          // charges in the two windows and total charge
1372   UShort_t nHits;                 // number of hits
1373   Int_t slope, offset;            // slope and offset of the tracklet
1374   Int_t sumX, sumY, sumXY, sumX2; // fit sums from fit registers
1375   //int32_t SumY2;                // not used in the current TRAP program
1376   FitReg_t *fit0, *fit1;          // pointers to relevant fit registers
1377   
1378 //  const uint32_t OneDivN[32] = {  // 2**31/N : exactly like in the TRAP, the simple division here gives the same result!
1379 //      0x00000000, 0x80000000, 0x40000000, 0x2AAAAAA0, 0x20000000, 0x19999990, 0x15555550, 0x12492490,
1380 //      0x10000000, 0x0E38E380, 0x0CCCCCC0, 0x0BA2E8B0, 0x0AAAAAA0, 0x09D89D80, 0x09249240, 0x08888880,
1381 //      0x08000000, 0x07878780, 0x071C71C0, 0x06BCA1A0, 0x06666660, 0x06186180, 0x05D17450, 0x0590B210,
1382 //      0x05555550, 0x051EB850, 0x04EC4EC0, 0x04BDA120, 0x04924920, 0x0469EE50, 0x04444440, 0x04210840};
1383
1384   for (Int_t cpu = 0; cpu < 4; cpu++) {
1385     if (fFitPtr[cpu] == 31)
1386     {
1387       fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker(); 
1388     }
1389     else
1390     {
1391       fit0 = &fFitReg[fFitPtr[cpu]  ];
1392       fit1 = &fFitReg[fFitPtr[cpu]+1]; // next channel
1393
1394       mult = 1;
1395       mult = mult << (32 + decPlaces);
1396       mult = -mult;
1397
1398       // Merging
1399       nHits   = fit0->fNhits + fit1->fNhits; // number of hits
1400       sumX    = fit0->fSumX  + fit1->fSumX;
1401       sumX2   = fit0->fSumX2 + fit1->fSumX2;
1402       denom   = nHits*sumX2 - sumX*sumX;
1403
1404       mult    = mult / denom; // exactly like in the TRAP program
1405       q0      = fit0->fQ0    + fit1->fQ0;
1406       q1      = fit0->fQ1    + fit1->fQ1;
1407       sumY    = fit0->fSumY  + fit1->fSumY  + 256*fit1->fNhits;
1408       sumXY   = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX;
1409
1410       slope   = nHits*sumXY - sumX * sumY;
1411       offset  = sumX2*sumY  - sumX * sumXY;
1412       temp    = mult * slope;
1413       slope   = temp >> 32; // take the upper 32 bits
1414       temp    = mult * offset;
1415       offset  = temp >> 32; // take the upper 32 bits
1416
1417       offset = offset + yoffs + (18 << (8 + decPlaces)); 
1418       slope  = slope * ndrift + deflCorr;
1419       offset = offset - (fFitPtr[cpu] << (8 + decPlaces));
1420       
1421       if ((slope < minslope) || (slope > maxslope))
1422       {
1423         fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1424       }
1425       else
1426       {
1427         temp    = slope;
1428         temp    = temp * scaleD;
1429         slope   = (temp >> 32);
1430         
1431         temp    = offset;
1432         temp    = temp * scaleY;
1433         offset  = (temp >> 32);
1434         
1435         // rounding, like in the TRAP
1436         slope   = (slope  + rndAdd) >> decPlaces;
1437         offset  = (offset + rndAdd) >> decPlaces;
1438
1439         if (slope > 0x3f || slope < -0x3f)
1440           AliWarning("Overflow in slope");
1441         slope   = slope  &   0x7F; // 7 bit
1442
1443         if (offset > 0xfff || offset < 0xfff)
1444           AliWarning("Overflow in offset");
1445         offset  = offset & 0x1FFF; // 13 bit
1446
1447         qTotal  = (q1 / nHits) >> 1;
1448         if (qTotal > 0xff)
1449           AliWarning("Overflow in charge");
1450         qTotal  = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program
1451
1452         // assemble and store the tracklet word
1453         fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
1454         new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
1455       }
1456     }
1457   }
1458 }
1459
1460 void AliTRDmcmSim::Tracklet()
1461 {
1462   // Run the tracklet calculation by calling sequentially:
1463   // CalcFitreg(); TrackletSelection(); FitTracklet()
1464   // and store the tracklets 
1465
1466   if (!fInitialized) {
1467     AliError("Called uninitialized! Nothing done!");
1468     return;
1469   }
1470
1471   fTrackletArray->Delete();
1472
1473   CalcFitreg();
1474   TrackletSelection();
1475   FitTracklet();
1476
1477   AliRunLoader *rl = AliRunLoader::Instance();
1478   AliDataLoader *dl = 0x0;
1479   if (rl)
1480     dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets");
1481   if (!dl) {
1482     AliError("Could not get the tracklets data loader!");
1483   }
1484   else {
1485     TTree *trackletTree = dl->Tree();
1486     if (!trackletTree) {
1487       dl->MakeTree();
1488       trackletTree = dl->Tree();
1489     }
1490
1491     AliTRDtrackletMCM *trkl = 0x0;
1492     TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
1493     if (!trkbranch)
1494       trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
1495
1496     for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
1497       trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]);
1498       trkbranch->SetAddress(&trkl);
1499 //      printf("filling tracklet 0x%08x\n", trkl->GetTrackletWord());
1500       trkbranch->Fill();
1501     }
1502     dl->WriteData("OVERWRITE");
1503   }
1504 }
1505
1506 void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
1507 {
1508   // write back the processed data configured by EBSF
1509   // EBSF = 1: unfiltered data; EBSF = 0: filtered data
1510   // zero-suppressed valued are written as -1 to digits
1511
1512   if (!fInitialized) {
1513     AliError("Called uninitialized! Nothing done!");
1514     return;
1515   }
1516
1517   Int_t firstAdc = 0;
1518   Int_t lastAdc = fNADC - 1;
1519
1520   while (GetCol(firstAdc) < 0)
1521     firstAdc++;
1522
1523   while (GetCol(lastAdc) < 0) 
1524     lastAdc--;
1525
1526   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
1527   {
1528     for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) {
1529       if (fZSM1Dim[iAdc] == 1) {
1530         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1531           digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, -1);
1532 //          printf("suppressed: %i, %i, %i, %i, now: %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin, 
1533 //                 digits->GetData(GetRow(), GetCol(iAdc), iTimeBin));
1534         }
1535       }
1536     }
1537   }
1538   else {
1539     for (Int_t iAdc = firstAdc; iAdc < lastAdc; iAdc++) {
1540       if (fZSM1Dim[iAdc] == 0) {
1541         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1542           digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, fADCF[iAdc][iTimeBin] >> fgkAddDigits);
1543         }
1544       }
1545       else {
1546         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1547           digits->SetData(GetRow(), GetCol(iAdc), iTimeBin, -1);
1548 //          printf("suppressed: %i, %i, %i, %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin);
1549         }
1550       }
1551     }
1552   }
1553 }
1554
1555 // help functions, to be cleaned up
1556
1557 UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
1558 {
1559   // 
1560   // This function adds a and b (unsigned) and clips to 
1561   // the specified number of bits. 
1562   //  
1563
1564   UInt_t sum = a + b;
1565   if (nbits < 32)
1566   {
1567     UInt_t maxv = (1 << nbits) - 1;;
1568     if (sum > maxv) 
1569       sum = maxv;
1570   }
1571   else
1572   {
1573     if ((sum < a) || (sum < b)) 
1574       sum = 0xFFFFFFFF;
1575   }
1576   return sum;
1577 }
1578
1579 void AliTRDmcmSim::Sort2(uint16_t  idx1i, uint16_t  idx2i, \
1580                             uint16_t  val1i, uint16_t  val2i, \
1581                             uint16_t *idx1o, uint16_t *idx2o, \
1582                             uint16_t *val1o, uint16_t *val2o) const
1583 {
1584   // sorting for tracklet selection
1585
1586     if (val1i > val2i)
1587     {
1588         *idx1o = idx1i;
1589         *idx2o = idx2i;
1590         *val1o = val1i;
1591         *val2o = val2i;
1592     }
1593     else
1594     {
1595         *idx1o = idx2i;
1596         *idx2o = idx1i;
1597         *val1o = val2i;
1598         *val2o = val1i;
1599     }
1600 }
1601
1602 void AliTRDmcmSim::Sort3(uint16_t  idx1i, uint16_t  idx2i, uint16_t  idx3i, \
1603                             uint16_t  val1i, uint16_t  val2i, uint16_t  val3i, \
1604                             uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, \
1605                             uint16_t *val1o, uint16_t *val2o, uint16_t *val3o)
1606 {
1607   // sorting for tracklet selection
1608
1609     int sel;
1610
1611
1612     if (val1i > val2i) sel=4; else sel=0;
1613     if (val2i > val3i) sel=sel + 2;
1614     if (val3i > val1i) sel=sel + 1;
1615     //printf("input channels %d %d %d, charges %d %d %d sel=%d\n",idx1i, idx2i, idx3i, val1i, val2i, val3i, sel);
1616     switch(sel)
1617     {
1618         case 6 : // 1 >  2  >  3            => 1 2 3
1619         case 0 : // 1 =  2  =  3            => 1 2 3 : in this case doesn't matter, but so is in hardware!
1620             *idx1o = idx1i;
1621             *idx2o = idx2i;
1622             *idx3o = idx3i;
1623             *val1o = val1i;
1624             *val2o = val2i;
1625             *val3o = val3i;
1626             break;
1627
1628         case 4 : // 1 >  2, 2 <= 3, 3 <= 1  => 1 3 2
1629             *idx1o = idx1i;
1630             *idx2o = idx3i;
1631             *idx3o = idx2i;
1632             *val1o = val1i;
1633             *val2o = val3i;
1634             *val3o = val2i;
1635             break;
1636
1637         case 2 : // 1 <= 2, 2 > 3, 3 <= 1   => 2 1 3
1638             *idx1o = idx2i;
1639             *idx2o = idx1i;
1640             *idx3o = idx3i;
1641             *val1o = val2i;
1642             *val2o = val1i;
1643             *val3o = val3i;
1644             break;
1645
1646         case 3 : // 1 <= 2, 2 > 3, 3  > 1   => 2 3 1
1647             *idx1o = idx2i;
1648             *idx2o = idx3i;
1649             *idx3o = idx1i;
1650             *val1o = val2i;
1651             *val2o = val3i;
1652             *val3o = val1i;
1653             break;
1654
1655         case 1 : // 1 <= 2, 2 <= 3, 3 > 1   => 3 2 1
1656             *idx1o = idx3i;
1657             *idx2o = idx2i;
1658             *idx3o = idx1i;
1659             *val1o = val3i;
1660             *val2o = val2i;
1661             *val3o = val1i;
1662         break;
1663
1664         case 5 : // 1 > 2, 2 <= 3, 3 >  1   => 3 1 2
1665             *idx1o = idx3i;
1666             *idx2o = idx1i;
1667             *idx3o = idx2i;
1668             *val1o = val3i;
1669             *val2o = val1i;
1670             *val3o = val2i;
1671         break;
1672
1673         default: // the rest should NEVER happen!
1674             printf("ERROR in Sort3!!!\n");
1675         break;
1676     }
1677 //    printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o);
1678 }
1679
1680 void AliTRDmcmSim::Sort6To4(uint16_t  idx1i, uint16_t  idx2i, uint16_t  idx3i, uint16_t  idx4i, uint16_t  idx5i, uint16_t  idx6i, \
1681                                uint16_t  val1i, uint16_t  val2i, uint16_t  val3i, uint16_t  val4i, uint16_t  val5i, uint16_t  val6i, \
1682                                uint16_t *idx1o, uint16_t *idx2o, uint16_t *idx3o, uint16_t *idx4o, \
1683                                uint16_t *val1o, uint16_t *val2o, uint16_t *val3o, uint16_t *val4o)
1684 {
1685   // sorting for tracklet selection
1686
1687     uint16_t idx21s, idx22s, idx23s, dummy;
1688     uint16_t val21s, val22s, val23s;
1689     uint16_t idx23as, idx23bs;
1690     uint16_t val23as, val23bs;
1691
1692     Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
1693                  idx1o, &idx21s, &idx23as,
1694                  val1o, &val21s, &val23as);
1695
1696     Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1697                  idx2o, &idx22s, &idx23bs,
1698                  val2o, &val22s, &val23bs);
1699
1700     Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, &dummy, &val23s, &dummy);
1701
1702     Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
1703                  idx3o, idx4o, &dummy,
1704                  val3o, val4o, &dummy);
1705
1706 }
1707
1708 void AliTRDmcmSim::Sort6To2Worst(uint16_t  idx1i, uint16_t  idx2i, uint16_t  idx3i, uint16_t  idx4i, uint16_t  idx5i, uint16_t  idx6i, \
1709                                     uint16_t  val1i, uint16_t  val2i, uint16_t  val3i, uint16_t  val4i, uint16_t  val5i, uint16_t  val6i, \
1710                                     uint16_t *idx5o, uint16_t *idx6o)
1711 {
1712   // sorting for tracklet selection
1713
1714     uint16_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5;
1715     uint16_t val21s, val22s, val23s;
1716     uint16_t idx23as, idx23bs;
1717     uint16_t val23as, val23bs;
1718
1719     Sort3(idx1i, idx2i,   idx3i, val1i, val2i, val3i,
1720                  &dummy1, &idx21s, &idx23as,
1721                  &dummy2, &val21s, &val23as);
1722
1723     Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
1724                  &dummy1, &idx22s, &idx23bs,
1725                  &dummy2, &val22s, &val23bs);
1726
1727     Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, idx5o, &val23s, &dummy1);
1728
1729     Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
1730                  &dummy1, &dummy2, idx6o,
1731                  &dummy3, &dummy4, &dummy5);
1732 //    printf("idx21s=%d, idx23as=%d, idx22s=%d, idx23bs=%d, idx5o=%d, idx6o=%d\n",
1733 //            idx21s,    idx23as,    idx22s,    idx23bs,    *idx5o,    *idx6o);
1734 }
1735
1736