Replace plane by layer and chamber by stack
[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 ) // 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
1259 // in order to use the Tracklets, please first 
1260 // -- set AliTRDfeeParam::fgkTracklet to kTRUE, in order to switch on Tracklet-calculation
1261 // -- set AliTRDfeeParam::fgkTFtype   to 3, in order to use the new tail cancellation filter
1262 // currently tracklets from filtered digits are only given when setting fgkTFtype (AliTRDfeeParam) to 3
1263
1264 // 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 outside 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-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 be in AliTRDfeeParam and 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   
1425   AliTRDtrapAlu inverseNAlu;
1426   inverseNAlu.Init(1,8);        // replaces the LUT for 1/N
1427   AliTRDtrapAlu MeanChargeAlu;  // mean charge in ADC counts
1428   MeanChargeAlu.Init(8,0);
1429   AliTRDtrapAlu TotalChargeAlu;
1430   TotalChargeAlu.Init(17,8);
1431   //nr of post comma bits should be the same for inverseN and TotalCharge
1432   
1433   
1434   SetPosLUT();                    // initialize the position correction LUT for this MCM;
1435
1436
1437   // fit-sums; remapping!; 0,1,2->0; 1,2,3->1; ... 18,19,20->18
1438   Int_t *X   = new Int_t[fNADC-2];
1439   Int_t *XX  = new Int_t[fNADC-2];
1440   Int_t *Y   = new Int_t[fNADC-2];
1441   Int_t *YY  = new Int_t[fNADC-2];
1442   Int_t *XY  = new Int_t[fNADC-2];
1443   Int_t *N   = new Int_t[fNADC-2];
1444   Int_t *QT0 = new Int_t[fNADC-2]; // accumulated charge
1445   Int_t *QT1 = new Int_t[fNADC-2]; // accumulated charge
1446   
1447   for (Int_t iCol = 0; iCol < fNADC-2; iCol++) { 
1448       
1449       // initialize fit-sums 
1450       X[iCol]   = 0;
1451       XX[iCol]  = 0;
1452       Y[iCol]   = 0;
1453       YY[iCol]  = 0;
1454       XY[iCol]  = 0;
1455       N[iCol]   = 0;
1456       QT0[iCol] = 0;
1457       QT1[iCol] = 0;
1458   }
1459   
1460
1461   filPed.Init(7,2);                         // convert filtered pedestal into 7+2Bits
1462   
1463   for(Int_t iT = 0; iT < fNTimeBin; iT++){
1464     
1465     if(iT<tFS || iT>=tFE) continue;         // linear fit yes/no? // !!**enable**!!
1466
1467     // reset
1468     Int_t portChannel[4]   = {-1,-1,-1,-1};   
1469     Int_t clusterCharge[4] = {0,0,0,0};
1470     Int_t leftCharge[4]    = {0,0,0,0};
1471     Int_t centerCharge[4]  = {0,0,0,0}; 
1472     Int_t rightCharge[4]   = {0,0,0,0};
1473     
1474     Int_t mark = 0;
1475     
1476     filPed.AssignFormatted(ieffped[iT]);   // no size-checking with AssignFormatted; ieffped>=0
1477     filPed = filPed;                       // this checks the size
1478     
1479     ieffped[iT] = filPed.GetValue();
1480         
1481     for(Int_t i = 0; i<7; i++){
1482       selection[i].next       = NULL;
1483       selection[i].iadc       =   -1;     // value of -1: invalid adc
1484       selection[i].value      =    0;
1485    
1486     }
1487     // selection[0] is starting list-element; just for pointing
1488
1489     // loop over inner adc's 
1490     for (Int_t iCol = 1; iCol < fNADC-1; iCol++) { 
1491       
1492       Int_t left   = fADCT[iCol-1][iT]; 
1493       Int_t center = fADCT[iCol][iT];
1494       Int_t right  = fADCT[iCol+1][iT];  
1495
1496       Int_t sum = left + center + right;            // cluster charge sum
1497       qsumAlu.AssignFormatted(sum);    
1498       qsumAlu = qsumAlu;                        // size-checking; redundant
1499  
1500       qsum[iCol] = qsumAlu.GetValue(); 
1501       
1502       //hit detection and masking
1503       if(center>=left){
1504         if(center>right){
1505           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
1506             mark |= 1;                              // marker
1507           }
1508         }
1509       }
1510       mark = mark<<1;                
1511     }
1512     mark = mark>>1;
1513
1514        
1515     // get selection of 6 adc's and sort,starting with greatest values
1516
1517     //read three from right side and sort (primitive sorting algorithm)
1518     Int_t i = 0; // adc number
1519     Int_t j = 1; // selection number
1520     while(i<fNADC-2 && j<=3){
1521       i = i + 1;
1522       if((mark>>(i-1)) & 1 == 1) {
1523         selection[j].iadc  = fNADC-1-i;
1524         selection[j].value = qsum[fNADC-1-i]>>6;   // for hit-selection only the first 8 out of the 14 Bits are used for comparison
1525         
1526         // insert into sorted list
1527         listLeft = &selection[0];
1528         list = listLeft->next;
1529         
1530         if(list!=NULL) {
1531           while((list->next != NULL) && (selection[j].value <= list->value)){
1532             listLeft = list;
1533             list = list->next;
1534           }
1535           
1536           if(selection[j].value<=list->value){
1537             selection[j].next = list->next;
1538             list->next = &selection[j];
1539           }
1540           else {
1541             listLeft->next = &selection[j];
1542             selection[j].next = list;
1543           }
1544         }
1545         else{
1546           listLeft->next = &selection[j];
1547           selection[j].next = list;
1548         }
1549         
1550         j = j + 1;
1551       }
1552     }
1553
1554
1555     // read three from left side
1556     Int_t k = fNADC-2;
1557     while(k>i && j<=6) {
1558       if((mark>>(k-1)) & 1 == 1) {
1559         selection[j].iadc  = fNADC-1-k;
1560         selection[j].value = qsum[fNADC-1-k]>>6;
1561         
1562         listLeft = &selection[0];
1563         list = listLeft->next;
1564         
1565         if(list!=NULL){
1566           while((list->next != NULL) && (selection[j].value <= list->value)){
1567             listLeft = list;
1568             list = list->next;
1569           }
1570         
1571           if(selection[j].value<=list->value){
1572             selection[j].next = list->next;
1573             list->next = &selection[j];
1574           }
1575           else {
1576             listLeft->next = &selection[j];
1577             selection[j].next = list;
1578           }
1579         }
1580         else{
1581           listLeft->next = &selection[j];
1582           selection[j].next = list;
1583         }
1584
1585         j = j + 1;
1586       }
1587       k = k - 1;
1588     }
1589
1590     // get the four with greatest charge-sum
1591     list = &selection[0];
1592     for(i = 0; i<4; i++){
1593       if(list->next == NULL) continue;
1594       list = list->next;
1595       if(list->iadc == -1) continue;
1596       Int_t adc = list->iadc;                              // channel number with selected hit
1597       
1598       // the following arrays contain the four chosen channels in 1 time-bin
1599       portChannel[i]   = adc; 
1600       clusterCharge[i] = qsum[adc];
1601       leftCharge[i]    = fADCT[adc-1][iT] - ieffped[iT]; // reduce by filtered pedestal (pedestal is part of the signal)
1602       centerCharge[i]  = fADCT[adc][iT] - ieffped[iT];           
1603       rightCharge[i]   = fADCT[adc+1][iT] - ieffped[iT];         
1604     }
1605
1606     // arithmetic unit
1607     
1608     // cluster verification
1609     if(!bVBY){
1610       for(i = 0; i<4; i++){
1611         Int_t lr = leftCharge[i]*rightCharge[i]*1024;
1612         Int_t cc = centerCharge[i]*centerCharge[i]*cQT;
1613         if (lr>=cc){
1614           portChannel[i]   = -1;                                 // set to invalid address 
1615           clusterCharge[i] = 0;
1616         }
1617       }
1618     }
1619
1620     // fit-sums of valid channels
1621     // local hit position
1622     for(i = 0; i<4; i++){
1623       if (centerCharge[i] ==  0) {
1624         portChannel[i] = -1; 
1625       }// prevent division by 0
1626       
1627       if (portChannel[i]  == -1) continue;
1628       
1629       Double_t dCOG = (Double_t)(rightCharge[i]-leftCharge[i])/centerCharge[i];
1630        
1631       Bool_t sign = (dCOG>=0.0) ? kFALSE : kTRUE;
1632       dCOG = (sign == kFALSE) ? dCOG : -dCOG;     // AssignDouble doesn't allow for signed doubles
1633       dCOGAlu.AssignDouble(dCOG);
1634       Int_t iLUTpos = dCOGAlu.GetValue();       // steers position in LUT
1635             
1636       dCOG = dCOG/2;
1637       yrawAlu.AssignDouble(dCOG);
1638       Int_t iCOG = yrawAlu.GetValue();
1639       Int_t y = iCOG + fPosLUT[iLUTpos % 128];    // local position in pad-units
1640       yrawAlu.AssignFormatted(y);               // 0<y<1           
1641       yAlu  = yrawAlu;                        // convert to 16 past-comma bits
1642       
1643       if(sign == kTRUE) yAlu.SetSign(-1);       // buffer width of 9 bits; sign on real (not estimated) position
1644       xAlu.AssignInt(iT);                       // buffer width of 5 bits 
1645       
1646
1647       xxAlu = xAlu * xAlu;                  // buffer width of 10 bits -> fulfilled by x*x       
1648       
1649       yyAlu = yAlu * yAlu;                  // buffer width of 16 bits
1650    
1651       xyAlu = xAlu * yAlu;                  // buffer width of 14 bits
1652                   
1653       Int_t adc = portChannel[i]-1;              // remapping! port-channel contains channel-nr. of inner adc's (1..19; mapped to 0..18)
1654
1655       // calculate fit-sums recursively
1656       // interpretation of their bit-length is given as comment
1657       
1658       // be aware that the accuracy of the result of a calculation is always determined by the accuracy of the less accurate value
1659
1660       XAlu.AssignFormatted(X[adc]);
1661       XAlu = XAlu + xAlu;                   // buffer width of 9 bits 
1662       X[adc] = XAlu.GetValue();
1663              
1664       XXAlu.AssignFormatted(XX[adc]);
1665       XXAlu = XXAlu + xxAlu;                // buffer width of 14 bits    
1666       XX[adc] = XXAlu.GetValue();
1667
1668       if (Y[adc] < 0) {
1669         YAlu.AssignFormatted(-Y[adc]);          // make sure that only positive values are assigned; sign-setting must be done by hand
1670         YAlu.SetSign(-1);
1671       }
1672       else {
1673         YAlu.AssignFormatted(Y[adc]);
1674         YAlu.SetSign(1);
1675       }
1676         
1677       YAlu = YAlu + yAlu;                   // buffer width of 14 bits (8 past-comma);     
1678       Y[adc] = YAlu.GetSignedValue();
1679             
1680       YYAlu.AssignFormatted(YY[adc]);
1681       YYAlu = YYAlu + yyAlu;                // buffer width of 21 bits (16 past-comma) 
1682       YY[adc] = YYAlu.GetValue();
1683            
1684       if (XY[adc] < 0) {
1685         XYAlu.AssignFormatted(-XY[adc]);
1686         XYAlu.SetSign(-1);
1687       }
1688       else {
1689         XYAlu.AssignFormatted(XY[adc]);
1690         XYAlu.SetSign(1);
1691       }
1692
1693       XYAlu = XYAlu + xyAlu;                // buffer allows 17 bits (8 past-comma) 
1694       XY[adc] = XYAlu.GetSignedValue();
1695             
1696       N[adc]  = N[adc] + 1;
1697    
1698
1699       // accumulated charge
1700       qsumAlu.AssignFormatted(qsum[adc+1]); // qsum was not remapped!
1701       qtruncAlu = qsumAlu;
1702
1703       if(iT>=tQS0 && iT<=tQE0){
1704         QT0Alu.AssignFormatted(QT0[adc]);
1705         QT0Alu = QT0Alu + qtruncAlu;
1706         QT0[adc] = QT0Alu.GetValue();
1707         //interpretation of QT0 as 12bit-value (all pre-comma); is this as it should be done?; buffer allows 15 Bit
1708       }
1709       
1710       if(iT>=tQS1 && iT<=tQE1){
1711         QT1Alu.AssignFormatted(QT1[adc]);
1712         QT1Alu = QT1Alu + qtruncAlu;
1713         QT1[adc] = QT1Alu.GetValue();
1714         //interpretation of QT1 as 12bit-value; buffer allows 16 Bit
1715       }
1716     }// i
1717       
1718     // remapping is done!!
1719      
1720   }//iT
1721  
1722   
1723     
1724   // tracklet-assembly
1725   
1726   // put into AliTRDfeeParam and take care that values are in proper range
1727   const Int_t cTCL = 1;      // left adc: number of hits; 8<=TCL<=31 (?? 1<=cTCL<+8 ??) 
1728   const Int_t cTCT = 8;      // joint number of hits;     8<=TCT<=31
1729   
1730   Int_t mPair   = 0;         // marker for possible tracklet pairs
1731   Int_t* hitSum = new Int_t[fNADC-3];
1732   // hitSum[0] means: hit sum of remapped channels 0 and 1; hitSum[17]: 17 and 18; 
1733   
1734   // check for all possible tracklet-pairs of adjacent channels (two are merged); mark the left channel of the chosen pairs
1735   for (Int_t iCol = 0; iCol < fNADC-3; iCol++) {
1736     hitSum[iCol] = N[iCol] + N[iCol+1];
1737     if ((N[iCol]>=cTCL) && (hitSum[iCol]>=cTCT)) {
1738         mPair |= 1;         // mark as possible channel-pair
1739      
1740     }
1741     mPair = mPair<<1;
1742   }
1743   mPair = mPair>>1;
1744   
1745   List_t* selectPair = new List_t[fNADC-2];      // list with 18 elements (0..18) containing the left channel-nr and hit sums
1746                                                  // selectPair[18] is starting list-element just for pointing
1747   for(Int_t k = 0; k<fNADC-2; k++){
1748       selectPair[k].next       = NULL;
1749       selectPair[k].iadc       =   -1;           // invalid adc
1750       selectPair[k].value      =    0;
1751    
1752     }
1753
1754  list = NULL;
1755  listLeft = NULL;
1756   
1757   // read marker and sort according to hit-sum
1758   
1759   Int_t adcL  = 0;            // left adc-channel-number (remapped)
1760   Int_t selNr = 0;            // current number in list
1761   
1762   // insert marked channels into list and sort according to hit-sum
1763   while(adcL < fNADC-3 && selNr < fNADC-3){
1764      
1765     if((mPair>>((fNADC-4)-(adcL))) & 1 == 1) {
1766       selectPair[selNr].iadc  = adcL;
1767       selectPair[selNr].value = hitSum[adcL];   
1768       
1769       listLeft = &selectPair[fNADC-3];
1770       list = listLeft->next;
1771         
1772       if(list!=NULL) {
1773         while((list->next != NULL) && (selectPair[selNr].value <= list->value)){
1774           listLeft = list;
1775           list = list->next;
1776         }
1777         
1778         if(selectPair[selNr].value <= list->value){
1779           selectPair[selNr].next = list->next;
1780           list->next = &selectPair[selNr];
1781         }
1782         else {
1783           listLeft->next = &selectPair[selNr];
1784           selectPair[selNr].next = list;
1785         }
1786         
1787       }
1788       else{
1789         listLeft->next = &selectPair[selNr];
1790         selectPair[selNr].next = list;
1791       }
1792       
1793       selNr = selNr + 1;
1794     }
1795     adcL = adcL + 1;
1796   }
1797   
1798   //select up to 4 channels with maximum number of hits
1799   Int_t lpairChannel[4] = {-1,-1,-1,-1}; // save the left channel-numbers of pairs with most hit-sum
1800   Int_t rpairChannel[4] = {-1,-1,-1,-1}; // save the right channel, too; needed for detecting double tracklets
1801   list = &selectPair[fNADC-3];
1802   
1803   for (Int_t i = 0; i<4; i++) {
1804     if(list->next == NULL) continue;
1805     list = list->next;
1806     if(list->iadc == -1) continue;
1807     lpairChannel[i] = list->iadc;        // channel number with selected hit
1808     rpairChannel[i] = lpairChannel[i]+1;
1809   }
1810   
1811   // avoid submission of double tracklets  
1812   for (Int_t i = 3; i>0; i--) {
1813     for (Int_t j = i-1; j>-1; j--) {
1814       if(lpairChannel[i] == rpairChannel[j]) {
1815         lpairChannel[i] = -1;
1816         rpairChannel[i] = -1;
1817         break;
1818       }
1819       if(rpairChannel[i] == lpairChannel[j]) {
1820         lpairChannel[i] = -1;
1821         rpairChannel[i] = -1;
1822         break;
1823       }
1824     }
1825   }
1826   
1827   // merging of the fit-sums of the remainig channels
1828   // assume same data-word-width as for fit-sums for 1 channel
1829   // relative scales!
1830   Int_t mADC[4];                      
1831   Int_t mN[4];
1832   Int_t mQT0[4];
1833   Int_t mQT1[4];
1834   Int_t mX[4];
1835   Int_t mXX[4];
1836   Int_t mY[4];
1837   Int_t mYY[4];
1838   Int_t mXY[4];
1839   Int_t mOffset[4];
1840   Int_t mSlope[4];
1841   Int_t mMeanCharge[4]; 
1842   Int_t inverseN = 0;
1843   Double_t invN = 0;
1844  
1845   for (Int_t i = 0; i<4; i++){
1846     mADC[i] = -1;                        // set to invalid number
1847     mN[i]   =  0;
1848     mQT0[i] =  0;
1849     mQT1[i] =  0;
1850     mX[i]   =  0;
1851     mXX[i]  =  0;
1852     mY[i]   =  0;
1853     mYY[i]  =  0;
1854     mXY[i]  =  0;
1855     mOffset[i] = 0;
1856     mSlope[i]  = 0;
1857     mMeanCharge[i] = 0;
1858   }
1859   
1860
1861   YAlu.AssignInt(1);
1862   Int_t wpad  = YAlu.GetValue();       // 1 with 8 past-comma bits
1863  
1864   for (Int_t i = 0; i<4; i++){
1865     
1866     mADC[i] = lpairChannel[i];          // mapping of merged sums to left channel nr. (0,1->0; 1,2->1; ... 17,18->17)
1867                                          // the adc and pad-mapping should now be one to one: adc i is linked to pad i; TRAP-numbering
1868     Int_t madc = mADC[i];
1869     if (madc == -1) continue;
1870     mN[i]    = hitSum[madc];
1871   
1872     // 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
1873     if (N[madc+1] == 0) {
1874         mQT0[i] = QT0[madc];
1875         mQT1[i] = QT1[madc];
1876         
1877     }
1878     else {
1879
1880         // is it ok to do the size-checking for the merged fit-sums with the same format as for single-channel fit-sums?
1881         
1882         mQT0[i]   = QT0[madc] + QT0[madc+1];
1883         QT0Alu.AssignFormatted(mQT0[i]);   
1884         QT0Alu  = QT0Alu;                // size-check
1885         mQT0[i]   = QT0Alu.GetValue();     // write back
1886         
1887         mQT1[i]   = QT1[madc] + QT1[madc+1];
1888         QT1Alu.AssignFormatted(mQT1[i]);
1889         QT1Alu  = QT1Alu;
1890         mQT1[i]   = QT1Alu.GetValue();
1891     }
1892     
1893     // calculate the mean charge in adc values; later to be replaced by electron likelihood
1894     mMeanCharge[i] = mQT0[i] + mQT1[i]; // total charge
1895     mMeanCharge[i] = mMeanCharge[i]>>2; // losing of accuracy; accounts for high mean charge
1896     // 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
1897     invN = 1.0/(mN[i]);
1898     inverseNAlu.AssignDouble(invN);
1899     inverseN = inverseNAlu.GetValue();
1900     mMeanCharge[i] = mMeanCharge[i] * inverseN;  // now to be interpreted with 8 past-comma bits
1901     TotalChargeAlu.AssignFormatted(mMeanCharge[i]);
1902     TotalChargeAlu = TotalChargeAlu;
1903     MeanChargeAlu = TotalChargeAlu;
1904     mMeanCharge[i] = MeanChargeAlu.GetValue();
1905     
1906     if (N[madc+1] == 0) {
1907         mX[i]     =   X[madc];
1908         mXX[i]    =  XX[madc];
1909         mY[i]     =   Y[madc];
1910         mXY[i]    =  XY[madc];
1911         mYY[i]    =  YY[madc];
1912     }
1913     else {
1914         
1915         mX[i]     =   X[madc] +  X[madc+1];
1916         XAlu.AssignFormatted(mX[i]);
1917         XAlu    = XAlu;
1918         mX[i]     = XAlu.GetValue();
1919         
1920         mXX[i]    =  XX[madc] + XX[madc+1];
1921         XXAlu.AssignFormatted(mXX[i]);
1922         XXAlu   = XXAlu;
1923         mXX[i]    = XXAlu.GetValue();
1924  
1925     
1926         mY[i]     =   Y[madc] + Y[madc+1] + wpad;
1927         if (mY[i] < 0) {
1928             YAlu.AssignFormatted(-mY[i]);
1929             YAlu.SetSign(-1);
1930         }
1931         else {
1932             YAlu.AssignFormatted(mY[i]);
1933             YAlu.SetSign(1);
1934         }
1935         YAlu    = YAlu;
1936         mY[i]     = YAlu.GetSignedValue();
1937         
1938         mXY[i]    = XY[madc] + XY[madc+1] + X[madc+1]*wpad;
1939         
1940         if (mXY[i] < 0) {
1941             XYAlu.AssignFormatted(-mXY[i]);
1942             XYAlu.SetSign(-1);
1943         }
1944         else {
1945             XYAlu.AssignFormatted(mXY[i]);
1946             XYAlu.SetSign(1);
1947         }
1948         XYAlu   = XYAlu;
1949         mXY[i]    = XYAlu.GetSignedValue();
1950     
1951         mYY[i]    = YY[madc] + YY[madc+1] + 2*Y[madc+1]*wpad + wpad*wpad;
1952         if (mYY[i] < 0) {
1953             YYAlu.AssignFormatted(-mYY[i]);
1954             YYAlu.SetSign(-1);
1955         }
1956         else {
1957             YYAlu.AssignFormatted(mYY[i]);
1958             YYAlu.SetSign(1);
1959         }
1960         
1961         YYAlu   = YYAlu;
1962         mYY[i]    = YYAlu.GetSignedValue();
1963     }
1964   
1965   }
1966     
1967   // calculation of offset and slope from the merged fit-sums; 
1968   // YY is needed for some error measure only; still to be done
1969   // be aware that all values are relative values (scale: timebin-width; pad-width) and are integer values on special scale
1970   
1971   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1972   // !!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   !!
1973   // !!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.    !!
1974   // !!This has to be taken into account by the GTU. Furthermore a Lorentz correction might have to be applied to the offset (see below).                                                                                             !!
1975   // !!In this implementation it is assumed that no miscalibration containing changing drift velocities in the amplification region is used.                                                                                          !!
1976   // !!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    !!
1977   // !!valid (in approximation) if tFS is close to the beginning of the drift region.                                                                                                                                                 !!
1978   // !!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                                                                    !!
1979   // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1980
1981   // which formats should be chosen?
1982   AliTRDtrapAlu denomAlu;
1983   denomAlu.Init(20,8);       
1984   AliTRDtrapAlu numAlu;
1985   numAlu.Init(20,8);     
1986   // is this enough pre-comma place? covers the range of the 13 bit-word of the transmitted offset
1987   // offset measured in coord. of left channel must be between -0.5 and 1.5; 14 pre-comma bits because numerator can be big
1988
1989   for (Int_t i = 0; i<4; i++) {
1990     if (mADC[i] == -1) continue;
1991       
1992     Int_t num0  = (mN[i]*mXX[i]-mX[i]*mX[i]);
1993     if (num0 < 0) {
1994       denomAlu.AssignInt(-num0);    // num0 does not have to be interpreted as having past-comma bits -> AssignInt
1995       denomAlu.SetSign(-1);
1996     }
1997     else {
1998       denomAlu.AssignInt(num0);
1999       denomAlu.SetSign(1);
2000     }
2001     
2002     Int_t num1  = mN[i]*mXY[i] - mX[i]*mY[i];
2003     if (num1 < 0) {
2004       numAlu.AssignFormatted(-num1); // value of num1 is already formatted to have 8 past-comma bits
2005       numAlu.SetSign(-1);
2006     }
2007     else {
2008       numAlu.AssignFormatted(num1);
2009       numAlu.SetSign(1);
2010     }
2011     numAlu    = numAlu/denomAlu;
2012     mSlope[i]   = numAlu.GetSignedValue();
2013    
2014     Int_t num2  = mXX[i]*mY[i] - mX[i]*mXY[i];
2015    
2016     if (num2 < 0) {
2017       numAlu.AssignFormatted(-num2);
2018       numAlu.SetSign(-1);
2019     }
2020     else {
2021       numAlu.AssignFormatted(num2);
2022       numAlu.SetSign(1);
2023     }
2024    
2025     numAlu    = numAlu/denomAlu;
2026    
2027     
2028     mOffset[i]  = numAlu.GetSignedValue();
2029     numAlu.SetSign(1);
2030     denomAlu.SetSign(1);
2031        
2032                                  
2033     //numAlu.AssignInt(mADC[i]+1);   // according to TRAP-manual but trafo not to middle of chamber (0.5 channels away)             
2034     numAlu.AssignDouble((Double_t)mADC[i] + 1.5);      // numAlu has enough pre-comma place for that; correct trafo, best values
2035     mOffset[i]  = mOffset[i] + numAlu.GetValue();      // transform offset to a coord.system relative to chip; +1 to avoid neg. values 
2036     
2037     // up to here: adc-mapping according to TRAP and in line with pad-col mapping
2038     // reverse adc-counting to be again in line with the online mapping
2039     mADC[i]     = fNADC - 4 - mADC[i];                 // fNADC-4-mADC[i]: 0..17; remapping necessary;
2040     mADC[i]     = mADC[i] + 2; 
2041     // +2: mapping onto original ADC-online-counting: inner adc's corresponding to a chip's pasa: number 2..19
2042   }
2043
2044   // adc-counting is corresponding to online mapping; use AliTRDfeeParam::GetPadColFromADC to get the pad to which adc is connected; 
2045   // pad-column mapping is reverse to adc-online mapping; TRAP adc-mapping is in line with pad-mapping (increase in same direction);
2046   
2047   // transform parameters to the local coordinate-system of a stack (used by GTU)
2048   AliTRDpadPlane* padPlane = fGeo->CreatePadPlane(fLayer,fStack);
2049   
2050   Double_t padWidthI = padPlane->GetWidthIPad()*10.0; // get values in cm; want them in mm
2051   //Double_t padWidthO = padPlane->GetWidthOPad()*10; // difference between outer pad-widths not included; in real TRAP??
2052   
2053   // difference between width of inner and outer pads of a row is not accounted for;
2054   
2055   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 
2056   Double_t eCharge  = 0.3;                            // unit charge in (GeV/c)/m*T
2057   Double_t ptMin   = 2.3;                            // minimum transverse momentum (GeV/c); to be adjusted(?)
2058   
2059   Double_t granularityOffset = 0.160;                // granularity for offset in mm
2060   Double_t granularitySlope  = 0.140;                // granularity for slope  in mm     
2061     
2062   // get the coordinates in SM-system; parameters: 
2063   
2064   Double_t zPos       =  (padPlane->GetRowPos(fRow))*10.0;  // z-position of the MCM; fRow is counted on a chamber; SM consists of 5 
2065   // zPos is position of pad-borders;
2066   Double_t zOffset = 0.0;
2067   if ( fRow == 0 || fRow == 15 ) {
2068       zOffset = padPlane->GetLengthOPad();
2069   }
2070   else {
2071       zOffset = padPlane->GetLengthIPad();
2072   }
2073   zOffset = (-1.0) * zOffset/2.0;
2074   // turn zPos to be z-coordinate at middle of pad-row
2075   zPos = zPos + zOffset;
2076
2077       
2078   Double_t xPos       =  0.0;                               // x-position of the upper border of the drift-chamber of actual layer
2079   Int_t    icol       =  0;                                 // column-number of adc-channel
2080   Double_t yPos[4];                                         // y-position of the pad to which ADC is connected
2081   Double_t dx         = 30.0;                               // height of drift-chamber in mm; maybe retrieve from AliTRDGeometry
2082   Double_t freqSample = fFeeParam->GetSamplingFrequency();  // retrieve the sampling frequency (10.019750 MHz)
2083   Double_t vdrift     = fCal->GetVdriftAverage(fChaId);     // averaged drift velocity for this detector (1.500000 cm/us)
2084   Int_t    nrOfDriftTimeBins = Int_t(dx/10.0*freqSample/vdrift); // the number of time bins in the drift region (20)
2085   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
2086   //Double_t lorAngle   =  7.0;                             // Lorentz-angle in degrees
2087   Double_t tiltAngle  = padPlane->GetTiltingAngle();        // sign-respecting tilting angle of pads in actual layer
2088   Double_t tiltTan    = TMath::Tan(TMath::Pi()/180.0 * tiltAngle);
2089   //Double_t lorTan     = TMath::Tan(TMath::Pi()/180.0 * lorAngle);
2090
2091   Double_t alphaMax[4];                            // maximum deflection from the direction to the primary vertex; granularity of hit pads
2092   Double_t slopeMin[4];                            // local limits for the deflection
2093   Double_t slopeMax[4];
2094   Int_t   mslopeMin[4];                            // in granularity units; to be compared to mSlope[i]
2095   Int_t   mslopeMax[4];
2096
2097
2098   //x coord. of upper side of drift chambers in local SM-system (in mm)
2099   //obtained by evaluating the x-range of the hits; should be crosschecked; only drift, not amplification region taken into account (30mm);
2100   //the y-deflection is given as difference of y between lower and upper side of drift-chamber, not pad-plane;
2101   switch(fLayer) 
2102     {
2103     case 0: 
2104       xPos = 3003.0;
2105       break;
2106     case 1:
2107       xPos = 3129.0;
2108       break;
2109     case 2:
2110       xPos = 3255.0;
2111       break;
2112     case 3:
2113       xPos = 3381.0;
2114       break;
2115     case 4:
2116       xPos = 3507.0;
2117       break;
2118     case 5:
2119       xPos = 3633.0;
2120       break;
2121     }
2122  
2123   // calculation of offset-correction n: 
2124
2125   Int_t nCorrectOffset = (fRobPos % 2 == 0) ? ((fMcmPos % 4)) : ( 4 + (fMcmPos % 4));  
2126  
2127   nCorrectOffset = (nCorrectOffset - 4)*18 - 1;
2128   if (nCorrectOffset < 0) {
2129     numAlu.AssignInt(-nCorrectOffset);
2130     numAlu.SetSign(-1);
2131   }
2132   else {
2133     numAlu.AssignInt(nCorrectOffset);
2134     numAlu.SetSign(1);
2135   }
2136   nCorrectOffset = numAlu.GetSignedValue();         
2137   Double_t mCorrectOffset = padWidthI/granularityOffset; // >= 0.0
2138  
2139   // calculation of slope-correction
2140
2141   // this is only true for tracks coming (approx.) from primary vertex
2142   // everything is evaluated for a tracklet covering the whole drift chamber
2143   Double_t cCorrectSlope = (-lorTan*dx + zPos/xPos*dx*tiltTan)/granularitySlope;
2144   // Double_t cCorrectSlope =  zPos/xPos*dx*tiltTan/granularitySlope;
2145   // zPos can be negative! for track from primary vertex: zOut-zIn > 0 <=> zPos > 0
2146   
2147   if (cCorrectSlope < 0) {
2148       numAlu.AssignDouble(-cCorrectSlope);
2149       numAlu.SetSign(-1);
2150   }
2151   else {
2152       numAlu.AssignDouble(cCorrectSlope);
2153       numAlu.SetSign(1);
2154   }
2155   cCorrectSlope = numAlu.GetSignedValue();
2156  
2157   // convert slope to deflection between upper and lower drift-chamber position (slope is given in pad-unit/time-bins)
2158   // different pad-width of outer pads of a pad-plane not taken into account
2159   // 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)
2160   // 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
2161   
2162   
2163   Double_t mCorrectSlope = (Double_t)(nrOfDriftTimeBins)*padWidthI/granularitySlope;  // >= 0.0
2164
2165   AliTRDtrapAlu correctAlu;
2166   correctAlu.Init(20,8);
2167   
2168   AliTRDtrapAlu offsetAlu;
2169   offsetAlu.Init(13,0,-0x1000,0x0FFF);          // 13 bit-word; 2-complement (1 sign-bit); asymmetric range
2170   
2171   AliTRDtrapAlu slopeAlu;
2172   slopeAlu.Init(7,0,-0x40,0x3F);                // 7 bit-word;  2-complement (1 sign-bit);
2173
2174   for (Int_t i = 0; i<4; i++) {
2175     
2176     if (mADC[i] == -1) continue;
2177     
2178     icol = fFeeParam->GetPadColFromADC(fRobPos,fMcmPos,mADC[i]); // be aware that mADC[i] contains the ADC-number according to online-mapping
2179     yPos[i]   = (padPlane->GetColPos(icol))*10.0;
2180     
2181     
2182     // offset:
2183     
2184     correctAlu.AssignDouble(mCorrectOffset);     // done because max. accuracy is 8 bit
2185     mCorrectOffset = correctAlu.GetValueWhole(); // cut offset correction to 8 past-comma bit accuracy
2186     mOffset[i]  = (Int_t)((mCorrectOffset)*(Double_t)(mOffset[i] + nCorrectOffset)); 
2187     mOffset[i]  = mOffset[i]*(-1);                   // adjust to direction of y-axes in online simulation
2188     
2189     if (mOffset[i] < 0) {
2190       numAlu.AssignFormatted(-mOffset[i]);
2191       numAlu.SetSign(-1);
2192     }
2193     else {
2194       numAlu.AssignFormatted(mOffset[i]);
2195       numAlu.SetSign(1);
2196     }
2197
2198     offsetAlu = numAlu; 
2199     mOffset[i]  = offsetAlu.GetSignedValue();  
2200
2201     
2202     // slope:
2203     
2204     correctAlu.AssignDouble(mCorrectSlope);
2205     mCorrectSlope = correctAlu.GetValueWhole();
2206     
2207     mSlope[i]   = (Int_t)((mCorrectSlope*(Double_t)mSlope[i]) + cCorrectSlope);
2208
2209     if (mSlope[i] < 0) {
2210       numAlu.AssignFormatted(-mSlope[i]);
2211       numAlu.SetSign(-1);
2212     }
2213     else {
2214       numAlu.AssignFormatted(mSlope[i]);
2215       numAlu.SetSign(1);
2216     }
2217
2218     slopeAlu  = numAlu;     // here all past-comma values are cut, not rounded; alternatively add +0.5 before cutting (means rounding)
2219     mSlope[i]   = slopeAlu.GetSignedValue(); 
2220        
2221     // local (LTU) limits for the deflection 
2222     // ATan returns angles in radian
2223     alphaMax[i]  = TMath::ASin(eCharge*magField/(2.0*ptMin)*(TMath::Sqrt(xPos*xPos + yPos[i]*yPos[i]))/1000.0); // /1000: mm->m
2224     slopeMin[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) - alphaMax[i]))/granularitySlope;
2225     slopeMax[i]  = dx*(TMath::Tan(TMath::ATan(yPos[i]/xPos) + alphaMax[i]))/granularitySlope;
2226     
2227     if (slopeMin[i] < 0) {
2228       slopeAlu.AssignDouble(-slopeMin[i]);
2229       slopeAlu.SetSign(-1);
2230     }
2231     else { 
2232       slopeAlu.AssignDouble(slopeMin[i]);
2233       slopeAlu.SetSign(1);
2234     }
2235     mslopeMin[i] = slopeAlu.GetSignedValue();  // the borders should lie inside the range of mSlope -> usage of slopeAlu again
2236    
2237     if (slopeMax[i] < 0) {
2238       slopeAlu.AssignDouble(-slopeMax[i]);
2239       slopeAlu.SetSign(-1);
2240     }
2241     else {
2242       slopeAlu.AssignDouble(slopeMax[i]);
2243       slopeAlu.SetSign(1);
2244     }
2245     mslopeMax[i] = slopeAlu.GetSignedValue();
2246   }
2247
2248   // suppress submission of tracks with low stiffness
2249   // put parameters in 32bit-word and submit (write to file as root-file; sort after SM, stack, layer, chamber) 
2250
2251   // sort tracklet-words in ascending y-order according to the offset (according to mADC would also be possible)
2252   // up to now they are sorted according to maximum hit sum
2253   // is the sorting really done in the TRAP-chip?
2254   
2255   Int_t order[4] = {-1,-1,-1,-1};
2256   Int_t wordnr = 0;   // number of tracklet-words
2257   
2258   for(Int_t j = 0; j < fMaxTracklets; j++) {
2259       if( mADC[j] == -1) continue; 
2260       //if( (mADC[j] == -1) || (mSlope[j] < mslopeMin[j]) || (mSlope[j] > mslopeMax[j])) continue; // this applies a pt-cut
2261       wordnr++;
2262       if( wordnr-1 == 0) {
2263           order[0] = j;
2264           continue;
2265       }
2266       // wordnr-1>0, wordnr-1<4
2267       order[wordnr-1] = j;
2268       for( Int_t k = 0; k < wordnr-1; k++) {
2269           if( mOffset[j] < mOffset[order[k]] ) {
2270               for( Int_t l = wordnr-1; l > k; l-- ) {
2271                   order[l] = order[l-1];
2272               }
2273               order[k] = j;
2274               break;
2275           }
2276           
2277       }
2278   }
2279         
2280   // fill the bit-words in ascending order and without gaps
2281   UInt_t bitWord[4] = {0,0,0,0};                 // attention: unsigned int to have real 32 bits (no 2-complement)
2282   for(Int_t j = 0; j < wordnr; j++) { // only "wordnr" tracklet-words
2283       //Bool_t rem1 = kTRUE;
2284     
2285     Int_t i = order[j];
2286     //bit-word is 2-complement and therefore without sign
2287     bitWord[j] =   1; // this is the starting 1 of the bit-word (at 33rd position); the 1 must be ignored
2288     //printf("\n");
2289     UInt_t shift  = 0;
2290     UInt_t shift2 = 0;
2291         
2292         
2293
2294
2295     /*printf("mean charge: %d\n",mMeanCharge[i]);
2296     printf("row: %d\n",fRow);
2297     printf("slope: %d\n",mSlope[i]);
2298     printf("pad position: %d\n",mOffset[i]);
2299     printf("channel: %d\n",mADC[i]);*/
2300
2301     // electron probability (currently not implemented; the mean charge is just scaled)
2302     shift = (UInt_t)mMeanCharge[i];
2303     for(Int_t iBit = 0; iBit < 8; iBit++) {
2304       bitWord[j]  = bitWord[j]<<1;
2305       bitWord[j] |= (shift>>(7-iBit))&1;               
2306       //printf("0");
2307     }
2308
2309     // pad row
2310     shift = (UInt_t)fRow;
2311     for(Int_t iBit = 0; iBit < 4; iBit++) {
2312       bitWord[j]  = bitWord[j]<<1;
2313       bitWord[j] |= (shift>>(3-iBit))&1;
2314       //printf("%d", (fRow>>(3-iBit))&1);
2315     }
2316     
2317     // deflection length
2318     if(mSlope[i] < 0) {
2319         shift = (UInt_t)(-mSlope[i]);
2320         // shift2 is 2-complement of shift
2321         shift2 = 1;
2322         for(Int_t iBit = 1; iBit < 7; iBit++) {
2323             shift2  = shift2<<1;
2324             shift2 |= (1-((shift)>>(6-iBit))&1);
2325             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2326         }
2327         shift2 = shift2 + 1;
2328         //printf("1");
2329         for(Int_t iBit = 0; iBit < 7; iBit++) {
2330             bitWord[j]  = bitWord[j]<<1;
2331             bitWord[j] |= (shift2>>(6-iBit))&1;
2332             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2333         }
2334     }
2335     else {
2336         shift = (UInt_t)(mSlope[i]);
2337         bitWord[j]  = bitWord[j]<<1;
2338         bitWord[j]   |= 0;
2339         //printf("0");
2340         for(Int_t iBit = 1; iBit < 7; iBit++) {
2341             bitWord[j]  = bitWord[j]<<1;
2342             bitWord[j] |= (shift>>(6-iBit))&1;
2343             //printf("%d",(mSlope[i]>>(6-iBit))&1);
2344         }
2345     }
2346
2347     // pad position
2348     if(mOffset[i] < 0) {
2349         shift = (UInt_t)(-mOffset[i]);
2350         shift2 = 1;
2351         for(Int_t iBit = 1; iBit < 13; iBit++) {
2352             shift2  = shift2<<1;
2353             shift2 |= (1-((shift)>>(12-iBit))&1);
2354             //printf("%d",(1-((-mOffset[i])>>(12-iBit))&1));
2355         }
2356         shift2 = shift2 + 1;
2357         //printf("1");
2358         for(Int_t iBit = 0; iBit < 13; iBit++) {
2359             bitWord[j]  = bitWord[j]<<1;
2360             bitWord[j] |= (shift2>>(12-iBit))&1;
2361             //printf("%d",(1-((-mSlope[i])>>(6-iBit))&1));
2362         }
2363     }
2364     else {
2365         shift = (UInt_t)mOffset[i];
2366         bitWord[j] = bitWord[j]<<1;
2367         bitWord[j]   |= 0; 
2368         //printf("0");
2369         for(Int_t iBit = 1; iBit < 13; iBit++) {
2370             bitWord[j]  = bitWord[j]<<1;
2371             bitWord[j] |= (shift>>(12-iBit))&1;
2372             //printf("%d",(mOffset[i]>>(12-iBit))&1);
2373         }
2374     }
2375
2376
2377         
2378     //printf("bitWord: %u\n",bitWord[j]);
2379     //printf("adc: %d\n",mADC[i]);
2380     fMCMT[j] = bitWord[j];
2381   }
2382     
2383   //printf("\n");
2384
2385   
2386   delete [] qsum;
2387   delete [] ieffped;
2388
2389   delete [] X;
2390   delete [] XX;
2391   delete [] Y;
2392   delete [] YY;
2393   delete [] XY;
2394   delete [] N;
2395   delete [] QT0;
2396   delete [] QT1;
2397
2398   delete [] hitSum;
2399   delete [] selectPair;
2400
2401   delete padPlane;
2402
2403
2404
2405  
2406
2407   // output-part; creates an independent output .root folder if uncommented;
2408   
2409   // structure: in the current directory a root-file called "TRD_readout_tree.root" is stored with subdirectories SMxx/sx (supermodule, stack);
2410   // in each of these subdirectories 6 trees according to layers are saved, called lx; 
2411   // 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;
2412   // branch-name: mcmxxxwd; 
2413   // another branch contains the channel-number (mcmxxxch)
2414   
2415   /*
2416   AliLog::SetClassDebugLevel("AliTRDmcmSim", 10);
2417   AliLog::SetFileOutput("../log/tracklet.log");
2418   
2419  
2420   UInt_t* trackletWord;
2421   Int_t*  adcChannel;
2422
2423   Int_t u = 0;
2424     
2425   // testing for wordnr in order to speed up the simulation
2426   
2427   if (wordnr == 0 && fNextEvent == 0) {
2428     return;
2429   }
2430
2431    
2432   Int_t mcmNr = fRobPos * (fGeo->MCMmax()) + fMcmPos;
2433   
2434   Char_t* SMName    = new Char_t[4];
2435   Char_t* stackName = new Char_t[2];
2436   Char_t* layerName = new Char_t[2];
2437
2438   Char_t* treeName  = new Char_t[2];
2439   Char_t* treeTitle = new Char_t[8];
2440   Char_t* branchNameWd = new Char_t[8]; // mcmxxxwd bit word
2441   Char_t* branchNameCh = new Char_t[8]; // mcmxxxch channel
2442    
2443   Char_t* dirName = NULL;
2444   Char_t* treeFile  = NULL; 
2445   Char_t* evFile = NULL;
2446   Char_t* curDir = new Char_t[255];
2447   
2448   for (Int_t i = 0; i<255; i++) {
2449     curDir[i]='n';
2450   }
2451   sprintf(curDir,"%s",gSystem->BaseName(gSystem->WorkingDirectory()));
2452   Int_t rawEvent = 0;
2453   Int_t nrPos = 3;
2454   
2455
2456   while(curDir[nrPos]!='n'){
2457     
2458    
2459     switch(curDir[nrPos]) {
2460     case '0':
2461       rawEvent = rawEvent*10 + 0;
2462       break;
2463     case '1':
2464       rawEvent = rawEvent*10 + 1;
2465       break;
2466     case '2':
2467       rawEvent = rawEvent*10 + 2;
2468       break;
2469     case '3':
2470       rawEvent = rawEvent*10 + 3;
2471       break;
2472     case '4':
2473       rawEvent = rawEvent*10 + 4;
2474       break;
2475     case '5':
2476       rawEvent = rawEvent*10 + 5;
2477       break;
2478     case '6':
2479       rawEvent = rawEvent*10 + 6;
2480       break;
2481     case '7':
2482       rawEvent = rawEvent*10 + 7;
2483       break;
2484     case '8':
2485       rawEvent = rawEvent*10 + 8;
2486       break;
2487     case '9':
2488       rawEvent = rawEvent*10 + 9;
2489       break;
2490    
2491     }
2492     nrPos = nrPos + 1; 
2493   }
2494   delete [] curDir;
2495     
2496   //if (!gSystem->ChangeDirectory("../TRD_Tracklet")) {
2497    // gSystem->MakeDirectory("../TRD_Tracklet");
2498     //gSystem->ChangeDirectory("../TRD_Tracklet");
2499     //}
2500
2501   gSystem->ChangeDirectory("..");
2502   
2503   
2504   TFile *f = new TFile("TRD_readout_tree.root","update");
2505   TTree *tree          = NULL;
2506   TBranch *branch      = NULL;
2507   TBranch *branchCh   = NULL; 
2508    
2509   Int_t iEventNr = 0; 
2510   Int_t dignr = 10; // nr of digits of a integer
2511   Int_t space = 1;  // additional char-space
2512   
2513   
2514   evFile = new Char_t[2+space];
2515   sprintf(evFile,"ev%d",iEventNr);
2516
2517   
2518   while(f->cd(evFile)){
2519     iEventNr = iEventNr + 1;
2520     if (iEventNr/dignr > 0) {
2521       dignr = dignr * 10;
2522       space = space + 1;
2523     }
2524     delete [] evFile;
2525     evFile = NULL;
2526     evFile = new Char_t[2+space];
2527     sprintf(evFile,"ev%d",iEventNr); 
2528   }
2529   
2530   if(iEventNr == rawEvent) { fNextEvent = 1; } // new event
2531    
2532   if (fNextEvent == 1) {
2533     fNextEvent = 0;
2534     // turn to head directory
2535     f->mkdir(evFile);
2536     f->cd(evFile);
2537     // create all subdirectories and trees in case of new event
2538    
2539      
2540     for (Int_t iSector = 0; iSector < 18; iSector++) {
2541   
2542       if (iSector < 10) {
2543         sprintf(SMName,"SM0%d",iSector);
2544       }
2545       else {
2546         sprintf(SMName,"SM%d",iSector);
2547       }
2548
2549       for (Int_t iStack = 0; iStack < 5; iStack++) {
2550         sprintf(stackName,"s%d",iStack);
2551         
2552         f->cd(evFile);
2553         if (iStack == 0) {
2554           gDirectory->mkdir(SMName);
2555         }
2556         gDirectory->cd(SMName);
2557         gDirectory->mkdir(stackName);
2558         gDirectory->cd(stackName);
2559         
2560         for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
2561           sprintf(layerName,"l%d",iLayer);
2562           sprintf(treeName,"%s",layerName);
2563           sprintf(treeTitle,"%s%s%s",SMName,stackName,layerName);
2564           tree = new TTree(treeName,treeTitle);
2565           tree->Write("",TObject::kOverwrite);
2566           delete tree;
2567           tree = NULL;
2568         }
2569       }
2570     }
2571    
2572     
2573   }
2574
2575     
2576   else {
2577     iEventNr = iEventNr - 1;
2578     dignr = dignr/10;
2579     if (iEventNr/dignr == 0) space = space - 1;
2580     delete [] evFile;
2581     evFile = NULL;
2582     evFile = new Char_t[2+space];
2583     sprintf(evFile,"ev%d",iEventNr);    
2584   }
2585   
2586   if (wordnr == 0) {
2587     delete [] SMName;
2588     delete [] stackName;
2589     delete [] layerName;
2590     delete [] treeName;
2591     delete [] treeTitle;
2592     delete [] branchNameWd;
2593     delete [] branchNameCh;
2594     delete [] evFile;
2595     f->Close();
2596     dirName   = new Char_t[6+space];
2597     //sprintf(dirName,"../raw%d",iEventNr);
2598     sprintf(dirName,"raw%d",iEventNr);
2599     gSystem->ChangeDirectory(dirName);
2600     delete [] dirName;
2601     return;
2602   }
2603   
2604   dirName   = new Char_t[6+space];
2605   //sprintf(dirName,"../raw%d",iEventNr);
2606   sprintf(dirName,"raw%d",iEventNr);
2607  
2608   f->cd(evFile);
2609  
2610   if (fSector < 10) {
2611     sprintf(SMName,"SM0%d",fSector);
2612   }
2613   else {
2614     sprintf(SMName,"SM%d",fSector);
2615   }
2616   sprintf(stackName,"s%d",fStack);
2617   sprintf(layerName,"l%d",fLayer);
2618   sprintf(treeName,"%s",layerName);
2619   sprintf(treeTitle,"%s%s%s",SMName,stackName,layerName);
2620  
2621   treeFile = new Char_t[13+space];
2622   sprintf(treeFile,"%s/%s/%s/%s",evFile,SMName,stackName,treeName);
2623   delete [] evFile;
2624   evFile = NULL;
2625   gDirectory->cd(SMName);
2626   gDirectory->cd(stackName);
2627   tree = (TTree*)f->Get(treeFile);
2628   delete [] treeFile;
2629   treeFile = NULL;
2630
2631
2632   //make branch with number of words and fill
2633
2634   if (mcmNr < 10) {
2635     sprintf(branchNameWd,"mcm00%dwd",mcmNr);
2636     sprintf(branchNameCh,"mcm00%dch",mcmNr);
2637   }
2638   else {
2639     if (mcmNr < 100) {
2640       sprintf(branchNameWd,"mcm0%dwd",mcmNr); 
2641       sprintf(branchNameCh,"mcm0%dch",mcmNr);
2642     }
2643     else {
2644       sprintf(branchNameWd,"mcm%dwd",mcmNr); 
2645       sprintf(branchNameCh,"mcm%dch",mcmNr);
2646     }
2647   }
2648
2649       
2650   
2651   // fill the tracklet word; here wordnr > 0
2652
2653   trackletWord = new UInt_t[fMaxTracklets];
2654   adcChannel   = new Int_t[fMaxTracklets];
2655
2656   for (Int_t j = 0; j < fMaxTracklets; j++) {
2657       Int_t i = order[j];
2658       trackletWord[j] = 0;
2659       adcChannel[j] = -1;
2660       if (bitWord[j]!=0) {
2661           trackletWord[u] = bitWord[j];
2662           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
2663           
2664           //fMCMT[u] = bitWord[j];
2665           u = u + 1;
2666       }
2667   }
2668   
2669   branch = tree->GetBranch(branchNameWd);
2670   if(!branch) {
2671     //make branch and fill
2672     branch = tree->Branch(branchNameWd,trackletWord,"trackletWord[4]/i"); // 32 bit unsigned integer
2673     branch->Fill();
2674   }
2675
2676   branchCh = tree->GetBranch(branchNameCh);
2677   if(!branchCh) {
2678     //make branch and fill
2679     branchCh = tree->Branch(branchNameCh,adcChannel,"adcChannel[4]/i"); // 32 bit unsigned integer
2680     branchCh->Fill();
2681   }
2682  
2683
2684   tree->Write("",TObject::kOverwrite);
2685  
2686
2687   delete [] SMName;
2688   delete [] stackName;
2689   delete [] layerName;
2690   delete [] treeName;
2691   delete [] treeTitle;
2692   delete [] branchNameWd;
2693   delete [] branchNameCh;
2694   delete [] trackletWord;
2695   delete [] adcChannel;
2696   
2697   f->Close();
2698   gSystem->ChangeDirectory(dirName);
2699   delete [] dirName;
2700   */
2701
2702
2703   // to be done:
2704   // error measure for quality of fit (not necessarily needed for the trigger)
2705   // cluster quality threshold (not yet set)
2706   // electron probability
2707   
2708
2709    
2710 }
2711
2712
2713
2714
2715
2716