- adjusted error parameterization
[u/mrichter/AliRoot.git] / TRD / AliTRDcluster.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 //                                                                           //
21 //  TRD cluster                                                              //
22 //                                                                           //
23 /////////////////////////////////////////////////////////////////////////////// 
24
25 #include "TMath.h"
26
27 #include "AliLog.h"
28 #include "AliTRDcluster.h"
29 #include "AliTRDgeometry.h"
30 #include "AliTRDCommonParam.h"
31
32 ClassImp(AliTRDcluster)
33
34 const Int_t AliTRDcluster::fgkNlut = 128;
35 Double_t *AliTRDcluster::fgLUT = 0x0;
36
37 //___________________________________________________________________________
38 AliTRDcluster::AliTRDcluster() 
39   :AliCluster() 
40   ,fPadCol(0)
41   ,fPadRow(0)
42   ,fPadTime(0)
43   ,fLocalTimeBin(0)
44   ,fNPads(0)
45   ,fClusterMasking(0)
46   ,fDetector(0)
47   ,fQ(0)
48   ,fCenter(0)
49
50   //
51   // Default constructor
52   //
53
54   for (Int_t i = 0; i < 7; i++) {
55     fSignals[i] = 0;
56   }
57   SetBit(kLUT);
58 }
59
60 //___________________________________________________________________________
61 AliTRDcluster::AliTRDcluster(Int_t det, UChar_t col, UChar_t row, UChar_t time, const Short_t *sig, UShort_t vid) 
62   :AliCluster() 
63   ,fPadCol(col)
64   ,fPadRow(row)
65   ,fPadTime(time)
66   ,fLocalTimeBin(0)
67   ,fNPads(0)
68   ,fClusterMasking(0)
69   ,fDetector(det)
70   ,fQ(0.)
71   ,fCenter(0.)
72
73   //
74   // Constructor for self constructing cluster. In this approach the information is inserted gradualy into the 
75   // cluster and all dependencies are (re)calculated inside the cluster itself.
76   //
77   // A.Bercuci <A.Bercuci@gsi.de>
78
79   memcpy(&fSignals, sig, 7*sizeof(Short_t));
80   fQ = fSignals[2]+fSignals[3]+fSignals[4];
81   SetVolumeId(vid);
82   SetBit(kLUT);
83 }
84
85 //___________________________________________________________________________
86 AliTRDcluster::AliTRDcluster(Int_t det, Float_t q
87                            , Float_t *pos, Float_t *sig
88                            , Int_t *tracks, Char_t npads, Short_t *signals
89                            , UChar_t col, UChar_t row, UChar_t time
90                            , Char_t timebin, Float_t center, UShort_t volid)
91   :AliCluster(volid,pos[0],pos[1],pos[2],sig[0],sig[1],0.0,0x0) 
92   ,fPadCol(col)
93   ,fPadRow(row)
94   ,fPadTime(time)
95   ,fLocalTimeBin(timebin)
96   ,fNPads(npads)
97   ,fClusterMasking(0)
98   ,fDetector(det)
99   ,fQ(q)
100   ,fCenter(center)
101
102   //
103   // Constructor
104   //
105
106   for (Int_t i = 0; i < 7; i++) {
107     fSignals[i] = signals[i];
108   }
109
110   if (tracks) {
111     AddTrackIndex(tracks);
112   }
113   SetBit(kLUT);
114 }
115
116 //_____________________________________________________________________________
117 AliTRDcluster::AliTRDcluster(const AliTRDcluster &c)
118   :AliCluster(c)
119   ,fPadCol(c.fPadCol)
120   ,fPadRow(c.fPadRow)
121   ,fPadTime(c.fPadTime)
122   ,fLocalTimeBin(c.fLocalTimeBin)
123   ,fNPads(c.fNPads)
124   ,fClusterMasking(c.fClusterMasking)
125   ,fDetector(c.fDetector)
126   ,fQ(c.fQ)
127   ,fCenter(c.fCenter)
128 {
129   //
130   // Copy constructor 
131   //
132
133   SetLabel(c.GetLabel(0),0);
134   SetLabel(c.GetLabel(1),1);
135   SetLabel(c.GetLabel(2),2);
136
137   SetY(c.GetY());
138   SetZ(c.GetZ());
139   AliCluster::SetSigmaY2(c.GetSigmaY2());
140   SetSigmaZ2(c.GetSigmaZ2());  
141
142   for (Int_t i = 0; i < 7; i++) {
143     fSignals[i] = c.fSignals[i];
144   }
145
146 }
147
148 //_____________________________________________________________________________
149 void AliTRDcluster::AddTrackIndex(Int_t *track)
150 {
151   //
152   // Adds track index. Currently assumed that track is an array of
153   // size 9, and up to 3 track indexes are stored in fTracks[3].
154   // Indexes are sorted according to:
155   //  1) index of max number of appearances is stored first
156   //  2) if two or more indexes appear equal number of times, the lowest
157   //     ones are stored first;
158   //
159
160   const Int_t kSize = 9;
161   Int_t  entries[kSize][2];
162
163   Int_t  i = 0;
164   Int_t  j = 0;
165   Int_t  k = 0;
166   Int_t  index;
167   Bool_t indexAdded;
168
169   for (i = 0; i < kSize; i++) {
170     entries[i][0] = -1;
171     entries[i][1] =  0;
172   }                                 
173
174   for (k = 0; k < kSize; k++) {
175
176     index      = track[k];
177     indexAdded = kFALSE; 
178
179     j = 0;
180     if (index >= 0) {
181       while ((!indexAdded) && (j < kSize)) {
182         if ((entries[j][0] == index) || 
183             (entries[j][1] ==     0)) {
184           entries[j][0] = index;
185           entries[j][1] = entries[j][1] + 1;
186           indexAdded    = kTRUE;
187         }
188         j++;
189       }
190     }
191
192   }
193
194   // Sort by number of appearances and index value
195   Int_t swap = 1;
196   Int_t tmp0;
197   Int_t tmp1;
198   while (swap > 0) {
199     swap = 0;
200     for (i = 0; i < (kSize - 1); i++) {
201       if ((entries[i][0]   >= 0) && 
202           (entries[i+1][0] >= 0)) {
203         if ((entries[i][1] < entries[i+1][1]) ||
204             ((entries[i][1] == entries[i+1][1]) &&
205              (entries[i][0] >  entries[i+1][0]))) {
206           tmp0            = entries[i][0];
207           tmp1            = entries[i][1];
208           entries[i][0]   = entries[i+1][0];
209           entries[i][1]   = entries[i+1][1];
210           entries[i+1][0] = tmp0;
211           entries[i+1][1] = tmp1;
212           swap++;
213         }
214       }
215     }
216   }               
217
218   // Set track indexes
219   for (i = 0; i < 3; i++) {
220     SetLabel(entries[i][0],i);
221   }
222
223   return;
224
225 }          
226
227 //_____________________________________________________________________________
228 void AliTRDcluster::Clear(Option_t *)
229 {
230   //
231   // Reset all member to the default value
232   //
233   fPadCol=0;
234   fPadRow=0;
235   fPadTime=0;
236   fLocalTimeBin=0;
237   fNPads=0;
238   fClusterMasking=0;
239   fDetector=0;
240   for (Int_t i=0; i < 7; i++) fSignals[i]=0;
241   fQ = 0;
242   fCenter = 0;
243   for (Int_t i = 0; i < 3; i++) SetLabel(0,i);
244   SetX(0.);
245   SetY(0.);
246   SetZ(0.);
247   AliCluster::SetSigmaY2(0.);
248   SetSigmaZ2(0.);
249   SetVolumeId(0);
250 }
251
252 //_____________________________________________________________________________
253 Float_t AliTRDcluster::GetSumS() const
254 {
255   //
256   // Returns the total charge from a not unfolded cluster
257   //
258
259   Float_t sum = 0.0;
260   for (Int_t i = 0; i < 7; i++) {
261     sum += fSignals[i];
262   }
263
264   return sum;
265
266 }
267
268 //___________________________________________________________________________
269 Double_t AliTRDcluster::GetSX(Int_t tb, Double_t z)
270 {
271 // Returns the error parameterization in the radial direction for TRD clusters as function of 
272 // the calibrated time bin (tb) and optionally distance to anode wire (z). By default (no z information) 
273 // the largest value over all cluster to wire values is chosen.
274 // 
275 // The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. 
276 //
277 //Begin_Html
278 //<img src="TRD/clusterXerrorDiff2D.gif">
279 //End_Html
280 //
281 // Author
282 // A.Bercuci <A.Bercuci@gsi.de>
283
284   if(tb<1 || tb>=24) return 10.; // return huge [10cm]
285   const Double_t sx[24][10]={
286     {0.000e+00, 9.352e-01, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 2.309e+00},
287     {8.387e-02, 8.718e-02, 8.816e-02, 9.444e-02, 9.993e-02, 1.083e-01, 1.161e-01, 1.280e-01, 1.417e-01, 1.406e-01},
288     {1.097e-01, 1.105e-01, 1.127e-01, 1.151e-01, 1.186e-01, 1.223e-01, 1.272e-01, 1.323e-01, 1.389e-01, 1.490e-01},
289     {1.407e-01, 1.404e-01, 1.414e-01, 1.430e-01, 1.429e-01, 1.449e-01, 1.476e-01, 1.494e-01, 1.515e-01, 1.589e-01},
290     {1.681e-01, 1.679e-01, 1.666e-01, 1.657e-01, 1.656e-01, 1.649e-01, 1.652e-01, 1.662e-01, 1.671e-01, 1.694e-01},
291     {1.745e-01, 1.737e-01, 1.707e-01, 1.690e-01, 1.643e-01, 1.610e-01, 1.612e-01, 1.628e-01, 1.638e-01, 1.659e-01},
292     {1.583e-01, 1.558e-01, 1.535e-01, 1.488e-01, 1.445e-01, 1.419e-01, 1.428e-01, 1.451e-01, 1.462e-01, 1.494e-01},
293     {1.414e-01, 1.391e-01, 1.368e-01, 1.300e-01, 1.256e-01, 1.259e-01, 1.285e-01, 1.326e-01, 1.358e-01, 1.406e-01},
294     {1.307e-01, 1.289e-01, 1.261e-01, 1.216e-01, 1.193e-01, 1.165e-01, 1.201e-01, 1.241e-01, 1.274e-01, 1.344e-01},
295     {1.251e-01, 1.227e-01, 1.208e-01, 1.155e-01, 1.110e-01, 1.116e-01, 1.133e-01, 1.187e-01, 1.229e-01, 1.308e-01},
296     {1.234e-01, 1.209e-01, 1.175e-01, 1.127e-01, 1.094e-01, 1.093e-01, 1.109e-01, 1.155e-01, 1.210e-01, 1.275e-01},
297     {1.215e-01, 1.187e-01, 1.156e-01, 1.108e-01, 1.070e-01, 1.065e-01, 1.090e-01, 1.134e-01, 1.196e-01, 1.251e-01},
298     {1.202e-01, 1.180e-01, 1.151e-01, 1.108e-01, 1.070e-01, 1.058e-01, 1.089e-01, 1.127e-01, 1.183e-01, 1.256e-01},
299     {1.207e-01, 1.176e-01, 1.142e-01, 1.109e-01, 1.072e-01, 1.069e-01, 1.088e-01, 1.122e-01, 1.182e-01, 1.252e-01},
300     {1.213e-01, 1.182e-01, 1.156e-01, 1.102e-01, 1.076e-01, 1.063e-01, 1.091e-01, 1.132e-01, 1.181e-01, 1.243e-01},
301     {1.205e-01, 1.180e-01, 1.150e-01, 1.104e-01, 1.072e-01, 1.063e-01, 1.083e-01, 1.132e-01, 1.183e-01, 1.243e-01},
302     {1.212e-01, 1.195e-01, 1.135e-01, 1.107e-01, 1.070e-01, 1.065e-01, 1.097e-01, 1.126e-01, 1.185e-01, 1.238e-01},
303     {1.201e-01, 1.184e-01, 1.155e-01, 1.111e-01, 1.088e-01, 1.075e-01, 1.089e-01, 1.131e-01, 1.189e-01, 1.237e-01},
304     {1.197e-01, 1.186e-01, 1.147e-01, 1.113e-01, 1.085e-01, 1.077e-01, 1.105e-01, 1.137e-01, 1.188e-01, 1.245e-01},
305     {1.213e-01, 1.194e-01, 1.154e-01, 1.114e-01, 1.091e-01, 1.082e-01, 1.098e-01, 1.140e-01, 1.194e-01, 1.247e-01},
306     {1.210e-01, 1.189e-01, 1.155e-01, 1.119e-01, 1.088e-01, 1.080e-01, 1.105e-01, 1.141e-01, 1.195e-01, 1.244e-01},
307     {1.196e-01, 1.189e-01, 1.145e-01, 1.105e-01, 1.095e-01, 1.083e-01, 1.087e-01, 1.121e-01, 1.173e-01, 1.208e-01},
308     {1.123e-01, 1.129e-01, 1.108e-01, 1.110e-01, 1.080e-01, 1.065e-01, 1.056e-01, 1.066e-01, 1.071e-01, 1.095e-01},
309     {1.136e-01, 1.135e-01, 1.130e-01, 1.122e-01, 1.113e-01, 1.071e-01, 1.041e-01, 1.025e-01, 1.014e-01, 9.973e-02}
310   };
311   if(z>=0. && z<.25) return sx[tb][Int_t(z/.025)];
312   
313   Double_t m = 0.; for(Int_t id=10; id--;) m+=sx[tb][id];
314   return m*.1;
315 }
316
317 //___________________________________________________________________________
318 Double_t AliTRDcluster::GetSYdrift(Int_t tb, Int_t ly, Double_t/* z*/)
319 {
320 // Returns the error parameterization for TRD clusters as function of the drift length (here calibrated time bin tb)
321 // and optionally distance to anode wire (z) for the LUT r-phi cluster shape. By default (no z information) the largest 
322 // value over all cluster to wire values is chosen.
323 // 
324 // For the LUT method the dependence of s_y with x and d is obtained via a fit to the cluster to MC 
325 // resolution. (see class AliTRDclusterResolution for more details). A normalization to the reference radial position
326 // x0 = 0.675 (tb=5 for ideal vd) is also applied (see GetSYprf()). The function is *NOT* calibration aware ! 
327 // The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. A comparison 
328 // with the GAUS parameterization is also given
329 //
330 // For the GAUS method the dependence of s_y with x is *analytic* and it is expressed by the relation.
331 // BEGIN_LATEX
332 // #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}}
333 // END_LATEX
334 // The result is displayed in the figure below as function of the drift time and compared with the LUT parameterization.
335 //Begin_Html
336 //<img src="TRD/clusterYerrorDiff2D.gif">
337 //<img src="TRD/clusterYerrorDiff1D.gif">
338 //End_Html
339 //
340 // Author
341 // A.Bercuci <A.Bercuci@gsi.de>
342
343
344   if(tb<1 || tb>=24) return 10.; // return huge [10cm]
345   const Float_t lSy[6][24] = {
346     {75.7561, 0.0325, 0.0175, 0.0174, 0.0206, 0.0232,
347      0.0253, 0.0262, 0.0265, 0.0264, 0.0266, 0.0257,
348      0.0258, 0.0261, 0.0259, 0.0253, 0.0257, 0.0261,
349      0.0255, 0.0250, 0.0259, 0.0266, 0.0278, 0.0319
350     },
351     {49.2252, 0.0371, 0.0204, 0.0189, 0.0230, 0.0261,
352      0.0281, 0.0290, 0.0292, 0.0286, 0.0277, 0.0279,
353      0.0285, 0.0281, 0.0291, 0.0281, 0.0281, 0.0282,
354      0.0272, 0.0282, 0.0282, 0.0284, 0.0310, 0.0334
355     },
356     {55.1674, 0.0388, 0.0212, 0.0200, 0.0239, 0.0271,
357      0.0288, 0.0299, 0.0306, 0.0300, 0.0296, 0.0303,
358      0.0293, 0.0290, 0.0291, 0.0294, 0.0295, 0.0290,
359      0.0293, 0.0292, 0.0292, 0.0293, 0.0316, 0.0358
360     },
361     {45.1004, 0.0411, 0.0225, 0.0215, 0.0249, 0.0281,
362      0.0301, 0.0315, 0.0320, 0.0308, 0.0318, 0.0321,
363      0.0312, 0.0311, 0.0316, 0.0315, 0.0310, 0.0308,
364      0.0313, 0.0303, 0.0314, 0.0314, 0.0324, 0.0369
365     },
366     {43.8614, 0.0420, 0.0239, 0.0224, 0.0268, 0.0296,
367      0.0322, 0.0336, 0.0333, 0.0326, 0.0321, 0.0325,
368      0.0329, 0.0326, 0.0323, 0.0322, 0.0326, 0.0320,
369      0.0329, 0.0319, 0.0314, 0.0329, 0.0341, 0.0373
370     },
371     {40.5440, 0.0434, 0.0246, 0.0236, 0.0275, 0.0311,
372      0.0332, 0.0345, 0.0347, 0.0347, 0.0340, 0.0336,
373      0.0339, 0.0344, 0.0339, 0.0341, 0.0341, 0.0342,
374      0.0345, 0.0328, 0.0341, 0.0332, 0.0356, 0.0398
375     },
376   };
377   // adjusted ...
378   return TMath::Max(lSy[ly][tb]-0.0150, 0.0010);
379
380 /*  const Double_t sy[24][10]={
381     {0.000e+00, 2.610e-01, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 4.680e-01},
382     {3.019e-02, 3.036e-02, 3.131e-02, 3.203e-02, 3.294e-02, 3.407e-02, 3.555e-02, 3.682e-02, 3.766e-02, 3.824e-02},
383     {1.773e-02, 1.778e-02, 1.772e-02, 1.790e-02, 1.807e-02, 1.833e-02, 1.873e-02, 1.905e-02, 1.958e-02, 2.029e-02},
384     {1.774e-02, 1.772e-02, 1.746e-02, 1.738e-02, 1.756e-02, 1.756e-02, 1.739e-02, 1.720e-02, 1.743e-02, 1.769e-02},
385     {2.064e-02, 2.078e-02, 2.069e-02, 2.060e-02, 2.033e-02, 2.024e-02, 2.022e-02, 1.961e-02, 1.922e-02, 1.901e-02},
386     {2.382e-02, 2.379e-02, 2.371e-02, 2.333e-02, 2.318e-02, 2.285e-02, 2.255e-02, 2.244e-02, 2.174e-02, 2.132e-02},
387     {2.615e-02, 2.589e-02, 2.539e-02, 2.493e-02, 2.420e-02, 2.396e-02, 2.362e-02, 2.342e-02, 2.321e-02, 2.330e-02},
388     {2.640e-02, 2.638e-02, 2.577e-02, 2.548e-02, 2.477e-02, 2.436e-02, 2.416e-02, 2.401e-02, 2.399e-02, 2.402e-02},
389     {2.647e-02, 2.632e-02, 2.587e-02, 2.546e-02, 2.465e-02, 2.447e-02, 2.429e-02, 2.415e-02, 2.429e-02, 2.475e-02},
390     {2.657e-02, 2.637e-02, 2.580e-02, 2.525e-02, 2.492e-02, 2.441e-02, 2.446e-02, 2.441e-02, 2.478e-02, 2.491e-02},
391     {2.640e-02, 2.608e-02, 2.583e-02, 2.539e-02, 2.478e-02, 2.440e-02, 2.456e-02, 2.464e-02, 2.486e-02, 2.533e-02},
392     {2.636e-02, 2.630e-02, 2.584e-02, 2.542e-02, 2.483e-02, 2.451e-02, 2.449e-02, 2.467e-02, 2.496e-02, 2.554e-02},
393     {2.634e-02, 2.629e-02, 2.583e-02, 2.526e-02, 2.480e-02, 2.460e-02, 2.458e-02, 2.472e-02, 2.518e-02, 2.549e-02},
394     {2.629e-02, 2.621e-02, 2.581e-02, 2.527e-02, 2.480e-02, 2.458e-02, 2.451e-02, 2.485e-02, 2.516e-02, 2.547e-02},
395     {2.629e-02, 2.607e-02, 2.573e-02, 2.543e-02, 2.485e-02, 2.464e-02, 2.452e-02, 2.476e-02, 2.505e-02, 2.550e-02},
396     {2.635e-02, 2.613e-02, 2.578e-02, 2.523e-02, 2.491e-02, 2.465e-02, 2.470e-02, 2.467e-02, 2.515e-02, 2.564e-02},
397     {2.613e-02, 2.602e-02, 2.587e-02, 2.526e-02, 2.507e-02, 2.482e-02, 2.456e-02, 2.486e-02, 2.509e-02, 2.572e-02},
398     {2.620e-02, 2.599e-02, 2.563e-02, 2.528e-02, 2.484e-02, 2.462e-02, 2.464e-02, 2.476e-02, 2.513e-02, 2.571e-02},
399     {2.634e-02, 2.596e-02, 2.565e-02, 2.519e-02, 2.497e-02, 2.457e-02, 2.450e-02, 2.481e-02, 2.511e-02, 2.540e-02},
400     {2.593e-02, 2.589e-02, 2.563e-02, 2.511e-02, 2.472e-02, 2.453e-02, 2.452e-02, 2.474e-02, 2.501e-02, 2.543e-02},
401     {2.576e-02, 2.582e-02, 2.526e-02, 2.505e-02, 2.462e-02, 2.446e-02, 2.445e-02, 2.466e-02, 2.486e-02, 2.510e-02},
402     {2.571e-02, 2.549e-02, 2.533e-02, 2.501e-02, 2.453e-02, 2.443e-02, 2.445e-02, 2.450e-02, 2.448e-02, 2.469e-02},
403     {2.812e-02, 2.786e-02, 2.776e-02, 2.723e-02, 2.695e-02, 2.650e-02, 2.642e-02, 2.617e-02, 2.612e-02, 2.610e-02},
404     {3.251e-02, 3.267e-02, 3.223e-02, 3.183e-02, 3.125e-02, 3.106e-02, 3.067e-02, 3.010e-02, 2.936e-02, 2.927e-02}
405   };
406   if(z>=0. && z<.25) return sy[tb][Int_t(z/.025)] -  sy[5][Int_t(z/.025)];
407
408   Double_t m = -1.e8; for(Int_t id=10; id--;) if((sy[tb][id] - sy[5][id])>m) m=sy[tb][id]-sy[5][id];
409
410   return m;*/
411 }
412
413 //___________________________________________________________________________
414 Double_t AliTRDcluster::GetSYcharge(Float_t q)
415 {
416 // Parameterization of the r-phi resolution component due to cluster charge.
417 // The value is the offset from the nominal cluster resolution defined as the
418 // cluster resolution at average cluster charge (q0).
419 // 
420 // BEGIN_LATEX
421 // #Delta #sigma_{y}(q) = a*(#frac{1}{q} - #frac{1}{q_{0}})
422 // q_{0} #approx 50
423 // END_LATEX
424 // The definition is *NOT* robust against gain fluctuations and thus two approaches are possible
425 // when residual miscalibration are available:
426 //   - determine parameterization with a resolution matching those of the gain
427 //   - define an analytic model which scales with the gain.
428 // 
429 // For more details please see AliTRDclusterResolution::ProcessCharge()
430 //
431 //Begin_Html
432 //<img src="TRD/clusterQerror.gif">
433 //End_Html
434 //
435 // Author
436 // A.Bercuci <A.Bercuci@gsi.de>
437
438   const Float_t sq0inv = 0.019962; // [1/q0]
439   const Float_t sqb    = 0.037328; // [cm]
440
441   return sqb*(1./q - sq0inv);
442 }
443
444 //___________________________________________________________________________
445 Double_t AliTRDcluster::GetSYprf(Int_t ly, Double_t center, Double_t s2)
446 {
447 // Parameterization of the cluster error in the r-phi direction due to charge sharing between 
448 // adiacent pads. Should be identical to what is provided in the OCDB as PRF [TODO]
449 //
450 // The parameterization is obtained from fitting cluster resolution at phi=exb and |x-0.675|<0.225. 
451 // For more details see AliTRDclusterResolution::ProcessCenter().
452 //
453 //Begin_Html
454 //<img src="TRD/clusterPRFerror.gif">
455 //End_Html
456 //
457 // Author
458 // A.Bercuci <A.Bercuci@gsi.de>
459
460 /*  const Float_t scy[AliTRDgeometry::kNlayer][4] = {
461     {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
462     {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
463     {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
464     {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
465     {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
466     {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
467   };*/
468   const Float_t lPRF[] = {0.438, 0.403, 0.393, 0.382, 0.376, 0.345};
469
470   return s2*TMath::Gaus(center, 0., lPRF[ly]);
471 }
472
473
474 //___________________________________________________________________________
475 Double_t AliTRDcluster::GetXcorr(Int_t tb, Double_t z)
476 {
477 // Drift length correction [cm]. Due to variation of mean drift velocity along the drift region
478 // from nominal vd at xd->infinity. For drift velocity determination based on tracking information 
479 // the correction should be negligible.
480 //
481 // TODO to be parametrized in term of drift velocity at infinite drift length
482 // A.Bercuci (Mar 28 2009)
483
484   if(tb<0 || tb>=24) return 0.;
485   const Int_t nd = 5;
486   const Double_t dx[24][nd]={
487     {+1.747e-01,+3.195e-01,+1.641e-01,+1.607e-01,+6.002e-01},
488     {+5.468e-02,+5.760e-02,+6.365e-02,+8.003e-02,+1.067e-01},
489     {-6.327e-02,-6.339e-02,-6.423e-02,-6.900e-02,-7.949e-02},
490     {-1.417e-01,-1.424e-01,-1.450e-01,-1.465e-01,-1.514e-01},
491     {-1.637e-01,-1.619e-01,-1.622e-01,-1.613e-01,-1.648e-01},
492     {-1.386e-01,-1.334e-01,-1.261e-01,-1.276e-01,-1.314e-01},
493     {-8.799e-02,-8.299e-02,-7.861e-02,-8.038e-02,-8.436e-02},
494     {-5.139e-02,-4.849e-02,-4.641e-02,-4.965e-02,-5.286e-02},
495     {-2.927e-02,-2.773e-02,-2.807e-02,-3.021e-02,-3.378e-02},
496     {-1.380e-02,-1.229e-02,-1.335e-02,-1.547e-02,-1.984e-02},
497     {-4.168e-03,-4.601e-03,-5.462e-03,-8.164e-03,-1.035e-02},
498     {+2.044e-03,+1.889e-03,+9.603e-04,-1.342e-03,-3.736e-03},
499     {+3.568e-03,+3.581e-03,+2.391e-03,+2.942e-05,-1.585e-03},
500     {+4.403e-03,+4.571e-03,+3.509e-03,+8.703e-04,-1.425e-03},
501     {+4.941e-03,+4.808e-03,+3.284e-03,+1.105e-03,-1.208e-03},
502     {+5.124e-03,+5.022e-03,+4.305e-03,+2.023e-03,-1.145e-03},
503     {+4.882e-03,+4.008e-03,+3.408e-03,+7.886e-04,-1.356e-03},
504     {+3.852e-03,+3.539e-03,+2.057e-03,+1.670e-04,-1.993e-03},
505     {+2.154e-03,+2.111e-03,+5.723e-04,-1.254e-03,-3.256e-03},
506     {+1.755e-03,+2.101e-03,+9.516e-04,-1.649e-03,-3.394e-03},
507     {+1.617e-03,+1.662e-03,+4.169e-04,-9.843e-04,-4.309e-03},
508     {-9.204e-03,-9.069e-03,-1.182e-02,-1.458e-02,-1.880e-02},
509     {-6.727e-02,-6.820e-02,-6.804e-02,-7.134e-02,-7.615e-02},
510     {-1.802e-01,-1.733e-01,-1.633e-01,-1.601e-01,-1.632e-01}
511   };
512 //   const Double_t dx[24][nd]={
513 //     {+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00,+0.000e+00},
514 //     {-2.763e-04,-2.380e-04,-6.286e-04,-9.424e-04,+1.046e-03,+1.932e-03,+1.620e-03,+1.951e-03,-1.321e-03,-1.115e-03},
515 //     {-1.825e-03,-9.245e-04,-1.012e-03,-8.215e-04,+2.703e-05,+1.403e-03,+2.340e-03,+2.577e-03,+2.017e-03,+8.006e-04},
516 //     {-3.070e-03,-8.563e-04,-1.257e-03,+8.491e-05,+4.503e-04,-2.467e-05,-1.793e-04,+5.085e-04,+1.321e-03,+4.056e-04},
517 //     {-3.637e-03,-2.857e-03,-3.098e-03,-2.304e-03,-1.467e-03,-1.755e-03,+4.585e-04,+2.757e-03,+3.184e-03,+3.525e-03},
518 //     {-9.884e-03,-7.695e-03,-7.290e-03,-3.990e-03,-9.982e-04,+2.226e-03,+3.375e-03,+6.419e-03,+7.209e-03,+6.891e-03},
519 //     {-6.844e-03,-5.807e-03,-4.012e-03,-1.566e-03,+5.035e-04,+2.024e-03,+3.225e-03,+3.918e-03,+5.942e-03,+6.024e-03},
520 //     {-2.628e-03,-2.201e-03,-4.562e-04,+9.832e-04,+3.411e-03,+2.062e-03,+1.526e-03,+9.350e-04,+8.842e-04,+1.007e-03},
521 //     {+6.603e-04,+1.545e-03,+1.681e-03,+1.918e-03,+2.165e-03,+1.825e-03,+1.691e-03,-1.923e-04,+1.835e-04,-1.284e-03},
522 //     {+1.895e-03,+1.586e-03,+2.000e-03,+3.537e-03,+2.526e-03,+1.316e-03,+8.229e-04,-7.671e-05,-2.175e-03,-3.528e-03},
523 //     {+2.927e-03,+3.369e-03,+3.603e-03,+2.675e-03,+2.737e-03,+1.133e-03,+4.318e-04,-1.215e-03,-2.443e-03,-3.116e-03},
524 //     {+3.072e-03,+3.564e-03,+3.612e-03,+3.149e-03,+2.768e-03,+1.186e-03,+3.083e-04,-1.447e-03,-2.480e-03,-3.263e-03},
525 //     {+2.697e-03,+3.565e-03,+3.759e-03,+2.855e-03,+2.909e-03,+6.564e-04,-5.224e-04,-3.309e-04,-1.636e-03,-3.739e-03},
526 //     {+3.633e-03,+3.232e-03,+3.727e-03,+3.024e-03,+3.365e-03,+1.598e-03,-6.903e-04,-1.039e-03,-3.176e-03,-4.472e-03},
527 //     {+2.999e-03,+3.942e-03,+3.322e-03,+3.162e-03,+1.978e-03,+1.657e-03,-4.760e-04,-8.343e-04,-2.346e-03,-3.281e-03},
528 //     {+3.734e-03,+3.098e-03,+3.435e-03,+2.512e-03,+2.651e-03,+1.745e-03,+9.424e-04,-1.404e-03,-3.177e-03,-4.444e-03},
529 //     {+3.204e-03,+4.003e-03,+3.068e-03,+2.697e-03,+3.187e-03,+3.878e-04,-1.124e-04,-1.855e-03,-2.584e-03,-3.807e-03},
530 //     {+2.653e-03,+3.631e-03,+2.327e-03,+3.460e-03,+1.810e-03,+1.244e-03,-3.651e-04,-2.664e-04,-2.307e-03,-3.642e-03},
531 //     {+2.538e-03,+3.208e-03,+2.390e-03,+3.519e-03,+1.763e-03,+1.330e-04,+1.669e-04,-1.422e-03,-1.685e-03,-3.519e-03},
532 //     {+2.605e-03,+2.465e-03,+2.771e-03,+2.966e-03,+2.361e-03,+6.029e-04,-4.435e-04,-1.876e-03,-1.694e-03,-3.757e-03},
533 //     {+2.866e-03,+3.315e-03,+3.146e-03,+2.117e-03,+1.933e-03,+9.339e-04,+9.556e-04,-1.314e-03,-3.615e-03,-3.558e-03},
534 //     {+4.002e-03,+3.543e-03,+3.631e-03,+4.127e-03,+1.919e-03,-2.852e-04,-9.484e-04,-2.060e-03,-4.477e-03,-5.491e-03},
535 //     {+6.029e-03,+5.147e-03,+4.286e-03,+2.215e-03,+9.240e-04,-1.554e-03,-2.366e-03,-3.635e-03,-5.372e-03,-6.467e-03},
536 //     {+3.941e-03,+3.995e-03,+5.638e-04,-3.332e-04,-2.539e-03,-3.764e-03,-3.647e-03,-4.900e-03,-5.414e-03,-5.202e-03}
537 //   };
538   if(z>=0. && z<.25) return dx[tb][Int_t(z/.025)];
539
540   Double_t m = 0.; for(Int_t id=nd; id--;) m+=dx[tb][id];
541   return m/nd;
542 }
543
544 //___________________________________________________________________________
545 Double_t AliTRDcluster::GetYcorr(Int_t ly, Float_t y)
546 {
547 // PRF correction for the LUT r-phi cluster shape.
548
549   const Float_t cy[AliTRDgeometry::kNlayer][3] = {
550     { 4.014e-04, 8.605e-03, -6.880e+00},
551     {-3.061e-04, 9.663e-03, -6.789e+00},
552     { 1.124e-03, 1.105e-02, -6.825e+00},
553     {-1.527e-03, 1.231e-02, -6.777e+00},
554     { 2.150e-03, 1.387e-02, -6.783e+00},
555     {-1.296e-03, 1.486e-02, -6.825e+00}
556   }; 
557
558   return cy[ly][0] + cy[ly][1] * TMath::Sin(cy[ly][2] * y);
559 }
560
561 //_____________________________________________________________________________
562 Float_t AliTRDcluster::GetXloc(Double_t t0, Double_t vd, Double_t *const /*q*/, Double_t *const /*xq*/, Double_t /*z*/)
563 {
564 //
565 // (Re)Calculate cluster position in the x direction in local chamber coordinates (with respect to the anode wire 
566 // position) using all available information from tracking.
567 // Input parameters:
568 //   t0 - calibration aware trigger delay [us]
569 //   vd - drift velocity in the region of the cluster [cm/us]
570 //   z  - distance to the anode wire [cm]. By default 0.2 !!
571 //   q & xq - array of charges and cluster positions from previous clusters in the tracklet [a.u.]
572 // Output values :
573 //   return x position of the cluster with respect to the 
574 //   anode wire using all tracking information
575 //
576 // The estimation of the radial position is based on calculating the drift time and the drift velocity at the point of 
577 // estimation. The drift time can be estimated according to the expression:
578 // BEGIN_LATEX
579 // t_{drift} = t_{bin} - t_{0} - t_{cause}(x) - t_{TC}(q_{i-1}, q_{i-2}, ...)
580 // END_LATEX
581 // where t_0 is the delay of the trigger signal. t_cause is the causality delay between ionisation electrons hitting 
582 // the anode and the registration of maximum signal by the electronics - it is due to the rising time of the TRF 
583 // convoluted with the diffusion width. t_TC is the residual charge from previous bins due to residual tails after tail 
584 // cancellation.
585 //
586 // The drift velocity is considered to vary linearly with the drift length (independent of the distance to the anode wire 
587 // in the z direction). Thus one can write the calculate iteratively the drift length from the expression:
588 // BEGIN_LATEX
589 // x = t_{drift}(x)*v_{drfit}(x)
590 // END_LATEX
591 //
592 // Authors
593 // Alex Bercuci <A.Bercuci@gsi.de>
594 //
595
596   AliTRDCommonParam *cp = AliTRDCommonParam::Instance(); 
597   Double_t fFreq = cp->GetSamplingFrequency();
598
599   //drift time corresponding to the center of the time bin
600   Double_t td = (fPadTime + .5)/fFreq; // [us] 
601   // correction for t0
602   td -= t0;
603   // time bin corrected for t0
604   // Bug in TMath::Nint().root-5.23.02
605   // TMath::Nint(3.5) = 4 and TMath::Nint(4.5) = 4
606   Double_t tmp = td*fFreq;
607   fLocalTimeBin = Char_t(TMath::Floor(tmp));
608   if(tmp-fLocalTimeBin > .5) fLocalTimeBin++;
609   if(td < .2) return 0.;
610   // TRF rising time (fitted)
611   // It should be absorbed by the t0. For the moment t0 is 0 for simulations.
612   // A.Bercuci (Mar 26 2009)
613   td -= 0.189;
614
615   // apply fitted correction 
616   Float_t x = td*vd + GetXcorr(fLocalTimeBin);
617   if(x>0.&&x<.5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()) SetInChamber();
618
619   return x;
620
621 /*
622   // calculate radial posion of clusters in the drift region
623
624   // invert drift time function
625   Double_t xM= AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght(),
626            x = vd*td + .5*AliTRDgeometry::CamHght(), 
627            t = cp->TimeStruct(vd, x, z), dx1=0.,dx2;
628   while(TMath::Abs(td-t)>1.e-4){ // convergence on 100ps
629     dx2 = vd*(td-t);
630     if(TMath::Abs(TMath::Abs(dx2)-TMath::Abs(dx1))<1.e-6){
631       x+=.5*dx2;
632       break;
633     } else x+=dx2;
634
635     if(x<0. || x>xM) return 0.;
636     t = cp->TimeStruct(vd, x, z);
637     dx1 = dx2;
638   }
639
640   return x-.5*AliTRDgeometry::CamHght();
641 */
642 }
643
644 //_____________________________________________________________________________
645 Float_t AliTRDcluster::GetYloc(Double_t y0, Double_t s2, Double_t W, Double_t *const y1, Double_t *const y2)
646 {
647
648   //printf("  s[%3d %3d %3d] w[%f %f] yr[%f %f]\n", fSignals[2], fSignals[3], fSignals[4], w1/(w1+w2), w2/(w1+w2), y1r*W, y2r*W);
649   if(IsRPhiMethod(kCOG)) GetDYcog();
650   else if(IsRPhiMethod(kLUT)) GetDYlut();
651   else if(IsRPhiMethod(kGAUS)) GetDYgauss(s2/W/W, y1, y2);
652   else return 0.;
653
654   if(y1) (*y1)*=W;
655   if(y2) (*y2)*=W;
656
657   return y0+fCenter*W+(IsRPhiMethod(kLUT)?GetYcorr(AliTRDgeometry::GetLayer(fDetector), fCenter):0.);
658 }
659
660 //___________________________________________________________________________
661 void AliTRDcluster::SetSigmaY2(Float_t s2, Float_t dt, Float_t exb, Float_t x, Float_t z, Float_t tgp)
662 {
663 // Set variance of TRD cluster in the r-phi direction for each method.
664 // Parameters :
665 //   - s2  - variance due to PRF width for the case of Gauss model. Replaced by parameterization in case of LUT.
666 //   - dt  - transversal diffusion coeficient 
667 //   - exb - tg of lorentz angle
668 //   - x   - drift length - with respect to the anode wire
669 //   - z   - offset from the anode wire
670 //   - tgp - local tangent of the track momentum azimuthal angle
671 //
672 // The ingredients from which the error is computed are:
673 //   - PRF (charge sharing on adjacent pads)  - see AliTRDcluster::GetSYprf()
674 //   - diffusion (dependence with drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSYdrift()
675 //   - charge of the cluster (complex dependence on gain and tail cancellation) - see AliTRDcluster::GetSYcharge()
676 //   - lorentz angle (dependence on the drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSX()
677 //   - track angle (superposition of charges on the anode wire) - see AliTRDseedV1::Fit()
678 //   - projection of radial(x) error on r-phi due to fixed value assumed in tracking for x - see AliTRDseedV1::Fit()
679 //
680 // The last 2 contributions to cluster error can be estimated only during tracking when the track angle 
681 // is known (tgp). For this reason the errors (and optional position) of TRD clusters are recalculated during 
682 // tracking and thus clusters attached to tracks might differ from bare clusters.
683 // 
684 // Taking into account all contributions one can write the the TRD cluster error parameterization as:
685 // BEGIN_LATEX
686 // #sigma_{y}^{2} = (#sigma_{diff}*Gauss(0, s_{ly}) + #delta_{#sigma}(q))^{2} + tg^{2}(#alpha_{L})*#sigma_{x}^{2} + tg^{2}(#phi-#alpha_{L})*#sigma_{x}^{2}+[tg(#phi-#alpha_{L})*tg(#alpha_{L})*x]^{2}/12
687 // END_LATEX
688 // From this formula one can deduce a that the simplest calibration method for PRF and diffusion contributions is 
689 // by measuring resolution at B=0T and phi=0. To disentangle further the two remaining contributions one has 
690 // to represent s2 as a function of drift length. 
691 // 
692 // In the gaussian model the diffusion contribution can be expressed as:
693 // BEGIN_LATEX
694 // #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}}
695 // END_LATEX
696 // thus resulting the PRF contribution. For the case of the LUT model both contributions have to be determined from 
697 // the fit (see AliTRDclusterResolution::ProcessCenter() for details).
698 // 
699 // Author:
700 // A.Bercuci <A.Bercuci@gsi.de>
701
702   Float_t sigmaY2 = 0.;
703   Int_t ly = AliTRDgeometry::GetLayer(fDetector);
704   if(IsRPhiMethod(kCOG)) sigmaY2 = 4.e-4;
705   else if(IsRPhiMethod(kLUT)){ 
706     Float_t sd = GetSYdrift(fLocalTimeBin, ly, z); //printf("drift[%6.2f] ", 1.e4*sd);
707     sigmaY2 = GetSYprf(ly, fCenter, sd); //printf("PRF[%6.2f] ", 1.e4*sigmaY2);
708     // add charge contribution TODO scale with respect to s2
709     sigmaY2+= GetSYcharge(TMath::Abs(fQ)); //printf("Q[%6.2f] ", 1.e4*sigmaY2);
710     sigmaY2 = TMath::Max(sigmaY2, Float_t(0.0010)); //!! protection 
711     sigmaY2*= sigmaY2;
712   } else if(IsRPhiMethod(kGAUS)){
713     // PRF contribution
714     sigmaY2 = s2;
715     // Diffusion contribution
716     Double_t sD2 = dt/(1.+exb); sD2 *= sD2; sD2 *= x;
717     sigmaY2+= sD2; 
718     // add charge contribution TODO scale with respect to s2
719     //sigmaY2+= GetSYcharge(TMath::Abs(fQ));
720   }
721
722   // store tg^2(phi-a_L) and tg^2(a_L)
723   Double_t tgg = (tgp-exb)/(1.+tgp*exb); tgg *= tgg;
724   Double_t exb2= exb*exb;
725
726   // Lorentz angle shift contribution 
727   Float_t sx = GetSX(fLocalTimeBin, z); sx*=sx;
728   sigmaY2+= exb2*sx; //printf("Al[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2));
729
730   // Radial contribution due to not measuring x in Kalman model 
731   sigmaY2+= tgg*sx; //printf("x[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2));
732
733   // Track angle contribution
734   sigmaY2+= tgg*x*x*exb2/12.; //printf("angle[%6.2f]\n", 1.e4*TMath::Sqrt(sigmaY2));
735
736   AliCluster::SetSigmaY2(sigmaY2);
737 }
738
739 //_____________________________________________________________________________
740 Bool_t AliTRDcluster::IsEqual(const TObject *o) const
741 {
742   //
743   // Compare relevant information of this cluster with another one
744   //
745   
746   const AliTRDcluster *inCluster = dynamic_cast<const AliTRDcluster*>(o);
747   if (!o || !inCluster) return kFALSE;
748
749   if ( AliCluster::GetX() != inCluster->GetX() ) return kFALSE;
750   if ( AliCluster::GetY() != inCluster->GetY() ) return kFALSE;
751   if ( AliCluster::GetZ() != inCluster->GetZ() ) return kFALSE;
752   if ( fQ != inCluster->fQ ) return kFALSE;
753   if ( fDetector != inCluster->fDetector ) return kFALSE;
754   if ( fPadCol != inCluster->fPadCol ) return kFALSE;
755   if ( fPadRow != inCluster->fPadRow ) return kFALSE;
756   if ( fPadTime != inCluster->fPadTime ) return kFALSE;
757   if ( fClusterMasking != inCluster->fClusterMasking ) return kFALSE;
758   if ( IsInChamber() != inCluster->IsInChamber() ) return kFALSE;
759   if ( IsShared() != inCluster->IsShared() ) return kFALSE;
760   if ( IsUsed() != inCluster->IsUsed() ) return kFALSE;
761   
762   return kTRUE;
763 }
764
765 //_____________________________________________________________________________
766 void AliTRDcluster::Print(Option_t *o) const
767 {
768   AliInfo(Form("Det[%3d] LTrC[%+6.2f %+6.2f %+6.2f] Q[%5.1f] FLAG[in(%c) use(%c) sh(%c)] Y[%s]", 
769     fDetector, GetX(), GetY(), GetZ(), fQ, 
770     IsInChamber() ? 'y' : 'n', 
771     IsUsed() ? 'y' : 'n', 
772     IsShared() ? 'y' : 'n',
773     IsRPhiMethod(kGAUS)?"GAUS":(IsRPhiMethod(kLUT)?"LUT":"COG")
774   ));
775
776   if(strcmp(o, "a")!=0) return;
777   AliInfo(Form("LChC[c(%3d) r(%2d) t(%2d)] t-t0[%2d] Npad[%d] cen[%5.3f] mask[%d]", fPadCol, fPadRow, fPadTime, fLocalTimeBin, fNPads, fCenter, fClusterMasking)); 
778   AliInfo(Form("Signals[%3d %3d %3d %3d %3d %3d %3d]", fSignals[0], fSignals[1], fSignals[2], fSignals[3], fSignals[4], fSignals[5], fSignals[6]));
779 }
780
781
782 //_____________________________________________________________________________
783 void AliTRDcluster::SetPadMaskedPosition(UChar_t position)
784 {
785   //
786   // store the pad corruption position code
787   // 
788   // Code: 1 = left cluster
789   //       2 = middle cluster;
790   //       4 = right cluster
791   //
792   for(Int_t ipos = 0; ipos < 3; ipos++)
793     if(TESTBIT(position, ipos))
794       SETBIT(fClusterMasking, ipos);
795 }
796
797 //_____________________________________________________________________________
798 void AliTRDcluster::SetPadMaskedStatus(UChar_t status)
799 {
800   //
801   // store the status of the corrupted pad
802   //
803   // Code: 2 = noisy
804   //       4 = Bridged Left
805   //       8 = Bridged Right
806   //      32 = Not Connected
807   for(Int_t ipos = 0; ipos < 5; ipos++)
808     if(TESTBIT(status, ipos))
809       SETBIT(fClusterMasking, ipos + 3);
810 }
811
812 //___________________________________________________________________________
813 Float_t AliTRDcluster::GetDYcog(Double_t *const, Double_t *const)
814 {
815 //
816 // Get COG position
817 // Used for clusters with more than 3 pads - where LUT not applicable
818 //
819   Double_t sum = fSignals[1]
820                 +fSignals[2]
821                 +fSignals[3] 
822                 +fSignals[4]
823                 +fSignals[5];
824
825   // ???????????? CBL
826   // Go to 3 pad COG ????
827   // ???????????? CBL
828   fCenter = (0.0 * (-fSignals[1] + fSignals[5])
829                       + (-fSignals[2] + fSignals[4])) / sum;
830
831   return fCenter;
832 }
833
834 //___________________________________________________________________________
835 Float_t AliTRDcluster::GetDYlut(Double_t *const, Double_t *const)
836 {
837   //
838   // Calculates the cluster position using the lookup table.
839   // Method provided by Bogdan Vulpescu.
840   //
841
842   if(!fgLUT) FillLUT();
843
844   Double_t ampL = fSignals[2],
845            ampC = fSignals[3],
846            ampR = fSignals[4];
847   Int_t ilayer = AliTRDgeometry::GetLayer(fDetector);
848
849   Double_t x    = 0.0;
850   Double_t xmin, xmax, xwid;
851
852   Int_t    side = 0;
853   Int_t    ix;
854
855   Double_t xMin[AliTRDgeometry::kNlayer] = { 
856     0.006492, 0.006377, 0.006258, 0.006144, 0.006030, 0.005980 
857   };
858   Double_t xMax[AliTRDgeometry::kNlayer] = { 
859     0.960351, 0.965870, 0.970445, 0.974352, 0.977667, 0.996101 
860   };
861
862   if      (ampL > ampR) {
863     x    = (ampL - ampR) / ampC;
864     side = -1;
865   } 
866   else if (ampL < ampR) {
867     x    = (ampR - ampL) / ampC;
868     side = +1;
869   }
870
871   if (ampL != ampR) {
872
873     xmin = xMin[ilayer] + 0.000005;
874     xmax = xMax[ilayer] - 0.000005;
875     xwid = (xmax - xmin) / 127.0;
876
877     if      (x < xmin) fCenter = 0.0000;
878     else if (x > xmax) fCenter = side * 0.5000;
879     else {
880       ix      = (Int_t) ((x - xmin) / xwid);
881       fCenter = side * fgLUT[ilayer*fgkNlut+ix];
882     }
883   } else fCenter = 0.0;
884
885   return fCenter;
886 }
887
888 //___________________________________________________________________________
889 Float_t AliTRDcluster::GetDYgauss(Double_t s2w, Double_t *const y1, Double_t *const y2)
890 {
891 //
892 // (Re)Calculate cluster position in the y direction in local chamber coordinates using all available information from tracking.
893 //
894 // Input parameters:
895 //   s2 - sigma of gaussian parameterization (see bellow for the exact parameterization)
896 //   W  - pad width
897 //   xd - drift length (with respect to the anode wire) [cm]
898 //   wt - omega*tau = tg(a_L)
899 // Output values :
900 //   y1 and y2 - partial positions based on 2 pads clusters
901 //   return y position of the cluster from all information
902 //
903 // Estimation of y coordinate is based on the gaussian approximation of the PRF. Thus one may
904 // calculate the y position knowing the signals q_i-1, q_i and q_i+1 in the 3 adiacent pads by:
905 // BEGIN_LATEX
906 // y = #frac{1}{w_{1}+w_{2}}#[]{w_{1}#(){y_{0}-#frac{W}{2}+#frac{s^{2}}{W}ln#frac{q_{i}}{q_{i-1}}}+w_{2}#(){y_{0}+ #frac{W}{2}+#frac{s^{2}}{W}ln#frac{q_{i+1}}{q_{i}}}}
907 // END_LATEX
908 // where W is the pad width, y_0 is the position of the center pad and s^2 is given by
909 // BEGIN_LATEX
910 // s^{2} = s^{2}_{0} + s^{2}_{diff} (x,B) + #frac{tg^{2}(#phi-#alpha_{L})*l^{2}}{12}
911 // END_LATEX
912 // with s_0 being the PRF for 0 drift and track incidence phi equal to the lorentz angle a_L and the diffusion term 
913 // being described by:
914 // BEGIN_LATEX
915 // s_{diff} (x,B) = #frac{D_{L}#sqrt{x}}{1+#(){#omega#tau}^{2}}
916 // END_LATEX
917 // with x being the drift length. The weights w_1 and w_2 are taken to be q_i-1^2 and q_i+1^2 respectively
918 // 
919 // Authors
920 // Alex Bercuci <A.Bercuci@gsi.de>
921 // Theodor Rascanu <trascanu@stud.uni-frankfurt.de>
922 //
923   Double_t w1 = fSignals[2]*fSignals[2];
924   Double_t w2 = fSignals[4]*fSignals[4];
925   Double_t w = w1+w2;
926   if(w<1.){
927     AliError("Missing side signals for cluster.");
928     Print("a");
929     return 0.;
930   }  
931
932   //Double_t s2w = s2/W/W;
933   Float_t y1r  = fSignals[2]>0 ? (-0.5 + s2w*TMath::Log(fSignals[3]/(Float_t)fSignals[2])) : 0.;
934   Float_t y2r  = fSignals[4]>0 ? (0.5 + s2w*TMath::Log(fSignals[4]/(Float_t)fSignals[3])) : 0.;
935
936   if(y1) (*y1) = y1r;
937   if(y2) (*y2) = y2r;
938
939   return fCenter      = (w1*y1r+w2*y2r)/w;
940 }
941
942
943
944 //_____________________________________________________________________________
945 void AliTRDcluster::FillLUT()
946 {
947   //
948   // Create the LUT
949   //
950
951   // The lookup table from Bogdan
952   Float_t lut[AliTRDgeometry::kNlayer][fgkNlut] = {  
953     {
954       0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611, 
955       0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187, 
956       0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704, 
957       0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160, 
958       0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557, 
959       0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901, 
960       0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207, 
961       0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483, 
962       0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727, 
963       0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952, 
964       0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157, 
965       0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345, 
966       0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520, 
967       0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683, 
968       0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824, 
969       0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978 
970     },
971     {
972       0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642, 
973       0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241, 
974       0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770, 
975       0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229, 
976       0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627, 
977       0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968, 
978       0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266, 
979       0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536, 
980       0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772, 
981       0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991, 
982       0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187, 
983       0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369, 
984       0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538, 
985       0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694, 
986       0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832, 
987       0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
988     },
989     {
990       0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674, 
991       0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296, 
992       0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839, 
993       0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299, 
994       0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693, 
995       0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033, 
996       0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323, 
997       0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586, 
998       0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815, 
999       0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027, 
1000       0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217, 
1001       0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393, 
1002       0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555, 
1003       0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705, 
1004       0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839, 
1005       0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1006     },
1007     {
1008       0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708, 
1009       0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356, 
1010       0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910, 
1011       0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371, 
1012       0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759, 
1013       0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095, 
1014       0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378, 
1015       0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633, 
1016       0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857, 
1017       0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061, 
1018       0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244, 
1019       0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414, 
1020       0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571, 
1021       0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715, 
1022       0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846, 
1023       0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1024     },
1025     {
1026       0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744, 
1027       0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419, 
1028       0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984, 
1029       0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444, 
1030       0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824, 
1031       0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152, 
1032       0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432, 
1033       0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677, 
1034       0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896, 
1035       0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094, 
1036       0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270, 
1037       0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435, 
1038       0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586, 
1039       0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725, 
1040       0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852, 
1041       0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1042     },
1043     {
1044       0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792, 
1045       0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505, 
1046       0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081, 
1047       0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541, 
1048       0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917, 
1049       0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235, 
1050       0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511, 
1051       0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748, 
1052       0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962, 
1053       0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152, 
1054       0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324, 
1055       0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483, 
1056       0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629, 
1057       0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759, 
1058       0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859, 
1059       0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1060     }
1061   }; 
1062
1063   if(!fgLUT) fgLUT = new Double_t[AliTRDgeometry::kNlayer*fgkNlut];
1064
1065   for (Int_t ilayer = 0; ilayer < AliTRDgeometry::kNlayer; ilayer++) {
1066     for (Int_t ilut  = 0; ilut  < fgkNlut; ilut++  ) {
1067       fgLUT[ilayer*fgkNlut+ilut] = lut[ilayer][ilut];
1068     }
1069   }
1070 }
1071