Cleaup the situation NTimeBin in simulation
[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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //  TRD MCM (Multi Chip Module) simulator                                    //
21 //  which simulated the TRAP processing after the AD-conversion              //
22 //  The relevant parameters (i.e. configuration registers of the TRAP        //
23 //  configuration are taken from AliTRDtrapConfig.                           //
24 //                                                                           //
25 ///////////////////////////////////////////////////////////////////////////////
26
27 #include <fstream>  // needed for raw data dump
28
29 #include <TCanvas.h>
30 #include <TH1F.h>
31 #include <TH2F.h>
32 #include <TGraph.h>
33 #include <TLine.h>
34 #include <TMath.h>
35 #include <TRandom.h>
36 #include <TClonesArray.h>
37
38 #include "AliLog.h"
39 #include "AliRun.h"
40 #include "AliRunLoader.h"
41 #include "AliLoader.h"
42 #include "AliTRDdigit.h"
43
44 #include "AliTRDfeeParam.h"
45 #include "AliTRDtrapConfig.h"
46 #include "AliTRDSimParam.h"
47 #include "AliTRDgeometry.h"
48 #include "AliTRDcalibDB.h"
49 #include "AliTRDdigitsManager.h"
50 #include "AliTRDarrayADC.h"
51 #include "AliTRDarrayDictionary.h"
52 #include "AliTRDpadPlane.h"
53 #include "AliTRDtrackletMCM.h"
54 #include "AliTRDmcmSim.h"
55
56 #include "AliMagF.h"
57 #include "TGeoGlobalMagField.h"
58
59 ClassImp(AliTRDmcmSim)
60
61 Bool_t AliTRDmcmSim::fgApplyCut = kTRUE;
62
63 Float_t AliTRDmcmSim::fgChargeNorm = 65000.;
64 Int_t AliTRDmcmSim::fgAddBaseline = 0;
65
66 Int_t AliTRDmcmSim::fgPidNBinsQ0 = 40;
67 Int_t AliTRDmcmSim::fgPidNBinsQ1 = 50;
68 Bool_t AliTRDmcmSim::fgPidLutDelete = kFALSE;
69 Int_t AliTRDmcmSim::fgPidLutDefault[40][50] = {
70   { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0  }, 
71   { 0, 9, 6, 12, 29, 53, 76, 94, 107, 116, 122, 126, 128, 129, 129, 129, 128, 127, 126, 124, 122, 120, 117, 115, 112, 109, 107, 104, 101, 99, 96, 94, 91, 89, 87, 85, 83, 81, 79, 78, 77, 75, 74, 73, 72, 72, 71, 71, 70, 70  }, 
72   { 0, 14, 8, 17, 37, 66, 94, 116, 131, 140, 146, 150, 152, 153, 153, 152, 150, 148, 145, 143, 139, 136, 132, 129, 125, 121, 118, 114, 110, 107, 104, 101, 98, 95, 93, 91, 89, 87, 85, 83, 82, 81, 80, 79, 78, 77, 77, 76, 76, 75  }, 
73   { 0, 33, 19, 34, 69, 112, 145, 167, 181, 189, 194, 196, 197, 197, 196, 194, 191, 188, 184, 180, 175, 170, 165, 159, 154, 148, 143, 137, 132, 127, 123, 118, 114, 111, 107, 104, 101, 99, 96, 94, 92, 91, 89, 88, 87, 86, 85, 85, 84, 84  }, 
74   { 0, 82, 52, 83, 136, 180, 205, 218, 226, 230, 232, 233, 233, 233, 232, 230, 228, 226, 223, 219, 215, 210, 205, 199, 193, 187, 180, 173, 167, 160, 154, 148, 142, 136, 131, 127, 122, 119, 115, 112, 109, 106, 104, 102, 100, 99, 97, 96, 95, 94  }, 
75   { 0, 132, 96, 136, 185, 216, 231, 238, 242, 244, 245, 245, 245, 245, 245, 244, 243, 242, 240, 238, 236, 233, 230, 226, 222, 217, 212, 206, 200, 193, 187, 180, 173, 167, 161, 155, 149, 144, 139, 134, 130, 126, 123, 120, 117, 114, 112, 110, 108, 107  }, 
76   { 0, 153, 120, 160, 203, 227, 238, 243, 246, 247, 248, 249, 249, 249, 248, 248, 247, 246, 245, 244, 243, 241, 239, 237, 234, 231, 228, 224, 219, 215, 209, 204, 198, 192, 186, 180, 174, 168, 163, 157, 152, 147, 143, 139, 135, 131, 128, 125, 123, 120  }, 
77   { 0, 156, 128, 166, 207, 229, 239, 244, 247, 248, 249, 249, 249, 249, 249, 249, 248, 247, 247, 246, 244, 243, 242, 240, 238, 236, 233, 230, 227, 224, 220, 216, 212, 207, 202, 197, 192, 187, 181, 176, 171, 166, 161, 156, 152, 148, 144, 140, 137, 134  }, 
78   { 0, 152, 128, 166, 206, 228, 239, 244, 246, 248, 249, 249, 249, 249, 249, 248, 248, 247, 246, 245, 244, 243, 241, 240, 238, 236, 234, 232, 229, 226, 224, 220, 217, 214, 210, 206, 202, 197, 193, 188, 184, 179, 174, 170, 166, 161, 157, 153, 150, 146  }, 
79   { 0, 146, 126, 164, 203, 226, 237, 243, 246, 247, 248, 248, 248, 248, 248, 247, 247, 246, 245, 244, 242, 241, 239, 238, 236, 234, 232, 230, 227, 225, 223, 220, 217, 215, 212, 209, 205, 202, 199, 195, 191, 187, 183, 179, 175, 171, 168, 164, 160, 156  }, 
80   { 0, 140, 123, 160, 200, 224, 235, 241, 244, 246, 247, 247, 247, 247, 247, 246, 245, 244, 243, 242, 240, 238, 237, 235, 233, 230, 228, 226, 224, 221, 219, 217, 215, 212, 210, 207, 205, 202, 200, 197, 194, 191, 188, 184, 181, 178, 174, 171, 168, 164  }, 
81   { 0, 133, 119, 156, 196, 220, 233, 239, 243, 245, 245, 246, 246, 246, 245, 244, 243, 242, 241, 239, 237, 235, 233, 231, 229, 226, 224, 221, 219, 216, 214, 212, 210, 208, 206, 204, 202, 199, 197, 195, 193, 191, 188, 186, 183, 181, 178, 175, 172, 169  }, 
82   { 0, 127, 115, 152, 192, 217, 230, 237, 241, 243, 244, 244, 244, 244, 243, 242, 241, 240, 238, 236, 234, 232, 229, 227, 224, 221, 218, 216, 213, 210, 208, 206, 203, 201, 200, 198, 196, 194, 193, 191, 190, 188, 186, 185, 183, 181, 179, 177, 174, 172  }, 
83   { 0, 121, 111, 147, 187, 213, 227, 235, 239, 241, 242, 243, 243, 242, 241, 240, 239, 237, 236, 233, 231, 228, 225, 222, 219, 216, 213, 210, 207, 204, 201, 199, 196, 194, 192, 191, 189, 188, 187, 185, 184, 183, 182, 181, 180, 178, 177, 176, 174, 172  }, 
84   { 0, 116, 107, 142, 181, 209, 224, 232, 237, 239, 240, 241, 241, 240, 239, 238, 237, 235, 233, 230, 227, 224, 221, 218, 214, 211, 207, 204, 200, 197, 194, 191, 189, 187, 185, 183, 182, 180, 179, 178, 178, 177, 176, 175, 175, 174, 173, 172, 172, 170  }, 
85   { 0, 112, 103, 136, 176, 204, 220, 229, 234, 237, 238, 239, 239, 238, 237, 236, 234, 232, 230, 227, 224, 221, 217, 213, 209, 205, 201, 198, 194, 190, 187, 184, 181, 179, 177, 175, 174, 172, 171, 171, 170, 169, 169, 169, 168, 168, 168, 168, 167, 167  }, 
86   { 0, 107, 99, 131, 170, 199, 216, 226, 231, 234, 236, 237, 237, 236, 235, 234, 232, 230, 227, 224, 221, 217, 213, 209, 205, 200, 196, 192, 188, 184, 180, 177, 174, 172, 169, 167, 166, 164, 163, 162, 162, 161, 161, 161, 161, 161, 161, 162, 162, 162  }, 
87   { 0, 104, 94, 125, 164, 193, 212, 222, 228, 232, 233, 234, 234, 234, 233, 231, 229, 227, 224, 221, 218, 214, 210, 205, 201, 196, 191, 187, 182, 178, 174, 171, 168, 165, 162, 160, 158, 157, 155, 154, 154, 153, 153, 153, 153, 154, 154, 154, 155, 155  }, 
88   { 0, 100, 90, 119, 157, 188, 207, 219, 225, 229, 231, 232, 232, 231, 230, 229, 227, 224, 222, 218, 215, 211, 206, 202, 197, 192, 187, 182, 178, 173, 169, 165, 162, 158, 156, 153, 151, 149, 148, 147, 146, 146, 145, 145, 145, 146, 146, 147, 148, 148  }, 
89   { 0, 97, 86, 113, 150, 182, 202, 215, 222, 226, 228, 229, 229, 229, 228, 226, 224, 222, 219, 216, 212, 208, 203, 199, 194, 188, 183, 178, 173, 169, 164, 160, 156, 153, 150, 147, 145, 143, 141, 140, 139, 138, 138, 138, 138, 138, 139, 139, 140, 141  }, 
90   { 0, 94, 82, 107, 144, 176, 197, 210, 218, 223, 225, 227, 227, 227, 226, 224, 222, 220, 217, 213, 209, 205, 201, 196, 191, 186, 180, 175, 170, 165, 160, 156, 152, 148, 145, 142, 139, 137, 135, 134, 132, 131, 131, 131, 131, 131, 131, 132, 132, 133  }, 
91   { 0, 92, 78, 102, 137, 169, 192, 206, 215, 220, 223, 224, 224, 224, 223, 222, 220, 217, 215, 211, 207, 203, 199, 194, 188, 183, 178, 172, 167, 162, 157, 152, 148, 144, 140, 137, 134, 132, 130, 128, 127, 125, 125, 124, 124, 124, 124, 125, 125, 126  }, 
92   { 0, 90, 75, 96, 131, 163, 187, 202, 211, 216, 220, 221, 222, 222, 221, 220, 218, 215, 212, 209, 205, 201, 197, 192, 187, 181, 176, 170, 165, 159, 154, 149, 145, 141, 137, 133, 130, 128, 125, 123, 122, 120, 119, 118, 118, 118, 118, 118, 119, 119  }, 
93   { 0, 88, 71, 91, 124, 157, 181, 197, 207, 213, 217, 219, 219, 219, 219, 217, 216, 213, 211, 207, 204, 200, 195, 190, 185, 180, 174, 169, 163, 158, 152, 147, 142, 138, 134, 130, 127, 124, 121, 119, 117, 116, 114, 114, 113, 112, 112, 112, 112, 113  }, 
94   { 0, 87, 68, 86, 118, 151, 176, 192, 203, 210, 214, 216, 217, 217, 217, 215, 214, 212, 209, 206, 202, 198, 194, 189, 184, 179, 173, 167, 162, 156, 151, 146, 141, 136, 132, 128, 124, 121, 118, 116, 114, 112, 110, 109, 108, 108, 107, 107, 107, 107  }, 
95   { 0, 85, 65, 81, 112, 144, 170, 188, 199, 206, 211, 213, 214, 215, 214, 213, 212, 210, 207, 204, 201, 197, 193, 188, 183, 178, 172, 167, 161, 155, 150, 145, 140, 135, 130, 126, 122, 119, 116, 113, 111, 109, 107, 106, 105, 104, 103, 103, 102, 102  }, 
96   { 0, 84, 62, 77, 106, 138, 165, 183, 195, 203, 208, 210, 212, 212, 212, 211, 210, 208, 206, 203, 200, 196, 192, 187, 183, 177, 172, 166, 161, 155, 150, 144, 139, 134, 129, 125, 121, 117, 114, 111, 109, 106, 104, 103, 101, 100, 99, 99, 98, 98  }, 
97   { 0, 84, 60, 73, 101, 133, 159, 178, 191, 199, 204, 208, 209, 210, 210, 209, 208, 206, 204, 202, 199, 195, 191, 187, 182, 177, 172, 166, 161, 155, 150, 144, 139, 134, 129, 124, 120, 116, 113, 110, 107, 104, 102, 100, 99, 98, 96, 96, 95, 95  }, 
98   { 0, 83, 58, 69, 96, 127, 154, 174, 187, 196, 201, 205, 207, 208, 208, 207, 206, 205, 203, 200, 197, 194, 190, 186, 182, 177, 172, 167, 161, 156, 150, 145, 139, 134, 129, 124, 120, 116, 112, 109, 106, 103, 101, 99, 97, 95, 94, 93, 92, 92  }, 
99   { 0, 82, 56, 66, 91, 121, 149, 169, 183, 192, 198, 202, 204, 206, 206, 206, 205, 203, 201, 199, 196, 193, 190, 186, 182, 177, 172, 167, 162, 156, 151, 145, 140, 135, 129, 125, 120, 116, 112, 108, 105, 102, 100, 97, 95, 94, 92, 91, 90, 89  }, 
100   { 0, 82, 54, 62, 86, 116, 144, 165, 179, 189, 195, 199, 202, 203, 204, 204, 203, 202, 200, 198, 196, 193, 189, 186, 182, 177, 173, 168, 163, 157, 152, 146, 141, 136, 130, 125, 121, 116, 112, 108, 105, 102, 99, 96, 94, 92, 91, 89, 88, 87  }, 
101   { 0, 82, 52, 59, 82, 111, 139, 160, 175, 185, 192, 197, 200, 201, 202, 202, 201, 200, 199, 197, 195, 192, 189, 186, 182, 178, 173, 168, 163, 158, 153, 148, 142, 137, 132, 127, 122, 117, 113, 109, 105, 102, 99, 96, 94, 92, 90, 88, 87, 85  }, 
102   { 0, 82, 50, 56, 78, 106, 134, 156, 171, 182, 189, 194, 197, 199, 200, 200, 200, 199, 198, 196, 194, 191, 188, 185, 182, 178, 174, 169, 164, 159, 154, 149, 144, 138, 133, 128, 123, 118, 114, 110, 106, 102, 99, 96, 93, 91, 89, 87, 86, 84  }, 
103   { 0, 82, 49, 54, 74, 102, 129, 151, 167, 179, 186, 191, 195, 197, 198, 198, 198, 197, 196, 195, 193, 191, 188, 185, 182, 178, 174, 170, 165, 161, 156, 151, 145, 140, 135, 130, 125, 120, 115, 111, 107, 103, 100, 97, 94, 91, 89, 87, 85, 83  }, 
104   { 0, 82, 47, 51, 70, 97, 124, 147, 164, 175, 183, 189, 192, 195, 196, 197, 197, 196, 195, 194, 192, 190, 188, 185, 182, 178, 175, 171, 166, 162, 157, 152, 147, 142, 137, 132, 127, 122, 117, 112, 108, 104, 101, 97, 94, 91, 89, 87, 85, 83  }, 
105   { 0, 83, 46, 49, 67, 93, 120, 143, 160, 172, 180, 186, 190, 192, 194, 195, 195, 195, 194, 193, 191, 189, 187, 185, 182, 179, 175, 172, 167, 163, 159, 154, 149, 144, 139, 134, 129, 124, 119, 114, 110, 106, 102, 98, 95, 92, 89, 87, 85, 83  }, 
106   { 0, 83, 45, 47, 64, 89, 116, 139, 156, 169, 177, 184, 188, 190, 192, 193, 193, 193, 193, 192, 190, 189, 187, 184, 182, 179, 176, 172, 168, 164, 160, 156, 151, 146, 141, 136, 131, 126, 121, 116, 112, 108, 104, 100, 96, 93, 90, 88, 85, 83  }, 
107   { 0, 84, 44, 45, 61, 85, 111, 134, 152, 165, 175, 181, 185, 188, 190, 191, 192, 192, 191, 191, 189, 188, 186, 184, 182, 179, 176, 173, 169, 166, 162, 157, 153, 148, 143, 138, 133, 128, 124, 119, 114, 110, 106, 102, 98, 95, 91, 89, 86, 84  }, 
108   { 0, 85, 43, 43, 58, 81, 107, 131, 149, 162, 172, 178, 183, 186, 188, 190, 190, 190, 190, 189, 188, 187, 186, 184, 182, 179, 176, 173, 170, 167, 163, 159, 155, 150, 145, 141, 136, 131, 126, 121, 117, 112, 108, 104, 100, 96, 93, 90, 87, 85  }, 
109   { 0, 85, 42, 41, 55, 78, 103, 127, 145, 159, 169, 176, 181, 184, 186, 188, 189, 189, 189, 188, 188, 186, 185, 183, 181, 179, 177, 174, 171, 168, 164, 160, 156, 152, 148, 143, 138, 134, 129, 124, 119, 115, 110, 106, 102, 98, 95, 91, 88, 86  } 
110 };
111
112 Int_t (*AliTRDmcmSim::fgPidLut) = *fgPidLutDefault;
113
114 //_____________________________________________________________________________
115 AliTRDmcmSim::AliTRDmcmSim() : TObject()
116   ,fInitialized(kFALSE)
117   ,fMaxTracklets(-1) 
118   ,fDetector(-1)
119   ,fRobPos(-1)
120   ,fMcmPos(-1)
121   ,fRow (-1)
122   ,fNADC(-1)
123   ,fNTimeBin(-1)
124   ,fADCR(NULL)
125   ,fADCF(NULL)
126   ,fMCMT(NULL)
127   ,fTrackletArray(NULL)      
128   ,fZSM(NULL)
129   ,fZSM1Dim(NULL)
130   ,fFeeParam(NULL)
131   ,fTrapConfig(NULL)
132   ,fSimParam(NULL)
133   ,fCommonParam(NULL)
134   ,fCal(NULL)
135   ,fGeo(NULL)
136   ,fDigitsManager(NULL)
137   ,fPedAcc(NULL)
138   ,fGainCounterA(NULL)
139   ,fGainCounterB(NULL)
140   ,fTailAmplLong(NULL)
141   ,fTailAmplShort(NULL)
142   ,fNHits(0)
143   ,fFitReg(NULL)
144 {
145   //
146   // AliTRDmcmSim default constructor
147   // By default, nothing is initialized.
148   // It is necessary to issue Init before use.
149 }
150
151 AliTRDmcmSim::~AliTRDmcmSim() 
152 {
153   //
154   // AliTRDmcmSim destructor
155   //
156
157   if(fInitialized) {
158     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
159       delete [] fADCR[iadc];
160       delete [] fADCF[iadc];
161       delete [] fZSM [iadc];
162     }
163     delete [] fADCR;
164     delete [] fADCF;
165     delete [] fZSM;
166     delete [] fZSM1Dim;
167     delete [] fMCMT;
168  
169     delete [] fPedAcc;
170     delete [] fGainCounterA;
171     delete [] fGainCounterB;
172     delete [] fTailAmplLong;
173     delete [] fTailAmplShort;
174     delete [] fFitReg;
175     
176     fTrackletArray->Delete();
177     delete fTrackletArray;
178     delete fGeo;
179   }
180 }
181
182 void AliTRDmcmSim::Init( Int_t det, Int_t robPos, Int_t mcmPos, Bool_t /* newEvent */ ) 
183 {
184   //
185   // Initialize the class with new geometry information
186   // fADC array will be reused with filled by zero
187   //
188    
189   if (!fInitialized) {
190     fFeeParam      = AliTRDfeeParam::Instance();
191     fTrapConfig    = AliTRDtrapConfig::Instance();
192     fSimParam      = AliTRDSimParam::Instance();
193     fCommonParam   = AliTRDCommonParam::Instance();
194     fCal           = AliTRDcalibDB::Instance();
195     fGeo           = new AliTRDgeometry();
196   }
197
198   fDetector      = det;
199   fRobPos        = robPos;
200   fMcmPos        = mcmPos;
201   fNADC          = fFeeParam->GetNadcMcm();
202   //fNTimeBin      = fCal->GetNumberOfTimeBins();
203   fNTimeBin      = AliTRDSimParam::Instance()->GetNTimeBins();
204   fRow           = fFeeParam->GetPadRowFromMCM( fRobPos, fMcmPos );
205   fMaxTracklets  = fFeeParam->GetMaxNrOfTracklets();  
206   
207   if (!fInitialized) {
208     fADCR    = new Int_t *[fNADC];
209     fADCF    = new Int_t *[fNADC];
210     fZSM     = new Int_t *[fNADC];
211     fZSM1Dim = new Int_t  [fNADC];
212     fGainCounterA = new UInt_t[fNADC];
213     fGainCounterB = new UInt_t[fNADC];
214     for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
215       fADCR[iadc] = new Int_t[fNTimeBin];
216       fADCF[iadc] = new Int_t[fNTimeBin];
217       fZSM [iadc] = new Int_t[fNTimeBin];
218     }
219     
220     // filter registers
221     fPedAcc = new UInt_t[fNADC]; // accumulator for pedestal filter
222     fTailAmplLong = new UShort_t[fNADC];
223     fTailAmplShort = new UShort_t[fNADC];
224     
225     // tracklet calculation
226     fFitReg = new FitReg_t[fNADC]; 
227     fTrackletArray = new TClonesArray("AliTRDtrackletMCM", fMaxTracklets);
228     
229     fMCMT = new UInt_t[fMaxTracklets];
230   }
231
232   fInitialized = kTRUE;
233
234   Reset();
235 }
236
237 void AliTRDmcmSim::Reset()
238 {
239   // Resets the data values and internal filter registers
240   // by re-initialising them
241
242   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
243     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
244       fADCR[iadc][it] = 0;
245       fADCF[iadc][it] = 0;
246       fZSM [iadc][it] = 1;   // Default unread = 1  simulator.SetMakeSDigits("TRD TOF PHOS HMPID EMCAL MUON FMD ZDC PMD T0 VZERO");
247
248     }
249     fZSM1Dim[iadc] = 1;      // Default unread = 1
250     fGainCounterA[iadc] = 0;
251     fGainCounterB[iadc] = 0;
252   }
253   
254   for(Int_t i = 0; i < fMaxTracklets; i++) {
255     fMCMT[i] = 0;
256   }
257   
258   FilterPedestalInit();
259   FilterGainInit();
260   FilterTailInit(fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP)); //??? not really correct if gain filter is active
261 }
262
263 void AliTRDmcmSim::SetNTimebins(Int_t ntimebins) 
264 {
265   fNTimeBin = ntimebins;
266   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
267     delete fADCR[iadc];
268     delete fADCF[iadc];
269     delete fZSM[iadc];
270     fADCR[iadc] = new Int_t[fNTimeBin];
271     fADCF[iadc] = new Int_t[fNTimeBin];
272     fZSM [iadc] = new Int_t[fNTimeBin];
273   }
274 }
275
276 Bool_t AliTRDmcmSim::LoadMCM(AliRunLoader* const runloader, Int_t det, Int_t rob, Int_t mcm) 
277 {
278   // loads the ADC data as obtained from the digitsManager for the specified MCM
279
280   Init(det, rob, mcm);
281
282   if (!runloader) {
283     AliError("No Runloader given");
284     return kFALSE;
285   }
286
287   AliLoader *trdLoader = runloader->GetLoader("TRDLoader");
288   if (!trdLoader) {
289     AliError("Could not get TRDLoader");
290     return kFALSE;
291   }
292
293   Bool_t retval = kTRUE;
294   trdLoader->LoadDigits();
295   fDigitsManager = 0x0;
296   AliTRDdigitsManager *digMgr = new AliTRDdigitsManager();
297   digMgr->SetSDigits(0);
298   digMgr->CreateArrays();
299   digMgr->ReadDigits(trdLoader->TreeD());
300   AliTRDarrayADC *digits = (AliTRDarrayADC*) digMgr->GetDigits(det);
301   if (digits->HasData()) {
302     digits->Expand();
303
304     if (fNTimeBin != digits->GetNtime()) {
305       SetNTimebins(digits->GetNtime());
306     }
307
308     Int_t padrow = fFeeParam->GetPadRowFromMCM(rob, mcm);
309     Int_t padcol = 0;
310     for (Int_t ch = 0; ch < fNADC; ch++) {
311       padcol = GetCol(ch);
312       fZSM1Dim[ch] = 1;
313       if (padcol < 0) {
314         fZSM1Dim[ch] = 0;
315         for (Int_t tb = 0; tb < fNTimeBin; tb++) {
316           fADCR[ch][tb] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
317           fADCF[ch][tb] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
318         }
319       }
320       else {
321         for (Int_t tb = 0; tb < fNTimeBin; tb++) {
322           if (digits->GetData(padrow,padcol, tb) < 0) {
323             fZSM1Dim[ch] = 0; 
324             fADCR[ch][tb] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
325             fADCF[ch][tb] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
326           }
327           else {
328             fADCR[ch][tb] = (digits->GetData(padrow, padcol, tb) + fgAddBaseline) << fgkAddDigits;
329             fADCF[ch][tb] = (digits->GetData(padrow, padcol, tb) + fgAddBaseline) << fgkAddDigits;
330           }
331         }
332       }
333     }
334   }
335   else 
336     retval = kFALSE;
337   
338   delete digMgr;
339   
340   return retval;
341 }
342
343 void AliTRDmcmSim::NoiseTest(Int_t nsamples, Int_t mean, Int_t sigma, Int_t inputGain, Int_t inputTail)
344 {
345   // This function can be used to test the filters. 
346   // It feeds nsamples of ADC values with a gaussian distribution specified by mean and sigma.
347   // The filter chain implemented here consists of:
348   // Pedestal -> Gain -> Tail
349   // With inputGain and inputTail the input to the gain and tail filter, respectively, 
350   // can be chosen where 
351   // 0: noise input
352   // 1: pedestal output
353   // 2: gain output
354   // The input has to be chosen from a stage before. 
355   // The filter behaviour is controlled by the TRAP parameters from AliTRDtrapConfig in the 
356   // same way as in normal simulation.
357   // The functions produces four histograms with the values at the different stages.
358
359   TH1F *h   = new TH1F("noise", "Gaussian Noise;sample;ADC count",
360                        nsamples, 0, nsamples);
361   TH1F *hfp = new TH1F("pedf", "Noise #rightarrow Pedestal filter;sample;ADC count", nsamples, 0, nsamples);
362   TH1F *hfg = new TH1F("pedg", "Pedestal #rightarrow Gain;sample;ADC count", nsamples, 0, nsamples);
363   TH1F *hft = new TH1F("pedt", "Gain #rightarrow Tail;sample;ADC count", nsamples, 0, nsamples);
364   h->SetStats(kFALSE);
365   hfp->SetStats(kFALSE);
366   hfg->SetStats(kFALSE);
367   hft->SetStats(kFALSE);
368   
369   Int_t value;  // ADC count with noise (10 bit)
370   Int_t valuep; // pedestal filter output (12 bit)
371   Int_t valueg; // gain filter output (12 bit)
372   Int_t valuet; // tail filter value (12 bit)
373   
374   for (Int_t i = 0; i < nsamples; i++) {
375     value = (Int_t) gRandom->Gaus(mean, sigma);  // generate noise with gaussian distribution 
376     h->SetBinContent(i, value);
377
378     valuep = FilterPedestalNextSample(1, 0, ((Int_t) value) << 2);
379     
380     if (inputGain == 0)
381       valueg = FilterGainNextSample(1, ((Int_t) value) << 2);
382     else 
383       valueg = FilterGainNextSample(1, valuep); 
384     
385     if (inputTail == 0)
386       valuet = FilterTailNextSample(1, ((Int_t) value) << 2);
387     else if (inputTail == 1)
388       valuet = FilterTailNextSample(1, valuep); 
389     else
390       valuet = FilterTailNextSample(1, valueg); 
391
392     hfp->SetBinContent(i, valuep >> 2);
393     hfg->SetBinContent(i, valueg >> 2);
394     hft->SetBinContent(i, valuet >> 2);
395   }
396
397   TCanvas *c = new TCanvas; 
398   c->Divide(2,2);
399   c->cd(1);
400   h->Draw();
401   c->cd(2);
402   hfp->Draw();
403   c->cd(3);
404   hfg->Draw();
405   c->cd(4);
406   hft->Draw();
407 }
408
409 Bool_t AliTRDmcmSim::CheckInitialized()
410 {
411   //
412   // Check whether object is initialized
413   //
414
415   if( ! fInitialized ) {
416     AliDebug(2, Form ("AliTRDmcmSim is not initialized but function other than Init() is called."));
417   }
418   return fInitialized;
419 }
420
421 void AliTRDmcmSim::Print(Option_t* const option) const
422 {
423   // Prints the data stored and/or calculated for this MCM.
424   // The output is controlled by option which can be a sequence of any of 
425   // the following characters:
426   // R - prints raw ADC data
427   // F - prints filtered data 
428   // H - prints detected hits
429   // T - prints found tracklets
430   // The later stages are only useful when the corresponding calculations 
431   // have been performed.
432
433   printf("MCM %i on ROB %i in detector %i\n", fMcmPos, fRobPos, fDetector);
434
435   TString opt = option;
436   if (opt.Contains("R")) {
437     printf("Raw ADC data (10 bit):\n");
438     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
439       for (Int_t iChannel = 0; iChannel < fNADC; iChannel++) {
440         printf("%5i", fADCR[iChannel][iTimeBin] >> fgkAddDigits);
441       }
442       printf("\n");
443     }
444   }
445
446   if (opt.Contains("F")) {
447     printf("Filtered data (12 bit):\n");
448     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
449       for (Int_t iChannel = 0; iChannel < fNADC; iChannel++) {
450         printf("%5i", fADCF[iChannel][iTimeBin]);
451       }
452       printf("\n");
453     }
454   }
455
456   if (opt.Contains("H")) {
457     printf("Found %i hits:\n", fNHits);
458     for (Int_t iHit = 0; iHit < fNHits; iHit++) {
459       printf("Hit %3i in timebin %2i, ADC %2i has charge %3i and position %3i\n",
460              iHit,  fHits[iHit].fTimebin, fHits[iHit].fChannel, fHits[iHit].fQtot, fHits[iHit].fYpos);
461     }
462   }
463
464   if (opt.Contains("T")) {
465     printf("Tracklets:\n");
466     for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntriesFast(); iTrkl++) {
467       printf("tracklet %i: 0x%08x\n", iTrkl, ((AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl])->GetTrackletWord());
468     }
469   }
470 }
471
472 void AliTRDmcmSim::Draw(Option_t* const option) 
473 {
474   // Plots the data stored in a 2-dim. timebin vs. ADC channel plot.
475   // The option selects what data is plotted and can be a sequence of 
476   // the following characters:
477   // R - plot raw data (default)
478   // F - plot filtered data (meaningless if R is specified)
479   // In addition to the ADC values:
480   // H - plot hits 
481   // T - plot tracklets
482
483   TString opt = option;
484
485   TH2F *hist = new TH2F("mcmdata", Form("Data of MCM %i on ROB %i in detector %i", \
486                                         fMcmPos, fRobPos, fDetector), \
487                         fNADC, -0.5, fNADC-.5, fNTimeBin, -.5, fNTimeBin-.5);
488   hist->GetXaxis()->SetTitle("ADC Channel");
489   hist->GetYaxis()->SetTitle("Timebin");
490   hist->SetStats(kFALSE);
491
492   if (opt.Contains("R")) {
493     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
494       for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
495         hist->SetBinContent(iAdc+1, iTimeBin+1, fADCR[iAdc][iTimeBin] >> fgkAddDigits);
496       }
497     }
498   }
499   else {
500     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
501       for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
502         hist->SetBinContent(iAdc+1, iTimeBin+1, fADCF[iAdc][iTimeBin] >> fgkAddDigits);
503       }
504     }
505   }
506   hist->Draw("colz");
507
508   if (opt.Contains("H")) {
509     TGraph *grHits = new TGraph();
510     for (Int_t iHit = 0; iHit < fNHits; iHit++) {
511       grHits->SetPoint(iHit, 
512                        fHits[iHit].fChannel + 1 + fHits[iHit].fYpos/256., 
513                        fHits[iHit].fTimebin);
514     }
515     grHits->Draw("*");
516   }
517
518   if (opt.Contains("T")) {
519     TLine *trklLines = new TLine[4];
520     for (Int_t iTrkl = 0; iTrkl < fTrackletArray->GetEntries(); iTrkl++) {
521       AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
522       AliTRDtrackletMCM *trkl = (AliTRDtrackletMCM*) (*fTrackletArray)[iTrkl];
523       Float_t offset = pp->GetColPos(fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 19)) + 19 * pp->GetWidthIPad();
524       trklLines[iTrkl].SetX1((offset -  trkl->GetY()) / pp->GetWidthIPad());
525       trklLines[iTrkl].SetY1(0);
526       trklLines[iTrkl].SetX2((offset - (trkl->GetY() + ((Float_t) trkl->GetdY())*140e-4)) / pp->GetWidthIPad());
527       trklLines[iTrkl].SetY2(fNTimeBin - 1);
528       trklLines[iTrkl].SetLineColor(2);
529       trklLines[iTrkl].SetLineWidth(2);
530       printf("Tracklet %i: y = %f, dy = %f, offset = %f\n", iTrkl, trkl->GetY(), (trkl->GetdY() * 140e-4), offset);
531       trklLines[iTrkl].Draw();
532     }
533   }
534 }
535
536 void AliTRDmcmSim::SetData( Int_t iadc, Int_t* const adc )
537 {
538   //
539   // Store ADC data into array of raw data
540   //
541
542   if( !CheckInitialized() ) return;
543
544   if( iadc < 0 || iadc >= fNADC ) {
545                         //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
546     return;
547   }
548
549   for( Int_t it = 0 ;  it < fNTimeBin ; it++ ) {
550     fADCR[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits;
551     fADCF[iadc][it] = (Int_t) (adc[it]) << fgkAddDigits;
552   }
553 }
554
555 void AliTRDmcmSim::SetData( Int_t iadc, Int_t it, Int_t adc )
556 {
557   //
558   // Store ADC data into array of raw data
559   //
560
561   if( !CheckInitialized() ) return;
562
563   if( iadc < 0 || iadc >= fNADC ) {
564     //Log (Form ("Error: iadc is out of range (should be 0 to %d).", fNADC-1));
565     return;
566   }
567
568   fADCR[iadc][it] = adc << fgkAddDigits;
569   fADCF[iadc][it] = adc << fgkAddDigits;
570 }
571
572 void AliTRDmcmSim::SetData(AliTRDarrayADC* const adcArray, AliTRDdigitsManager *digitsManager)
573 {
574   // Set the ADC data from an AliTRDarrayADC
575
576   if (!fInitialized) {
577     AliError("Called uninitialized! Nothing done!");
578     return;
579   }
580
581   fDigitsManager = digitsManager;
582
583   if (fNTimeBin != adcArray->GetNtime()) {
584     SetNTimebins(adcArray->GetNtime());
585   }
586
587   Int_t offset = (fMcmPos % 4) * 21 + (fRobPos % 2) * 84;
588
589 //  Int_t firstAdc = 0;
590 //  Int_t lastAdc = fNADC-1;
591 //
592 //  while (GetCol(firstAdc) < 0) {
593 //    for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
594 //      fADCR[firstAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
595 //      fADCF[firstAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
596 //    }
597 //    firstAdc++;
598 //  }
599 //
600 //  while (GetCol(lastAdc) < 0) {
601 //    for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
602 //      fADCR[lastAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
603 //      fADCF[lastAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
604 //    }
605 //    lastAdc--;
606 //  }
607
608   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
609     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
610       Int_t value = adcArray->GetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin);
611       if (value < 0 || (20-iAdc + offset < 1) || (20-iAdc + offset > 165)) {
612         fADCR[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
613         fADCF[iAdc][iTimeBin] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
614       }
615       else {
616         fADCR[iAdc][iTimeBin] = (adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin)  + fgAddBaseline) << fgkAddDigits;
617         fADCF[iAdc][iTimeBin] = (adcArray->GetData(GetRow(), GetCol(iAdc), iTimeBin)  + fgAddBaseline) << fgkAddDigits;
618       }
619     }
620   }
621 }
622
623 void AliTRDmcmSim::SetDataPedestal( Int_t iadc )
624 {
625   //
626   // Store ADC data into array of raw data
627   //
628
629   if( !CheckInitialized() ) return;
630
631   if( iadc < 0 || iadc >= fNADC ) {
632     return;
633   }
634
635   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
636     fADCR[iadc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
637     fADCF[iadc][it] = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP) + (fgAddBaseline << fgkAddDigits);
638   }
639 }
640
641 Int_t AliTRDmcmSim::GetCol( Int_t iadc )
642 {
643   //
644   // Return column id of the pad for the given ADC channel
645   //
646
647   if( !CheckInitialized() ) 
648     return -1;
649
650   Int_t col = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, iadc);
651   if (col < 0 || col >= fFeeParam->GetNcol()) 
652     return -1;
653   else 
654     return col;
655 }
656
657 Int_t AliTRDmcmSim::ProduceRawStream( UInt_t *buf, Int_t maxSize, UInt_t iEv)
658 {
659   //
660   // Produce raw data stream from this MCM and put in buf
661   // Returns number of words filled, or negative value 
662   // with -1 * number of overflowed words
663   //
664
665   UInt_t  x;
666   Int_t   nw  = 0;  // Number of written words
667   Int_t   of  = 0;  // Number of overflowed words
668   Int_t   rawVer   = fFeeParam->GetRAWversion();
669   Int_t **adc;
670   Int_t   nActiveADC = 0;       // number of activated ADC bits in a word
671
672   if( !CheckInitialized() ) return 0;
673
674   if( fFeeParam->GetRAWstoreRaw() ) {
675     adc = fADCR;
676   } else {
677     adc = fADCF;
678   }
679
680   // Produce MCM header
681   x = (1<<31) | (fRobPos << 28) | (fMcmPos << 24) | ((iEv % 0x100000) << 4) | 0xC;
682
683   if (nw < maxSize) {
684     buf[nw++] = x;
685     //printf("\nMCM header: %X ",x);
686   }
687   else {
688     of++;
689   }
690
691   // Produce ADC mask : nncc cccm mmmm mmmm mmmm mmmm mmmm 1100
692   //                            n : unused , c : ADC count, m : selected ADCs
693   if( rawVer >= 3 ) {
694     x = 0;
695     for( Int_t iAdc = 0 ; iAdc < fNADC ; iAdc++ ) {
696       if( fZSM1Dim[iAdc] == 0 ) { //  0 means not suppressed
697                 x = x | (1 << (iAdc+4) );       // last 4 digit reserved for 1100=0xc
698                 nActiveADC++;           // number of 1 in mmm....m
699       }
700     }
701         x = x | (1 << 30) | ( ( 0x3FFFFFFC ) & (~(nActiveADC) << 25) ) | 0xC;   // nn = 01, ccccc are inverted, 0xc=1100
702         //printf("nActiveADC=%d=%08X, inverted=%X ",nActiveADC,nActiveADC,x );
703
704     if (nw < maxSize) {
705       buf[nw++] = x;
706       //printf("ADC mask: %X nMask=%d ADC data: ",x,nActiveADC);
707     }
708     else {
709       of++;
710     }
711   }
712
713   // Produce ADC data. 3 timebins are packed into one 32 bits word
714   // In this version, different ADC channel will NOT share the same word
715
716   UInt_t aa=0, a1=0, a2=0, a3=0;
717
718   for (Int_t iAdc = 0; iAdc < 21; iAdc++ ) {
719     if( rawVer>= 3 && fZSM1Dim[iAdc] != 0 ) continue; // Zero Suppression, 0 means not suppressed
720     aa = !(iAdc & 1) + 2;
721     for (Int_t iT = 0; iT < fNTimeBin; iT+=3 ) {
722       a1 = ((iT    ) < fNTimeBin ) ? adc[iAdc][iT  ] >> fgkAddDigits : 0;
723       a2 = ((iT + 1) < fNTimeBin ) ? adc[iAdc][iT+1] >> fgkAddDigits : 0;
724       a3 = ((iT + 2) < fNTimeBin ) ? adc[iAdc][iT+2] >> fgkAddDigits : 0;
725       x = (a3 << 22) | (a2 << 12) | (a1 << 2) | aa;
726       if (nw < maxSize) {
727         buf[nw++] = x;
728         //printf("%08X ",x);
729       }
730       else {
731         of++;
732       }
733     }
734   }
735
736   if( of != 0 ) return -of; else return nw;
737 }
738
739 Int_t AliTRDmcmSim::ProduceTrackletStream( UInt_t *buf, Int_t maxSize )
740 {
741   //
742   // Produce tracklet data stream from this MCM and put in buf
743   // Returns number of words filled, or negative value 
744   // with -1 * number of overflowed words
745   //
746
747   Int_t   nw  = 0;  // Number of written words
748   Int_t   of  = 0;  // Number of overflowed words
749     
750   if( !CheckInitialized() ) return 0;
751
752   // Produce tracklet data. A maximum of four 32 Bit words will be written per MCM 
753   // fMCMT is filled continuously until no more tracklet words available
754
755   for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
756     if (nw < maxSize) 
757       buf[nw++] = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet])->GetTrackletWord();
758     else 
759       of++;
760   }
761   
762   if( of != 0 ) return -of; else return nw;
763 }
764
765 void AliTRDmcmSim::Filter()
766 {
767   //
768   // Filter the raw ADC values. The active filter stages and their
769   // parameters are taken from AliTRDtrapConfig.
770   // The raw data is stored separate from the filtered data. Thus, 
771   // it is possible to run the filters on a set of raw values 
772   // sequentially for parameter tuning.
773   //
774
775   if( !CheckInitialized() ) {
776     AliError("got called before initialization! Nothing done!");
777     return;
778   }
779
780   // Apply filters sequentially. Bypass is handled by filters
781   // since counters and internal registers may be updated even 
782   // if the filter is bypassed.
783   // The first filter takes the data from fADCR and 
784   // outputs to fADCF. 
785   
786   // Non-linearity filter not implemented.
787   FilterPedestal();
788   FilterGain();
789   FilterTail();
790   // Crosstalk filter not implemented.
791 }
792
793 void AliTRDmcmSim::FilterPedestalInit() 
794 {
795   // Initializes the pedestal filter assuming that the input has 
796   // been constant for a long time (compared to the time constant).
797
798 //  UShort_t    fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
799   UShort_t    fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
800   UShort_t    shifts[4] = {11, 14, 17, 21}; //??? where to take shifts from?
801
802   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++)
803     fPedAcc[iAdc] = (fSimParam->GetADCbaseline() << 2) * (1<<shifts[fptc]);
804 }
805
806 UShort_t AliTRDmcmSim::FilterPedestalNextSample(Int_t adc, Int_t timebin, UShort_t value)
807 {
808   // Returns the output of the pedestal filter given the input value.
809   // The output depends on the internal registers and, thus, the 
810   // history of the filter.
811
812   UShort_t    fpnp = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); // 0..511 -> 0..127.75, pedestal at the output
813   UShort_t    fptc = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPTC); // 0..3, 0 - fastest, 3 - slowest
814   UShort_t    fpby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPBY); // 0..1 the bypass, active low
815   UShort_t    shifts[4] = {11, 14, 17, 21}; //??? where to come from
816
817   UShort_t accumulatorShifted;
818   Int_t correction;
819   UShort_t inpAdd;
820   
821   inpAdd = value + fpnp;
822
823   if (fpby == 0) //??? before or after update of accumulator
824     return value;
825
826   accumulatorShifted = (fPedAcc[adc] >> shifts[fptc]) & 0x3FF;   // 10 bits
827   if (timebin == 0) // the accumulator is disabled in the drift time
828   {
829     correction = (value & 0x3FF) - accumulatorShifted;
830     fPedAcc[adc] = (fPedAcc[adc] + correction) & 0x7FFFFFFF;             // 31 bits
831   }
832   
833   if (inpAdd <= accumulatorShifted)
834     return 0;
835   else
836   {
837     inpAdd = inpAdd - accumulatorShifted;
838     if (inpAdd > 0xFFF) 
839       return 0xFFF;
840     else 
841       return inpAdd;
842   }
843 }
844
845 void AliTRDmcmSim::FilterPedestal()
846 {
847   //
848   // Apply pedestal filter
849   //
850   // As the first filter in the chain it reads data from fADCR 
851   // and outputs to fADCF. 
852   // It has only an effect if previous samples have been fed to 
853   // find the pedestal. Currently, the simulation assumes that 
854   // the input has been stable for a sufficiently long time.
855
856   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
857     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
858       fADCF[iAdc][iTimeBin] = FilterPedestalNextSample(iAdc, iTimeBin, fADCR[iAdc][iTimeBin]);
859     }
860   }
861 }
862
863 void AliTRDmcmSim::FilterGainInit()
864 {
865   // Initializes the gain filter. In this case, only threshold 
866   // counters are reset.
867
868   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
869     // these are counters which in hardware continue 
870     // until maximum or reset
871     fGainCounterA[iAdc] = 0;
872     fGainCounterB[iAdc] = 0;
873   }
874 }
875
876 UShort_t AliTRDmcmSim::FilterGainNextSample(Int_t adc, UShort_t value)
877 {
878   // Apply the gain filter to the given value.
879   // BEGIN_LATEX O_{i}(t) = #gamma_{i} * I_{i}(t) + a_{i} END_LATEX
880   // The output depends on the internal registers and, thus, the 
881   // history of the filter.
882
883   UShort_t    fgby = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGBY); // bypass, active low
884   UShort_t    fgf  = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGF0 + adc)); // 0x700 + (0 & 0x1ff);
885   UShort_t    fga  = fTrapConfig->GetTrapReg(AliTRDtrapConfig::TrapReg_t(AliTRDtrapConfig::kFGA0 + adc)); // 40;
886   UShort_t    fgta = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTA); // 20;
887   UShort_t    fgtb = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFGTB); // 2060;
888
889   UInt_t tmp;
890
891   value &= 0xFFF;
892   tmp = (value * fgf) >> 11;
893   if (tmp > 0xFFF) tmp = 0xFFF;
894
895   if (fgby == 1)
896     value = AddUintClipping(tmp, fga, 12);
897
898   // Update threshold counters 
899   // not really useful as they are cleared with every new event
900   if ((fGainCounterA[adc] == 0x3FFFFFF) || (fGainCounterB[adc] == 0x3FFFFFF))
901   {
902     if (value >= fgtb) 
903       fGainCounterB[adc]++;
904     else if (value >= fgta) 
905       fGainCounterA[adc]++;
906   }
907
908   return value;
909 }
910
911 void AliTRDmcmSim::FilterGain()
912 {
913   // Read data from fADCF and apply gain filter.
914
915   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
916     for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
917         fADCF[iAdc][iTimeBin] = FilterGainNextSample(iAdc, fADCF[iAdc][iTimeBin]);
918     }
919   }
920 }
921
922 void AliTRDmcmSim::FilterTailInit(Int_t baseline)
923 {
924   // Initializes the tail filter assuming that the input has 
925   // been at the baseline value (configured by FTFP) for a 
926   // sufficiently long time.
927
928   // exponents and weight calculated from configuration
929   UShort_t    alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
930   UShort_t    lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
931   UShort_t    lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
932
933   Float_t lambdaL = lambdaLong  * 1.0 / (1 << 11);
934   Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
935   Float_t alphaL  = alphaLong   * 1.0 / (1 << 11);
936   Float_t qup, qdn;
937   qup = (1 - lambdaL) * (1 - lambdaS);
938   qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
939   Float_t kdc = qup/qdn;
940
941   Float_t kt, ql, qs;
942   UShort_t aout;
943   
944   kt = kdc * baseline;
945   aout = baseline - (UShort_t) kt;
946   ql = lambdaL * (1 - lambdaS) *      alphaL;
947   qs = lambdaS * (1 - lambdaL) * (1 - alphaL);
948
949   for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
950     fTailAmplLong[iAdc]  = (UShort_t) (aout * ql / (ql + qs));
951     fTailAmplShort[iAdc] = (UShort_t) (aout * qs / (ql + qs));
952   }
953 }
954
955 UShort_t AliTRDmcmSim::FilterTailNextSample(Int_t adc, UShort_t value)
956 {
957   // Returns the output of the tail filter for the given input value. 
958   // The output depends on the internal registers and, thus, the 
959   // history of the filter.
960
961   // exponents and weight calculated from configuration
962   UShort_t    alphaLong = 0x3ff & fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTAL); // the weight of the long component
963   UShort_t    lambdaLong = (1 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLL) & 0x1FF); // the multiplier
964   UShort_t    lambdaShort = (0 << 10) | (1 << 9) | (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTLS) & 0x1FF); // the multiplier
965
966   Float_t lambdaL = lambdaLong  * 1.0 / (1 << 11);
967   Float_t lambdaS = lambdaShort * 1.0 / (1 << 11);
968   Float_t alphaL  = alphaLong   * 1.0 / (1 << 11);
969   Float_t qup, qdn;
970   qup = (1 - lambdaL) * (1 - lambdaS);
971   qdn = 1 - lambdaS * alphaL - lambdaL * (1 - alphaL);
972 //  Float_t kdc = qup/qdn;
973
974   UInt_t aDiff;
975   UInt_t alInpv;
976   UShort_t aQ;
977   UInt_t tmp;
978   
979   UShort_t inpVolt = value & 0xFFF;    // 12 bits
980       
981   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFTBY) == 0) // bypass mode, active low
982     return value;
983   else
984   {   
985     // add the present generator outputs
986     aQ = AddUintClipping(fTailAmplLong[adc], fTailAmplShort[adc], 12);
987
988     // calculate the difference between the input the generated signal
989     if (inpVolt > aQ) 
990       aDiff = inpVolt - aQ;
991     else                
992       aDiff = 0;
993
994     // the inputs to the two generators, weighted
995     alInpv = (aDiff * alphaLong) >> 11;
996
997     // the new values of the registers, used next time
998     // long component
999     tmp = AddUintClipping(fTailAmplLong[adc], alInpv, 12);
1000     tmp =  (tmp * lambdaLong) >> 11;
1001     fTailAmplLong[adc] = tmp & 0xFFF;
1002     // short component
1003     tmp = AddUintClipping(fTailAmplShort[adc], aDiff - alInpv, 12);
1004     tmp =  (tmp * lambdaShort) >> 11;
1005     fTailAmplShort[adc] = tmp & 0xFFF;
1006
1007     // the output of the filter
1008     return aDiff;
1009   }
1010 }
1011
1012 void AliTRDmcmSim::FilterTail()
1013 {
1014   // Apply tail filter
1015
1016   for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1017     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
1018       fADCF[iAdc][iTimeBin] = FilterTailNextSample(iAdc, fADCF[iAdc][iTimeBin]);
1019     }
1020   }
1021 }
1022
1023 void AliTRDmcmSim::ZSMapping()
1024 {
1025   //
1026   // Zero Suppression Mapping implemented in TRAP chip
1027   //
1028   // See detail TRAP manual "Data Indication" section:
1029   // http://www.kip.uni-heidelberg.de/ti/TRD/doc/trap/TRAP-UserManual.pdf
1030   //
1031
1032   //??? values should come from TRAPconfig
1033   Int_t eBIS = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIS); // TRAP default = 0x4  (Tis=4)
1034   Int_t eBIT = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIT); // TRAP default = 0x28 (Tit=40)
1035   Int_t eBIL = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIL); // TRAP default = 0xf0
1036                                                                  // (lookup table accept (I2,I1,I0)=(111)
1037                                                                  // or (110) or (101) or (100))
1038   Int_t eBIN = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBIN); // TRAP default = 1 (no neighbor sensitivity)
1039   Int_t ep   = 0; // fTrapConfig->GetTrapReg(AliTRDtrapConfig::kFPNP); //??? really subtracted here
1040
1041   Int_t **adc = fADCF;
1042
1043   if( !CheckInitialized() ) {
1044     AliError("got called uninitialized! Nothing done!");    
1045     return;
1046   }
1047
1048   for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1049     for( Int_t iadc = 1 ; iadc < fNADC-1; iadc++ ) {
1050
1051       // Get ADC data currently in filter buffer
1052       Int_t ap = adc[iadc-1][it] - ep; // previous
1053       Int_t ac = adc[iadc  ][it] - ep; // current
1054       Int_t an = adc[iadc+1][it] - ep; // next
1055
1056       // evaluate three conditions
1057       Int_t i0 = ( ac >=  ap && ac >=  an ) ? 0 : 1; // peak center detection
1058       Int_t i1 = ( ap + ac + an > eBIT )    ? 0 : 1; // cluster
1059       Int_t i2 = ( ac > eBIS )              ? 0 : 1; // absolute large peak
1060
1061       Int_t i = i2 * 4 + i1 * 2 + i0;    // Bit position in lookup table
1062       Int_t d = (eBIL >> i) & 1;         // Looking up  (here d=0 means true
1063                                          // and d=1 means false according to TRAP manual)
1064
1065       fZSM[iadc][it] &= d;
1066       if( eBIN == 0 ) {  // turn on neighboring ADCs
1067         fZSM[iadc-1][it] &= d;
1068         fZSM[iadc+1][it] &= d;
1069       }
1070     }
1071   }
1072
1073   // do 1 dim projection
1074   for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1075     for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1076       fZSM1Dim[iadc] &= fZSM[iadc][it];
1077     }
1078   }
1079 }
1080
1081 void AliTRDmcmSim::DumpData( const char * const f, const char * const target )
1082 {
1083   //
1084   // Dump data stored (for debugging).
1085   // target should contain one or multiple of the following characters
1086   //   R   for raw data
1087   //   F   for filtered data
1088   //   Z   for zero suppression map
1089   //   S   Raw dat astream
1090   // other characters are simply ignored
1091   //
1092
1093   UInt_t tempbuf[1024];
1094
1095   if( !CheckInitialized() ) return;
1096
1097   std::ofstream of( f, std::ios::out | std::ios::app );
1098   of << Form("AliTRDmcmSim::DumpData det=%03d sm=%02d stack=%d layer=%d rob=%d mcm=%02d\n",
1099              fDetector, fGeo->GetSector(fDetector), fGeo->GetStack(fDetector), 
1100              fGeo->GetSector(fDetector), fRobPos, fMcmPos );
1101
1102   for( Int_t t=0 ; target[t] != 0 ; t++ ) {
1103     switch( target[t] ) {
1104     case 'R' :
1105     case 'r' :
1106       of << Form("fADCR (raw ADC data)\n");
1107       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1108         of << Form("  ADC %02d: ", iadc);
1109         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1110           of << Form("% 4d",  fADCR[iadc][it]);
1111         }
1112         of << Form("\n");
1113       }
1114       break;
1115     case 'F' :
1116     case 'f' :
1117       of << Form("fADCF (filtered ADC data)\n");
1118       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1119         of << Form("  ADC %02d: ", iadc);
1120         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1121           of << Form("% 4d",  fADCF[iadc][it]);
1122         }
1123         of << Form("\n");
1124       }
1125       break;
1126     case 'Z' :
1127     case 'z' :
1128       of << Form("fZSM and fZSM1Dim (Zero Suppression Map)\n");
1129       for( Int_t iadc = 0 ; iadc < fNADC; iadc++ ) {
1130         of << Form("  ADC %02d: ", iadc);
1131         if( fZSM1Dim[iadc] == 0 ) { of << " R   " ; } else { of << " .   "; } // R:read .:suppressed
1132         for( Int_t it = 0 ; it < fNTimeBin ; it++ ) {
1133           if( fZSM[iadc][it] == 0 ) { of << " R"; } else { of << " ."; } // R:read .:suppressed
1134         }
1135         of << Form("\n");
1136       }
1137       break;
1138     case 'S' :
1139     case 's' :
1140       Int_t s = ProduceRawStream( tempbuf, 1024 ); 
1141       of << Form("Stream for Raw Simulation size=%d rawver=%d\n", s, fFeeParam->GetRAWversion());
1142       of << Form("  address  data\n");
1143       for( Int_t i = 0 ; i < s ; i++ ) {
1144         of << Form("  %04x     %08x\n", i, tempbuf[i]);
1145       }
1146     }
1147   }
1148 }
1149
1150 void AliTRDmcmSim::AddHitToFitreg(Int_t adc, UShort_t timebin, UShort_t qtot, Short_t ypos, Int_t label) 
1151 {
1152   // Add the given hit to the fit register which is lateron used for 
1153   // the tracklet calculation. 
1154   // In addition to the fit sums in the fit register MC information 
1155   // is stored.
1156
1157   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0)) && 
1158       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0)))
1159     fFitReg[adc].fQ0 += qtot;
1160   
1161   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1)) && 
1162       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1)))
1163     fFitReg[adc].fQ1 += qtot;
1164   
1165   if ((timebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS) ) && 
1166       (timebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE)))
1167   {
1168     fFitReg[adc].fSumX  += timebin;
1169     fFitReg[adc].fSumX2 += timebin*timebin;
1170     fFitReg[adc].fNhits++;
1171     fFitReg[adc].fSumY  += ypos;
1172     fFitReg[adc].fSumY2 += ypos*ypos;
1173     fFitReg[adc].fSumXY += timebin*ypos;
1174   }
1175
1176   // register hits (MC info)
1177   fHits[fNHits].fChannel = adc;
1178   fHits[fNHits].fQtot = qtot;
1179   fHits[fNHits].fYpos = ypos;
1180   fHits[fNHits].fTimebin = timebin;
1181   fHits[fNHits].fLabel = label;
1182   fNHits++;
1183 }
1184
1185 void AliTRDmcmSim::CalcFitreg() 
1186 {
1187   // Preprocessing.
1188   // Detect the hits and fill the fit registers.
1189   // Requires 12-bit data from fADCF which means Filter() 
1190   // has to be called before even if all filters are bypassed.
1191
1192   //???
1193   // TRAP parameters:
1194   const UShort_t lutPos[128] = {   // move later to some other file
1195     0,  1,  1,  2,  2,  3,  3,  4,  4,  5,  5,  6,  6,  7,  7,  8,  8,  9,  9, 10, 10, 11, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15,
1196     16, 16, 16, 17, 17, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 24, 25, 25, 25, 26, 26, 26, 26,
1197     27, 27, 27, 27, 27, 27, 27, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 27, 27, 27, 27, 26,
1198     26, 26, 26, 25, 25, 25, 24, 24, 23, 23, 22, 22, 21, 21, 20, 20, 19, 18, 18, 17, 17, 16, 15, 14, 13, 12, 11, 10,  9,  8,  7,  7};
1199   
1200   //??? to be clarified:
1201   UInt_t adcMask = 0xffffffff;
1202   
1203   UShort_t timebin, adcch, adcLeft, adcCentral, adcRight, hitQual, timebin1, timebin2, qtotTemp;
1204   Short_t ypos, fromLeft, fromRight, found;
1205   UShort_t qTotal[19]; // the last is dummy
1206   UShort_t marked[6], qMarked[6], worse1, worse2;
1207   
1208   timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFS); 
1209   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0) 
1210       < timebin1)
1211     timebin1 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0);
1212   timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFE); 
1213   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1) 
1214       > timebin2)
1215     timebin2 = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1);
1216
1217   // reset the fit registers
1218   fNHits = 0; 
1219   for (adcch = 0; adcch < fNADC-2; adcch++) // due to border channels
1220   {
1221     fFitReg[adcch].fNhits = 0;
1222     fFitReg[adcch].fQ0    = 0;
1223     fFitReg[adcch].fQ1    = 0;
1224     fFitReg[adcch].fSumX  = 0;
1225     fFitReg[adcch].fSumY  = 0;
1226     fFitReg[adcch].fSumX2 = 0;
1227     fFitReg[adcch].fSumY2 = 0;
1228     fFitReg[adcch].fSumXY = 0;
1229   }
1230   
1231   for (timebin = timebin1; timebin < timebin2; timebin++)
1232   {
1233     // first find the hit candidates and store the total cluster charge in qTotal array
1234     // in case of not hit store 0 there.
1235     for (adcch = 0; adcch < fNADC-2; adcch++) {
1236       if ( ( (adcMask >> adcch) & 7) == 7) //??? all 3 channels are present in case of ZS
1237       {
1238         adcLeft  = fADCF[adcch  ][timebin];
1239         adcCentral  = fADCF[adcch+1][timebin];
1240         adcRight = fADCF[adcch+2][timebin];
1241         if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVBY) == 1) 
1242           hitQual = ( (adcLeft * adcRight) < 
1243                        (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPVT) * adcCentral) );
1244         else            
1245           hitQual = 1;
1246         // The accumulated charge is with the pedestal!!!
1247         qtotTemp = adcLeft + adcCentral + adcRight;
1248         if ( (hitQual) &&
1249              (qtotTemp >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)) &&
1250              (adcLeft <= adcCentral) &&
1251              (adcCentral > adcRight) )
1252           qTotal[adcch] = qtotTemp;
1253         else
1254           qTotal[adcch] = 0;
1255       }
1256       else
1257         qTotal[adcch] = 0; //jkl
1258       AliDebug(10,Form("ch %2d   qTotal %5d",adcch, qTotal[adcch]));
1259     }
1260
1261     fromLeft = -1;
1262     adcch = 0;
1263     found = 0;
1264     marked[4] = 19; // invalid channel
1265     marked[5] = 19; // invalid channel
1266     qTotal[19] = 0;
1267     while ((adcch < 16) && (found < 3))
1268     {
1269       if (qTotal[adcch] > 0)
1270       {
1271         fromLeft = adcch;
1272         marked[2*found+1]=adcch;
1273         found++;
1274       }
1275       adcch++;
1276     }
1277     
1278     fromRight = -1;
1279     adcch = 18;
1280     found = 0;
1281     while ((adcch > 2) && (found < 3))
1282     {
1283       if (qTotal[adcch] > 0)
1284       {
1285         marked[2*found]=adcch;
1286         found++;
1287         fromRight = adcch;
1288       }
1289       adcch--;
1290     }
1291
1292     AliDebug(10,Form("Fromleft=%d, Fromright=%d",fromLeft, fromRight));
1293     // here mask the hit candidates in the middle, if any
1294     if ((fromLeft >= 0) && (fromRight >= 0) && (fromLeft < fromRight))
1295       for (adcch = fromLeft+1; adcch < fromRight; adcch++)
1296         qTotal[adcch] = 0;
1297     
1298     found = 0;
1299     for (adcch = 0; adcch < 19; adcch++)
1300       if (qTotal[adcch] > 0) found++;
1301     // NOT READY
1302
1303     if (found > 4) // sorting like in the TRAP in case of 5 or 6 candidates!
1304     {
1305       if (marked[4] == marked[5]) marked[5] = 19;
1306       for (found=0; found<6; found++)
1307       {
1308         qMarked[found] = qTotal[marked[found]] >> 4;
1309         AliDebug(10,Form("ch_%d qTotal %d qTotals %d",marked[found],qTotal[marked[found]],qMarked[found]));
1310       }
1311       
1312       Sort6To2Worst(marked[0], marked[3], marked[4], marked[1], marked[2], marked[5],
1313                     qMarked[0],
1314                     qMarked[3],
1315                     qMarked[4],
1316                     qMarked[1],
1317                     qMarked[2],
1318                     qMarked[5],
1319                     &worse1, &worse2);
1320       // Now mask the two channels with the smallest charge
1321       if (worse1 < 19)
1322       {
1323         qTotal[worse1] = 0;
1324         AliDebug(10,Form("Kill ch %d\n",worse1));
1325       }
1326       if (worse2 < 19)
1327       {
1328         qTotal[worse2] = 0;
1329         AliDebug(10,Form("Kill ch %d\n",worse2));
1330       }
1331     }
1332     
1333     for (adcch = 0; adcch < 19; adcch++) {
1334       if (qTotal[adcch] > 0) // the channel is marked for processing
1335       {
1336         adcLeft  = fADCF[adcch  ][timebin];
1337         adcCentral  = fADCF[adcch+1][timebin];
1338         adcRight = fADCF[adcch+2][timebin];
1339         // hit detected, in TRAP we have 4 units and a hit-selection, here we proceed all channels!
1340         // subtract the pedestal TPFP, clipping instead of wrapping
1341         
1342         Int_t regTPFP = fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPFP);
1343         AliDebug(10, Form("Hit found, time=%d, adcch=%d/%d/%d, adc values=%d/%d/%d, regTPFP=%d, TPHT=%d\n",
1344                timebin, adcch, adcch+1, adcch+2, adcLeft, adcCentral, adcRight, regTPFP, 
1345                fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPHT)));
1346
1347         if (adcLeft  < regTPFP) adcLeft  = 0; else adcLeft  -= regTPFP;
1348         if (adcCentral  < regTPFP) adcCentral  = 0; else adcCentral  -= regTPFP;
1349         if (adcRight < regTPFP) adcRight = 0; else adcRight -= regTPFP;
1350
1351         // Calculate the center of gravity
1352         // checking for adcCentral != 0 (in case of "bad" configuration)
1353         if (adcCentral == 0)
1354           continue;
1355         ypos = 128*(adcLeft - adcRight) / adcCentral;
1356         if (ypos < 0) ypos = -ypos;
1357         // make the correction using the LUT
1358         ypos = ypos + lutPos[ypos & 0x7F];
1359         if (adcLeft > adcRight) ypos = -ypos;
1360
1361         // label calculation
1362         Int_t mcLabel = -1;
1363         if (fDigitsManager) {
1364           Int_t label[9] = { 0 }; // up to 9 different labels possible
1365           Int_t count[9] = { 0 };
1366           Int_t maxIdx = -1;
1367           Int_t maxCount = 0;
1368           Int_t nLabels = 0;
1369           Int_t padcol[3]; 
1370           padcol[0] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch);
1371           padcol[1] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch+1);
1372           padcol[2] = fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, adcch+2);
1373           Int_t padrow = fFeeParam->GetPadRowFromMCM(fRobPos, fMcmPos);
1374           for (Int_t iDict = 0; iDict < 3; iDict++) {
1375             if (!fDigitsManager->UsesDictionaries() || fDigitsManager->GetDictionary(fDetector, iDict) == 0) {
1376               AliError("Cannot get dictionary");
1377               continue;
1378             }
1379             AliTRDarrayDictionary *dict = (AliTRDarrayDictionary*) fDigitsManager->GetDictionary(fDetector, iDict);
1380             if (dict->GetDim() == 0) {
1381               AliError(Form("Dictionary %i of det. %i has dim. 0", fDetector, iDict));
1382               continue;
1383             }
1384             dict->Expand();
1385             for (Int_t iPad = 0; iPad < 3; iPad++) {
1386               if (padcol[iPad] < 0) 
1387                 continue;
1388               Int_t currLabel = dict->GetData(padrow, padcol[iPad], timebin); //fDigitsManager->GetTrack(iDict, padrow, padcol, timebin, fDetector);
1389               AliDebug(10, Form("Read label: %4i for det: %3i, row: %i, col: %i, tb: %i\n", currLabel, fDetector, padrow, padcol[iPad], timebin));
1390               for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
1391                 if (currLabel == label[iLabel]) {
1392                   count[iLabel]++;
1393                   if (count[iLabel] > maxCount) {
1394                     maxCount = count[iLabel];
1395                     maxIdx = iLabel;
1396                   }
1397                   currLabel = 0;
1398                   break;
1399                 }
1400               } 
1401               if (currLabel > 0) {
1402                 label[nLabels++] = currLabel;
1403               }
1404             }
1405           }
1406           if (maxIdx >= 0)
1407             mcLabel = label[maxIdx];
1408         }
1409
1410         // add the hit to the fitregister
1411         AddHitToFitreg(adcch, timebin, qTotal[adcch], ypos, mcLabel);
1412       }
1413     }
1414   }
1415 }
1416
1417 void AliTRDmcmSim::TrackletSelection() 
1418 {
1419   // Select up to 4 tracklet candidates from the fit registers  
1420   // and assign them to the CPUs.
1421
1422   UShort_t adcIdx, i, j, ntracks, tmp;
1423   UShort_t trackletCand[18][2]; // store the adcch[0] and number of hits[1] for all tracklet candidates
1424
1425   ntracks = 0;
1426   for (adcIdx = 0; adcIdx < 18; adcIdx++) // ADCs
1427     if ( (fFitReg[adcIdx].fNhits 
1428           >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCL)) &&
1429          (fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits
1430           >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPCT)))
1431     {
1432       trackletCand[ntracks][0] = adcIdx;
1433       trackletCand[ntracks][1] = fFitReg[adcIdx].fNhits+fFitReg[adcIdx+1].fNhits;
1434       AliDebug(10,Form("%d  %2d %4d\n", ntracks, trackletCand[ntracks][0], trackletCand[ntracks][1]));
1435       ntracks++;
1436     };
1437
1438   for (i=0; i<ntracks;i++) 
1439     AliDebug(10,Form("%d %d %d\n",i,trackletCand[i][0], trackletCand[i][1]));
1440
1441   if (ntracks > 4)
1442   {
1443     // primitive sorting according to the number of hits
1444     for (j = 0; j < (ntracks-1); j++)
1445     {
1446       for (i = j+1; i < ntracks; i++)
1447       {
1448         if ( (trackletCand[j][1]  < trackletCand[i][1]) ||
1449              ( (trackletCand[j][1] == trackletCand[i][1]) && (trackletCand[j][0] < trackletCand[i][0]) ) )
1450         {
1451           // swap j & i
1452           tmp = trackletCand[j][1];
1453           trackletCand[j][1] = trackletCand[i][1];
1454           trackletCand[i][1] = tmp;
1455           tmp = trackletCand[j][0];
1456           trackletCand[j][0] = trackletCand[i][0];
1457           trackletCand[i][0] = tmp;
1458         }
1459       }
1460     }
1461     ntracks = 4; // cut the rest, 4 is the max
1462   }
1463   // else is not necessary to sort
1464   
1465   // now sort, so that the first tracklet going to CPU0 corresponds to the highest adc channel - as in the TRAP
1466   for (j = 0; j < (ntracks-1); j++)
1467   {
1468     for (i = j+1; i < ntracks; i++)
1469     {
1470       if (trackletCand[j][0] < trackletCand[i][0])
1471       {
1472         // swap j & i
1473         tmp = trackletCand[j][1];
1474         trackletCand[j][1] = trackletCand[i][1];
1475         trackletCand[i][1] = tmp;
1476         tmp = trackletCand[j][0];
1477         trackletCand[j][0] = trackletCand[i][0];
1478         trackletCand[i][0] = tmp;
1479       }
1480     }
1481   }
1482   for (i = 0; i < ntracks; i++)  // CPUs with tracklets.
1483     fFitPtr[i] = trackletCand[i][0]; // pointer to the left channel with tracklet for CPU[i]
1484   for (i = ntracks; i < 4; i++)  // CPUs without tracklets
1485     fFitPtr[i] = 31;            // pointer to the left channel with tracklet for CPU[i] = 31 (invalid)
1486   AliDebug(10,Form("found %i tracklet candidates\n", ntracks));
1487   for (i = 0; i < 4; i++)
1488     AliDebug(10,Form("fitPtr[%i]: %i\n", i, fFitPtr[i]));
1489 }
1490
1491 void AliTRDmcmSim::FitTracklet()
1492 {
1493   // Perform the actual tracklet fit based on the fit sums 
1494   // which have been filled in the fit registers. 
1495
1496   // parameters in fitred.asm (fit program)
1497   Int_t decPlaces = 5;
1498   Int_t rndAdd = 0;
1499   if (decPlaces >  1) 
1500     rndAdd = (1 << (decPlaces-1)) + 1;
1501   else if (decPlaces == 1)
1502     rndAdd = 1;
1503   Int_t ndriftDp = 5;  // decimal places for drift time
1504   Long64_t shift = ((Long64_t) 1 << 32);
1505
1506
1507   // calculated in fitred.asm
1508   Int_t padrow = ((fRobPos >> 1) << 2) | (fMcmPos >> 2);
1509   Int_t yoffs = (((((fRobPos & 0x1) << 2) + (fMcmPos & 0x3)) * 18) << 8) - 
1510     ((18*4*2 - 18*2 - 1) << 7);
1511   yoffs = yoffs << decPlaces; // holds position of ADC channel 1
1512   Int_t layer = fDetector % 6;
1513   UInt_t scaleY = (UInt_t) ((0.635 + 0.03 * layer)/(256.0 * 160.0e-4) * shift);
1514   UInt_t scaleD = (UInt_t) ((0.635 + 0.03 * layer)/(256.0 * 140.0e-4) * shift);
1515   // previously taken from geometry:
1516   // UInt_t scaleYold = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 160e-4)));
1517   // UInt_t scaleDold = (UInt_t) (shift * (pp->GetWidthIPad() / (256 * 140e-4)));
1518
1519
1520   // should come from trapConfig (DMEM) 
1521   AliTRDpadPlane *pp = fGeo->GetPadPlane(fDetector);
1522   Float_t scaleSlope = (256 / pp->GetWidthIPad()) * (1 << decPlaces); // only used for calculation of corrections and cut
1523   Int_t ndrift   = 20 << ndriftDp; //??? value in simulation?
1524   Int_t deflCorr = (Int_t) (TMath::Tan(fCommonParam->GetOmegaTau(fCal->GetVdriftAverage(fDetector))) * fGeo->CdrHght() * scaleSlope); // -370;
1525   Int_t tiltCorr = (Int_t) (pp->GetRowPos(padrow) / fGeo->GetTime0(fDetector % 6) * fGeo->CdrHght() * scaleSlope * 
1526                             TMath::Tan(pp->GetTiltingAngle() / 180. * TMath::Pi()));
1527 //  printf("vdrift av.: %f\n", fCal->GetVdriftAverage(fDetector));
1528 //  printf("chamber height: %f\n", fGeo->CdrHght());
1529 //  printf("omega tau: %f\n", fCommonParam->GetOmegaTau(fCal->GetVdriftAverage(fDetector)));
1530 //  printf("deflection correction: %i\n", deflCorr);
1531   Float_t ptcut = 2.3;
1532   AliMagF* fld = (AliMagF *) TGeoGlobalMagField::Instance()->GetField();
1533   Double_t bz = 0;
1534   if (fld) {
1535     bz       = 0.1 * fld->SolenoidField();   // kGauss -> Tesla
1536   }
1537 //  printf("Bz: %f\n", bz);
1538   Float_t x0 = fGeo->GetTime0(fDetector % 6);
1539   Float_t y0 = pp->GetColPos(fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 10));
1540   Float_t alphaMax = TMath::ASin( (TMath::Sqrt(TMath::Power(x0/100., 2) + TMath::Power(y0/100., 2)) * 
1541                                    0.3 * TMath::Abs(bz) ) / (2 * ptcut));
1542 //  printf("alpha max: %f\n", alphaMax * 180/TMath::Pi());
1543   Int_t minslope = -1 * (Int_t) (fGeo->CdrHght() * TMath::Tan(TMath::ATan(y0/x0) + alphaMax) / 140.e-4);
1544   Int_t maxslope = -1 * (Int_t) (fGeo->CdrHght() * TMath::Tan(TMath::ATan(y0/x0) - alphaMax) / 140.e-4);
1545
1546
1547   // local variables for calculation
1548   Long64_t mult, temp, denom; //???
1549   UInt_t q0, q1, qTotal;          // charges in the two windows and total charge
1550   UShort_t nHits;                 // number of hits
1551   Int_t slope, offset;            // slope and offset of the tracklet
1552   Int_t sumX, sumY, sumXY, sumX2; // fit sums from fit registers
1553   //int32_t SumY2;                // not used in the current TRAP program
1554   FitReg_t *fit0, *fit1;          // pointers to relevant fit registers
1555   
1556 //  const uint32_t OneDivN[32] = {  // 2**31/N : exactly like in the TRAP, the simple division here gives the same result!
1557 //      0x00000000, 0x80000000, 0x40000000, 0x2AAAAAA0, 0x20000000, 0x19999990, 0x15555550, 0x12492490,
1558 //      0x10000000, 0x0E38E380, 0x0CCCCCC0, 0x0BA2E8B0, 0x0AAAAAA0, 0x09D89D80, 0x09249240, 0x08888880,
1559 //      0x08000000, 0x07878780, 0x071C71C0, 0x06BCA1A0, 0x06666660, 0x06186180, 0x05D17450, 0x0590B210,
1560 //      0x05555550, 0x051EB850, 0x04EC4EC0, 0x04BDA120, 0x04924920, 0x0469EE50, 0x04444440, 0x04210840};
1561
1562   for (Int_t cpu = 0; cpu < 4; cpu++) {
1563     if (fFitPtr[cpu] == 31)
1564     {
1565       fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker(); 
1566     }
1567     else
1568     {
1569       fit0 = &fFitReg[fFitPtr[cpu]  ];
1570       fit1 = &fFitReg[fFitPtr[cpu]+1]; // next channel
1571
1572       mult = 1;
1573       mult = mult << (32 + decPlaces);
1574       mult = -mult;
1575
1576       // Merging
1577       nHits   = fit0->fNhits + fit1->fNhits; // number of hits
1578       sumX    = fit0->fSumX  + fit1->fSumX;
1579       sumX2   = fit0->fSumX2 + fit1->fSumX2;
1580       denom   = nHits*sumX2 - sumX*sumX;
1581
1582       mult    = mult / denom; // exactly like in the TRAP program
1583       q0      = fit0->fQ0    + fit1->fQ0;
1584       q1      = fit0->fQ1    + fit1->fQ1;
1585       sumY    = fit0->fSumY  + fit1->fSumY  + 256*fit1->fNhits;
1586       sumXY   = fit0->fSumXY + fit1->fSumXY + 256*fit1->fSumX;
1587
1588       slope   = nHits*sumXY - sumX * sumY;
1589       AliDebug(5, Form("slope from fitreg: %i", slope));
1590       offset  = sumX2*sumY  - sumX * sumXY;
1591       temp    = mult * slope;
1592       slope   = temp >> 32; // take the upper 32 bits
1593       slope   = -slope;
1594       temp    = mult * offset;
1595       offset  = temp >> 32; // take the upper 32 bits
1596
1597       offset = offset + yoffs;
1598       AliDebug(5, Form("slope: %i, slope * ndrift: %i, deflCorr: %i, tiltCorr: %i", slope, slope * ndrift, deflCorr, tiltCorr));
1599       slope  = ((slope * ndrift) >> ndriftDp) + deflCorr + tiltCorr;
1600       offset = offset - (fFitPtr[cpu] << (8 + decPlaces));
1601       
1602       AliDebug(5, Form("Det: %3i, ROB: %i, MCM: %2i: deflection: %i, min: %i, max: %i", fDetector, fRobPos, fMcmPos, slope, minslope, maxslope));
1603       temp    = slope;
1604       temp    = temp * scaleD;
1605       slope   = (temp >> 32);
1606       AliDebug(5, Form("slope after scaling: %i", slope));
1607
1608       temp    = offset;
1609       temp    = temp * scaleY;
1610       offset  = (temp >> 32);
1611         
1612       // rounding, like in the TRAP
1613       slope   = (slope  + rndAdd) >> decPlaces;
1614       AliDebug(5, Form("slope after shifting: %i", slope));
1615       offset  = (offset + rndAdd) >> decPlaces;
1616
1617       Bool_t rejected = kFALSE;
1618       if ((slope < minslope) || (slope > maxslope))
1619         rejected = kTRUE;
1620
1621       if (rejected && GetApplyCut())
1622       {
1623         fMCMT[cpu] = 0x10001000; //??? AliTRDfeeParam::GetTrackletEndmarker();
1624       }
1625       else
1626       {
1627         if (slope > 63 || slope < -64) { // wrapping in TRAP!
1628           AliError(Form("Overflow in slope: %i, tracklet discarded!", slope));
1629           fMCMT[cpu] = 0x10001000;
1630           continue;
1631         }
1632
1633         slope   = slope  &   0x7F; // 7 bit
1634         
1635         if (offset > 0xfff || offset < -0xfff) 
1636           AliWarning("Overflow in offset");
1637         offset  = offset & 0x1FFF; // 13 bit
1638
1639         Float_t length = TMath::Sqrt(1 + (pp->GetRowPos(padrow) * pp->GetRowPos(padrow) +
1640                                           (fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 10) * pp->GetWidthIPad() *
1641                                            fFeeParam->GetPadColFromADC(fRobPos, fMcmPos, 10) * pp->GetWidthIPad())) /
1642                                      (fGeo->GetTime0(fDetector % 6)*fGeo->GetTime0(fDetector % 6)));
1643
1644         //        qTotal  = (q1 / nHits) >> 1;
1645         qTotal = GetPID(q0/length/fgChargeNorm, q1/length/fgChargeNorm);
1646         if (qTotal > 0xff)
1647           AliWarning("Overflow in charge");
1648         qTotal  = qTotal & 0xFF; // 8 bit, exactly like in the TRAP program
1649         
1650         // assemble and store the tracklet word
1651         fMCMT[cpu] = (qTotal << 24) | (padrow << 20) | (slope << 13) | offset;
1652
1653         // calculate MC label
1654         Int_t mcLabel = -1;
1655         Int_t nHits0 = 0;
1656         Int_t nHits1 = 0;
1657         if (fDigitsManager) {
1658           Int_t label[30] = {0}; // up to 30 different labels possible
1659           Int_t count[30] = {0};
1660           Int_t maxIdx = -1;
1661           Int_t maxCount = 0;
1662           Int_t nLabels = 0;
1663           for (Int_t iHit = 0; iHit < fNHits; iHit++) {
1664             if ((fHits[iHit].fChannel - fFitPtr[cpu] < 0) ||
1665                 (fHits[iHit].fChannel - fFitPtr[cpu] > 1))
1666               continue;
1667
1668             // counting contributing hits
1669             if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS0) &&
1670                 fHits[iHit].fTimebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE0))
1671               nHits0++;
1672             if (fHits[iHit].fTimebin >= fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQS1) &&
1673                 fHits[iHit].fTimebin <  fTrapConfig->GetTrapReg(AliTRDtrapConfig::kTPQE1))
1674               nHits1++;
1675
1676             Int_t currLabel = fHits[iHit].fLabel;
1677             for (Int_t iLabel = 0; iLabel < nLabels; iLabel++) {
1678               if (currLabel == label[iLabel]) {
1679                 count[iLabel]++;
1680                 if (count[iLabel] > maxCount) {
1681                   maxCount = count[iLabel];
1682                   maxIdx = iLabel;
1683                 }
1684                 currLabel = 0;
1685                 break;
1686               }
1687             }
1688             if (currLabel > 0) {
1689               label[nLabels++] = currLabel;
1690             }
1691           }
1692           if (maxIdx >= 0)
1693             mcLabel = label[maxIdx];
1694         }
1695         new ((*fTrackletArray)[fTrackletArray->GetEntriesFast()]) AliTRDtrackletMCM((UInt_t) fMCMT[cpu], fDetector*2 + fRobPos%2, fRobPos, fMcmPos);
1696         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetLabel(mcLabel);
1697
1698        
1699         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits(fit0->fNhits + fit1->fNhits);
1700         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits0(nHits0);
1701         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetNHits1(nHits1);
1702         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ0(q0);
1703         ((AliTRDtrackletMCM*) (*fTrackletArray)[fTrackletArray->GetEntriesFast()-1])->SetQ1(q1);
1704       }
1705     }
1706   }
1707 }
1708
1709 Int_t AliTRDmcmSim::GetPID(Float_t q0, Float_t q1) 
1710 {
1711   // get PID from accumulated charges q0 and q1
1712
1713   Int_t binQ0 = (Int_t) (q0 * fgPidNBinsQ0) + 1;
1714   Int_t binQ1 = (Int_t) (q1 * fgPidNBinsQ1) + 1;
1715   binQ0 = binQ0 >= fgPidNBinsQ0 ? fgPidNBinsQ0-1 : binQ0;
1716   binQ1 = binQ1 >= fgPidNBinsQ0 ? fgPidNBinsQ0-1 : binQ1;
1717
1718   return fgPidLut[binQ0*fgPidNBinsQ1+binQ1];
1719 }
1720
1721 void AliTRDmcmSim::SetPIDlut(Int_t *lut, Int_t nbinsq0, Int_t nbinsq1)
1722 {
1723   // set a user-defined PID LUT
1724
1725   if (fgPidLutDelete)
1726     delete [] fgPidLut;
1727
1728   fgPidLutDelete = kFALSE;
1729   fgPidLut = lut;
1730   fgPidNBinsQ0 = nbinsq0;
1731   fgPidNBinsQ1 = nbinsq1;
1732 }
1733
1734 void AliTRDmcmSim::SetPIDlut(TH2F *lut)
1735 {
1736   // set a user-defined PID LUT from a 2D histogram
1737
1738   if (fgPidLutDelete)
1739     delete [] fgPidLut;
1740
1741   fgPidNBinsQ0 = lut->GetNbinsX();
1742   fgPidNBinsQ1 = lut->GetNbinsY();
1743
1744   fgPidLut = new Int_t[fgPidNBinsQ0*fgPidNBinsQ1];
1745
1746   for (Int_t ix = 0; ix < fgPidNBinsQ0; ix++) {
1747     for (Int_t iy = 0; iy < fgPidNBinsQ1; iy++) {
1748       fgPidLut[ix*fgPidNBinsQ1 + iy] = (Int_t) (256. * lut->GetBinContent(ix, iy));
1749     }
1750   }
1751
1752   fgPidLutDelete = kTRUE;
1753 }
1754
1755 void AliTRDmcmSim::SetPIDlutDefault()
1756 {
1757   // use the default PID LUT
1758
1759   if (fgPidLutDelete )
1760     delete [] fgPidLut;
1761
1762   fgPidLutDelete = kFALSE;
1763   fgPidLut = *fgPidLutDefault;
1764   fgPidNBinsQ0 = 40;
1765   fgPidNBinsQ1 = 50;
1766 }
1767
1768 void AliTRDmcmSim::Tracklet()
1769 {
1770   // Run the tracklet calculation by calling sequentially:
1771   // CalcFitreg(); TrackletSelection(); FitTracklet()
1772   // and store the tracklets 
1773
1774   if (!fInitialized) {
1775     AliError("Called uninitialized! Nothing done!");
1776     return;
1777   }
1778
1779   fTrackletArray->Delete();
1780
1781   CalcFitreg();
1782   if (fNHits == 0)
1783     return;
1784   TrackletSelection();
1785   FitTracklet();
1786 }
1787
1788 Bool_t AliTRDmcmSim::StoreTracklets() 
1789 {
1790   // store the found tracklets via the loader
1791
1792   if (fTrackletArray->GetEntriesFast() == 0) 
1793     return kTRUE;
1794
1795   AliRunLoader *rl = AliRunLoader::Instance();
1796   AliDataLoader *dl = 0x0;
1797   if (rl)
1798     dl = rl->GetLoader("TRDLoader")->GetDataLoader("tracklets");
1799   if (!dl) {
1800     AliError("Could not get the tracklets data loader!");
1801     return kFALSE;
1802   }
1803
1804   TTree *trackletTree = dl->Tree();
1805   if (!trackletTree) {
1806     dl->MakeTree();
1807     trackletTree = dl->Tree();
1808   }
1809   
1810   AliTRDtrackletMCM *trkl = 0x0;
1811   TBranch *trkbranch = trackletTree->GetBranch("mcmtrklbranch");
1812   if (!trkbranch)
1813     trkbranch = trackletTree->Branch("mcmtrklbranch", "AliTRDtrackletMCM", &trkl, 32000);
1814   
1815   for (Int_t iTracklet = 0; iTracklet < fTrackletArray->GetEntriesFast(); iTracklet++) {
1816     trkl = ((AliTRDtrackletMCM*) (*fTrackletArray)[iTracklet]);
1817     trkbranch->SetAddress(&trkl);
1818 //      printf("filling tracklet 0x%08x\n", trkl->GetTrackletWord());
1819     trkbranch->Fill();
1820   }
1821   dl->WriteData("OVERWRITE");
1822
1823   return kTRUE;
1824 }
1825
1826 void AliTRDmcmSim::WriteData(AliTRDarrayADC *digits)
1827 {
1828   // write back the processed data configured by EBSF
1829   // EBSF = 1: unfiltered data; EBSF = 0: filtered data
1830   // zero-suppressed valued are written as -1 to digits
1831
1832   if (!fInitialized) {
1833     AliError("Called uninitialized! Nothing done!");
1834     return;
1835   }
1836
1837 //  Int_t firstAdc = 0;
1838 //  Int_t lastAdc = fNADC - 1;
1839 //
1840 //  while (GetCol(firstAdc) < 0)
1841 //    firstAdc++;
1842 //
1843 //  while (GetCol(lastAdc) < 0) 
1844 //    lastAdc--;
1845
1846   Int_t offset = (fMcmPos % 4) * 21 + (fRobPos % 2) * 84;
1847
1848   if (fTrapConfig->GetTrapReg(AliTRDtrapConfig::kEBSF) != 0) // store unfiltered data
1849   {
1850     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
1851       if (fZSM1Dim[iAdc] == 1) {
1852         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1853           digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, -1);
1854 //          printf("suppressed: %i, %i, %i, %i, now: %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin, 
1855 //                 digits->GetData(GetRow(), GetCol(iAdc), iTimeBin));
1856         }
1857       }
1858     }
1859   }
1860   else {
1861     for (Int_t iAdc = 0; iAdc < fNADC; iAdc++) {
1862       if (fZSM1Dim[iAdc] == 0) {
1863         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1864           digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, (fADCF[iAdc][iTimeBin] >> fgkAddDigits) - fgAddBaseline);
1865         }
1866       }
1867       else {
1868         for (Int_t iTimeBin = 0; iTimeBin < fNTimeBin; iTimeBin++) {
1869           digits->SetDataByAdcCol(GetRow(), 20-iAdc + offset, iTimeBin, -1);
1870 //          printf("suppressed: %i, %i, %i, %i\n", fDetector, GetRow(), GetCol(iAdc), iTimeBin);
1871         }
1872       }
1873     }
1874   }
1875 }
1876
1877 // help functions, to be cleaned up
1878
1879 UInt_t AliTRDmcmSim::AddUintClipping(UInt_t a, UInt_t b, UInt_t nbits) const
1880 {
1881   // 
1882   // This function adds a and b (unsigned) and clips to 
1883   // the specified number of bits. 
1884   //  
1885
1886   UInt_t sum = a + b;
1887   if (nbits < 32)
1888   {
1889     UInt_t maxv = (1 << nbits) - 1;;
1890     if (sum > maxv) 
1891       sum = maxv;
1892   }
1893   else
1894   {
1895     if ((sum < a) || (sum < b)) 
1896       sum = 0xFFFFFFFF;
1897   }
1898   return sum;
1899 }
1900
1901 void AliTRDmcmSim::Sort2(UShort_t  idx1i, UShort_t  idx2i, \
1902                             UShort_t  val1i, UShort_t  val2i, \
1903                             UShort_t *idx1o, UShort_t *idx2o, \
1904                             UShort_t *val1o, UShort_t *val2o) const
1905 {
1906   // sorting for tracklet selection
1907
1908     if (val1i > val2i)
1909     {
1910         *idx1o = idx1i;
1911         *idx2o = idx2i;
1912         *val1o = val1i;
1913         *val2o = val2i;
1914     }
1915     else
1916     {
1917         *idx1o = idx2i;
1918         *idx2o = idx1i;
1919         *val1o = val2i;
1920         *val2o = val1i;
1921     }
1922 }
1923
1924 void AliTRDmcmSim::Sort3(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, \
1925                             UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, \
1926                             UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, \
1927                             UShort_t *val1o, UShort_t *val2o, UShort_t *val3o)
1928 {
1929   // sorting for tracklet selection
1930
1931     Int_t sel;
1932
1933
1934     if (val1i > val2i) sel=4; else sel=0;
1935     if (val2i > val3i) sel=sel + 2;
1936     if (val3i > val1i) sel=sel + 1;
1937     //printf("input channels %d %d %d, charges %d %d %d sel=%d\n",idx1i, idx2i, idx3i, val1i, val2i, val3i, sel);
1938     switch(sel)
1939     {
1940         case 6 : // 1 >  2  >  3            => 1 2 3
1941         case 0 : // 1 =  2  =  3            => 1 2 3 : in this case doesn't matter, but so is in hardware!
1942             *idx1o = idx1i;
1943             *idx2o = idx2i;
1944             *idx3o = idx3i;
1945             *val1o = val1i;
1946             *val2o = val2i;
1947             *val3o = val3i;
1948             break;
1949
1950         case 4 : // 1 >  2, 2 <= 3, 3 <= 1  => 1 3 2
1951             *idx1o = idx1i;
1952             *idx2o = idx3i;
1953             *idx3o = idx2i;
1954             *val1o = val1i;
1955             *val2o = val3i;
1956             *val3o = val2i;
1957             break;
1958
1959         case 2 : // 1 <= 2, 2 > 3, 3 <= 1   => 2 1 3
1960             *idx1o = idx2i;
1961             *idx2o = idx1i;
1962             *idx3o = idx3i;
1963             *val1o = val2i;
1964             *val2o = val1i;
1965             *val3o = val3i;
1966             break;
1967
1968         case 3 : // 1 <= 2, 2 > 3, 3  > 1   => 2 3 1
1969             *idx1o = idx2i;
1970             *idx2o = idx3i;
1971             *idx3o = idx1i;
1972             *val1o = val2i;
1973             *val2o = val3i;
1974             *val3o = val1i;
1975             break;
1976
1977         case 1 : // 1 <= 2, 2 <= 3, 3 > 1   => 3 2 1
1978             *idx1o = idx3i;
1979             *idx2o = idx2i;
1980             *idx3o = idx1i;
1981             *val1o = val3i;
1982             *val2o = val2i;
1983             *val3o = val1i;
1984         break;
1985
1986         case 5 : // 1 > 2, 2 <= 3, 3 >  1   => 3 1 2
1987             *idx1o = idx3i;
1988             *idx2o = idx1i;
1989             *idx3o = idx2i;
1990             *val1o = val3i;
1991             *val2o = val1i;
1992             *val3o = val2i;
1993         break;
1994
1995         default: // the rest should NEVER happen!
1996             AliError("ERROR in Sort3!!!\n");
1997         break;
1998     }
1999 //    printf("output channels %d %d %d, charges %d %d %d \n",*idx1o, *idx2o, *idx3o, *val1o, *val2o, *val3o);
2000 }
2001
2002 void AliTRDmcmSim::Sort6To4(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, UShort_t  idx4i, UShort_t  idx5i, UShort_t  idx6i, \
2003                                UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, UShort_t  val4i, UShort_t  val5i, UShort_t  val6i, \
2004                                UShort_t *idx1o, UShort_t *idx2o, UShort_t *idx3o, UShort_t *idx4o, \
2005                                UShort_t *val1o, UShort_t *val2o, UShort_t *val3o, UShort_t *val4o)
2006 {
2007   // sorting for tracklet selection
2008
2009     UShort_t idx21s, idx22s, idx23s, dummy;
2010     UShort_t val21s, val22s, val23s;
2011     UShort_t idx23as, idx23bs;
2012     UShort_t val23as, val23bs;
2013
2014     Sort3(idx1i, idx2i, idx3i, val1i, val2i, val3i,
2015                  idx1o, &idx21s, &idx23as,
2016                  val1o, &val21s, &val23as);
2017
2018     Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
2019                  idx2o, &idx22s, &idx23bs,
2020                  val2o, &val22s, &val23bs);
2021
2022     Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, &dummy, &val23s, &dummy);
2023
2024     Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
2025                  idx3o, idx4o, &dummy,
2026                  val3o, val4o, &dummy);
2027
2028 }
2029
2030 void AliTRDmcmSim::Sort6To2Worst(UShort_t  idx1i, UShort_t  idx2i, UShort_t  idx3i, UShort_t  idx4i, UShort_t  idx5i, UShort_t  idx6i, \
2031                                     UShort_t  val1i, UShort_t  val2i, UShort_t  val3i, UShort_t  val4i, UShort_t  val5i, UShort_t  val6i, \
2032                                     UShort_t *idx5o, UShort_t *idx6o)
2033 {
2034   // sorting for tracklet selection
2035
2036     UShort_t idx21s, idx22s, idx23s, dummy1, dummy2, dummy3, dummy4, dummy5;
2037     UShort_t val21s, val22s, val23s;
2038     UShort_t idx23as, idx23bs;
2039     UShort_t val23as, val23bs;
2040
2041     Sort3(idx1i, idx2i,   idx3i, val1i, val2i, val3i,
2042                  &dummy1, &idx21s, &idx23as,
2043                  &dummy2, &val21s, &val23as);
2044
2045     Sort3(idx4i, idx5i, idx6i, val4i, val5i, val6i,
2046                  &dummy1, &idx22s, &idx23bs,
2047                  &dummy2, &val22s, &val23bs);
2048
2049     Sort2(idx23as, idx23bs, val23as, val23bs, &idx23s, idx5o, &val23s, &dummy1);
2050
2051     Sort3(idx21s, idx22s, idx23s, val21s, val22s, val23s,
2052                  &dummy1, &dummy2, idx6o,
2053                  &dummy3, &dummy4, &dummy5);
2054 //    printf("idx21s=%d, idx23as=%d, idx22s=%d, idx23bs=%d, idx5o=%d, idx6o=%d\n",
2055 //            idx21s,    idx23as,    idx22s,    idx23bs,    *idx5o,    *idx6o);
2056 }
2057
2058