b096c82b3aa6909b491fbb17e56c311d3a8a0ccb
[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
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 void AliTRDmcmSim::ZSMapping()
846 {
847   //
848   // Zero Suppression Mapping implemented in TRAP chip
849   //
850   // See detail TRAP manual "Data Indication" section:
851   // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
852   //
853
854   Int_t eBIS = fFeeParam->GetEBsglIndThr();       // TRAP default = 0x4  (Tis=4)
855   Int_t eBIT = fFeeParam->GetEBsumIndThr();       // TRAP default = 0x28 (Tit=40)
856   Int_t eBIL = fFeeParam->GetEBindLUT();          // TRAP default = 0xf0
857                                                   // (lookup table accept (I2,I1,I0)=(111)
858                                                   // or (110) or (101) or (100))
859   Int_t eBIN = fFeeParam->GetEBignoreNeighbour(); // TRAP default = 1 (no neighbor sensitivity)
860   Int_t ep   = AliTRDfeeParam::GetPFeffectPedestal();
861
862   if( !CheckInitialized() ) return;
863
864   for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
865     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
866
867       // Get ADC data currently in filter buffer
868       Int_t ap = fADCF[iadc-1][it] - ep; // previous
869       Int_t ac = fADCF[iadc  ][it] - ep; // current
870       Int_t an = fADCF[iadc+1][it] - ep; // next
871
872       // evaluate three conditions
873       Int_t i0 = ( ac >=  ap && ac >=  an ) ? 0 : 1; // peak center detection
874       Int_t i1 = ( ap + ac + an > eBIT )    ? 0 : 1; // cluster
875       Int_t i2 = ( ac > eBIS )              ? 0 : 1; // absolute large peak
876
877       Int_t i = i2 * 4 + i1 * 2 + i0;    // Bit position in lookup table
878       Int_t d = (eBIL >> i) & 1;         // Looking up  (here d=0 means true
879                                          // and d=1 means false according to TRAP manual)
880
881       fZSM[iadc][it] &= d;
882       if( eBIN == 0 ) {  // turn on neighboring ADCs
883         fZSM[iadc-1][it] &= d;
884         fZSM[iadc+1][it] &= d;
885       }
886
887     }
888   }
889
890   // do 1 dim projection
891   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
892     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
893       fZSM1Dim[iadc] &= fZSM[iadc][it];
894     }
895   }
896
897 }
898
899 //_____________________________________________________________________________
900 void AliTRDmcmSim::DumpData( char *f, char *target )
901 {
902   //
903   // Dump data stored (for debugging).
904   // target should contain one or multiple of the following characters
905   //   R   for raw data
906   //   F   for filtered data
907   //   Z   for zero suppression map
908   //   S   Raw dat astream
909   // other characters are simply ignored
910   //
911
912   UInt_t tempbuf[1024];
913
914   if( !CheckInitialized() ) return;
915
916   std::ofstream of( f, std::ios::out | std::ios::app );
917   of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
918              fChaId, fSector, fStack, fLayer, fRobPos, fMcmPos );
919
920   for( int t=0 ; target[t] != 0 ; t++ ) {
921     switch( target[t] ) {
922     case 'R' :
923     case 'r' :
924       of << Form("fADCR (raw ADC data)\n");
925       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
926         of << Form("  ADC %02d: ", iadc);
927         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
928           of << Form("% 4d",  fADCR[iadc][it]);
929         }
930         of << Form("\n");
931       }
932       break;
933     case 'F' :
934     case 'f' :
935       of << Form("fADCF (filtered ADC data)\n");
936       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
937         of << Form("  ADC %02d: ", iadc);
938         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
939           of << Form("% 4d",  fADCF[iadc][it]);
940         }
941         of << Form("\n");
942       }
943       break;
944     case 'Z' :
945     case 'z' :
946       of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
947       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
948         of << Form("  ADC %02d: ", iadc);
949         if( fZSM1Dim[iadc] == 0 ) { of << " R   " ; } else { of << " .   "; } // R:read .:suppressed
950         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
951           if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
952         }
953         of << Form("\n");
954       }
955       break;
956     case 'S' :
957     case 's' :
958       Int_t s = ProduceRawStream( tempbuf, 1024 ); 
959       of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
960       of << Form("  address  data\n");
961       for( int i = 0 ; i < s ; i++ ) {
962         of << Form("  %04x     %08x\n", i, tempbuf[i]);
963       }
964     }
965   }
966 }
967
968 //_____________________________________________________________________________
969 void AliTRDmcmSim::FilterSimDeConvExpA(Int_t *source, Double_t *target
970                                      , Int_t n, Int_t nexp) 
971 {
972   //
973   // Exponential filter "analog"
974   // source will not be changed
975   //
976
977   Int_t    i = 0;
978   Int_t    k = 0;
979   Double_t reminder[2];
980   Double_t correction;
981   Double_t result;
982   Double_t rates[2];
983   Double_t coefficients[2];
984
985   // Initialize (coefficient = alpha, rates = lambda)
986   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
987
988   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
989   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
990   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
991   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
992   
993   coefficients[0] = c1;
994   coefficients[1] = c2;
995
996   Double_t dt = 0.1;
997   rates[0] = TMath::Exp(-dt/(r1));
998   rates[1] = TMath::Exp(-dt/(r2));
999
1000   // Attention: computation order is important
1001   correction = 0.0;
1002   for (k = 0; k < nexp; k++) {
1003     reminder[k] = 0.0;
1004   }
1005     
1006   for (i = 0; i < n; i++) {
1007
1008     result    = ((Double_t)source[i] - correction);    // no rescaling
1009     target[i] = result;
1010     
1011     for (k = 0; k < nexp; k++) {
1012       reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
1013     }
1014       
1015     correction = 0.0;
1016     for (k = 0; k < nexp; k++) {
1017       correction += reminder[k];
1018     }
1019   }
1020 }
1021
1022 //_____________________________________________________________________________
1023 void AliTRDmcmSim::FilterSimDeConvExpD(Int_t *source, Int_t *target, Int_t n
1024                                      , Int_t nexp) 
1025 {
1026   //
1027   // Exponential filter "digital"
1028   // source will not be changed
1029   //
1030
1031   Int_t i        = 0;
1032   Int_t fAlphaL  = 0;
1033   Int_t fAlphaS  = 0;
1034   Int_t fTailPed = 0;
1035   Int_t iAlphaL  = 0;
1036   Int_t iAlphaS  = 0;
1037
1038   // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
1039   // initialize (coefficient = alpha, rates = lambda)
1040
1041   Double_t dt = 0.1;
1042   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
1043   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
1044   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
1045   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
1046
1047   Int_t fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
1048   Int_t fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
1049   Int_t iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; //  9 bit paramter + fixed bits
1050   Int_t iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; //  9 bit paramter + fixed bits
1051
1052   if (nexp == 1) {
1053     fAlphaL = (Int_t) (c1 * 2048.0);
1054     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
1055   }
1056   if (nexp == 2) {
1057     fAlphaL = (Int_t) (c1 * 2048.0);
1058     fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
1059     iAlphaL = fAlphaL & 0x03FF;                         // 10 bit paramter
1060     iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400;              // 10 bit paramter + fixed bits
1061   }
1062   
1063   Double_t iAl = iAlphaL  / 2048.0;            // alpha  L: correspondence to floating point numbers
1064   Double_t iAs = iAlphaS  / 2048.0;            // alpha  S: correspondence to floating point numbers
1065   Double_t iLl = iLambdaL / 2048.0;            // lambda L: correspondence to floating point numbers
1066   Double_t iLs = iLambdaS / 2048.0;            // lambda S: correspondence to floating point numbers
1067
1068   Int_t h1;
1069   Int_t h2;
1070   Int_t rem1;
1071   Int_t rem2;
1072   Int_t correction;
1073   Int_t result;
1074   Int_t iFactor = ((Int_t) fFeeParam->GetPFeffectPedestal() ) << 2;
1075
1076   Double_t xi = 1 - (iLl*iAs + iLs*iAl);             // Calculation of equilibrium values of the
1077   rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
1078   rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
1079   
1080   // further initialization
1081   if ((rem1 + rem2) > 0x0FFF) {
1082     correction = 0x0FFF;
1083   } 
1084   else {
1085     correction = (rem1 + rem2) & 0x0FFF;
1086   }
1087
1088   fTailPed = iFactor - correction;
1089
1090   for (i = 0; i < n; i++) {
1091
1092     result = (source[i]  - correction);
1093     if (result < 0) { // Too much undershoot
1094       result = 0;
1095     }
1096
1097     target[i] = result;
1098                                                         
1099     h1 = (rem1 + ((iAlphaL * result) >> 11));
1100     if (h1 > 0x0FFF) {
1101       h1 = 0x0FFF;
1102     } 
1103     else {
1104       h1 &= 0x0FFF;
1105     }
1106
1107     h2 = (rem2 + ((iAlphaS * result) >> 11));
1108     if (h2 > 0x0FFF) {
1109       h2 = 0x0FFF;
1110     } 
1111     else {
1112       h2 &= 0x0FFF;
1113     }
1114   
1115     rem1 = (iLambdaL * h1 ) >> 11;
1116     rem2 = (iLambdaS * h2 ) >> 11;
1117     
1118     if ((rem1 + rem2) > 0x0FFF) {
1119       correction = 0x0FFF;
1120     } 
1121     else {
1122       correction = (rem1 + rem2) & 0x0FFF;
1123     }
1124
1125   }
1126
1127 }
1128
1129 //_____________________________________________________________________________
1130 void AliTRDmcmSim::FilterSimDeConvExpMI(Int_t *source, Double_t *target
1131                                       , Int_t n) 
1132 {
1133   //
1134   // Exponential filter (M. Ivanov)
1135   // source will not be changed
1136   //
1137
1138   Int_t i = 0;
1139   Double_t sig1[100];
1140   Double_t sig2[100];
1141   Double_t sig3[100];
1142
1143   for (i = 0; i < n; i++) {
1144     sig1[i] = (Double_t)source[i];
1145   }
1146
1147   Float_t dt      = 0.1;
1148   Float_t lambda0 = (1.0 / fFeeParam->GetTFr2()) * dt;
1149   Float_t lambda1 = (1.0 / fFeeParam->GetTFr1()) * dt;
1150
1151   FilterSimTailMakerSpline( sig1, sig2, lambda0, n);
1152   FilterSimTailCancelationMI( sig2, sig3, 0.7, lambda1, n);
1153
1154   for (i = 0; i < n; i++) {
1155     target[i] = sig3[i];
1156   }
1157
1158 }
1159
1160 //______________________________________________________________________________
1161 void AliTRDmcmSim::FilterSimTailMakerSpline(Double_t *ampin, Double_t *ampout
1162                                           , Double_t lambda, Int_t n) 
1163 {
1164   //
1165   // Special filter (M. Ivanov)
1166   //
1167
1168   Int_t    i = 0;
1169   Double_t l = TMath::Exp(-lambda*0.5);
1170   Double_t in[1000];
1171   Double_t out[1000];
1172
1173   // Initialize in[] and out[] goes 0 ... 2*n+19
1174   for (i = 0; i < n*2+20; i++) {
1175     in[i] = out[i] = 0;
1176   }
1177
1178   // in[] goes 0, 1
1179   in[0] = ampin[0];
1180   in[1] = (ampin[0] + ampin[1]) * 0.5;
1181    
1182   // Add charge to the end
1183   for (i = 0; i < 22; i++) {
1184     in[2*(n-1)+i] = ampin[n-1]; // in[] goes 2*n-2, 2*n-1, ... , 2*n+19 
1185   }
1186
1187   // Use arithmetic mean
1188   for (i = 1; i < n-1; i++) {
1189     in[2*i]   = ampin[i];    // in[] goes 2, 3, ... , 2*n-4, 2*n-3
1190     in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
1191   }
1192
1193   Double_t temp;
1194   out[2*n]    = in[2*n];
1195   temp        = 0;
1196   for (i = 2*n; i >= 0; i--) {
1197     out[i]    = in[i] + temp;
1198     temp      = l*(temp+in[i]);
1199   }
1200
1201   for (i = 0; i < n; i++){
1202     //ampout[i] = out[2*i+1];  // org
1203     ampout[i] = out[2*i];
1204   }
1205
1206 }
1207
1208 //______________________________________________________________________________
1209 void AliTRDmcmSim::FilterSimTailCancelationMI(Double_t *ampin, Double_t *ampout
1210                                             , Double_t norm, Double_t lambda
1211                                             , Int_t n) 
1212 {
1213   //
1214   // Special filter (M. Ivanov)
1215   //
1216
1217   Int_t    i = 0;
1218
1219   Double_t l = TMath::Exp(-lambda*0.5);
1220   Double_t k = l*(1.0 - norm*lambda*0.5);
1221   Double_t in[1000];
1222   Double_t out[1000];
1223
1224   // Initialize in[] and out[] goes 0 ... 2*n+19
1225   for (i = 0; i < n*2+20; i++) {
1226     in[i] = out[i] = 0;
1227   }
1228
1229   // in[] goes 0, 1
1230   in[0] = ampin[0];
1231   in[1] = (ampin[0]+ampin[1])*0.5;
1232
1233   // Add charge to the end
1234   for (i =-2; i < 22; i++) {
1235     // in[] goes 2*n-4, 2*n-3, ... , 2*n+19 
1236     in[2*(n-1)+i] = ampin[n-1];
1237   }
1238
1239   for (i = 1; i < n-2; i++) {
1240     // in[] goes 2, 3, ... , 2*n-6, 2*n-5
1241     in[2*i]    = ampin[i];
1242     in[2*i+1]  = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
1243     //in[2*i+1]  = ((ampin[i]+ampin[i+1]))/2.0;
1244   }
1245
1246   Double_t temp;
1247   out[0] = in[0];
1248   temp   = in[0];
1249   for (i = 1; i <= 2*n; i++) {
1250     out[i] = in[i] + (k-l)*temp;
1251     temp   = in[i] +  k   *temp;
1252   }
1253
1254   for (i = 0; i < n; i++) {
1255     //ampout[i] = out[2*i+1];  // org
1256     //ampout[i] = TMath::Max(out[2*i+1],0.0);  // org
1257     ampout[i] = TMath::Max(out[2*i],0.0);
1258   }
1259 }
1260
1261
1262 //_____________________________________________________________________________________
1263 //the following filter uses AliTRDtrapAlu-class
1264
1265 void AliTRDmcmSim::FilterSimDeConvExpEl(Int_t *source, Int_t *target, Int_t n, Int_t nexp) {
1266   //static Int_t count = 0;
1267  
1268   Double_t dt = 0.1;
1269   Double_t r1 = (Double_t)fFeeParam->GetTFr1();
1270   Double_t r2 = (Double_t)fFeeParam->GetTFr2();
1271   Double_t c1 = (Double_t)fFeeParam->GetTFc1();
1272   Double_t c2 = (Double_t)fFeeParam->GetTFc2();
1273   
1274   nexp = 1;
1275
1276   //it is assumed that r1,r2,c1,c2 are given such, that the configuration values are in the ranges according to TRAP-manual
1277   //parameters need to be adjusted
1278   AliTRDtrapAlu lambdaL;
1279   AliTRDtrapAlu lambdaS;
1280   AliTRDtrapAlu alphaL;
1281   AliTRDtrapAlu alphaS;
1282   
1283   AliTRDtrapAlu correction;
1284   AliTRDtrapAlu result;
1285   AliTRDtrapAlu bufL;
1286   AliTRDtrapAlu bufS;
1287  
1288   AliTRDtrapAlu bSource;
1289   
1290   lambdaL.Init(1,11);
1291   lambdaS.Init(1,11);
1292   alphaL.Init(1,11);
1293   alphaS.Init(1,11);
1294   
1295   //count=count+1;
1296
1297   lambdaL.AssignDouble(TMath::Exp(-dt/r1));
1298   lambdaS.AssignDouble(TMath::Exp(-dt/r2));
1299   alphaL.AssignDouble(c1); // in AliTRDfeeParam the number of exponentials is set and also the according time constants
1300   alphaS.AssignDouble(c2); // later it should be: alphaS=1-alphaL
1301   
1302   //data is enlarged to 12 bits, including 2 bits after the comma; class AliTRDtrapAlu is used to handle arithmetics correctly
1303   correction.Init(10,2);
1304   result.Init(10,2);
1305   bufL.Init(10,2);
1306   bufS.Init(10,2);
1307   bSource.Init(10,2);
1308   
1309   for(Int_t i = 0; i < n; i++) {
1310     bSource.AssignInt(source[i]);
1311     result = bSource - correction; // subtraction can produce an underflow
1312     if(result.GetSign() == kTRUE) {
1313       result.AssignInt(0);
1314     }
1315     
1316     //target[i] = result.GetValuePre();  // later, target and source should become AliTRDtrapAlu,too in order to simulate the 10+2Bits through the filter properly
1317     
1318     target[i] = result.GetValue(); // 12 bit-value; to get the corresponding integer value, target must be shifted: target>>2 
1319
1320     //printf("target-Wert zur Zeit %d : %d",i,target[i]);
1321     //printf("\n");
1322     
1323     bufL  =  bufL + (result * alphaL);
1324     bufL  =  bufL * lambdaL; 
1325     
1326     bufS  =  bufS + (result * alphaS);
1327     bufS  =  bufS * lambdaS;  // eventually this should look like:
1328     // bufS = (bufS + (result - result * alphaL)) * lambdaS // alphaS=1-alphaL; then alphaS-variable is not needed any more
1329
1330     correction = bufL + bufS; //check for overflow intrinsic; if overflowed, correction is set to 0x03FF
1331   }
1332  
1333   
1334 }
1335
1336
1337
1338
1339
1340
1341
1342 //__________________________________________________________________________________
1343
1344
1345 // in order to use the Tracklets, please first 
1346 // -- set AliTRDfeeParam::fgkTracklet to kTRUE, in order to switch on Tracklet-calculation
1347 // -- set AliTRDfeeParam::fgkTFtype   to 3, in order to use the new tail cancellation filter
1348 //    currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3
1349 // -- set AliTRDfeeParam::fgkMCTrackletOutput to kTRUE, if you want to use the Tracklet output container with           information about the Tracklet position (MCM, channel number)
1350
1351 // 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
1352
1353 void AliTRDmcmSim::Tracklet(){
1354     // tracklet calculation
1355     // if you use this code after a simulation, please make sure the same filter-settings as in the simulation are set in AliTRDfeeParam
1356
1357   if(!CheckInitialized()){ return; }
1358   
1359   Bool_t filtered = kTRUE;
1360   
1361   
1362   
1363   AliTRDtrapAlu data;
1364   data.Init(10,2);
1365   if(fADCT[0][0]==-1){                      // check if filter was applied
1366     filtered = kFALSE;
1367     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1368       for( Int_t iT = 0 ; iT < fNTimeBin ; iT++ ) {
1369         data.AssignInt(fADCR[iadc][iT]);
1370         fADCT[iadc][iT] = data.GetValue(); // all incoming values are positive 10+2 bit values; if el.filter was called, this is done correctly
1371       }
1372     }
1373    
1374   }
1375   
1376   // the online ordering of mcm's is reverse to the TRAP-manual-ordering! reverse fADCT (to be consistent to TRAP), then do all calculations
1377   // reverse fADCT:
1378   Int_t** rev0 = new Int_t *[fNADC];
1379   Int_t** rev1 = new Int_t *[fNADC];
1380   
1381   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1382     rev0[iadc] = new Int_t[fNTimeBin];
1383     rev1[iadc] = new Int_t[fNTimeBin];
1384     for( Int_t iT = 0; iT < fNTimeBin; iT++) {
1385       if( iadc <= fNADC-iadc-1 ) {
1386         rev0[iadc][iT]  = fADCT[fNADC-iadc-1][iT];
1387         rev1[iadc][iT]  = fADCT[iadc][iT];
1388         fADCT[iadc][iT] = rev0[iadc][iT];
1389       }
1390       else {
1391         rev0[iadc][iT]  = rev1[fNADC-iadc-1][iT];
1392         fADCT[iadc][iT] = rev0[iadc][iT];
1393       }
1394     }
1395   }
1396   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1397     delete[] rev0[iadc];
1398     delete[] rev1[iadc];
1399   }
1400   
1401   delete[] rev0;
1402   delete[] rev1;
1403   
1404   rev0 = NULL;
1405   rev1 = NULL;
1406     
1407   // get the filtered pedestal; supports only electronic tail-cancellation filter
1408   AliTRDtrapAlu filPed;
1409   Int_t ep = 0;
1410   Int_t *ieffped = new Int_t[fNTimeBin];
1411   for(Int_t iT = 0; iT < fNTimeBin; iT++){
1412     ieffped[iT] = ep; 
1413   }
1414  
1415   if( filtered == kTRUE ) {
1416     if( fFeeParam->IsPFon() ){
1417       ep = fFeeParam->GetPFeffectPedestal();
1418     }
1419     Int_t      nexp  = fFeeParam->GetTFnExp();
1420     Int_t  *isource  = new Int_t[fNTimeBin];
1421     filPed.Init(10,2);
1422     filPed.AssignInt(ep);           
1423     Int_t epf = filPed.GetValue();  
1424     for(Int_t iT = 0; iT < fNTimeBin; iT++){
1425       isource[iT] = ep;                  
1426       ieffped[iT] = epf;
1427     }
1428  
1429     if( fFeeParam->IsTFon() ) {
1430       FilterSimDeConvExpEl( isource, ieffped, fNTimeBin, nexp);
1431     }
1432   
1433     delete[] isource;
1434   }
1435   
1436   //the following values should go to AliTRDfeeParam once they are defined; then they have to be read in properly
1437   //naming follows conventions in TRAP-manual
1438   
1439   
1440   Bool_t bVBY = kTRUE;                         // cluster-verification bypass
1441
1442   Double_t cQTParam = 0;                      // cluster quality threshold; granularity 2^-10; range: 0<=cQT/2^-10<=2^-4 - 2^-10
1443   AliTRDtrapAlu cQTAlu; 
1444   cQTAlu.Init(1,10,0,63);
1445   cQTAlu.AssignDouble(cQTParam);
1446   Int_t cQT = cQTAlu.GetValue();
1447
1448   // linear fit 
1449   Int_t tFS = fFeeParam->GetLinearFitStart();  // linear fit start
1450   Int_t tFE = fFeeParam->GetLinearFitEnd();    // linear fit stop
1451    
1452   // charge accumulators
1453   Int_t tQS0 = fFeeParam->GetQacc0Start();     // start-time for charge-accumulator 0
1454   Int_t tQE0 = fFeeParam->GetQacc0End();       // stop-time for charge-accumulator 0
1455   Int_t tQS1 = fFeeParam->GetQacc1Start();     // start-time for charge-accumulator 1 
1456   Int_t tQE1 = fFeeParam->GetQacc1End();       // stop-time for charge-accumulator 1
1457   // values set such that tQS0=tFS; tQE0=tQS1-1; tFE=tQE1; want to do (QS0+QS1)/N
1458  
1459   Double_t cTHParam = (Double_t)fFeeParam->GetMinClusterCharge(); // cluster charge threshold
1460   AliTRDtrapAlu cTHAlu;  
1461   cTHAlu.Init(12,2);
1462   cTHAlu.AssignDouble(cTHParam);
1463   Int_t cTH = cTHAlu.GetValue();                                 // cTH used for comparison
1464
1465   struct List_t {
1466     List_t *next;
1467     Int_t iadc;
1468     Int_t value;
1469   };
1470   
1471   List_t selection[7];            // list with 7 elements
1472   List_t *list = NULL;
1473   List_t *listLeft = NULL;
1474     
1475   Int_t* qsum = new Int_t[fNADC];
1476    
1477   // fit sums
1478   AliTRDtrapAlu qsumAlu;
1479   qsumAlu.Init(12,2);           // charge sum will be 12+2 bits
1480   AliTRDtrapAlu dCOGAlu; 
1481   dCOGAlu.Init(1,7,0,127);      // COG will be 1+7 Bits; maximum 1 - 2^-7 for LUT
1482   AliTRDtrapAlu yrawAlu;
1483   yrawAlu.Init(1,8,-1,255);
1484   AliTRDtrapAlu yAlu;
1485   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 -
1486   AliTRDtrapAlu xAlu;
1487   xAlu.Init(5,8);               // 8 past-comma bits because value will be added/multiplied to another value with this accuracy
1488   AliTRDtrapAlu xxAlu;
1489   xxAlu.Init(10,0);            
1490   AliTRDtrapAlu yyAlu;
1491   yyAlu.Init(1,16,0,0xFFFF);    // maximum is 2^16-1; 16Bit for past-commas
1492   AliTRDtrapAlu xyAlu;
1493   xyAlu.Init(6,8);
1494   AliTRDtrapAlu XAlu;
1495   XAlu.Init(9,0);
1496   AliTRDtrapAlu XXAlu;
1497   XXAlu.Init(14,0);
1498   AliTRDtrapAlu YAlu;
1499   YAlu.Init(5,8);               // 14 bit, 1 is sign-bit; therefore only 13 bit 
1500   AliTRDtrapAlu YYAlu;
1501   YYAlu.Init(5,16);
1502   AliTRDtrapAlu XYAlu;
1503   XYAlu.Init(8,8);              // 17 bit, 1 is sign-bit; therefore only 16 bit        
1504   AliTRDtrapAlu qtruncAlu;
1505   qtruncAlu.Init(12,0);
1506   AliTRDtrapAlu QT0Alu;
1507   QT0Alu.Init(15,0);
1508   AliTRDtrapAlu QT1Alu;
1509   QT1Alu.Init(16,0);
1510
1511   AliTRDtrapAlu oneAlu;
1512   oneAlu.Init(1,8);
1513  
1514   
1515   AliTRDtrapAlu inverseNAlu;
1516   inverseNAlu.Init(1,8);        // simulates the LUT for 1/N
1517   AliTRDtrapAlu MeanChargeAlu;  // mean charge in ADC counts
1518   MeanChargeAlu.Init(8,0);
1519   AliTRDtrapAlu TotalChargeAlu;
1520   TotalChargeAlu.Init(17,8);
1521   //nr of post comma bits should be the same for inverseN and TotalCharge
1522   
1523   
1524   SetPosLUT();                    // initialize the position correction LUT for this MCM;
1525
1526
1527   // fit-sums; remapping!; 0,1,2->0; 1,2,3->1; ... 18,19,20->18
1528   Int_t *X   = new Int_t[fNADC-2];
1529   Int_t *XX  = new Int_t[fNADC-2];
1530   Int_t *Y   = new Int_t[fNADC-2];
1531   Int_t *YY  = new Int_t[fNADC-2];
1532   Int_t *XY  = new Int_t[fNADC-2];
1533   Int_t *N   = new Int_t[fNADC-2];
1534   Int_t *QT0 = new Int_t[fNADC-2]; // accumulated charge
1535   Int_t *QT1 = new Int_t[fNADC-2]; // accumulated charge
1536   
1537   for (Int_t iCol = 0; iCol < fNADC-2; iCol++) { 
1538       
1539       // initialize fit-sums 
1540       X[iCol]   = 0;
1541       XX[iCol]  = 0;
1542       Y[iCol]   = 0;
1543       YY[iCol]  = 0;
1544       XY[iCol]  = 0;
1545       N[iCol]   = 0;
1546       QT0[iCol] = 0;
1547       QT1[iCol] = 0;
1548   }
1549   
1550
1551   filPed.Init(7,2);                         // convert filtered pedestal into 7+2Bits
1552   
1553   for(Int_t iT = 0; iT < fNTimeBin; iT++){
1554     
1555     if(iT<tFS || iT>=tFE) continue;         // linear fit yes/no? 
1556
1557     // reset
1558     Int_t portChannel[4]   = {-1,-1,-1,-1};   
1559     Int_t clusterCharge[4] = {0,0,0,0};
1560     Int_t leftCharge[4]    = {0,0,0,0};
1561     Int_t centerCharge[4]  = {0,0,0,0}; 
1562     Int_t rightCharge[4]   = {0,0,0,0};
1563     
1564     Int_t mark = 0;
1565     
1566     filPed.AssignFormatted(ieffped[iT]);   // no size-checking when using AssignFormatted; ieffped>=0
1567     filPed = filPed;                       // this checks the size
1568     
1569     ieffped[iT] = filPed.GetValue();
1570         
1571     for(Int_t i = 0; i<7; i++){
1572       selection[i].next       = NULL;
1573       selection[i].iadc       =   -1;     // value of -1: invalid adc
1574       selection[i].value      =    0;
1575    
1576     }
1577     // selection[0] is starting list-element; just for pointing
1578
1579     // loop over inner adc's 
1580     for (Int_t iCol = 1; iCol < fNADC-1; iCol++) { 
1581       
1582       Int_t left   = fADCT[iCol-1][iT]; 
1583       Int_t center = fADCT[iCol][iT];
1584       Int_t right  = fADCT[iCol+1][iT];  
1585
1586       Int_t sum = left + center + right;            // cluster charge sum
1587       qsumAlu.AssignFormatted(sum);    
1588       qsumAlu = qsumAlu;                        // size-checking; redundant
1589  
1590       qsum[iCol] = qsumAlu.GetValue(); 
1591       
1592       //hit detection and masking
1593       if(center>=left){
1594         if(center>right){
1595           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
1596             mark |= 1;                              // marker
1597           }
1598         }
1599       }
1600       mark = mark<<1;                
1601     }
1602     mark = mark>>1;
1603
1604        
1605     // get selection of 6 adc's and sort,starting with greatest values
1606
1607     //read three from right side and sort (primitive sorting algorithm)
1608     Int_t i = 0; // adc number
1609     Int_t j = 1; // selection number
1610     while(i<fNADC-2 && j<=3){
1611       i = i + 1;
1612       if( ((mark>>(i-1)) & 1) == 1) {
1613         selection[j].iadc  = fNADC-1-i;
1614         selection[j].value = qsum[fNADC-1-i]>>6;   // for hit-selection only the first 8 out of the 14 Bits are used for comparison
1615         
1616         // insert into sorted list
1617         listLeft = &selection[0];
1618         list = listLeft->next;
1619         
1620         if(list!=NULL) {
1621           while((list->next != NULL) && (selection[j].value <= list->value)){
1622             listLeft = list;
1623             list = list->next;
1624           }
1625           
1626           if(selection[j].value<=list->value){
1627             selection[j].next = list->next;
1628             list->next = &selection[j];
1629           }
1630           else {
1631             listLeft->next = &selection[j];
1632             selection[j].next = list;
1633           }
1634         }
1635         else{
1636           listLeft->next = &selection[j];
1637           selection[j].next = list;
1638         }
1639         
1640         j = j + 1;
1641       }
1642     }
1643
1644
1645     // read three from left side
1646     Int_t k = fNADC-2;
1647     while(k>i && j<=6) {
1648       if( ((mark>>(k-1)) & 1) == 1) {
1649         selection[j].iadc  = fNADC-1-k;
1650         selection[j].value = qsum[fNADC-1-k]>>6;
1651         
1652         listLeft = &selection[0];
1653         list = listLeft->next;
1654         
1655         if(list!=NULL){
1656           while((list->next != NULL) && (selection[j].value <= list->value)){
1657             listLeft = list;
1658             list = list->next;
1659           }
1660         
1661           if(selection[j].value<=list->value){
1662             selection[j].next = list->next;
1663             list->next = &selection[j];
1664           }
1665           else {
1666             listLeft->next = &selection[j];
1667             selection[j].next = list;
1668           }
1669         }
1670         else{
1671           listLeft->next = &selection[j];
1672           selection[j].next = list;
1673         }
1674
1675         j = j + 1;
1676       }
1677       k = k - 1;
1678     }
1679
1680     // get the four with greatest charge-sum
1681     list = &selection[0];
1682     for(i = 0; i<4; i++){
1683       if(list->next == NULL) continue;
1684       list = list->next;
1685       if(list->iadc == -1) continue;
1686       Int_t adc = list->iadc;                              // channel number with selected hit
1687       
1688       // the following arrays contain the four chosen channels in 1 time-bin
1689       portChannel[i]   = adc; 
1690       clusterCharge[i] = qsum[adc];
1691       leftCharge[i]    = fADCT[adc-1][iT] - ieffped[iT]; // reduce by filtered pedestal (pedestal is part of the signal)
1692       centerCharge[i]  = fADCT[adc][iT] - ieffped[iT];           
1693       rightCharge[i]   = fADCT[adc+1][iT] - ieffped[iT];         
1694     }
1695
1696     // arithmetic unit
1697     
1698     // cluster verification
1699     if(!bVBY){
1700       for(i = 0; i<4; i++){
1701         Int_t lr = leftCharge[i]*rightCharge[i]*1024;
1702         Int_t cc = centerCharge[i]*centerCharge[i]*cQT;
1703         if (lr>=cc){
1704           portChannel[i]   = -1;                                 // set to invalid address 
1705           clusterCharge[i] = 0;
1706         }
1707       }
1708     }
1709
1710     // fit-sums of valid channels
1711     // local hit position
1712     for(i = 0; i<4; i++){
1713       if (centerCharge[i] ==  0) {
1714         portChannel[i] = -1; 
1715       }// prevent division by 0
1716       
1717       if (portChannel[i]  == -1) continue;
1718       
1719       Double_t dCOG = (Double_t)(rightCharge[i]-leftCharge[i])/centerCharge[i];
1720        
1721       Bool_t sign = (dCOG>=0.0) ? kFALSE : kTRUE;
1722       dCOG = (sign == kFALSE) ? dCOG : -dCOG;     // AssignDouble doesn't allow for signed doubles
1723       dCOGAlu.AssignDouble(dCOG);
1724       Int_t iLUTpos = dCOGAlu.GetValue();       // steers position in LUT
1725             
1726       dCOG = dCOG/2;
1727       yrawAlu.AssignDouble(dCOG);
1728       Int_t iCOG = yrawAlu.GetValue();
1729       Int_t y = iCOG + fPosLUT[iLUTpos % 128];    // local position in pad-units
1730       yrawAlu.AssignFormatted(y);               // 0<y<1           
1731       yAlu  = yrawAlu;                        // convert to 16 past-comma bits
1732       
1733       if(sign == kTRUE) yAlu.SetSign(-1);       // buffer width of 9 bits; sign on real (not estimated) position
1734       xAlu.AssignInt(iT);                       // buffer width of 5 bits 
1735       
1736
1737       xxAlu = xAlu * xAlu;                  // buffer width of 10 bits -> fulfilled by x*x       
1738       
1739       yyAlu = yAlu * yAlu;                  // buffer width of 16 bits
1740    
1741       xyAlu = xAlu * yAlu;                  // buffer width of 14 bits
1742                   
1743       Int_t adc = portChannel[i]-1;              // remapping! port-channel contains channel-nr. of inner adc's (1..19; mapped to 0..18)
1744
1745       // calculate fit-sums recursively
1746       // interpretation of their bit-length is given as comment
1747       
1748       // be aware that the accuracy of the result of a calculation is always determined by the accuracy of the less accurate value
1749
1750       XAlu.AssignFormatted(X[adc]);
1751       XAlu = XAlu + xAlu;                   // buffer width of 9 bits 
1752       X[adc] = XAlu.GetValue();
1753              
1754       XXAlu.AssignFormatted(XX[adc]);
1755       XXAlu = XXAlu + xxAlu;                // buffer width of 14 bits    
1756       XX[adc] = XXAlu.GetValue();
1757
1758       if (Y[adc] < 0) {
1759         YAlu.AssignFormatted(-Y[adc]);          // make sure that only positive values are assigned; sign-setting must be done by hand
1760         YAlu.SetSign(-1);
1761       }
1762       else {
1763         YAlu.AssignFormatted(Y[adc]);
1764         YAlu.SetSign(1);
1765       }
1766         
1767       YAlu = YAlu + yAlu;                   // buffer width of 14 bits (8 past-comma);     
1768       Y[adc] = YAlu.GetSignedValue();
1769             
1770       YYAlu.AssignFormatted(YY[adc]);
1771       YYAlu = YYAlu + yyAlu;                // buffer width of 21 bits (16 past-comma) 
1772       YY[adc] = YYAlu.GetValue();
1773            
1774       if (XY[adc] < 0) {
1775         XYAlu.AssignFormatted(-XY[adc]);
1776         XYAlu.SetSign(-1);
1777       }
1778       else {
1779         XYAlu.AssignFormatted(XY[adc]);
1780         XYAlu.SetSign(1);
1781       }
1782
1783       XYAlu = XYAlu + xyAlu;                // buffer allows 17 bits (8 past-comma) 
1784       XY[adc] = XYAlu.GetSignedValue();
1785             
1786       N[adc]  = N[adc] + 1;
1787    
1788
1789       // accumulated charge
1790       qsumAlu.AssignFormatted(qsum[adc+1]); // qsum was not remapped!
1791       qtruncAlu = qsumAlu;
1792
1793       if(iT>=tQS0 && iT<=tQE0){
1794         QT0Alu.AssignFormatted(QT0[adc]);
1795         QT0Alu = QT0Alu + qtruncAlu;
1796         QT0[adc] = QT0Alu.GetValue();
1797         //interpretation of QT0 as 12bit-value (all pre-comma); is this as it should be done?; buffer allows 15 Bit
1798       }
1799       
1800       if(iT>=tQS1 && iT<=tQE1){
1801         QT1Alu.AssignFormatted(QT1[adc]);
1802         QT1Alu = QT1Alu + qtruncAlu;
1803         QT1[adc] = QT1Alu.GetValue();
1804         //interpretation of QT1 as 12bit-value; buffer allows 16 Bit
1805       }
1806     }// i
1807       
1808     // remapping is done!!
1809      
1810   }//iT
1811  
1812   
1813     
1814   // tracklet-assembly
1815   
1816   // put into AliTRDfeeParam and take care that values are in proper range
1817   const Int_t cTCL = 1;      // left adc: number of hits; 8<=TCL<=31 (?? 1<=cTCL<+8 ??) 
1818   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%)
1819   
1820   Int_t mPair   = 0;         // marker for possible tracklet pairs
1821   Int_t* hitSum = new Int_t[fNADC-3];
1822   // hitSum[0] means: hit sum of remapped channels 0 and 1; hitSum[17]: 17 and 18; 
1823   
1824   // check for all possible tracklet-pairs of adjacent channels (two are merged); mark the left channel of the chosen pairs
1825   for (Int_t iCol = 0; iCol < fNADC-3; iCol++) {
1826     hitSum[iCol] = N[iCol] + N[iCol+1];
1827     if ((N[iCol]>=cTCL) && (hitSum[iCol]>=cTCT)) {
1828         mPair |= 1;         // mark as possible channel-pair
1829      
1830     }
1831     mPair = mPair<<1;
1832   }
1833   mPair = mPair>>1;
1834   
1835   List_t* selectPair = new List_t[fNADC-2];      // list with 18 elements (0..18) containing the left channel-nr and hit sums
1836                                                  // selectPair[18] is starting list-element just for pointing
1837   for(Int_t k = 0; k<fNADC-2; k++){
1838       selectPair[k].next       = NULL;
1839       selectPair[k].iadc       =   -1;           // invalid adc
1840       selectPair[k].value      =    0;
1841    
1842     }
1843
1844  list = NULL;
1845  listLeft = NULL;
1846   
1847   // read marker and sort according to hit-sum
1848   
1849   Int_t adcL  = 0;            // left adc-channel-number (remapped)
1850   Int_t selNr = 0;            // current number in list
1851   
1852   // insert marked channels into list and sort according to hit-sum
1853   while(adcL < fNADC-3 && selNr < fNADC-3){
1854      
1855     if( ((mPair>>((fNADC-4)-(adcL))) & 1) == 1) {
1856       selectPair[selNr].iadc  = adcL;
1857       selectPair[selNr].value = hitSum[adcL];   
1858       
1859       listLeft = &selectPair[fNADC-3];
1860       list = listLeft->next;
1861         
1862       if(list!=NULL) {
1863         while((list->next != NULL) && (selectPair[selNr].value <= list->value)){
1864           listLeft = list;
1865           list = list->next;
1866         }
1867         
1868         if(selectPair[selNr].value <= list->value){
1869           selectPair[selNr].next = list->next;
1870           list->next = &selectPair[selNr];
1871         }
1872         else {
1873           listLeft->next = &selectPair[selNr];
1874           selectPair[selNr].next = list;
1875         }
1876         
1877       }
1878       else{
1879         listLeft->next = &selectPair[selNr];
1880         selectPair[selNr].next = list;
1881       }
1882       
1883       selNr = selNr + 1;
1884     }
1885     adcL = adcL + 1;
1886   }
1887   
1888   //select up to 4 channels with maximum number of hits
1889   Int_t lpairChannel[4] = {-1,-1,-1,-1}; // save the left channel-numbers of pairs with most hit-sum
1890   Int_t rpairChannel[4] = {-1,-1,-1,-1}; // save the right channel, too; needed for detecting double tracklets
1891   list = &selectPair[fNADC-3];
1892   
1893   for (Int_t i = 0; i<4; i++) {
1894     if(list->next == NULL) continue;
1895     list = list->next;
1896     if(list->iadc == -1) continue;
1897     lpairChannel[i] = list->iadc;        // channel number with selected hit
1898     rpairChannel[i] = lpairChannel[i]+1;
1899   }
1900   
1901   // avoid submission of double tracklets  
1902   for (Int_t i = 3; i>0; i--) {
1903     for (Int_t j = i-1; j>-1; j--) {
1904       if(lpairChannel[i] == rpairChannel[j]) {
1905         lpairChannel[i] = -1;
1906         rpairChannel[i] = -1;
1907         break;
1908       }
1909       /* if(rpairChannel[i] == lpairChannel[j]) {
1910         lpairChannel[i] = -1;
1911         rpairChannel[i] = -1;
1912         break;
1913         }*/
1914     }
1915   }
1916   
1917   // merging of the fit-sums of the remainig channels
1918   // assume same data-word-width as for fit-sums for 1 channel
1919   // relative scales!
1920   Int_t mADC[4];                      
1921   Int_t mN[4];
1922   Int_t mQT0[4];
1923   Int_t mQT1[4];
1924   Int_t mX[4];
1925   Int_t mXX[4];
1926   Int_t mY[4];
1927   Int_t mYY[4];
1928   Int_t mXY[4];
1929   Int_t mOffset[4];
1930   Int_t mSlope[4];
1931   Int_t mMeanCharge[4]; 
1932   Int_t inverseN = 0;
1933   Double_t invN = 0;
1934   Int_t one = 0;
1935
1936   for (Int_t i = 0; i<4; i++){
1937     mADC[i] = -1;                        // set to invalid number
1938     mN[i]   =  0;
1939     mQT0[i] =  0;
1940     mQT1[i] =  0;
1941     mX[i]   =  0;
1942     mXX[i]  =  0;
1943     mY[i]   =  0;
1944     mYY[i]  =  0;
1945     mXY[i]  =  0;
1946     mOffset[i] = 0;
1947     mSlope[i]  = 0;
1948     mMeanCharge[i] = 0;
1949   }
1950   
1951   oneAlu.AssignInt(1);
1952   one = oneAlu.GetValue();              // one with 8 past comma bits
1953  
1954   for (Int_t i = 0; i<4; i++){
1955           
1956
1957     mADC[i] = lpairChannel[i];          // mapping of merged sums to left channel nr. (0,1->0; 1,2->1; ... 17,18->17)
1958                                         // the adc and pad-mapping should now be one to one: adc i is linked to pad i; TRAP-numbering
1959     Int_t madc = mADC[i];
1960     if (madc == -1) continue;
1961     
1962     YAlu.AssignInt(N[rpairChannel[i]]);
1963     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)
1964
1965     mN[i]    = hitSum[madc];
1966   
1967     // 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
1968     if (N[madc+1] == 0) {
1969         mQT0[i] = QT0[madc];
1970         mQT1[i] = QT1[madc];
1971         
1972     }
1973     else {
1974
1975         // is it ok to do the size-checking for the merged fit-sums with the same format as for single-channel fit-sums?
1976         
1977         mQT0[i]   = QT0[madc] + QT0[madc+1];
1978         QT0Alu.AssignFormatted(mQT0[i]);   
1979         QT0Alu  = QT0Alu;                // size-check
1980         mQT0[i]   = QT0Alu.GetValue();     // write back
1981         
1982         mQT1[i]   = QT1[madc] + QT1[madc+1];
1983         QT1Alu.AssignFormatted(mQT1[i]);
1984         QT1Alu  = QT1Alu;
1985         mQT1[i]   = QT1Alu.GetValue();
1986     }
1987     
1988     // calculate the mean charge in adc values; later to be replaced by electron likelihood
1989     mMeanCharge[i] = mQT0[i] + mQT1[i]; // total charge
1990     mMeanCharge[i] = mMeanCharge[i]>>2; // losing of accuracy; accounts for high mean charge
1991     // 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
1992     invN = 1.0/(mN[i]);
1993     inverseNAlu.AssignDouble(invN);
1994     inverseN = inverseNAlu.GetValue();
1995     mMeanCharge[i] = mMeanCharge[i] * inverseN;  // now to be interpreted with 8 past-comma bits
1996     TotalChargeAlu.AssignFormatted(mMeanCharge[i]);
1997     TotalChargeAlu = TotalChargeAlu;
1998     MeanChargeAlu = TotalChargeAlu;
1999     mMeanCharge[i] = MeanChargeAlu.GetValue();
2000     
2001     // this check is not necessary; it is just for efficiency reasons
2002     if (N[madc+1] == 0) {
2003         mX[i]     =   X[madc];
2004         mXX[i]    =  XX[madc];
2005         mY[i]     =   Y[madc];
2006         mXY[i]    =  XY[madc];
2007         mYY[i]    =  YY[madc];
2008     }
2009     else {
2010         
2011         mX[i]     =   X[madc] +  X[madc+1];
2012         XAlu.AssignFormatted(mX[i]);
2013         XAlu      = XAlu;
2014         mX[i]     = XAlu.GetValue();
2015         
2016         mXX[i]    =  XX[madc] + XX[madc+1];
2017         XXAlu.AssignFormatted(mXX[i]);
2018         XXAlu     = XXAlu;
2019         mXX[i]    = XXAlu.GetValue();
2020  
2021     
2022         mY[i]     =   Y[madc] + Y[madc+1] + wpad;
2023         if (mY[i] < 0) {
2024             YAlu.AssignFormatted(-mY[i]);
2025             YAlu.SetSign(-1);
2026         }
2027         else {
2028             YAlu.AssignFormatted(mY[i]);
2029             YAlu.SetSign(1);
2030         }
2031         YAlu    = YAlu;
2032         mY[i]     = YAlu.GetSignedValue();
2033         
2034         mXY[i]    = XY[madc] + XY[madc+1] + X[madc+1]*one;    // multiplication by one to maintain the data format
2035         
2036         if (mXY[i] < 0) {
2037             XYAlu.AssignFormatted(-mXY[i]);
2038             XYAlu.SetSign(-1);
2039         }
2040         else {
2041             XYAlu.AssignFormatted(mXY[i]);
2042             XYAlu.SetSign(1);
2043         }
2044         XYAlu   = XYAlu;
2045         mXY[i]    = XYAlu.GetSignedValue();
2046     
2047         mYY[i]    = YY[madc] + YY[madc+1] + 2*Y[madc+1]*one+ wpad*one;
2048         if (mYY[i] < 0) {
2049             YYAlu.AssignFormatted(-mYY[i]);
2050             YYAlu.SetSign(-1);
2051         }
2052         else {
2053             YYAlu.AssignFormatted(mYY[i]);
2054             YYAlu.SetSign(1);
2055         }
2056         
2057         YYAlu   = YYAlu;
2058         mYY[i]    = YYAlu.GetSignedValue();
2059     }
2060   
2061   }
2062     
2063   // calculation of offset and slope from the merged fit-sums; 
2064   // YY is needed for some error measure only; still to be done
2065   // be aware that all values are relative values (scale: timebin-width; pad-width) and are integer values on special scale
2066   
2067   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2068   // !!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   !!
2069   // !!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.    !!
2070   // !!This has to be taken into account by the GTU. Furthermore a Lorentz correction might have to be applied to the offset (see below).                                                                                             !!
2071   // !!In this implementation it is assumed that no miscalibration containing changing drift velocities in the amplification region is used.                                                                                          !!
2072   // !!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    !!
2073   // !!valid (in approximation) if tFS is close to the beginning of the drift region.                                                                                                                                                 !!
2074   // !!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                                                                    !!
2075   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2076
2077   // which formats should be chosen?
2078   AliTRDtrapAlu denomAlu;
2079   denomAlu.Init(20,8);       
2080   AliTRDtrapAlu numAlu;
2081   numAlu.Init(20,8);     
2082   // is this enough pre-comma place? covers the range of the 13 bit-word of the transmitted offset
2083   // offset measured in coord. of left channel must be between -0.5 and 1.5; 14 pre-comma bits because numerator can be big
2084
2085   for (Int_t i = 0; i<4; i++) {
2086     if (mADC[i] == -1) continue;
2087       
2088     Int_t num0  = (mN[i]*mXX[i]-mX[i]*mX[i]);
2089     if (num0 < 0) {
2090       denomAlu.AssignInt(-num0);    // num0 does not have to be interpreted as having past-comma bits -> AssignInt
2091       denomAlu.SetSign(-1);
2092     }
2093     else {
2094       denomAlu.AssignInt(num0);
2095       denomAlu.SetSign(1);
2096     }
2097     
2098     Int_t num1  = mN[i]*mXY[i] - mX[i]*mY[i];
2099     if (num1 < 0) {
2100       numAlu.AssignFormatted(-num1); // value of num1 is already formatted to have 8 past-comma bits
2101       numAlu.SetSign(-1);
2102     }
2103     else {
2104       numAlu.AssignFormatted(num1);
2105       numAlu.SetSign(1);
2106     }
2107     numAlu    = numAlu/denomAlu;
2108     mSlope[i]   = numAlu.GetSignedValue();
2109    
2110     Int_t num2  = mXX[i]*mY[i] - mX[i]*mXY[i];
2111    
2112     if (num2 < 0) {
2113       numAlu.AssignFormatted(-num2);
2114       numAlu.SetSign(-1);
2115     }
2116     else {
2117       numAlu.AssignFormatted(num2);
2118       numAlu.SetSign(1);
2119     }
2120    
2121     numAlu    = numAlu/denomAlu;
2122    
2123     
2124     mOffset[i]  = numAlu.GetSignedValue();
2125     numAlu.SetSign(1);
2126     denomAlu.SetSign(1);
2127        
2128                                  
2129     //numAlu.AssignInt(mADC[i]+1);   // according to TRAP-manual but trafo not to middle of chamber (0.5 channels away)             
2130     numAlu.AssignDouble((Double_t)mADC[i] + 1.5);      // numAlu has enough pre-comma place for that; correct trafo, best values
2131     mOffset[i]  = mOffset[i] + numAlu.GetValue();      // transform offset to a coord.system relative to chip; +1 to avoid neg. values 
2132     
2133     // up to here: adc-mapping according to TRAP manual and in line with pad-col mapping
2134     // reverse adc-counting to be again in line with the online mapping
2135     mADC[i]     = fNADC - 4 - mADC[i];                 // fNADC-4-mADC[i]: 0..17; remapping necessary;
2136     mADC[i]     = mADC[i] + 2; 
2137     // +2: mapping onto original ADC-online-counting: inner adc's corresponding to a chip's pasa: number 2..19
2138   }
2139
2140   // adc-counting is corresponding to online mapping; use AliTRDfeeParam::GetPadColFromADC to get the pad to which adc is connected; 
2141   // pad-column mapping is reverse to adc-online mapping; TRAP adc-mapping is in line with pad-mapping (increase in same direction);
2142   
2143   // transform parameters to the local coordinate-system of a stack (used by GTU)
2144   AliTRDpadPlane* padPlane = fGeo->CreatePadPlane(fLayer,fStack);
2145   
2146   Double_t padWidthI = padPlane->GetWidthIPad()*10.0; // get values in cm; want them in mm
2147   //Double_t padWidthO = padPlane->GetWidthOPad()*10; // difference between outer pad-widths not included; in real TRAP??
2148   
2149   // difference between width of inner and outer pads of a row is not accounted for;
2150   
2151   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 
2152   Double_t eCharge  = 0.3;                           // unit charge in (GeV/c)/m*T
2153   Double_t ptMin   = 2.3;                            // minimum transverse momentum (GeV/c); to be adjusted(?)
2154   
2155   Double_t granularityOffset = 0.160;                // granularity for offset in mm
2156   Double_t granularitySlope  = 0.140;                // granularity for slope  in mm     
2157     
2158   // get the coordinates in SM-system; parameters: 
2159   
2160   Double_t zPos       =  (padPlane->GetRowPos(fRow))*10.0;  // z-position of the MCM; fRow is counted on a chamber; SM consists of 5 
2161   // zPos is position of pad-borders;
2162   Double_t zOffset = 0.0;
2163   if ( fRow == 0 || fRow == 15 ) {
2164       zOffset = padPlane->GetLengthOPad();
2165   }
2166   else {
2167       zOffset = padPlane->GetLengthIPad();
2168   }
2169   zOffset = (-1.0) * zOffset/2.0;
2170   // turn zPos to be z-coordinate at middle of pad-row
2171   zPos = zPos + zOffset;
2172
2173       
2174   Double_t xPos       =  0.0;                               // x-position of the upper border of the drift-chamber of actual layer
2175   Int_t    icol       =  0;                                 // column-number of adc-channel
2176   Double_t yPos[4];                                         // y-position of the pad to which ADC is connected
2177   Double_t dx         = 30.0;                               // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry
2178   Double_t freqSample = fFeeParam->GetSamplingFrequency();  // retrieve the sampling frequency (10.019750 MHz)
2179   Double_t vdrift     = fCal->GetVdriftAverage(fChaId);     // averaged drift velocity for this detector (1.500000 cm/us)
2180   Int_t    nrOfDriftTimeBins = Int_t(dx/10.0*freqSample/vdrift); // the number of time bins in the drift region (20)
2181   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))
2182   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
2183   if(nrOfOffsetCorrTimeBins < 0) nrOfOffsetCorrTimeBins = 0;// don't apply offset correction if no drift time bins before tFS can be assumed 
2184   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
2185   //Double_t lorAngle   =  7.0;                             // Lorentz-angle in degrees
2186   Double_t tiltAngle  = padPlane->GetTiltingAngle();        // sign-respecting tilting angle of pads in actual layer
2187   Double_t tiltTan    = TMath::Tan(TMath::Pi()/180.0 * tiltAngle);
2188   //Double_t lorTan     = TMath::Tan(TMath::Pi()/180.0 * lorAngle);
2189
2190   Double_t alphaMax[4];                            // maximum deflection from the direction to the primary vertex; granularity of hit pads
2191   Double_t slopeMin[4];                            // local limits for the deflection
2192   Double_t slopeMax[4];
2193   Int_t   mslopeMin[4];                            // in granularity units; to be compared to mSlope[i]
2194   Int_t   mslopeMax[4];
2195
2196
2197   // x coord. of upper side of drift chambers in local SM-system (in mm)
2198   // obtained by evaluating the x-range of the hits; should be crosschecked; only drift, not amplification region taken into account (30mm);
2199   // the y-deflection is given as difference of y between lower and upper side of drift-chamber, not pad-plane;
2200   switch(fLayer) 
2201     {
2202     case 0: 
2203       xPos = 3003.0;
2204       break;
2205     case 1:
2206       xPos = 3129.0;
2207       break;
2208     case 2:
2209       xPos = 3255.0;
2210       break;
2211     case 3:
2212       xPos = 3381.0;
2213       break;
2214     case 4:
2215       xPos = 3507.0;
2216       break;
2217     case 5:
2218       xPos = 3633.0;
2219       break;
2220     }
2221  
2222   // calculation of offset-correction n: 
2223
2224   Int_t nCorrectOffset = (fRobPos % 2 == 0) ? ((fMcmPos % 4)) : ( 4 + (fMcmPos % 4));  
2225  
2226   nCorrectOffset = (nCorrectOffset - 4)*18 - 1;
2227   if (nCorrectOffset < 0) {
2228     numAlu.AssignInt(-nCorrectOffset);
2229     numAlu.SetSign(-1);
2230   }
2231   else {
2232     numAlu.AssignInt(nCorrectOffset);
2233     numAlu.SetSign(1);
2234   }
2235   nCorrectOffset = numAlu.GetSignedValue();   
2236
2237   // the Lorentz correction to the offset
2238   Double_t lorCorrectOffset = lorTan *(Double_t)nrOfOffsetCorrTimeBins*vdrift*10.0/freqSample; // Lorentz offset correction in mm
2239   
2240
2241   lorCorrectOffset = lorCorrectOffset/padWidthI; // Lorentz correction in pad width units
2242   
2243   if(lorCorrectOffset < 0) {
2244       numAlu.AssignDouble(-lorCorrectOffset);
2245       numAlu.SetSign(-1);
2246   }
2247   else{
2248       numAlu.AssignDouble(lorCorrectOffset);
2249       numAlu.SetSign(1);
2250   }
2251   
2252   Int_t mlorCorrectOffset = numAlu.GetSignedValue();
2253   
2254   
2255   Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0
2256  
2257   // calculation of slope-correction
2258
2259   // this is only true for tracks coming (approx.) from primary vertex
2260   // everything is evaluated for a tracklet covering the whole drift chamber
2261   Double_t cCorrectSlope = (-lorTan*dx + zPos/xPos*dx*tiltTan)/granularitySlope;
2262   // Double_t cCorrectSlope =  zPos/xPos*dx*tiltTan/granularitySlope;
2263   // zPos can be negative! for track from primary vertex: zOut-zIn > 0 <=> zPos > 0
2264   
2265   if (cCorrectSlope < 0) {
2266       numAlu.AssignDouble(-cCorrectSlope);
2267       numAlu.SetSign(-1);
2268   }
2269   else {
2270       numAlu.AssignDouble(cCorrectSlope);
2271       numAlu.SetSign(1);
2272   }
2273   cCorrectSlope = numAlu.GetSignedValue();
2274  
2275   // convert slope to deflection between upper and lower drift-chamber position (slope is given in pad-unit/time-bins)
2276   // different pad-width of outer pads of a pad-plane not taken into account
2277   // 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)
2278   // 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
2279   
2280   
2281   Double_t mCorrectSlope = (Double_t)(nrOfDriftTimeBins)*padWidthI/granularitySlope;  // >= 0.0
2282
2283   AliTRDtrapAlu correctAlu;
2284   correctAlu.Init(20,8);
2285   
2286   AliTRDtrapAlu offsetAlu;
2287   offsetAlu.Init(13,0,-0x1000,0x0FFF);          // 13 bit-word; 2-complement (1 sign-bit); asymmetric range
2288   
2289   AliTRDtrapAlu slopeAlu;
2290   slopeAlu.Init(7,0,-0x40,0x3F);                // 7 bit-word;  2-complement (1 sign-bit);
2291
2292   for (Int_t i = 0; i<4; i++) {
2293     
2294     if (mADC[i] == -1) continue;
2295     
2296     icol = fFeeParam->GetPadColFromADC(fRobPos,fMcmPos,mADC[i]); // be aware that mADC[i] contains the ADC-number according to online-mapping
2297     yPos[i]   = (padPlane->GetColPos(icol))*10.0;
2298     
2299     
2300     // offset:
2301     
2302     correctAlu.AssignDouble(mCorrectOffset);     // done because max. accuracy is 8 bit
2303     mCorrectOffset = correctAlu.GetValueWhole(); // cut offset correction to 8 past-comma bit accuracy
2304     mOffset[i]  = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset - mlorCorrectOffset)); 
2305     //mOffset[i]  = mOffset[i]*(-1);                   // adjust to direction of y-axes in online simulation
2306     
2307     if (mOffset[i] < 0) {
2308       numAlu.AssignFormatted(-mOffset[i]);
2309       numAlu.SetSign(-1);
2310     }
2311     else {
2312       numAlu.AssignFormatted(mOffset[i]);
2313       numAlu.SetSign(1);
2314     }
2315
2316     offsetAlu = numAlu; 
2317     mOffset[i]  = offsetAlu.GetSignedValue();  
2318
2319     
2320     // slope:
2321     
2322     correctAlu.AssignDouble(mCorrectSlope);
2323     mCorrectSlope = correctAlu.GetValueWhole();
2324     
2325     mSlope[i]   = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope);
2326
2327     if (mSlope[i] < 0) {
2328       numAlu.AssignFormatted(-mSlope[i]);
2329       numAlu.SetSign(-1);
2330     }
2331     else {
2332       numAlu.AssignFormatted(mSlope[i]);
2333       numAlu.SetSign(1);
2334     }
2335
2336     slopeAlu  = numAlu;     // here all past-comma values are cut, not rounded; alternatively add +0.5 before cutting (means rounding)
2337     mSlope[i]   = slopeAlu.GetSignedValue(); 
2338        
2339     // local (LTU) limits for the deflection 
2340     // ATan returns angles in radian
2341     alphaMax[i]  = TMath::ASin(eCharge*magField/(2.0*ptMin)*(TMath::Sqrt(xPos*xPos + yPos[i]*yPos[i]))/1000.0); // /1000: mm->m
2342     slopeMin[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) - alphaMax[i]))/granularitySlope;
2343     slopeMax[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) + alphaMax[i]))/granularitySlope;
2344     
2345     if (slopeMin[i] < 0) {
2346       slopeAlu.AssignDouble(-slopeMin[i]);
2347       slopeAlu.SetSign(-1);
2348     }
2349     else { 
2350       slopeAlu.AssignDouble(slopeMin[i]);
2351       slopeAlu.SetSign(1);
2352     }
2353     mslopeMin[i] = slopeAlu.GetSignedValue();  // the borders should lie inside the range of mSlope -> usage of slopeAlu again
2354    
2355     if (slopeMax[i] < 0) {
2356       slopeAlu.AssignDouble(-slopeMax[i]);
2357       slopeAlu.SetSign(-1);
2358     }
2359     else {
2360       slopeAlu.AssignDouble(slopeMax[i]);
2361       slopeAlu.SetSign(1);
2362     }
2363     mslopeMax[i] = slopeAlu.GetSignedValue();
2364   }
2365
2366   // suppress submission of tracks with low stiffness
2367   // put parameters in 32bit-word and submit (write to file as root-file; sort after SM, stack, layer, chamber) 
2368
2369   // sort tracklet-words in ascending y-order according to the offset (according to mADC would also be possible)
2370   // up to now they are sorted according to maximum hit sum
2371   // is the sorting really done in the TRAP-chip?
2372   
2373   Int_t order[4] = {-1,-1,-1,-1};
2374   Int_t wordnr = 0;   // number of tracklet-words
2375   
2376   for(Int_t j = 0; j < fMaxTracklets; j++) {
2377       //if( mADC[j] == -1) continue; 
2378       if( (mADC[j] == -1) || (mSlope[j] < mslopeMin[j]) || (mSlope[j] > mslopeMax[j])) continue; // this applies a pt-cut
2379       wordnr++;
2380       if( wordnr-1 == 0) {
2381           order[0] = j;
2382           continue;
2383       }
2384       // wordnr-1>0, wordnr-1<4
2385       order[wordnr-1] = j;
2386       for( Int_t k = 0; k < wordnr-1; k++) {
2387           if( mOffset[j] < mOffset[order[k]] ) {
2388               for( Int_t l = wordnr-1; l > k; l-- ) {
2389                   order[l] = order[l-1];
2390               }
2391               order[k] = j;
2392               break;
2393           }
2394           
2395       }
2396   }
2397         
2398   // fill the bit-words in ascending order and without gaps
2399   UInt_t bitWord[4] = {0,0,0,0};                 // attention: unsigned int to have real 32 bits (no 2-complement)
2400   for(Int_t j = 0; j < wordnr; j++) { // only "wordnr" tracklet-words
2401       //Bool_t rem1 = kTRUE;
2402     
2403     Int_t i = order[j];
2404     //bit-word is 2-complement and therefore without sign
2405     bitWord[j] =   1; // this is the starting 1 of the bit-word (at 33rd position); the 1 must be ignored
2406     //printf("\n");
2407     UInt_t shift  = 0;
2408     UInt_t shift2 = 0;
2409         
2410         
2411
2412
2413     /*printf("mean charge: %d\n",mMeanCharge[i]);
2414     printf("row: %d\n",fRow);
2415     printf("slope: %d\n",mSlope[i]);
2416     printf("pad position: %d\n",mOffset[i]);
2417     printf("channel: %d\n",mADC[i]);*/
2418
2419     // electron probability (currently not implemented; the mean charge is just scaled)
2420     shift = (UInt_t)mMeanCharge[i];
2421     for(Int_t iBit = 0; iBit < 8; iBit++) {
2422       bitWord[j]  = bitWord[j]<<1;
2423       bitWord[j] |= (shift>>(7-iBit))&1;               
2424       //printf("0");
2425     }
2426
2427     // pad row
2428     shift = (UInt_t)fRow;
2429     for(Int_t iBit = 0; iBit < 4; iBit++) {
2430       bitWord[j]  = bitWord[j]<<1;
2431       bitWord[j] |= (shift>>(3-iBit))&1;
2432       //printf("%d", (fRow>>(3-iBit))&1);
2433     }
2434     
2435     // deflection length
2436     if(mSlope[i] < 0) {
2437         shift = (UInt_t)(-mSlope[i]);
2438         // shift2 is 2-complement of shift
2439         shift2 = 1;
2440         for(Int_t iBit = 1; iBit < 7; iBit++) {
2441             shift2  = shift2<<1;
2442             shift2 |= (1- (((shift)>>(6-iBit))&1) );
2443             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2444         }
2445         shift2 = shift2 + 1;
2446         //printf("1");
2447         for(Int_t iBit = 0; iBit < 7; iBit++) {
2448             bitWord[j]  = bitWord[j]<<1;
2449             bitWord[j] |= (shift2>>(6-iBit))&1;
2450             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2451         }
2452     }
2453     else {
2454         shift = (UInt_t)(mSlope[i]);
2455         bitWord[j]  = bitWord[j]<<1;
2456         bitWord[j]   |= 0;
2457         //printf("0");
2458         for(Int_t iBit = 1; iBit < 7; iBit++) {
2459             bitWord[j]  = bitWord[j]<<1;
2460             bitWord[j] |= (shift>>(6-iBit))&1;
2461             //printf("%d",(mSlope[i]>>(6-iBit))&1);
2462         }
2463     }
2464
2465     // pad position
2466     if(mOffset[i] < 0) {
2467         shift = (UInt_t)(-mOffset[i]);
2468         shift2 = 1;
2469         for(Int_t iBit = 1; iBit < 13; iBit++) {
2470             shift2  = shift2<<1;
2471             shift2 |= (1-(((shift)>>(12-iBit))&1));
2472             //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
2473         }
2474         shift2 = shift2 + 1;
2475         //printf("1");
2476         for(Int_t iBit = 0; iBit < 13; iBit++) {
2477             bitWord[j]  = bitWord[j]<<1;
2478             bitWord[j] |= (shift2>>(12-iBit))&1;
2479             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2480         }
2481     }
2482     else {
2483         shift = (UInt_t)mOffset[i];
2484         bitWord[j] = bitWord[j]<<1;
2485         bitWord[j]   |= 0; 
2486         //printf("0");
2487         for(Int_t iBit = 1; iBit < 13; iBit++) {
2488             bitWord[j]  = bitWord[j]<<1;
2489             bitWord[j] |= (shift>>(12-iBit))&1;
2490             //printf("%d",(mOffset[i]>>(12-iBit))&1);
2491         }
2492     }
2493
2494
2495         
2496     //printf("bitWord: %u\n",bitWord[j]);
2497     //printf("adc: %d\n",mADC[i]);
2498     fMCMT[j] = bitWord[j];
2499   }
2500     
2501   //printf("\n");
2502
2503   
2504   delete [] qsum;
2505   delete [] ieffped;
2506
2507   delete [] X;
2508   delete [] XX;
2509   delete [] Y;
2510   delete [] YY;
2511   delete [] XY;
2512   delete [] N;
2513   delete [] QT0;
2514   delete [] QT1;
2515
2516   delete [] hitSum;
2517   delete [] selectPair;
2518
2519   delete padPlane;
2520
2521 //if you want to activate the MC tracklet output, set fgkMCTrackletOutput=kTRUE in AliTRDfeeParam
2522         
2523   if (!fFeeParam->GetMCTrackletOutput()) 
2524       return;
2525  
2526   AliLog::SetClassDebugLevel("AliTRDmcmSim", 10);
2527   AliLog::SetFileOutput("../log/tracklet.log");
2528   
2529   // testing for wordnr in order to speed up the simulation
2530   if (wordnr == 0) 
2531     return;
2532    
2533   UInt_t        *trackletWord = new UInt_t[fMaxTracklets];
2534   Int_t         *adcChannel   = new Int_t[fMaxTracklets];
2535   Int_t         *trackRef     = new Int_t[fMaxTracklets];
2536
2537   Int_t u = 0;
2538
2539   AliTRDdigitsManager *digman = new AliTRDdigitsManager();
2540   digman->ReadDigits(gAlice->GetRunLoader()->GetLoader("TRDLoader")->TreeD());
2541   digman->SetUseDictionaries(kTRUE);
2542   AliTRDfeeParam *feeParam = AliTRDfeeParam::Instance();
2543
2544   for (Int_t j = 0; j < fMaxTracklets; j++) {
2545       Int_t i = order[j];
2546       trackletWord[j] = 0;
2547       adcChannel[j] = -1;
2548       if (bitWord[j]!=0) {
2549           trackletWord[u] = bitWord[j];
2550           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
2551
2552 // Finding label of MC track
2553           TH1F *hTrkRef = new TH1F("trackref", "trackref", 100000, 0, 100000);
2554           Int_t track[3];
2555           Int_t padcol = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u]);
2556           Int_t padcol_ngb = feeParam->GetPadColFromADC(fRobPos, fMcmPos, adcChannel[u] - 1);
2557           Int_t padrow = 4 * (fRobPos / 2) + fMcmPos / 4;
2558           Int_t det = 30 * fSector + 6 * fStack + fLayer;
2559           for(Int_t iTimebin = feeParam->GetLinearFitStart(); iTimebin < feeParam->GetLinearFitEnd(); iTimebin++) {
2560               track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2561               track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2562               track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2563               hTrkRef->Fill(track[0]);
2564               if (track[1] != track[0] && track[1] != -1)
2565                   hTrkRef->Fill(track[1]);
2566               if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2567                   hTrkRef->Fill(track[2]);
2568               if (padcol_ngb >= 0) {
2569                   track[0] = digman->GetTrack(0, padrow, padcol, iTimebin, det);
2570                   track[1] = digman->GetTrack(1, padrow, padcol, iTimebin, det);
2571                   track[2] = digman->GetTrack(2, padrow, padcol, iTimebin, det);
2572                   hTrkRef->Fill(track[0]);
2573                   if (track[1] != track[0] && track[1] != -1)
2574                       hTrkRef->Fill(track[1]);
2575                   if (track[2] != track[0] && track[2] != track[1] && track[2] != -1)
2576                       hTrkRef->Fill(track[2]);
2577               }
2578           }
2579           trackRef[u] = hTrkRef->GetMaximumBin() - 1;
2580           delete hTrkRef;
2581           u = u + 1;
2582       }
2583   }
2584
2585  AliDataLoader *dl = gAlice->GetRunLoader()->GetLoader("TRDLoader")->GetDataLoader("tracklets");
2586   if (!dl) {
2587     AliError("Could not get the tracklets data loader!");
2588   }
2589   else {
2590     TTree *trackletTree = dl->Tree();
2591     if (!trackletTree)
2592       dl->MakeTree();
2593     trackletTree = dl->Tree();
2594
2595    AliTRDtrackletMCM *trkl = new AliTRDtrackletMCM(); 
2596    TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
2597    if (!trkbranch)
2598        trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
2599     trkbranch->SetAddress(&trkl);
2600
2601     for (Int_t iTracklet = 0; iTracklet < fMaxTracklets; iTracklet++) {
2602         if (trackletWord[iTracklet] == 0)
2603             continue;
2604         trkl->SetTrackletWord(trackletWord[iTracklet]);
2605         trkl->SetDetector(30*fSector + 6*fStack + fLayer);
2606         trkl->SetROB(fRobPos);
2607         trkl->SetMCM(fMcmPos);
2608         trkl->SetLabel(trackRef[iTracklet]);
2609         trackletTree->Fill();
2610     }
2611     delete trkl;
2612     dl->WriteData("OVERWRITE");
2613   }
2614
2615   delete [] trackletWord;
2616   delete [] adcChannel; 
2617   delete [] trackRef;
2618   delete digman;
2619
2620   // to be done:
2621   // error measure for quality of fit (not necessarily needed for the trigger)
2622   // cluster quality threshold (not yet set)
2623   // electron probability
2624 }
2625