First cleanup of Loaders
[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 //                                                                           //
20 ///////////////////////////////////////////////////////////////////////////////
21
22 /* $Id$ */
23
24 /*
25
26   New release on 2007/08/17
27
28 AliTRDmcmSim is now stably working and zero suppression function seems ok.
29 From now, the default version of raw data is set to 3 in AliTRDfeeParam.
30
31 The following internal parameters were abolished because it is useless and
32 made trouble:
33
34    fColOfADCbeg
35    fColOfADCend
36
37 GetCol member was modified accordingly. 
38
39 New member function DumpData was prepared for diagnostics.
40
41 ZSMapping member function was debugged. It was causing crash due to
42 wrong indexing in 1 dimensional numbering. Also code was shaped up better.
43
44 */
45
46 /*Semi-final version of TRD raw data simulation code with zero suppression (ZS)
47 similar to TRD FEE. ZS is realized by the class group:
48
49   AliTRDfeeParam
50   AliTRDmcmSim
51   AliTRDrawData
52
53 AliTRDfeeParam has been modified to have more parameters like raw data
54 production version and so on. AliTRDmcmSim is new class and this is the core
55 of MCM (PASA+TRAP) simulator. It has still very simple function and it will be
56 another project to improve this to make it closer to the reall FEE.
57 AliTRDrawData has been modified to use new class AliTRDmcmSim.
58
59 These modifications were tested on Aug. 02 HEAD version that code itself
60 compiles. I'm sure there must be still bugs and we need testing by as many as
61 possible persons now. Especially it seems HLT part is impacted by problems
62 because some parameters were moved from AliTRDrawData to AliTRDfeeParam (like
63 fRawVersion disappeared from AliTRDrawData).
64
65 In TRD definition, we have now 4 raw data versions.
66
67   0 very old offline version (by Bogdan)
68   1 test version (no zero suppression)
69   2 final version (no zero suppression)
70   3 test version (with zero suppression)
71
72 The default is still set to 2 in AliTRDfeeParam::fgkRAWversion and it uses
73 previously existing codes. If you set this to 3, AliTRDrawData changes behavior
74 to use AliTRDmcmSim with ZS.
75
76 Plan is after we make sure it works stably, we delete AliTRDmcm which is obsolete.
77 However it still take time because tracklet part is not yet touched.
78 The default raw version is 2.
79
80                                                                  Ken Oyama
81 */
82
83 // if no histo is drawn, these are obsolete
84 #include <TH1.h>
85 #include <TCanvas.h>
86
87 // only needed if I/O of tracklets is activated
88 #include <TObject.h>
89 #include <TFile.h>
90 #include <TTree.h>
91 #include <TSystem.h>
92
93 #include <fstream>
94
95 #include <TMath.h>
96
97 #include "AliLog.h"
98
99 #include "AliTRDmcmSim.h"
100 #include "AliTRDfeeParam.h"
101 #include "AliTRDSimParam.h"
102 #include "AliTRDgeometry.h"
103 #include "AliTRDcalibDB.h"
104 #include "AliTRDdigitsManager.h"
105 #include "AliTRDarrayADC.h"
106 // additional for new tail filter and/or tracklet
107 #include "AliTRDtrapAlu.h"
108 #include "AliTRDpadPlane.h"
109 #include "AliTRDtrackletMCM.h"
110
111 #include "AliRun.h"
112 #include "AliLoader.h"
113
114 ClassImp(AliTRDmcmSim)
115
116 //_____________________________________________________________________________
117 AliTRDmcmSim::AliTRDmcmSim() :TObject()
118   ,fInitialized(kFALSE)
119   ,fNextEvent(-1)    
120   ,fMaxTracklets(-1) 
121   ,fChaId(-1)
122   ,fSector(-1)
123   ,fStack(-1)
124   ,fLayer(-1)
125   ,fRobPos(-1)
126   ,fMcmPos(-1)
127   ,fNADC(-1)
128   ,fNTimeBin(-1)
129   ,fRow (-1)
130   ,fADCR(NULL)
131   ,fADCF(NULL)
132   ,fADCT(NULL)     
133   ,fPosLUT(NULL)    
134   ,fMCMT(NULL)      
135   ,fZSM(NULL)
136   ,fZSM1Dim(NULL)
137   ,fFeeParam(NULL)
138   ,fSimParam(NULL)
139   ,fCal(NULL)
140   ,fGeo(NULL)
141 {
142   //
143   // AliTRDmcmSim default constructor
144   //
145
146   // By default, nothing is initialized.
147   // It is necessary to issue Init before use.
148 }
149
150 //_____________________________________________________________________________
151 AliTRDmcmSim::AliTRDmcmSim(const AliTRDmcmSim &m) 
152   :TObject(m)
153   ,fInitialized(kFALSE) 
154   ,fNextEvent(-1)    
155   ,fMaxTracklets(-1) 
156   ,fChaId(-1)
157   ,fSector(-1)
158   ,fStack(-1)
159   ,fLayer(-1)
160   ,fRobPos(-1)
161   ,fMcmPos(-1)
162   ,fNADC(-1)
163   ,fNTimeBin(-1)
164   ,fRow(-1)
165   ,fADCR(NULL)
166   ,fADCF(NULL)
167   ,fADCT(NULL)      
168   ,fPosLUT(NULL)    
169   ,fMCMT(NULL)      
170   ,fZSM(NULL)
171   ,fZSM1Dim(NULL)
172   ,fFeeParam(NULL)
173   ,fSimParam(NULL)
174   ,fCal(NULL)
175   ,fGeo(NULL)
176   
177 {
178   //
179   // AliTRDmcmSim copy constructor
180   //
181
182   // By default, nothing is initialized.
183   // It is necessary to issue Init before use.
184 }
185
186 //_____________________________________________________________________________
187 AliTRDmcmSim::~AliTRDmcmSim() 
188 {
189   //
190   // AliTRDmcmSim destructor
191   //
192
193   if( fADCR != NULL ) {
194     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
195       delete [] fADCR[iadc];
196       delete [] fADCF[iadc];
197       delete [] fADCT[iadc];
198       delete [] fZSM [iadc];
199     }
200     delete [] fADCR;
201     delete [] fADCF;
202     delete [] fADCT;
203     delete [] fZSM;
204     delete [] fZSM1Dim;
205   }
206  
207   if(fInitialized){
208     delete [] fPosLUT;
209     delete [] fMCMT;
210   }
211
212   delete fGeo;
213
214 }
215
216 //_____________________________________________________________________________
217 AliTRDmcmSim &AliTRDmcmSim::operator=(const AliTRDmcmSim &m)
218 {
219   //
220   // Assignment operator
221   //
222
223   if (this != &m) {
224     ((AliTRDmcmSim &) m).Copy(*this);
225   }
226   return *this;
227
228 }
229
230 //_____________________________________________________________________________
231 void AliTRDmcmSim::Copy(TObject &m) const
232 {
233   //
234   // Copy function
235   //
236   ((AliTRDmcmSim &) m).fNextEvent     = 0; //new
237   ((AliTRDmcmSim &) m).fMaxTracklets  = 0; //new
238   ((AliTRDmcmSim &) m).fInitialized   = 0;
239   ((AliTRDmcmSim &) m).fChaId         = 0;
240   ((AliTRDmcmSim &) m).fSector        = 0;
241   ((AliTRDmcmSim &) m).fStack         = 0;
242   ((AliTRDmcmSim &) m).fLayer         = 0;
243   ((AliTRDmcmSim &) m).fRobPos        = 0;
244   ((AliTRDmcmSim &) m).fMcmPos        = 0;
245   ((AliTRDmcmSim &) m).fNADC          = 0;
246   ((AliTRDmcmSim &) m).fNTimeBin      = 0;
247   ((AliTRDmcmSim &) m).fRow           = 0;
248   ((AliTRDmcmSim &) m).fADCR          = 0;
249   ((AliTRDmcmSim &) m).fADCF          = 0;
250   ((AliTRDmcmSim &) m).fADCT          = 0; //new
251   ((AliTRDmcmSim &) m).fPosLUT        = 0; //new
252   ((AliTRDmcmSim &) m).fMCMT          = 0; //new
253   ((AliTRDmcmSim &) m).fZSM           = 0;
254   ((AliTRDmcmSim &) m).fZSM1Dim       = 0;
255   ((AliTRDmcmSim &) m).fFeeParam      = 0;
256   ((AliTRDmcmSim &) m).fSimParam      = 0;
257   ((AliTRDmcmSim &) m).fCal           = 0;
258   ((AliTRDmcmSim &) m).fGeo           = 0;
259
260 }
261
262 //_____________________________________________________________________________
263
264 //void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos ) 
265 void AliTRDmcmSim::Init( Int_t chaId, Int_t robPos, Int_t mcmPos, Bool_t newEvent = kFALSE ) // only for readout tree (new event)
266 {
267   //
268   // Initialize the class with new geometry information
269   // fADC array will be reused with filled by zero
270   //
271    
272   fNextEvent     = 0; 
273   fFeeParam      = AliTRDfeeParam::Instance();
274   fSimParam      = AliTRDSimParam::Instance();
275   fCal           = AliTRDcalibDB::Instance();
276   fGeo           = new AliTRDgeometry();
277   fChaId         = chaId;
278   fSector        = fGeo->GetSector( fChaId );
279   fStack         = fGeo->GetStack( fChaId );
280   fLayer         = fGeo->GetLayer( fChaId );
281   fRobPos        = robPos;
282   fMcmPos        = mcmPos;
283   fNADC          = fFeeParam->GetNadcMcm();
284   fNTimeBin      = fCal->GetNumberOfTimeBins();
285   fRow           = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
286
287   fMaxTracklets  = fFeeParam->GetMaxNrOfTracklets();
288
289  
290
291   
292   if (newEvent == kTRUE) {
293       fNextEvent = 1;
294   }
295
296
297
298   // Allocate ADC data memory if not yet done
299   if( fADCR == NULL ) {
300     fADCR    = new Int_t *[fNADC];
301     fADCF    = new Int_t *[fNADC];
302     fADCT    = new Int_t *[fNADC]; //new
303     fZSM     = new Int_t *[fNADC];
304     fZSM1Dim = new Int_t  [fNADC];
305     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
306       fADCR[iadc] = new Int_t[fNTimeBin];
307       fADCF[iadc] = new Int_t[fNTimeBin];
308       fADCT[iadc] = new Int_t[fNTimeBin]; //new
309       fZSM [iadc] = new Int_t[fNTimeBin];
310     }
311   }
312
313   // Initialize ADC data
314   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
315     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
316       fADCR[iadc][it] = 0;
317       fADCF[iadc][it] = 0;
318       fADCT[iadc][it] = -1;  //new
319       fZSM [iadc][it] = 1;   // Default unread = 1
320     }
321     fZSM1Dim[iadc] = 1;      // Default unread = 1
322   }
323   
324   //new:
325   fPosLUT = new Int_t[128];
326   for(Int_t i = 0; i<128; i++){
327     fPosLUT[i] = 0;
328   }
329   
330   fMCMT = new UInt_t[fMaxTracklets];
331   for(Int_t i = 0; i < fMaxTracklets; i++) {
332     fMCMT[i] = 0;
333   }
334
335
336   fInitialized = kTRUE;
337 }
338
339 //_____________________________________________________________________________
340 Bool_t AliTRDmcmSim::CheckInitialized()
341 {
342   //
343   // Check whether object is initialized
344   //
345
346   if( ! fInitialized ) {
347     AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
348   }
349   return fInitialized;
350 }
351
352 //_____________________________________________________________________________
353
354
355 void AliTRDmcmSim::SetPosLUT() {
356   Double_t iHi  = (Double_t)fCal->GetPRFhi();
357   Double_t iLo  = (Double_t)fCal->GetPRFlo();
358   Int_t   nBin  = fCal->GetPRFbin();
359   Int_t   iOff  = fLayer * nBin;
360   Int_t kNlayer = fGeo->Nlayer();
361
362   Float_t  *sPRFsmp   = new Float_t[nBin*kNlayer];
363   Double_t *sPRFlayer = new Double_t[nBin];
364   
365   
366   for(Int_t i = 0; i<nBin*kNlayer; i++){
367     
368     //printf("%f\n",fCal->GetSampledPRF()[i]);
369     sPRFsmp[i] = fCal->GetSampledPRF()[i]; 
370   
371   }
372
373   Double_t sWidth = (iHi-iLo)/((Double_t) nBin);
374   Int_t   sPad    = (Int_t) (1.0/sWidth);
375   
376   // get the PRF for actual layer (interpolated to ibin data-points; 61 measured)
377   for(Int_t iBin = 0; iBin < nBin; iBin++){
378     sPRFlayer[iBin] = (Double_t)sPRFsmp[iOff+iBin];
379   }
380
381   Int_t bin0 = (Int_t)(-iLo / sWidth - 0.5);                           // bin-nr. for pad-position 0
382   
383   Int_t bin1 = (Int_t)((Double_t)(0.5 - iLo) / sWidth - 0.5);          // bin-nr. for pad-position 0.5
384   bin1 = bin1 + 1;
385   bin0 = bin0 + 1;  //avoid negative values in aYest (start right of symmetry center)
386   while (bin0-sPad<0) {
387     bin0 = bin0 + 1;
388   }
389   while (bin1+sPad>=nBin) {
390     bin1 = bin1 - 1;
391   }
392   
393   Double_t* aYest = new Double_t[bin1-bin0+1];
394
395   /*TH1F* hist1 = new TH1F("h1","yest(y)",128,0,0.5);
396   TH1F* hist2 = new TH1F("h2","y(yest)",128,0,0.5);
397   TH1F* hist3 = new TH1F("h3","y(yest)-yest",128,0,0.5);
398   TH1F* hist4 = new TH1F("h4","y(yest)-yest,discrete",128,0,0.5);
399  
400   TCanvas *c1 = new TCanvas("c1","c1",800,1000);
401   hist1->Draw();
402   TCanvas *c2 = new TCanvas("c2","c2",800,1000);
403   hist2->Draw();
404   TCanvas *c3 = new TCanvas("c3","c3",800,1000);
405   hist3->Draw();
406   TCanvas *c4 = new TCanvas("c4","c4",800,1000);
407   hist4->Draw();*/
408   
409   for(Int_t iBin = bin0; iBin <= bin1; iBin++){
410     aYest[iBin-bin0] = 0.5*(sPRFlayer[iBin-sPad] - sPRFlayer[iBin+sPad])/(sPRFlayer[iBin]); // estimated position from PRF; between 0 and 1
411     //Double_t position = ((Double_t)(iBin)+0.5)*sWidth+iLo;
412     //  hist1->Fill(position,aYest[iBin-bin0]);
413   }
414   
415
416
417   Double_t aY[128]; // reversed function
418
419   AliTRDtrapAlu a;
420   a.Init(1,8,0,31);
421   
422   for(Int_t j = 0; j<128; j++) { // loop over all Yest; LUT has 128 entries; 
423     Double_t yest = ((Double_t)j)/256; 
424     
425     Int_t iBin = 0;
426     while (yest>aYest[iBin] && iBin<(bin1-bin0)) {
427       iBin = iBin+1;
428     }
429     if((iBin == bin1 - bin0)&&(yest>aYest[iBin])) {
430       aY[j] = 0.5;                      // yest too big
431       //hist2->Fill(yest,aY[j]);
432       
433     }
434     else {
435       Int_t bin_d = iBin + bin0 - 1;
436       Int_t bin_u = iBin + bin0;
437       Double_t y_d = ((Double_t)bin_d + 0.5)*sWidth + iLo; // lower y
438       Double_t y_u = ((Double_t)bin_u + 0.5)*sWidth + iLo; // upper y
439       Double_t yest_d = aYest[iBin-1];                     // lower estimated y
440       Double_t yest_u = aYest[iBin];                       // upper estimated y
441       
442       aY[j] = ((yest-yest_d)/(yest_u-yest_d))*(y_u-y_d) + y_d;
443       //hist2->Fill(yest,aY[j]);
444      
445     }
446     aY[j] = aY[j] - yest;
447     //hist3->Fill(yest,aY[j]);
448     // formatting
449     a.AssignDouble(aY[j]);
450     //a.WriteWord();
451     fPosLUT[j] = a.GetValue(); // 1+8Bit value;128 entries;LUT is steered by abs(Q(i+1)-Q(i-1))/Q(i)=COG and gives the correction to COG/2
452     //hist4->Fill(yest,fPosLUT[j]);
453     
454   }
455  
456    
457   
458   delete [] sPRFsmp;
459   delete [] sPRFlayer;
460   delete [] aYest;
461   
462 }
463
464
465 //_____________________________________________________________________________
466 Int_t* AliTRDmcmSim::GetPosLUT(){
467   return fPosLUT;
468 }
469
470
471
472 void AliTRDmcmSim::SetData( Int_t iadc, Int_t *adc )
473 {
474   //
475   // Store ADC data into array of raw data
476   //
477
478   if( !CheckInitialized() ) return;
479
480   if( iadc < 0 || iadc >= fNADC ) {
481     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
482     return;
483   }
484
485   for( int it = 0 ;  it < fNTimeBin ; it++ ) {
486     fADCR[iadc][it] = (Int_t)(adc[it]);
487   }
488 }
489
490 //_____________________________________________________________________________
491 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
492 {
493   //
494   // Store ADC data into array of raw data
495   //
496
497   if( !CheckInitialized() ) return;
498
499   if( iadc < 0 || iadc >= fNADC ) {
500     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
501     return;
502   }
503
504   fADCR[iadc][it] = adc;
505 }
506
507 //_____________________________________________________________________________
508 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
509 {
510   //
511   // Store ADC data into array of raw data
512   //
513
514   if( !CheckInitialized() ) return;
515
516   if( iadc < 0 || iadc >= fNADC ) {
517     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
518     return;
519   }
520
521   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
522     fADCR[iadc][it] = fSimParam->GetADCbaseline();
523   }
524 }
525
526 //_____________________________________________________________________________
527 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
528 {
529   //
530   // Return column id of the pad for the given ADC channel
531   //
532
533   if( !CheckInitialized() ) return -1;
534
535   return fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
536 }
537
538 //_____________________________________________________________________________
539 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize )
540 {
541   //
542   // Produce raw data stream from this MCM and put in buf
543   // Returns number of words filled, or negative value 
544   // with -1 * number of overflowed words
545   //
546
547   UInt_t  x;
548   UInt_t  iEv = 0;
549   Int_t   nw  = 0;  // Number of written words
550   Int_t   of  = 0;  // Number of overflowed words
551   Int_t   rawVer   = fFeeParam->GetRAWversion();
552   Int_t **adc;
553
554   if( !CheckInitialized() ) return 0;
555
556   if( fFeeParam->GetRAWstoreRaw() ) {
557     adc = fADCR;
558   } else {
559     adc = fADCF;
560   }
561
562   // Produce MCM header
563   x = ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
564   if (nw < maxSize) {
565     buf[nw++] = x;
566   }
567   else {
568     of++;
569   }
570
571   // Produce ADC mask
572   if( rawVer >= 3 ) {
573     x = 0;
574     for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
575       if( fZSM1Dim[iAdc] == 0 ) { //  0 means not suppressed
576         x = x | (1 << iAdc);
577       }
578     }
579     if (nw < maxSize) {
580       buf[nw++] = x;
581     }
582     else {
583       of++;
584     }
585   }
586
587   // Produce ADC data. 3 timebins are packed into one 32 bits word
588   // In this version, different ADC channel will NOT share the same word
589
590   UInt_t aa=0, a1=0, a2=0, a3=0;
591
592   for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
593     if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // suppressed
594     aa = !(iAdc & 1) + 2;
595     for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
596       a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] : 0;
597       a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
598       a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
599       x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
600       if (nw < maxSize) {
601         buf[nw++] = x;
602       }
603       else {
604         of++;
605       }
606     }
607   }
608
609   if( of != 0 ) return -of; else return nw;
610 }
611
612 //_____________________________________________________________________________
613 Int_t AliTRDmcmSim::ProduceRawStreamV2( UInt_t *buf, Int_t maxSize )
614 {
615   //
616   // Produce raw data stream from this MCM and put in buf
617   // Returns number of words filled, or negative value 
618   // with -1 * number of overflowed words
619   //
620
621   UInt_t  x;
622   UInt_t  iEv = 0;
623   Int_t   nw  = 0;  // Number of written words
624   Int_t   of  = 0;  // Number of overflowed words
625   Int_t   rawVer   = fFeeParam->GetRAWversion();
626   Int_t **adc;
627   Int_t   nActiveADC = 0;       // number of activated ADC bits in a word
628
629   if( !CheckInitialized() ) return 0;
630
631   if( fFeeParam->GetRAWstoreRaw() ) {
632     adc = fADCR;
633   } else {
634     adc = fADCF;
635   }
636
637   // Produce MCM header
638   x = (1<<31) | ((fRobPos * fFeeParam->GetNmcmRob() + fMcmPos) << 24) | ((iEv % 0x100000) << 4) | 0xC;
639   if (nw < maxSize) {
640     buf[nw++] = x;
641         //printf("\nMCM header: %X ",x);
642   }
643   else {
644     of++;
645   }
646
647   // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
648   //                            n : unused , c : ADC count, m : selected ADCs
649   if( rawVer >= 3 ) {
650     x = 0;
651     for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
652       if( fZSM1Dim[iAdc] == 0 ) { //  0 means not suppressed
653                 x = x | (1 << (iAdc+4) );       // last 4 digit reserved for 1100=0xc
654                 nActiveADC++;           // number of 1 in mmm....m
655       }
656     }
657         x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC;   // nn = 01, ccccc are inverted, 0xc=1100
658         //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x );
659
660     if (nw < maxSize) {
661       buf[nw++] = x;
662           //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
663     }
664     else {
665       of++;
666     }
667   }
668
669   // Produce ADC data. 3 timebins are packed into one 32 bits word
670   // In this version, different ADC channel will NOT share the same word
671
672   UInt_t aa=0, a1=0, a2=0, a3=0;
673
674   for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
675     if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
676     aa = !(iAdc & 1) + 2;
677     for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
678       a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] : 0;
679       a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] : 0;
680       a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] : 0;
681       x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
682       if (nw < maxSize) {
683         buf[nw++] = x;
684         //printf("%08X ",x);
685       }
686       else {
687         of++;
688       }
689     }
690   }
691
692   if( of != 0 ) return -of; else return nw;
693 }
694
695 //_____________________________________________________________________________
696 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
697 {
698   //
699   // Produce tracklet data stream from this MCM and put in buf
700   // Returns number of words filled, or negative value 
701   // with -1 * number of overflowed words
702   //
703
704   UInt_t  x;
705   Int_t   nw  = 0;  // Number of written words
706   Int_t   of  = 0;  // Number of overflowed words
707     
708   if( !CheckInitialized() ) return 0;
709
710   // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM 
711   // fMCMT is filled continuously until no more tracklet words available
712
713   Int_t wd = 0;
714   while ( (wd < fMaxTracklets) && (fMCMT[wd] > 0) ){
715       x = fMCMT[wd];
716       if (nw < maxSize) {
717         buf[nw++] = x;
718       }
719       else {
720         of++;
721       }
722       wd++;
723   }
724   
725   if( of != 0 ) return -of; else return nw;
726 }
727
728
729 //_____________________________________________________________________________
730 void AliTRDmcmSim::Filter()
731 {
732   //
733   // Apply digital filter
734   //
735
736   if( !CheckInitialized() ) return;
737
738   // Initialize filtered data array with raw data
739   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
740     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
741       fADCF[iadc][it] = fADCR[iadc][it]; 
742     }
743   }
744
745   // Then apply fileters one by one to filtered data array
746   if( fFeeParam->IsPFon() ) FilterPedestal();
747   if( fFeeParam->IsGFon() ) FilterGain();
748   if( fFeeParam->IsTFon() ) FilterTail();
749 }
750
751 //_____________________________________________________________________________
752 void AliTRDmcmSim::FilterPedestal()
753 {
754
755      
756   //
757   // Apply pedestal filter
758   //
759
760   Int_t ap = fSimParam->GetADCbaseline();      // ADC instrinsic pedestal
761   Int_t ep = fFeeParam->GetPFeffectPedestal(); // effective pedestal
762   //Int_t tc = fFeeParam->GetPFtimeConstant();   // this makes no sense yet
763
764   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
765     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
766       fADCF[iadc][it] = fADCF[iadc][it] - ap + ep;
767     }
768   }
769 }
770
771 //_____________________________________________________________________________
772 void AliTRDmcmSim::FilterGain()
773 {
774   //
775   // Apply gain filter (not implemented)
776   // Later it will be implemented because gain digital filiter will
777   // increase noise level.
778   //
779
780 }
781
782 //_____________________________________________________________________________
783 void AliTRDmcmSim::FilterTail()
784 {
785   //
786   // Apply exponential tail filter (Bogdan's version)
787   //
788
789   Double_t *dtarg  = new Double_t[fNTimeBin];
790   Int_t    *itarg  = new Int_t[fNTimeBin];
791   Int_t     nexp   = fFeeParam->GetTFnExp();
792   Int_t     tftype = fFeeParam->GetTFtype();
793
794   switch( tftype ) {
795     
796   case 0: // Exponential Filter Analog Bogdan
797     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
798       FilterSimDeConvExpA( fADCF[iCol], dtarg, fNTimeBin, nexp);
799       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
800         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
801       }
802     }
803     break;
804
805   case 1: // Exponential filter digital Bogdan
806     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
807       FilterSimDeConvExpD( fADCF[iCol], itarg, fNTimeBin, nexp);
808       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
809         fADCF[iCol][iTime] = itarg[iTime];
810       }
811     }
812     break;
813     
814   case 2: // Exponential filter Marian special
815     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
816       FilterSimDeConvExpMI( fADCF[iCol], dtarg, fNTimeBin);
817       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
818         fADCF[iCol][iTime] = (Int_t) TMath::Max(0.0,dtarg[iTime]);
819       }
820     }
821     break;
822
823     //new
824   case 3: // Exponential filter using AliTRDtrapAlu class
825     for (Int_t iCol = 0; iCol < fNADC; iCol++) {
826       FilterSimDeConvExpEl( fADCF[iCol], itarg, fNTimeBin, nexp);
827       for (Int_t iTime = 0; iTime < fNTimeBin; iTime++) {
828         fADCF[iCol][iTime] = itarg[iTime]>>2; // to be used for raw-data
829         fADCT[iCol][iTime] = itarg[iTime];    // 12bits; to be used for tracklet; tracklet will have own container; 
830       }
831     }
832     break;
833
834     
835   default:
836     AliError(Form("Invalid filter type %d ! \n", tftype ));
837     break;
838   }
839
840   delete [] dtarg;
841   delete [] itarg;
842
843 }
844
845 //_____________________________________________________________________________
846 void AliTRDmcmSim::ZSMapping()
847 {
848   //
849   // Zero Suppression Mapping implemented in TRAP chip
850   //
851   // See detail TRAP manual "Data Indication" section:
852   // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
853   //
854
855   Int_t eBIS = fFeeParam->GetEBsglIndThr();       // TRAP default = 0x4  (Tis=4)
856   Int_t eBIT = fFeeParam->GetEBsumIndThr();       // TRAP default = 0x28 (Tit=40)
857   Int_t eBIL = fFeeParam->GetEBindLUT();          // TRAP default = 0xf0
858                                                   // (lookup table accept (I2,I1,I0)=(111)
859                                                   // or (110) or (101) or (100))
860   Int_t eBIN = fFeeParam->GetEBignoreNeighbour(); // TRAP default = 1 (no neighbor sensitivity)
861   Int_t ep   = AliTRDfeeParam::GetPFeffectPedestal();
862
863   if( !CheckInitialized() ) return;
864
865   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
866     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
867
868       // Get ADC data currently in filter buffer
869       Int_t ap = fADCF[iadc-1][it] - ep; // previous
870       Int_t ac = fADCF[iadc  ][it] - ep; // current
871       Int_t an = fADCF[iadc+1][it] - ep; // next
872
873       // evaluate three conditions
874       Int_t i0 = ( ac >=  ap && ac >=  an ) ? 0 : 1; // peak center detection
875       Int_t i1 = ( ap + ac + an > eBIT )    ? 0 : 1; // cluster
876       Int_t i2 = ( ac > eBIS )              ? 0 : 1; // absolute large peak
877
878       Int_t i = i2 * 4 + i1 * 2 + i0;    // Bit position in lookup table
879       Int_t d = (eBIL >> i) & 1;         // Looking up  (here d=0 means true
880                                          // and d=1 means false according to TRAP manual)
881
882       fZSM[iadc][it] &= d;
883       if( eBIN == 0 ) {  // turn on neighboring ADCs
884         fZSM[iadc-1][it] &= d;
885         fZSM[iadc+1][it] &= d;
886       }
887
888     }
889   }
890
891   // do 1 dim projection
892   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
893     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
894       fZSM1Dim[iadc] &= fZSM[iadc][it];
895     }
896   }
897
898 }
899
900 //_____________________________________________________________________________
901 void AliTRDmcmSim::DumpData( char *f, char *target )
902 {
903   //
904   // Dump data stored (for debugging).
905   // target should contain one or multiple of the following characters
906   //   R   for raw data
907   //   F   for filtered data
908   //   Z   for zero suppression map
909   //   S   Raw dat astream
910   // other characters are simply ignored
911   //
912
913   UInt_t tempbuf[1024];
914
915   if( !CheckInitialized() ) return;
916
917   std::ofstream of( f, std::ios::out | std::ios::app );
918   of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
919              fChaId, fSector, fStack, fLayer, fRobPos, fMcmPos );
920
921   for( int t=0 ; target[t] != 0 ; t++ ) {
922     switch( target[t] ) {
923     case 'R' :
924     case 'r' :
925       of << Form("fADCR (raw ADC data)\n");
926       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
927         of << Form("  ADC %02d: ", iadc);
928         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
929           of << Form("% 4d",  fADCR[iadc][it]);
930         }
931         of << Form("\n");
932       }
933       break;
934     case 'F' :
935     case 'f' :
936       of << Form("fADCF (filtered ADC data)\n");
937       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
938         of << Form("  ADC %02d: ", iadc);
939         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
940           of << Form("% 4d",  fADCF[iadc][it]);
941         }
942         of << Form("\n");
943       }
944       break;
945     case 'Z' :
946     case 'z' :
947       of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
948       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
949         of << Form("  ADC %02d: ", iadc);
950         if( fZSM1Dim[iadc] == 0 ) { of << " R   " ; } else { of << " .   "; } // R:read .:suppressed
951         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
952           if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
953         }
954         of << Form("\n");
955       }
956       break;
957     case 'S' :
958     case 's' :
959       Int_t s = ProduceRawStream( tempbuf, 1024 ); 
960       of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
961       of << Form("  address  data\n");
962       for( int i = 0 ; i < s ; i++ ) {
963         of << Form("  %04x     %08x\n", i, tempbuf[i]);
964       }
965     }
966   }
967 }
968
969 //_____________________________________________________________________________
970 void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target
971                                      , Int_t n, Int_t nexp) 
972 {
973   //
974   // Exponential filter "analog"
975   // source will not be changed
976   //
977
978   Int_t    i = 0;
979   Int_t    k = 0;
980   Double_t reminder[2];
981   Double_t correction;
982   Double_t result;
983   Double_t rates[2];
984   Double_t coefficients[2];
985
986   // Initialize (coefficient = alpha, rates = lambda)
987   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
988
989   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
990   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
991   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
992   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
993   
994   coefficients[0] = c1;
995   coefficients[1] = c2;
996
997   Double_t dt = 0.1;
998   rates[0] = TMath::Exp(-dt/(r1));
999   rates[1] = TMath::Exp(-dt/(r2));
1000
1001   // Attention: computation order is important
1002   correction = 0.0;
1003   for (k = 0; k < nexp; k++) {
1004     reminder[k] = 0.0;
1005   }
1006     
1007   for (i = 0; i < n; i++) {
1008
1009     result    = ((Double_t)source[i] - correction);    // no rescaling
1010     target[i] = result;
1011     
1012     for (k = 0; k < nexp; k++) {
1013       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1014     }
1015       
1016     correction = 0.0;
1017     for (k = 0; k < nexp; k++) {
1018       correction += reminder[k];
1019     }
1020   }
1021 }
1022
1023 //_____________________________________________________________________________
1024 void AliTRDmcmSim::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n
1025                                      , Int_t nexp) 
1026 {
1027   //
1028   // Exponential filter "digital"
1029   // source will not be changed
1030   //
1031
1032   Int_t i        = 0;
1033   Int_t fAlphaL  = 0;
1034   Int_t fAlphaS  = 0;
1035   Int_t fTailPed = 0;
1036   Int_t iAlphaL  = 0;
1037   Int_t iAlphaS  = 0;
1038
1039   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
1040   // initialize (coefficient = alpha, rates = lambda)
1041
1042   Double_t dt = 0.1;
1043   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
1044   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
1045   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
1046   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
1047
1048   Int_t fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
1049   Int_t fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
1050   Int_t iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; //  9 bit paramter + fixed bits
1051   Int_t iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; //  9 bit paramter + fixed bits
1052
1053   if (nexp == 1) {
1054     fAlphaL = (Int_t) (c1 * 2048.0);
1055     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
1056   }
1057   if (nexp == 2) {
1058     fAlphaL = (Int_t) (c1 * 2048.0);
1059     fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
1060     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
1061     iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400;              // 10 bit paramter + fixed bits
1062   }
1063   
1064   Double_t iAl = iAlphaL  / 2048.0;            // alpha  L: correspondence to floating point numbers
1065   Double_t iAs = iAlphaS  / 2048.0;            // alpha  S: correspondence to floating point numbers
1066   Double_t iLl = iLambdaL / 2048.0;            // lambda L: correspondence to floating point numbers
1067   Double_t iLs = iLambdaS / 2048.0;            // lambda S: correspondence to floating point numbers
1068
1069   Int_t h1;
1070   Int_t h2;
1071   Int_t rem1;
1072   Int_t rem2;
1073   Int_t correction;
1074   Int_t result;
1075   Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2;
1076
1077   Double_t xi = 1 - (iLl*iAs + iLs*iAl);             // Calculation of equilibrium values of the
1078   rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
1079   rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
1080   
1081   // further initialization
1082   if ((rem1 + rem2) > 0x0FFF) {
1083     correction = 0x0FFF;
1084   } 
1085   else {
1086     correction = (rem1 + rem2) & 0x0FFF;
1087   }
1088
1089   fTailPed = iFactor - correction;
1090
1091   for (i = 0; i < n; i++) {
1092
1093     result = (source[i]  - correction);
1094     if (result < 0) { // Too much undershoot
1095       result = 0;
1096     }
1097
1098     target[i] = result;
1099                                                         
1100     h1 = (rem1 + ((iAlphaL * result) >> 11));
1101     if (h1 > 0x0FFF) {
1102       h1 = 0x0FFF;
1103     } 
1104     else {
1105       h1 &= 0x0FFF;
1106     }
1107
1108     h2 = (rem2 + ((iAlphaS * result) >> 11));
1109     if (h2 > 0x0FFF) {
1110       h2 = 0x0FFF;
1111     } 
1112     else {
1113       h2 &= 0x0FFF;
1114     }
1115   
1116     rem1 = (iLambdaL * h1 ) >> 11;
1117     rem2 = (iLambdaS * h2 ) >> 11;
1118     
1119     if ((rem1 + rem2) > 0x0FFF) {
1120       correction = 0x0FFF;
1121     } 
1122     else {
1123       correction = (rem1 + rem2) & 0x0FFF;
1124     }
1125
1126   }
1127
1128 }
1129
1130 //_____________________________________________________________________________
1131 void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target
1132                                       , Int_t n) 
1133 {
1134   //
1135   // Exponential filter (M. Ivanov)
1136   // source will not be changed
1137   //
1138
1139   Int_t i = 0;
1140   Double_t sig1[100];
1141   Double_t sig2[100];
1142   Double_t sig3[100];
1143
1144   for (i = 0; i < n; i++) {
1145     sig1[i] = (Double_t)source[i];
1146   }
1147
1148   Float_t dt      = 0.1;
1149   Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt;
1150   Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt;
1151
1152   FilterSimTailMakerSpline( sig1, sig2, lambda0, n);
1153   FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n);
1154
1155   for (i = 0; i < n; i++) {
1156     target[i] = sig3[i];
1157   }
1158
1159 }
1160
1161 //______________________________________________________________________________
1162 void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout
1163                                           , Double_t lambda, Int_t n) 
1164 {
1165   //
1166   // Special filter (M. Ivanov)
1167   //
1168
1169   Int_t    i = 0;
1170   Double_t l = TMath::Exp(-lambda*0.5);
1171   Double_t in[1000];
1172   Double_t out[1000];
1173
1174   // Initialize in[] and out[] goes 0 ... 2*n+19
1175   for (i = 0; i < n*2+20; i++) {
1176     in[i] = out[i] = 0;
1177   }
1178
1179   // in[] goes 0, 1
1180   in[0] = ampin[0];
1181   in[1] = (ampin[0] + ampin[1]) * 0.5;
1182    
1183   // Add charge to the end
1184   for (i = 0; i < 22; i++) {
1185     in[2*(n-1)+i] = ampin[n-1]; // in[] goes 2*n-2, 2*n-1, ... , 2*n+19 
1186   }
1187
1188   // Use arithmetic mean
1189   for (i = 1; i < n-1; i++) {
1190     in[2*i]   = ampin[i];    // in[] goes 2, 3, ... , 2*n-4, 2*n-3
1191     in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
1192   }
1193
1194   Double_t temp;
1195   out[2*n]    = in[2*n];
1196   temp        = 0;
1197   for (i = 2*n; i >= 0; i--) {
1198     out[i]    = in[i] + temp;
1199     temp      = l*(temp+in[i]);
1200   }
1201
1202   for (i = 0; i < n; i++){
1203     //ampout[i] = out[2*i+1];  // org
1204     ampout[i] = out[2*i];
1205   }
1206
1207 }
1208
1209 //______________________________________________________________________________
1210 void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout
1211                                             , Double_t norm, Double_t lambda
1212                                             , Int_t n) 
1213 {
1214   //
1215   // Special filter (M. Ivanov)
1216   //
1217
1218   Int_t    i = 0;
1219
1220   Double_t l = TMath::Exp(-lambda*0.5);
1221   Double_t k = l*(1.0 - norm*lambda*0.5);
1222   Double_t in[1000];
1223   Double_t out[1000];
1224
1225   // Initialize in[] and out[] goes 0 ... 2*n+19
1226   for (i = 0; i < n*2+20; i++) {
1227     in[i] = out[i] = 0;
1228   }
1229
1230   // in[] goes 0, 1
1231   in[0] = ampin[0];
1232   in[1] = (ampin[0]+ampin[1])*0.5;
1233
1234   // Add charge to the end
1235   for (i =-2; i < 22; i++) {
1236     // in[] goes 2*n-4, 2*n-3, ... , 2*n+19 
1237     in[2*(n-1)+i] = ampin[n-1];
1238   }
1239
1240   for (i = 1; i < n-2; i++) {
1241     // in[] goes 2, 3, ... , 2*n-6, 2*n-5
1242     in[2*i]    = ampin[i];
1243     in[2*i+1]  = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
1244     //in[2*i+1]  = ((ampin[i]+ampin[i+1]))/2.0;
1245   }
1246
1247   Double_t temp;
1248   out[0] = in[0];
1249   temp   = in[0];
1250   for (i = 1; i <= 2*n; i++) {
1251     out[i] = in[i] + (k-l)*temp;
1252     temp   = in[i] +  k   *temp;
1253   }
1254
1255   for (i = 0; i < n; i++) {
1256     //ampout[i] = out[2*i+1];  // org
1257     //ampout[i] = TMath::Max(out[2*i+1],0.0);  // org
1258     ampout[i] = TMath::Max(out[2*i],0.0);
1259   }
1260 }
1261
1262
1263 //_____________________________________________________________________________________
1264 //the following filter uses AliTRDtrapAlu-class
1265
1266 void AliTRDmcmSim::FilterSimDeConvExpEl(Int_t *source, Int_t *target, Int_t n, Int_t nexp) {
1267   //static Int_t count = 0;
1268  
1269   Double_t dt = 0.1;
1270   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
1271   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
1272   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
1273   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
1274   
1275   nexp = 1;
1276
1277   //it is assumed that r1,r2,c1,c2 are given such, that the configuration values are in the ranges according to TRAP-manual
1278   //parameters need to be adjusted
1279   AliTRDtrapAlu lambdaL;
1280   AliTRDtrapAlu lambdaS;
1281   AliTRDtrapAlu alphaL;
1282   AliTRDtrapAlu alphaS;
1283   
1284   AliTRDtrapAlu correction;
1285   AliTRDtrapAlu result;
1286   AliTRDtrapAlu bufL;
1287   AliTRDtrapAlu bufS;
1288  
1289   AliTRDtrapAlu bSource;
1290   
1291   lambdaL.Init(1,11);
1292   lambdaS.Init(1,11);
1293   alphaL.Init(1,11);
1294   alphaS.Init(1,11);
1295   
1296   //count=count+1;
1297
1298   lambdaL.AssignDouble(TMath::Exp(-dt/r1));
1299   lambdaS.AssignDouble(TMath::Exp(-dt/r2));
1300   alphaL.AssignDouble(c1); // in AliTRDfeeParam the number of exponentials is set and also the according time constants
1301   alphaS.AssignDouble(c2); // later it should be: alphaS=1-alphaL
1302   
1303   //data is enlarged to 12 bits, including 2 bits after the comma; class AliTRDtrapAlu is used to handle arithmetics correctly
1304   correction.Init(10,2);
1305   result.Init(10,2);
1306   bufL.Init(10,2);
1307   bufS.Init(10,2);
1308   bSource.Init(10,2);
1309   
1310   for(Int_t i = 0; i < n; i++) {
1311     bSource.AssignInt(source[i]);
1312     result = bSource - correction; // subtraction can produce an underflow
1313     if(result.GetSign() == kTRUE) {
1314       result.AssignInt(0);
1315     }
1316     
1317     //target[i] = result.GetValuePre();  // later, target and source should become AliTRDtrapAlu,too in order to simulate the 10+2Bits through the filter properly
1318     
1319     target[i] = result.GetValue(); // 12 bit-value; to get the corresponding integer value, target must be shifted: target>>2 
1320
1321     //printf("target-Wert zur Zeit %d : %d",i,target[i]);
1322     //printf("\n");
1323     
1324     bufL  =  bufL + (result * alphaL);
1325     bufL  =  bufL * lambdaL; 
1326     
1327     bufS  =  bufS + (result * alphaS);
1328     bufS  =  bufS * lambdaS;  // eventually this should look like:
1329     // bufS = (bufS + (result - result * alphaL)) * lambdaS // alphaS=1-alphaL; then alphaS-variable is not needed any more
1330
1331     correction = bufL + bufS; //check for overflow intrinsic; if overflowed, correction is set to 0x03FF
1332   }
1333  
1334   
1335 }
1336
1337
1338
1339
1340
1341
1342
1343 //__________________________________________________________________________________
1344
1345
1346 // in order to use the Tracklets, please first 
1347 // -- set AliTRDfeeParam::fgkTracklet to kTRUE, in order to switch on Tracklet-calculation
1348 // -- set AliTRDfeeParam::fgkTFtype   to 3, in order to use the new tail cancellation filter
1349 //    currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3
1350 // -- set AliTRDfeeParam::fgkMCTrackletOutput to kTRUE, if you want to use the Tracklet output container with           information about the Tracklet position (MCM, channel number)
1351
1352 // The code is designed such that the less possible calculations with AliTRDtrapAlu class-objects are performed; whenever possible calculations are done with doubles or integers and the results are transformed into the right format
1353
1354 void AliTRDmcmSim::Tracklet(){
1355     // tracklet calculation
1356     // if you use this code after a simulation, please make sure the same filter-settings as in the simulation are set in AliTRDfeeParam
1357
1358   if(!CheckInitialized()){ return; }
1359   
1360   Bool_t filtered = kTRUE;
1361   
1362   
1363   
1364   AliTRDtrapAlu data;
1365   data.Init(10,2);
1366   if(fADCT[0][0]==-1){                      // check if filter was applied
1367     filtered = kFALSE;
1368     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1369       for( Int_t iT = 0 ; iT < fNTimeBin ; iT++ ) {
1370         data.AssignInt(fADCR[iadc][iT]);
1371         fADCT[iadc][iT] = data.GetValue(); // all incoming values are positive 10+2 bit values; if el.filter was called, this is done correctly
1372       }
1373     }
1374    
1375   }
1376   
1377   // the online ordering of mcm's is reverse to the TRAP-manual-ordering! reverse fADCT (to be consistent to TRAP), then do all calculations
1378   // reverse fADCT:
1379   Int_t** rev0 = new Int_t *[fNADC];
1380   Int_t** rev1 = new Int_t *[fNADC];
1381   
1382   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1383     rev0[iadc] = new Int_t[fNTimeBin];
1384     rev1[iadc] = new Int_t[fNTimeBin];
1385     for( Int_t iT = 0; iT < fNTimeBin; iT++) {
1386       if( iadc <= fNADC-iadc-1 ) {
1387         rev0[iadc][iT]  = fADCT[fNADC-iadc-1][iT];
1388         rev1[iadc][iT]  = fADCT[iadc][iT];
1389         fADCT[iadc][iT] = rev0[iadc][iT];
1390       }
1391       else {
1392         rev0[iadc][iT]  = rev1[fNADC-iadc-1][iT];
1393         fADCT[iadc][iT] = rev0[iadc][iT];
1394       }
1395     }
1396   }
1397   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1398     delete[] rev0[iadc];
1399     delete[] rev1[iadc];
1400   }
1401   
1402   delete[] rev0;
1403   delete[] rev1;
1404   
1405   rev0 = NULL;
1406   rev1 = NULL;
1407     
1408   // get the filtered pedestal; supports only electronic tail-cancellation filter
1409   AliTRDtrapAlu filPed;
1410   Int_t ep = 0;
1411   Int_t *ieffped = new Int_t[fNTimeBin];
1412   for(Int_t iT = 0; iT < fNTimeBin; iT++){
1413     ieffped[iT] = ep; 
1414   }
1415  
1416   if( filtered == kTRUE ) {
1417     if( fFeeParam->IsPFon() ){
1418       ep = fFeeParam->GetPFeffectPedestal();
1419     }
1420     Int_t      nexp  = fFeeParam->GetTFnExp();
1421     Int_t  *isource  = new Int_t[fNTimeBin];
1422     filPed.Init(10,2);
1423     filPed.AssignInt(ep);           
1424     Int_t epf = filPed.GetValue();  
1425     for(Int_t iT = 0; iT < fNTimeBin; iT++){
1426       isource[iT] = ep;                  
1427       ieffped[iT] = epf;
1428     }
1429  
1430     if( fFeeParam->IsTFon() ) {
1431       FilterSimDeConvExpEl( isource, ieffped, fNTimeBin, nexp);
1432     }
1433   
1434     delete[] isource;
1435   }
1436   
1437   //the following values should go to AliTRDfeeParam once they are defined; then they have to be read in properly
1438   //naming follows conventions in TRAP-manual
1439   
1440   
1441   Bool_t bVBY = kTRUE;                         // cluster-verification bypass
1442
1443   Double_t cQTParam = 0;                      // cluster quality threshold; granularity 2^-10; range: 0<=cQT/2^-10<=2^-4 - 2^-10
1444   AliTRDtrapAlu cQTAlu; 
1445   cQTAlu.Init(1,10,0,63);
1446   cQTAlu.AssignDouble(cQTParam);
1447   Int_t cQT = cQTAlu.GetValue();
1448
1449   // linear fit 
1450   Int_t tFS = fFeeParam->GetLinearFitStart();  // linear fit start
1451   Int_t tFE = fFeeParam->GetLinearFitEnd();    // linear fit stop
1452    
1453   // charge accumulators
1454   Int_t tQS0 = fFeeParam->GetQacc0Start();     // start-time for charge-accumulator 0
1455   Int_t tQE0 = fFeeParam->GetQacc0End();       // stop-time for charge-accumulator 0
1456   Int_t tQS1 = fFeeParam->GetQacc1Start();     // start-time for charge-accumulator 1 
1457   Int_t tQE1 = fFeeParam->GetQacc1End();       // stop-time for charge-accumulator 1
1458   // values set such that tQS0=tFS; tQE0=tQS1-1; tFE=tQE1; want to do (QS0+QS1)/N
1459  
1460   Double_t cTHParam = (Double_t)fFeeParam->GetMinClusterCharge(); // cluster charge threshold
1461   AliTRDtrapAlu cTHAlu;  
1462   cTHAlu.Init(12,2);
1463   cTHAlu.AssignDouble(cTHParam);
1464   Int_t cTH = cTHAlu.GetValue();                                 // cTH used for comparison
1465
1466   struct List_t {
1467     List_t *next;
1468     Int_t iadc;
1469     Int_t value;
1470   };
1471   
1472   List_t selection[7];            // list with 7 elements
1473   List_t *list = NULL;
1474   List_t *listLeft = NULL;
1475     
1476   Int_t* qsum = new Int_t[fNADC];
1477    
1478   // fit sums
1479   AliTRDtrapAlu qsumAlu;
1480   qsumAlu.Init(12,2);           // charge sum will be 12+2 bits
1481   AliTRDtrapAlu dCOGAlu; 
1482   dCOGAlu.Init(1,7,0,127);      // COG will be 1+7 Bits; maximum 1 - 2^-7 for LUT
1483   AliTRDtrapAlu yrawAlu;
1484   yrawAlu.Init(1,8,-1,255);
1485   AliTRDtrapAlu yAlu;
1486   yAlu.Init(1,16,-1,0xFF00);    // only first 8 past-comma bits filled;additional 8 bits for accuracy;maximum 1 - 2^-8; sign is given by + or -
1487   AliTRDtrapAlu xAlu;
1488   xAlu.Init(5,8);               // 8 past-comma bits because value will be added/multiplied to another value with this accuracy
1489   AliTRDtrapAlu xxAlu;
1490   xxAlu.Init(10,0);            
1491   AliTRDtrapAlu yyAlu;
1492   yyAlu.Init(1,16,0,0xFFFF);    // maximum is 2^16-1; 16Bit for past-commas
1493   AliTRDtrapAlu xyAlu;
1494   xyAlu.Init(6,8);
1495   AliTRDtrapAlu XAlu;
1496   XAlu.Init(9,0);
1497   AliTRDtrapAlu XXAlu;
1498   XXAlu.Init(14,0);
1499   AliTRDtrapAlu YAlu;
1500   YAlu.Init(5,8);               // 14 bit, 1 is sign-bit; therefore only 13 bit 
1501   AliTRDtrapAlu YYAlu;
1502   YYAlu.Init(5,16);
1503   AliTRDtrapAlu XYAlu;
1504   XYAlu.Init(8,8);              // 17 bit, 1 is sign-bit; therefore only 16 bit        
1505   AliTRDtrapAlu qtruncAlu;
1506   qtruncAlu.Init(12,0);
1507   AliTRDtrapAlu QT0Alu;
1508   QT0Alu.Init(15,0);
1509   AliTRDtrapAlu QT1Alu;
1510   QT1Alu.Init(16,0);
1511
1512   AliTRDtrapAlu oneAlu;
1513   oneAlu.Init(1,8);
1514  
1515   
1516   AliTRDtrapAlu inverseNAlu;
1517   inverseNAlu.Init(1,8);        // simulates the LUT for 1/N
1518   AliTRDtrapAlu MeanChargeAlu;  // mean charge in ADC counts
1519   MeanChargeAlu.Init(8,0);
1520   AliTRDtrapAlu TotalChargeAlu;
1521   TotalChargeAlu.Init(17,8);
1522   //nr of post comma bits should be the same for inverseN and TotalCharge
1523   
1524   
1525   SetPosLUT();                    // initialize the position correction LUT for this MCM;
1526
1527
1528   // fit-sums; remapping!; 0,1,2->0; 1,2,3->1; ... 18,19,20->18
1529   Int_t *X   = new Int_t[fNADC-2];
1530   Int_t *XX  = new Int_t[fNADC-2];
1531   Int_t *Y   = new Int_t[fNADC-2];
1532   Int_t *YY  = new Int_t[fNADC-2];
1533   Int_t *XY  = new Int_t[fNADC-2];
1534   Int_t *N   = new Int_t[fNADC-2];
1535   Int_t *QT0 = new Int_t[fNADC-2]; // accumulated charge
1536   Int_t *QT1 = new Int_t[fNADC-2]; // accumulated charge
1537   
1538   for (Int_t iCol = 0; iCol < fNADC-2; iCol++) { 
1539       
1540       // initialize fit-sums 
1541       X[iCol]   = 0;
1542       XX[iCol]  = 0;
1543       Y[iCol]   = 0;
1544       YY[iCol]  = 0;
1545       XY[iCol]  = 0;
1546       N[iCol]   = 0;
1547       QT0[iCol] = 0;
1548       QT1[iCol] = 0;
1549   }
1550   
1551
1552   filPed.Init(7,2);                         // convert filtered pedestal into 7+2Bits
1553   
1554   for(Int_t iT = 0; iT < fNTimeBin; iT++){
1555     
1556     if(iT<tFS || iT>=tFE) continue;         // linear fit yes/no? 
1557
1558     // reset
1559     Int_t portChannel[4]   = {-1,-1,-1,-1};   
1560     Int_t clusterCharge[4] = {0,0,0,0};
1561     Int_t leftCharge[4]    = {0,0,0,0};
1562     Int_t centerCharge[4]  = {0,0,0,0}; 
1563     Int_t rightCharge[4]   = {0,0,0,0};
1564     
1565     Int_t mark = 0;
1566     
1567     filPed.AssignFormatted(ieffped[iT]);   // no size-checking when using AssignFormatted; ieffped>=0
1568     filPed = filPed;                       // this checks the size
1569     
1570     ieffped[iT] = filPed.GetValue();
1571         
1572     for(Int_t i = 0; i<7; i++){
1573       selection[i].next       = NULL;
1574       selection[i].iadc       =   -1;     // value of -1: invalid adc
1575       selection[i].value      =    0;
1576    
1577     }
1578     // selection[0] is starting list-element; just for pointing
1579
1580     // loop over inner adc's 
1581     for (Int_t iCol = 1; iCol < fNADC-1; iCol++) { 
1582       
1583       Int_t left   = fADCT[iCol-1][iT]; 
1584       Int_t center = fADCT[iCol][iT];
1585       Int_t right  = fADCT[iCol+1][iT];  
1586
1587       Int_t sum = left + center + right;            // cluster charge sum
1588       qsumAlu.AssignFormatted(sum);    
1589       qsumAlu = qsumAlu;                        // size-checking; redundant
1590  
1591       qsum[iCol] = qsumAlu.GetValue(); 
1592       
1593       //hit detection and masking
1594       if(center>=left){
1595         if(center>right){
1596           if(qsum[iCol]>=(cTH + 3*ieffped[iT])){    // effective pedestal of all three channels must be added to cTH(+20); this is not parallel to TRAP manual; maybe cTH has to be adjusted in fFeeParam; therefore channels are not yet reduced by their pedestal
1597             mark |= 1;                              // marker
1598           }
1599         }
1600       }
1601       mark = mark<<1;                
1602     }
1603     mark = mark>>1;
1604
1605        
1606     // get selection of 6 adc's and sort,starting with greatest values
1607
1608     //read three from right side and sort (primitive sorting algorithm)
1609     Int_t i = 0; // adc number
1610     Int_t j = 1; // selection number
1611     while(i<fNADC-2 && j<=3){
1612       i = i + 1;
1613       if( ((mark>>(i-1)) & 1) == 1) {
1614         selection[j].iadc  = fNADC-1-i;
1615         selection[j].value = qsum[fNADC-1-i]>>6;   // for hit-selection only the first 8 out of the 14 Bits are used for comparison
1616         
1617         // insert into sorted list
1618         listLeft = &selection[0];
1619         list = listLeft->next;
1620         
1621         if(list!=NULL) {
1622           while((list->next != NULL) && (selection[j].value <= list->value)){
1623             listLeft = list;
1624             list = list->next;
1625           }
1626           
1627           if(selection[j].value<=list->value){
1628             selection[j].next = list->next;
1629             list->next = &selection[j];
1630           }
1631           else {
1632             listLeft->next = &selection[j];
1633             selection[j].next = list;
1634           }
1635         }
1636         else{
1637           listLeft->next = &selection[j];
1638           selection[j].next = list;
1639         }
1640         
1641         j = j + 1;
1642       }
1643     }
1644
1645
1646     // read three from left side
1647     Int_t k = fNADC-2;
1648     while(k>i && j<=6) {
1649       if( ((mark>>(k-1)) & 1) == 1) {
1650         selection[j].iadc  = fNADC-1-k;
1651         selection[j].value = qsum[fNADC-1-k]>>6;
1652         
1653         listLeft = &selection[0];
1654         list = listLeft->next;
1655         
1656         if(list!=NULL){
1657           while((list->next != NULL) && (selection[j].value <= list->value)){
1658             listLeft = list;
1659             list = list->next;
1660           }
1661         
1662           if(selection[j].value<=list->value){
1663             selection[j].next = list->next;
1664             list->next = &selection[j];
1665           }
1666           else {
1667             listLeft->next = &selection[j];
1668             selection[j].next = list;
1669           }
1670         }
1671         else{
1672           listLeft->next = &selection[j];
1673           selection[j].next = list;
1674         }
1675
1676         j = j + 1;
1677       }
1678       k = k - 1;
1679     }
1680
1681     // get the four with greatest charge-sum
1682     list = &selection[0];
1683     for(i = 0; i<4; i++){
1684       if(list->next == NULL) continue;
1685       list = list->next;
1686       if(list->iadc == -1) continue;
1687       Int_t adc = list->iadc;                              // channel number with selected hit
1688       
1689       // the following arrays contain the four chosen channels in 1 time-bin
1690       portChannel[i]   = adc; 
1691       clusterCharge[i] = qsum[adc];
1692       leftCharge[i]    = fADCT[adc-1][iT] - ieffped[iT]; // reduce by filtered pedestal (pedestal is part of the signal)
1693       centerCharge[i]  = fADCT[adc][iT] - ieffped[iT];           
1694       rightCharge[i]   = fADCT[adc+1][iT] - ieffped[iT];         
1695     }
1696
1697     // arithmetic unit
1698     
1699     // cluster verification
1700     if(!bVBY){
1701       for(i = 0; i<4; i++){
1702         Int_t lr = leftCharge[i]*rightCharge[i]*1024;
1703         Int_t cc = centerCharge[i]*centerCharge[i]*cQT;
1704         if (lr>=cc){
1705           portChannel[i]   = -1;                                 // set to invalid address 
1706           clusterCharge[i] = 0;
1707         }
1708       }
1709     }
1710
1711     // fit-sums of valid channels
1712     // local hit position
1713     for(i = 0; i<4; i++){
1714       if (centerCharge[i] ==  0) {
1715         portChannel[i] = -1; 
1716       }// prevent division by 0
1717       
1718       if (portChannel[i]  == -1) continue;
1719       
1720       Double_t dCOG = (Double_t)(rightCharge[i]-leftCharge[i])/centerCharge[i];
1721        
1722       Bool_t sign = (dCOG>=0.0) ? kFALSE : kTRUE;
1723       dCOG = (sign == kFALSE) ? dCOG : -dCOG;     // AssignDouble doesn't allow for signed doubles
1724       dCOGAlu.AssignDouble(dCOG);
1725       Int_t iLUTpos = dCOGAlu.GetValue();       // steers position in LUT
1726             
1727       dCOG = dCOG/2;
1728       yrawAlu.AssignDouble(dCOG);
1729       Int_t iCOG = yrawAlu.GetValue();
1730       Int_t y = iCOG + fPosLUT[iLUTpos % 128];    // local position in pad-units
1731       yrawAlu.AssignFormatted(y);               // 0<y<1           
1732       yAlu  = yrawAlu;                        // convert to 16 past-comma bits
1733       
1734       if(sign == kTRUE) yAlu.SetSign(-1);       // buffer width of 9 bits; sign on real (not estimated) position
1735       xAlu.AssignInt(iT);                       // buffer width of 5 bits 
1736       
1737
1738       xxAlu = xAlu * xAlu;                  // buffer width of 10 bits -> fulfilled by x*x       
1739       
1740       yyAlu = yAlu * yAlu;                  // buffer width of 16 bits
1741    
1742       xyAlu = xAlu * yAlu;                  // buffer width of 14 bits
1743                   
1744       Int_t adc = portChannel[i]-1;              // remapping! port-channel contains channel-nr. of inner adc's (1..19; mapped to 0..18)
1745
1746       // calculate fit-sums recursively
1747       // interpretation of their bit-length is given as comment
1748       
1749       // be aware that the accuracy of the result of a calculation is always determined by the accuracy of the less accurate value
1750
1751       XAlu.AssignFormatted(X[adc]);
1752       XAlu = XAlu + xAlu;                   // buffer width of 9 bits 
1753       X[adc] = XAlu.GetValue();
1754              
1755       XXAlu.AssignFormatted(XX[adc]);
1756       XXAlu = XXAlu + xxAlu;                // buffer width of 14 bits    
1757       XX[adc] = XXAlu.GetValue();
1758
1759       if (Y[adc] < 0) {
1760         YAlu.AssignFormatted(-Y[adc]);          // make sure that only positive values are assigned; sign-setting must be done by hand
1761         YAlu.SetSign(-1);
1762       }
1763       else {
1764         YAlu.AssignFormatted(Y[adc]);
1765         YAlu.SetSign(1);
1766       }
1767         
1768       YAlu = YAlu + yAlu;                   // buffer width of 14 bits (8 past-comma);     
1769       Y[adc] = YAlu.GetSignedValue();
1770             
1771       YYAlu.AssignFormatted(YY[adc]);
1772       YYAlu = YYAlu + yyAlu;                // buffer width of 21 bits (16 past-comma) 
1773       YY[adc] = YYAlu.GetValue();
1774            
1775       if (XY[adc] < 0) {
1776         XYAlu.AssignFormatted(-XY[adc]);
1777         XYAlu.SetSign(-1);
1778       }
1779       else {
1780         XYAlu.AssignFormatted(XY[adc]);
1781         XYAlu.SetSign(1);
1782       }
1783
1784       XYAlu = XYAlu + xyAlu;                // buffer allows 17 bits (8 past-comma) 
1785       XY[adc] = XYAlu.GetSignedValue();
1786             
1787       N[adc]  = N[adc] + 1;
1788    
1789
1790       // accumulated charge
1791       qsumAlu.AssignFormatted(qsum[adc+1]); // qsum was not remapped!
1792       qtruncAlu = qsumAlu;
1793
1794       if(iT>=tQS0 && iT<=tQE0){
1795         QT0Alu.AssignFormatted(QT0[adc]);
1796         QT0Alu = QT0Alu + qtruncAlu;
1797         QT0[adc] = QT0Alu.GetValue();
1798         //interpretation of QT0 as 12bit-value (all pre-comma); is this as it should be done?; buffer allows 15 Bit
1799       }
1800       
1801       if(iT>=tQS1 && iT<=tQE1){
1802         QT1Alu.AssignFormatted(QT1[adc]);
1803         QT1Alu = QT1Alu + qtruncAlu;
1804         QT1[adc] = QT1Alu.GetValue();
1805         //interpretation of QT1 as 12bit-value; buffer allows 16 Bit
1806       }
1807     }// i
1808       
1809     // remapping is done!!
1810      
1811   }//iT
1812  
1813   
1814     
1815   // tracklet-assembly
1816   
1817   // put into AliTRDfeeParam and take care that values are in proper range
1818   const Int_t cTCL = 1;      // left adc: number of hits; 8<=TCL<=31 (?? 1<=cTCL<+8 ??) 
1819   const Int_t cTCT = 8;      // joint number of hits;     8<=TCT<=31; note that according to TRAP manual this number cannot be lower than 8; however it should be adjustable to the number of hits in the fit time range (40%)
1820   
1821   Int_t mPair   = 0;         // marker for possible tracklet pairs
1822   Int_t* hitSum = new Int_t[fNADC-3];
1823   // hitSum[0] means: hit sum of remapped channels 0 and 1; hitSum[17]: 17 and 18; 
1824   
1825   // check for all possible tracklet-pairs of adjacent channels (two are merged); mark the left channel of the chosen pairs
1826   for (Int_t iCol = 0; iCol < fNADC-3; iCol++) {
1827     hitSum[iCol] = N[iCol] + N[iCol+1];
1828     if ((N[iCol]>=cTCL) && (hitSum[iCol]>=cTCT)) {
1829         mPair |= 1;         // mark as possible channel-pair
1830      
1831     }
1832     mPair = mPair<<1;
1833   }
1834   mPair = mPair>>1;
1835   
1836   List_t* selectPair = new List_t[fNADC-2];      // list with 18 elements (0..18) containing the left channel-nr and hit sums
1837                                                  // selectPair[18] is starting list-element just for pointing
1838   for(Int_t k = 0; k<fNADC-2; k++){
1839       selectPair[k].next       = NULL;
1840       selectPair[k].iadc       =   -1;           // invalid adc
1841       selectPair[k].value      =    0;
1842    
1843     }
1844
1845  list = NULL;
1846  listLeft = NULL;
1847   
1848   // read marker and sort according to hit-sum
1849   
1850   Int_t adcL  = 0;            // left adc-channel-number (remapped)
1851   Int_t selNr = 0;            // current number in list
1852   
1853   // insert marked channels into list and sort according to hit-sum
1854   while(adcL < fNADC-3 && selNr < fNADC-3){
1855      
1856     if( ((mPair>>((fNADC-4)-(adcL))) & 1) == 1) {
1857       selectPair[selNr].iadc  = adcL;
1858       selectPair[selNr].value = hitSum[adcL];   
1859       
1860       listLeft = &selectPair[fNADC-3];
1861       list = listLeft->next;
1862         
1863       if(list!=NULL) {
1864         while((list->next != NULL) && (selectPair[selNr].value <= list->value)){
1865           listLeft = list;
1866           list = list->next;
1867         }
1868         
1869         if(selectPair[selNr].value <= list->value){
1870           selectPair[selNr].next = list->next;
1871           list->next = &selectPair[selNr];
1872         }
1873         else {
1874           listLeft->next = &selectPair[selNr];
1875           selectPair[selNr].next = list;
1876         }
1877         
1878       }
1879       else{
1880         listLeft->next = &selectPair[selNr];
1881         selectPair[selNr].next = list;
1882       }
1883       
1884       selNr = selNr + 1;
1885     }
1886     adcL = adcL + 1;
1887   }
1888   
1889   //select up to 4 channels with maximum number of hits
1890   Int_t lpairChannel[4] = {-1,-1,-1,-1}; // save the left channel-numbers of pairs with most hit-sum
1891   Int_t rpairChannel[4] = {-1,-1,-1,-1}; // save the right channel, too; needed for detecting double tracklets
1892   list = &selectPair[fNADC-3];
1893   
1894   for (Int_t i = 0; i<4; i++) {
1895     if(list->next == NULL) continue;
1896     list = list->next;
1897     if(list->iadc == -1) continue;
1898     lpairChannel[i] = list->iadc;        // channel number with selected hit
1899     rpairChannel[i] = lpairChannel[i]+1;
1900   }
1901   
1902   // avoid submission of double tracklets  
1903   for (Int_t i = 3; i>0; i--) {
1904     for (Int_t j = i-1; j>-1; j--) {
1905       if(lpairChannel[i] == rpairChannel[j]) {
1906         lpairChannel[i] = -1;
1907         rpairChannel[i] = -1;
1908         break;
1909       }
1910       /* if(rpairChannel[i] == lpairChannel[j]) {
1911         lpairChannel[i] = -1;
1912         rpairChannel[i] = -1;
1913         break;
1914         }*/
1915     }
1916   }
1917   
1918   // merging of the fit-sums of the remainig channels
1919   // assume same data-word-width as for fit-sums for 1 channel
1920   // relative scales!
1921   Int_t mADC[4];                      
1922   Int_t mN[4];
1923   Int_t mQT0[4];
1924   Int_t mQT1[4];
1925   Int_t mX[4];
1926   Int_t mXX[4];
1927   Int_t mY[4];
1928   Int_t mYY[4];
1929   Int_t mXY[4];
1930   Int_t mOffset[4];
1931   Int_t mSlope[4];
1932   Int_t mMeanCharge[4]; 
1933   Int_t inverseN = 0;
1934   Double_t invN = 0;
1935   Int_t one = 0;
1936
1937   for (Int_t i = 0; i<4; i++){
1938     mADC[i] = -1;                        // set to invalid number
1939     mN[i]   =  0;
1940     mQT0[i] =  0;
1941     mQT1[i] =  0;
1942     mX[i]   =  0;
1943     mXX[i]  =  0;
1944     mY[i]   =  0;
1945     mYY[i]  =  0;
1946     mXY[i]  =  0;
1947     mOffset[i] = 0;
1948     mSlope[i]  = 0;
1949     mMeanCharge[i] = 0;
1950   }
1951   
1952   oneAlu.AssignInt(1);
1953   one = oneAlu.GetValue();              // one with 8 past comma bits
1954  
1955   for (Int_t i = 0; i<4; i++){
1956           
1957
1958     mADC[i] = lpairChannel[i];          // mapping of merged sums to left channel nr. (0,1->0; 1,2->1; ... 17,18->17)
1959                                         // the adc and pad-mapping should now be one to one: adc i is linked to pad i; TRAP-numbering
1960     Int_t madc = mADC[i];
1961     if (madc == -1) continue;
1962     
1963     YAlu.AssignInt(N[rpairChannel[i]]);
1964     Int_t wpad  = YAlu.GetValue();       // enlarge hit counter of right channel by 8 past-comma bits; YAlu can have 5 pre-comma bits (values up to 63); hit counter<=nr of time bins (24)
1965
1966     mN[i]    = hitSum[madc];
1967   
1968     // don't merge fit sums in case of a stand-alone tracklet (consisting of only 1 channel); in that case only left channel makes up the fit sums
1969     if (N[madc+1] == 0) {
1970         mQT0[i] = QT0[madc];
1971         mQT1[i] = QT1[madc];
1972         
1973     }
1974     else {
1975
1976         // is it ok to do the size-checking for the merged fit-sums with the same format as for single-channel fit-sums?
1977         
1978         mQT0[i]   = QT0[madc] + QT0[madc+1];
1979         QT0Alu.AssignFormatted(mQT0[i]);   
1980         QT0Alu  = QT0Alu;                // size-check
1981         mQT0[i]   = QT0Alu.GetValue();     // write back
1982         
1983         mQT1[i]   = QT1[madc] + QT1[madc+1];
1984         QT1Alu.AssignFormatted(mQT1[i]);
1985         QT1Alu  = QT1Alu;
1986         mQT1[i]   = QT1Alu.GetValue();
1987     }
1988     
1989     // calculate the mean charge in adc values; later to be replaced by electron likelihood
1990     mMeanCharge[i] = mQT0[i] + mQT1[i]; // total charge
1991     mMeanCharge[i] = mMeanCharge[i]>>2; // losing of accuracy; accounts for high mean charge
1992     // simulate LUT for 1/N; LUT is fed with the double-accurate pre-calculated value of 1/N; accuracy of entries has to be adjusted to real TRAP
1993     invN = 1.0/(mN[i]);
1994     inverseNAlu.AssignDouble(invN);
1995     inverseN = inverseNAlu.GetValue();
1996     mMeanCharge[i] = mMeanCharge[i] * inverseN;  // now to be interpreted with 8 past-comma bits
1997     TotalChargeAlu.AssignFormatted(mMeanCharge[i]);
1998     TotalChargeAlu = TotalChargeAlu;
1999     MeanChargeAlu = TotalChargeAlu;
2000     mMeanCharge[i] = MeanChargeAlu.GetValue();
2001     
2002     // this check is not necessary; it is just for efficiency reasons
2003     if (N[madc+1] == 0) {
2004         mX[i]     =   X[madc];
2005         mXX[i]    =  XX[madc];
2006         mY[i]     =   Y[madc];
2007         mXY[i]    =  XY[madc];
2008         mYY[i]    =  YY[madc];
2009     }
2010     else {
2011         
2012         mX[i]     =   X[madc] +  X[madc+1];
2013         XAlu.AssignFormatted(mX[i]);
2014         XAlu      = XAlu;
2015         mX[i]     = XAlu.GetValue();
2016         
2017         mXX[i]    =  XX[madc] + XX[madc+1];
2018         XXAlu.AssignFormatted(mXX[i]);
2019         XXAlu     = XXAlu;
2020         mXX[i]    = XXAlu.GetValue();
2021  
2022     
2023         mY[i]     =   Y[madc] + Y[madc+1] + wpad;
2024         if (mY[i] < 0) {
2025             YAlu.AssignFormatted(-mY[i]);
2026             YAlu.SetSign(-1);
2027         }
2028         else {
2029             YAlu.AssignFormatted(mY[i]);
2030             YAlu.SetSign(1);
2031         }
2032         YAlu    = YAlu;
2033         mY[i]     = YAlu.GetSignedValue();
2034         
2035         mXY[i]    = XY[madc] + XY[madc+1] + X[madc+1]*one;    // multiplication by one to maintain the data format
2036         
2037         if (mXY[i] < 0) {
2038             XYAlu.AssignFormatted(-mXY[i]);
2039             XYAlu.SetSign(-1);
2040         }
2041         else {
2042             XYAlu.AssignFormatted(mXY[i]);
2043             XYAlu.SetSign(1);
2044         }
2045         XYAlu   = XYAlu;
2046         mXY[i]    = XYAlu.GetSignedValue();
2047     
2048         mYY[i]    = YY[madc] + YY[madc+1] + 2*Y[madc+1]*one+ wpad*one;
2049         if (mYY[i] < 0) {
2050             YYAlu.AssignFormatted(-mYY[i]);
2051             YYAlu.SetSign(-1);
2052         }
2053         else {
2054             YYAlu.AssignFormatted(mYY[i]);
2055             YYAlu.SetSign(1);
2056         }
2057         
2058         YYAlu   = YYAlu;
2059         mYY[i]    = YYAlu.GetSignedValue();
2060     }
2061   
2062   }
2063     
2064   // calculation of offset and slope from the merged fit-sums; 
2065   // YY is needed for some error measure only; still to be done
2066   // be aware that all values are relative values (scale: timebin-width; pad-width) and are integer values on special scale
2067   
2068   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2069   // !!important note: the offset is calculated from hits in the time bin range between tFS and tFE; it corresponds to the value at the height of the time bin tFS which does NOT need to correspond to the upper side of the drift   !!
2070   // !!volume (cathode wire plane). The offset cannot be rescaled as long as it is unknown which is the first time bin that contains hits from the drift region and thus to which distance from the cathode plane tFS corresponds.    !!
2071   // !!This has to be taken into account by the GTU. Furthermore a Lorentz correction might have to be applied to the offset (see below).                                                                                             !!
2072   // !!In this implementation it is assumed that no miscalibration containing changing drift velocities in the amplification region is used.                                                                                          !!
2073   // !!The corrections to the offset (e.g. no ExB correction applied as offset is supposed to be on top of drift region; however not at anode wire, so some inclination of drifting clusters due to Lorentz angle exists) are only    !!
2074   // !!valid (in approximation) if tFS is close to the beginning of the drift region.                                                                                                                                                 !!
2075   // !!The slope however can be converted to a deflection length between electrode and cathode wire plane as it is clear that the drift region is sampled 20 times                                                                    !!
2076   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2077
2078   // which formats should be chosen?
2079   AliTRDtrapAlu denomAlu;
2080   denomAlu.Init(20,8);       
2081   AliTRDtrapAlu numAlu;
2082   numAlu.Init(20,8);     
2083   // is this enough pre-comma place? covers the range of the 13 bit-word of the transmitted offset
2084   // offset measured in coord. of left channel must be between -0.5 and 1.5; 14 pre-comma bits because numerator can be big
2085
2086   for (Int_t i = 0; i<4; i++) {
2087     if (mADC[i] == -1) continue;
2088       
2089     Int_t num0  = (mN[i]*mXX[i]-mX[i]*mX[i]);
2090     if (num0 < 0) {
2091       denomAlu.AssignInt(-num0);    // num0 does not have to be interpreted as having past-comma bits -> AssignInt
2092       denomAlu.SetSign(-1);
2093     }
2094     else {
2095       denomAlu.AssignInt(num0);
2096       denomAlu.SetSign(1);
2097     }
2098     
2099     Int_t num1  = mN[i]*mXY[i] - mX[i]*mY[i];
2100     if (num1 < 0) {
2101       numAlu.AssignFormatted(-num1); // value of num1 is already formatted to have 8 past-comma bits
2102       numAlu.SetSign(-1);
2103     }
2104     else {
2105       numAlu.AssignFormatted(num1);
2106       numAlu.SetSign(1);
2107     }
2108     numAlu    = numAlu/denomAlu;
2109     mSlope[i]   = numAlu.GetSignedValue();
2110    
2111     Int_t num2  = mXX[i]*mY[i] - mX[i]*mXY[i];
2112    
2113     if (num2 < 0) {
2114       numAlu.AssignFormatted(-num2);
2115       numAlu.SetSign(-1);
2116     }
2117     else {
2118       numAlu.AssignFormatted(num2);
2119       numAlu.SetSign(1);
2120     }
2121    
2122     numAlu    = numAlu/denomAlu;
2123    
2124     
2125     mOffset[i]  = numAlu.GetSignedValue();
2126     numAlu.SetSign(1);
2127     denomAlu.SetSign(1);
2128        
2129                                  
2130     //numAlu.AssignInt(mADC[i]+1);   // according to TRAP-manual but trafo not to middle of chamber (0.5 channels away)             
2131     numAlu.AssignDouble((Double_t)mADC[i] + 1.5);      // numAlu has enough pre-comma place for that; correct trafo, best values
2132     mOffset[i]  = mOffset[i] + numAlu.GetValue();      // transform offset to a coord.system relative to chip; +1 to avoid neg. values 
2133     
2134     // up to here: adc-mapping according to TRAP manual and in line with pad-col mapping
2135     // reverse adc-counting to be again in line with the online mapping
2136     mADC[i]     = fNADC - 4 - mADC[i];                 // fNADC-4-mADC[i]: 0..17; remapping necessary;
2137     mADC[i]     = mADC[i] + 2; 
2138     // +2: mapping onto original ADC-online-counting: inner adc's corresponding to a chip's pasa: number 2..19
2139   }
2140
2141   // adc-counting is corresponding to online mapping; use AliTRDfeeParam::GetPadColFromADC to get the pad to which adc is connected; 
2142   // pad-column mapping is reverse to adc-online mapping; TRAP adc-mapping is in line with pad-mapping (increase in same direction);
2143   
2144   // transform parameters to the local coordinate-system of a stack (used by GTU)
2145   AliTRDpadPlane* padPlane = fGeo->CreatePadPlane(fLayer,fStack);
2146   
2147   Double_t padWidthI = padPlane->GetWidthIPad()*10.0; // get values in cm; want them in mm
2148   //Double_t padWidthO = padPlane->GetWidthOPad()*10; // difference between outer pad-widths not included; in real TRAP??
2149   
2150   // difference between width of inner and outer pads of a row is not accounted for;
2151   
2152   Double_t magField = 0.4;                           // z-component of magnetic field in Tesla; adjust to current simulation!!; magnetic field can hardly be evaluated for the position of each mcm 
2153   Double_t eCharge  = 0.3;                           // unit charge in (GeV/c)/m*T
2154   Double_t ptMin   = 2.3;                            // minimum transverse momentum (GeV/c); to be adjusted(?)
2155   
2156   Double_t granularityOffset = 0.160;                // granularity for offset in mm
2157   Double_t granularitySlope  = 0.140;                // granularity for slope  in mm     
2158     
2159   // get the coordinates in SM-system; parameters: 
2160   
2161   Double_t zPos       =  (padPlane->GetRowPos(fRow))*10.0;  // z-position of the MCM; fRow is counted on a chamber; SM consists of 5 
2162   // zPos is position of pad-borders;
2163   Double_t zOffset = 0.0;
2164   if ( fRow == 0 || fRow == 15 ) {
2165       zOffset = padPlane->GetLengthOPad();
2166   }
2167   else {
2168       zOffset = padPlane->GetLengthIPad();
2169   }
2170   zOffset = (-1.0) * zOffset/2.0;
2171   // turn zPos to be z-coordinate at middle of pad-row
2172   zPos = zPos + zOffset;
2173
2174       
2175   Double_t xPos       =  0.0;                               // x-position of the upper border of the drift-chamber of actual layer
2176   Int_t    icol       =  0;                                 // column-number of adc-channel
2177   Double_t yPos[4];                                         // y-position of the pad to which ADC is connected
2178   Double_t dx         = 30.0;                               // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry
2179   Double_t freqSample = fFeeParam->GetSamplingFrequency();  // retrieve the sampling frequency (10.019750 MHz)
2180   Double_t vdrift     = fCal->GetVdriftAverage(fChaId);     // averaged drift velocity for this detector (1.500000 cm/us)
2181   Int_t    nrOfDriftTimeBins = Int_t(dx/10.0*freqSample/vdrift); // the number of time bins in the drift region (20)
2182   Int_t    nrOfAmplTimeBins  = 2;                           // the number of time bins between anode wire and cathode wires in ampl.region (3.5mm)(guess)(suppose v_drift+3.5cm/us there=>all clusters arrive at anode wire within one time bin (100ns))
2183   Int_t    nrOfOffsetCorrTimeBins = tFS - nrOfAmplTimeBins - 1; // -1 is  to be conservative; offset correction will not remove the shift but is supposed to improve it; if tFS = 5, 2 drift time bins before tFS are assumed
2184   if(nrOfOffsetCorrTimeBins < 0) nrOfOffsetCorrTimeBins = 0;// don't apply offset correction if no drift time bins before tFS can be assumed 
2185   Double_t lorTan     = fCal->GetOmegaTau(vdrift,magField); // tan of the Lorentz-angle for this detector; could be evaluated and set as a parameter for each mcm
2186   //Double_t lorAngle   =  7.0;                             // Lorentz-angle in degrees
2187   Double_t tiltAngle  = padPlane->GetTiltingAngle();        // sign-respecting tilting angle of pads in actual layer
2188   Double_t tiltTan    = TMath::Tan(TMath::Pi()/180.0 * tiltAngle);
2189   //Double_t lorTan     = TMath::Tan(TMath::Pi()/180.0 * lorAngle);
2190
2191   Double_t alphaMax[4];                            // maximum deflection from the direction to the primary vertex; granularity of hit pads
2192   Double_t slopeMin[4];                            // local limits for the deflection
2193   Double_t slopeMax[4];
2194   Int_t   mslopeMin[4];                            // in granularity units; to be compared to mSlope[i]
2195   Int_t   mslopeMax[4];
2196
2197
2198   // x coord. of upper side of drift chambers in local SM-system (in mm)
2199   // obtained by evaluating the x-range of the hits; should be crosschecked; only drift, not amplification region taken into account (30mm);
2200   // the y-deflection is given as difference of y between lower and upper side of drift-chamber, not pad-plane;
2201   switch(fLayer) 
2202     {
2203     case 0: 
2204       xPos = 3003.0;
2205       break;
2206     case 1:
2207       xPos = 3129.0;
2208       break;
2209     case 2:
2210       xPos = 3255.0;
2211       break;
2212     case 3:
2213       xPos = 3381.0;
2214       break;
2215     case 4:
2216       xPos = 3507.0;
2217       break;
2218     case 5:
2219       xPos = 3633.0;
2220       break;
2221     }
2222  
2223   // calculation of offset-correction n: 
2224
2225   Int_t nCorrectOffset = (fRobPos % 2 == 0) ? ((fMcmPos % 4)) : ( 4 + (fMcmPos % 4));  
2226  
2227   nCorrectOffset = (nCorrectOffset - 4)*18 - 1;
2228   if (nCorrectOffset < 0) {
2229     numAlu.AssignInt(-nCorrectOffset);
2230     numAlu.SetSign(-1);
2231   }
2232   else {
2233     numAlu.AssignInt(nCorrectOffset);
2234     numAlu.SetSign(1);
2235   }
2236   nCorrectOffset = numAlu.GetSignedValue();   
2237
2238   // the Lorentz correction to the offset
2239   Double_t lorCorrectOffset = lorTan *(Double_t)nrOfOffsetCorrTimeBins*vdrift*10.0/freqSample; // Lorentz offset correction in mm
2240   
2241
2242   lorCorrectOffset = lorCorrectOffset/padWidthI; // Lorentz correction in pad width units
2243   
2244   if(lorCorrectOffset < 0) {
2245       numAlu.AssignDouble(-lorCorrectOffset);
2246       numAlu.SetSign(-1);
2247   }
2248   else{
2249       numAlu.AssignDouble(lorCorrectOffset);
2250       numAlu.SetSign(1);
2251   }
2252   
2253   Int_t mlorCorrectOffset = numAlu.GetSignedValue();
2254   
2255   
2256   Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0
2257  
2258   // calculation of slope-correction
2259
2260   // this is only true for tracks coming (approx.) from primary vertex
2261   // everything is evaluated for a tracklet covering the whole drift chamber
2262   Double_t cCorrectSlope = (-lorTan*dx + zPos/xPos*dx*tiltTan)/granularitySlope;
2263   // Double_t cCorrectSlope =  zPos/xPos*dx*tiltTan/granularitySlope;
2264   // zPos can be negative! for track from primary vertex: zOut-zIn > 0 <=> zPos > 0
2265   
2266   if (cCorrectSlope < 0) {
2267       numAlu.AssignDouble(-cCorrectSlope);
2268       numAlu.SetSign(-1);
2269   }
2270   else {
2271       numAlu.AssignDouble(cCorrectSlope);
2272       numAlu.SetSign(1);
2273   }
2274   cCorrectSlope = numAlu.GetSignedValue();
2275  
2276   // convert slope to deflection between upper and lower drift-chamber position (slope is given in pad-unit/time-bins)
2277   // different pad-width of outer pads of a pad-plane not taken into account
2278   // note that the fit was only done in the range tFS to tFE, however this range does not need to cover the whole drift region (neither start nor end of it)
2279   // however the tracklets are supposed to be a fit in the drift region thus the linear function is stretched to fit the drift region of 30 mm
2280   
2281   
2282   Double_t mCorrectSlope = (Double_t)(nrOfDriftTimeBins)*padWidthI/granularitySlope;  // >= 0.0
2283
2284   AliTRDtrapAlu correctAlu;
2285   correctAlu.Init(20,8);
2286   
2287   AliTRDtrapAlu offsetAlu;
2288   offsetAlu.Init(13,0,-0x1000,0x0FFF);          // 13 bit-word; 2-complement (1 sign-bit); asymmetric range
2289   
2290   AliTRDtrapAlu slopeAlu;
2291   slopeAlu.Init(7,0,-0x40,0x3F);                // 7 bit-word;  2-complement (1 sign-bit);
2292
2293   for (Int_t i = 0; i<4; i++) {
2294     
2295     if (mADC[i] == -1) continue;
2296     
2297     icol = fFeeParam->GetPadColFromADC(fRobPos,fMcmPos,mADC[i]); // be aware that mADC[i] contains the ADC-number according to online-mapping
2298     yPos[i]   = (padPlane->GetColPos(icol))*10.0;
2299     
2300     
2301     // offset:
2302     
2303     correctAlu.AssignDouble(mCorrectOffset);     // done because max. accuracy is 8 bit
2304     mCorrectOffset = correctAlu.GetValueWhole(); // cut offset correction to 8 past-comma bit accuracy
2305     mOffset[i]  = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset - mlorCorrectOffset)); 
2306     //mOffset[i]  = mOffset[i]*(-1);                   // adjust to direction of y-axes in online simulation
2307     
2308     if (mOffset[i] < 0) {
2309       numAlu.AssignFormatted(-mOffset[i]);
2310       numAlu.SetSign(-1);
2311     }
2312     else {
2313       numAlu.AssignFormatted(mOffset[i]);
2314       numAlu.SetSign(1);
2315     }
2316
2317     offsetAlu = numAlu; 
2318     mOffset[i]  = offsetAlu.GetSignedValue();  
2319
2320     
2321     // slope:
2322     
2323     correctAlu.AssignDouble(mCorrectSlope);
2324     mCorrectSlope = correctAlu.GetValueWhole();
2325     
2326     mSlope[i]   = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope);
2327
2328     if (mSlope[i] < 0) {
2329       numAlu.AssignFormatted(-mSlope[i]);
2330       numAlu.SetSign(-1);
2331     }
2332     else {
2333       numAlu.AssignFormatted(mSlope[i]);
2334       numAlu.SetSign(1);
2335     }
2336
2337     slopeAlu  = numAlu;     // here all past-comma values are cut, not rounded; alternatively add +0.5 before cutting (means rounding)
2338     mSlope[i]   = slopeAlu.GetSignedValue(); 
2339        
2340     // local (LTU) limits for the deflection 
2341     // ATan returns angles in radian
2342     alphaMax[i]  = TMath::ASin(eCharge*magField/(2.0*ptMin)*(TMath::Sqrt(xPos*xPos + yPos[i]*yPos[i]))/1000.0); // /1000: mm->m
2343     slopeMin[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) - alphaMax[i]))/granularitySlope;
2344     slopeMax[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) + alphaMax[i]))/granularitySlope;
2345     
2346     if (slopeMin[i] < 0) {
2347       slopeAlu.AssignDouble(-slopeMin[i]);
2348       slopeAlu.SetSign(-1);
2349     }
2350     else { 
2351       slopeAlu.AssignDouble(slopeMin[i]);
2352       slopeAlu.SetSign(1);
2353     }
2354     mslopeMin[i] = slopeAlu.GetSignedValue();  // the borders should lie inside the range of mSlope -> usage of slopeAlu again
2355    
2356     if (slopeMax[i] < 0) {
2357       slopeAlu.AssignDouble(-slopeMax[i]);
2358       slopeAlu.SetSign(-1);
2359     }
2360     else {
2361       slopeAlu.AssignDouble(slopeMax[i]);
2362       slopeAlu.SetSign(1);
2363     }
2364     mslopeMax[i] = slopeAlu.GetSignedValue();
2365   }
2366
2367   // suppress submission of tracks with low stiffness
2368   // put parameters in 32bit-word and submit (write to file as root-file; sort after SM, stack, layer, chamber) 
2369
2370   // sort tracklet-words in ascending y-order according to the offset (according to mADC would also be possible)
2371   // up to now they are sorted according to maximum hit sum
2372   // is the sorting really done in the TRAP-chip?
2373   
2374   Int_t order[4] = {-1,-1,-1,-1};
2375   Int_t wordnr = 0;   // number of tracklet-words
2376   
2377   for(Int_t j = 0; j < fMaxTracklets; j++) {
2378       //if( mADC[j] == -1) continue; 
2379       if( (mADC[j] == -1) || (mSlope[j] < mslopeMin[j]) || (mSlope[j] > mslopeMax[j])) continue; // this applies a pt-cut
2380       wordnr++;
2381       if( wordnr-1 == 0) {
2382           order[0] = j;
2383           continue;
2384       }
2385       // wordnr-1>0, wordnr-1<4
2386       order[wordnr-1] = j;
2387       for( Int_t k = 0; k < wordnr-1; k++) {
2388           if( mOffset[j] < mOffset[order[k]] ) {
2389               for( Int_t l = wordnr-1; l > k; l-- ) {
2390                   order[l] = order[l-1];
2391               }
2392               order[k] = j;
2393               break;
2394           }
2395           
2396       }
2397   }
2398         
2399   // fill the bit-words in ascending order and without gaps
2400   UInt_t bitWord[4] = {0,0,0,0};                 // attention: unsigned int to have real 32 bits (no 2-complement)
2401   for(Int_t j = 0; j < wordnr; j++) { // only "wordnr" tracklet-words
2402       //Bool_t rem1 = kTRUE;
2403     
2404     Int_t i = order[j];
2405     //bit-word is 2-complement and therefore without sign
2406     bitWord[j] =   1; // this is the starting 1 of the bit-word (at 33rd position); the 1 must be ignored
2407     //printf("\n");
2408     UInt_t shift  = 0;
2409     UInt_t shift2 = 0;
2410         
2411         
2412
2413
2414     /*printf("mean charge: %d\n",mMeanCharge[i]);
2415     printf("row: %d\n",fRow);
2416     printf("slope: %d\n",mSlope[i]);
2417     printf("pad position: %d\n",mOffset[i]);
2418     printf("channel: %d\n",mADC[i]);*/
2419
2420     // electron probability (currently not implemented; the mean charge is just scaled)
2421     shift = (UInt_t)mMeanCharge[i];
2422     for(Int_t iBit = 0; iBit < 8; iBit++) {
2423       bitWord[j]  = bitWord[j]<<1;
2424       bitWord[j] |= (shift>>(7-iBit))&1;               
2425       //printf("0");
2426     }
2427
2428     // pad row
2429     shift = (UInt_t)fRow;
2430     for(Int_t iBit = 0; iBit < 4; iBit++) {
2431       bitWord[j]  = bitWord[j]<<1;
2432       bitWord[j] |= (shift>>(3-iBit))&1;
2433       //printf("%d", (fRow>>(3-iBit))&1);
2434     }
2435     
2436     // deflection length
2437     if(mSlope[i] < 0) {
2438         shift = (UInt_t)(-mSlope[i]);
2439         // shift2 is 2-complement of shift
2440         shift2 = 1;
2441         for(Int_t iBit = 1; iBit < 7; iBit++) {
2442             shift2  = shift2<<1;
2443             shift2 |= (1- (((shift)>>(6-iBit))&1) );
2444             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2445         }
2446         shift2 = shift2 + 1;
2447         //printf("1");
2448         for(Int_t iBit = 0; iBit < 7; iBit++) {
2449             bitWord[j]  = bitWord[j]<<1;
2450             bitWord[j] |= (shift2>>(6-iBit))&1;
2451             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2452         }
2453     }
2454     else {
2455         shift = (UInt_t)(mSlope[i]);
2456         bitWord[j]  = bitWord[j]<<1;
2457         bitWord[j]   |= 0;
2458         //printf("0");
2459         for(Int_t iBit = 1; iBit < 7; iBit++) {
2460             bitWord[j]  = bitWord[j]<<1;
2461             bitWord[j] |= (shift>>(6-iBit))&1;
2462             //printf("%d",(mSlope[i]>>(6-iBit))&1);
2463         }
2464     }
2465
2466     // pad position
2467     if(mOffset[i] < 0) {
2468         shift = (UInt_t)(-mOffset[i]);
2469         shift2 = 1;
2470         for(Int_t iBit = 1; iBit < 13; iBit++) {
2471             shift2  = shift2<<1;
2472             shift2 |= (1-(((shift)>>(12-iBit))&1));
2473             //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
2474         }
2475         shift2 = shift2 + 1;
2476         //printf("1");
2477         for(Int_t iBit = 0; iBit < 13; iBit++) {
2478             bitWord[j]  = bitWord[j]<<1;
2479             bitWord[j] |= (shift2>>(12-iBit))&1;
2480             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2481         }
2482     }
2483     else {
2484         shift = (UInt_t)mOffset[i];
2485         bitWord[j] = bitWord[j]<<1;
2486         bitWord[j]   |= 0; 
2487         //printf("0");
2488         for(Int_t iBit = 1; iBit < 13; iBit++) {
2489             bitWord[j]  = bitWord[j]<<1;
2490             bitWord[j] |= (shift>>(12-iBit))&1;
2491             //printf("%d",(mOffset[i]>>(12-iBit))&1);
2492         }
2493     }
2494
2495
2496         
2497     //printf("bitWord: %u\n",bitWord[j]);
2498     //printf("adc: %d\n",mADC[i]);
2499     fMCMT[j] = bitWord[j];
2500   }
2501     
2502   //printf("\n");
2503
2504   
2505   delete [] qsum;
2506   delete [] ieffped;
2507
2508   delete [] X;
2509   delete [] XX;
2510   delete [] Y;
2511   delete [] YY;
2512   delete [] XY;
2513   delete [] N;
2514   delete [] QT0;
2515   delete [] QT1;
2516
2517   delete [] hitSum;
2518   delete [] selectPair;
2519
2520   delete padPlane;
2521
2522 //if you want to activate the MC tracklet output, set fgkMCTrackletOutput=kTRUE in AliTRDfeeParam
2523         
2524   if (!fFeeParam->GetMCTrackletOutput()) 
2525       return;
2526  
2527   AliLog::SetClassDebugLevel("AliTRDmcmSim", 10);
2528   AliLog::SetFileOutput("../log/tracklet.log");
2529   
2530   // testing for wordnr in order to speed up the simulation
2531   if (wordnr == 0) 
2532     return;
2533    
2534   UInt_t        *trackletWord = new UInt_t[fMaxTracklets];
2535   Int_t         *adcChannel   = new Int_t[fMaxTracklets];
2536   Int_t         *trackRef     = new Int_t[fMaxTracklets];
2537
2538   Int_t u = 0;
2539
2540   AliTRDdigitsManager *digman = new AliTRDdigitsManager();
2541   digman->ReadDigits(AliRunLoader::GetRunLoader()->GetLoader("TRDLoader")->TreeD());
2542   digman->SetUseDictionaries(kTRUE);
2543   AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
2544
2545   for (Int_t j = 0; j < fMaxTracklets; j++) {
2546       Int_t i = order[j];
2547       trackletWord[j] = 0;
2548       adcChannel[j] = -1;
2549       if (bitWord[j]!=0) {
2550           trackletWord[u] = bitWord[j];
2551           adcChannel[u]   = mADC[i];   // mapping onto the original adc-array to be in line with the digits-adc-ordering (21 channels in total on 1 mcm, 18 belonging to pads); mADC[i] should be >-1 in case bitWord[i]>0
2552
2553 // Finding label of MC track
2554           TH1F *hTrkRef = new TH1F("trackref", "trackref", 100000, 0, 100000);
2555           Int_t track[3];
2556           Int_t padcol = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u]);
2557           Int_t padcol_ngb = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u] - 1);
2558           Int_t padrow = 4 * (fRobPos / 2) + fMcmPos / 4;
2559           Int_t det = 30 * fSector + 6 * fStack + fLayer;
2560           for(Int_t iTimebin = feeParam->GetLinearFitStart(); iTimebin < feeParam->GetLinearFitEnd(); iTimebin++) {
2561               track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2562               track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2563               track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2564               hTrkRef->Fill(track[0]);
2565               if (track[1] != track[0] && track[1] != -1)
2566                   hTrkRef->Fill(track[1]);
2567               if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2568                   hTrkRef->Fill(track[2]);
2569               if (padcol_ngb >= 0) {
2570                   track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2571                   track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2572                   track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2573                   hTrkRef->Fill(track[0]);
2574                   if (track[1] != track[0] && track[1] != -1)
2575                       hTrkRef->Fill(track[1]);
2576                   if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2577                       hTrkRef->Fill(track[2]);
2578               }
2579           }
2580           trackRef[u] = hTrkRef->GetMaximumBin() - 1;
2581           delete hTrkRef;
2582           u = u + 1;
2583       }
2584   }
2585
2586   AliDataLoader *dl = AliRunLoader::GetRunLoader()->GetLoader("TRDLoader")->GetDataLoader("tracklets");
2587   if (!dl) {
2588     AliError("Could not get the tracklets data loader!");
2589   }
2590   else {
2591     TTree *trackletTree = dl->Tree();
2592     if (!trackletTree)
2593       dl->MakeTree();
2594     trackletTree = dl->Tree();
2595
2596    AliTRDtrackletMCM *trkl = new AliTRDtrackletMCM(); 
2597    TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
2598    if (!trkbranch)
2599        trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
2600     trkbranch->SetAddress(&trkl);
2601
2602     for (Int_t iTracklet = 0; iTracklet < fMaxTracklets; iTracklet++) {
2603         if (trackletWord[iTracklet] == 0)
2604             continue;
2605         trkl->SetTrackletWord(trackletWord[iTracklet]);
2606         trkl->SetDetector(30*fSector + 6*fStack + fLayer);
2607         trkl->SetROB(fRobPos);
2608         trkl->SetMCM(fMcmPos);
2609         trkl->SetLabel(trackRef[iTracklet]);
2610         trackletTree->Fill();
2611     }
2612     delete trkl;
2613     dl->WriteData("OVERWRITE");
2614   }
2615
2616   delete [] trackletWord;
2617   delete [] adcChannel; 
2618   delete [] trackRef;
2619   delete digman;
2620
2621   // to be done:
2622   // error measure for quality of fit (not necessarily needed for the trigger)
2623   // cluster quality threshold (not yet set)
2624   // electron probability
2625 }
2626 //_____________________________________________________________________________________
2627 void AliTRDmcmSim::GeneratefZSM1Dim()
2628 {
2629   //
2630   // Generate the array fZSM1Dim necessary
2631   // for the method ProduceRawStream
2632   //
2633
2634   // Fill the mapping
2635   // Supressed zeros indicated by -1 in digits array
2636   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) 
2637     {
2638       for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
2639         {
2640           
2641           if(fADCF[iadc][it]==-1)  // If is a supressed value
2642             {
2643               fZSM[iadc][it]=1;
2644             }
2645           else                    // Not suppressed
2646             {
2647               fZSM[iadc][it]=0;
2648             }
2649         }
2650     }
2651
2652   // Make the 1 dim projection
2653   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) 
2654     {
2655       for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
2656         {
2657           fZSM1Dim[iadc] &= fZSM[iadc][it];
2658         }
2659     }
2660 }
2661 //_______________________________________________________________________________________
2662 void AliTRDmcmSim::CopyArrays()
2663 {
2664   //
2665   // Initialize filtered data array with raw data
2666   // Method added for internal consistency
2667   //
2668
2669   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) 
2670     {
2671       for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
2672         {
2673           fADCF[iadc][it] = fADCR[iadc][it]; 
2674         }
2675     }
2676 }
2677 //_______________________________________________________________________________________
2678 void AliTRDmcmSim::StartfastZS(Int_t pads, Int_t timebins)
2679 {
2680   //
2681   // Initialize just the necessary elements to perform
2682   // the zero suppression in the digitizer
2683   //
2684    
2685   fFeeParam  = AliTRDfeeParam::Instance();
2686   fSimParam  = AliTRDSimParam::Instance();
2687   fNADC      = pads;      
2688   fNTimeBin  = timebins; 
2689
2690   if( fADCR == NULL ) 
2691     {
2692       fADCR    = new Int_t *[fNADC];
2693       fADCF    = new Int_t *[fNADC];
2694       fADCT    = new Int_t *[fNADC]; 
2695       fZSM     = new Int_t *[fNADC];
2696       fZSM1Dim = new Int_t  [fNADC];
2697     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ )
2698       {
2699         fADCR[iadc] = new Int_t[fNTimeBin];
2700         fADCF[iadc] = new Int_t[fNTimeBin];
2701         fADCT[iadc] = new Int_t[fNTimeBin]; 
2702         fZSM [iadc] = new Int_t[fNTimeBin];
2703       }
2704     }
2705
2706   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) 
2707     {
2708       for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
2709         {
2710           fADCR[iadc][it] =  0;
2711           fADCF[iadc][it] =  0;
2712           fADCT[iadc][it] = -1;  
2713           fZSM [iadc][it] =  1;   
2714         }
2715       fZSM1Dim[iadc] = 1;      
2716     }
2717   
2718   fInitialized = kTRUE;
2719 }
2720 //_______________________________________________________________________________________
2721 void AliTRDmcmSim::FlagDigitsArray(AliTRDarrayADC *tempdigs, Int_t valrow)
2722 {
2723   //
2724   // Modify the digits array to flag suppressed values
2725   //
2726
2727   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) 
2728     {
2729       for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
2730         {
2731           if(fZSM[iadc][it]==1)
2732             {
2733               tempdigs->SetData(valrow,iadc,it,-1);
2734             }
2735         }
2736     }
2737 }
2738 //_______________________________________________________________________________________
2739 void AliTRDmcmSim::RestoreZeros()
2740 {
2741   //
2742   // Restore the zero-suppressed values (set as -1) to the value 0
2743   //
2744
2745   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) 
2746     {
2747       for( Int_t it = 0 ; it < fNTimeBin ; it++ ) 
2748         {
2749           
2750           if(fADCF[iadc][it]==-1)  //if is a supressed zero, reset to zero
2751             {
2752               fADCF[iadc][it]=0;
2753               fADCR[iadc][it]=0;
2754             }     
2755         }
2756     }
2757
2758 }
2759