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