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