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 #include "AliTRDCommonParam.h"
  49
  50 #include "Cal/AliTRDCalPID.h"
  51 #include "Cal/AliTRDCalROC.h"
  52 #include "Cal/AliTRDCalDet.h"
  53
  54 ClassImp(AliTRDseedV1)
  55
  56 //____________________________________________________________________
  57 AliTRDseedV1::AliTRDseedV1(Int_t det)
  58   :AliTRDtrackletBase()
  59   ,fReconstructor(0x0)
  60   ,fClusterIter(0x0)
  61   ,fExB(0.)
  62   ,fVD(0.)
  63   ,fT0(0.)
  64   ,fS2PRF(0.)
  65   ,fDiffL(0.)
  66   ,fDiffT(0.)
  67   ,fClusterIdx(0)
  68   ,fN(0)
  69   ,fDet(det)
  70   ,fTilt(0.)
  71   ,fPadLength(0.)
  72   ,fMom(0.)
  73   ,fdX(0.)
  74   ,fX0(0.)
  75   ,fX(0.)
  76   ,fY(0.)
  77   ,fZ(0.)
  78   ,fS2Y(0.)
  79   ,fS2Z(0.)
  80   ,fC(0.)
  81   ,fChi2(0.)
  82 {
  83   //
  84   // Constructor
  85   //
  86   for(Int_t ic=kNclusters; ic--;) fIndexes[ic] = -1;
  87   memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
  88   fYref[0] = 0.; fYref[1] = 0.;
  89   fZref[0] = 0.; fZref[1] = 0.;
  90   fYfit[0] = 0.; fYfit[1] = 0.;
  91   fZfit[0] = 0.; fZfit[1] = 0.;
  92   memset(fdEdx, 0, kNslices*sizeof(Float_t));
  93   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
  94   fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
  95   fLabels[2]=0;  // number of different labels for tracklet
  96   fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
  97   // covariance matrix [diagonal]
  98   // default sy = 200um and sz = 2.3 cm
  99   fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
 100   SetStandAlone(kFALSE);
 101 }
 102
 103 //____________________________________________________________________
 104 AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
 105   :AliTRDtrackletBase((AliTRDtrackletBase&)ref)
 106   ,fReconstructor(0x0)
 107   ,fClusterIter(0x0)
 108   ,fExB(0.)
 109   ,fVD(0.)
 110   ,fT0(0.)
 111   ,fS2PRF(0.)
 112   ,fDiffL(0.)
 113   ,fDiffT(0.)
 114   ,fClusterIdx(0)
 115   ,fN(0)
 116   ,fDet(-1)
 117   ,fTilt(0.)
 118   ,fPadLength(0.)
 119   ,fMom(0.)
 120   ,fdX(0.)
 121   ,fX0(0.)
 122   ,fX(0.)
 123   ,fY(0.)
 124   ,fZ(0.)
 125   ,fS2Y(0.)
 126   ,fS2Z(0.)
 127   ,fC(0.)
 128   ,fChi2(0.)
 129 {
 130   //
 131   // Copy Constructor performing a deep copy
 132   //
 133   if(this != &ref){
 134     ref.Copy(*this);
 135   }
 136   SetBit(kOwner, kFALSE);
 137   SetStandAlone(ref.IsStandAlone());
 138 }
 139
 140
 141 //____________________________________________________________________
 142 AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
 143 {
 144   //
 145   // Assignment Operator using the copy function
 146   //
 147
 148   if(this != &ref){
 149     ref.Copy(*this);
 150   }
 151   SetBit(kOwner, kFALSE);
 152
 153   return *this;
 154 }
 155
 156 //____________________________________________________________________
 157 AliTRDseedV1::~AliTRDseedV1()
 158 {
 159   //
 160   // Destructor. The RecoParam object belongs to the underlying tracker.
 161   //
 162
 163   //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
 164
 165   if(IsOwner()) {
 166     for(int itb=0; itb<kNclusters; itb++){
 167       if(!fClusters[itb]) continue;
 168       //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
 169       delete fClusters[itb];
 170       fClusters[itb] = 0x0;
 171     }
 172   }
 173 }
 174
 175 //____________________________________________________________________
 176 void AliTRDseedV1::Copy(TObject &ref) const
 177 {
 178   //
 179   // Copy function
 180   //
 181
 182   //AliInfo("");
 183   AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
 184
 185   target.fReconstructor = fReconstructor;
 186   target.fClusterIter   = 0x0;
 187   target.fExB           = fExB;
 188   target.fVD            = fVD;
 189   target.fT0            = fT0;
 190   target.fS2PRF         = fS2PRF;
 191   target.fDiffL         = fDiffL;
 192   target.fDiffT         = fDiffT;
 193   target.fClusterIdx    = 0;
 194   target.fN             = fN;
 195   target.fDet           = fDet;
 196   target.fTilt          = fTilt;
 197   target.fPadLength     = fPadLength;
 198   target.fMom           = fMom;
 199   target.fdX            = fdX;
 200   target.fX0            = fX0;
 201   target.fX             = fX;
 202   target.fY             = fY;
 203   target.fZ             = fZ;
 204   target.fS2Y           = fS2Y;
 205   target.fS2Z           = fS2Z;
 206   target.fC             = fC;
 207   target.fChi2          = fChi2;
 208
 209   memcpy(target.fIndexes, fIndexes, kNclusters*sizeof(Int_t));
 210   memcpy(target.fClusters, fClusters, kNclusters*sizeof(AliTRDcluster*));
 211   target.fYref[0] = fYref[0]; target.fYref[1] = fYref[1];
 212   target.fZref[0] = fZref[0]; target.fZref[1] = fZref[1];
 213   target.fYfit[0] = fYfit[0]; target.fYfit[1] = fYfit[1];
 214   target.fZfit[0] = fZfit[0]; target.fZfit[1] = fZfit[1];
 215   memcpy(target.fdEdx, fdEdx, kNslices*sizeof(Float_t));
 216   memcpy(target.fProb, fProb, AliPID::kSPECIES*sizeof(Float_t));
 217   memcpy(target.fLabels, fLabels, 3*sizeof(Int_t));
 218   memcpy(target.fRefCov, fRefCov, 3*sizeof(Double_t));
 219   memcpy(target.fCov, fCov, 3*sizeof(Double_t));
 220
 221   TObject::Copy(ref);
 222 }
 223
 224
 225 //____________________________________________________________
 226 Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
 227 {
 228 // Initialize this tracklet using the track information
 229 //
 230 // Parameters:
 231 //   track - the TRD track used to initialize the tracklet
 232 //
 233 // Detailed description
 234 // The function sets the starting point and direction of the
 235 // tracklet according to the information from the TRD track.
 236 //
 237 // Caution
 238 // The TRD track has to be propagated to the beginning of the
 239 // chamber where the tracklet will be constructed
 240 //
 241
 242   Double_t y, z;
 243   if(!track->GetProlongation(fX0, y, z)) return kFALSE;
 244   UpDate(track);
 245   return kTRUE;
 246 }
 247
 248
 249 //_____________________________________________________________________________
 250 void AliTRDseedV1::Reset()
 251 {
 252   //
 253   // Reset seed
 254   //
 255   fExB=0.;fVD=0.;fT0=0.;fS2PRF=0.;
 256   fDiffL=0.;fDiffT=0.;
 257   fClusterIdx=0;
 258   fN=0;
 259   fDet=-1;fTilt=0.;fPadLength=0.;
 260   fMom=0.;
 261   fdX=0.;fX0=0.; fX=0.; fY=0.; fZ=0.;
 262   fS2Y=0.; fS2Z=0.;
 263   fC=0.; fChi2 = 0.;
 264
 265   for(Int_t ic=kNclusters; ic--;) fIndexes[ic] = -1;
 266   memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
 267   fYref[0] = 0.; fYref[1] = 0.;
 268   fZref[0] = 0.; fZref[1] = 0.;
 269   fYfit[0] = 0.; fYfit[1] = 0.;
 270   fZfit[0] = 0.; fZfit[1] = 0.;
 271   memset(fdEdx, 0, kNslices*sizeof(Float_t));
 272   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
 273   fLabels[0]=-1; fLabels[1]=-1; // most freq MC labels
 274   fLabels[2]=0;  // number of different labels for tracklet
 275   fRefCov[0] = 1.; fRefCov[1] = 0.; fRefCov[2] = 1.;
 276   // covariance matrix [diagonal]
 277   // default sy = 200um and sz = 2.3 cm
 278   fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
 279 }
 280
 281 //____________________________________________________________________
 282 void AliTRDseedV1::UpDate(const AliTRDtrackV1 *trk)
 283 {
 284   // update tracklet reference position from the TRD track
 285   // Funny name to avoid the clash with the function AliTRDseed::Update() (has to be made obselete)
 286
 287   Double_t fSnp = trk->GetSnp();
 288   Double_t fTgl = trk->GetTgl();
 289   fMom = trk->GetP();
 290   fYref[1] = fSnp/(1. - fSnp*fSnp);
 291   fZref[1] = fTgl;
 292   SetCovRef(trk->GetCovariance());
 293
 294   Double_t dx = trk->GetX() - fX0;
 295   fYref[0] = trk->GetY() - dx*fYref[1];
 296   fZref[0] = trk->GetZ() - dx*fZref[1];
 297 }
 298
 299 //_____________________________________________________________________________
 300 void AliTRDseedV1::UpdateUsed()
 301 {
 302   //
 303   // Calculate number of used clusers in the tracklet
 304   //
 305
 306   Int_t nused = 0, nshared = 0;
 307   for (Int_t i = kNclusters; i--; ) {
 308     if (!fClusters[i]) continue;
 309     if(fClusters[i]->IsUsed()){
 310       nused++;
 311     } else if(fClusters[i]->IsShared()){
 312       if(IsStandAlone()) nused++;
 313       else nshared++;
 314     }
 315   }
 316   SetNUsed(nused);
 317   SetNShared(nshared);
 318 }
 319
 320 //_____________________________________________________________________________
 321 void AliTRDseedV1::UseClusters()
 322 {
 323   //
 324   // Use clusters
 325   //
 326   // In stand alone mode:
 327   // Clusters which are marked as used or shared from another track are
 328   // removed from the tracklet
 329   //
 330   // In barrel mode:
 331   // - Clusters which are used by another track become shared
 332   // - Clusters which are attached to a kink track become shared
 333   //
 334   AliTRDcluster **c = &fClusters[0];
 335   for (Int_t ic=kNclusters; ic--; c++) {
 336     if(!(*c)) continue;
 337     if(IsStandAlone()){
 338       if((*c)->IsShared() || (*c)->IsUsed()){
 339         if((*c)->IsShared()) SetNShared(GetNShared()-1);
 340         else SetNUsed(GetNUsed()-1);
 341         (*c) = 0x0;
 342         fIndexes[ic] = -1;
 343         SetN(GetN()-1);
 344         continue;
 345       }
 346     } else {
 347       if((*c)->IsUsed() || IsKink()){
 348         (*c)->SetShared();
 349         continue;
 350       }
 351     }
 352     (*c)->Use();
 353   }
 354 }
 355
 356
 357
 358 //____________________________________________________________________
 359 void AliTRDseedV1::CookdEdx(Int_t nslices)
 360 {
 361 // Calculates average dE/dx for all slices and store them in the internal array fdEdx.
 362 //
 363 // Parameters:
 364 //  nslices : number of slices for which dE/dx should be calculated
 365 // Output:
 366 //  store results in the internal array fdEdx. This can be accessed with the method
 367 //  AliTRDseedV1::GetdEdx()
 368 //
 369 // Detailed description
 370 // Calculates average dE/dx for all slices. Depending on the PID methode
 371 // the number of slices can be 3 (LQ) or 8(NN).
 372 // The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t))
 373 //
 374 // The following effects are included in the calculation:
 375 // 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
 376 // 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
 377 // 3. cluster size
 378 //
 379
 380   Int_t nclusters[kNslices];
 381   memset(nclusters, 0, kNslices*sizeof(Int_t));
 382   memset(fdEdx, 0, kNslices*sizeof(Float_t));
 383
 384   const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
 385
 386   AliTRDcluster *c = 0x0;
 387   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
 388     if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
 389     Float_t dx = TMath::Abs(fX0 - c->GetX());
 390
 391     // Filter clusters for dE/dx calculation
 392
 393     // 1.consider calibration effects for slice determination
 394     Int_t slice;
 395     if(dx<kDriftLength){ // TODO should be replaced by c->IsInChamber()
 396       slice = Int_t(dx * nslices / kDriftLength);
 397     } else slice = c->GetX() < fX0 ? nslices-1 : 0;
 398
 399
 400     // 2. take sharing into account
 401     Float_t w = /*c->IsShared() ? .5 :*/ 1.;
 402
 403     // 3. take into account large clusters TODO
 404     //w *= c->GetNPads() > 3 ? .8 : 1.;
 405
 406     //CHECK !!!
 407     fdEdx[slice]   += w * GetdQdl(ic); //fdQdl[ic];
 408     nclusters[slice]++;
 409   } // End of loop over clusters
 410
 411   //if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
 412   if(nslices == AliTRDpidUtil::kLQslices){
 413   // calculate mean charge per slice (only LQ PID)
 414     for(int is=0; is<nslices; is++){
 415       if(nclusters[is]) fdEdx[is] /= nclusters[is];
 416     }
 417   }
 418 }
 419
 420 //_____________________________________________________________________________
 421 void AliTRDseedV1::CookLabels()
 422 {
 423   //
 424   // Cook 2 labels for seed
 425   //
 426
 427   Int_t labels[200];
 428   Int_t out[200];
 429   Int_t nlab = 0;
 430   for (Int_t i = 0; i < kNclusters; i++) {
 431     if (!fClusters[i]) continue;
 432     for (Int_t ilab = 0; ilab < 3; ilab++) {
 433       if (fClusters[i]->GetLabel(ilab) >= 0) {
 434         labels[nlab] = fClusters[i]->GetLabel(ilab);
 435         nlab++;
 436       }
 437     }
 438   }
 439
 440   fLabels[2] = AliMathBase::Freq(nlab,labels,out,kTRUE);
 441   fLabels[0] = out[0];
 442   if ((fLabels[2]  > 1) && (out[3] > 1)) fLabels[1] = out[2];
 443 }
 444
 445
 446 //____________________________________________________________________
 447 void AliTRDseedV1::GetClusterXY(const AliTRDcluster *c, Double_t &x, Double_t &y)
 448 {
 449 // Return corrected position of the cluster taking into
 450 // account variation of the drift velocity with drift length.
 451
 452
 453   // drift velocity correction TODO to be moved to the clusterizer
 454   const Float_t cx[] = {
 455     -9.6280e-02, 1.3091e-01,-1.7415e-02,-9.9221e-02,-1.2040e-01,-9.5493e-02,
 456     -5.0041e-02,-1.6726e-02, 3.5756e-03, 1.8611e-02, 2.6378e-02, 3.3823e-02,
 457      3.4811e-02, 3.5282e-02, 3.5386e-02, 3.6047e-02, 3.5201e-02, 3.4384e-02,
 458      3.2864e-02, 3.1932e-02, 3.2051e-02, 2.2539e-02,-2.5154e-02,-1.2050e-01,
 459     -1.2050e-01
 460   };
 461
 462   // PRF correction TODO to be replaced by the gaussian
 463   // approximation with full error parametrization and // moved to the clusterizer
 464   const Float_t cy[AliTRDgeometry::kNlayer][3] = {
 465     { 4.014e-04, 8.605e-03, -6.880e+00},
 466     {-3.061e-04, 9.663e-03, -6.789e+00},
 467     { 1.124e-03, 1.105e-02, -6.825e+00},
 468     {-1.527e-03, 1.231e-02, -6.777e+00},
 469     { 2.150e-03, 1.387e-02, -6.783e+00},
 470     {-1.296e-03, 1.486e-02, -6.825e+00}
 471   };
 472
 473   Int_t ily = AliTRDgeometry::GetLayer(c->GetDetector());
 474   x = c->GetX() - cx[c->GetLocalTimeBin()];
 475   y = c->GetY() + cy[ily][0] + cy[ily][1] * TMath::Sin(cy[ily][2] * c->GetCenter());
 476   return;
 477 }
 478
 479 //____________________________________________________________________
 480 Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
 481 {
 482 // Using the linear approximation of the track inside one TRD chamber (TRD tracklet)
 483 // the charge per unit length can be written as:
 484 // BEGIN_LATEX
 485 // #frac{dq}{dl} = #frac{q_{c}}{dx * #sqrt{1 + #(){#frac{dy}{dx}}^{2}_{fit} + #(){#frac{dy}{dx}}^{2}_{ref}}}
 486 // END_LATEX
 487 // where qc is the total charge collected in the current time bin and dx is the length
 488 // of the time bin. For the moment (Jan 20 2009) only pad row cross corrections are
 489 // considered for the charge but none are applied for drift velocity variations along
 490 // the drift region or assymetry of the TRF
 491 //
 492 // Author : Alex Bercuci <A.Bercuci@gsi.de>
 493 //
 494   Float_t dq = 0.;
 495   if(fClusters[ic]) dq += TMath::Abs(fClusters[ic]->GetQ());
 496   if(fClusters[ic+kNtb]) dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
 497   if(dq<1.e-3 || fdX < 1.e-3) return 0.;
 498
 499   return dq/fdX/TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
 500 }
 501
 502 //____________________________________________________________________
 503 Float_t* AliTRDseedV1::GetProbability(Bool_t force)
 504 {
 505   if(!force) return &fProb[0];
 506   if(!CookPID()) return 0x0;
 507   return &fProb[0];
 508 }
 509
 510 //____________________________________________________________
 511 Bool_t AliTRDseedV1::CookPID()
 512 {
 513 // Fill probability array for tracklet from the DB.
 514 //
 515 // Parameters
 516 //
 517 // Output
 518 //   returns pointer to the probability array and 0x0 if missing DB access
 519 //
 520 // Detailed description
 521
 522
 523   // retrive calibration db
 524   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 525   if (!calibration) {
 526     AliError("No access to calibration data");
 527     return kFALSE;
 528   }
 529
 530   if (!fReconstructor) {
 531     AliError("Reconstructor not set.");
 532     return kFALSE;
 533   }
 534
 535   // Retrieve the CDB container class with the parametric detector response
 536   const AliTRDCalPID *pd = calibration->GetPIDObject(fReconstructor->GetPIDMethod());
 537   if (!pd) {
 538     AliError("No access to AliTRDCalPID object");
 539     return kFALSE;
 540   }
 541   //AliInfo(Form("Method[%d] : %s", fReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
 542
 543   // calculate tracklet length TO DO
 544   Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
 545   /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
 546
 547   //calculate dE/dx
 548   CookdEdx(fReconstructor->GetNdEdxSlices());
 549
 550   // Sets the a priori probabilities
 551   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
 552     fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, GetPlane());
 553   }
 554
 555   return kTRUE;
 556 }
 557
 558 //____________________________________________________________________
 559 Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
 560 {
 561   //
 562   // Returns a quality measurement of the current seed
 563   //
 564
 565   Float_t zcorr = kZcorr ? fTilt * (fZfit[0] - fZref[0]) : 0.;
 566   return
 567       .5 * TMath::Abs(18.0 - GetN())
 568     + 10.* TMath::Abs(fYfit[1] - fYref[1])
 569     + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
 570     + 2. * TMath::Abs(fZfit[0] - fZref[0]) / fPadLength;
 571 }
 572
 573 //____________________________________________________________________
 574 void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
 575 {
 576 // Computes covariance in the y-z plane at radial point x (in tracking coordinates)
 577 // and returns the results in the preallocated array cov[3] as :
 578 //   cov[0] = Var(y)
 579 //   cov[1] = Cov(yz)
 580 //   cov[2] = Var(z)
 581 //
 582 // Details
 583 //
 584 // For the linear transformation
 585 // BEGIN_LATEX
 586 // Y = T_{x} X^{T}
 587 // END_LATEX
 588 //   The error propagation has the general form
 589 // BEGIN_LATEX
 590 // C_{Y} = T_{x} C_{X} T_{x}^{T}
 591 // END_LATEX
 592 //  We apply this formula 2 times. First to calculate the covariance of the tracklet
 593 // at point x we consider:
 594 // BEGIN_LATEX
 595 // T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
 596 // END_LATEX
 597 // and secondly to take into account the tilt angle
 598 // BEGIN_LATEX
 599 // T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y)    0}{0   Var(z)}}
 600 // END_LATEX
 601 //
 602 // using simple trigonometrics one can write for this last case
 603 // BEGIN_LATEX
 604 // 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})}}
 605 // END_LATEX
 606 // which can be aproximated for small alphas (2 deg) with
 607 // BEGIN_LATEX
 608 // 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}}}
 609 // END_LATEX
 610 //
 611 // before applying the tilt rotation we also apply systematic uncertainties to the tracklet
 612 // position which can be tunned from outside via the AliTRDrecoParam::SetSysCovMatrix(). They might
 613 // account for extra misalignment/miscalibration uncertainties.
 614 //
 615 // Author :
 616 // Alex Bercuci <A.Bercuci@gsi.de>
 617 // Date : Jan 8th 2009
 618 //
 619
 620
 621   Double_t xr     = fX0-x;
 622   Double_t sy2    = fCov[0] +2.*xr*fCov[1] + xr*xr*fCov[2];
 623   Double_t sz2    = fPadLength*fPadLength/12.;
 624
 625   // insert systematic uncertainties
 626   Double_t sys[15];
 627   fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
 628   sy2 += sys[0];
 629   sz2 += sys[1];
 630
 631   // rotate covariance matrix
 632   Double_t t2 = fTilt*fTilt;
 633   Double_t correction = 1./(1. + t2);
 634   cov[0] = (sy2+t2*sz2)*correction;
 635   cov[1] = fTilt*(sz2 - sy2)*correction;
 636   cov[2] = (t2*sy2+sz2)*correction;
 637 }
 638
 639
 640 //____________________________________________________________________
 641 void AliTRDseedV1::Calibrate()
 642 {
 643 // Retrieve calibration and position parameters from OCDB.
 644 // The following information are used
 645 //  - detector index
 646 //  - column and row position of first attached cluster. If no clusters are attached
 647 // to the tracklet a random central chamber position (c=70, r=7) will be used.
 648 //
 649 // The following information is cached in the tracklet
 650 //   t0 (trigger delay)
 651 //   drift velocity
 652 //   PRF width
 653 //   omega*tau = tg(a_L)
 654 //   diffusion coefficients (longitudinal and transversal)
 655 //
 656 // Author :
 657 // Alex Bercuci <A.Bercuci@gsi.de>
 658 // Date : Jan 8th 2009
 659 //
 660
 661   AliCDBManager *cdb = AliCDBManager::Instance();
 662   if(cdb->GetRun() < 0){
 663     AliError("OCDB manager not properly initialized");
 664     return;
 665   }
 666
 667   AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
 668   AliTRDCalROC  *vdROC = calib->GetVdriftROC(fDet),
 669                 *t0ROC = calib->GetT0ROC(fDet);;
 670   const AliTRDCalDet *vdDet = calib->GetVdriftDet();
 671   const AliTRDCalDet *t0Det = calib->GetT0Det();
 672
 673   Int_t col = 70, row = 7;
 674   AliTRDcluster **c = &fClusters[0];
 675   if(GetN()){
 676     Int_t ic = 0;
 677     while (ic<kNclusters && !(*c)){ic++; c++;}
 678     if(*c){
 679       col = (*c)->GetPadCol();
 680       row = (*c)->GetPadRow();
 681     }
 682   }
 683
 684   fT0    = t0Det->GetValue(fDet) + t0ROC->GetValue(col,row);
 685   fVD    = vdDet->GetValue(fDet) * vdROC->GetValue(col, row);
 686   fS2PRF = calib->GetPRFWidth(fDet, col, row); fS2PRF *= fS2PRF;
 687   fExB   = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
 688   AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
 689   fDiffT, fVD);
 690   SetBit(kCalib, kTRUE);
 691 }
 692
 693 //____________________________________________________________________
 694 void AliTRDseedV1::SetOwner()
 695 {
 696   //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
 697
 698   if(TestBit(kOwner)) return;
 699   for(int ic=0; ic<kNclusters; ic++){
 700     if(!fClusters[ic]) continue;
 701     fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
 702   }
 703   SetBit(kOwner);
 704 }
 705
 706 // //____________________________________________________________________
 707 // Bool_t       AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr, AliTRDcluster *c)
 708 // {
 709 //   //
 710 //   // Iterative process to register clusters to the seed.
 711 //   // In iteration 0 we try only one pad-row and if quality not
 712 //   // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
 713 //   //
 714 //   // debug level 7
 715 //   //
 716 //
 717 //   if(!fReconstructor->GetRecoParam() ){
 718 //     AliError("Seed can not be used without a valid RecoParam.");
 719 //     return kFALSE;
 720 //   }
 721 //
 722 //   AliTRDchamberTimeBin *layer = 0x0;
 723 //   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7){
 724 //     AliTRDtrackingChamber ch(*chamber);
 725 //     ch.SetOwner();
 726 //     TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
 727 //     cstreamer << "AttachClustersIter"
 728 //       << "chamber.="   << &ch
 729 //       << "tracklet.="  << this
 730 //       << "\n";
 731 //   }
 732 //
 733 //   Float_t  tquality;
 734 //   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
 735 //   Double_t kroadz = fPadLength * .5 + 1.;
 736 //
 737 //   // initialize configuration parameters
 738 //   Float_t zcorr = kZcorr ? fTilt * (fZfit[0] - fZref[0]) : 0.;
 739 //   Int_t   niter = kZcorr ? 1 : 2;
 740 //
 741 //   Double_t yexp, zexp;
 742 //   Int_t ncl = 0;
 743 //   // start seed update
 744 //   for (Int_t iter = 0; iter < niter; iter++) {
 745 //     ncl = 0;
 746 //     for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 747 //       if(!(layer = chamber->GetTB(iTime))) continue;
 748 //       if(!Int_t(*layer)) continue;
 749 //
 750 //       // define searching configuration
 751 //       Double_t dxlayer = layer->GetX() - fX0;
 752 //       if(c){
 753 //         zexp = c->GetZ();
 754 //         //Try 2 pad-rows in second iteration
 755 //         if (iter > 0) {
 756 //           zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
 757 //           if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
 758 //           if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
 759 //         }
 760 //       } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
 761 //       yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
 762 //
 763 //       // Get and register cluster
 764 //       Int_t    index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
 765 //       if (index < 0) continue;
 766 //       AliTRDcluster *cl = (*layer)[index];
 767 //
 768 //       fIndexes[iTime]  = layer->GetGlobalIndex(index);
 769 //       fClusters[iTime] = cl;
 770 // //       fY[iTime]        = cl->GetY();
 771 // //       fZ[iTime]        = cl->GetZ();
 772 //       ncl++;
 773 //     }
 774 //     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fDet, ncl));
 775 //
 776 //     if(ncl>1){
 777 //       // calculate length of the time bin (calibration aware)
 778 //       Int_t irp = 0; Float_t x[2]={0., 0.}; Int_t tb[2] = {0,0};
 779 //       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 780 //         if(!fClusters[iTime]) continue;
 781 //         x[irp]  = fClusters[iTime]->GetX();
 782 //         tb[irp] = iTime;
 783 //         irp++;
 784 //         if(irp==2) break;
 785 //       }
 786 //       Int_t dtb = tb[1] - tb[0];
 787 //       fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
 788 //
 789 //       // update X0 from the clusters (calibration/alignment aware)
 790 //       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 791 //         if(!(layer = chamber->GetTB(iTime))) continue;
 792 //         if(!layer->IsT0()) continue;
 793 //         if(fClusters[iTime]){
 794 //           fX0 = fClusters[iTime]->GetX();
 795 //           break;
 796 //         } else { // we have to infere the position of the anode wire from the other clusters
 797 //           for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
 798 //             if(!fClusters[jTime]) continue;
 799 //             fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
 800 //             break;
 801 //           }
 802 //         }
 803 //       }
 804 //
 805 //       // update YZ reference point
 806 //       // TODO
 807 //
 808 //       // update x reference positions (calibration/alignment aware)
 809 // //       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 810 // //         if(!fClusters[iTime]) continue;
 811 // //         fX[iTime] = fX0 - fClusters[iTime]->GetX();
 812 // //       }
 813 //
 814 //       FitMI();
 815 //     }
 816 //     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fDet, fN2));
 817 //
 818 //     if(IsOK()){
 819 //       tquality = GetQuality(kZcorr);
 820 //       if(tquality < quality) break;
 821 //       else quality = tquality;
 822 //     }
 823 //     kroadz *= 2.;
 824 //   } // Loop: iter
 825 //   if (!IsOK()) return kFALSE;
 826 //
 827 //   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=1) CookLabels();
 828 //
 829 //   // load calibration params
 830 //   Calibrate();
 831 //   UpdateUsed();
 832 //   return kTRUE;
 833 // }
 834
 835 //____________________________________________________________________
 836 Bool_t  AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt)
 837 {
 838   //
 839   // Projective algorithm to attach clusters to seeding tracklets
 840   //
 841   // Parameters
 842   //
 843   // Output
 844   //
 845   // Detailed description
 846   // 1. Collapse x coordinate for the full detector plane
 847   // 2. truncated mean on y (r-phi) direction
 848   // 3. purge clusters
 849   // 4. truncated mean on z direction
 850   // 5. purge clusters
 851   // 6. fit tracklet
 852   //
 853   Bool_t kPRINT = kFALSE;
 854   if(!fReconstructor->GetRecoParam() ){
 855     AliError("Seed can not be used without a valid RecoParam.");
 856     return kFALSE;
 857   }
 858   // Initialize reco params for this tracklet
 859   // 1. first time bin in the drift region
 860   Int_t t0 = 4;
 861   Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
 862
 863   Double_t syRef  = TMath::Sqrt(fRefCov[0]);
 864   //define roads
 865   Double_t kroady = 1.;
 866   //fReconstructor->GetRecoParam() ->GetRoad1y();
 867   Double_t kroadz = fPadLength * 1.5 + 1.;
 868   if(kPRINT) printf("AttachClusters() sy[%f] road[%f]\n", syRef, kroady);
 869
 870   // working variables
 871   const Int_t kNrows = 16;
 872   AliTRDcluster *clst[kNrows][kNclusters];
 873   Double_t cond[4], dx, dy, yt, zt,
 874     yres[kNrows][kNclusters];
 875   Int_t idxs[kNrows][kNclusters], ncl[kNrows], ncls = 0;
 876   memset(ncl, 0, kNrows*sizeof(Int_t));
 877   memset(clst, 0, kNrows*kNclusters*sizeof(AliTRDcluster*));
 878
 879   // Do cluster projection
 880   AliTRDcluster *c = 0x0;
 881   AliTRDchamberTimeBin *layer = 0x0;
 882   Bool_t kBUFFER = kFALSE;
 883   for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
 884     if(!(layer = chamber->GetTB(it))) continue;
 885     if(!Int_t(*layer)) continue;
 886
 887     dx   = fX0 - layer->GetX();
 888     yt = fYref[0] - fYref[1] * dx;
 889     zt = fZref[0] - fZref[1] * dx;
 890     if(kPRINT) printf("\t%2d dx[%f] yt[%f] zt[%f]\n", it, dx, yt, zt);
 891
 892     // select clusters on a 5 sigmaKalman level
 893     cond[0] = yt; cond[2] = kroady;
 894     cond[1] = zt; cond[3] = kroadz;
 895     Int_t n=0, idx[6];
 896     layer->GetClusters(cond, idx, n, 6);
 897     for(Int_t ic = n; ic--;){
 898       c  = (*layer)[idx[ic]];
 899       dy = yt - c->GetY();
 900       dy += tilt ? fTilt * (c->GetZ() - zt) : 0.;
 901       // select clusters on a 3 sigmaKalman level
 902 /*      if(tilt && TMath::Abs(dy) > 3.*syRef){
 903         printf("too large !!!\n");
 904         continue;
 905       }*/
 906       Int_t r = c->GetPadRow();
 907       if(kPRINT) printf("\t\t%d dy[%f] yc[%f] r[%d]\n", ic, TMath::Abs(dy), c->GetY(), r);
 908       clst[r][ncl[r]] = c;
 909       idxs[r][ncl[r]] = idx[ic];
 910       yres[r][ncl[r]] = dy;
 911       ncl[r]++; ncls++;
 912
 913       if(ncl[r] >= kNclusters) {
 914         AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kNclusters));
 915         kBUFFER = kTRUE;
 916         break;
 917       }
 918     }
 919     if(kBUFFER) break;
 920   }
 921   if(kPRINT) printf("Found %d clusters\n", ncls);
 922   if(ncls<kClmin) return kFALSE;
 923
 924   // analyze each row individualy
 925   Double_t mean, syDis;
 926   Int_t nrow[] = {0, 0, 0}, nr = 0, lr=-1;
 927   for(Int_t ir=kNrows; ir--;){
 928     if(!(ncl[ir])) continue;
 929     if(lr>0 && lr-ir != 1){
 930       if(kPRINT) printf("W - gap in rows attached !!\n");
 931     }
 932     if(kPRINT) printf("\tir[%d] lr[%d] n[%d]\n", ir, lr, ncl[ir]);
 933     // Evaluate truncated mean on the y direction
 934     if(ncl[ir] > 3) AliMathBase::EvaluateUni(ncl[ir], yres[ir], mean, syDis, Int_t(ncl[ir]*.8));
 935     else {
 936       mean = 0.; syDis = 0.;
 937     }
 938
 939     // TODO check mean and sigma agains cluster resolution !!
 940     if(kPRINT) printf("\tr[%2d] m[%f %5.3fsigma] s[%f]\n", ir, mean, TMath::Abs(mean/syRef), syDis);
 941     // select clusters on a 3 sigmaDistr level
 942     Bool_t kFOUND = kFALSE;
 943     for(Int_t ic = ncl[ir]; ic--;){
 944       if(yres[ir][ic] - mean > 3. * syDis){
 945         clst[ir][ic] = 0x0; continue;
 946       }
 947       nrow[nr]++; kFOUND = kTRUE;
 948     }
 949     // exit loop
 950     if(kFOUND) nr++;
 951     lr = ir; if(nr>=3) break;
 952   }
 953   if(kPRINT) printf("lr[%d] nr[%d] nrow[0]=%d nrow[1]=%d nrow[2]=%d\n", lr, nr, nrow[0], nrow[1], nrow[2]);
 954
 955   // classify cluster rows
 956   Int_t row = -1;
 957   switch(nr){
 958   case 1:
 959     row = lr;
 960     break;
 961   case 2:
 962     SetBit(kRowCross, kTRUE); // mark pad row crossing
 963     if(nrow[0] > nrow[1]){ row = lr+1; lr = -1;}
 964     else{
 965       row = lr; lr = 1;
 966       nrow[2] = nrow[1];
 967       nrow[1] = nrow[0];
 968       nrow[0] = nrow[2];
 969     }
 970     break;
 971   case 3:
 972     SetBit(kRowCross, kTRUE); // mark pad row crossing
 973     break;
 974   }
 975   if(kPRINT) printf("\trow[%d] n[%d]\n\n", row, nrow[0]);
 976   if(row<0) return kFALSE;
 977
 978   // Select and store clusters
 979   // We should consider here :
 980   //  1. How far is the chamber boundary
 981   //  2. How big is the mean
 982   Int_t n = 0;
 983   for (Int_t ir = 0; ir < nr; ir++) {
 984     Int_t jr = row + ir*lr;
 985     if(kPRINT) printf("\tattach %d clusters for row %d\n", ncl[jr], jr);
 986     for (Int_t ic = 0; ic < ncl[jr]; ic++) {
 987       if(!(c = clst[jr][ic])) continue;
 988       Int_t it = c->GetPadTime();
 989       // TODO proper indexing of clusters !!
 990       fIndexes[it+kNtb*ir]  = chamber->GetTB(it)->GetGlobalIndex(idxs[jr][ic]);
 991       fClusters[it+kNtb*ir] = c;
 992
 993       //printf("\tid[%2d] it[%d] idx[%d]\n", ic, it, fIndexes[it]);
 994
 995       n++;
 996     }
 997   }
 998
 999   // number of minimum numbers of clusters expected for the tracklet
