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 //___________________________________________________________________________
  35 AliTRDcluster::AliTRDcluster()
  36   :AliCluster()
  37   ,fPadCol(0)
  38   ,fPadRow(0)
  39   ,fPadTime(0)
  40   ,fLocalTimeBin(0)
  41   ,fNPads(0)
  42   ,fClusterMasking(0)
  43   ,fDetector(0)
  44   ,fQ(0)
  45   ,fCenter(0)
  46 {
  47   //
  48   // Default constructor
  49   //
  50
  51   for (Int_t i = 0; i < 7; i++) {
  52     fSignals[i] = 0;
  53   }
  54
  55 }
  56
  57 //___________________________________________________________________________
  58 AliTRDcluster::AliTRDcluster(Int_t det, Float_t q
  59                            , Float_t *pos, Float_t *sig
  60                            , Int_t *tracks, Char_t npads, Short_t *signals
  61                            , UChar_t col, UChar_t row, UChar_t time
  62                            , Char_t timebin, Float_t center, UShort_t volid)
  63   :AliCluster(volid,pos[0],pos[1],pos[2],sig[0],sig[1],0.0,0x0)
  64   ,fPadCol(col)
  65   ,fPadRow(row)
  66   ,fPadTime(time)
  67   ,fLocalTimeBin(timebin)
  68   ,fNPads(npads)
  69   ,fClusterMasking(0)
  70   ,fDetector(det)
  71   ,fQ(q)
  72   ,fCenter(center)
  73 {
  74   //
  75   // Constructor
  76   //
  77
  78   for (Int_t i = 0; i < 7; i++) {
  79     fSignals[i] = signals[i];
  80   }
  81
  82   if (tracks) {
  83     AddTrackIndex(tracks);
  84   }
  85
  86 }
  87
  88 //_____________________________________________________________________________
  89 AliTRDcluster::AliTRDcluster(const AliTRDcluster &c)
  90   :AliCluster(c)
  91   ,fPadCol(c.fPadCol)
  92   ,fPadRow(c.fPadRow)
  93   ,fPadTime(c.fPadTime)
  94   ,fLocalTimeBin(c.fLocalTimeBin)
  95   ,fNPads(c.fNPads)
  96   ,fClusterMasking(c.fClusterMasking)
  97   ,fDetector(c.fDetector)
  98   ,fQ(c.fQ)
  99   ,fCenter(c.fCenter)
 100 {
 101   //
 102   // Copy constructor
 103   //
 104
 105   SetBit(kInChamber, c.IsInChamber());
 106   SetLabel(c.GetLabel(0),0);
 107   SetLabel(c.GetLabel(1),1);
 108   SetLabel(c.GetLabel(2),2);
 109
 110   SetY(c.GetY());
 111   SetZ(c.GetZ());
 112   SetSigmaY2(c.GetSigmaY2());
 113   SetSigmaZ2(c.GetSigmaZ2());
 114
 115   for (Int_t i = 0; i < 7; i++) {
 116     fSignals[i] = c.fSignals[i];
 117   }
 118
 119 }
 120
 121 //_____________________________________________________________________________
 122 void AliTRDcluster::AddTrackIndex(Int_t *track)
 123 {
 124   //
 125   // Adds track index. Currently assumed that track is an array of
 126   // size 9, and up to 3 track indexes are stored in fTracks[3].
 127   // Indexes are sorted according to:
 128   //  1) index of max number of appearances is stored first
 129   //  2) if two or more indexes appear equal number of times, the lowest
 130   //     ones are stored first;
 131   //
 132
 133   const Int_t kSize = 9;
 134   Int_t  entries[kSize][2];
 135
 136   Int_t  i = 0;
 137   Int_t  j = 0;
 138   Int_t  k = 0;
 139   Int_t  index;
 140   Bool_t indexAdded;
 141
 142   for (i = 0; i < kSize; i++) {
 143     entries[i][0] = -1;
 144     entries[i][1] =  0;
 145   }
 146
 147   for (k = 0; k < kSize; k++) {
 148
 149     index      = track[k];
 150     indexAdded = kFALSE;
 151
 152     j = 0;
 153     if (index >= 0) {
 154       while ((!indexAdded) && (j < kSize)) {
 155         if ((entries[j][0] == index) ||
 156             (entries[j][1] ==     0)) {
 157           entries[j][0] = index;
 158           entries[j][1] = entries[j][1] + 1;
 159           indexAdded    = kTRUE;
 160         }
 161         j++;
 162       }
 163     }
 164
 165   }
 166
 167   // Sort by number of appearances and index value
 168   Int_t swap = 1;
 169   Int_t tmp0;
 170   Int_t tmp1;
 171   while (swap > 0) {
 172     swap = 0;
 173     for (i = 0; i < (kSize - 1); i++) {
 174       if ((entries[i][0]   >= 0) &&
 175           (entries[i+1][0] >= 0)) {
 176         if ((entries[i][1] < entries[i+1][1]) ||
 177             ((entries[i][1] == entries[i+1][1]) &&
 178              (entries[i][0] >  entries[i+1][0]))) {
 179           tmp0            = entries[i][0];
 180           tmp1            = entries[i][1];
 181           entries[i][0]   = entries[i+1][0];
 182           entries[i][1]   = entries[i+1][1];
 183           entries[i+1][0] = tmp0;
 184           entries[i+1][1] = tmp1;
 185           swap++;
 186         }
 187       }
 188     }
 189   }
 190
 191   // Set track indexes
 192   for (i = 0; i < 3; i++) {
 193     SetLabel(entries[i][0],i);
 194   }
 195
 196   return;
 197
 198 }
 199
 200 //_____________________________________________________________________________
 201 void AliTRDcluster::Clear(Option_t *)
 202 {
 203   //
 204   // Reset all member to the default value
 205   //
 206   fPadCol=0;
 207   fPadRow=0;
 208   fPadTime=0;
 209   fLocalTimeBin=0;
 210   fNPads=0;
 211   fClusterMasking=0;
 212   fDetector=0;
 213   for (Int_t i=0; i < 7; i++) fSignals[i]=0;
 214   fQ = 0;
 215   fCenter = 0;
 216   for (Int_t i = 0; i < 3; i++) SetLabel(0,i);
 217   SetX(0);
 218   SetY(0);
 219   SetZ(0);
 220   SetSigmaY2(0);
 221   SetSigmaZ2(0);
 222   SetVolumeId(0);
 223 }
 224
 225 //_____________________________________________________________________________
 226 Float_t AliTRDcluster::GetSumS() const
 227 {
 228   //
 229   // Returns the total charge from a not unfolded cluster
 230   //
 231
 232   Float_t sum = 0.0;
 233   for (Int_t i = 0; i < 7; i++) {
 234     sum += fSignals[i];
 235   }
 236
 237   return sum;
 238
 239 }
 240
 241 //___________________________________________________________________________
 242 Float_t AliTRDcluster::GetXcorr(Int_t tb)
 243 {
 244   // drift length correction [cm]
 245   // TODO to be parametrized in term of drift velocity
 246   // A.Bercuci (Mar 28 2009)
 247
 248   if(tb<0 || tb>=24) return 0.;
 249   const Float_t cx[] = {
 250     1.6402e-01, 7.2917e-02,-6.7848e-02,-1.4529e-01,-1.6279e-01,-1.3116e-01,
 251    -8.2712e-02,-4.9453e-02,-2.9501e-02,-1.4543e-02,-6.1749e-03, 3.9221e-04,
 252     1.9711e-03, 2.7388e-03, 2.9070e-03, 3.4183e-03, 2.8014e-03, 1.9351e-03,
 253     4.9252e-04, 4.5742e-04, 1.2263e-04,-1.2219e-02,-6.9658e-02,-1.6681e-01};
 254   return cx[tb];
 255 }
 256
 257 //___________________________________________________________________________
 258 Float_t AliTRDcluster::GetYcorr(Int_t ly, Float_t y)
 259 {
 260   // PRF correction TODO to be replaced by the gaussian
 261   // approximation with full error parametrization and // moved to the clusterizer
 262   const Float_t cy[AliTRDgeometry::kNlayer][3] = {
 263     { 4.014e-04, 8.605e-03, -6.880e+00},
 264     {-3.061e-04, 9.663e-03, -6.789e+00},
 265     { 1.124e-03, 1.105e-02, -6.825e+00},
 266     {-1.527e-03, 1.231e-02, -6.777e+00},
 267     { 2.150e-03, 1.387e-02, -6.783e+00},
 268     {-1.296e-03, 1.486e-02, -6.825e+00}
 269   };
 270
 271   return cy[ly][0] + cy[ly][1] * TMath::Sin(cy[ly][2] * y);
 272 }
 273
 274 //_____________________________________________________________________________
 275 Float_t AliTRDcluster::GetXloc(Double_t t0, Double_t vd, Double_t *const /*q*/, Double_t *const /*xq*/, Double_t z)
 276 {
 277 //
 278 // (Re)Calculate cluster position in the x direction in local chamber coordinates (with respect to the anode wire
 279 // position) using all available information from tracking.
 280 // Input parameters:
 281 //   t0 - calibration aware trigger delay [us]
 282 //   vd - drift velocity in the region of the cluster [cm/us]
 283 //   z  - distance to the anode wire [cm]. By default 0.2 !!
 284 //   q & xq - array of charges and cluster positions from previous clusters in the tracklet [a.u.]
 285 // Output values :
 286 //   return x position of the cluster with respect to the
 287 //   anode wire using all tracking information
 288 //
 289 // The estimation of the radial position is based on calculating the drift time and the drift velocity at the point of
 290 // estimation. The drift time can be estimated according to the expression:
 291 // BEGIN_LATEX
 292 // t_{drift} = t_{bin} - t_{0} - t_{cause}(x) - t_{TC}(q_{i-1}, q_{i-2}, ...)
 293 // END_LATEX
 294 // where t_0 is the delay of the trigger signal. t_cause is the causality delay between ionisation electrons hitting
 295 // the anode and the registration of maximum signal by the electronics - it is due to the rising time of the TRF
 296 // convoluted with the diffusion width. t_TC is the residual charge from previous bins due to residual tails after tail
 297 // cancellation.
 298 //
 299 // The drift velocity is considered to vary linearly with the drift length (independent of the distance to the anode wire
 300 // in the z direction). Thus one can write the calculate iteratively the drift length from the expression:
 301 // BEGIN_LATEX
 302 // x = t_{drift}(x)*v_{drfit}(x)
 303 // END_LATEX
 304 //
 305 // Authors
 306 // Alex Bercuci <A.Bercuci@gsi.de>
 307 //
 308
 309   AliTRDCommonParam *cp = AliTRDCommonParam::Instance();
 310   Double_t fFreq = cp->GetSamplingFrequency();
 311   //drift time corresponding to the center of the time bin
 312   Double_t td = (fPadTime + .5)/fFreq; // [us]
 313   // correction for t0
 314   td -= t0;
 315   // calculate radial posion of clusters in the drift region
 316   if(td < .2 || td > 2.4) return 0.;
 317   // correction for TRF rising time 0.2us
 318   td -= 0.189;
 319
 320   // invert drift time function
 321   Double_t xM= AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght(),
 322            x = vd*td + .5*AliTRDgeometry::CamHght(),
 323            t = cp->TimeStruct(vd, x, z), dx1=0.,dx2;
 324   while(TMath::Abs(td-t)>1.e-4){ // convergence on 100ps
 325     dx2 = vd*(td-t);
 326     if(TMath::Abs(TMath::Abs(dx2)-TMath::Abs(dx1))<1.e-6){
 327       x+=.5*dx2;
 328       break;
 329     } else x+=dx2;
 330
 331     if(x<0. || x>xM) return 0.;
 332     t = cp->TimeStruct(vd, x, z);
 333     dx1 = dx2;
 334   }
 335
 336   return x-.5*AliTRDgeometry::CamHght();
 337 }
 338
 339 //_____________________________________________________________________________
 340 Float_t AliTRDcluster::GetYloc(Double_t s2, Double_t W, Double_t xd, Double_t wt, Double_t *const y1, Double_t *const y2)
 341 {
 342 //
 343 // (Re)Calculate cluster position in the y direction in local chamber coordinates using all available information from tracking.
 344 //
 345 // Input parameters:
 346 //   s2 - sigma of gaussian parameterization (see bellow for the exact parameterization)
 347 //   W  - pad width
 348 //   xd - drift length (with respect to the anode wire) [cm]
 349 //   wt - omega*tau = tg(a_L)
 350 // Output values :
 351 //   y1 and y2 - partial positions based on 2 pads clusters
 352 //   return y position of the cluster from all information
 353 //
 354 // Estimation of y coordinate is based on the gaussian approximation of the PRF. Thus one may
 355 // calculate the y position knowing the signals q_i-1, q_i and q_i+1 in the 3 adiacent pads by:
 356 // BEGIN_LATEX
 357 // 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}}}}
 358 // END_LATEX
 359 // where W is the pad width, y_0 is the position of the center pad and s^2 is given by
 360 // BEGIN_LATEX
 361 // s^{2} = s^{2}_{0} + s^{2}_{diff} (x,B) + #frac{tg^{2}(#phi-#alpha_{L})*l^{2}}{12}
 362 // END_LATEX
 363 // with s_0 being the PRF for 0 drift and track incidence phi equal to the lorentz angle a_L and the diffusion term
 364 // being described by:
 365 // BEGIN_LATEX
 366 // s_{diff} (x,B) = #frac{D_{L}#sqrt{x}}{1+#(){#omega#tau}^{2}}
 367 // END_LATEX
 368 // 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
 369 //
 370 // Authors
 371 // Alex Bercuci <A.Bercuci@gsi.de>
 372 // Theodor Rascanu <trascanu@stud.uni-frankfurt.de>
 373 //
 374   Float_t y0 = GetY()-W*fCenter;
 375   Double_t w1 = fSignals[2]*fSignals[2];
 376   Double_t w2 = fSignals[4]*fSignals[4];
 377   Float_t y1r = fSignals[2]>0 ? (y0 - .5*W + s2*TMath::Log(fSignals[3]/(Float_t)fSignals[2])/W) : 0.;
 378   Float_t y2r = fSignals[4]>0 ? (y0 + .5*W + s2*TMath::Log(fSignals[4]/(Float_t)fSignals[3])/W) : 0.;
 379
 380   if(y1) (*y1) = y1r;
 381   if(y2) (*y2) = y2r;
 382
 383   Double_t ld  = TMath::Max(xd - 0.*AliTRDgeometry::CamHght(), 0.);
 384
 385   return (w1*y1r+w2*y2r)/(w1+w2) - ld*wt;
 386 }
 387
 388 //_____________________________________________________________________________
 389 Bool_t AliTRDcluster::IsEqual(const TObject *o) const
 390 {
 391   //
 392   // Compare relevant information of this cluster with another one
 393   //
 394
 395   const AliTRDcluster *inCluster = dynamic_cast<const AliTRDcluster*>(o);
 396   if (!o || !inCluster) return kFALSE;
 397
 398   if ( AliCluster::GetX() != inCluster->GetX() ) return kFALSE;
 399   if ( AliCluster::GetY() != inCluster->GetY() ) return kFALSE;
 400   if ( AliCluster::GetZ() != inCluster->GetZ() ) return kFALSE;
 401   if ( fQ != inCluster->fQ ) return kFALSE;
 402   if ( fDetector != inCluster->fDetector ) return kFALSE;
 403   if ( fPadCol != inCluster->fPadCol ) return kFALSE;
 404   if ( fPadRow != inCluster->fPadRow ) return kFALSE;
 405   if ( fPadTime != inCluster->fPadTime ) return kFALSE;
 406   if ( fClusterMasking != inCluster->fClusterMasking ) return kFALSE;
 407   if ( IsInChamber() != inCluster->IsInChamber() ) return kFALSE;
 408   if ( IsShared() != inCluster->IsShared() ) return kFALSE;
 409   if ( IsUsed() != inCluster->IsUsed() ) return kFALSE;
 410
 411   return kTRUE;
 412 }
 413
 414 //_____________________________________________________________________________
 415 void AliTRDcluster::Print(Option_t *o) const
 416 {
 417   AliInfo(Form("Det[%3d] LTrC[%7.2f %7.2f %7.2f] Q[%f] Stat[in(%c) use(%c) sh(%c)]",
 418     fDetector, GetX(), GetY(), GetZ(), fQ,
 419     IsInChamber() ? 'y' : 'n', IsUsed() ? 'y' : 'n', IsShared() ? 'y' : 'n'));
 420
 421   if(strcmp(o, "a")!=0) return;
 422   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));
 423   AliInfo(Form("Signals[%3d %3d %3d %3d %3d %3d %3d]", fSignals[0], fSignals[1], fSignals[2], fSignals[3], fSignals[4], fSignals[5], fSignals[6]));
 424 }
 425
 426
 427 //_____________________________________________________________________________
 428 void AliTRDcluster::SetPadMaskedPosition(UChar_t position)
 429 {
 430   //
 431   // store the pad corruption position code
 432   //
 433   // Code: 1 = left cluster
 434   //       2 = middle cluster;
 435   //       4 = right cluster
 436   //
 437   for(Int_t ipos = 0; ipos < 3; ipos++)
 438     if(TESTBIT(position, ipos))
 439       SETBIT(fClusterMasking, ipos);
 440 }
 441
 442 //_____________________________________________________________________________
 443 void AliTRDcluster::SetPadMaskedStatus(UChar_t status)
 444 {
 445   //
 446   // store the status of the corrupted pad
 447   //
 448   // Code: 2 = noisy
 449   //       4 = Bridged Left
 450   //       8 = Bridged Right
 451   //      32 = Not Connected
 452   for(Int_t ipos = 0; ipos < 5; ipos++)
 453     if(TESTBIT(status, ipos))
 454       SETBIT(fClusterMasking, ipos + 3);
 455 }