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