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