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