1000   if (n < kClmin){
1001     //AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", n, kClmin));
1002     return kFALSE;
1003   }
1004   SetN(n);
1005
1006   // Load calibration parameters for this tracklet
1007   Calibrate();
1008
1009   // calculate dx for time bins in the drift region (calibration aware)
1010   Int_t irp = 0; Float_t x[2] = {0.,0.}; Int_t tb[2]={0,0};
1011   for (Int_t it = t0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
1012     if(!fClusters[it]) continue;
1013     x[irp]  = fClusters[it]->GetX();
1014     tb[irp] = it;
1015     irp++;
1016     if(irp==2) break;
1017   }
1018   Int_t dtb = tb[1] - tb[0];
1019   fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
1020
1021   // update X0 from the clusters (calibration/alignment aware) TODO remove dependence on x0 !!
1022   for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
1023     if(!(layer = chamber->GetTB(it))) continue;
1024     if(!layer->IsT0()) continue;
1025     if(fClusters[it]){
1026       fX0 = fClusters[it]->GetX();
1027       break;
1028     } else { // we have to infere the position of the anode wire from the other clusters
1029       for (Int_t jt = it+1; jt < AliTRDtrackerV1::GetNTimeBins(); jt++) {
1030         if(!fClusters[jt]) continue;
1031         fX0 = fClusters[jt]->GetX() + fdX * (jt - it);
1032         break;
1033       }
1034     }
1035   }
1036
1037   return kTRUE;
1038 }
1039
1040 //____________________________________________________________
1041 void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
1042 {
1043 //   Fill in all derived information. It has to be called after recovery from file or HLT.
1044 //   The primitive data are
1045 //   - list of clusters
1046 //   - detector (as the detector will be removed from clusters)
1047 //   - position of anode wire (fX0) - temporary
1048 //   - track reference position and direction
1049 //   - momentum of the track
1050 //   - time bin length [cm]
1051 //
1052 //   A.Bercuci <A.Bercuci@gsi.de> Oct 30th 2008
1053 //
1054   fReconstructor = rec;
1055   AliTRDgeometry g;
1056   AliTRDpadPlane *pp = g.GetPadPlane(fDet);
1057   fTilt      = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
1058   fPadLength = pp->GetLengthIPad();
1059   //fSnp = fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);
1060   //fTgl = fZref[1];
1061   Int_t n = 0, nshare = 0, nused = 0;
1062   AliTRDcluster **cit = &fClusters[0];
1063   for(Int_t ic = kNclusters; ic--; cit++){
1064     if(!(*cit)) return;
1065     n++;
1066     if((*cit)->IsShared()) nshare++;
1067     if((*cit)->IsUsed()) nused++;
1068   }
1069   SetN(n); SetNUsed(nused); SetNShared(nshare);
1070   Fit();
1071   CookLabels();
1072   GetProbability();
1073 }
1074
1075
1076 //____________________________________________________________________
1077 Bool_t AliTRDseedV1::Fit(Bool_t tilt, Int_t errors)
1078 {
1079   //
1080   // Linear fit of the tracklet
1081   //
1082   // Parameters :
1083   //
1084   // Output :
1085   //  True if successful
1086   //
1087   // Detailed description
1088   // 2. Check if tracklet crosses pad row boundary
1089   // 1. Calculate residuals in the y (r-phi) direction
1090   // 3. Do a Least Square Fit to the data
1091   //
1092
1093   if(!IsCalibrated()){
1094     AliWarning("Tracklet fit failed. Call Calibrate().");
1095     return kFALSE;
1096   }
1097
1098   const Int_t kClmin = 8;
1099
1100
1101   // cluster error parametrization parameters
1102   // 1. sy total charge
1103   const Float_t sq0inv = 0.019962; // [1/q0]
1104   const Float_t sqb    = 1.0281564;    //[cm]
1105   // 2. sy for the PRF
1106   const Float_t scy[AliTRDgeometry::kNlayer][4] = {
1107     {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02},
1108     {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02},
1109     {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02},
1110     {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02},
1111     {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02},
1112     {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02},
1113   };
1114   // 3. sy parallel to the track
1115   const Float_t sy0 =  2.649e-02; // [cm]
1116   const Float_t sya = -8.864e-04; // [cm]
1117   const Float_t syb = -2.435e-01; // [cm]
1118
1119   // 4. sx parallel to the track
1120   const Float_t sxgc = 5.427e-02;
1121   const Float_t sxgm = 7.783e-01;
1122   const Float_t sxgs = 2.743e-01;
1123   const Float_t sxe0 =-2.065e+00;
1124   const Float_t sxe1 =-2.978e-02;
1125
1126   // 5. sx perpendicular to the track
1127 //   const Float_t sxd0 = 1.881e-02;
1128 //   const Float_t sxd1 =-4.101e-01;
1129 //   const Float_t sxd2 = 1.572e+00;
1130
1131   // get track direction
1132   Double_t y0   = fYref[0];
1133   Double_t dydx = fYref[1];
1134   Double_t z0   = fZref[0];
1135   Double_t dzdx = fZref[1];
1136   Double_t yt, zt;
1137
1138   // calculation of tg^2(phi - a_L) and tg^2(a_L)
1139   Double_t tgg = (dydx-fExB)/(1.+dydx*fExB); tgg *= tgg;
1140   //Double_t exb2= fExB*fExB;
1141
1142   //AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
1143   TLinearFitter  fitterY(1, "pol1");
1144   // convertion factor from square to gauss distribution for sigma
1145   //Double_t convert = 1./TMath::Sqrt(12.);
1146
1147   // book cluster information
1148   Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
1149
1150   Int_t ily = AliTRDgeometry::GetLayer(fDet);
1151   Int_t fN = 0;
1152   AliTRDcluster *c=0x0, **jc = &fClusters[0];
1153   for (Int_t ic=0; ic<kNtb; ic++, ++jc) {
1154     //zRow[ic] = -1;
1155     xc[ic]  = -1.;
1156     yc[ic]  = 999.;
1157     zc[ic]  = 999.;
1158     sy[ic]  = 0.;
1159     //sz[ic]  = 0.;
1160     if(!(c = (*jc))) continue;
1161     if(!c->IsInChamber()) continue;
1162
1163     Float_t w = 1.;
1164     if(c->GetNPads()>4) w = .5;
1165     if(c->GetNPads()>5) w = .2;
1166
1167     //zRow[fN] = c->GetPadRow();
1168     qc[fN]   = TMath::Abs(c->GetQ());
1169     // correct cluster position for PRF and v drift
1170     //Int_t jc = TMath::Max(fN-3, 0);
1171     //xc[fN]   = c->GetXloc(fT0, fVD, &qc[jc], &xc[jc]/*, z0 - c->GetX()*dzdx*/);
1172     //Double_t s2 = fS2PRF + fDiffL*fDiffL*xc[fN]/(1.+2.*exb2)+tgg*xc[fN]*xc[fN]*exb2/12.;
1173     //yc[fN]   = c->GetYloc(s2, fPadLength, xc[fN], fExB);
1174
1175     // uncalibrated cluster correction
1176     // TODO remove
1177     Double_t x, y; GetClusterXY(c, x, y);
1178     xc[fN]   = fX0 - x;
1179     yc[fN]   = y;
1180     zc[fN]   = c->GetZ();
1181
1182     // extrapolated y value for the track
1183     yt = y0 - xc[fN]*dydx;
1184     // extrapolated z value for the track
1185     zt = z0 - xc[fN]*dzdx;
1186     // tilt correction
1187     if(tilt) yc[fN] -= fTilt*(zc[fN] - zt);
1188
1189     // ELABORATE CLUSTER ERROR
1190     // TODO to be moved to AliTRDcluster
1191     // basic y error (|| to track).
1192     sy[fN]  = xc[fN] < AliTRDgeometry::CamHght() ? 2. : sy0 + sya*TMath::Exp(1./(xc[fN]+syb));
1193     //printf("cluster[%d]\n\tsy[0] = %5.3e [um]\n", fN,  sy[fN]*1.e4);
1194     // y error due to total charge
1195     sy[fN] += sqb*(1./qc[fN] - sq0inv);
1196     //printf("\tsy[1] = %5.3e [um]\n", sy[fN]*1.e4);
1197     // y error due to PRF
1198     sy[fN] += scy[ily][0]*TMath::Gaus(c->GetCenter(), scy[ily][1], scy[ily][2]) - scy[ily][3];
1199     //printf("\tsy[2] = %5.3e [um]\n", sy[fN]*1.e4);
1200
1201     sy[fN] *= sy[fN];
1202
1203     // ADD ERROR ON x
1204     // error of drift length parallel to the track
1205     Double_t sx = sxgc*TMath::Gaus(xc[fN], sxgm, sxgs) + TMath::Exp(sxe0+sxe1*xc[fN]); // [cm]
1206     //printf("\tsx[0] = %5.3e [um]\n", sx*1.e4);
1207     // error of drift length perpendicular to the track
1208     //sx += sxd0 + sxd1*d + sxd2*d*d;
1209     sx *= sx; // square sx
1210
1211     // add error from ExB
1212     if(errors>0) sy[fN] += fExB*fExB*sx;
1213     //printf("\tsy[3] = %5.3e [um^2]\n", sy[fN]*1.e8);
1214
1215     // global radial error due to misalignment/miscalibration
1216     Double_t sx0  = 0.; sx0 *= sx0;
1217     // add sx contribution to sy due to track angle
1218     if(errors>1) sy[fN] += tgg*(sx+sx0);
1219     // TODO we should add tilt pad correction here
1220     //printf("\tsy[4] = %5.3e [um^2]\n", sy[fN]*1.e8);
1221     c->SetSigmaY2(sy[fN]);
1222
1223     sy[fN]  = TMath::Sqrt(sy[fN]);
1224     fitterY.AddPoint(&xc[fN], yc[fN], sy[fN]);
1225     fN++;
1226   }
1227   // to few clusters
1228   if (fN < kClmin) return kFALSE;
1229
1230   // fit XY
1231   fitterY.Eval();
1232   fYfit[0] = fitterY.GetParameter(0);
1233   fYfit[1] = -fitterY.GetParameter(1);
1234   // store covariance
1235   Double_t *p = fitterY.GetCovarianceMatrix();
1236   fCov[0] = p[0]; // variance of y0
1237   fCov[1] = p[1]; // covariance of y0, dydx
1238   fCov[2] = p[3]; // variance of dydx
1239   // the ref radial position is set at the minimum of
1240   // the y variance of the tracklet
1241   fX   = -fCov[1]/fCov[2]; //fXref = fX0 - fXref;
1242   fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
1243
1244   // fit XZ
1245   if(IsRowCross()){
1246     // TODO pad row cross position estimation !!!
1247     //AliInfo(Form("Padrow cross in detector %d", fDet));
1248     fZfit[0] = .5*(zc[0]+zc[fN-1]); fZfit[1] = 0.;
1249     fS2Z     = 0.02+1.55*fZref[1]; fS2Z *= fS2Z;
1250   } else {
1251     fZfit[0] = zc[0]; fZfit[1] = 0.;
1252     fS2Z     = fPadLength*fPadLength/12.;
1253   }
1254
1255
1256 //   // determine z offset of the fit
1257 //   Float_t zslope = 0.;
1258 //   Int_t nchanges = 0, nCross = 0;
1259 //   if(nz==2){ // tracklet is crossing pad row
1260 //     // Find the break time allowing one chage on pad-rows
1261 //     // with maximal number of accepted clusters
1262 //     Int_t padRef = zRow[0];
1263 //     for (Int_t ic=1; ic<fN; ic++) {
1264 //       if(zRow[ic] == padRef) continue;
1265 //
1266 //       // debug
1267 //       if(zRow[ic-1] == zRow[ic]){
1268 //         printf("ERROR in pad row change!!!\n");
1269 //       }
1270 //
1271 //       // evaluate parameters of the crossing point
1272 //       Float_t sx = (xc[ic-1] - xc[ic])*convert;
1273 //       fCross[0] = .5 * (xc[ic-1] + xc[ic]);
1274 //       fCross[2] = .5 * (zc[ic-1] + zc[ic]);
1275 //       fCross[3] = TMath::Max(dzdx * sx, .01);
1276 //       zslope    = zc[ic-1] > zc[ic] ? 1. : -1.;
1277 //       padRef    = zRow[ic];
1278 //       nCross    = ic;
1279 //       nchanges++;
1280 //     }
1281 //   }
1282 //
1283 //   // condition on nCross and reset nchanges TODO
1284 //
1285 //   if(nchanges==1){
1286 //     if(dzdx * zslope < 0.){
1287 //       AliInfo("Tracklet-Track mismatch in dzdx. TODO.");
1288 //     }
1289 //
1290 //
1291 //     //zc[nc] = fitterZ.GetFunctionParameter(0);
1292 //     fCross[1] = fYfit[0] - fCross[0] * fYfit[1];
1293 //     fCross[0] = fX0 - fCross[0];
1294 //   }
1295
1296   return kTRUE;
1297 }
1298
1299
1300 /*
1301 //_____________________________________________________________________________
1302 void AliTRDseedV1::FitMI()
1303 {
1304 //
1305 // Fit the seed.
1306 // Marian Ivanov's version
1307 //
1308 // linear fit on the y direction with respect to the reference direction.
1309 // The residuals for each x (x = xc - x0) are deduced from:
1310 // dy = y - yt             (1)
1311 // the tilting correction is written :
1312 // y = yc + h*(zc-zt)      (2)
1313 // yt = y0+dy/dx*x         (3)
1314 // zt = z0+dz/dx*x         (4)
1315 // from (1),(2),(3) and (4)
1316 // dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0)
1317 // the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this:
1318 // 1. use tilting correction for calculating the y
1319 // 2. neglect tilting correction here and account for it in the error parametrization of the tracklet.
1320   const Float_t kRatio  = 0.8;
1321   const Int_t   kClmin  = 5;
1322   const Float_t kmaxtan = 2;
1323
1324   if (TMath::Abs(fYref[1]) > kmaxtan){
1325                 //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
1326                 return;              // Track inclined too much
1327         }
1328
1329   Float_t  sigmaexp  = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
1330   Float_t  ycrosscor = fPadLength * fTilt * 0.5;           // Y correction for crossing
1331   Int_t fNChange = 0;
1332
1333   Double_t sumw;
1334   Double_t sumwx;
1335   Double_t sumwx2;
1336   Double_t sumwy;
1337   Double_t sumwxy;
1338   Double_t sumwz;
1339   Double_t sumwxz;
1340
1341         // Buffering: Leave it constant fot Performance issues
1342   Int_t    zints[kNtb];            // Histograming of the z coordinate
1343                                          // Get 1 and second max probable coodinates in z
1344   Int_t    zouts[2*kNtb];
1345   Float_t  allowedz[kNtb];         // Allowed z for given time bin
1346   Float_t  yres[kNtb];             // Residuals from reference
1347   //Float_t  anglecor = fTilt * fZref[1];  // Correction to the angle
1348
1349   Float_t pos[3*kNtb]; memset(pos, 0, 3*kNtb*sizeof(Float_t));
1350   Float_t *fX = &pos[0], *fY = &pos[kNtb], *fZ = &pos[2*kNtb];
1351
1352   Int_t fN  = 0; AliTRDcluster *c = 0x0;
1353   fN2 = 0;
1354   for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1355     yres[i] = 10000.0;
1356     if (!(c = fClusters[i])) continue;
1357     if(!c->IsInChamber()) continue;
1358     // Residual y
1359     //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + fTilt*(fZ[i] - fZref[0]);
1360     fX[i] = fX0 - c->GetX();
1361     fY[i] = c->GetY();
1362     fZ[i] = c->GetZ();
1363     yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
1364     zints[fN] = Int_t(fZ[i]);
1365     fN++;
1366   }
1367
1368   if (fN < kClmin){
1369     //printf("Exit fN < kClmin: fN = %d\n", fN);
1370     return;
1371   }
1372   Int_t nz = AliTRDtrackerV1::Freq(fN, zints, zouts, kFALSE);
1373   Float_t fZProb   = zouts[0];
1374   if (nz <= 1) zouts[3] = 0;
1375   if (zouts[1] + zouts[3] < kClmin) {
1376     //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
1377     return;
1378   }
1379
1380   // Z distance bigger than pad - length
1381   if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
1382
1383   Int_t  breaktime = -1;
1384   Bool_t mbefore   = kFALSE;
1385   Int_t  cumul[kNtb][2];
1386   Int_t  counts[2] = { 0, 0 };
1387
1388   if (zouts[3] >= 3) {
1389
1390     //
1391     // Find the break time allowing one chage on pad-rows
1392     // with maximal number of accepted clusters
1393     //
1394     fNChange = 1;
1395     for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1396       cumul[i][0] = counts[0];
1397       cumul[i][1] = counts[1];
1398       if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
1399       if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
1400     }
1401     Int_t  maxcount = 0;
1402     for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
1403       Int_t after  = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
1404       Int_t before = cumul[i][1];
1405       if (after + before > maxcount) {
1406         maxcount  = after + before;
1407         breaktime = i;
1408         mbefore   = kFALSE;
1409       }
1410       after  = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
1411       before = cumul[i][0];
1412       if (after + before > maxcount) {
1413         maxcount  = after + before;
1414         breaktime = i;
1415         mbefore   = kTRUE;
1416       }
1417     }
1418     breaktime -= 1;
1419   }
1420
1421   for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1422     if (i >  breaktime) allowedz[i] =   mbefore  ? zouts[2] : zouts[0];
1423     if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
1424   }
1425
1426   if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
1427       ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
1428     //
1429     // Tracklet z-direction not in correspondance with track z direction
1430     //
1431     fNChange = 0;
1432     for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1433       allowedz[i] = zouts[0];  // Only longest taken
1434     }
1435   }
1436
1437   if (fNChange > 0) {
1438     //
1439     // Cross pad -row tracklet  - take the step change into account
1440     //
1441     for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1442       if (!fClusters[i]) continue;
1443       if(!fClusters[i]->IsInChamber()) continue;
1444       if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1445       // Residual y
1446       //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + fTilt*(fZ[i] - fZref[0]);
1447       yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
1448 //       if (TMath::Abs(fZ[i] - fZProb) > 2) {
1449 //         if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength;
1450 //         if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength;
1451       }
1452     }
1453   }
1454
1455   Double_t yres2[kNtb];
1456   Double_t mean;
1457   Double_t sigma;
1458   for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1459     if (!fClusters[i]) continue;
1460     if(!fClusters[i]->IsInChamber()) continue;
1461     if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
1462     yres2[fN2] = yres[i];
1463     fN2++;
1464   }
1465   if (fN2 < kClmin) {
1466                 //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
1467     fN2 = 0;
1468     return;
1469   }
1470   AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
1471   if (sigma < sigmaexp * 0.8) {
1472     sigma = sigmaexp;
1473   }
1474   //Float_t fSigmaY = sigma;
1475
1476   // Reset sums
1477   sumw   = 0;
1478   sumwx  = 0;
1479   sumwx2 = 0;
1480   sumwy  = 0;
1481   sumwxy = 0;
1482   sumwz  = 0;
1483   sumwxz = 0;
1484
1485   fN2    = 0;
1486   Float_t fMeanz = 0;
1487   Float_t fMPads = 0;
1488   fUsable = 0;
1489   for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1490     if (!fClusters[i]) continue;
1491     if (!fClusters[i]->IsInChamber()) continue;
1492     if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
1493     if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0;  continue;}
1494     SETBIT(fUsable,i);
1495     fN2++;
1496     fMPads += fClusters[i]->GetNPads();
1497     Float_t weight = 1.0;
1498     if (fClusters[i]->GetNPads() > 4) weight = 0.5;
1499     if (fClusters[i]->GetNPads() > 5) weight = 0.2;
1500
1501
1502     Double_t x = fX[i];
1503     //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
1504
1505     sumw   += weight;
1506     sumwx  += x * weight;
1507     sumwx2 += x*x * weight;
1508     sumwy  += weight * yres[i];
1509     sumwxy += weight * (yres[i]) * x;
1510     sumwz  += weight * fZ[i];
1511     sumwxz += weight * fZ[i] * x;
1512
1513   }
1514
1515   if (fN2 < kClmin){
1516                 //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
1517     fN2 = 0;
1518     return;
1519   }
1520   fMeanz = sumwz / sumw;
1521   Float_t correction = 0;
1522   if (fNChange > 0) {
1523     // Tracklet on boundary
1524     if (fMeanz < fZProb) correction =  ycrosscor;
1525     if (fMeanz > fZProb) correction = -ycrosscor;
1526   }
1527
1528   Double_t det = sumw * sumwx2 - sumwx * sumwx;
1529   fYfit[0]    = (sumwx2 * sumwy  - sumwx * sumwxy) / det;
1530   fYfit[1]    = (sumw   * sumwxy - sumwx * sumwy)  / det;
1531
1532   fS2Y = 0;
1533   for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
1534     if (!TESTBIT(fUsable,i)) continue;
1535     Float_t delta = yres[i] - fYfit[0] - fYfit[1] * fX[i];
1536     fS2Y += delta*delta;
1537   }
1538   fS2Y = TMath::Sqrt(fS2Y / Float_t(fN2-2));
1539         // TEMPORARY UNTIL covariance properly calculated
1540         fS2Y = TMath::Max(fS2Y, Float_t(.1));
1541
1542   fZfit[0]   = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
1543   fZfit[1]   = (sumw   * sumwxz - sumwx * sumwz)  / det;
1544 //   fYfitR[0] += fYref[0] + correction;
1545 //   fYfitR[1] += fYref[1];
1546 //  fYfit[0]   = fYfitR[0];
1547   fYfit[1]   = -fYfit[1];
1548
1549   UpdateUsed();
1550 }*/
1551
1552 //___________________________________________________________________
1553 void AliTRDseedV1::Print(Option_t *o) const
1554 {
1555   //
1556   // Printing the seedstatus
1557   //
1558
1559   AliInfo(Form("Det[%3d] Tilt[%+6.2f] Pad[%5.2f]", fDet, fTilt, fPadLength));
1560   AliInfo(Form("N[%2d] Nused[%2d] Nshared[%2d]", GetN(), GetNUsed(), GetNShared()));
1561   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]));
1562   AliInfo(Form("Ref        y[%7.2f] z[%7.2f] dydx[%5.2f] dzdx[%5.2f]", fYref[0], fZref[0], fYref[1], fZref[1]))
1563
1564
1565   if(strcmp(o, "a")!=0) return;
1566
1567   AliTRDcluster* const* jc = &fClusters[0];
1568   for(int ic=0; ic<kNclusters; ic++, jc++) {
1569     if(!(*jc)) continue;
1570     (*jc)->Print(o);
1571   }
1572 }
1573
1574
1575 //___________________________________________________________________
1576 Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
1577 {
1578   // Checks if current instance of the class has the same essential members
1579   // as the given one
1580
1581   if(!o) return kFALSE;
1582   const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
1583   if(!inTracklet) return kFALSE;
1584
1585   for (Int_t i = 0; i < 2; i++){
1586     if ( fYref[i] != inTracklet->fYref[i] ) return kFALSE;
1587     if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
1588   }
1589
1590   if ( fS2Y != inTracklet->fS2Y ) return kFALSE;
1591   if ( fTilt != inTracklet->fTilt ) return kFALSE;
1592   if ( fPadLength != inTracklet->fPadLength ) return kFALSE;
1593
1594   for (Int_t i = 0; i < kNclusters; i++){
1595 //     if ( fX[i] != inTracklet->GetX(i) ) return kFALSE;
1596 //     if ( fY[i] != inTracklet->GetY(i) ) return kFALSE;
1597 //     if ( fZ[i] != inTracklet->GetZ(i) ) return kFALSE;
1598     if ( fIndexes[i] != inTracklet->fIndexes[i] ) return kFALSE;
1599   }
1600 //   if ( fUsable != inTracklet->fUsable ) return kFALSE;
1601
1602   for (Int_t i=0; i < 2; i++){
1603     if ( fYfit[i] != inTracklet->fYfit[i] ) return kFALSE;
1604     if ( fZfit[i] != inTracklet->fZfit[i] ) return kFALSE;
1605     if ( fLabels[i] != inTracklet->fLabels[i] ) return kFALSE;
1606   }
1607
1608 /*  if ( fMeanz != inTracklet->GetMeanz() ) return kFALSE;
1609   if ( fZProb != inTracklet->GetZProb() ) return kFALSE;*/
1610   if ( fN != inTracklet->fN ) return kFALSE;
1611   //if ( fNUsed != inTracklet->fNUsed ) return kFALSE;
1612   //if ( fFreq != inTracklet->GetFreq() ) return kFALSE;
1613   //if ( fNChange != inTracklet->GetNChange() ) return kFALSE;
1614
1615   if ( fC != inTracklet->fC ) return kFALSE;
1616   //if ( fCC != inTracklet->GetCC() ) return kFALSE;
1617   if ( fChi2 != inTracklet->fChi2 ) return kFALSE;
1618   //  if ( fChi2Z != inTracklet->GetChi2Z() ) return kFALSE;
1619
1620   if ( fDet != inTracklet->fDet ) return kFALSE;
1621   if ( fMom != inTracklet->fMom ) return kFALSE;
1622   if ( fdX != inTracklet->fdX ) return kFALSE;
1623
1624   for (Int_t iCluster = 0; iCluster < kNclusters; iCluster++){
1625     AliTRDcluster *curCluster = fClusters[iCluster];
1626     AliTRDcluster *inCluster = inTracklet->fClusters[iCluster];
1627     if (curCluster && inCluster){
1628       if (! curCluster->IsEqual(inCluster) ) {
1629         curCluster->Print();
1630         inCluster->Print();
1631         return kFALSE;
1632       }
1633     } else {
1634       // if one cluster exists, and corresponding
1635       // in other tracklet doesn't - return kFALSE
1636       if(curCluster || inCluster) return kFALSE;
1637     }
1638   }
1639   return kTRUE;
1640 }