]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDcluster.cxx
coverity fix
[u/mrichter/AliRoot.git] / TRD / AliTRDcluster.cxx
index 58bc599f91abeb30b833ff57362424c342df15d1..4c612bfcbb0d5a5b09651a90a9e5e63b9b458f1b 100644 (file)
@@ -14,7 +14,6 @@
  **************************************************************************/
 
 /* $Id$ */
-
  
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
 //                                                                           //
 /////////////////////////////////////////////////////////////////////////////// 
 
-#include "TMath.h"
+#include <TMath.h>
 
 #include "AliLog.h"
+
 #include "AliTRDcluster.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDCommonParam.h"
+#include "AliTRDtrackletWord.h"
 
 ClassImp(AliTRDcluster)
 
@@ -54,11 +55,12 @@ AliTRDcluster::AliTRDcluster()
   for (Int_t i = 0; i < 7; i++) {
     fSignals[i] = 0;
   }
-
+  SetBit(kLUT);
 }
 
 //___________________________________________________________________________
-AliTRDcluster::AliTRDcluster(Int_t det, UChar_t col, UChar_t row, UChar_t time, const Short_t *sig, UShort_t vid) 
+AliTRDcluster::AliTRDcluster(Int_t det, UChar_t col, UChar_t row, UChar_t time
+                           , const Short_t *sig, UShort_t vid) 
   :AliCluster() 
   ,fPadCol(col)
   ,fPadRow(row)
@@ -79,12 +81,13 @@ AliTRDcluster::AliTRDcluster(Int_t det, UChar_t col, UChar_t row, UChar_t time,
   memcpy(&fSignals, sig, 7*sizeof(Short_t));
   fQ = fSignals[2]+fSignals[3]+fSignals[4];
   SetVolumeId(vid);
+  SetBit(kLUT);
 }
 
 //___________________________________________________________________________
 AliTRDcluster::AliTRDcluster(Int_t det, Float_t q
                            , Float_t *pos, Float_t *sig
-                           , Int_t *tracks, Char_t npads, Short_t *signals
+                           , Int_t *tracks, Char_t npads, Short_t * const signals
                            , UChar_t col, UChar_t row, UChar_t time
                            , Char_t timebin, Float_t center, UShort_t volid)
   :AliCluster(volid,pos[0],pos[1],pos[2],sig[0],sig[1],0.0,0x0) 
@@ -109,6 +112,29 @@ AliTRDcluster::AliTRDcluster(Int_t det, Float_t q
   if (tracks) {
     AddTrackIndex(tracks);
   }
+  SetBit(kLUT);
+}
+
+//_____________________________________________________________________________
+AliTRDcluster::AliTRDcluster(const AliTRDtrackletWord *const tracklet, Int_t det, UShort_t volid)
+  :AliCluster(volid,tracklet->GetX(),tracklet->GetY(),tracklet->GetZ(),0,0,0)
+  ,fPadCol(0)
+  ,fPadRow(0)
+  ,fPadTime(0)
+  ,fLocalTimeBin(0)
+  ,fNPads(0)
+  ,fClusterMasking(0)
+  ,fDetector(det)
+  ,fQ(0.)
+  ,fCenter(0.)
+{
+  //
+  // Constructor from online tracklet 
+  //
+
+  for (Int_t i = 0; i < 7; i++) {
+    fSignals[i] = 0;
+  }
 
 }
 
@@ -129,14 +155,13 @@ AliTRDcluster::AliTRDcluster(const AliTRDcluster &c)
   // Copy constructor 
   //
 
-  SetBit(kInChamber, c.IsInChamber());
   SetLabel(c.GetLabel(0),0);
   SetLabel(c.GetLabel(1),1);
   SetLabel(c.GetLabel(2),2);
 
   SetY(c.GetY());
   SetZ(c.GetZ());
-  SetSigmaY2(c.GetSigmaY2());
+  AliCluster::SetSigmaY2(c.GetSigmaY2());
   SetSigmaZ2(c.GetSigmaZ2());  
 
   for (Int_t i = 0; i < 7; i++) {
@@ -146,7 +171,47 @@ AliTRDcluster::AliTRDcluster(const AliTRDcluster &c)
 }
 
 //_____________________________________________________________________________
