b02d1cd986a2ed88ca0d28d03e421dfc3f24b10c
[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 //  TRD cluster                                                              //
21 //                                                                           //
22 /////////////////////////////////////////////////////////////////////////////// 
23
24 #include "TMath.h"
25
26 #include "AliLog.h"
27 #include "AliTRDcluster.h"
28 #include "AliTRDgeometry.h"
29 #include "AliTRDCommonParam.h"
30 #include "AliTRDtrackletWord.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 AliTRDtrackletWord *const tracklet, Int_t det, UShort_t volid)
118   :AliCluster(volid,tracklet->GetX(),tracklet->GetY(),tracklet->GetZ(),0,0,0)
119   ,fPadCol(0)
120   ,fPadRow(0)
121   ,fPadTime(0)
122   ,fLocalTimeBin(0)
123   ,fNPads(0)
124   ,fClusterMasking(0)
125   ,fDetector(det)
126   ,fQ(0.)
127   ,fCenter(0.)
128 {
129   //
130   // Constructor from online tracklet 
131   //
132 }
133
134 //_____________________________________________________________________________
135 AliTRDcluster::AliTRDcluster(const AliTRDcluster &c)
136   :AliCluster(c)
137   ,fPadCol(c.fPadCol)
138   ,fPadRow(c.fPadRow)
139   ,fPadTime(c.fPadTime)
140   ,fLocalTimeBin(c.fLocalTimeBin)
141   ,fNPads(c.fNPads)
142   ,fClusterMasking(c.fClusterMasking)
143   ,fDetector(c.fDetector)
144   ,fQ(c.fQ)
145   ,fCenter(c.fCenter)
146 {
147   //
148   // Copy constructor 
149   //
150
151   SetLabel(c.GetLabel(0),0);
152   SetLabel(c.GetLabel(1),1);
153   SetLabel(c.GetLabel(2),2);
154
155   SetY(c.GetY());
156   SetZ(c.GetZ());
157   AliCluster::SetSigmaY2(c.GetSigmaY2());
158   SetSigmaZ2(c.GetSigmaZ2());  
159
160   for (Int_t i = 0; i < 7; i++) {
161     fSignals[i] = c.fSignals[i];
162   }
163
164 }
165
166 //_____________________________________________________________________________
167 void AliTRDcluster::AddTrackIndex(Int_t *track)
168 {
169   //
170   // Adds track index. Currently assumed that track is an array of
171   // size 9, and up to 3 track indexes are stored in fTracks[3].
172   // Indexes are sorted according to:
173   //  1) index of max number of appearances is stored first
174   //  2) if two or more indexes appear equal number of times, the lowest
175   //     ones are stored first;
176   //
177
178   const Int_t kSize = 9;
179   Int_t  entries[kSize][2];
180
181   Int_t  i = 0;
182   Int_t  j = 0;
183   Int_t  k = 0;
184   Int_t  index;
185   Bool_t indexAdded;
186
187   for (i = 0; i < kSize; i++) {
188     entries[i][0] = -1;
189     entries[i][1] =  0;
190   }                                 
191
192   for (k = 0; k < kSize; k++) {
193
194     index      = track[k];
195     indexAdded = kFALSE; 
196
197     j = 0;
198     if (index >= 0) {
199       while ((!indexAdded) && (j < kSize)) {
200         if ((entries[j][0] == index) || 
201             (entries[j][1] ==     0)) {
202           entries[j][0] = index;
203           entries[j][1] = entries[j][1] + 1;
204           indexAdded    = kTRUE;
205         }
206         j++;
207       }
208     }
209
210   }
211
212   // Sort by number of appearances and index value
213   Int_t swap = 1;
214   Int_t tmp0;
215   Int_t tmp1;
216   while (swap > 0) {
217     swap = 0;
218     for (i = 0; i < (kSize - 1); i++) {
219       if ((entries[i][0]   >= 0) && 
220           (entries[i+1][0] >= 0)) {
221         if ((entries[i][1] < entries[i+1][1]) ||
222             ((entries[i][1] == entries[i+1][1]) &&
223              (entries[i][0] >  entries[i+1][0]))) {
224           tmp0            = entries[i][0];
225           tmp1            = entries[i][1];
226           entries[i][0]   = entries[i+1][0];
227           entries[i][1]   = entries[i+1][1];
228           entries[i+1][0] = tmp0;
229           entries[i+1][1] = tmp1;
230           swap++;
231         }
232       }
233     }
234   }               
235
236   // Set track indexes
237   for (i = 0; i < 3; i++) {
238     SetLabel(entries[i][0],i);
239   }
240
241   return;
242
243 }          
244
245 //_____________________________________________________________________________
246 void AliTRDcluster::Clear(Option_t *)
247 {
248   //
249   // Reset all member to the default value
250   //
251   fPadCol=0;
252   fPadRow=0;
253   fPadTime=0;
254   fLocalTimeBin=0;
255   fNPads=0;
256   fClusterMasking=0;
257   fDetector=0;
258   for (Int_t i=0; i < 7; i++) fSignals[i]=0;
259   fQ = 0;
260   fCenter = 0;
261   for (Int_t i = 0; i < 3; i++) SetLabel(0,i);
262   SetX(0.);
263   SetY(0.);
264   SetZ(0.);
265   AliCluster::SetSigmaY2(0.);
266   SetSigmaZ2(0.);
267   SetVolumeId(0);
268 }
269
270 //_____________________________________________________________________________
271 Float_t AliTRDcluster::GetSumS() const
272 {
273   //
274   // Returns the total charge from a not unfolded cluster
275   //
276
277   Float_t sum = 0.0;
278   for (Int_t i = 0; i < 7; i++) {
279     sum += fSignals[i];
280   }
281
282   return sum;
283
284 }
285
286 //___________________________________________________________________________
287 Double_t AliTRDcluster::GetSX(Int_t tb, Double_t z)
288 {
289 // Returns the error parameterization in the radial direction for TRD clusters as function of 
290 // the calibrated time bin (tb) and optionally distance to anode wire (z). By default (no z information) 
291 // the mean value over all cluster to wire distance is chosen.
292 // 
293 // There are several contributions which are entering in the definition of the radial errors of the clusters. 
294 // Although an analytic defition should be possible for the moment this is not yet available but instead a 
295 // numerical parameterization is provided (see AliTRDclusterResolution::ProcessSigma() for the calibration 
296 // method). The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. 
297 //
298 //Begin_Html
299 //<img src="TRD/clusterXerrorDiff2D.gif">
300 //End_Html
301 //
302 // Here is a list of uncertainty components:
303 // - Time Response Function (TRF) - the major contribution. since TRF is also not symmetric (even if tail is 
304 //   cancelled) it also creates a systematic shift dependent on the charge distribution before and after the cluster.
305 // - longitudinal diffusion - increase the width of TRF and scales with square root of drift length
306 // - variation in the drift velocity within the drift cell 
307 //
308 // Author
309 // A.Bercuci <A.Bercuci@gsi.de>
310
311   if(tb<1 || tb>=24) return 10.; // return huge [10cm]
312   const Double_t sx[24][10]={
313     {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},
314     {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},
315     {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},
316     {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},
317     {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},
318     {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},
319     {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},
320     {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},
321     {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},
322     {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},
323     {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},
324     {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},
325     {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},
326     {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},
327     {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},
328     {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},
329     {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},
330     {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},
331     {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},
332     {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},
333     {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},
334     {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},
335     {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},
336     {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}
337   };
338   if(z>=0. && z<.25) return sx[tb][Int_t(z/.025)];
339   
340   Double_t m = 0.; for(Int_t id=10; id--;) m+=sx[tb][id];
341   return m*.1;
342 }
343
344 //___________________________________________________________________________
345 Double_t AliTRDcluster::GetSYdrift(Int_t tb, Int_t ly, Double_t/* z*/)
346 {
347 // Returns the error parameterization for TRD clusters as function of the drift length (here calibrated time bin tb)
348 // and optionally distance to anode wire (z) for the LUT r-phi cluster shape. By default (no z information) the largest 
349 // value over all cluster to wire values is chosen.
350 // 
351 // For the LUT method the dependence of s_y with x and d is obtained via a fit to the cluster to MC 
352 // resolution. (see class AliTRDclusterResolution for more details). A normalization to the reference radial position
353 // x0 = 0.675 (tb=5 for ideal vd) is also applied (see GetSYprf()). The function is *NOT* calibration aware ! 
354 // The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. A comparison 
355 // with the GAUS parameterization is also given
356 //
357 // For the GAUS method the dependence of s_y with x is *analytic* and it is expressed by the relation.
358 // BEGIN_LATEX
359 // #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}}
360 // END_LATEX
361 // The result is displayed in the figure below as function of the drift time and compared with the LUT parameterization.
362 //Begin_Html
363 //<img src="TRD/clusterYerrorDiff2D.gif">
364 //<img src="TRD/clusterYerrorDiff1D.gif">
365 //End_Html
366 //
367 // Author
368 // A.Bercuci <A.Bercuci@gsi.de>
369
370
371   if(tb<1 || tb>=24) return 10.; // return huge [10cm]
372   const Float_t lSy[6][24] = {
373     {75.7561, 0.0325, 0.0175, 0.0174, 0.0206, 0.0232,
374      0.0253, 0.0262, 0.0265, 0.0264, 0.0266, 0.0257,
375      0.0258, 0.0261, 0.0259, 0.0253, 0.0257, 0.0261,
376      0.0255, 0.0250, 0.0259, 0.0266, 0.0278, 0.0319
377     },
378     {49.2252, 0.0371, 0.0204, 0.0189, 0.0230, 0.0261,
379      0.0281, 0.0290, 0.0292, 0.0286, 0.0277, 0.0279,
380      0.0285, 0.0281, 0.0291, 0.0281, 0.0281, 0.0282,
381      0.0272, 0.0282, 0.0282, 0.0284, 0.0310, 0.0334
382     },
383     {55.1674, 0.0388, 0.0212, 0.0200, 0.0239, 0.0271,
384      0.0288, 0.0299, 0.0306, 0.0300, 0.0296, 0.0303,
385      0.0293, 0.0290, 0.0291, 0.0294, 0.0295, 0.0290,
386      0.0293, 0.0292, 0.0292, 0.0293, 0.0316, 0.0358
387     },
388     {45.1004, 0.0411, 0.0225, 0.0215, 0.0249, 0.0281,
389      0.0301, 0.0315, 0.0320, 0.0308, 0.0318, 0.0321,
390      0.0312, 0.0311, 0.0316, 0.0315, 0.0310, 0.0308,
391      0.0313, 0.0303, 0.0314, 0.0314, 0.0324, 0.0369
392     },
393     {43.8614, 0.0420, 0.0239, 0.0224, 0.0268, 0.0296,
394      0.0322, 0.0336, 0.0333, 0.0326, 0.0321, 0.0325,
395      0.0329, 0.0326, 0.0323, 0.0322, 0.0326, 0.0320,
396      0.0329, 0.0319, 0.0314, 0.0329, 0.0341, 0.0373
397     },
398     {40.5440, 0.0434, 0.0246, 0.0236, 0.0275, 0.0311,
399      0.0332, 0.0345, 0.0347, 0.0347, 0.0340, 0.0336,
400      0.0339, 0.0344, 0.0339, 0.0341, 0.0341, 0.0342,
401      0.0345, 0.0328, 0.0341, 0.0332, 0.0356, 0.0398
402     },
403   };
404   // adjusted ...
405   return TMath::Max(lSy[ly][tb]-0.0150, 0.0010);
406
407 /*  const Double_t sy[24][10]={
408     {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},
409     {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},
410     {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},
411     {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},
412     {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},
413     {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},
414     {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},
415     {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},
416     {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},
417     {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},
418     {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},
419     {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},
420     {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},
421     {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},
422     {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},
423     {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},
424     {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},
425     {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},
426     {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},
427     {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},
428     {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},
429     {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},
430     {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},
431     {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}
432   };
433   if(z>=0. && z<.25) return sy[tb][Int_t(z/.025)] -  sy[5][Int_t(z/.025)];
434
435   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];
436
437   return m;*/
438 }
439
440 //___________________________________________________________________________
441 Double_t AliTRDcluster::GetSYcharge(Float_t q)
442 {
443 // Parameterization of the r-phi resolution component due to cluster charge.
444 // The value is the offset from the nominal cluster resolution defined as the
445 // cluster resolution at average cluster charge (q0).
446 // 
447 // BEGIN_LATEX
448 // #Delta #sigma_{y}(q) = a*(#frac{1}{q} - #frac{1}{q_{0}})
449 // q_{0} #approx 50
450 // END_LATEX
451 // The definition is *NOT* robust against gain fluctuations and thus two approaches are possible
452 // when residual miscalibration are available:
453 //   - determine parameterization with a resolution matching those of the gain
454 //   - define an analytic model which scales with the gain.
455 // 
456 // For more details please see AliTRDclusterResolution::ProcessCharge()
457 //
458 //Begin_Html
459 //<img src="TRD/clusterQerror.gif">
460 //End_Html
461 //
462 // Author
463 // A.Bercuci <A.Bercuci@gsi.de>
464
465   const Float_t sq0inv = 0.019962; // [1/q0]
466   const Float_t sqb    = 0.037328; // [cm]
467
468   return sqb*(1./q - sq0inv);
469 }
470
471 //___________________________________________________________________________
472 Double_t AliTRDcluster::GetSYprf(Int_t ly, Double_t center, Double_t s2)
473 {
474 // Parameterization of the cluster error in the r-phi direction due to charge sharing between 
475 // adiacent pads. Should be identical to what is provided in the OCDB as PRF [TODO]
476 //
477 // The parameterization is obtained from fitting cluster resolution at phi=exb and |x-0.675|<0.225. 
478 // For more details see AliTRDclusterResolution::ProcessCenter().
479 //
480 //Begin_Html
481 //<img src="TRD/clusterPRFerror.gif">
482 //End_Html
483 //
484 // Author
485 // A.Bercuci <A.Bercuci@gsi.de>
486
487 /*  const Float_t scy[AliTRDgeometry::kNlayer][4] = {
488     {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
489     {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
490     {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
491     {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
492     {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
493     {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
494   };*/
495   const Float_t lPRF[] = {0.438, 0.403, 0.393, 0.382, 0.376, 0.345};
496
497   return s2*TMath::Gaus(center, 0., lPRF[ly]);
498 }
499
500
501 //___________________________________________________________________________
502 Double_t AliTRDcluster::GetXcorr(Int_t tb, Double_t z)
503 {
504 // Drift length correction [cm]. Due to variation of mean drift velocity along the drift region
505 // from nominal vd at xd->infinity. For drift velocity determination based on tracking information 
506 // the correction should be negligible.
507 //Begin_Html
508 //<img src="TRD/clusterXcorr.gif">
509 //End_Html
510 // TODO to be parametrized in term of drift velocity at infinite drift length
511 // A.Bercuci (Mar 28 2009)
512
513   if(tb<0 || tb>=24) return 0.;
514   const Int_t nd = 5;
515   const Double_t dx[24][nd]={
516     {+1.747e-01,+3.195e-01,+1.641e-01,+1.607e-01,+6.002e-01},
517     {+5.468e-02,+5.760e-02,+6.365e-02,+8.003e-02,+1.067e-01},
518     {-6.327e-02,-6.339e-02,-6.423e-02,-6.900e-02,-7.949e-02},
519     {-1.417e-01,-1.424e-01,-1.450e-01,-1.465e-01,-1.514e-01},
520     {-1.637e-01,-1.619e-01,-1.622e-01,-1.613e-01,-1.648e-01},
521     {-1.386e-01,-1.334e-01,-1.261e-01,-1.276e-01,-1.314e-01},
522     {-8.799e-02,-8.299e-02,-7.861e-02,-8.038e-02,-8.436e-02},
523     {-5.139e-02,-4.849e-02,-4.641e-02,-4.965e-02,-5.286e-02},
524     {-2.927e-02,-2.773e-02,-2.807e-02,-3.021e-02,-3.378e-02},
525     {-1.380e-02,-1.229e-02,-1.335e-02,-1.547e-02,-1.984e-02},
526     {-4.168e-03,-4.601e-03,-5.462e-03,-8.164e-03,-1.035e-02},
527     {+2.044e-03,+1.889e-03,+9.603e-04,-1.342e-03,-3.736e-03},
528     {+3.568e-03,+3.581e-03,+2.391e-03,+2.942e-05,-1.585e-03},
529     {+4.403e-03,+4.571e-03,+3.509e-03,+8.703e-04,-1.425e-03},
530     {+4.941e-03,+4.808e-03,+3.284e-03,+1.105e-03,-1.208e-03},
531     {+5.124e-03,+5.022e-03,+4.305e-03,+2.023e-03,-1.145e-03},
532     {+4.882e-03,+4.008e-03,+3.408e-03,+7.886e-04,-1.356e-03},
533     {+3.852e-03,+3.539e-03,+2.057e-03,+1.670e-04,-1.993e-03},
534     {+2.154e-03,+2.111e-03,+5.723e-04,-1.254e-03,-3.256e-03},
535     {+1.755e-03,+2.101e-03,+9.516e-04,-1.649e-03,-3.394e-03},
536     {+1.617e-03,+1.662e-03,+4.169e-04,-9.843e-04,-4.309e-03},
537     {-9.204e-03,-9.069e-03,-1.182e-02,-1.458e-02,-1.880e-02},
538     {-6.727e-02,-6.820e-02,-6.804e-02,-7.134e-02,-7.615e-02},
539     {-1.802e-01,-1.733e-01,-1.633e-01,-1.601e-01,-1.632e-01}
540   };
541 //   const Double_t dx[24][nd]={
542 //     {+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},
543 //     {-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},
544 //     {-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},
545 //     {-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},
546 //     {-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},
547 //     {-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},
548 //     {-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},
549 //     {-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},
550 //     {+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},
551 //     {+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},
552 //     {+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},
553 //     {+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},
554 //     {+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},
555 //     {+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},
556 //     {+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},
557 //     {+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},
558 //     {+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},
559 //     {+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},
560 //     {+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},
561 //     {+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},
562 //     {+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},
563 //     {+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},
564 //     {+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},
565 //     {+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}
566 //   };
567   if(z>=0. && z<.25) return dx[tb][Int_t(z/.025)];
568
569   Double_t m = 0.; for(Int_t id=nd; id--;) m+=dx[tb][id];
570   return m/nd;
571 }
572
573 //___________________________________________________________________________
574 Double_t AliTRDcluster::GetYcorr(Int_t ly, Float_t y)
575 {
576 // PRF correction for the LUT r-phi cluster shape.
577 //Begin_Html
578 //<img src="TRD/clusterYcorr.gif">
579 //End_Html
580
581   const Float_t cy[AliTRDgeometry::kNlayer][3] = {
582     { 4.014e-04, 8.605e-03, -6.880e+00},
583     {-3.061e-04, 9.663e-03, -6.789e+00},
584     { 1.124e-03, 1.105e-02, -6.825e+00},
585     {-1.527e-03, 1.231e-02, -6.777e+00},
586     { 2.150e-03, 1.387e-02, -6.783e+00},
587     {-1.296e-03, 1.486e-02, -6.825e+00}
588   }; 
589
590   return cy[ly][0] + cy[ly][1] * TMath::Sin(cy[ly][2] * y);
591 }
592
593 //_____________________________________________________________________________
594 Float_t AliTRDcluster::GetXloc(Double_t t0, Double_t vd, Double_t *const /*q*/, Double_t *const /*xq*/, Double_t /*z*/)
595 {
596 //
597 // (Re)Calculate cluster position in the x direction in local chamber coordinates (with respect to the anode wire 
598 // position) using all available information from tracking.
599 // Input parameters:
600 //   t0 - calibration aware trigger delay [us]
601 //   vd - drift velocity in the region of the cluster [cm/us]
602 //   z  - distance to the anode wire [cm]. By default average over the drift cell width.
603 //   q & xq - array of charges and cluster positions from previous clusters in the tracklet [a.u.]
604 // Output values :
605 //   return x position of the cluster with respect to the 
606 //   anode wire using all tracking information
607 //
608 // The estimation of the radial position is based on calculating the drift time and the drift velocity at the point of 
609 // estimation. The drift time can be estimated according to the expression:
610 // BEGIN_LATEX
611 // t_{drift} = t_{bin} - t_{0} - t_{cause}(x) - t_{TC}(q_{i-1}, q_{i-2}, ...)
612 // END_LATEX
613 // where t_0 is the delay of the trigger signal. t_cause is the causality delay between ionisation electrons hitting 
614 // the anode and the registration of maximum signal by the electronics - it is due to the rising time of the TRF 
615 // A second order correction here comes from the fact that the time spreading of charge at anode is the convolution of
616 // TRF with the diffusion and thus cross-talk between clusters before and after local clusters changes with drift length. 
617 // t_TC is the residual charge from previous (in time) clusters due to residual tails after tail cancellation. 
618 // This tends to push cluster forward and depends on the magnitude of their charge.
619 //
620 // The drift velocity varies with the drift length (and distance to anode wire) as described by cell structure simulation. 
621 // Thus one, in principle, can calculate iteratively the drift length from the expression:
622 // BEGIN_LATEX
623 // x = t_{drift}(x)*v_{drift}(x)
624 // END_LATEX
625 // In practice we use a numerical approach (AliTRDcluster::GetXcorr()) to correct for anisochronity obtained from MC 
626 // comparison (see AliTRDclusterResolution::ProcessSigma()). Also the calibration of 0 approximation (no x dependence)
627 // for t_cause is obtained from MC comparisons and impossible to disentangle in real life from trigger delay.
628 //
629 // Author
630 // Alex Bercuci <A.Bercuci@gsi.de>
631 //
632
633   AliTRDCommonParam *cp = AliTRDCommonParam::Instance(); 
634   Double_t fFreq = cp->GetSamplingFrequency();
635
636   //drift time corresponding to the center of the time bin
637   Double_t td = (fPadTime + .5)/fFreq; // [us] 
638   // correction for t0
639   td -= t0;
640   // time bin corrected for t0
641   // Bug in TMath::Nint().root-5.23.02
642   // TMath::Nint(3.5) = 4 and TMath::Nint(4.5) = 4
643   Double_t tmp = td*fFreq;
644   fLocalTimeBin = Char_t(TMath::Floor(tmp));
645   if(tmp-fLocalTimeBin > .5) fLocalTimeBin++;
646   if(td < .2) return 0.;
647   // TRF rising time (fitted)
648   // It should be absorbed by the t0. For the moment t0 is 0 for simulations.
649   // A.Bercuci (Mar 26 2009)
650   td -= 0.189;
651
652   // apply fitted correction 
653   Float_t x = td*vd + GetXcorr(fLocalTimeBin);
654   if(x>0.&&x<.5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()) SetInChamber();
655
656   return x;
657
658 /*
659   // calculate radial posion of clusters in the drift region
660
661   // invert drift time function
662   Double_t xM= AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght(),
663            x = vd*td + .5*AliTRDgeometry::CamHght(), 
664            t = cp->TimeStruct(vd, x, z), dx1=0.,dx2;
665   while(TMath::Abs(td-t)>1.e-4){ // convergence on 100ps
666     dx2 = vd*(td-t);
667     if(TMath::Abs(TMath::Abs(dx2)-TMath::Abs(dx1))<1.e-6){
668       x+=.5*dx2;
669       break;
670     } else x+=dx2;
671
672     if(x<0. || x>xM) return 0.;
673     t = cp->TimeStruct(vd, x, z);
674     dx1 = dx2;
675   }
676
677   return x-.5*AliTRDgeometry::CamHght();
678 */
679 }
680
681 //_____________________________________________________________________________
682 Float_t AliTRDcluster::GetYloc(Double_t y0, Double_t s2, Double_t W, Double_t *const y1, Double_t *const y2)
683 {
684 // Calculate, in tracking cooordinate system, the r-phi offset the cluster from the middle of the center pad. Three possible methods are implemented:
685 //   - Center of Gravity (COG) see AliTRDcluster::GetDYcog()
686 //   - Look-up Table (LUT) see AliTRDcluster::GetDYlut()
687 //   - Gauss shape (GAUS) see AliTRDcluster::GetDYgauss()
688 // In addition for the case of LUT method position corrections are also applied (see AliTRDcluster::GetYcorr())
689
690   if(IsRPhiMethod(kCOG)) GetDYcog();
691   else if(IsRPhiMethod(kLUT)) GetDYlut();
692   else if(IsRPhiMethod(kGAUS)) GetDYgauss(s2/W/W, y1, y2);
693   else return 0.;
694
695   if(y1) (*y1)*=W;
696   if(y2) (*y2)*=W;
697
698   return y0+fCenter*W+(IsRPhiMethod(kLUT)?GetYcorr(AliTRDgeometry::GetLayer(fDetector), fCenter):0.);
699 }
700
701 //___________________________________________________________________________
702 void AliTRDcluster::SetSigmaY2(Float_t s2, Float_t dt, Float_t exb, Float_t x, Float_t z, Float_t tgp)
703 {
704 // Set variance of TRD cluster in the r-phi direction for each method.
705 // Parameters :
706 //   - s2  - variance due to PRF width for the case of Gauss model. Replaced by parameterization in case of LUT.
707 //   - dt  - transversal diffusion coeficient 
708 //   - exb - tg of lorentz angle
709 //   - x   - drift length - with respect to the anode wire
710 //   - z   - offset from the anode wire
711 //   - tgp - local tangent of the track momentum azimuthal angle
712 //
713 // The ingredients from which the error is computed are:
714 //   - PRF (charge sharing on adjacent pads)  - see AliTRDcluster::GetSYprf()
715 //   - diffusion (dependence with drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSYdrift()
716 //   - charge of the cluster (complex dependence on gain and tail cancellation) - see AliTRDcluster::GetSYcharge()
717 //   - lorentz angle (dependence on the drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSX()
718 //   - track angle (superposition of charges on the anode wire) - see AliTRDseedV1::Fit()
719 //   - projection of radial(x) error on r-phi due to fixed value assumed in tracking for x - see AliTRDseedV1::Fit()
720 //
721 // The last 2 contributions to cluster error can be estimated only during tracking when the track angle 
722 // is known (tgp). For this reason the errors (and optional position) of TRD clusters are recalculated during 
723 // tracking and thus clusters attached to tracks might differ from bare clusters.
724 // 
725 // Taking into account all contributions one can write the the TRD cluster error parameterization as:
726 // BEGIN_LATEX
727 // #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
728 // END_LATEX
729 // From this formula one can deduce a that the simplest calibration method for PRF and diffusion contributions is 
730 // by measuring resolution at B=0T and phi=0. To disentangle further the two remaining contributions one has 
731 // to represent s2 as a function of drift length. 
732 // 
733 // In the gaussian model the diffusion contribution can be expressed as:
734 // BEGIN_LATEX
735 // #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}}
736 // END_LATEX
737 // thus resulting the PRF contribution. For the case of the LUT model both contributions have to be determined from 
738 // the fit (see AliTRDclusterResolution::ProcessCenter() for details).
739 // 
740 // Author:
741 // A.Bercuci <A.Bercuci@gsi.de>
742
743   Float_t sigmaY2 = 0.;
744   Int_t ly = AliTRDgeometry::GetLayer(fDetector);
745   if(IsRPhiMethod(kCOG)) sigmaY2 = 4.e-4;
746   else if(IsRPhiMethod(kLUT)){ 
747     Float_t sd = GetSYdrift(fLocalTimeBin, ly, z); //printf("drift[%6.2f] ", 1.e4*sd);
748     sigmaY2 = GetSYprf(ly, fCenter, sd); //printf("PRF[%6.2f] ", 1.e4*sigmaY2);
749     // add charge contribution TODO scale with respect to s2
750     sigmaY2+= GetSYcharge(TMath::Abs(fQ)); //printf("Q[%6.2f] ", 1.e4*sigmaY2);
751     sigmaY2 = TMath::Max(sigmaY2, Float_t(0.0010)); //!! protection 
752     sigmaY2*= sigmaY2;
753   } else if(IsRPhiMethod(kGAUS)){
754     // PRF contribution
755     sigmaY2 = s2;
756     // Diffusion contribution
757     Double_t sD2 = dt/(1.+exb); sD2 *= sD2; sD2 *= x;
758     sigmaY2+= sD2; 
759     // add charge contribution TODO scale with respect to s2
760     //sigmaY2+= GetSYcharge(TMath::Abs(fQ));
761   }
762
763   // store tg^2(phi-a_L) and tg^2(a_L)
764   Double_t tgg = (tgp-exb)/(1.+tgp*exb); tgg *= tgg;
765   Double_t exb2= exb*exb;
766
767   // Lorentz angle shift contribution 
768   Float_t sx = GetSX(fLocalTimeBin, z); sx*=sx;
769   sigmaY2+= exb2*sx; //printf("Al[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2));
770
771   // Radial contribution due to not measuring x in Kalman model 
772   sigmaY2+= tgg*sx; //printf("x[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2));
773
774   // Track angle contribution
775   sigmaY2+= tgg*x*x*exb2/12.; //printf("angle[%6.2f]\n", 1.e4*TMath::Sqrt(sigmaY2));
776
777   AliCluster::SetSigmaY2(sigmaY2);
778 }
779
780 //_____________________________________________________________________________
781 Bool_t AliTRDcluster::IsEqual(const TObject *o) const
782 {
783   //
784   // Compare relevant information of this cluster with another one
785   //
786   
787   const AliTRDcluster *inCluster = dynamic_cast<const AliTRDcluster*>(o);
788   if (!o || !inCluster) return kFALSE;
789
790   if ( AliCluster::GetX() != inCluster->GetX() ) return kFALSE;
791   if ( AliCluster::GetY() != inCluster->GetY() ) return kFALSE;
792   if ( AliCluster::GetZ() != inCluster->GetZ() ) return kFALSE;
793   if ( fQ != inCluster->fQ ) return kFALSE;
794   if ( fDetector != inCluster->fDetector ) return kFALSE;
795   if ( fPadCol != inCluster->fPadCol ) return kFALSE;
796   if ( fPadRow != inCluster->fPadRow ) return kFALSE;
797   if ( fPadTime != inCluster->fPadTime ) return kFALSE;
798   if ( fClusterMasking != inCluster->fClusterMasking ) return kFALSE;
799   if ( IsInChamber() != inCluster->IsInChamber() ) return kFALSE;
800   if ( IsShared() != inCluster->IsShared() ) return kFALSE;
801   if ( IsUsed() != inCluster->IsUsed() ) return kFALSE;
802   
803   return kTRUE;
804 }
805
806 //_____________________________________________________________________________
807 void AliTRDcluster::Print(Option_t *o) const
808 {
809   AliInfo(Form("Det[%3d] LTrC[%+6.2f %+6.2f %+6.2f] Q[%5.1f] FLAG[in(%c) use(%c) sh(%c)] Y[%s]", 
810     fDetector, GetX(), GetY(), GetZ(), fQ, 
811     IsInChamber() ? 'y' : 'n', 
812     IsUsed() ? 'y' : 'n', 
813     IsShared() ? 'y' : 'n',
814     IsRPhiMethod(kGAUS)?"GAUS":(IsRPhiMethod(kLUT)?"LUT":"COG")
815   ));
816
817   if(strcmp(o, "a")!=0) return;
818   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)); 
819   AliInfo(Form("Signals[%3d %3d %3d %3d %3d %3d %3d]", fSignals[0], fSignals[1], fSignals[2], fSignals[3], fSignals[4], fSignals[5], fSignals[6]));
820 }
821
822
823 //_____________________________________________________________________________
824 void AliTRDcluster::SetPadMaskedPosition(UChar_t position)
825 {
826   //
827   // store the pad corruption position code
828   // 
829   // Code: 1 = left cluster
830   //       2 = middle cluster;
831   //       4 = right cluster
832   //
833   for(Int_t ipos = 0; ipos < 3; ipos++)
834     if(TESTBIT(position, ipos))
835       SETBIT(fClusterMasking, ipos);
836 }
837
838 //_____________________________________________________________________________
839 void AliTRDcluster::SetPadMaskedStatus(UChar_t status)
840 {
841   //
842   // store the status of the corrupted pad
843   //
844   // Code: 2 = noisy
845   //       4 = Bridged Left
846   //       8 = Bridged Right
847   //      32 = Not Connected
848   for(Int_t ipos = 0; ipos < 5; ipos++)
849     if(TESTBIT(status, ipos))
850       SETBIT(fClusterMasking, ipos + 3);
851 }
852
853 //___________________________________________________________________________
854 Float_t AliTRDcluster::GetDYcog(Double_t *const, Double_t *const)
855 {
856 //
857 // Get COG position
858 // Used for clusters with more than 3 pads - where LUT not applicable
859 //
860   Double_t sum = fSignals[1]
861                 +fSignals[2]
862                 +fSignals[3] 
863                 +fSignals[4]
864                 +fSignals[5];
865
866   // ???????????? CBL
867   // Go to 3 pad COG ????
868   // ???????????? CBL
869   fCenter = (0.0 * (-fSignals[1] + fSignals[5])
870                       + (-fSignals[2] + fSignals[4])) / sum;
871
872   return fCenter;
873 }
874
875 //___________________________________________________________________________
876 Float_t AliTRDcluster::GetDYlut(Double_t *const, Double_t *const)
877 {
878   //
879   // Calculates the cluster position using the lookup table.
880   // Method provided by Bogdan Vulpescu.
881   //
882
883   if(!fgLUT) FillLUT();
884
885   Double_t ampL = fSignals[2],
886            ampC = fSignals[3],
887            ampR = fSignals[4];
888   Int_t ilayer = AliTRDgeometry::GetLayer(fDetector);
889
890   Double_t x    = 0.0;
891   Double_t xmin, xmax, xwid;
892
893   Int_t    side = 0;
894   Int_t    ix;
895
896   Double_t xMin[AliTRDgeometry::kNlayer] = { 
897     0.006492, 0.006377, 0.006258, 0.006144, 0.006030, 0.005980 
898   };
899   Double_t xMax[AliTRDgeometry::kNlayer] = { 
900     0.960351, 0.965870, 0.970445, 0.974352, 0.977667, 0.996101 
901   };
902
903   if      (ampL > ampR) {
904     x    = (ampL - ampR) / ampC;
905     side = -1;
906   } 
907   else if (ampL < ampR) {
908     x    = (ampR - ampL) / ampC;
909     side = +1;
910   }
911
912   if (ampL != ampR) {
913
914     xmin = xMin[ilayer] + 0.000005;
915     xmax = xMax[ilayer] - 0.000005;
916     xwid = (xmax - xmin) / 127.0;
917
918     if      (x < xmin) fCenter = 0.0000;
919     else if (x > xmax) fCenter = side * 0.5000;
920     else {
921       ix      = (Int_t) ((x - xmin) / xwid);
922       fCenter = side * fgLUT[ilayer*fgkNlut+ix];
923     }
924   } else fCenter = 0.0;
925
926   return fCenter;
927 }
928
929 //___________________________________________________________________________
930 Float_t AliTRDcluster::GetDYgauss(Double_t s2w, Double_t *const y1, Double_t *const y2)
931 {
932 //
933 // (Re)Calculate cluster position in the y direction in local chamber coordinates using all available information from tracking.
934 //
935 // Input parameters:
936 //   s2 - sigma of gaussian parameterization (see bellow for the exact parameterization)
937 //   W  - pad width
938 //   xd - drift length (with respect to the anode wire) [cm]
939 //   wt - omega*tau = tg(a_L)
940 // Output values :
941 //   y1 and y2 - partial positions based on 2 pads clusters
942 //   return y position of the cluster from all information
943 //
944 // Estimation of y coordinate is based on the gaussian approximation of the PRF. Thus one may
945 // calculate the y position knowing the signals q_i-1, q_i and q_i+1 in the 3 adiacent pads by:
946 // BEGIN_LATEX
947 // 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}}}}
948 // END_LATEX
949 // where W is the pad width, y_0 is the position of the center pad and s^2 is given by
950 // BEGIN_LATEX
951 // s^{2} = s^{2}_{0} + s^{2}_{diff} (x,B) + #frac{tg^{2}(#phi-#alpha_{L})*l^{2}}{12}
952 // END_LATEX
953 // with s_0 being the PRF for 0 drift and track incidence phi equal to the lorentz angle a_L and the diffusion term 
954 // being described by:
955 // BEGIN_LATEX
956 // s_{diff} (x,B) = #frac{D_{L}#sqrt{x}}{1+#(){#omega#tau}^{2}}
957 // END_LATEX
958 // 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
959 // 
960 // Authors
961 // Alex Bercuci <A.Bercuci@gsi.de>
962 // Theodor Rascanu <trascanu@stud.uni-frankfurt.de>
963 //
964   Double_t w1 = fSignals[2]*fSignals[2];
965   Double_t w2 = fSignals[4]*fSignals[4];
966   Double_t w = w1+w2;
967   if(w<1.){
968     AliError("Missing side signals for cluster.");
969     Print("a");
970     return 0.;
971   }  
972
973   //Double_t s2w = s2/W/W;
974   Float_t y1r  = fSignals[2]>0 ? (-0.5 + s2w*TMath::Log(fSignals[3]/(Float_t)fSignals[2])) : 0.;
975   Float_t y2r  = fSignals[4]>0 ? (0.5 + s2w*TMath::Log(fSignals[4]/(Float_t)fSignals[3])) : 0.;
976
977   if(y1) (*y1) = y1r;
978   if(y2) (*y2) = y2r;
979
980   return fCenter      = (w1*y1r+w2*y2r)/w;
981 }
982
983
984
985 //_____________________________________________________________________________
986 void AliTRDcluster::FillLUT()
987 {
988   //
989   // Create the LUT
990   //
991
992   // The lookup table from Bogdan
993   Float_t lut[AliTRDgeometry::kNlayer][fgkNlut] = {  
994     {
995       0.0070, 0.0150, 0.0224, 0.0298, 0.0374, 0.0454, 0.0533, 0.0611, 
996       0.0684, 0.0755, 0.0827, 0.0900, 0.0975, 0.1049, 0.1120, 0.1187, 
997       0.1253, 0.1318, 0.1385, 0.1453, 0.1519, 0.1584, 0.1646, 0.1704, 
998       0.1762, 0.1821, 0.1879, 0.1938, 0.1996, 0.2053, 0.2108, 0.2160, 
999       0.2210, 0.2260, 0.2310, 0.2361, 0.2411, 0.2461, 0.2509, 0.2557, 
1000       0.2602, 0.2646, 0.2689, 0.2732, 0.2774, 0.2816, 0.2859, 0.2901, 
1001       0.2942, 0.2983, 0.3022, 0.3061, 0.3099, 0.3136, 0.3172, 0.3207, 
1002       0.3242, 0.3278, 0.3312, 0.3347, 0.3382, 0.3416, 0.3450, 0.3483, 
1003       0.3515, 0.3547, 0.3579, 0.3609, 0.3639, 0.3669, 0.3698, 0.3727, 
1004       0.3756, 0.3785, 0.3813, 0.3842, 0.3870, 0.3898, 0.3926, 0.3952, 
1005       0.3979, 0.4005, 0.4032, 0.4057, 0.4082, 0.4108, 0.4132, 0.4157, 
1006       0.4181, 0.4205, 0.4228, 0.4252, 0.4275, 0.4299, 0.4322, 0.4345, 
1007       0.4367, 0.4390, 0.4412, 0.4434, 0.4456, 0.4478, 0.4499, 0.4520, 
1008       0.4541, 0.4562, 0.4583, 0.4603, 0.4623, 0.4643, 0.4663, 0.4683, 
1009       0.4702, 0.4722, 0.4741, 0.4758, 0.4774, 0.4790, 0.4805, 0.4824, 
1010       0.4844, 0.4863, 0.4883, 0.4902, 0.4921, 0.4940, 0.4959, 0.4978 
1011     },
1012     {
1013       0.0072, 0.0156, 0.0235, 0.0313, 0.0394, 0.0478, 0.0561, 0.0642, 
1014       0.0718, 0.0792, 0.0868, 0.0947, 0.1025, 0.1101, 0.1172, 0.1241, 
1015       0.1309, 0.1378, 0.1449, 0.1518, 0.1586, 0.1650, 0.1710, 0.1770, 
1016       0.1830, 0.1891, 0.1952, 0.2011, 0.2070, 0.2125, 0.2177, 0.2229, 
1017       0.2280, 0.2332, 0.2383, 0.2435, 0.2484, 0.2533, 0.2581, 0.2627, 
1018       0.2670, 0.2714, 0.2757, 0.2799, 0.2842, 0.2884, 0.2927, 0.2968, 
1019       0.3008, 0.3048, 0.3086, 0.3123, 0.3159, 0.3195, 0.3231, 0.3266, 
1020       0.3301, 0.3335, 0.3370, 0.3404, 0.3438, 0.3471, 0.3504, 0.3536, 
1021       0.3567, 0.3598, 0.3628, 0.3657, 0.3686, 0.3715, 0.3744, 0.3772, 
1022       0.3800, 0.3828, 0.3856, 0.3884, 0.3911, 0.3938, 0.3965, 0.3991, 
1023       0.4016, 0.4042, 0.4067, 0.4092, 0.4116, 0.4140, 0.4164, 0.4187, 
1024       0.4211, 0.4234, 0.4257, 0.4280, 0.4302, 0.4325, 0.4347, 0.4369, 
1025       0.4391, 0.4413, 0.4434, 0.4456, 0.4477, 0.4497, 0.4518, 0.4538, 
1026       0.4558, 0.4578, 0.4598, 0.4618, 0.4637, 0.4656, 0.4675, 0.4694, 
1027       0.4713, 0.4732, 0.4750, 0.4766, 0.4781, 0.4797, 0.4813, 0.4832, 
1028       0.4851, 0.4870, 0.4888, 0.4906, 0.4925, 0.4942, 0.4960, 0.4978
1029     },
1030     {
1031       0.0075, 0.0163, 0.0246, 0.0328, 0.0415, 0.0504, 0.0592, 0.0674, 
1032       0.0753, 0.0832, 0.0914, 0.0996, 0.1077, 0.1154, 0.1225, 0.1296, 
1033       0.1369, 0.1442, 0.1515, 0.1585, 0.1652, 0.1714, 0.1776, 0.1839, 
1034       0.1902, 0.1965, 0.2025, 0.2085, 0.2141, 0.2194, 0.2247, 0.2299, 
1035       0.2352, 0.2405, 0.2457, 0.2507, 0.2557, 0.2604, 0.2649, 0.2693, 
1036       0.2737, 0.2780, 0.2823, 0.2867, 0.2909, 0.2951, 0.2992, 0.3033, 
1037       0.3072, 0.3110, 0.3146, 0.3182, 0.3218, 0.3253, 0.3288, 0.3323, 
1038       0.3357, 0.3392, 0.3426, 0.3459, 0.3492, 0.3524, 0.3555, 0.3586, 
1039       0.3616, 0.3645, 0.3674, 0.3703, 0.3731, 0.3759, 0.3787, 0.3815, 
1040       0.3843, 0.3870, 0.3897, 0.3925, 0.3950, 0.3976, 0.4002, 0.4027, 
1041       0.4052, 0.4076, 0.4101, 0.4124, 0.4148, 0.4171, 0.4194, 0.4217, 
1042       0.4239, 0.4262, 0.4284, 0.4306, 0.4328, 0.4350, 0.4371, 0.4393, 
1043       0.4414, 0.4435, 0.4455, 0.4476, 0.4496, 0.4516, 0.4536, 0.4555, 
1044       0.4575, 0.4594, 0.4613, 0.4632, 0.4650, 0.4669, 0.4687, 0.4705, 
1045       0.4723, 0.4741, 0.4758, 0.4773, 0.4789, 0.4804, 0.4821, 0.4839, 
1046       0.4857, 0.4875, 0.4893, 0.4910, 0.4928, 0.4945, 0.4961, 0.4978
1047     },
1048     {
1049       0.0078, 0.0171, 0.0258, 0.0345, 0.0438, 0.0532, 0.0624, 0.0708, 
1050       0.0791, 0.0875, 0.0962, 0.1048, 0.1130, 0.1206, 0.1281, 0.1356, 
1051       0.1432, 0.1508, 0.1582, 0.1651, 0.1716, 0.1780, 0.1845, 0.1910, 
1052       0.1975, 0.2038, 0.2099, 0.2155, 0.2210, 0.2263, 0.2317, 0.2371, 
1053       0.2425, 0.2477, 0.2528, 0.2578, 0.2626, 0.2671, 0.2715, 0.2759, 
1054       0.2803, 0.2846, 0.2890, 0.2933, 0.2975, 0.3016, 0.3056, 0.3095, 
1055       0.3132, 0.3168, 0.3204, 0.3239, 0.3274, 0.3309, 0.3344, 0.3378, 
1056       0.3412, 0.3446, 0.3479, 0.3511, 0.3543, 0.3574, 0.3603, 0.3633, 
1057       0.3662, 0.3690, 0.3718, 0.3747, 0.3774, 0.3802, 0.3829, 0.3857, 
1058       0.3883, 0.3910, 0.3936, 0.3962, 0.3987, 0.4012, 0.4037, 0.4061, 
1059       0.4085, 0.4109, 0.4132, 0.4155, 0.4177, 0.4200, 0.4222, 0.4244, 
1060       0.4266, 0.4288, 0.4309, 0.4331, 0.4352, 0.4373, 0.4394, 0.4414, 
1061       0.4435, 0.4455, 0.4475, 0.4494, 0.4514, 0.4533, 0.4552, 0.4571, 
1062       0.4590, 0.4608, 0.4626, 0.4645, 0.4662, 0.4680, 0.4698, 0.4715, 
1063       0.4733, 0.4750, 0.4766, 0.4781, 0.4796, 0.4812, 0.4829, 0.4846, 
1064       0.4863, 0.4880, 0.4897, 0.4914, 0.4930, 0.4946, 0.4963, 0.4979
1065     },
1066     {
1067       0.0081, 0.0178, 0.0270, 0.0364, 0.0463, 0.0562, 0.0656, 0.0744, 
1068       0.0831, 0.0921, 0.1013, 0.1102, 0.1183, 0.1261, 0.1339, 0.1419, 
1069       0.1499, 0.1576, 0.1648, 0.1715, 0.1782, 0.1849, 0.1917, 0.1984, 
1070       0.2048, 0.2110, 0.2167, 0.2223, 0.2278, 0.2333, 0.2389, 0.2444, 
1071       0.2497, 0.2548, 0.2598, 0.2645, 0.2691, 0.2735, 0.2780, 0.2824, 
1072       0.2868, 0.2912, 0.2955, 0.2997, 0.3038, 0.3078, 0.3116, 0.3152, 
1073       0.3188, 0.3224, 0.3259, 0.3294, 0.3329, 0.3364, 0.3398, 0.3432, 
1074       0.3465, 0.3497, 0.3529, 0.3561, 0.3591, 0.3620, 0.3649, 0.3677, 
1075       0.3705, 0.3733, 0.3761, 0.3788, 0.3816, 0.3843, 0.3869, 0.3896, 
1076       0.3922, 0.3948, 0.3973, 0.3998, 0.4022, 0.4047, 0.4070, 0.4094, 
1077       0.4117, 0.4139, 0.4162, 0.4184, 0.4206, 0.4227, 0.4249, 0.4270, 
1078       0.4291, 0.4313, 0.4334, 0.4354, 0.4375, 0.4395, 0.4415, 0.4435, 
1079       0.4455, 0.4474, 0.4493, 0.4512, 0.4531, 0.4550, 0.4568, 0.4586, 
1080       0.4604, 0.4622, 0.4639, 0.4657, 0.4674, 0.4691, 0.4708, 0.4725, 
1081       0.4742, 0.4758, 0.4773, 0.4788, 0.4803, 0.4819, 0.4836, 0.4852, 
1082       0.4869, 0.4885, 0.4901, 0.4917, 0.4933, 0.4948, 0.4964, 0.4979
1083     },
1084     {
1085       0.0085, 0.0189, 0.0288, 0.0389, 0.0497, 0.0603, 0.0699, 0.0792, 
1086       0.0887, 0.0985, 0.1082, 0.1170, 0.1253, 0.1336, 0.1421, 0.1505, 
1087       0.1587, 0.1662, 0.1733, 0.1803, 0.1874, 0.1945, 0.2014, 0.2081, 
1088       0.2143, 0.2201, 0.2259, 0.2316, 0.2374, 0.2431, 0.2487, 0.2541, 
1089       0.2593, 0.2642, 0.2689, 0.2735, 0.2781, 0.2826, 0.2872, 0.2917, 
1090       0.2961, 0.3003, 0.3045, 0.3086, 0.3125, 0.3162, 0.3198, 0.3235, 
1091       0.3270, 0.3306, 0.3342, 0.3377, 0.3411, 0.3446, 0.3479, 0.3511, 
1092       0.3543, 0.3575, 0.3605, 0.3634, 0.3663, 0.3691, 0.3720, 0.3748, 
1093       0.3775, 0.3803, 0.3830, 0.3857, 0.3884, 0.3911, 0.3937, 0.3962, 
1094       0.3987, 0.4012, 0.4036, 0.4060, 0.4084, 0.4107, 0.4129, 0.4152, 
1095       0.4174, 0.4196, 0.4218, 0.4239, 0.4261, 0.4282, 0.4303, 0.4324, 
1096       0.4344, 0.4365, 0.4385, 0.4405, 0.4425, 0.4445, 0.4464, 0.4483, 
1097       0.4502, 0.4521, 0.4539, 0.4558, 0.4576, 0.4593, 0.4611, 0.4629, 
1098       0.4646, 0.4663, 0.4680, 0.4697, 0.4714, 0.4730, 0.4747, 0.4759, 
1099       0.4769, 0.4780, 0.4790, 0.4800, 0.4811, 0.4827, 0.4843, 0.4859, 
1100       0.4874, 0.4889, 0.4905, 0.4920, 0.4935, 0.4950, 0.4965, 0.4979
1101     }
1102   }; 
1103
1104   if(!fgLUT) fgLUT = new Double_t[AliTRDgeometry::kNlayer*fgkNlut];
1105
1106   for (Int_t ilayer = 0; ilayer < AliTRDgeometry::kNlayer; ilayer++) {
1107     for (Int_t ilut  = 0; ilut  < fgkNlut; ilut++  ) {
1108       fgLUT[ilayer*fgkNlut+ilut] = lut[ilayer][ilut];
1109     }
1110   }
1111 }
1112