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
  36 #include "AliTRDcluster.h"
  37 #include "AliTRDseedV1.h"
  38 #include "AliTRDtrackV1.h"
  39 #include "AliTRDcalibDB.h"
  40 #include "AliTRDchamberTimeBin.h"
  41 #include "AliTRDtrackingChamber.h"
  42 #include "AliTRDtrackerV1.h"
  43 #include "AliTRDReconstructor.h"
  44 #include "AliTRDrecoParam.h"
  45 #include "AliTRDgeometry.h"
  46 #include "Cal/AliTRDCalPID.h"
  47
  48 ClassImp(AliTRDseedV1)
  49
  50 //____________________________________________________________________
  51 AliTRDseedV1::AliTRDseedV1(Int_t plane)
  52   :AliTRDseed()
  53   ,fReconstructor(0x0)
  54   ,fPlane(plane)
  55   ,fMom(0.)
  56   ,fSnp(0.)
  57   ,fTgl(0.)
  58   ,fdX(0.)
  59 {
  60   //
  61   // Constructor
  62   //
  63   //printf("AliTRDseedV1::AliTRDseedV1()\n");
  64
  65   for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = 0.;
  66   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
  67 }
  68
  69 //____________________________________________________________________
  70 AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
  71   :AliTRDseed((AliTRDseed&)ref)
  72   ,fReconstructor(ref.fReconstructor)
  73   ,fPlane(ref.fPlane)
  74   ,fMom(ref.fMom)
  75   ,fSnp(ref.fSnp)
  76   ,fTgl(ref.fTgl)
  77   ,fdX(ref.fdX)
  78 {
  79   //
  80   // Copy Constructor performing a deep copy
  81   //
  82
  83   //printf("AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &)\n");
  84   SetBit(kOwner, kFALSE);
  85   for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice];
  86   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec];
  87 }
  88
  89
  90 //____________________________________________________________________
  91 AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
  92 {
  93   //
  94   // Assignment Operator using the copy function
  95   //
  96
  97   if(this != &ref){
  98     ref.Copy(*this);
  99   }
 100   return *this;
 101
 102 }
 103
 104 //____________________________________________________________________
 105 AliTRDseedV1::~AliTRDseedV1()
 106 {
 107   //
 108   // Destructor. The RecoParam object belongs to the underlying tracker.
 109   //
 110
 111   //printf("I-AliTRDseedV1::~AliTRDseedV1() : Owner[%s]\n", IsOwner()?"YES":"NO");
 112
 113   if(IsOwner())
 114     for(int itb=0; itb<knTimebins; itb++){
 115       if(!fClusters[itb]) continue;
 116       //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
 117       delete fClusters[itb];
 118       fClusters[itb] = 0x0;
 119     }
 120 }
 121
 122 //____________________________________________________________________
 123 void AliTRDseedV1::Copy(TObject &ref) const
 124 {
 125   //
 126   // Copy function
 127   //
 128
 129   //AliInfo("");
 130   AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
 131
 132   target.fPlane         = fPlane;
 133   target.fMom           = fMom;
 134   target.fSnp           = fSnp;
 135   target.fTgl           = fTgl;
 136   target.fdX            = fdX;
 137   target.fReconstructor = fReconstructor;
 138
 139   for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice];
 140   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec];
 141
 142   AliTRDseed::Copy(target);
 143 }
 144
 145
 146 //____________________________________________________________
 147 Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
 148 {
 149 // Initialize this tracklet using the track information
 150 //
 151 // Parameters:
 152 //   track - the TRD track used to initialize the tracklet
 153 //
 154 // Detailed description
 155 // The function sets the starting point and direction of the
 156 // tracklet according to the information from the TRD track.
 157 //
 158 // Caution
 159 // The TRD track has to be propagated to the beginning of the
 160 // chamber where the tracklet will be constructed
 161 //
 162
 163   Double_t y, z;
 164   if(!track->GetProlongation(fX0, y, z)) return kFALSE;
 165   fYref[0] = y;
 166   fYref[1] = track->GetSnp()/(1. - track->GetSnp()*track->GetSnp());
 167   fZref[0] = z;
 168   fZref[1] = track->GetTgl();
 169
 170   //printf("Tracklet ref x[%7.3f] y[%7.3f] z[%7.3f], snp[%f] tgl[%f]\n", fX0, fYref[0], fZref[0], track->GetSnp(), track->GetTgl());
 171   return kTRUE;
 172 }
 173
 174
 175 //____________________________________________________________________
 176 void AliTRDseedV1::CookdEdx(Int_t nslices)
 177 {
 178 // Calculates average dE/dx for all slices and store them in the internal array fdEdx.
 179 //
 180 // Parameters:
 181 //  nslices : number of slices for which dE/dx should be calculated
 182 // Output:
 183 //  store results in the internal array fdEdx. This can be accessed with the method
 184 //  AliTRDseedV1::GetdEdx()
 185 //
 186 // Detailed description
 187 // Calculates average dE/dx for all slices. Depending on the PID methode
 188 // the number of slices can be 3 (LQ) or 8(NN).
 189 // The calculation of dQ/dl are done using the tracklet fit results (see AliTRDseedV1::GetdQdl(Int_t)) i.e.
 190 //
 191 // dQ/dl = qc/(dx * sqrt(1 + dy/dx^2 + dz/dx^2))
 192 //
 193 // The following effects are included in the calculation:
 194 // 1. calibration values for t0 and vdrift (using x coordinate to calculate slice)
 195 // 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
 196 // 3. cluster size
 197 //
 198
 199   Int_t nclusters[knSlices];
 200   for(int i=0; i<knSlices; i++){
 201     fdEdx[i]     = 0.;
 202     nclusters[i] = 0;
 203   }
 204   Float_t clength = (/*.5 * */AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
 205
 206   AliTRDcluster *cluster = 0x0;
 207   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
 208     if(!(cluster = fClusters[ic])) continue;
 209     Float_t x = cluster->GetX();
 210
 211     // Filter clusters for dE/dx calculation
 212
 213     // 1.consider calibration effects for slice determination
 214     Int_t slice;
 215     if(cluster->IsInChamber()) slice = Int_t(TMath::Abs(fX0 - x) * nslices / clength);
 216     else slice = x < fX0 ? 0 : nslices-1;
 217
 218     // 2. take sharing into account
 219     Float_t w = cluster->IsShared() ? .5 : 1.;
 220
 221     // 3. take into account large clusters TODO
 222     //w *= c->GetNPads() > 3 ? .8 : 1.;
 223
 224     //CHECK !!!
 225     fdEdx[slice]   += w * GetdQdl(ic); //fdQdl[ic];
 226     nclusters[slice]++;
 227   } // End of loop over clusters
 228
 229   if(fReconstructor->GetPIDMethod() == AliTRDReconstructor::kLQPID){
 230   // calculate mean charge per slice (only LQ PID)
 231     for(int is=0; is<nslices; is++){
 232       if(nclusters[is]) fdEdx[is] /= nclusters[is];
 233     }
 234   }
 235 }
 236
 237
 238 //____________________________________________________________________
 239 Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
 240 {
 241   return fClusters[ic] ? TMath::Abs(fClusters[ic]->GetQ()) /fdX / TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]) : 0.;
 242 }
 243
 244 //____________________________________________________________________
 245 Double_t* AliTRDseedV1::GetProbability()
 246 {
 247 // Fill probability array for tracklet from the DB.
 248 //
 249 // Parameters
 250 //
 251 // Output
 252 //   returns pointer to the probability array and 0x0 if missing DB access
 253 //
 254 // Detailed description
 255
 256
 257   // retrive calibration db
 258   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 259   if (!calibration) {
 260     AliError("No access to calibration data");
 261     return 0x0;
 262   }
 263
 264   if (!fReconstructor) {
 265     AliError("Reconstructor not set.");
 266     return 0x0;
 267   }
 268
 269   // Retrieve the CDB container class with the parametric detector response
 270   const AliTRDCalPID *pd = calibration->GetPIDObject(fReconstructor->GetPIDMethod());
 271   if (!pd) {
 272     AliError("No access to AliTRDCalPID object");
 273     return 0x0;
 274   }
 275   //AliInfo(Form("Method[%d] : %s", fReconstructor->GetRecoParam() ->GetPIDMethod(), pd->IsA()->GetName()));
 276
 277   // calculate tracklet length TO DO
 278   Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
 279   /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
 280
 281   //calculate dE/dx
 282   CookdEdx(fReconstructor->GetNdEdxSlices());
 283
 284   // Sets the a priori probabilities
 285   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
 286     fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, fPlane);
 287   }
 288
 289   return &fProb[0];
 290 }
 291
 292 //____________________________________________________________________
 293 Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
 294 {
 295   //
 296   // Returns a quality measurement of the current seed
 297   //
 298
 299   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 300   return
 301       .5 * TMath::Abs(18.0 - fN2)
 302     + 10.* TMath::Abs(fYfit[1] - fYref[1])
 303     + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
 304     + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
 305 }
 306
 307 //____________________________________________________________________
 308 void AliTRDseedV1::GetCovAt(Double_t /*x*/, Double_t *cov) const
 309 {
 310 // Computes covariance in the y-z plane at radial point x
 311
 312   Int_t ic = 0; while (!fClusters[ic]) ic++;
 313   AliTRDcalibDB *fCalib = AliTRDcalibDB::Instance();
 314   Double_t exB         = fCalib->GetOmegaTau(fCalib->GetVdriftAverage(fClusters[ic]->GetDetector()), -AliTracker::GetBz()*0.1);
 315
 316   Double_t sy2    = fSigmaY2*fSigmaY2 + .2*(fYfit[1]-exB)*(fYfit[1]-exB);
 317   Double_t sz2    = fPadLength/12.;
 318
 319
 320   //printf("Yfit[1] %f sy20 %f SigmaY2 %f\n", fYfit[1], sy20, fSigmaY2);
 321
 322   cov[0] = sy2;
 323   cov[1] = fTilt*(sy2-sz2);
 324   cov[2] = sz2;
 325
 326   // insert systematic uncertainties calibration and misalignment
 327   Double_t sys[15];
 328   fReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
 329   cov[0] += (sys[0]*sys[0]);
 330   cov[2] += (sys[1]*sys[1]);
 331 }
 332
 333
 334 //____________________________________________________________________
 335 void AliTRDseedV1::SetOwner()
 336 {
 337   //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
 338
 339   if(TestBit(kOwner)) return;
 340   for(int ic=0; ic<knTimebins; ic++){
 341     if(!fClusters[ic]) continue;
 342     fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
 343   }
 344   SetBit(kOwner);
 345 }
 346
 347 //____________________________________________________________________
 348 Bool_t  AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr, AliTRDcluster *c)
 349 {
 350   //
 351   // Iterative process to register clusters to the seed.
 352   // In iteration 0 we try only one pad-row and if quality not
 353   // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
 354   //
 355   // debug level 7
 356   //
 357
 358   if(!fReconstructor->GetRecoParam() ){
 359     AliError("Seed can not be used without a valid RecoParam.");
 360     return kFALSE;
 361   }
 362
 363   AliTRDchamberTimeBin *layer = 0x0;
 364   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7 && c){
 365     TClonesArray clusters("AliTRDcluster", 24);
 366     clusters.SetOwner(kTRUE);
 367     AliTRDcluster *cc = 0x0;
 368     Int_t det=-1, ncl, ncls = 0;
 369     for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 370       if(!(layer = chamber->GetTB(iTime))) continue;
 371       if(!(ncl = Int_t(*layer))) continue;
 372       for(int ic=0; ic<ncl; ic++){
 373         cc = (*layer)[ic];
 374         det = cc->GetDetector();
 375         new(clusters[ncls++]) AliTRDcluster(*cc);
 376       }
 377     }
 378     AliInfo(Form("N clusters[%d] = %d", fPlane, ncls));
 379
 380     Int_t ref = c ? 1 : 0;
 381     TTreeSRedirector &cstreamer = *AliTRDtrackerV1::DebugStreamer();
 382     cstreamer << "AttachClustersIter"
 383       << "det="        << det
 384       << "ref="        << ref
 385       << "clusters.="  << &clusters
 386       << "tracklet.="  << this
 387       << "cl.="        << c
 388       << "\n";
 389   }
 390
 391   Float_t  tquality;
 392   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
 393   Double_t kroadz = fPadLength * .5 + 1.;
 394
 395   // initialize configuration parameters
 396   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 397   Int_t   niter = kZcorr ? 1 : 2;
 398
 399   Double_t yexp, zexp;
 400   Int_t ncl = 0;
 401   // start seed update
 402   for (Int_t iter = 0; iter < niter; iter++) {
 403     ncl = 0;
 404     for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 405       if(!(layer = chamber->GetTB(iTime))) continue;
 406       if(!Int_t(*layer)) continue;
 407
 408       // define searching configuration
 409       Double_t dxlayer = layer->GetX() - fX0;
 410       if(c){
 411         zexp = c->GetZ();
 412         //Try 2 pad-rows in second iteration
 413         if (iter > 0) {
 414           zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
 415           if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
 416           if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
 417         }
 418       } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
 419       yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
 420
 421       // Get and register cluster
 422       Int_t    index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
 423       if (index < 0) continue;
 424       AliTRDcluster *cl = (*layer)[index];
 425
 426       fIndexes[iTime]  = layer->GetGlobalIndex(index);
 427       fClusters[iTime] = cl;
 428       fY[iTime]        = cl->GetY();
 429       fZ[iTime]        = cl->GetZ();
 430       ncl++;
 431     }
 432     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fPlane, ncl));
 433
 434     if(ncl>1){
 435       // calculate length of the time bin (calibration aware)
 436       Int_t irp = 0; Float_t x[2]; Int_t tb[2];
 437       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 438         if(!fClusters[iTime]) continue;
 439         x[irp]  = fClusters[iTime]->GetX();
 440         tb[irp] = iTime;
 441         irp++;
 442         if(irp==2) break;
 443       }
 444       fdX = (x[1] - x[0]) / (tb[0] - tb[1]);
 445
 446       // update X0 from the clusters (calibration/alignment aware)
 447       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 448         if(!(layer = chamber->GetTB(iTime))) continue;
 449         if(!layer->IsT0()) continue;
 450         if(fClusters[iTime]){
 451           fX0 = fClusters[iTime]->GetX();
 452           break;
 453         } else { // we have to infere the position of the anode wire from the other clusters
 454           for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
 455             if(!fClusters[jTime]) continue;
 456             fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
 457           }
 458           break;
 459         }
 460       }
 461
 462       // update YZ reference point
 463       // TODO
 464
 465       // update x reference positions (calibration/alignment aware)
 466       for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 467         if(!fClusters[iTime]) continue;
 468         fX[iTime] = fClusters[iTime]->GetX() - fX0;
 469       }
 470
 471       AliTRDseed::Update();
 472     }
 473     if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker)>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fPlane, fN2));
 474
 475     if(IsOK()){
 476       tquality = GetQuality(kZcorr);
 477       if(tquality < quality) break;
 478       else quality = tquality;
 479     }
 480     kroadz *= 2.;
 481   } // Loop: iter
 482   if (!IsOK()) return kFALSE;
 483
 484   CookLabels();
 485   UpdateUsed();
 486   return kTRUE;
 487 }
 488
 489 //____________________________________________________________________
 490 Bool_t  AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber
 491                                       ,Bool_t kZcorr)
 492 {
 493   //
 494   // Projective algorithm to attach clusters to seeding tracklets
 495   //
 496   // Parameters
 497   //
 498   // Output
 499   //
 500   // Detailed description
 501   // 1. Collapse x coordinate for the full detector plane
 502   // 2. truncated mean on y (r-phi) direction
 503   // 3. purge clusters
 504   // 4. truncated mean on z direction
 505   // 5. purge clusters
 506   // 6. fit tracklet
 507   //
 508
 509   if(!fReconstructor->GetRecoParam() ){
 510     AliError("Seed can not be used without a valid RecoParam.");
 511     return kFALSE;
 512   }
 513
 514   const Int_t kClusterCandidates = 2 * knTimebins;
 515
 516   //define roads
 517   Double_t kroady = fReconstructor->GetRecoParam() ->GetRoad1y();
 518   Double_t kroadz = fPadLength * 1.5 + 1.;
 519   // correction to y for the tilting angle
 520   Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
 521
 522   // working variables
 523   AliTRDcluster *clusters[kClusterCandidates];
 524   Double_t cond[4], yexp[knTimebins], zexp[knTimebins],
 525     yres[kClusterCandidates], zres[kClusterCandidates];
 526   Int_t ncl, *index = 0x0, tboundary[knTimebins];
 527
 528   // Do cluster projection
 529   AliTRDchamberTimeBin *layer = 0x0;
 530   Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
 531   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 532     if(!(layer = chamber->GetTB(iTime))) continue;
 533     if(!Int_t(*layer)) continue;
 534
 535     fX[iTime] = layer->GetX() - fX0;
 536     zexp[iTime] = fZref[0] + fZref[1] * fX[iTime];
 537     yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr;
 538
 539     // build condition and process clusters
 540     cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady;
 541     cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz;
 542     layer->GetClusters(cond, index, ncl);
 543     for(Int_t ic = 0; ic<ncl; ic++){
 544       AliTRDcluster *c = layer->GetCluster(index[ic]);
 545       clusters[nYclusters] = c;
 546       yres[nYclusters++] = c->GetY() - yexp[iTime];
 547       if(nYclusters >= kClusterCandidates) {
 548         AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates));
 549         kEXIT = kTRUE;
 550         break;
 551       }
 552     }
 553     tboundary[iTime] = nYclusters;
 554     if(kEXIT) break;
 555   }
 556
 557   // Evaluate truncated mean on the y direction
 558   Double_t mean, sigma;
 559   AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2);
 560   // purge cluster candidates
 561   Int_t nZclusters = 0;
 562   for(Int_t ic = 0; ic<nYclusters; ic++){
 563     if(yres[ic] - mean > 4. * sigma){
 564       clusters[ic] = 0x0;
 565       continue;
 566     }
 567     zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()];
 568   }
 569
 570   // Evaluate truncated mean on the z direction
 571   AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2);
 572   // purge cluster candidates
 573   for(Int_t ic = 0; ic<nZclusters; ic++){
 574     if(zres[ic] - mean > 4. * sigma){
 575       clusters[ic] = 0x0;
 576       continue;
 577     }
 578   }
 579
 580
 581   // Select only one cluster/TimeBin
 582   Int_t lastCluster = 0;
 583   fN2 = 0;
 584   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 585     ncl = tboundary[iTime] - lastCluster;
 586     if(!ncl) continue;
 587     Int_t iptr = lastCluster;
 588     if(ncl > 1){
 589       Float_t dold = 9999.;
 590       for(int ic=lastCluster; ic<tboundary[iTime]; ic++){
 591         if(!clusters[ic]) continue;
 592         Float_t y = yexp[iTime] - clusters[ic]->GetY();
 593         Float_t z = zexp[iTime] - clusters[ic]->GetZ();
 594         Float_t d = y * y + z * z;
 595         if(d > dold) continue;
 596         dold = d;
 597         iptr = ic;
 598       }
 599     }
 600     fIndexes[iTime]  = chamber->GetTB(iTime)->GetGlobalIndex(iptr);
 601     fClusters[iTime] = clusters[iptr];
 602     fY[iTime]        = clusters[iptr]->GetY();
 603     fZ[iTime]        = clusters[iptr]->GetZ();
 604     lastCluster      = tboundary[iTime];
 605     fN2++;
 606   }
 607
 608   // number of minimum numbers of clusters expected for the tracklet
 609   Int_t kClmin = Int_t(fReconstructor->GetRecoParam() ->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
 610   if (fN2 < kClmin){
 611     AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
 612     fN2 = 0;
 613     return kFALSE;
 614   }
 615
 616   // update used clusters
 617   fNUsed = 0;
 618   for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
 619     if(!fClusters[iTime]) continue;
 620     if((fClusters[iTime]->IsUsed())) fNUsed++;
 621   }
 622
 623   if (fN2-fNUsed < kClmin){
 624     AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
 625     fN2 = 0;
 626     return kFALSE;
 627   }
 628
 629   return kTRUE;
 630 }
 631
 632 //____________________________________________________________________
 633 Bool_t AliTRDseedV1::Fit()
 634 {
 635   //
 636   // Linear fit of the tracklet
 637   //
 638   // Parameters :
 639   //
 640   // Output :
 641   //  True if successful
 642   //
 643   // Detailed description
 644   // 2. Check if tracklet crosses pad row boundary
 645   // 1. Calculate residuals in the y (r-phi) direction
 646   // 3. Do a Least Square Fit to the data
 647   //
 648
 649   const Int_t kClmin = 8;
 650   const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
 651   AliTRDtrackerV1::AliTRDLeastSquare fitterY, fitterZ;
 652
 653   // convertion factor from square to gauss distribution for sigma
 654   Double_t convert = 1./TMath::Sqrt(12.);
 655
 656   // book cluster information
 657   Double_t xc[knTimebins+1], yc[knTimebins], zc[knTimebins+1], sy[knTimebins], sz[knTimebins+1];
 658   Int_t zRow[knTimebins];
 659   AliTRDcluster *c = 0x0;
 660   Int_t nc = 0;
 661   for (Int_t ic=0; ic<kNtb; ic++) {
 662     zRow[ic] = -1;
 663     xc[ic]  = -1.;
 664     yc[ic]  = 999.;
 665     zc[ic]  = 999.;
 666     sy[ic]  = 0.;
 667     sz[ic]  = 0.;
 668     if(!(c = fClusters[ic])) continue;
 669     if(!c->IsInChamber()) continue;
 670     Float_t w = 1.;
 671     if(c->GetNPads()>4) w = .5;
 672     if(c->GetNPads()>5) w = .2;
 673     zRow[nc] = c->GetPadRow();
 674     xc[nc]   = fX0 - c->GetX();
 675     yc[nc]   = c->GetY();
 676     zc[nc]   = c->GetZ();
 677     sy[nc]   = w; // all clusters have the same sigma
 678     sz[nc]   = fPadLength*convert;
 679     fitterZ.AddPoint(&xc[nc], zc[nc], sz[nc]);
 680     nc++;
 681   }
 682   // to few clusters
 683   if (nc < kClmin) return kFALSE;
 684
 685
 686   Int_t zN[2*35];
 687   Int_t nz = AliTRDtrackerV1::Freq(nc, zRow, zN, kFALSE);
 688   // more than one pad row crossing
 689   if(nz>2) return kFALSE;
 690
 691   // estimate reference parameter at average x
 692   Double_t y0 = fYref[0];
 693   Double_t dydx = fYref[1];
 694   Double_t dzdx = fZref[1];
 695   zc[nc]  = fZref[0];
 696
 697   // determine z offset of the fit
 698   Int_t nchanges = 0, nCross = 0;
 699   if(nz==2){ // tracklet is crossing pad row
 700     // Find the break time allowing one chage on pad-rows
 701     // with maximal number of accepted clusters
 702     Int_t padRef = zRow[0];
 703     for (Int_t ic=1; ic<nc; ic++) {
 704       if(zRow[ic] == padRef) continue;
 705
 706       // debug
 707       if(zRow[ic-1] == zRow[ic]){
 708         printf("ERROR in pad row change!!!\n");
 709       }
 710
 711       // evaluate parameters of the crossing point
 712       Float_t sx = (xc[ic-1] - xc[ic])*convert;
 713       xc[nc] = .5 * (xc[ic-1] + xc[ic]);
 714       zc[nc] = .5 * (zc[ic-1] + zc[ic]);
 715       sz[nc] = TMath::Max(dzdx * sx, .01);
 716       dzdx   = zc[ic-1] > zc[ic] ? 1. : -1.;
 717       padRef = zRow[ic];
 718       nCross = ic;
 719       nchanges++;
 720     }
 721   }
 722
 723   // condition on nCross and reset nchanges TODO
 724
 725   if(nchanges==1){
 726     if(dzdx * fZref[1] < 0.){
 727       AliInfo("tracklet direction does not correspond to the track direction. TODO.");
 728     }
 729     SetBit(kRowCross, kTRUE); // mark pad row crossing
 730     fCross[0] = xc[nc]; fCross[2] = zc[nc]; fCross[3] = sz[nc];
 731     fitterZ.AddPoint(&xc[nc], zc[nc], sz[nc]);
 732     fitterZ.Eval();
 733     dzdx = fZref[1]; // we don't trust Parameter[1] ??;
 734     zc[nc] = fitterZ.GetFunctionParameter(0);
 735   } else if(nchanges > 1){ // debug
 736     AliInfo("ERROR in n changes!!!");
 737     return kFALSE;
 738   }
 739
 740
 741   // estimate deviation from reference direction
 742   dzdx *= fTilt;
 743   for (Int_t ic=0; ic<nc; ic++) {
 744     yc[ic] -= y0 + xc[ic]*(dydx + dzdx) + fTilt * (zc[ic] - zc[nc]);
 745     fitterY.AddPoint(&xc[ic], yc[ic], sy[ic]);
 746   }
 747   fitterY.Eval();
 748   fYfit[0] = y0+fitterY.GetFunctionParameter(0);
 749   fYfit[1] = dydx+fitterY.GetFunctionParameter(1);
 750   if(nchanges) fCross[1] = fYfit[0] + fCross[0] * fYfit[1];
 751
 752 //      printf("\nnz = %d\n", nz);
 753 //      for(int ic=0; ic<35; ic++) printf("%d row[%d]\n", ic, zRow[ic]);
 754 //
 755 //      for(int ic=0; ic<nz; ic++) printf("%d n[%d]\n", ic, zN[ic]);
 756
 757   return kTRUE;
 758 }
 759
 760 //___________________________________________________________________
 761 void AliTRDseedV1::Draw(Option_t*)
 762 {
 763 }
 764
 765 //___________________________________________________________________
 766 void AliTRDseedV1::Print(Option_t*) const
 767 {
 768   //
 769   // Printing the seedstatus
 770   //
 771
 772   printf("Seed status :\n");
 773   printf("  fTilt      = %f\n", fTilt);
 774   printf("  fPadLength = %f\n", fPadLength);
 775   printf("  fX0        = %f\n", fX0);
 776   for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++) {
 777           const Char_t *isUsable = fUsable[ic]?"Yes":"No";
 778     printf("  %d X[%f] Y[%f] Z[%f] Indexes[%d] clusters[%p] usable[%s]\n"
 779                 , ic
 780                 , fX[ic]
 781                 , fY[ic]
 782                 , fZ[ic]
 783                 , fIndexes[ic]
 784                 , ((void*) fClusters[ic])
 785                 , isUsable);
 786         }
 787
 788   printf("  fYref[0] =%f fYref[1] =%f\n", fYref[0], fYref[1]);
 789   printf("  fZref[0] =%f fZref[1] =%f\n", fZref[0], fZref[1]);
 790   printf("  fYfit[0] =%f fYfit[1] =%f\n", fYfit[0], fYfit[1]);
 791   printf("  fYfitR[0]=%f fYfitR[1]=%f\n", fYfitR[0], fYfitR[1]);
 792   printf("  fZfit[0] =%f fZfit[1] =%f\n", fZfit[0], fZfit[1]);
 793   printf("  fZfitR[0]=%f fZfitR[1]=%f\n", fZfitR[0], fZfitR[1]);
 794   printf("  fSigmaY =%f\n", fSigmaY);
 795   printf("  fSigmaY2=%f\n", fSigmaY2);
 796   printf("  fMeanz  =%f\n", fMeanz);
 797   printf("  fZProb  =%f\n", fZProb);
 798   printf("  fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]);
 799   printf("  fN      =%d\n", fN);
 800   printf("  fN2     =%d (>8 isOK)\n",fN2);
 801   printf("  fNUsed  =%d\n", fNUsed);
 802   printf("  fFreq   =%d\n", fFreq);
 803   printf("  fNChange=%d\n",  fNChange);
 804   printf("  fMPads  =%f\n", fMPads);
 805
 806   printf("  fC      =%f\n", fC);
 807   printf("  fCC     =%f\n",fCC);
 808   printf("  fChi2   =%f\n", fChi2);
 809   printf("  fChi2Z  =%f\n", fChi2Z);
 810 }
 811