-void AliTRDcluster::AddTrackIndex(Int_t *track)
+AliTRDcluster &AliTRDcluster::operator=(const AliTRDcluster &c)
+{
+  //
+  // Assignment operator
+  //
+
+  if (&c == this) {
+     return *this;
+  }
+
+  // Call the assignment operator of the base class
+  AliCluster::operator=(c);
+
+  fPadCol         = c.fPadCol;
+  fPadRow         = c.fPadRow;
+  fPadTime        = c.fPadTime;
+  fLocalTimeBin   = c.fLocalTimeBin;
+  fNPads          = c.fNPads;
+  fClusterMasking = c.fClusterMasking;
+  fDetector       = c.fDetector;
+  fQ              = c.fQ;
+  fCenter         = c.fCenter;
+
+  SetLabel(c.GetLabel(0),0);
+  SetLabel(c.GetLabel(1),1);
+  SetLabel(c.GetLabel(2),2);
+
+  SetY(c.GetY());
+  SetZ(c.GetZ());
+  SetSigmaZ2(c.GetSigmaZ2());
+
+  for (Int_t i = 0; i < 7; i++) {
+    fSignals[i] = c.fSignals[i];
+  }
+
+  return *this;
+
+} 
+
+//_____________________________________________________________________________
+void AliTRDcluster::AddTrackIndex(const Int_t * const track)
 {
   //
   // Adds track index. Currently assumed that track is an array of
@@ -241,11 +306,11 @@ void AliTRDcluster::Clear(Option_t *)
   fQ = 0;
   fCenter = 0;
   for (Int_t i = 0; i < 3; i++) SetLabel(0,i);
-  SetX(0);
-  SetY(0);
-  SetZ(0);
-  SetSigmaY2(0);
-  SetSigmaZ2(0);
+  SetX(0.);
+  SetY(0.);
+  SetZ(0.);
+  AliCluster::SetSigmaY2(0.);
+  SetSigmaZ2(0.);
   SetVolumeId(0);
 }
 
