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