TRD/AliTRDseedV1.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 //  The TRD track seed                                                    //
  21 //                                                                        //
  22 //  Authors:                                                              //
  23 //    Alex Bercuci <A.Bercuci@gsi.de>                                     //
  24 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  25 //                                                                        //
  26 ////////////////////////////////////////////////////////////////////////////
  27
  28 #include "TMath.h"
  29 #include "TLinearFitter.h"
  30 #include "TClonesArray.h" // tmp
  31 #include <TTreeStream.h>
  32
  33 #include "AliLog.h"
  34 #include "AliMathBase.h"
  35 #include "AliCDBManager.h"
  36 #include "AliTracker.h"
  37
  38 #include "AliTRDpadPlane.h"
  39 #include "AliTRDcluster.h"
  40 #include "AliTRDseedV1.h"
  41 #include "AliTRDtrackV1.h"
  42 #include "AliTRDcalibDB.h"
  43 #include "AliTRDchamberTimeBin.h"
  44 #include "AliTRDtrackingChamber.h"
  45 #include "AliTRDtrackerV1.h"
  46 #include "AliTRDReconstructor.h"
  47 #include "AliTRDrecoParam.h"
  48
  49 #include "Cal/AliTRDCalPID.h"
  50 #include "Cal/AliTRDCalROC.h"
  51 #include "Cal/AliTRDCalDet.h"
  52
  53 ClassImp(AliTRDseedV1)
  54
  55 //____________________________________________________________________
  56 AliTRDseedV1::AliTRDseedV1(Int_t det)
  57   :AliTRDseed()
  58   ,fReconstructor(0x0)
  59   ,fClusterIter(0x0)
  60   ,fClusterIdx(0)
  61   ,fDet(det)
  62   ,fMom(0.)
  63   ,fSnp(0.)
  64   ,fTgl(0.)
  65   ,fdX(0.)
  66   ,fXref(0.)
  67   ,fExB(0.)
  68 {
  69   //
  70   // Constructor
  71   //
  72   for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = 0.;
  73   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
  74   fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
  75   // covariance matrix [diagonal]
  76   // default sy = 200um and sz = 2.3 cm
  77   fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
  78 }
  79
  80 //____________________________________________________________________
  81 AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
  82   :AliTRDseed((AliTRDseed&)ref)
  83   ,fReconstructor(ref.fReconstructor)
  84   ,fClusterIter(0x0)
  85   ,fClusterIdx(0)
  86   ,fDet(ref.fDet)
  87   ,fMom(ref.fMom)
  88   ,fSnp(ref.fSnp)
  89   ,fTgl(ref.fTgl)
  90   ,fdX(ref.fdX)
  91   ,fXref(ref.fXref)
  92   ,fExB(ref.fExB)
  93 {
  94   //
  95   // Copy Constructor performing a deep copy
  96   //
  97
  98   //printf("AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &)\n");
  99   SetBit(kOwner, kFALSE);
 100   for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice];
 101   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec];
 102   memcpy(fRefCov, ref.fRefCov, 3*sizeof(Double_t));
 103   memcpy(fCov, ref.fCov, 3*sizeof(Double_t));
 104 }
 105
 106
 107 //____________________________________________________________________
 108 AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
 109 {
 110   //
 111   // Assignment Operator using the copy function
 112   //
 113
 114   if(this != &ref){
 115     ref.Copy(*this);
 116   }
 117   SetBit(kOwner, kFALSE);
 118
 119   return *this;
 120
 121 }
 122
 123 //____________________________________________________________________
 124 AliTRDseedV1::~AliTRDseedV1()
 125 {
 126   //
 127   // Destructor. The RecoParam object belongs to the underlying tracker.
 128   //
 129
 130   //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
 131
 132   if(IsOwner())
 133     for(int itb=0; itb<knTimebins; itb++){
 134       if(!fClusters[itb]) continue;
 135       //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
 136       delete fClusters[itb];
 137       fClusters[itb] = 0x0;
 138     }
 139 }
 140
 141 //____________________________________________________________________
 142 void AliTRDseedV1::Copy(TObject &ref) const
 143 {
 144   //
 145   // Copy function
 146   //
 147
 148   //AliInfo("");
 149   AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
 150
 151   target.fClusterIter   = 0x0;
 152   target.fClusterIdx    = 0;
 153   target.fDet           = fDet;
 154   target.fMom           = fMom;
 155   target.fSnp           = fSnp;
 156   target.fTgl           = fTgl;
 157   target.fdX            = fdX;
 158   target.fXref          = fXref;
 159   target.fExB           = fExB;
 160   target.fReconstructor = fReconstructor;
 161
 162   for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice];
 163   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec];
 164   memcpy(target.fRefCov, fRefCov, 3*sizeof(Double_t));
 165   memcpy(target.fCov, fCov, 3*sizeof(Double_t));
 166
 167   AliTRDseed::Copy(target);
 168 }
 169
 170
 171 //____________________________________________________________
 172 Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
 173 {
 174 // Initialize this tracklet using the track information
 175 //
 176 // Parameters:
 177 //   track - the TRD track used to initialize the tracklet
 178 //
 179 // Detailed description
 180 // The function sets the starting point and direction of the
 181 // tracklet according to the information from the TRD track.
 182 //
 183 // Caution
 184 // The TRD track has to be propagated to the beginning of the
 185 // chamber where the tracklet will be constructed
 186 //
 187
 188   Double_t y, z;
 189   if(!track->GetProlongation(fX0, y, z)) return kFALSE;
 190   fYref[0] = y;
 191   fYref[1] = track->GetSnp()/(1. - track->GetSnp()*track->GetSnp());
 192   fZref[0] = z;
 193   fZref[1] = track->GetTgl();
 194
 195   const Double_t *cov = track->GetCovariance();
 196   fRefCov[0] = cov[0]; // Var(y)
 197   fRefCov[1] = cov[1]; // Cov(yz)
 198   fRefCov[2] = cov[5]; // Var(z)
 199
 200   //printf("Tracklet ref x[%7.3f] y[%7.3f] z[%7.3f], snp[%f] tgl[%f]\n", fX0, fYref[0], fZref[0], track->GetSnp(), track->GetTgl());
 201   return kTRUE;
 202 }
 203
 204
 205 //____________________________________________________________________
 206 void AliTRDseedV1::CookdEdx(Int_t nslices)
 207 {
 208 // Calculates average dE/dx for all slices and store them in the internal array fdEdx.
 209 //
 210 // Parameters:
 211 //  nslices : number of slices for which dE/dx should be calculated
 212 // Output:
 213 //  store results in the internal array fdEdx. This can be accessed with the method
 214 //  AliTRDseedV1::GetdEdx()
 215 //
 216 // Detailed description
 217 // Calculates average dE/dx for all slices. Depending on the PID methode
 218 // the number of slices can be 3 (LQ) or 8(NN).
 219 // The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t)) i.e.
 220 //
 221 // dQ/dl = qc/(dx * sqrt(1 + dy/dx^2 + dz/dx^2))
 222 //
 223 // The following effects are included in the calculation:
 224 // 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
 225 // 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
 226 // 3. cluster size
 227 //
 228
 229   Int_t nclusters[knSlices];
 230   for(int i=0; i<knSlices; i++){
 231     fdEdx[i]     = 0.;
 232     nclusters[i] = 0;
 233   }
 234   Float_t clength = (/*.5 * */AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
 235
 236   AliTRDcluster *cluster = 0x0;
 237   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
 238     if(!(cluster = fClusters[ic])) continue;
 239     Float_t x = cluster->GetX();
 240
 241     // Filter clusters for dE/dx calculation
 242
 243     // 1.consider calibration effects for slice determination
 244     Int_t slice;
 245     if(cluster->IsInChamber()) slice = Int_t(TMath::Abs(fX0 - x) * nslices / clength);
 246     else slice = x < fX0 ? 0 : nslices-1;
 247
 248     // 2. take sharing into account
 249     Float_t w = cluster->IsShared() ? .5 : 1.;
 250
 251     // 3. take into account large clusters TODO
 252     //w *= c->GetNPads() > 3 ? .8 : 1.;
 253
 254     //CHECK !!!
 255     fdEdx[slice]   += w * GetdQdl(ic); //fdQdl[ic];
 256     nclusters[slice]++;
 257   } // End of loop over clusters
 258
 259   //if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
 260   if(nslices == AliTRDReconstructor::kLQslices){
 261   // calculate mean charge per slice (only LQ PID)
 262     for(int is=0; is<nslices; is++){
 263       if(nclusters[is]) fdEdx[is] /= nclusters[is];
 264     }
 265   }
 266 }
 267
 268 //____________________________________________________________________
 269 void AliTRDseedV1::GetClusterXY(const AliTRDcluster *c, Double_t &x, Double_t &y)
 270 {
 271 // Return corrected position of the cluster taking into
 272 // account variation of the drift velocity with drift length.
 273
 274
 275   // drift velocity correction TODO to be moved to the clusterizer
 276   const Float_t cx[] = {
 277     -9.6280e-02, 1.3091e-01,-1.7415e-02,-9.9221e-02,-1.2040e-01,-9.5493e-02,
 278     -5.0041e-02,-1.6726e-02, 3.5756e-03, 1.8611e-02, 2.6378e-02, 3.3823e-02,
 279      3.4811e-02, 3.5282e-02, 3.5386e-02, 3.6047e-02, 3.5201e-02, 3.4384e-02,
 280      3.2864e-02, 3.1932e-02, 3.2051e-02, 2.2539e-02,-2.5154e-02,-1.2050e-01,
 281     -1.2050e-01
 282   };
 283
 284   // PRF correction TODO to be replaced by the gaussian
 285   // approximation with full error parametrization and // moved to the clusterizer
 286   const Float_t cy[AliTRDgeometry::kNlayer][3] = {
 287     { 4.014e-04, 8.605e-03, -6.880e+00},
 288     {-3.061e-04, 9.663e-03, -6.789e+00},
 289     { 1.124e-03, 1.105e-02, -6.825e+00},
 290     {-1.527e-03, 1.231e-02, -6.777e+00},
 291     { 2.150e-03, 1.387e-02, -6.783e+00},
 292     {-1.296e-03, 1.486e-02, -6.825e+00}
 293   };
 294
 295   Int_t ily = AliTRDgeometry::GetLayer(c->GetDetector());
 296   x = c->GetX() - cx[c->GetLocalTimeBin()];
 297   y = c->GetY() + cy[ily][0] + cy[ily][1] * TMath::Sin(cy[ily][2] * c->GetCenter());
 298   return;
 299 }
 300
 301 //____________________________________________________________________
 302 Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
 303 {
 304   return fClusters[ic] ? TMath::Abs(fClusters[ic]->GetQ()) /fdX / TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]) : 0.;
 305 }
 306
 307 //____________________________________________________________________
 308 Double_t* AliTRDseedV1::GetProbability()
 309 {
 310 // Fill probability array for tracklet from the DB.
 311 //
 312 // Parameters
 313 //
 314 // Output
 315 //   returns pointer to the probability array and 0x0 if missing DB access
 316 //
 317 // Detailed description
 318
 319
 320   // retrive calibration db
 321   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 322   if (!calibration) {
 323     AliError("No access to calibration data");
 324     return 0x0;
 325   }
 326
 327   if (!fReconstructor) {
 328     AliError("Reconstructor not set.");
 329     return 0x0;
 330   }
 331
 332   // Retrieve the CDB container class with the parametric detector response
 333   const AliTRDCalPID *pd = calibration->GetPIDObject(fReconstructor->GetPIDMethod());
 334   if (!pd) {
 335     AliError("No access to AliTRDCalPID object");
 336     return 0x0;
 337   }
 338   //AliInfo(Form("Method[%d] : %s", fReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
 339
 340   // calculate tracklet length TO DO
 341   Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
 342   /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
 343
 344   //calculate dE/dx
 345   CookdEdx(fReconstructor->GetNdEdxSlices());
 346
 347   // Sets the a priori probabilities
 348   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
 349     fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, GetPlane());
 350   }
 351
 352   return &fProb[0];
 353 }
 354
 355 //____________________________________________________________________
 356 Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
 357 {
 358   //
 359   // Returns a quality measurement of the current seed
 360   //
 361
 362   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 363   return
 364       .5 * TMath::Abs(18.0 - fN2)
 365     + 10.* TMath::Abs(fYfit[1] - fYref[1])
 366     + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
 367     + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
 368 }
 369
 370 //____________________________________________________________________
 371 void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
 372 {
 373 // Computes covariance in the y-z plane at radial point x (in tracking coordinates)
 374 // and returns the results in the preallocated array cov[3] as :
 375 //   cov[0] = Var(y)
 376 //   cov[1] = Cov(yz)
 377 //   cov[2] = Var(z)
 378 //
 379 // Details
 380 //
 381 // For the linear transformation
 382 // BEGIN_LATEX
 383 // Y = T_{x} X^{T}
 384 // END_LATEX
 385 //   The error propagation has the general form
 386 // BEGIN_LATEX
 387 // C_{Y} = T_{x} C_{X} T_{x}^{T}
 388 // END_LATEX
 389 //  We apply this formula 2 times. First to calculate the covariance of the tracklet
 390 // at point x we consider:
 391 // BEGIN_LATEX
 392 // T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
 393 // END_LATEX
 394 // and secondly to take into account the tilt angle
 395 // BEGIN_LATEX
 396 // T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y)    0}{0   Var(z)}}
 397 // END_LATEX
 398 //
 399 // using simple trigonometrics one can write for this last case
 400 // BEGIN_LATEX
 401 // C_{Y}=#frac{1}{1+tg^{2}#alpha} #(){#splitline{(#sigma_{y}^{2}+tg^{2}#alpha#sigma_{z}^{2}) __ tg#alpha(#sigma_{z}^{2}-#sigma_{y}^{2})}{tg#alpha(#sigma_{z}^{2}-#sigma_{y}^{2}) __ (#sigma_{z}^{2}+tg^{2}#alpha#sigma_{y}^{2})}}
 402 // END_LATEX
 403 // which can be aproximated for small alphas (2 deg) with
 404 // BEGIN_LATEX
 405 // C_{Y}=#(){#splitline{#sigma_{y}^{2} __ (#sigma_{z}^{2}-#sigma_{y}^{2})tg#alpha}{((#sigma_{z}^{2}-#sigma_{y}^{2})tg#alpha __ #sigma_{z}^{2}}}
 406 // END_LATEX
 407 //
 408 // before applying the tilt rotation we also apply systematic uncertainties to the tracklet
 409 // position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might
 410 // account for extra misalignment/miscalibration uncertainties.
 411 //
 412 // Author :
 413 // Alex Bercuci <A.Bercuci@gsi.de>
 414 // Date : Jan 8th 2009
 415 //
 416   Double_t xr     = fX0-x;
 417   Double_t sy2    = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
 418   Double_t sz2    = fPadLength*fPadLength/12.;
 419
 420   // insert systematic uncertainties
 421   Double_t sys[15];
 422   fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
 423   sy2 += sys[0];
 424   sz2 += sys[1];
 425
 426   // rotate covariance matrix
 427   Double_t t2 = fTilt*fTilt;
 428   Double_t correction = 1./(1. + t2);
 429   cov[0] = (sy2+t2*sz2)*correction;
 430   cov[1] = fTilt*(sz2 - sy2)*correction;
 431   cov[2] = (t2*sy2+sz2)*correction;
 432 }
 433
 434
 435 //____________________________________________________________________
 436 void AliTRDseedV1::SetExB()
 437 {
 438 // Retrive the tg(a_L) from OCDB. The following information are used
 439 //  - detector index
 440 //  - column and row position of first attached cluster.
 441 //
 442 // If no clusters are attached to the tracklet a random central chamber
 443 // position (c=70, r=7) will be used to retrieve the drift velocity.
 444 //
 445 // Author :
 446 // Alex Bercuci <A.Bercuci@gsi.de>
 447 // Date : Jan 8th 2009
 448 //
 449
 450   AliCDBManager *cdb = AliCDBManager::Instance();
 451   if(cdb->GetRun() < 0){
 452     AliError("OCDB manager not properly initialized");
 453     return;
 454   }
 455
 456   AliTRDcalibDB *fCalib = AliTRDcalibDB::Instance();
 457   AliTRDCalROC  *fCalVdriftROC = fCalib->GetVdriftROC(fDet);
 458   const AliTRDCalDet  *fCalVdriftDet = fCalib->GetVdriftDet();
 459
 460   Int_t col = 70, row = 7;
 461   AliTRDcluster **c = &fClusters[0];
 462   if(fN){
 463     Int_t ic = 0;
 464     while (ic<AliTRDseed::knTimebins && !(*c)){ic++; c++;}
 465     if(*c){
 466       col = (*c)->GetPadCol();
 467       row = (*c)->GetPadRow();
 468     }
 469   }
 470
 471   Double_t vd = fCalVdriftDet->GetValue(fDet) * fCalVdriftROC->GetValue(col, row);
 472   fExB   = fCalib->GetOmegaTau(vd, -0.1*AliTracker::GetBz());
 473 }
 474
 475 //____________________________________________________________________
 476 void AliTRDseedV1::SetOwner()
 477 {
 478   //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
 479
 480   if(TestBit(kOwner)) return;
 481   for(int ic=0; ic<knTimebins; ic++){
 482     if(!fClusters[ic]) continue;
 483     fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
 484   }
 485   SetBit(kOwner);
 486 }
 487
 488 //____________________________________________________________________
 489 Bool_t  AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr, AliTRDcluster *c)
 490 {
 491   //
 492   // Iterative process to register clusters to the seed.
 493   // In iteration 0 we try only one pad-row and if quality not
 494   // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
 495   //
 496   // debug level 7
 497   //
 498
 499   if(!fReconstructor->GetRecoParam() ){
 500     AliError("Seed can not be used without a valid RecoParam.");
 501     return kFALSE;
 502   }
 503
 504   AliTRDchamberTimeBin *layer = 0x0;
 505   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7){
 506     AliTRDtrackingChamber ch(*chamber);
 507     ch.SetOwner();
 508     TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
 509     cstreamer << "AttachClustersIter"
 510       << "chamber.="   << &ch
 511       << "tracklet.="  << this
 512       << "\n";
 513   }
 514
 515   Float_t  tquality;
 516   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
 517   Double_t kroadz = fPadLength * .5 + 1.;
 518
 519   // initialize configuration parameters
 520   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 521   Int_t   niter = kZcorr ? 1 : 2;
 522
 523   Double_t yexp, zexp;
 524   Int_t ncl = 0;
 525   // start seed update
 526   for (Int_t iter = 0; iter < niter; iter++) {
 527     ncl = 0;
 528     for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 529       if(!(layer = chamber->GetTB(iTime))) continue;
 530       if(!Int_t(*layer)) continue;
 531
 532       // define searching configuration
 533       Double_t dxlayer = layer->GetX() - fX0;
 534       if(c){
 535         zexp = c->GetZ();
 536         //Try 2 pad-rows in second iteration
 537         if (iter > 0) {
 538           zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
 539           if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
 540           if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
 541         }
 542       } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
 543       yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
 544
 545       // Get and register cluster
 546       Int_t    index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
 547       if (index < 0) continue;
 548       AliTRDcluster *cl = (*layer)[index];
 549
 550       fIndexes[iTime]  = layer->GetGlobalIndex(index);
 551       fClusters[iTime] = cl;
 552       fY[iTime]        = cl->GetY();
 553       fZ[iTime]        = cl->GetZ();
 554       ncl++;
 555     }
 556     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
 557
 558     if(ncl>1){
 559       // calculate length of the time bin (calibration aware)
 560       Int_t irp = 0; Float_t x[2]; Int_t tb[2];
 561       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 562         if(!fClusters[iTime]) continue;
 563         x[irp]  = fClusters[iTime]->GetX();
 564         tb[irp] = iTime;
 565         irp++;
 566         if(irp==2) break;
 567       }
 568       fdX = (x[1] - x[0]) / (tb[0] - tb[1]);
 569
 570       // update X0 from the clusters (calibration/alignment aware)
 571       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 572         if(!(layer = chamber->GetTB(iTime))) continue;
 573         if(!layer->IsT0()) continue;
 574         if(fClusters[iTime]){
 575           fX0 = fClusters[iTime]->GetX();
 576           break;
 577         } else { // we have to infere the position of the anode wire from the other clusters
 578           for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
 579             if(!fClusters[jTime]) continue;
 580             fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
 581             break;
 582           }
 583         }
 584       }
 585
 586       // update YZ reference point
 587       // TODO
 588
 589       // update x reference positions (calibration/alignment aware)
 590       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 591         if(!fClusters[iTime]) continue;
 592         fX[iTime] = fX0 - fClusters[iTime]->GetX();
 593       }
 594
 595       AliTRDseed::Update();
 596     }
 597     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
 598
 599     if(IsOK()){
 600       tquality = GetQuality(kZcorr);
 601       if(tquality < quality) break;
 602       else quality = tquality;
 603     }
 604     kroadz *= 2.;
 605   } // Loop: iter
 606   if (!IsOK()) return kFALSE;
 607
 608   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
 609
 610   // refit tracklet with errors
 611   //SetExB(); Fit(kFALSE, 2);
 612
 613   UpdateUsed();
 614   return kTRUE;
 615 }
 616
 617 //____________________________________________________________________
 618 Bool_t  AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber
 619                                       ,Bool_t kZcorr)
 620 {
 621   //
 622   // Projective algorithm to attach clusters to seeding tracklets
 623   //
 624   // Parameters
 625   //
 626   // Output
 627   //
 628   // Detailed description
 629   // 1. Collapse x coordinate for the full detector plane
 630   // 2. truncated mean on y (r-phi) direction
 631   // 3. purge clusters
 632   // 4. truncated mean on z direction
 633   // 5. purge clusters
 634   // 6. fit tracklet
 635   //
 636
 637   if(!fReconstructor->GetRecoParam() ){
 638     AliError("Seed can not be used without a valid RecoParam.");
 639     return kFALSE;
 640   }
 641
 642   const Int_t kClusterCandidates = 2 * knTimebins;
 643
 644   //define roads
 645   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
 646   Double_t kroadz = fPadLength * 1.5 + 1.;
 647   // correction to y for the tilting angle
 648   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 649
 650   // working variables
 651   AliTRDcluster *clusters[kClusterCandidates];
 652   Double_t cond[4], yexp[knTimebins], zexp[knTimebins],
 653     yres[kClusterCandidates], zres[kClusterCandidates];
 654   Int_t ncl, *index = 0x0, tboundary[knTimebins];
 655
 656   // Do cluster projection
 657   AliTRDchamberTimeBin *layer = 0x0;
 658   Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
 659   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 660     if(!(layer = chamber->GetTB(iTime))) continue;
 661     if(!Int_t(*layer)) continue;
 662
 663     fX[iTime] = layer->GetX() - fX0;
 664     zexp[iTime] = fZref[0] + fZref[1] * fX[iTime];
 665     yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr;
 666
 667     // build condition and process clusters
 668     cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady;
 669     cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz;
 670     layer->GetClusters(cond, index, ncl);
 671     for(Int_t ic = 0; ic<ncl; ic++){
 672       AliTRDcluster *c = layer->GetCluster(index[ic]);
 673       clusters[nYclusters] = c;
 674       yres[nYclusters++] = c->GetY() - yexp[iTime];
 675       if(nYclusters >= kClusterCandidates) {
 676         AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates));
 677         kEXIT = kTRUE;
 678         break;
 679       }
 680     }
 681     tboundary[iTime] = nYclusters;
 682     if(kEXIT) break;
 683   }
 684
 685   // Evaluate truncated mean on the y direction
 686   Double_t mean, sigma;
 687   AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2);
 688   // purge cluster candidates
 689   Int_t nZclusters = 0;
 690   for(Int_t ic = 0; ic<nYclusters; ic++){
 691     if(yres[ic] - mean > 4. * sigma){
 692       clusters[ic] = 0x0;
 693       continue;
 694     }
 695     zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()];
 696   }
 697
 698   // Evaluate truncated mean on the z direction
 699   AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2);
 700   // purge cluster candidates
 701   for(Int_t ic = 0; ic<nZclusters; ic++){
 702     if(zres[ic] - mean > 4. * sigma){
 703       clusters[ic] = 0x0;
 704       continue;
 705     }
 706   }
 707
 708
 709   // Select only one cluster/TimeBin
 710   Int_t lastCluster = 0;
 711   fN2 = 0;
 712   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 713     ncl = tboundary[iTime] - lastCluster;
 714     if(!ncl) continue;
 715     Int_t iptr = lastCluster;
 716     if(ncl > 1){
 717       Float_t dold = 9999.;
 718       for(int ic=lastCluster; ic<tboundary[iTime]; ic++){
 719         if(!clusters[ic]) continue;
 720         Float_t y = yexp[iTime] - clusters[ic]->GetY();
 721         Float_t z = zexp[iTime] - clusters[ic]->GetZ();
 722         Float_t d = y * y + z * z;
 723         if(d > dold) continue;
 724         dold = d;
 725         iptr = ic;
 726       }
 727     }
 728     fIndexes[iTime]  = chamber->GetTB(iTime)->GetGlobalIndex(iptr);
 729     fClusters[iTime] = clusters[iptr];
 730     fY[iTime]        = clusters[iptr]->GetY();
 731     fZ[iTime]        = clusters[iptr]->GetZ();
 732     lastCluster      = tboundary[iTime];
 733     fN2++;
 734   }
 735
 736   // number of minimum numbers of clusters expected for the tracklet
 737   Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
 738   if (fN2 < kClmin){
 739     AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
 740     fN2 = 0;
 741     return kFALSE;
 742   }
 743
 744   // update used clusters
 745   fNUsed = 0;
 746   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 747     if(!fClusters[iTime]) continue;
 748     if((fClusters[iTime]->IsUsed())) fNUsed++;
 749   }
 750
 751   if (fN2-fNUsed < kClmin){
 752     AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
 753     fN2 = 0;
 754     return kFALSE;
 755   }
 756
 757   return kTRUE;
 758 }
 759
 760 //____________________________________________________________
 761 void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
 762 {
 763 //   Fill in all derived information. It has to be called after recovery from file or HLT.
 764 //   The primitive data are
 765 //   - list of clusters
 766 //   - detector (as the detector will be removed from clusters)
 767 //   - position of anode wire (fX0) - temporary
 768 //   - track reference position and direction
 769 //   - momentum of the track
 770 //   - time bin length [cm]
 771 //
 772 //   A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
 773 //
 774   fReconstructor = rec;
 775   AliTRDgeometry g;
 776   AliTRDpadPlane *pp = g.GetPadPlane(fDet);
 777   fTilt      = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
 778   fPadLength = pp->GetLengthIPad();
 779   fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
 780   fTgl = fZref[1];
 781   fN = 0; fN2 = 0; fMPads = 0.;
 782   AliTRDcluster **cit = &fClusters[0];
 783   for(Int_t ic = knTimebins; ic--; cit++){
 784     if(!(*cit)) return;
 785     fN++; fN2++;
 786     fX[ic] = (*cit)->GetX() - fX0;
 787     fY[ic] = (*cit)->GetY();
 788     fZ[ic] = (*cit)->GetZ();
 789   }
 790   Update(); // Fit();
 791   CookLabels();
 792   GetProbability();
 793 }
 794
 795
 796 //____________________________________________________________________
 797 Bool_t AliTRDseedV1::Fit(Bool_t tilt, Int_t errors)
 798 {
 799   //
 800   // Linear fit of the tracklet
 801   //
 802   // Parameters :
 803   //
 804   // Output :
 805   //  True if successful
 806   //
 807   // Detailed description
 808   // 2. Check if tracklet crosses pad row boundary
 809   // 1. Calculate residuals in the y (r-phi) direction
 810   // 3. Do a Least Square Fit to the data
 811   //
 812
 813   const Int_t kClmin = 8;
 814
 815
 816   // cluster error parametrization parameters
 817   // 1. sy total charge
 818   const Float_t sq0inv = 0.019962; // [1/q0]
 819   const Float_t sqb    = 1.0281564;    //[cm]
 820   // 2. sy for the PRF
 821   const Float_t scy[AliTRDgeometry::kNlayer][4] = {
 822     {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
 823     {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
 824     {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
 825     {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
 826     {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
 827     {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
 828   };
 829   // 3. sy parallel to the track
 830   const Float_t sy0 =  2.649e-02; // [cm]
 831   const Float_t sya = -8.864e-04; // [cm]
 832   const Float_t syb = -2.435e-01; // [cm]
 833
 834   // 4. sx parallel to the track
 835   const Float_t sxgc = 5.427e-02;
 836   const Float_t sxgm = 7.783e-01;
 837   const Float_t sxgs = 2.743e-01;
 838   const Float_t sxe0 =-2.065e+00;
 839   const Float_t sxe1 =-2.978e-02;
 840
 841   // 5. sx perpendicular to the track
 842 //   const Float_t sxd0 = 1.881e-02;
 843 //   const Float_t sxd1 =-4.101e-01;
 844 //   const Float_t sxd2 = 1.572e+00;
 845
 846   // get track direction
 847   Double_t y0   = fYref[0];
 848   Double_t dydx = fYref[1];
 849   Double_t z0   = fZref[0];
 850   Double_t dzdx = fZref[1];
 851   Double_t yt, zt;
 852
 853   const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
 854   AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
 855   TLinearFitter  fitterY(1, "pol1");
 856   // convertion factor from square to gauss distribution for sigma
 857   Double_t convert = 1./TMath::Sqrt(12.);
 858
 859   // book cluster information
 860   Double_t q, xc[knTimebins], yc[knTimebins], zc[knTimebins], sy[knTimebins], sz[knTimebins];
 861   Int_t zRow[knTimebins];
 862
 863   Int_t ily = AliTRDgeometry::GetLayer(fDet);
 864   fN = 0; fXref = 0.; Double_t ssx = 0.;
 865   AliTRDcluster *c=0x0, **jc = &fClusters[0];
 866   for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
 867     zRow[ic] = -1;
 868     xc[ic]  = -1.;
 869     yc[ic]  = 999.;
 870     zc[ic]  = 999.;
 871     sy[ic]  = 0.;
 872     sz[ic]  = 0.;
 873     if(!(c = (*jc))) continue;
 874     if(!c->IsInChamber()) continue;
 875
 876     Float_t w = 1.;
 877     if(c->GetNPads()>4) w = .5;
 878     if(c->GetNPads()>5) w = .2;
 879
 880     zRow[fN] = c->GetPadRow();
 881     // correct cluster position for PRF and v drift
 882     Double_t x, y; GetClusterXY(c, x, y);
 883     xc[fN]   = fX0 - x;
 884     yc[fN]   = y;
 885     zc[fN]   = c->GetZ();
 886
 887     // extrapolated y value for the track
 888     yt = y0 - xc[fN]*dydx;
 889     // extrapolated z value for the track
 890     zt = z0 - xc[fN]*dzdx;
 891     // tilt correction
 892     if(tilt) yc[fN] -= fTilt*(zc[fN] - zt);
 893
 894     // ELABORATE CLUSTER ERROR
 895     // TODO to be moved to AliTRDcluster
 896     q = TMath::Abs(c->GetQ());
 897     Double_t tgg = (dydx-fExB)/(1.+dydx*fExB); tgg *= tgg;
 898     // basic y error (|| to track).
 899     sy[fN]  = xc[fN] < AliTRDgeometry::CamHght() ? 2. : sy0 + sya*TMath::Exp(1./(xc[fN]+syb));
 900     //printf("cluster[%d]\n\tsy[0] = %5.3e [um]\n", fN,  sy[fN]*1.e4);
 901     // y error due to total charge
 902     sy[fN] += sqb*(1./q - sq0inv);
 903     //printf("\tsy[1] = %5.3e [um]\n", sy[fN]*1.e4);
 904     // y error due to PRF
 905     sy[fN] += scy[ily][0]*TMath::Gaus(c->GetCenter(), scy[ily][1], scy[ily][2]) - scy[ily][3];
 906     //printf("\tsy[2] = %5.3e [um]\n", sy[fN]*1.e4);
 907
 908     sy[fN] *= sy[fN];
 909
 910     // ADD ERROR ON x
 911     // error of drift length parallel to the track
 912     Double_t sx = sxgc*TMath::Gaus(xc[fN], sxgm, sxgs) + TMath::Exp(sxe0+sxe1*xc[fN]); // [cm]
 913     //printf("\tsx[0] = %5.3e [um]\n", sx*1.e4);
 914     // error of drift length perpendicular to the track
 915     //sx += sxd0 + sxd1*d + sxd2*d*d;
 916     sx *= sx; // square sx
 917     // update xref
 918     fXref += xc[fN]/sx; ssx+=1./sx;
 919
 920     // add error from ExB
 921     if(errors>0) sy[fN] += fExB*fExB*sx;
 922     //printf("\tsy[3] = %5.3e [um^2]\n", sy[fN]*1.e8);
 923
 924     // global radial error due to misalignment/miscalibration
 925     Double_t sx0  = 0.; sx0 *= sx0;
 926     // add sx contribution to sy due to track angle
 927     if(errors>1) sy[fN] += tgg*(sx+sx0);
 928     // TODO we should add tilt pad correction here
 929     //printf("\tsy[4] = %5.3e [um^2]\n", sy[fN]*1.e8);
 930     c->SetSigmaY2(sy[fN]);
 931
 932     sy[fN]  = TMath::Sqrt(sy[fN]);
 933     fitterY.AddPoint(&xc[fN], yc[fN]/*-yt*/, sy[fN]);
 934
 935     sz[fN]   = fPadLength*convert;
 936     fitterZ.AddPoint(&xc[fN], zc[fN], sz[fN]);
 937     fN++;
 938   }
 939   // to few clusters
 940   if (fN < kClmin) return kFALSE;
 941
 942   // fit XY
 943   fitterY.Eval();
 944   fYfit[0] = fitterY.GetParameter(0);
 945   fYfit[1] = -fitterY.GetParameter(1);
 946   // store covariance
 947   Double_t *p = fitterY.GetCovarianceMatrix();
 948   fCov[0] = p[0]; // variance of y0
 949   fCov[1] = p[1]; // covariance of y0, dydx
 950   fCov[2] = p[3]; // variance of dydx
 951   // store ref radial position.
 952   fXref /= ssx; fXref = fX0 - fXref;
 953
 954   // check par row crossing
 955   Int_t zN[2*AliTRDseed::knTimebins];
 956   Int_t nz = AliTRDtrackerV1::Freq(fN, zRow, zN, kFALSE);
 957   // more than one pad row crossing
 958   if(nz>2) return kFALSE;
 959
 960
 961   // determine z offset of the fit
 962   Float_t zslope = 0.;
 963   Int_t nchanges = 0, nCross = 0;
 964   if(nz==2){ // tracklet is crossing pad row
 965     // Find the break time allowing one chage on pad-rows
 966     // with maximal number of accepted clusters
 967     Int_t padRef = zRow[0];
 968     for (Int_t ic=1; ic<fN; ic++) {
 969       if(zRow[ic] == padRef) continue;
 970
 971       // debug
 972       if(zRow[ic-1] == zRow[ic]){
 973         printf("ERROR in pad row change!!!\n");
 974       }
 975
 976       // evaluate parameters of the crossing point
 977       Float_t sx = (xc[ic-1] - xc[ic])*convert;
 978       fCross[0] = .5 * (xc[ic-1] + xc[ic]);
 979       fCross[2] = .5 * (zc[ic-1] + zc[ic]);
 980       fCross[3] = TMath::Max(dzdx * sx, .01);
 981       zslope    = zc[ic-1] > zc[ic] ? 1. : -1.;
 982       padRef    = zRow[ic];
 983       nCross    = ic;
 984       nchanges++;
 985     }
 986   }
 987
 988   // condition on nCross and reset nchanges TODO
 989
 990   if(nchanges==1){
 991     if(dzdx * zslope < 0.){
 992       AliInfo("tracklet direction does not correspond to the track direction. TODO.");
 993     }
 994     SetBit(kRowCross, kTRUE); // mark pad row crossing
 995     fitterZ.AddPoint(&fCross[0], fCross[2], fCross[3]);
 996     fitterZ.Eval();
 997     //zc[nc] = fitterZ.GetFunctionParameter(0);
 998     fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
 999     fCross[0] = fX0 - fCross[0];
1000   } else if(nchanges > 1){ // debug
1001     AliError("N pad row crossing > 1.");
1002     return kFALSE;
1003   }
1004
1005   UpdateUsed();
1006
1007   return kTRUE;
1008 }
1009
1010
1011 //___________________________________________________________________
1012 void AliTRDseedV1::Print(Option_t *o) const
1013 {
1014   //
1015   // Printing the seedstatus
1016   //
1017
1018   AliInfo(Form("Det[%3d] Tilt[%+6.2f] Pad[%5.2f]", fDet, fTilt, fPadLength));
1019   AliInfo(Form("Nattach[%2d] Nfit[%2d] Nuse[%2d] pads[%f]", fN, fN2, fNUsed, fMPads));
1020   AliInfo(Form("x[%7.2f] y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fX0, fYfit[0], fZfit[0], fYfit[1], fZfit[1]));
1021   AliInfo(Form("Ref        y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fYref[0], fZref[0], fYref[1], fZref[1]))
1022
1023
1024   if(strcmp(o, "a")!=0) return;
1025
1026   AliTRDcluster* const* jc = &fClusters[0];
1027   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++, jc++) {
1028     if(!(*jc)) continue;
1029     (*jc)->Print(o);
1030   }
1031 }
1032
1033
1034 //___________________________________________________________________
1035 Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1036 {
1037   // Checks if current instance of the class has the same essential members
1038   // as the given one
1039
1040   if(!o) return kFALSE;
1041   const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1042   if(!inTracklet) return kFALSE;
1043
1044   for (Int_t i = 0; i < 2; i++){
1045     if ( fYref[i] != inTracklet->GetYref(i) ) return kFALSE;
1046     if ( fZref[i] != inTracklet->GetZref(i) ) return kFALSE;
1047   }
1048
1049   if ( fSigmaY != inTracklet->GetSigmaY() ) return kFALSE;
1050   if ( fSigmaY2 != inTracklet->GetSigmaY2() ) return kFALSE;
1051   if ( fTilt != inTracklet->GetTilt() ) return kFALSE;
1052   if ( fPadLength != inTracklet->GetPadLength() ) return kFALSE;
1053
1054   for (Int_t i = 0; i < knTimebins; i++){
1055     if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1056     if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1057     if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1058     if ( fIndexes[i] != inTracklet->GetIndexes(i) ) return kFALSE;
1059     if ( fUsable[i] != inTracklet->IsUsable(i) ) return kFALSE;
1060   }
1061
1062   for (Int_t i=0; i < 2; i++){
1063     if ( fYfit[i] != inTracklet->GetYfit(i) ) return kFALSE;
1064     if ( fZfit[i] != inTracklet->GetZfit(i) ) return kFALSE;
1065     if ( fYfitR[i] != inTracklet->GetYfitR(i) ) return kFALSE;
1066     if ( fZfitR[i] != inTracklet->GetZfitR(i) ) return kFALSE;
1067     if ( fLabels[i] != inTracklet->GetLabels(i) ) return kFALSE;
1068   }
1069
1070   if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1071   if ( fZProb != inTracklet->GetZProb() ) return kFALSE;
1072   if ( fN2 != inTracklet->GetN2() ) return kFALSE;
1073   if ( fNUsed != inTracklet->GetNUsed() ) return kFALSE;
1074   if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1075   if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
1076   if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
1077
1078   if ( fC != inTracklet->GetC() ) return kFALSE;
1079   if ( fCC != inTracklet->GetCC() ) return kFALSE;
1080   if ( fChi2 != inTracklet->GetChi2() ) return kFALSE;
1081   //  if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1082
1083   if ( fDet != inTracklet->GetDetector() ) return kFALSE;
1084   if ( fMom != inTracklet->GetMomentum() ) return kFALSE;
1085   if ( fdX != inTracklet->GetdX() ) return kFALSE;
1086
1087   for (Int_t iCluster = 0; iCluster < knTimebins; iCluster++){
1088     AliTRDcluster *curCluster = fClusters[iCluster];
1089     AliTRDcluster *inCluster = inTracklet->GetClusters(iCluster);
1090     if (curCluster && inCluster){
1091       if (! curCluster->IsEqual(inCluster) ) {
1092         curCluster->Print();
1093         inCluster->Print();
1094         return kFALSE;
1095       }
1096     } else {
1097       // if one cluster exists, and corresponding
1098       // in other tracklet doesn't - return kFALSE
1099       if(curCluster || inCluster) return kFALSE;
1100     }
1101   }
1102   return kTRUE;
1103 }