@@ -268,8 +333,32 @@ Float_t AliTRDcluster::GetSumS() const
 //___________________________________________________________________________
 Double_t AliTRDcluster::GetSX(Int_t tb, Double_t z)
 {
+  //
+  // Returns the error parameterization in the radial direction for TRD clusters as function of 
+  // the calibrated time bin (tb) and optionally distance to anode wire (z). By default (no z information) 
+  // the mean value over all cluster to wire distance is chosen.
+  // 
+  // There are several contributions which are entering in the definition of the radial errors of the clusters. 
+  // Although an analytic defition should be possible for the moment this is not yet available but instead a 
+  // numerical parameterization is provided (see AliTRDclusterResolution::ProcessSigma() for the calibration 
+  // method). The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. 
+  //
+  //Begin_Html
+  //<img src="TRD/clusterXerrorDiff2D.gif">
+  //End_Html
+  //
+  // Here is a list of uncertainty components:
+  // - Time Response Function (TRF) - the major contribution. since TRF is also not symmetric (even if tail is 
+  //   cancelled) it also creates a systematic shift dependent on the charge distribution before and after the cluster.
+  // - longitudinal diffusion - increase the width of TRF and scales with square root of drift length
+  // - variation in the drift velocity within the drift cell 
+  //
+  // Author
+  // A.Bercuci <A.Bercuci@gsi.de>
+  //
+
   if(tb<1 || tb>=24) return 10.; // return huge [10cm]
-  const Double_t sx[24][10]={
+  static const Double_t sx[24][10]={
     {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},
     {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},
     {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},
@@ -297,15 +386,76 @@ Double_t AliTRDcluster::GetSX(Int_t tb, Double_t z)
   };
   if(z>=0. && z<.25) return sx[tb][Int_t(z/.025)];
   
-  Double_t m = 1.e-8; for(Int_t id=10; id--;) if(sx[tb][id]>m) m=sx[tb][id];
-  return m;
+  Double_t m = 0.; for(Int_t id=10; id--;) m+=sx[tb][id];
+  return m*.1;
+
 }
 
 //___________________________________________________________________________
-Double_t AliTRDcluster::GetSY(Int_t tb, Double_t z)
+Double_t AliTRDcluster::GetSYdrift(Int_t tb, Int_t ly, Double_t/* z*/)
 {
+  //
+  // Returns the error parameterization for TRD clusters as function of the drift length (here calibrated time bin tb)
+  // and optionally distance to anode wire (z) for the LUT r-phi cluster shape. By default (no z information) the largest 
+  // value over all cluster to wire values is chosen.
+  // 
+  // For the LUT method the dependence of s_y with x and d is obtained via a fit to the cluster to MC 
+  // resolution. (see class AliTRDclusterResolution for more details). A normalization to the reference radial position
+  // x0 = 0.675 (tb=5 for ideal vd) is also applied (see GetSYprf()). The function is *NOT* calibration aware ! 
+  // The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. A comparison 
+  // with the GAUS parameterization is also given
+  //
+  // For the GAUS method the dependence of s_y with x is *analytic* and it is expressed by the relation.
+  // BEGIN_LATEX
+  // #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}}
+  // END_LATEX
+  // The result is displayed in the figure below as function of the drift time and compared with the LUT parameterization.
+  //Begin_Html
+  //<img src="TRD/clusterYerrorDiff2D.gif">
+  //<img src="TRD/clusterYerrorDiff1D.gif">
+  //End_Html
+  //
+  // Author
+  // A.Bercuci <A.Bercuci@gsi.de>
+  //
+
   if(tb<1 || tb>=24) return 10.; // return huge [10cm]
-  const Double_t sy[24][10]={
+  static const Float_t lSy[6][24] = {
+    {75.7561, 0.0325, 0.0175, 0.0174, 0.0206, 0.0232,
+     0.0253, 0.0262, 0.0265, 0.0264, 0.0266, 0.0257,
+     0.0258, 0.0261, 0.0259, 0.0253, 0.0257, 0.0261,
+     0.0255, 0.0250, 0.0259, 0.0266, 0.0278, 0.0319
+    },
+    {49.2252, 0.0371, 0.0204, 0.0189, 0.0230, 0.0261,
+     0.0281, 0.0290, 0.0292, 0.0286, 0.0277, 0.0279,
+     0.0285, 0.0281, 0.0291, 0.0281, 0.0281, 0.0282,
+     0.0272, 0.0282, 0.0282, 0.0284, 0.0310, 0.0334
+    },
+    {55.1674, 0.0388, 0.0212, 0.0200, 0.0239, 0.0271,
+     0.0288, 0.0299, 0.0306, 0.0300, 0.0296, 0.0303,
+     0.0293, 0.0290, 0.0291, 0.0294, 0.0295, 0.0290,
+     0.0293, 0.0292, 0.0292, 0.0293, 0.0316, 0.0358
+    },
+    {45.1004, 0.0411, 0.0225, 0.0215, 0.0249, 0.0281,
+     0.0301, 0.0315, 0.0320, 0.0308, 0.0318, 0.0321,
+     0.0312, 0.0311, 0.0316, 0.0315, 0.0310, 0.0308,
+     0.0313, 0.0303, 0.0314, 0.0314, 0.0324, 0.0369
+    },
+    {43.8614, 0.0420, 0.0239, 0.0224, 0.0268, 0.0296,
+     0.0322, 0.0336, 0.0333, 0.0326, 0.0321, 0.0325,
+     0.0329, 0.0326, 0.0323, 0.0322, 0.0326, 0.0320,
+     0.0329, 0.0319, 0.0314, 0.0329, 0.0341, 0.0373
+    },
+    {40.5440, 0.0434, 0.0246, 0.0236, 0.0275, 0.0311,
+     0.0332, 0.0345, 0.0347, 0.0347, 0.0340, 0.0336,
+     0.0339, 0.0344, 0.0339, 0.0341, 0.0341, 0.0342,
+     0.0345, 0.0328, 0.0341, 0.0332, 0.0356, 0.0398
+    },
+  };
+  // adjusted ...
+  return TMath::Max(lSy[ly][tb]-0.0150, 0.0010);
+
+/*  const Double_t sy[24][10]={
     {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},
     {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},
     {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},
@@ -331,23 +481,95 @@ Double_t AliTRDcluster::GetSY(Int_t tb, Double_t z)
     {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},
     {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}
   };
-  if(z>=0. && z<.25) return sy[tb][Int_t(z/.025)];
+  if(z>=0. && z<.25) return sy[tb][Int_t(z/.025)] -  sy[5][Int_t(z/.025)];
 
-  Double_t m = 1.e-8; for(Int_t id=10; id--;) if(sy[tb][id]>m) m=sy[tb][id];
+  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];
 
-  return m;
+  return m;*/
 }
 
+//___________________________________________________________________________
+Double_t AliTRDcluster::GetSYcharge(Float_t q)
+{
+  //
+  // Parameterization of the r-phi resolution component due to cluster charge.
+  // The value is the offset from the nominal cluster resolution defined as the
+  // cluster resolution at average cluster charge (q0).
+  // 
+  // BEGIN_LATEX
+  // #Delta #sigma_{y}(q) = a*(#frac{1}{q} - #frac{1}{q_{0}})
+  // q_{0} #approx 50
+  // END_LATEX
+  // The definition is *NOT* robust against gain fluctuations and thus two approaches are possible
+  // when residual miscalibration are available:
+  //   - determine parameterization with a resolution matching those of the gain
+  //   - define an analytic model which scales with the gain.
+  // 
+  // For more details please see AliTRDclusterResolution::ProcessCharge()
+  //
+  //Begin_Html
+  //<img src="TRD/clusterQerror.gif">
+  //End_Html
+  //
+  // Author
+  // A.Bercuci <A.Bercuci@gsi.de>
+  //
+
+  const Float_t sq0inv = 0.019962; // [1/q0]
+  const Float_t sqb    = 0.037328; // [cm]
+
+  return sqb*(1./q - sq0inv);
+}
+
+//___________________________________________________________________________
+Double_t AliTRDcluster::GetSYprf(Int_t ly, Double_t center, Double_t s2)
+{
+  //
+  // Parameterization of the cluster error in the r-phi direction due to charge sharing between 
+  // adiacent pads. Should be identical to what is provided in the OCDB as PRF [TODO]
+  //
+  // The parameterization is obtained from fitting cluster resolution at phi=exb and |x-0.675|<0.225. 
+  // For more details see AliTRDclusterResolution::ProcessCenter().
+  //
+  //Begin_Html
+  //<img src="TRD/clusterPRFerror.gif">
+  //End_Html
+  //
+  // Author
+  // A.Bercuci <A.Bercuci@gsi.de>
+  //
+
+/*  const Float_t scy[AliTRDgeometry::kNlayer][4] = {
+    {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
+    {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
+    {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
+    {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
+    {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
+    {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
+  };*/
+  const Float_t lPRF[] = {0.438, 0.403, 0.393, 0.382, 0.376, 0.345};
+
+  return s2*TMath::Gaus(center, 0., lPRF[ly]);
+}
+
+
 //___________________________________________________________________________
 Double_t AliTRDcluster::GetXcorr(Int_t tb, Double_t z)
 {
-  // drift length correction [cm]
-  // TODO to be parametrized in term of drift velocity
+  //
+  // Drift length correction [cm]. Due to variation of mean drift velocity along the drift region
+  // from nominal vd at xd->infinity. For drift velocity determination based on tracking information 
+  // the correction should be negligible.
+  //Begin_Html
+  //<img src="TRD/clusterXcorr.gif">
+  //End_Html
+  // TODO to be parametrized in term of drift velocity at infinite drift length
   // A.Bercuci (Mar 28 2009)
+  //
 
   if(tb<0 || tb>=24) return 0.;
   const Int_t nd = 5;
-  const Double_t dx[24][nd]={
+  static const Double_t dx[24][nd]={
     {+1.747e-01,+3.195e-01,+1.641e-01,+1.607e-01,+6.002e-01},
     {+5.468e-02,+5.760e-02,+6.365e-02,+8.003e-02,+1.067e-01},
     {-6.327e-02,-6.339e-02,-6.423e-02,-6.900e-02,-7.949e-02},
@@ -408,9 +630,14 @@ Double_t AliTRDcluster::GetXcorr(Int_t tb, Double_t z)
 //___________________________________________________________________________
 Double_t AliTRDcluster::GetYcorr(Int_t ly, Float_t y)
 {
-// PRF correction TODO to be replaced by the gaussian 
-// approximation with full error parametrization
-  const Float_t cy[AliTRDgeometry::kNlayer][3] = {
+  //
+  // PRF correction for the LUT r-phi cluster shape.
+  //Begin_Html
+  //<img src="TRD/clusterYcorr.gif">
+  //End_Html
+  //
+
+  static const Float_t cy[AliTRDgeometry::kNlayer][3] = {
     { 4.014e-04, 8.605e-03, -6.880e+00},
     {-3.061e-04, 9.663e-03, -6.789e+00},
     { 1.124e-03, 1.105e-02, -6.825e+00},
@@ -423,51 +650,61 @@ Double_t AliTRDcluster::GetYcorr(Int_t ly, Float_t y)
 }
 
 //_____________________________________________________________________________
-Float_t AliTRDcluster::GetXloc(Double_t t0, Double_t vd, Double_t *const /*q*/, Double_t *const /*xq*/, Double_t /*z*/)
+Float_t AliTRDcluster::GetXloc(Double_t t0, Double_t vd
+                             , const Double_t *const /*q*/
+                             , const Double_t *const /*xq*/
+                             , Double_t /*z*/)
 {
-//
-// (Re)Calculate cluster position in the x direction in local chamber coordinates (with respect to the anode wire 
-// position) using all available information from tracking.
-// Input parameters:
-//   t0 - calibration aware trigger delay [us]
-//   vd - drift velocity in the region of the cluster [cm/us]
-//   z  - distance to the anode wire [cm]. By default 0.2 !!
-//   q & xq - array of charges and cluster positions from previous clusters in the tracklet [a.u.]
-// Output values :
-//   return x position of the cluster with respect to the 
-//   anode wire using all tracking information
-//
-// The estimation of the radial position is based on calculating the drift time and the drift velocity at the point of 
-// estimation. The drift time can be estimated according to the expression:
-// BEGIN_LATEX
-// t_{drift} = t_{bin} - t_{0} - t_{cause}(x) - t_{TC}(q_{i-1}, q_{i-2}, ...)
-// END_LATEX
-// where t_0 is the delay of the trigger signal. t_cause is the causality delay between ionisation electrons hitting 
-// the anode and the registration of maximum signal by the electronics - it is due to the rising time of the TRF 
-// convoluted with the diffusion width. t_TC is the residual charge from previous bins due to residual tails after tail 
-// cancellation.
-//
-// The drift velocity is considered to vary linearly with the drift length (independent of the distance to the anode wire 
-// in the z direction). Thus one can write the calculate iteratively the drift length from the expression:
-// BEGIN_LATEX
-// x = t_{drift}(x)*v_{drfit}(x)
-// END_LATEX
-//
-// Authors
-// Alex Bercuci <A.Bercuci@gsi.de>
-//
+  //
+  // (Re)Calculate cluster position in the x direction in local chamber coordinates (with respect to the anode wire 
+  // position) using all available information from tracking.
+  // Input parameters:
+  //   t0 - calibration aware trigger delay [us]
+  //   vd - drift velocity in the region of the cluster [cm/us]
+  //   z  - distance to the anode wire [cm]. By default average over the drift cell width.
+  //   q & xq - array of charges and cluster positions from previous clusters in the tracklet [a.u.]
+  // Output values :
+  //   return x position of the cluster with respect to the 
+  //   anode wire using all tracking information
+  //
+  // The estimation of the radial position is based on calculating the drift time and the drift velocity at the point of 
+  // estimation. The drift time can be estimated according to the expression:
+  // BEGIN_LATEX
+  // t_{drift} = t_{bin} - t_{0} - t_{cause}(x) - t_{TC}(q_{i-1}, q_{i-2}, ...)
+  // END_LATEX
+  // where t_0 is the delay of the trigger signal. t_cause is the causality delay between ionisation electrons hitting 
+  // the anode and the registration of maximum signal by the electronics - it is due to the rising time of the TRF 
+  // A second order correction here comes from the fact that the time spreading of charge at anode is the convolution of
+  // TRF with the diffusion and thus cross-talk between clusters before and after local clusters changes with drift length. 
+  // t_TC is the residual charge from previous (in time) clusters due to residual tails after tail cancellation. 
+  // This tends to push cluster forward and depends on the magnitude of their charge.
+  //
+  // The drift velocity varies with the drift length (and distance to anode wire) as described by cell structure simulation. 
+  // Thus one, in principle, can calculate iteratively the drift length from the expression:
+  // BEGIN_LATEX
+  // x = t_{drift}(x)*v_{drift}(x)
+  // END_LATEX
+  // In practice we use a numerical approach (AliTRDcluster::GetXcorr()) to correct for anisochronity obtained from MC 
+  // comparison (see AliTRDclusterResolution::ProcessSigma()). Also the calibration of 0 approximation (no x dependence)
+  // for t_cause is obtained from MC comparisons and impossible to disentangle in real life from trigger delay.
+  //
+  // Author
+  // Alex Bercuci <A.Bercuci@gsi.de>
+  //
 
   AliTRDCommonParam *cp = AliTRDCommonParam::Instance(); 
   Double_t fFreq = cp->GetSamplingFrequency();
 
-  // calculate t0 corrected time bin
-  Double_t td = fPadTime - t0;
-  fLocalTimeBin = TMath::Nint(td);
   //drift time corresponding to the center of the time bin
-  td = (td + .5)/fFreq; // [us] 
+  Double_t td = (fPadTime + .5)/fFreq; // [us] 
   // correction for t0
   td -= t0;
-  // calculate radial posion of clusters in the drift region
+  // time bin corrected for t0
+  // Bug in TMath::Nint().root-5.23.02
+  // TMath::Nint(3.5) = 4 and TMath::Nint(4.5) = 4
+  Double_t tmp = td*fFreq;
+  fLocalTimeBin = Char_t(TMath::Floor(tmp));
+  if(tmp-fLocalTimeBin > .5) fLocalTimeBin++;
   if(td < .2) return 0.;
   // TRF rising time (fitted)
   // It should be absorbed by the t0. For the moment t0 is 0 for simulations.
@@ -475,45 +712,119 @@ Float_t AliTRDcluster::GetXloc(Double_t t0, Double_t vd, Double_t *const /*q*/,
   td -= 0.189;
 
   // apply fitted correction 
-  Float_t x = td*vd + GetXcorr(fLocalTimeBin);
+  Float_t x = td*vd + (HasXcorr() ? GetXcorr(fLocalTimeBin) : 0.);
   if(x>0.&&x<.5*AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()) SetInChamber();
 
   return x;
-
-/*
-  // invert drift time function
-  Double_t xM= AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght(),
-           x = vd*td + .5*AliTRDgeometry::CamHght(), 
-           t = cp->TimeStruct(vd, x, z), dx1=0.,dx2;
-  while(TMath::Abs(td-t)>1.e-4){ // convergence on 100ps
-    dx2 = vd*(td-t);
-    if(TMath::Abs(TMath::Abs(dx2)-TMath::Abs(dx1))<1.e-6){
-      x+=.5*dx2;
-      break;
-    } else x+=dx2;
-
-    if(x<0. || x>xM) return 0.;
-    t = cp->TimeStruct(vd, x, z);
-    dx1 = dx2;
-  }
-
-  return x-.5*AliTRDgeometry::CamHght();
-*/
 }
 
 //_____________________________________________________________________________
-Float_t AliTRDcluster::GetYloc(Double_t y0, Double_t s2, Double_t W, Double_t *const y1, Double_t *const y2)
+Float_t AliTRDcluster::GetYloc(Double_t y0, Double_t s2, Double_t W
+                             , Double_t *const y1, Double_t *const y2)
 {
+  //
+  // Calculate, in tracking cooordinate system, the r-phi offset the cluster
+  // from the middle of the center pad. Three possible methods are implemented:
+  //   - Center of Gravity (COG) see AliTRDcluster::GetDYcog()
+  //   - Look-up Table (LUT) see AliTRDcluster::GetDYlut()
+  //   - Gauss shape (GAUS) see AliTRDcluster::GetDYgauss()
+  // In addition for the case of LUT method position corrections are also
+  // applied (see AliTRDcluster::GetYcorr())
+  //
 
-  //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);
   if(IsRPhiMethod(kCOG)) GetDYcog();
   else if(IsRPhiMethod(kLUT)) GetDYlut();
-  if(IsRPhiMethod(kGAUS)) GetDYgauss(s2/W/W, y1, y2);
+  else if(IsRPhiMethod(kGAUS)) GetDYgauss(s2/W/W, y1, y2);
+  else return 0.;
 
   if(y1) (*y1)*=W;
   if(y2) (*y2)*=W;
 
-  return y0-fCenter*W;
+  return y0+fCenter*W+(IsRPhiMethod(kLUT)?GetYcorr(AliTRDgeometry::GetLayer(fDetector), fCenter):0.);
+}
+
+//___________________________________________________________________________
+void AliTRDcluster::SetSigmaY2(Float_t s2, Float_t dt, Float_t exb, Float_t x, Float_t z, Float_t tgp)
+{
+  //
+  // Set variance of TRD cluster in the r-phi direction for each method.
+  // Parameters :
+  //   - s2  - variance due to PRF width for the case of Gauss model. Replaced by parameterization in case of LUT.
+  //   - dt  - transversal diffusion coeficient 
+  //   - exb - tg of lorentz angle
+  //   - x   - drift length - with respect to the anode wire
+  //   - z   - offset from the anode wire
+  //   - tgp - local tangent of the track momentum azimuthal angle
+  //
+  // The ingredients from which the error is computed are:
+  //   - PRF (charge sharing on adjacent pads)  - see AliTRDcluster::GetSYprf()
+  //   - diffusion (dependence with drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSYdrift()
+  //   - charge of the cluster (complex dependence on gain and tail cancellation) - see AliTRDcluster::GetSYcharge()
+  //   - lorentz angle (dependence on the drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSX()
+  //   - track angle (superposition of charges on the anode wire) - see AliTRDseedV1::Fit()
+  //   - projection of radial(x) error on r-phi due to fixed value assumed in tracking for x - see AliTRDseedV1::Fit()
+  //
+  // The last 2 contributions to cluster error can be estimated only during tracking when the track angle 
+  // is known (tgp). For this reason the errors (and optional position) of TRD clusters are recalculated during 
+  // tracking and thus clusters attached to tracks might differ from bare clusters.
+  // 
+  // Taking into account all contributions one can write the the TRD cluster error parameterization as:
+  // BEGIN_LATEX
+  // #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
+  // END_LATEX
+  // From this formula one can deduce a that the simplest calibration method for PRF and diffusion contributions is 
+  // by measuring resolution at B=0T and phi=0. To disentangle further the two remaining contributions one has 
+  // to represent s2 as a function of drift length. 
+  // 
+  // In the gaussian model the diffusion contribution can be expressed as:
+  // BEGIN_LATEX
+  // #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}}
+  // END_LATEX
+  // thus resulting the PRF contribution. For the case of the LUT model both contributions have to be determined from 
+  // the fit (see AliTRDclusterResolution::ProcessCenter() for details).
+  // 
+  // Author:
+  // A.Bercuci <A.Bercuci@gsi.de>
+  //
+
+  Float_t sigmaY2 = 0.;
+  Int_t ly = AliTRDgeometry::GetLayer(fDetector);
+  if(IsRPhiMethod(kCOG)) sigmaY2 = 4.e-4;
+  else if(IsRPhiMethod(kLUT)){ 
+    Float_t sd = GetSYdrift(fLocalTimeBin, ly, z); //printf("drift[%6.2f] ", 1.e4*sd);
+    sigmaY2 = GetSYprf(ly, fCenter, sd); //printf("PRF[%6.2f] ", 1.e4*sigmaY2);
+    // add charge contribution TODO scale with respect to s2
+    sigmaY2+= GetSYcharge(TMath::Abs(fQ)); //printf("Q[%6.2f] ", 1.e4*sigmaY2);
+    sigmaY2 = TMath::Max(sigmaY2, Float_t(0.0010)); //!! protection 
+    sigmaY2*= sigmaY2;
+  } else if(IsRPhiMethod(kGAUS)){
+    // PRF contribution
+    sigmaY2 = s2;
+    // Diffusion contribution
+    Double_t sD2 = dt/(1.+exb); sD2 *= sD2; sD2 *= x;
+    sigmaY2+= sD2; 
+    // add charge contribution TODO scale with respect to s2
+    //sigmaY2+= GetSYcharge(TMath::Abs(fQ));
+  }
+
+  // store tg^2(phi-a_L) and tg^2(a_L)
+  // limit parametrization to a maximum angle of 25 deg
+  if(TMath::Abs(tgp)>0.466) tgp = (tgp>0.)?0.466:-0.466;
+  Double_t tgg = (tgp-exb)/(1.+tgp*exb); tgg *= tgg;
+  Double_t exb2= exb*exb;
+
+  // Lorentz angle shift contribution 
+  Float_t sx = GetSX(fLocalTimeBin, z); sx*=sx;
+  sigmaY2+= exb2*sx; //printf("Al[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2));
+
+  // Radial contribution due to not measuring x in Kalman model 
+  sigmaY2+= tgg*sx; //printf("x[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2));
+
+  // Track angle contribution
+  sigmaY2+= tgg*x*x*exb2/12.; //printf("angle[%6.2f]\n", 1.e4*TMath::Sqrt(sigmaY2));
+
+  AliCluster::SetSigmaY2(sigmaY2);
+
 }
 
 //_____________________________________________________________________________
@@ -540,58 +851,84 @@ Bool_t AliTRDcluster::IsEqual(const TObject *o) const
   if ( IsUsed() != inCluster->IsUsed() ) return kFALSE;
   
   return kTRUE;
+
 }
 
 //_____________________________________________________________________________
 void AliTRDcluster::Print(Option_t *o) const
 {
-  AliInfo(Form("Det[%3d] LTrC[%+6.2f %+6.2f %+6.2f] Q[%5.1f] Stat[in(%c) use(%c) sh(%c)]", 
-    fDetector, GetX(), GetY(), GetZ(), fQ, 
-    IsInChamber() ? 'y' : 'n', IsUsed() ? 'y' : 'n', IsShared() ? 'y' : 'n'));
+  //
+  // Print cluster information
+  //
 
-  if(strcmp(o, "a")!=0) return;
-  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)); 
-  AliInfo(Form("Signals[%3d %3d %3d %3d %3d %3d %3d]", fSignals[0], fSignals[1], fSignals[2], fSignals[3], fSignals[4], fSignals[5], fSignals[6]));
+  if(strcmp(o, "a")==0) {
+    AliInfo(Form(
+    "\nDet[%3d] LTrC[%+6.2f %+6.2f %+6.2f] Q[%5.1f] FLAG[in(%c) use(%c) sh(%c)] Y[%s]"
+    "\n         LChC[c(%3d) r(%2d) t(%2d)] t-t0[%2d] Npad[%d] cen[%5.3f] mask[%d]"
+    "\n         QS[%3d %3d %3d %3d %3d %3d %3d] S2[%e %e]"
+    , fDetector, GetX(), GetY(), GetZ(), fQ, 
+    IsInChamber() ? 'y' : 'n', 
+    IsUsed() ? 'y' : 'n', 
+    IsShared() ? 'y' : 'n',
+    IsRPhiMethod(kGAUS)?"GAUS":(IsRPhiMethod(kLUT)?"LUT":"COG")
+    , fPadCol, fPadRow, fPadTime, fLocalTimeBin, fNPads, fCenter, fClusterMasking
+    , fSignals[0], fSignals[1], fSignals[2], fSignals[3]
+    , fSignals[4], fSignals[5], fSignals[6]
+    , GetSigmaY2(), GetSigmaZ2()));
+  } else { 
+    AliInfo(Form("Det[%3d] LTrC[%+6.2f %+6.2f %+6.2f] Q[%5.1f] FLAG[in(%c) use(%c) sh(%c)] Y[%s]", 
+    fDetector, GetX(), GetY(), GetZ(), fQ, 
+    IsInChamber() ? 'y' : 'n', 
+    IsUsed() ? 'y' : 'n', 
+    IsShared() ? 'y' : 'n',
+    IsRPhiMethod(kGAUS)?"GAUS":(IsRPhiMethod(kLUT)?"LUT":"COG")
+    ));
+  }
 }
 
-
 //_____________________________________________________________________________
 void AliTRDcluster::SetPadMaskedPosition(UChar_t position)
 {
   //
-  // store the pad corruption position code
+  // Store the pad corruption position code
   // 
   // Code: 1 = left cluster
   //       2 = middle cluster;
   //       4 = right cluster
   //
-  for(Int_t ipos = 0; ipos < 3; ipos++)
-    if(TESTBIT(position, ipos))
+
+  for (Int_t ipos = 0; ipos < 3; ipos++) {
+    if (TESTBIT(position, ipos))
       SETBIT(fClusterMasking, ipos);
+  }
 }
 
 //_____________________________________________________________________________
 void AliTRDcluster::SetPadMaskedStatus(UChar_t status)
 {
   //
-  // store the status of the corrupted pad
+  // Store the status of the corrupted pad
   //
   // Code: 2 = noisy
   //       4 = Bridged Left
   //       8 = Bridged Right
   //      32 = Not Connected
-  for(Int_t ipos = 0; ipos < 5; ipos++)
+  //
+
+  for (Int_t ipos = 0; ipos < 5; ipos++) {
     if(TESTBIT(status, ipos))
       SETBIT(fClusterMasking, ipos + 3);
+  }
+
 }
 
 //___________________________________________________________________________
-Float_t AliTRDcluster::GetDYcog(Double_t *const, Double_t *const)
+Float_t AliTRDcluster::GetDYcog(const Double_t *const, const Double_t *const)
 {
-//
-// Get COG position
-// Used for clusters with more than 3 pads - where LUT not applicable
-//
+  //
+  // Get COG position
+  // Used for clusters with more than 3 pads - where LUT not applicable
+  //
 
   Double_t sum = fSignals[1]
                 +fSignals[2]
@@ -603,13 +940,13 @@ Float_t AliTRDcluster::GetDYcog(Double_t *const, Double_t *const)
   // Go to 3 pad COG ????
   // ???????????? CBL
   fCenter = (0.0 * (-fSignals[1] + fSignals[5])
-                      + (-fSignals[2] + fSignals[4])) / sum;
+                    + (-fSignals[2] + fSignals[4])) / sum;
 
   return fCenter;
 }
 
 //___________________________________________________________________________
-Float_t AliTRDcluster::GetDYlut(Double_t *const, Double_t *const)
+Float_t AliTRDcluster::GetDYlut(const Double_t *const, const Double_t *const)
 {
   //
   // Calculates the cluster position using the lookup table.
@@ -697,8 +1034,15 @@ Float_t AliTRDcluster::GetDYgauss(Double_t s2w, Double_t *const y1, Double_t *co
 // Alex Bercuci <A.Bercuci@gsi.de>
 // Theodor Rascanu <trascanu@stud.uni-frankfurt.de>
 //
-  Double_t w1  = fSignals[2]*fSignals[2];
-  Double_t w2  = fSignals[4]*fSignals[4];
+  Double_t w1 = fSignals[2]*fSignals[2];
+  Double_t w2 = fSignals[4]*fSignals[4];
+  Double_t w = w1+w2;
+  if(w<1.){
+    AliError("Missing side signals for cluster.");
+    Print("a");
+    return 0.;
+  }  
+
   //Double_t s2w = s2/W/W;
   Float_t y1r  = fSignals[2]>0 ? (-0.5 + s2w*TMath::Log(fSignals[3]/(Float_t)fSignals[2])) : 0.;
   Float_t y2r  = fSignals[4]>0 ? (0.5 + s2w*TMath::Log(fSignals[4]/(Float_t)fSignals[3])) : 0.;
@@ -706,7 +1050,7 @@ Float_t AliTRDcluster::GetDYgauss(Double_t s2w, Double_t *const y1, Double_t *co
   if(y1) (*y1) = y1r;
   if(y2) (*y2) = y2r;
 
-  return fCenter      = (w1*y1r+w2*y2r)/(w1+w2);
+  return fCenter      = (w1*y1r+w2*y2r)/w;
 }