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