TRD/AliTRDtracker.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 standard TRD tracker                                                 //
  21 //  Based on Kalman filltering approach                                      //
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  25 #include <Riostream.h>
  26 #include <TFile.h>
  27 #include <TBranch.h>
  28 #include <TTree.h>
  29 #include <TObjArray.h>
  30
  31 #include "AliTRDgeometry.h"
  32 #include "AliTRDpadPlane.h"
  33 #include "AliTRDgeometryFull.h"
  34 #include "AliTRDcluster.h"
  35 #include "AliTRDtrack.h"
  36 #include "AliESD.h"
  37
  38 #include "AliTRDcalibDB.h"
  39 #include "AliTRDCommonParam.h"
  40
  41 #include "TTreeStream.h"
  42 #include "TGraph.h"
  43 #include "AliTRDtracker.h"
  44 #include "TLinearFitter.h"
  45 #include "AliRieman.h"
  46 #include "AliTrackPointArray.h"
  47 #include "AliAlignObj.h"
  48 #include "AliTRDReconstructor.h"
  49 //
  50
  51 ClassImp(AliTRDtracker)
  52 ClassImp(AliTRDseed)
  53
  54
  55
  56   const  Float_t     AliTRDtracker::fgkMinClustersInTrack = 0.5;
  57   const  Float_t     AliTRDtracker::fgkLabelFraction      = 0.8;
  58   const  Double_t    AliTRDtracker::fgkMaxChi2            = 12.;
  59   const    Double_t    AliTRDtracker::fgkMaxSnp           = 0.95;  // correspond to tan = 3
  60   const    Double_t    AliTRDtracker::fgkMaxStep          = 2.;     // maximal step size in propagation
  61
  62
  63 //
  64
  65
  66
  67
  68 //____________________________________________________________________
  69 AliTRDtracker::AliTRDtracker():AliTracker(),
  70                                fGeom(0),
  71                                fNclusters(0),
  72                                fClusters(0),
  73                                fNseeds(0),
  74                                fSeeds(0),
  75                                fNtracks(0),
  76                                fTracks(0),
  77                                fTimeBinsPerPlane(0),
  78                                fAddTRDseeds(kFALSE),
  79                                fNoTilt(kFALSE)
  80 {
  81   // Default constructor
  82
  83   for(Int_t i=0;i<kTrackingSectors;i++) fTrSec[i]=0;
  84   for(Int_t j=0;j<5;j++)
  85     for(Int_t k=0;k<18;k++) fHoles[j][k]=kFALSE;
  86   fDebugStreamer = 0;
  87 }
  88 //____________________________________________________________________
  89 AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
  90 {
  91   //
  92   //  Main constructor
  93   //
  94
  95   fAddTRDseeds = kFALSE;
  96   fGeom = NULL;
  97   fNoTilt = kFALSE;
  98
  99   TDirectory *savedir=gDirectory;
 100   TFile *in=(TFile*)geomfile;
 101   if (!in->IsOpen()) {
 102     printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n");
 103     printf("    FULL TRD geometry and DEFAULT TRD parameter will be used\n");
 104   }
 105   else {
 106     in->cd();
 107     fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
 108   }
 109
 110   if(fGeom) {
 111     //    printf("Found geometry version %d on file \n", fGeom->IsVersion());
 112   }
 113   else {
 114     printf("AliTRDtracker::AliTRDtracker(): can't find TRD geometry!\n");
 115     fGeom = new AliTRDgeometryFull();
 116     fGeom->SetPHOShole();
 117     fGeom->SetRICHhole();
 118   }
 119   fGeom->ReadGeoMatrices();
 120
 121   savedir->cd();
 122
 123
 124   fNclusters = 0;
 125   fClusters  = new TObjArray(2000);
 126   fNseeds    = 0;
 127   fSeeds     = new TObjArray(2000);
 128   fNtracks   = 0;
 129   fTracks    = new TObjArray(1000);
 130
 131   for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
 132     Int_t trS = CookSectorIndex(geomS);
 133     fTrSec[trS] = new AliTRDtrackingSector(fGeom, geomS);
 134     for (Int_t icham=0;icham<AliTRDgeometry::kNcham; icham++){
 135       fHoles[icham][trS]=fGeom->IsHole(0,icham,geomS);
 136     }
 137   }
 138   AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
 139   Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
 140   if(tiltAngle < 0.1) {
 141     fNoTilt = kTRUE;
 142   }
 143
 144   fTimeBinsPerPlane =  AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
 145
 146   fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
 147
 148   savedir->cd();
 149 }
 150
 151 //___________________________________________________________________
 152 AliTRDtracker::~AliTRDtracker()
 153 {
 154   //
 155   // Destructor of AliTRDtracker
 156   //
 157
 158   if (fClusters) {
 159     fClusters->Delete();
 160     delete fClusters;
 161   }
 162   if (fTracks) {
 163     fTracks->Delete();
 164     delete fTracks;
 165   }
 166   if (fSeeds) {
 167     fSeeds->Delete();
 168     delete fSeeds;
 169   }
 170   delete fGeom;
 171
 172   for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
 173     delete fTrSec[geomS];
 174   }
 175   if (fDebugStreamer) {
 176     //fDebugStreamer->Close();
 177     delete fDebugStreamer;
 178   }
 179 }
 180
 181 //_____________________________________________________________________
 182
 183
 184 Int_t  AliTRDtracker::LocalToGlobalID(Int_t lid){
 185   //
 186   // transform internal TRD ID to global detector ID
 187   //
 188   Int_t  isector = fGeom->GetSector(lid);
 189   Int_t  ichamber= fGeom->GetChamber(lid);
 190   Int_t  iplan   = fGeom->GetPlane(lid);
 191   //
 192   AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
 193   switch (iplan) {
 194   case 0:
 195     iLayer = AliAlignObj::kTRD1;
 196     break;
 197   case 1:
 198     iLayer = AliAlignObj::kTRD2;
 199     break;
 200   case 2:
 201     iLayer = AliAlignObj::kTRD3;
 202     break;
 203   case 3:
 204     iLayer = AliAlignObj::kTRD4;
 205     break;
 206   case 4:
 207     iLayer = AliAlignObj::kTRD5;
 208     break;
 209   case 5:
 210     iLayer = AliAlignObj::kTRD6;
 211     break;
 212   };
 213   Int_t modId = isector*fGeom->Ncham()+ichamber;
 214   UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
 215   return volid;
 216 }
 217
 218 Int_t  AliTRDtracker::GlobalToLocalID(Int_t gid){
 219   //
 220   // transform global detector ID to local detector ID
 221   //
 222   Int_t modId=0;
 223   AliAlignObj::ELayerID  layerId  = AliAlignObj::VolUIDToLayer(gid, modId);
 224   Int_t     isector  = modId/fGeom->Ncham();
 225   Int_t     ichamber = modId%fGeom->Ncham();
 226   Int_t     iLayer    = -1;
 227   switch (layerId) {
 228   case AliAlignObj::kTRD1:
 229     iLayer = 0;
 230     break;
 231   case AliAlignObj::kTRD2:
 232     iLayer = 1;
 233     break;
 234   case AliAlignObj::kTRD3:
 235     iLayer = 2;
 236     break;
 237   case AliAlignObj::kTRD4:
 238     iLayer = 3;
 239     break;
 240   case AliAlignObj::kTRD5:
 241     iLayer = 4;
 242     break;
 243   case AliAlignObj::kTRD6:
 244     iLayer = 5;
 245     break;
 246   default:
 247     iLayer =-1;
 248   }
 249   if (iLayer<0) return -1;
 250   Int_t lid = fGeom->GetDetector(iLayer,ichamber,isector);
 251   return lid;
 252 }
 253
 254
 255 Bool_t  AliTRDtracker::Transform(AliTRDcluster * cluster){
 256   //
 257   //
 258   //
 259   const Double_t kX0shift           = 2.52;    // magic constants for geo manager transformation
 260   const Double_t kX0shift5          = 3.05;    //
 261   //
 262   //
 263   // apply alignment and calibration to transform cluster
 264   //
 265   //
 266   Int_t detector = cluster->GetDetector();
 267   Int_t plane   = fGeom->GetPlane(cluster->GetDetector());
 268   Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
 269   Int_t sector  = fGeom->GetSector(cluster->GetDetector());
 270
 271   Double_t dxAmp  = (Double_t) fGeom->CamHght();          // Amplification region
 272   Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.);  // drift distance
 273   //
 274   // ExB correction
 275   //
 276   Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
 277   Double_t exB =   AliTRDcalibDB::Instance()->GetOmegaTau(vdrift);
 278   //
 279   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
 280   AliTRDpadPlane * padPlane = commonParam->GetPadPlane(plane,chamber);
 281   Double_t zshiftIdeal  = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
 282   Double_t localPos[3], localPosTracker[3];
 283   localPos[0] = -cluster->GetX();
 284   localPos[1] =  cluster->GetY() - driftX*exB;
 285   localPos[2] =  cluster->GetZ() -zshiftIdeal;
 286   //
 287   cluster->SetY(cluster->GetY() - driftX*exB);
 288   Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
 289   cluster->SetX(xplane- cluster->GetX());
 290   //
 291   TGeoHMatrix * matrix =  fGeom->GetCorrectionMatrix(cluster->GetDetector());
 292   if (!matrix){
 293     // no matrix found - if somebody used geometry with holes
 294     AliError("Invalid Geometry - Default Geometry used\n");
 295     return kTRUE;
 296   }
 297   matrix->LocalToMaster(localPos, localPosTracker);
 298   //
 299   //
 300   //
 301   if (AliTRDReconstructor::StreamLevel()>1){
 302     (*fDebugStreamer)<<"Transform"<<
 303       "Cl.="<<cluster<<
 304       "matrix.="<<matrix<<
 305       "Detector="<<detector<<
 306       "Sector="<<sector<<
 307       "Plane="<<plane<<
 308       "Chamber="<<chamber<<
 309       "lx0="<<localPosTracker[0]<<
 310       "ly0="<<localPosTracker[1]<<
 311       "lz0="<<localPosTracker[2]<<
 312       "\n";
 313   }
 314   //
 315   if (plane==5)
 316      cluster->SetX(localPosTracker[0]+kX0shift5);
 317   else
 318     cluster->SetX(localPosTracker[0]+kX0shift);
 319
 320   cluster->SetY(localPosTracker[1]);
 321   cluster->SetZ(localPosTracker[2]);
 322   return kTRUE;
 323 }
 324
 325 // Bool_t  AliTRDtracker::Transform(AliTRDcluster * cluster){
 326 //   //
 327 //   //
 328 //   const Double_t kDriftCorrection  = 1.01;                 // drift coeficient correction
 329 //   const Double_t kTime0Cor         = 0.32;                 // time0 correction
 330 //   //
 331 //   const Double_t kX0shift           = 2.52;
 332 //   const Double_t kX0shift5          = 3.05;
 333
 334 //   //
 335 //   // apply alignment and calibration to transform cluster
 336 //   //
 337 //   //
 338 //   Int_t detector = cluster->GetDetector();
 339 //   Int_t plane   = fGeom->GetPlane(cluster->GetDetector());
 340 //   Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
 341 //   Int_t sector  = fGeom->GetSector(cluster->GetDetector());
 342
 343 //   Double_t dxAmp  = (Double_t) fGeom->CamHght();          // Amplification region
 344 //   Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.);  // drift distance
 345 //   //
 346 //   // ExB correction
 347 //   //
 348 //   Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
 349 //   Double_t exB =   AliTRDcalibDB::Instance()->GetOmegaTau(vdrift);
 350 //   //
 351
 352 //   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
 353 //   AliTRDpadPlane * padPlane = commonParam->GetPadPlane(plane,chamber);
 354 //   Double_t zshiftIdeal  = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
 355 //   Double_t localPos[3], globalPos[3], localPosTracker[3], localPosTracker2[3];
 356 //   localPos[2] = -cluster->GetX();
 357 //   localPos[0] =  cluster->GetY() - driftX*exB;
 358 //   localPos[1] =  cluster->GetZ() -zshiftIdeal;
 359 //   TGeoHMatrix * matrix =  fGeom->GetGeoMatrix(cluster->GetDetector());
 360 //   matrix->LocalToMaster(localPos, globalPos);
 361
 362 //   Double_t sectorAngle = 20.*(sector%18)+10;
 363 //   TGeoHMatrix  rotSector;
 364 //   rotSector.RotateZ(sectorAngle);
 365 //   rotSector.LocalToMaster(globalPos, localPosTracker);
 366 //   //
 367 //   //
 368 //   TGeoHMatrix  matrix2(*matrix);
 369 //   matrix2.MultiplyLeft(&rotSector);
 370 //   matrix2.LocalToMaster(localPos,localPosTracker2);
 371 //   //
 372 //   //
 373 //   //
 374 //   cluster->SetY(cluster->GetY() - driftX*exB);
 375 //   Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
 376 //   cluster->SetX(xplane- kDriftCorrection*(cluster->GetX()-kTime0Cor));
 377 //   (*fDebugStreamer)<<"Transform"<<
 378 //     "Cl.="<<cluster<<
 379 //     "matrix.="<<matrix<<
 380 //     "matrix2.="<<&matrix2<<
 381 //     "Detector="<<detector<<
 382 //     "Sector="<<sector<<
 383 //     "Plane="<<plane<<
 384 //     "Chamber="<<chamber<<
 385 //     "lx0="<<localPosTracker[0]<<
 386 //     "ly0="<<localPosTracker[1]<<
 387 //     "lz0="<<localPosTracker[2]<<
 388 //     "lx2="<<localPosTracker2[0]<<
 389 //     "ly2="<<localPosTracker2[1]<<
 390 //     "lz2="<<localPosTracker2[2]<<
 391 //     "\n";
 392 //   //
 393 //   if (plane==5)
 394 //      cluster->SetX(localPosTracker[0]+kX0shift5);
 395 //   else
 396 //     cluster->SetX(localPosTracker[0]+kX0shift);
 397
 398 //   cluster->SetY(localPosTracker[1]);
 399 //   cluster->SetZ(localPosTracker[2]);
 400 //   return kTRUE;
 401 // }
 402
 403 Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) {
 404   //
 405   // Rotates the track when necessary
 406   //
 407
 408   Double_t alpha = AliTRDgeometry::GetAlpha();
 409   Double_t y = track->GetY();
 410   Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
 411
 412   //Int_t ns = AliTRDgeometry::kNsect;
 413   //Int_t s=Int_t(track->GetAlpha()/alpha)%ns;
 414
 415   if (y > ymax) {
 416     //s = (s+1) % ns;
 417     if (!track->Rotate(alpha)) return kFALSE;
 418   } else if (y <-ymax) {
 419     //s = (s-1+ns) % ns;
 420     if (!track->Rotate(-alpha)) return kFALSE;
 421   }
 422
 423   return kTRUE;
 424 }
 425
 426
 427 AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin, UInt_t &index){
 428   //
 429   //try to find cluster in the backup list
 430   //
 431   AliTRDcluster * cl =0;
 432   Int_t *indexes = track->GetBackupIndexes();
 433   for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
 434     if (indexes[i]==0) break;
 435     AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
 436     if (!cli) break;
 437     if (cli->GetLocalTimeBin()!=timebin) continue;
 438     Int_t iplane = fGeom->GetPlane(cli->GetDetector());
 439     if (iplane==plane) {
 440       cl = cli;
 441       index = indexes[i];
 442       break;
 443     }
 444   }
 445   return cl;
 446 }
 447
 448
 449 Int_t  AliTRDtracker::GetLastPlane(AliTRDtrack * track){
 450   //
 451   //return last updated plane
 452   Int_t lastplane=0;
 453   Int_t *indexes = track->GetBackupIndexes();
 454   for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
 455     AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
 456     if (!cli) break;
 457     Int_t iplane = fGeom->GetPlane(cli->GetDetector());
 458     if (iplane>lastplane) {
 459       lastplane = iplane;
 460     }
 461   }
 462   return lastplane;
 463 }
 464 //___________________________________________________________________
 465 Int_t AliTRDtracker::Clusters2Tracks(AliESD* event)
 466 {
 467   //
 468   // Finds tracks within the TRD. The ESD event is expected to contain seeds
 469   // at the outer part of the TRD. The seeds
 470   // are found within the TRD if fAddTRDseeds is TRUE.
 471   // The tracks are propagated to the innermost time bin
 472   // of the TRD and the ESD event is updated
 473   //
 474
 475   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 476   Float_t foundMin = fgkMinClustersInTrack * timeBins;
 477   Int_t nseed = 0;
 478   Int_t found = 0;
 479   //  Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
 480
 481   Int_t n = event->GetNumberOfTracks();
 482   for (Int_t i=0; i<n; i++) {
 483     AliESDtrack* seed=event->GetTrack(i);
 484     ULong_t status=seed->GetStatus();
 485     if ( (status & AliESDtrack::kTRDout ) == 0 ) continue;
 486     if ( (status & AliESDtrack::kTRDin) != 0 ) continue;
 487     nseed++;
 488
 489     AliTRDtrack* seed2 = new AliTRDtrack(*seed);
 490     //seed2->ResetCovariance();
 491     AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
 492     AliTRDtrack &t=*pt;
 493     FollowProlongation(t);
 494     if (t.GetNumberOfClusters() >= foundMin) {
 495       UseClusters(&t);
 496       CookLabel(pt, 1-fgkLabelFraction);
 497       //      t.CookdEdx();
 498     }
 499     found++;
 500 //    cout<<found<<'\r';
 501
 502     Double_t xTPC = 250;
 503     if (PropagateToX(t,xTPC,fgkMaxStep)) {
 504       seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
 505     }
 506     delete seed2;
 507     delete pt;
 508   }
 509
 510   cout<<"Number of loaded seeds: "<<nseed<<endl;
 511   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
 512
 513   // after tracks from loaded seeds are found and the corresponding
 514   // clusters are used, look for additional seeds from TRD
 515
 516
 517   cout<<"Total number of found tracks: "<<found<<endl;
 518
 519   return 0;
 520 }
 521
 522
 523
 524 //_____________________________________________________________________________
 525 Int_t AliTRDtracker::PropagateBack(AliESD* event) {
 526   //
 527   // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
 528   // backpropagated by the TPC tracker. Each seed is first propagated
 529   // to the TRD, and then its prolongation is searched in the TRD.
 530   // If sufficiently long continuation of the track is found in the TRD
 531   // the track is updated, otherwise it's stored as originaly defined
 532   // by the TPC tracker.
 533   //
 534
 535   Int_t found=0;
 536   Float_t foundMin = 20;
 537   Int_t n = event->GetNumberOfTracks();
 538   //
 539   //Sort tracks
 540   Float_t *quality =new Float_t[n];
 541   Int_t *index   =new Int_t[n];
 542   for (Int_t i=0; i<n; i++) {
 543     AliESDtrack* seed=event->GetTrack(i);
 544     Double_t covariance[15];
 545     seed->GetExternalCovariance(covariance);
 546     quality[i] = covariance[0]+covariance[2];
 547   }
 548   TMath::Sort(n,quality,index,kFALSE);
 549   //
 550   for (Int_t i=0; i<n; i++) {
 551     //    AliESDtrack* seed=event->GetTrack(i);
 552     AliESDtrack* seed=event->GetTrack(index[i]);
 553
 554     ULong_t status=seed->GetStatus();
 555     if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
 556     if ( (status & AliESDtrack::kTRDout) != 0 ) continue;
 557
 558     Int_t lbl = seed->GetLabel();
 559     AliTRDtrack *track = new AliTRDtrack(*seed);
 560     track->SetSeedLabel(lbl);
 561     seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); //make backup
 562     fNseeds++;
 563     Float_t p4     = track->GetC();
 564     //
 565     Int_t expectedClr = FollowBackProlongation(*track);
 566     if (TMath::Abs(track->GetC()-p4)/TMath::Abs(p4)<0.2 || TMath::Abs(track->GetPt())>0.8 ) {
 567       //
 568       //make backup for back propagation
 569       //
 570       Int_t foundClr = track->GetNumberOfClusters();
 571       if (foundClr >= foundMin) {
 572         track->CookdEdx();
 573         CookdEdxTimBin(*track);
 574         CookLabel(track, 1-fgkLabelFraction);
 575         if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
 576         if(track->GetChi2()/track->GetNumberOfClusters()<4) {   // sign only gold tracks
 577           if (seed->GetKinkIndex(0)==0&&TMath::Abs(track->GetPt())<1.5 ) UseClusters(track);
 578         }
 579         Bool_t isGold = kFALSE;
 580
 581         if (track->GetChi2()/track->GetNumberOfClusters()<5) {  //full gold track
 582           // seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
 583            if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
 584           isGold = kTRUE;
 585         }
 586         if (!isGold && track->GetNCross()==0&&track->GetChi2()/track->GetNumberOfClusters()<7){ //almost gold track
 587           //      seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
 588           if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
 589           isGold = kTRUE;
 590         }
 591         if (!isGold && track->GetBackupTrack()){
 592           if (track->GetBackupTrack()->GetNumberOfClusters()>foundMin&&
 593               (track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1))<7){
 594             seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
 595             isGold = kTRUE;
 596           }
 597         }
 598         if (track->StatusForTOF()>0 &&track->fNCross==0 && Float_t(track->fN)/Float_t(track->fNExpected)>0.4){
 599           //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
 600         }
 601       }
 602     }
 603     // Debug part of tracking
 604     TTreeSRedirector& cstream = *fDebugStreamer;
 605     Int_t eventNr = event->GetEventNumber();
 606     if (AliTRDReconstructor::StreamLevel()>0){
 607       if (track->GetBackupTrack()){
 608         cstream<<"Tracks"<<
 609           "EventNr="<<eventNr<<
 610           "ESD.="<<seed<<
 611           "trd.="<<track<<
 612           "trdback.="<<track->GetBackupTrack()<<
 613           "\n";
 614       }else{
 615         cstream<<"Tracks"<<
 616           "EventNr="<<eventNr<<
 617           "ESD.="<<seed<<
 618           "trd.="<<track<<
 619           "trdback.="<<track<<
 620           "\n";
 621       }
 622     }
 623     //
 624     //Propagation to the TOF (I.Belikov)
 625     if (track->GetStop()==kFALSE){
 626
 627       Double_t xtof=371.;
 628       Double_t c2=track->GetC()*xtof - track->GetEta();
 629       if (TMath::Abs(c2)>=0.99) {
 630         delete track;
 631         continue;
 632       }
 633       Double_t xTOF0 = 370. ;
 634       PropagateToX(*track,xTOF0,fgkMaxStep);
 635       //
 636       //energy losses taken to the account - check one more time
 637       c2=track->GetC()*xtof - track->GetEta();
 638       if (TMath::Abs(c2)>=0.99) {
 639         delete track;
 640         continue;
 641       }
 642
 643       //
 644       Double_t ymax=xtof*TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
 645       Double_t y=track->GetYat(xtof);
 646       if (y > ymax) {
 647         if (!track->Rotate(AliTRDgeometry::GetAlpha())) {
 648           delete track;
 649           continue;
 650         }
 651       } else if (y <-ymax) {
 652         if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
 653           delete track;
 654           continue;
 655         }
 656       }
 657
 658       if (track->PropagateTo(xtof)) {
 659         seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
 660         for (Int_t i=0;i<kNPlane;i++) {
 661            seed->SetTRDsignals(track->GetPIDsignals(i),i);
 662            seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
 663         }
 664         //      seed->SetTRDtrack(new AliTRDtrack(*track));
 665         if (track->GetNumberOfClusters()>foundMin) found++;
 666       }
 667     }else{
 668       if (track->GetNumberOfClusters()>15&&track->GetNumberOfClusters()>0.5*expectedClr){
 669         seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
 670         //seed->SetStatus(AliESDtrack::kTRDStop);
 671         for (Int_t i=0;i<kNPlane;i++) {
 672            seed->SetTRDsignals(track->GetPIDsignals(i),i);
 673            seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
 674         }
 675         //seed->SetTRDtrack(new AliTRDtrack(*track));
 676         found++;
 677       }
 678     }
 679     seed->SetTRDQuality(track->StatusForTOF());
 680     seed->SetTRDBudget(track->fBudget[0]);
 681
 682     delete track;
 683     //
 684     //End of propagation to the TOF
 685     //if (foundClr>foundMin)
 686     //  seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
 687
 688
 689   }
 690
 691   cerr<<"Number of seeds: "<<fNseeds<<endl;
 692   cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
 693
 694   if (AliTRDReconstructor::SeedingOn()) MakeSeedsMI(3,5,event); //new seeding
 695
 696   fSeeds->Clear(); fNseeds=0;
 697   delete [] index;
 698   delete [] quality;
 699
 700   return 0;
 701
 702 }
 703
 704 //_____________________________________________________________________________
 705 Int_t AliTRDtracker::RefitInward(AliESD* event)
 706 {
 707   //
 708   // Refits tracks within the TRD. The ESD event is expected to contain seeds
 709   // at the outer part of the TRD.
 710   // The tracks are propagated to the innermost time bin
 711   // of the TRD and the ESD event is updated
 712   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
 713   //
 714
 715   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 716   Float_t foundMin = fgkMinClustersInTrack * timeBins;
 717   Int_t nseed = 0;
 718   Int_t found = 0;
 719   //  Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
 720   AliTRDtrack seed2;
 721
 722   Int_t n = event->GetNumberOfTracks();
 723   for (Int_t i=0; i<n; i++) {
 724     AliESDtrack* seed=event->GetTrack(i);
 725     new(&seed2) AliTRDtrack(*seed);
 726     if (seed2.GetX()<270){
 727       seed->UpdateTrackParams(&seed2, AliESDtrack::kTRDbackup); // backup TPC track - only update
 728       continue;
 729     }
 730
 731     ULong_t status=seed->GetStatus();
 732     if ( (status & AliESDtrack::kTRDout ) == 0 ) {
 733       continue;
 734     }
 735     if ( (status & AliESDtrack::kTRDin) != 0 ) {
 736       continue;
 737     }
 738     nseed++;
 739 //     if (1/seed2.Get1Pt()>1.5&& seed2.GetX()>260.) {
 740 //       Double_t oldx = seed2.GetX();
 741 //       seed2.PropagateTo(500.);
 742 //       seed2.ResetCovariance(1.);
 743 //       seed2.PropagateTo(oldx);
 744 //     }
 745 //     else{
 746 //       seed2.ResetCovariance(5.);
 747 //     }
 748
 749     AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
 750     Int_t * indexes2 = seed2.GetIndexes();
 751     for (Int_t i=0;i<kNPlane;i++) {
 752       pt->SetPIDsignals(seed2.GetPIDsignals(i),i);
 753       pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i);
 754     }
 755
 756     Int_t * indexes3 = pt->GetBackupIndexes();
 757     for (Int_t i=0;i<200;i++) {
 758       if (indexes2[i]==0) break;
 759       indexes3[i] = indexes2[i];
 760     }
 761     //AliTRDtrack *pt = seed2;
 762     AliTRDtrack &t=*pt;
 763     FollowProlongation(t);
 764     if (t.GetNumberOfClusters() >= foundMin) {
 765       //      UseClusters(&t);
 766       //CookLabel(pt, 1-fgkLabelFraction);
 767       t.CookdEdx();
 768       CookdEdxTimBin(t);
 769     }
 770     found++;
 771 //    cout<<found<<'\r';
 772     Double_t xTPC = 250;
 773     if(PropagateToX(t,xTPC,fgkMaxStep)) {
 774       seed->UpdateTrackParams(pt, AliESDtrack::kTRDrefit);
 775       for (Int_t i=0;i<kNPlane;i++) {
 776         seed->SetTRDsignals(pt->GetPIDsignals(i),i);
 777         seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
 778       }
 779     }else{
 780       //if not prolongation to TPC - propagate without update
 781       AliTRDtrack* seed2 = new AliTRDtrack(*seed);
 782       seed2->ResetCovariance(5.);
 783       AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha());
 784       delete seed2;
 785       if (PropagateToX(*pt2,xTPC,fgkMaxStep)) {
 786         //pt2->CookdEdx(0.,1.);
 787         pt2->CookdEdx( ); // Modification by PS
 788         CookdEdxTimBin(*pt2);
 789         seed->UpdateTrackParams(pt2, AliESDtrack::kTRDrefit);
 790         for (Int_t i=0;i<kNPlane;i++) {
 791           seed->SetTRDsignals(pt2->GetPIDsignals(i),i);
 792           seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
 793         }
 794       }
 795       delete pt2;
 796     }
 797     delete pt;
 798   }
 799
 800   cout<<"Number of loaded seeds: "<<nseed<<endl;
 801   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
 802
 803   return 0;
 804
 805 }
 806
 807
 808
 809
 810
 811 //---------------------------------------------------------------------------
 812 Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t)
 813 {
 814   // Starting from current position on track=t this function tries
 815   // to extrapolate the track up to timeBin=0 and to confirm prolongation
 816   // if a close cluster is found. Returns the number of clusters
 817   // expected to be found in sensitive layers
 818   // GeoManager used to estimate mean density
 819   Int_t sector;
 820   Int_t lastplane = GetLastPlane(&t);
 821   Double_t radLength = 0.0;
 822   Double_t rho = 0.0;
 823   Int_t expectedNumberOfClusters = 0;
 824   //
 825   //
 826   //
 827   for (Int_t iplane = lastplane; iplane>=0; iplane--){
 828     //
 829     Int_t    row0 = GetGlobalTimeBin(0, iplane,GetTimeBinsPerPlane()-1);
 830     Int_t    rowlast = GetGlobalTimeBin(0, iplane,0);
 831     //
 832     // propagate track close to the plane if neccessary
 833     //
 834     Double_t currentx  = fTrSec[0]->GetLayer(rowlast)->GetX();
 835     if (currentx < -fgkMaxStep +t.GetX()){
 836       //propagate closer to chamber - safety space fgkMaxStep
 837       if (!PropagateToX(t, currentx+fgkMaxStep, fgkMaxStep)) break;
 838     }
 839     if (!AdjustSector(&t)) break;
 840     //
 841     // get material budget
 842     //
 843     Double_t xyz0[3],xyz1[3],param[7],x,y,z;
 844     t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);   //starting global position
 845     // end global position
 846     x = fTrSec[0]->GetLayer(row0)->GetX();
 847     if (!t.GetProlongation(x,y,z)) break;
 848     xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
 849     xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
 850     xyz1[2] = z;
 851     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 852     rho = param[0];
 853     radLength = param[1];   // get mean propagation parameters
 854     //
 855     // propagate nad update
 856     //
 857     sector = t.GetSector();
 858     //    for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
 859     for (Int_t itime=0 ;itime<GetTimeBinsPerPlane();itime++) {
 860       Int_t    ilayer = GetGlobalTimeBin(0, iplane,itime);
 861       expectedNumberOfClusters++;
 862       t.fNExpected++;
 863       if (t.fX>345) t.fNExpectedLast++;
 864       AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(ilayer));
 865       AliTRDcluster *cl=0;
 866       UInt_t index=0;
 867       Double_t maxChi2=fgkMaxChi2;
 868       x = timeBin.GetX();
 869       if (timeBin) {
 870         AliTRDcluster * cl0 = timeBin[0];
 871         if (!cl0) continue;         // no clusters in given time bin
 872         Int_t plane = fGeom->GetPlane(cl0->GetDetector());
 873         if (plane>lastplane) continue;
 874         Int_t timebin = cl0->GetLocalTimeBin();
 875         AliTRDcluster * cl2= GetCluster(&t,plane, timebin,index);
 876         //
 877         if (cl2) {
 878           cl =cl2;
 879           Double_t h01 = GetTiltFactor(cl);
 880           maxChi2=t.GetPredictedChi2(cl,h01);
 881         }
 882         if (cl) {
 883           //      if (cl->GetNPads()<5)
 884           Double_t dxsample = timeBin.GetdX();
 885           t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
 886           Double_t h01 = GetTiltFactor(cl);
 887           Int_t det = cl->GetDetector();
 888           Int_t plane = fGeom->GetPlane(det);
 889           if (t.fX>345){
 890             t.fNLast++;
 891             t.fChi2Last+=maxChi2;
 892           }
 893           Double_t xcluster = cl->GetX();
 894           t.PropagateTo(xcluster,radLength,rho);
 895           if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
 896           }
 897         }
 898       }
 899     }
 900   }
 901   return expectedNumberOfClusters;
 902 }
 903
 904
 905
 906
 907
 908 //___________________________________________________________________
 909 Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
 910 {
 911
 912   // Starting from current radial position of track <t> this function
 913   // extrapolates the track up to outer timebin and in the sensitive
 914   // layers confirms prolongation if a close cluster is found.
 915   // Returns the number of clusters expected to be found in sensitive layers
 916   // Use GEO manager for material Description
 917
 918   Int_t sector;
 919   Int_t clusters[1000];
 920   for (Int_t i=0;i<1000;i++) clusters[i]=-1;
 921   Double_t radLength = 0.0;
 922   Double_t rho = 0.0;
 923   Int_t expectedNumberOfClusters = 0;
 924   Float_t ratio0=0;
 925   AliTRDtracklet tracklet;
 926   //
 927   //
 928   for (Int_t iplane = 0; iplane<kNPlane; iplane++){
 929     Int_t    row0    = GetGlobalTimeBin(0, iplane,GetTimeBinsPerPlane()-1);
 930     Int_t    rowlast = GetGlobalTimeBin(0, iplane,0);
 931     //
 932     Double_t currentx  = fTrSec[0]->GetLayer(row0)->GetX();
 933     if (currentx<t.GetX()) continue;
 934     //
 935     //       propagate closer to chamber if neccessary
 936     //
 937     if (currentx > fgkMaxStep +t.GetX()){
 938       if (!PropagateToX(t, currentx-fgkMaxStep, fgkMaxStep)) break;
 939     }
 940     if (!AdjustSector(&t)) break;
 941     if (TMath::Abs(t.GetSnp())>fgkMaxSnp) break;
 942     //
 943     // get material budget inside of chamber
 944     //
 945     Double_t xyz0[3],xyz1[3],param[7],x,y,z;
 946     t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);   //starting global position
 947     // end global position
 948     x = fTrSec[0]->GetLayer(rowlast)->GetX();
 949     if (!t.GetProlongation(x,y,z)) break;
 950     xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
 951     xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
 952     xyz1[2] = z;
 953     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 954     rho = param[0];
 955     radLength = param[1];   // get mean propagation parameters
 956     //
 957     // Find clusters
 958     //
 959     sector = t.GetSector();
 960     Float_t  ncl   = FindClusters(sector,row0,rowlast,&t,clusters,tracklet);
 961     if (tracklet.GetN()<GetTimeBinsPerPlane()/3) continue;
 962     //
 963     // Propagate and update track
 964     //
 965     for (Int_t itime= GetTimeBinsPerPlane()-1;itime>=0;itime--) {
 966       Int_t    ilayer = GetGlobalTimeBin(0, iplane,itime);
 967       expectedNumberOfClusters++;
 968       t.fNExpected++;
 969       if (t.fX>345) t.fNExpectedLast++;
 970       AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(ilayer));
 971       AliTRDcluster *cl=0;
 972       UInt_t index=0;
 973       Double_t maxChi2=fgkMaxChi2;
 974       x = timeBin.GetX();
 975       //
 976       if (timeBin) {
 977         if (clusters[ilayer]>0) {
 978           index = clusters[ilayer];
 979           cl    = (AliTRDcluster*)GetCluster(index);
 980           Double_t h01 = GetTiltFactor(cl);
 981           maxChi2=t.GetPredictedChi2(cl,h01);
 982         }
 983
 984         if (cl) {
 985           //      if (cl->GetNPads()<5)
 986           Double_t dxsample = timeBin.GetdX();
 987           t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
 988           Double_t h01 = GetTiltFactor(cl);
 989           Int_t det = cl->GetDetector();
 990           Int_t plane = fGeom->GetPlane(det);
 991           if (t.fX>345){
 992             t.fNLast++;
 993             t.fChi2Last+=maxChi2;
 994           }
 995           Double_t xcluster = cl->GetX();
 996           t.PropagateTo(xcluster,radLength,rho);
 997           if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
 998             if(!t.Update(cl,maxChi2,index,h01)) {
 999             }
1000           }
1001           //
1002           // reset material budget if 2 consecutive gold
1003           if (plane>0)
1004             if (t.fTracklets[plane].GetN()+t.fTracklets[plane-1].GetN()>20){
1005               t.fBudget[2] = 0;
1006             }
1007         }
1008       }
1009     }
1010     ratio0 = ncl/Float_t(fTimeBinsPerPlane);
1011     Float_t  ratio1 = Float_t(t.fN+1)/Float_t(t.fNExpected+1.);
1012     if (tracklet.GetChi2()<18.&&ratio0>0.8 && ratio1>0.6 && ratio0+ratio1>1.5 && t.GetNCross()==0 && TMath::Abs(t.GetSnp())<0.85&&t.fN>20){
1013       t.MakeBackupTrack();                            // make backup of the track until is gold
1014     }
1015
1016   }
1017   //
1018   return expectedNumberOfClusters;
1019 }
1020
1021
1022
1023 Int_t  AliTRDtracker::PropagateToX(AliTRDtrack& t, Double_t xToGo, Double_t maxStep)
1024 {
1025   // Starting from current radial position of track <t> this function
1026   // extrapolates the track up to radial position <xToGo>.
1027   // Returns 1 if track reaches the plane, and 0 otherwise
1028   const Double_t kEpsilon = 0.00001;
1029   //  Double_t tanmax = TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
1030   Double_t xpos     = t.GetX();
1031   Double_t dir      = (xpos<xToGo) ? 1.:-1.;
1032   //
1033   while ( (xToGo-xpos)*dir > kEpsilon){
1034     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
1035     //
1036     Double_t xyz0[3],xyz1[3],param[7],x,y,z;
1037     t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);   //starting global position
1038     x    = xpos+step;
1039     //
1040     if (!t.GetProlongation(x,y,z)) return 0;   // no prolongation
1041     //
1042     xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
1043     xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
1044     xyz1[2] = z;
1045     //
1046     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1047     if (!t.PropagateTo(x,param[1],param[0])) return 0;
1048     AdjustSector(&t);
1049     xpos = t.GetX();
1050   }
1051   return 1;
1052
1053 }
1054
1055
1056
1057 //_____________________________________________________________________________
1058 Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1059 {
1060   // Fills clusters into TRD tracking_sectors
1061   // Note that the numbering scheme for the TRD tracking_sectors
1062   // differs from that of TRD sectors
1063   cout<<"\n Read Sectors  clusters"<<endl;
1064   if (ReadClusters(fClusters,cTree)) {
1065      Error("LoadClusters","Problem with reading the clusters !");
1066      return 1;
1067   }
1068   Int_t ncl=fClusters->GetEntriesFast();
1069   fNclusters=ncl;
1070   cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1071
1072   UInt_t index;
1073   for (Int_t ichamber=0;ichamber<5;ichamber++)
1074     for (Int_t isector=0;isector<18;isector++){
1075       fHoles[ichamber][isector]=kTRUE;
1076     }
1077
1078
1079   while (ncl--) {
1080 //    printf("\r %d left  ",ncl);
1081     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
1082     Int_t detector=c->GetDetector();
1083     Int_t localTimeBin=c->GetLocalTimeBin();
1084     Int_t sector=fGeom->GetSector(detector);
1085     Int_t plane=fGeom->GetPlane(detector);
1086
1087     Int_t trackingSector = CookSectorIndex(sector);
1088     if (c->GetLabel(0)>0){
1089       Int_t chamber = fGeom->GetChamber(detector);
1090       fHoles[chamber][trackingSector]=kFALSE;
1091     }
1092
1093     Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
1094     if(gtb < 0) continue;
1095     Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1096
1097     index=ncl;
1098     //
1099     // apply pos correction
1100     Transform(c);
1101     fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
1102   }
1103   return 0;
1104 }
1105
1106 //_____________________________________________________________________________
1107 void AliTRDtracker::UnloadClusters()
1108 {
1109   //
1110   // Clears the arrays of clusters and tracks. Resets sectors and timebins
1111   //
1112
1113   Int_t i, nentr;
1114
1115   nentr = fClusters->GetEntriesFast();
1116   for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i);
1117   fNclusters = 0;
1118
1119   nentr = fSeeds->GetEntriesFast();
1120   for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i);
1121
1122   nentr = fTracks->GetEntriesFast();
1123   for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i);
1124
1125   Int_t nsec = AliTRDgeometry::kNsect;
1126
1127   for (i = 0; i < nsec; i++) {
1128     for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1129       fTrSec[i]->GetLayer(pl)->Clear();
1130     }
1131   }
1132
1133 }
1134
1135 //__________________________________________________________________________
1136 void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD * esd)
1137 {
1138   //
1139   // Creates  seeds using clusters between  position inner plane  and outer plane
1140   //
1141   const Double_t kMaxTheta = 1;
1142   const Double_t kMaxPhi   = 2.0;
1143   //
1144   const Double_t kRoad0y  =  6;     // road for middle cluster
1145   const Double_t kRoad0z  =  8.5;   // road for middle cluster
1146   //
1147   const Double_t kRoad1y  =  2;    // road in y for seeded cluster
1148   const Double_t kRoad1z  =  20;    // road in z for seeded cluster
1149   //
1150   const Double_t kRoad2y  =  3;    // road in y for extrapolated cluster
1151   const Double_t kRoad2z  =  20;   // road in z for extrapolated cluster
1152   const Int_t    kMaxSeed  = 3000;
1153   Int_t maxSec=AliTRDgeometry::kNsect;
1154
1155   //
1156   // linear fitters in planes
1157   TLinearFitter fitterTC(2,"hyp2");  // fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
1158   TLinearFitter fitterT2(4,"hyp4");  // fitting with tilting pads - kz not fixed
1159   fitterTC.StoreData(kTRUE);
1160   fitterT2.StoreData(kTRUE);
1161   AliRieman rieman(1000);   // rieman fitter
1162   AliRieman rieman2(1000);   // rieman fitter
1163   //
1164   // find the maximal and minimal layer for the planes
1165   //
1166   Int_t layers[6][2];
1167   AliTRDpropagationLayer* reflayers[6];
1168   for (Int_t i=0;i<6;i++){layers[i][0]=10000; layers[i][1]=0;}
1169   for (Int_t ns=0;ns<maxSec;ns++){
1170     for (Int_t ilayer=0;ilayer<fTrSec[ns]->GetNumberOfLayers();ilayer++){
1171       AliTRDpropagationLayer& layer=*(fTrSec[ns]->GetLayer(ilayer));
1172       if (layer==0) continue;
1173       Int_t det   = layer[0]->GetDetector();
1174       Int_t plane = fGeom->GetPlane(det);
1175       if (ilayer<layers[plane][0]) layers[plane][0] = ilayer;
1176       if (ilayer>layers[plane][1]) layers[plane][1] = ilayer;
1177     }
1178   }
1179   //
1180   AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
1181   Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
1182   Double_t hL[6];         // tilting angle
1183   Double_t xcl[6];        // x - position of reference cluster
1184   Double_t ycl[6];        // y - position of reference cluster
1185   Double_t zcl[6];        // z - position of reference cluster
1186   AliTRDcluster *cl[6]={0,0,0,0,0,0};    // seeding clusters
1187   Float_t padlength[6]={10,10,10,10,10,10};   //current pad-length
1188   Double_t chi2R =0, chi2Z=0;
1189   Double_t chi2RF =0, chi2ZF=0;
1190   //
1191   Int_t nclusters;     // total number of clusters
1192   for (Int_t i=0;i<6;i++) {hL[i]=h01; if (i%2==1) hL[i]*=-1.;}
1193   //
1194   //
1195   //         registered seed
1196   AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
1197   AliTRDseed *seed[kMaxSeed];
1198   for (Int_t iseed=0;iseed<kMaxSeed;iseed++) seed[iseed]= &pseed[iseed*6];
1199   AliTRDseed *cseed = seed[0];
1200   //
1201   Double_t   seedquality[kMaxSeed];
1202   Double_t   seedquality2[kMaxSeed];
1203   Double_t   seedparams[kMaxSeed][7];
1204   Int_t      seedlayer[kMaxSeed];
1205   Int_t      registered =0;
1206   Int_t      sort[kMaxSeed];
1207   //
1208   // seeding part
1209   //
1210   for (Int_t ns = 0; ns<maxSec; ns++){         //loop over sectors
1211   //for (Int_t ns = 0; ns<5; ns++){         //loop over sectors
1212     registered = 0;   // reset registerd seed counter
1213     cseed      = seed[registered];
1214     Float_t iter=0;
1215     for (Int_t sLayer=2; sLayer>=0;sLayer--){
1216       //for (Int_t dseed=5;dseed<15; dseed+=3){  //loop over central seeding time bins
1217       iter+=1.;
1218       Int_t dseed = 5+Int_t(iter)*3;
1219       // Initialize seeding layers
1220       for (Int_t ilayer=0;ilayer<6;ilayer++){
1221         reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
1222         xcl[ilayer]       = reflayers[ilayer]->GetX();
1223       }
1224       //
1225       Double_t xref                 = (xcl[sLayer+1] + xcl[sLayer+2])*0.5;
1226       AliTRDpropagationLayer& layer0=*reflayers[sLayer+0];
1227       AliTRDpropagationLayer& layer1=*reflayers[sLayer+1];
1228       AliTRDpropagationLayer& layer2=*reflayers[sLayer+2];
1229       AliTRDpropagationLayer& layer3=*reflayers[sLayer+3];
1230       //
1231       Int_t maxn3  = layer3;
1232       for (Int_t icl3=0;icl3<maxn3;icl3++){
1233         AliTRDcluster *cl3 = layer3[icl3];
1234         if (!cl3) continue;
1235         padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2()*12.);
1236         ycl[sLayer+3] = cl3->GetY();
1237         zcl[sLayer+3] = cl3->GetZ();
1238         Float_t yymin0 = ycl[sLayer+3] - 1- kMaxPhi *(xcl[sLayer+3]-xcl[sLayer+0]);
1239         Float_t yymax0 = ycl[sLayer+3] + 1+ kMaxPhi *(xcl[sLayer+3]-xcl[sLayer+0]);
1240         Int_t   maxn0 = layer0;  //
1241         for (Int_t icl0=layer0.Find(yymin0);icl0<maxn0;icl0++){
1242           AliTRDcluster *cl0 = layer0[icl0];
1243           if (!cl0) continue;
1244           if (cl3->IsUsed()&&cl0->IsUsed()) continue;
1245           ycl[sLayer+0] = cl0->GetY();
1246           zcl[sLayer+0] = cl0->GetZ();
1247           if ( ycl[sLayer+0]>yymax0) break;
1248           Double_t tanphi   = (ycl[sLayer+3]-ycl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
1249           if (TMath::Abs(tanphi)>kMaxPhi) continue;
1250           Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
1251           if (TMath::Abs(tantheta)>kMaxTheta) continue;
1252           padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2()*12.);
1253           //
1254           // expected position in 1 layer
1255           Double_t y1exp = ycl[sLayer+0]+(tanphi)  *(xcl[sLayer+1]-xcl[sLayer+0]);
1256           Double_t z1exp = zcl[sLayer+0]+(tantheta)*(xcl[sLayer+1]-xcl[sLayer+0]);
1257           Float_t yymin1 = y1exp - kRoad0y-tanphi;
1258           Float_t yymax1 = y1exp + kRoad0y+tanphi;
1259           Int_t   maxn1  = layer1;  //
1260           //
1261           for (Int_t icl1=layer1.Find(yymin1);icl1<maxn1;icl1++){
1262             AliTRDcluster *cl1 = layer1[icl1];
1263             if (!cl1) continue;
1264             Int_t nusedCl = 0;
1265             if (cl3->IsUsed()) nusedCl++;
1266             if (cl0->IsUsed()) nusedCl++;
1267             if (cl1->IsUsed()) nusedCl++;
1268             if (nusedCl>1) continue;
1269             ycl[sLayer+1] = cl1->GetY();
1270             zcl[sLayer+1] = cl1->GetZ();
1271             if ( ycl[sLayer+1]>yymax1) break;
1272             if (TMath::Abs(ycl[sLayer+1]-y1exp)>kRoad0y+tanphi) continue;
1273             if (TMath::Abs(zcl[sLayer+1]-z1exp)>kRoad0z)        continue;
1274             padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2()*12.);
1275             //
1276             Double_t y2exp  = ycl[sLayer+0]+(tanphi)  *(xcl[sLayer+2]-xcl[sLayer+0])+(ycl[sLayer+1]-y1exp);
1277             Double_t z2exp  = zcl[sLayer+0]+(tantheta)*(xcl[sLayer+2]-xcl[sLayer+0]);
1278             Int_t    index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y,  kRoad1z);
1279             if (index2<=0) continue;
1280             AliTRDcluster *cl2 = (AliTRDcluster*)GetCluster(index2);
1281             padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2()*12.);
1282             ycl[sLayer+2] = cl2->GetY();
1283             zcl[sLayer+2] = cl2->GetZ();
1284             if (TMath::Abs(cl2->GetZ()-z2exp)>kRoad0z)        continue;
1285             //
1286             rieman.Reset();
1287             rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
1288             rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
1289             rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
1290             rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
1291             rieman.Update();
1292             //
1293             // reset fitter
1294             for (Int_t iLayer=0;iLayer<6;iLayer++){
1295               cseed[iLayer].Reset();
1296             }
1297             chi2Z =0.; chi2R=0.;
1298             for (Int_t iLayer=0;iLayer<4;iLayer++){
1299               cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
1300               chi2Z += (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer])*
1301                 (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer]);
1302               cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
1303               cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
1304               chi2R += (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer])*
1305                 (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer]);
1306               cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
1307             }
1308             if (TMath::Sqrt(chi2R)>1./iter) continue;
1309             if (TMath::Sqrt(chi2Z)>7./iter) continue;
1310             //
1311             //
1312             //
1313             Float_t minmax[2]={-100,100};
1314             for (Int_t iLayer=0;iLayer<4;iLayer++){
1315               Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer]*0.5+1 -cseed[sLayer+iLayer].fZref[0];
1316               if (max<minmax[1]) minmax[1]=max;
1317               Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer]*0.5-1 -cseed[sLayer+iLayer].fZref[0];
1318               if (min>minmax[0]) minmax[0]=min;
1319             }
1320             Bool_t isFake = kFALSE;
1321             if (cl0->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
1322             if (cl1->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
1323             if (cl2->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
1324             if (AliTRDReconstructor::StreamLevel()>0){
1325               if ((!isFake) || (icl3%10)==0 ){  //debugging print
1326                 TTreeSRedirector& cstream = *fDebugStreamer;
1327                 cstream<<"Seeds0"<<
1328                   "isFake="<<isFake<<
1329                   "Cl0.="<<cl0<<
1330                   "Cl1.="<<cl1<<
1331                   "Cl2.="<<cl2<<
1332                   "Cl3.="<<cl3<<
1333                   "Xref="<<xref<<
1334                   "X0="<<xcl[sLayer+0]<<
1335                   "X1="<<xcl[sLayer+1]<<
1336                   "X2="<<xcl[sLayer+2]<<
1337                   "X3="<<xcl[sLayer+3]<<
1338                   "Y2exp="<<y2exp<<
1339                   "Z2exp="<<z2exp<<
1340                   "Chi2R="<<chi2R<<
1341                   "Chi2Z="<<chi2Z<<
1342                   "Seed0.="<<&cseed[sLayer+0]<<
1343                   "Seed1.="<<&cseed[sLayer+1]<<
1344                   "Seed2.="<<&cseed[sLayer+2]<<
1345                   "Seed3.="<<&cseed[sLayer+3]<<
1346                   "Zmin="<<minmax[0]<<
1347                   "Zmax="<<minmax[1]<<
1348                   "\n";
1349               }
1350             }
1351
1352             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1353             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1354             //<<<<<<<<<<<<<<<<<<    FIT SEEDING PART                  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1355             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1356             cl[sLayer+0] = cl0;
1357             cl[sLayer+1] = cl1;
1358             cl[sLayer+2] = cl2;
1359             cl[sLayer+3] = cl3;
1360             Bool_t isOK=kTRUE;
1361             for (Int_t jLayer=0;jLayer<4;jLayer++){
1362               cseed[sLayer+jLayer].fTilt = hL[sLayer+jLayer];
1363               cseed[sLayer+jLayer].fPadLength = padlength[sLayer+jLayer];
1364               cseed[sLayer+jLayer].fX0   = xcl[sLayer+jLayer];
1365               for (Int_t iter=0; iter<2; iter++){
1366                 //
1367                 // in iteration 0 we try only one pad-row
1368                 // if quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
1369                 //
1370                 AliTRDseed tseed = cseed[sLayer+jLayer];
1371                 Float_t    roadz  = padlength[sLayer+jLayer]*0.5;
1372                 if (iter>0) roadz = padlength[sLayer+jLayer];
1373                 //
1374                 Float_t quality =10000;
1375                 for (Int_t iTime=2;iTime<20;iTime++){
1376                   AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
1377                   Double_t dxlayer= layer.GetX()-xcl[sLayer+jLayer];
1378                   Double_t zexp   = cl[sLayer+jLayer]->GetZ() ;
1379                   if (iter>0){
1380                     // try 2 pad-rows in second iteration
1381                     zexp  = tseed.fZref[0]+ tseed.fZref[1]*dxlayer;
1382                     if (zexp>cl[sLayer+jLayer]->GetZ()) zexp = cl[sLayer+jLayer]->GetZ()+padlength[sLayer+jLayer]*0.5;
1383                     if (zexp<cl[sLayer+jLayer]->GetZ()) zexp = cl[sLayer+jLayer]->GetZ()-padlength[sLayer+jLayer]*0.5;
1384                   }
1385                   //
1386                   Double_t yexp  =  tseed.fYref[0]+
1387                     tseed.fYref[1]*dxlayer;
1388                   Int_t    index = layer.FindNearestCluster(yexp,zexp,kRoad1y, roadz);
1389                   if (index<=0) continue;
1390                   AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
1391                   //
1392                   tseed.fIndexes[iTime]  = index;
1393                   tseed.fClusters[iTime] = cl;   // register cluster
1394                   tseed.fX[iTime] = dxlayer;     // register cluster
1395                   tseed.fY[iTime] = cl->GetY();  // register cluster
1396                   tseed.fZ[iTime] = cl->GetZ();  // register cluster
1397                 }
1398                 tseed.Update();
1399                 //count the number of clusters and distortions into quality
1400                 Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1401                 Float_t tquality   = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1402                   TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
1403                   2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1404                 if (iter==0 && tseed.IsOK()) {
1405                   cseed[sLayer+jLayer] = tseed;
1406                   quality = tquality;
1407                   if (tquality<5) break;
1408                 }
1409                 if (tseed.IsOK() && tquality<quality)
1410                   cseed[sLayer+jLayer] = tseed;
1411               }
1412               if (!cseed[sLayer+jLayer].IsOK()){
1413                 isOK = kFALSE;
1414                 break;
1415               }
1416               cseed[sLayer+jLayer].CookLabels();
1417               cseed[sLayer+jLayer].UpdateUsed();
1418               nusedCl+= cseed[sLayer+jLayer].fNUsed;
1419               if (nusedCl>25){
1420                 isOK = kFALSE;
1421                 break;
1422               }
1423             }
1424             //
1425             if (!isOK) continue;
1426             nclusters=0;
1427             for (Int_t iLayer=0;iLayer<4;iLayer++){
1428               if (cseed[sLayer+iLayer].IsOK()){
1429                 nclusters+=cseed[sLayer+iLayer].fN2;
1430               }
1431             }
1432             //
1433             // iteration 0
1434             rieman.Reset();
1435             for (Int_t iLayer=0;iLayer<4;iLayer++){
1436               rieman.AddPoint(xcl[sLayer+iLayer],cseed[sLayer+iLayer].fYfitR[0],
1437                               cseed[sLayer+iLayer].fZProb,1,10);
1438             }
1439             rieman.Update();
1440             //
1441             //
1442             chi2R =0; chi2Z=0;
1443             for (Int_t iLayer=0;iLayer<4;iLayer++){
1444               cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
1445               chi2R += (cseed[sLayer+iLayer].fYref[0]-cseed[sLayer+iLayer].fYfitR[0])*
1446                 (cseed[sLayer+iLayer].fYref[0]-cseed[sLayer+iLayer].fYfitR[0]);
1447               cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
1448               cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
1449               chi2Z += (cseed[sLayer+iLayer].fZref[0]- cseed[sLayer+iLayer].fMeanz)*
1450                 (cseed[sLayer+iLayer].fZref[0]- cseed[sLayer+iLayer].fMeanz);
1451               cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
1452             }
1453             Double_t curv = rieman.GetC();
1454             //
1455             // likelihoods
1456             //
1457             Double_t sumda =
1458               TMath::Abs(cseed[sLayer+0].fYfitR[1]- cseed[sLayer+0].fYref[1])+
1459               TMath::Abs(cseed[sLayer+1].fYfitR[1]- cseed[sLayer+1].fYref[1])+
1460               TMath::Abs(cseed[sLayer+2].fYfitR[1]- cseed[sLayer+2].fYref[1])+
1461               TMath::Abs(cseed[sLayer+3].fYfitR[1]- cseed[sLayer+3].fYref[1]);
1462             Double_t likea = TMath::Exp(-sumda*10.6);
1463             Double_t likechi2 = 0.0000000001;
1464             if (chi2R<0.5) likechi2+=TMath::Exp(-TMath::Sqrt(chi2R)*7.73);
1465             Double_t likechi2z = TMath::Exp(-chi2Z*0.088)/TMath::Exp(-chi2Z*0.019);
1466             Double_t likeN    = TMath::Exp(-(72-nclusters)*0.19);
1467             Double_t like     = likea*likechi2*likechi2z*likeN;
1468             //
1469             Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fYref[1]-130*curv)*1.9);
1470             Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fZref[1]-
1471                                                         cseed[sLayer+0].fZref[0]/xcl[sLayer+0])*5.9);
1472             Double_t likePrim  = TMath::Max(likePrimY*likePrimZ,0.0005);
1473
1474             seedquality[registered]  = like;
1475             seedlayer[registered]    = sLayer;
1476             if (TMath::Log(0.000000000000001+like)<-15) continue;
1477             AliTRDseed seedb[6];
1478             for (Int_t iLayer=0;iLayer<6;iLayer++){
1479               seedb[iLayer] = cseed[iLayer];
1480             }
1481             //
1482             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1483             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1484             //<<<<<<<<<<<<<<<   FULL TRACK FIT PART         <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1485             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1486             //
1487             Int_t nlayers            = 0;
1488             Int_t nusedf             = 0;
1489             Int_t findable           = 0;
1490             //
1491             // add new layers  - avoid long extrapolation
1492             //
1493             Int_t tLayer[2]={0,0};
1494             if (sLayer==2) {tLayer[0]=1; tLayer[1]=0;}
1495             if (sLayer==1) {tLayer[0]=5; tLayer[1]=0;}
1496             if (sLayer==0) {tLayer[0]=4; tLayer[1]=5;}
1497             //
1498             for (Int_t iLayer=0;iLayer<2;iLayer++){
1499               Int_t jLayer = tLayer[iLayer];      // set tracking layer
1500               cseed[jLayer].Reset();
1501               cseed[jLayer].fTilt    = hL[jLayer];
1502               cseed[jLayer].fPadLength = padlength[jLayer];
1503               cseed[jLayer].fX0      = xcl[jLayer];
1504               // get pad length and rough cluster
1505               Int_t    indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].fYref[0],
1506                                                                           cseed[jLayer].fZref[0],kRoad2y,kRoad2z);
1507               if (indexdummy<=0) continue;
1508               AliTRDcluster *cldummy = (AliTRDcluster*)GetCluster(indexdummy);
1509               padlength[jLayer]      = TMath::Sqrt(cldummy->GetSigmaZ2()*12.);
1510             }
1511             AliTRDseed::FitRiemanTilt(cseed, kTRUE);
1512             //
1513             for (Int_t iLayer=0;iLayer<2;iLayer++){
1514               Int_t jLayer = tLayer[iLayer];      // set tracking layer
1515               if ( (jLayer==0) && !(cseed[1].IsOK())) continue;  // break not allowed
1516               if ( (jLayer==5) && !(cseed[4].IsOK())) continue;  // break not allowed
1517               Float_t  zexp  = cseed[jLayer].fZref[0];
1518               Double_t zroad =  padlength[jLayer]*0.5+1.;
1519               //
1520               //
1521               for (Int_t iter=0;iter<2;iter++){
1522                 AliTRDseed tseed = cseed[jLayer];
1523                 Float_t quality = 10000;
1524                 for (Int_t iTime=2;iTime<20;iTime++){
1525                   AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
1526                   Double_t dxlayer     = layer.GetX()-xcl[jLayer];
1527                   Double_t yexp        = tseed.fYref[0]+tseed.fYref[1]*dxlayer;
1528                   Float_t  yroad       = kRoad1y;
1529                   Int_t    index = layer.FindNearestCluster(yexp,zexp, yroad, zroad);
1530                   if (index<=0) continue;
1531                   AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
1532                   //
1533                   tseed.fIndexes[iTime]  = index;
1534                   tseed.fClusters[iTime] = cl;   // register cluster
1535                   tseed.fX[iTime] = dxlayer;     // register cluster
1536                   tseed.fY[iTime] = cl->GetY();  // register cluster
1537                   tseed.fZ[iTime] = cl->GetZ();  // register cluster
1538                 }
1539                 tseed.Update();
1540                 if (tseed.IsOK()){
1541                   Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1542                   Float_t tquality   = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1543                     TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
1544                     2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1545                   //
1546                   if (tquality<quality){
1547                     cseed[jLayer]=tseed;
1548                     quality = tquality;
1549                   }
1550                 }
1551                 zroad*=2.;
1552               }
1553               if ( cseed[jLayer].IsOK()){
1554                 cseed[jLayer].CookLabels();
1555                 cseed[jLayer].UpdateUsed();
1556                 nusedf+= cseed[jLayer].fNUsed;
1557                 AliTRDseed::FitRiemanTilt(cseed, kTRUE);
1558               }
1559             }
1560             //
1561             //
1562             // make copy
1563             AliTRDseed bseed[6];
1564             for (Int_t jLayer=0;jLayer<6;jLayer++){
1565               bseed[jLayer] = cseed[jLayer];
1566             }
1567             Float_t lastquality = 10000;
1568             Float_t lastchi2    = 10000;
1569             Float_t chi2        = 1000;
1570
1571             //
1572             for (Int_t iter =0; iter<4;iter++){
1573               //
1574               // sort tracklets according "quality", try to "improve" 4 worst
1575               //
1576               Float_t sumquality = 0;
1577               Float_t squality[6];
1578               Int_t   sortindexes[6];
1579               for (Int_t jLayer=0;jLayer<6;jLayer++){
1580                 if (bseed[jLayer].IsOK()){
1581                   AliTRDseed &tseed = bseed[jLayer];
1582                   Double_t zcor  =  tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
1583                   Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1584                   Float_t tquality  = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1585                     TMath::Abs(tseed.fYfit[0]-(tseed.fYref[0]-zcor))/0.2+
1586                     2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1587                   squality[jLayer] = tquality;
1588                 }
1589                 else  squality[jLayer]=-1;
1590                 sumquality +=squality[jLayer];
1591               }
1592
1593               if (sumquality>=lastquality ||  chi2>lastchi2) break;
1594               lastquality = sumquality;
1595               lastchi2    = chi2;
1596               if (iter>0){
1597                 for (Int_t jLayer=0;jLayer<6;jLayer++){
1598                   cseed[jLayer] = bseed[jLayer];
1599                 }
1600               }
1601               TMath::Sort(6,squality,sortindexes,kFALSE);
1602               //
1603               //
1604               for (Int_t jLayer=5;jLayer>1;jLayer--){
1605                 Int_t bLayer = sortindexes[jLayer];
1606                 AliTRDseed tseed = bseed[bLayer];
1607                 for (Int_t iTime=2;iTime<20;iTime++){
1608                   AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
1609                   Double_t dxlayer= layer.GetX()-xcl[bLayer];
1610                   //
1611                   Double_t zexp  =  tseed.fZref[0];
1612                   Double_t zcor  =  tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
1613                   //
1614                   Float_t  roadz = padlength[bLayer]+1;
1615                   if (TMath::Abs(tseed.fZProb-zexp)> padlength[bLayer]*0.5) {roadz = padlength[bLayer]*0.5;}
1616                   if (tseed.fZfit[1]*tseed.fZref[1]<0) {roadz = padlength[bLayer]*0.5;}
1617                   if (TMath::Abs(tseed.fZProb-zexp)<0.1*padlength[bLayer]) {
1618                     zexp = tseed.fZProb;
1619                     roadz = padlength[bLayer]*0.5;
1620                   }
1621                   //
1622                   Double_t yexp  =  tseed.fYref[0]+
1623                     tseed.fYref[1]*dxlayer-zcor;
1624                   Int_t    index = layer.FindNearestCluster(yexp,zexp,kRoad1y, roadz);
1625                   if (index<=0) continue;
1626                   AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
1627                   //
1628                   tseed.fIndexes[iTime]  = index;
1629                   tseed.fClusters[iTime] = cl;   // register cluster
1630                   tseed.fX[iTime] = dxlayer;     // register cluster
1631                   tseed.fY[iTime] = cl->GetY();  // register cluster
1632                   tseed.fZ[iTime] = cl->GetZ();  // register cluster
1633                 }
1634                 tseed.Update();
1635                 if (tseed.IsOK()) {
1636                   Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
1637                   Double_t zcor  =  tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
1638                   //
1639                   Float_t tquality   = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
1640                     TMath::Abs(tseed.fYfit[0]-(tseed.fYref[0]-zcor))/0.2+
1641                     2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
1642                   //
1643                   if (tquality<squality[bLayer])
1644                     bseed[bLayer] = tseed;
1645                 }
1646               }
1647               chi2 = AliTRDseed::FitRiemanTilt(bseed, kTRUE);
1648             }
1649             //
1650             //
1651             //
1652             nclusters  = 0;
1653             nlayers    = 0;
1654             findable   = 0;
1655             for (Int_t iLayer=0;iLayer<6;iLayer++) {
1656               if (TMath::Abs(cseed[iLayer].fYref[0]/cseed[iLayer].fX0)<0.15)
1657                 findable++;
1658               if (cseed[iLayer].IsOK()){
1659                 nclusters+=cseed[iLayer].fN2;
1660                 nlayers++;
1661               }
1662             }
1663             if (nlayers<3) continue;
1664             rieman.Reset();
1665             for (Int_t iLayer=0;iLayer<6;iLayer++){
1666               if (cseed[iLayer].IsOK()) rieman.AddPoint(xcl[iLayer],cseed[iLayer].fYfitR[0],
1667                                                                    cseed[iLayer].fZProb,1,10);
1668             }
1669             rieman.Update();
1670             //
1671             chi2RF =0;
1672             chi2ZF =0;
1673             for (Int_t iLayer=0;iLayer<6;iLayer++){
1674               if (cseed[iLayer].IsOK()){
1675                 cseed[iLayer].fYref[0] = rieman.GetYat(xcl[iLayer]);
1676                 chi2RF += (cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0])*
1677                   (cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0]);
1678                 cseed[iLayer].fYref[1] = rieman.GetDYat(xcl[iLayer]);
1679                 cseed[iLayer].fZref[0] = rieman.GetZat(xcl[iLayer]);
1680                 chi2ZF += (cseed[iLayer].fZref[0]- cseed[iLayer].fMeanz)*
1681                   (cseed[iLayer].fZref[0]- cseed[iLayer].fMeanz);
1682                 cseed[iLayer].fZref[1] = rieman.GetDZat(xcl[iLayer]);
1683               }
1684             }
1685             chi2RF/=TMath::Max((nlayers-3.),1.);
1686             chi2ZF/=TMath::Max((nlayers-3.),1.);
1687             curv = rieman.GetC();
1688
1689             //
1690
1691             Double_t xref2    = (xcl[2]+xcl[3])*0.5;  // middle of the chamber
1692             Double_t dzmf     = rieman.GetDZat(xref2);
1693             Double_t zmf      = rieman.GetZat(xref2);
1694             //
1695             // fit hyperplane
1696             //
1697             Int_t npointsT =0;
1698             fitterTC.ClearPoints();
1699             fitterT2.ClearPoints();
1700             rieman2.Reset();
1701             for (Int_t iLayer=0; iLayer<6;iLayer++){
1702               if (!cseed[iLayer].IsOK()) continue;
1703               for (Int_t itime=0;itime<25;itime++){
1704                 if (!cseed[iLayer].fUsable[itime]) continue;
1705                 Double_t x   = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2;  // x relative to the midle chamber
1706                 Double_t y   = cseed[iLayer].fY[itime];
1707                 Double_t z   = cseed[iLayer].fZ[itime];
1708                 // ExB correction to the correction
1709                 // tilted rieman
1710                 //
1711                 Double_t uvt[6];
1712                 Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0;      // global x
1713                 //
1714                 Double_t t = 1./(x2*x2+y*y);
1715                 uvt[1]  = t;    // t
1716                 uvt[0]  = 2.*x2*uvt[1];      // u
1717                 //
1718                 uvt[2]  = 2.0*hL[iLayer]*uvt[1];
1719                 uvt[3]  = 2.0*hL[iLayer]*x*uvt[1];
1720                 uvt[4]  = 2.0*(y+hL[iLayer]*z)*uvt[1];
1721                 //
1722                 Double_t error = 2*0.2*uvt[1];
1723                 fitterT2.AddPoint(uvt,uvt[4],error);
1724                 //
1725                 // constrained rieman
1726                 //
1727                 z =cseed[iLayer].fZ[itime];
1728                 uvt[0]  = 2.*x2*t;           // u
1729                 uvt[1]  = 2*hL[iLayer]*x2*uvt[1];
1730                 uvt[2]  = 2*(y+hL[iLayer]*(z-GetZ()))*t;
1731                 fitterTC.AddPoint(uvt,uvt[2],error);
1732                 //
1733                 rieman2.AddPoint(x2,y,z,1,10);
1734                 npointsT++;
1735               }
1736             }
1737             rieman2.Update();
1738             fitterTC.Eval();
1739             fitterT2.Eval();
1740             Double_t rpolz0 = fitterT2.GetParameter(3);
1741             Double_t rpolz1 = fitterT2.GetParameter(4);
1742             //
1743             // linear fitter  - not possible to make boundaries
1744             // non accept non possible z and dzdx combination
1745             //
1746             Bool_t   acceptablez =kTRUE;
1747             for (Int_t iLayer=0; iLayer<6;iLayer++){
1748               if (cseed[iLayer].IsOK()){
1749                 Double_t zT2 =  rpolz0+rpolz1*(xcl[iLayer] - xref2);
1750                 if (TMath::Abs(cseed[iLayer].fZProb-zT2)>padlength[iLayer]*0.5+1)
1751                   acceptablez = kFALSE;
1752               }
1753             }
1754             if (!acceptablez){
1755               fitterT2.FixParameter(3,zmf);
1756               fitterT2.FixParameter(4,dzmf);
1757               fitterT2.Eval();
1758               fitterT2.ReleaseParameter(3);
1759               fitterT2.ReleaseParameter(4);
1760               rpolz0 = fitterT2.GetParameter(3);
1761               rpolz1 = fitterT2.GetParameter(4);
1762             }
1763             //
1764             Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
1765             Double_t chi2TC = fitterTC.GetChisquare()/Float_t(npointsT);
1766             //
1767             Double_t polz1c = fitterTC.GetParameter(2);
1768             Double_t polz0c = polz1c*xref2;
1769             //
1770             Double_t aC     =  fitterTC.GetParameter(0);
1771             Double_t bC     =  fitterTC.GetParameter(1);
1772             Double_t cC     =  aC/TMath::Sqrt(bC*bC+1.);     // curvature
1773             //
1774             Double_t aR     =  fitterT2.GetParameter(0);
1775             Double_t bR     =  fitterT2.GetParameter(1);
1776             Double_t dR     =  fitterT2.GetParameter(2);
1777             Double_t cR     =  1+bR*bR-dR*aR;
1778             Double_t dca    =  0.;
1779             if (cR>0){
1780               dca = -dR/(TMath::Sqrt(1+bR*bR-dR*aR)+TMath::Sqrt(1+bR*bR));
1781               cR  = aR/TMath::Sqrt(cR);
1782             }
1783             //
1784             Double_t chi2ZT2=0, chi2ZTC=0;
1785             for (Int_t iLayer=0; iLayer<6;iLayer++){
1786               if (cseed[iLayer].IsOK()){
1787                 Double_t zT2 =  rpolz0+rpolz1*(xcl[iLayer] - xref2);
1788                 Double_t zTC =  polz0c+polz1c*(xcl[iLayer] - xref2);
1789                 chi2ZT2 += TMath::Abs(cseed[iLayer].fMeanz-zT2);
1790                 chi2ZTC += TMath::Abs(cseed[iLayer].fMeanz-zTC);
1791               }
1792             }
1793             chi2ZT2/=TMath::Max((nlayers-3.),1.);
1794             chi2ZTC/=TMath::Max((nlayers-3.),1.);
1795             //
1796             //
1797             //
1798             AliTRDseed::FitRiemanTilt(cseed, kTRUE);
1799             Float_t sumdaf = 0;
1800             for (Int_t iLayer=0;iLayer<6;iLayer++){
1801               if (cseed[iLayer].IsOK())
1802                 sumdaf += TMath::Abs((cseed[iLayer].fYfit[1]-cseed[iLayer].fYref[1])/cseed[iLayer].fSigmaY2);
1803             }
1804             sumdaf /= Float_t (nlayers-2.);
1805             //
1806             // likelihoods for full track
1807             //
1808             Double_t likezf      = TMath::Exp(-chi2ZF*0.14);
1809             Double_t likechi2C   = TMath::Exp(-chi2TC*0.677);
1810             Double_t likechi2TR  = TMath::Exp(-chi2TR*0.78);
1811             Double_t likeaf      = TMath::Exp(-sumdaf*3.23);
1812             seedquality2[registered] = likezf*likechi2TR*likeaf;
1813 //          Bool_t isGold = kFALSE;
1814 //
1815 //          if (nlayers == 6        && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE;   // gold
1816 //          if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE;   // gold
1817 //          if (isGold &&nusedf<10){
1818 //            for (Int_t jLayer=0;jLayer<6;jLayer++){
1819 //              if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
1820 //                seed[index][jLayer].UseClusters();  //sign gold
1821 //            }
1822 //          }
1823             //
1824             //
1825             //
1826             Int_t index0=0;
1827             if (!cseed[0].IsOK()){
1828               index0 = 1;
1829               if (!cseed[1].IsOK()) index0 = 2;
1830             }
1831             seedparams[registered][0] = cseed[index0].fX0;
1832             seedparams[registered][1] = cseed[index0].fYref[0];
1833             seedparams[registered][2] = cseed[index0].fZref[0];
1834             seedparams[registered][5] = cR;
1835             seedparams[registered][3] = cseed[index0].fX0*cR - TMath::Sin(TMath::ATan(cseed[0].fYref[1]));
1836             seedparams[registered][4] = cseed[index0].fZref[1]/
1837               TMath::Sqrt(1+cseed[index0].fYref[1]*cseed[index0].fYref[1]);
1838             seedparams[registered][6] = ns;
1839             //
1840             //
1841             Int_t labels[12], outlab[24];
1842             Int_t nlab=0;
1843             for (Int_t iLayer=0;iLayer<6;iLayer++){
1844               if (!cseed[iLayer].IsOK()) continue;
1845               if (cseed[iLayer].fLabels[0]>=0) {
1846                 labels[nlab] = cseed[iLayer].fLabels[0];
1847                 nlab++;
1848               }
1849               if (cseed[iLayer].fLabels[1]>=0) {
1850                 labels[nlab] = cseed[iLayer].fLabels[1];
1851                 nlab++;
1852               }
1853             }
1854             Freq(nlab,labels,outlab,kFALSE);
1855             Int_t label = outlab[0];
1856             Int_t frequency  = outlab[1];
1857             for (Int_t iLayer=0;iLayer<6;iLayer++){
1858               cseed[iLayer].fFreq  = frequency;
1859               cseed[iLayer].fC     = cR;
1860               cseed[iLayer].fCC     = cC;
1861               cseed[iLayer].fChi2  = chi2TR;
1862               cseed[iLayer].fChi2Z = chi2ZF;
1863             }
1864             //
1865             if (1||(!isFake)){  //debugging print
1866               Float_t zvertex = GetZ();
1867               TTreeSRedirector& cstream = *fDebugStreamer;
1868               if (AliTRDReconstructor::StreamLevel()>0)
1869               cstream<<"Seeds1"<<
1870                 "isFake="<<isFake<<
1871                 "Vertex="<<zvertex<<
1872                 "Rieman2.="<<&rieman2<<
1873                 "Rieman.="<<&rieman<<
1874                 "Xref="<<xref<<
1875                 "X0="<<xcl[0]<<
1876                 "X1="<<xcl[1]<<
1877                 "X2="<<xcl[2]<<
1878                 "X3="<<xcl[3]<<
1879                 "X4="<<xcl[4]<<
1880                 "X5="<<xcl[5]<<
1881                 "Chi2R="<<chi2R<<
1882                 "Chi2Z="<<chi2Z<<
1883                 "Chi2RF="<<chi2RF<<                          //chi2 of trackletes on full track
1884                 "Chi2ZF="<<chi2ZF<<                          //chi2 z on tracklets on full track
1885                 "Chi2ZT2="<<chi2ZT2<<                        //chi2 z on tracklets on full track  - rieman tilt
1886                 "Chi2ZTC="<<chi2ZTC<<                        //chi2 z on tracklets on full track  - rieman tilt const
1887                 //
1888                 "Chi2TR="<<chi2TR<<                           //chi2 without vertex constrain
1889                 "Chi2TC="<<chi2TC<<                           //chi2 with    vertex constrain
1890                 "C="<<curv<<                                  // non constrained - no tilt correction
1891                 "DR="<<dR<<                                   // DR parameter          - tilt correction
1892                 "DCA="<<dca<<                                 // DCA                   - tilt correction
1893                 "CR="<<cR<<                                   // non constrained curvature - tilt correction
1894                 "CC="<<cC<<                                   // constrained curvature
1895                 "Polz0="<<polz0c<<
1896                 "Polz1="<<polz1c<<
1897                 "RPolz0="<<rpolz0<<
1898                 "RPolz1="<<rpolz1<<
1899                 "Ncl="<<nclusters<<
1900                 "Nlayers="<<nlayers<<
1901                 "NUsedS="<<nusedCl<<
1902                 "NUsed="<<nusedf<<
1903                 "Findable="<<findable<<
1904                 "Like="<<like<<
1905                 "LikePrim="<<likePrim<<
1906                 "Likechi2C="<<likechi2C<<
1907                 "Likechi2TR="<<likechi2TR<<
1908                 "Likezf="<<likezf<<
1909                 "LikeF="<<seedquality2[registered]<<
1910                 "S0.="<<&cseed[0]<<
1911                 "S1.="<<&cseed[1]<<
1912                 "S2.="<<&cseed[2]<<
1913                 "S3.="<<&cseed[3]<<
1914                 "S4.="<<&cseed[4]<<
1915                 "S5.="<<&cseed[5]<<
1916                 "SB0.="<<&seedb[0]<<
1917                 "SB1.="<<&seedb[1]<<
1918                 "SB2.="<<&seedb[2]<<
1919                 "SB3.="<<&seedb[3]<<
1920                 "SB4.="<<&seedb[4]<<
1921                 "SB5.="<<&seedb[5]<<
1922                 "Label="<<label<<
1923                 "Freq="<<frequency<<
1924                 "sLayer="<<sLayer<<
1925                 "\n";
1926             }
1927             if (registered<kMaxSeed-1) {
1928               registered++;
1929               cseed = seed[registered];
1930             }
1931           }// end of loop over layer 1
1932         }  // end of loop over layer 0
1933       }    // end of loop over layer 3
1934     }      // end of loop over seeding time bins
1935     //
1936     // choos best
1937     //
1938     TMath::Sort(registered,seedquality2,sort,kTRUE);
1939     Bool_t signedseed[kMaxSeed];
1940     for (Int_t i=0;i<registered;i++){
1941       signedseed[i]= kFALSE;
1942     }
1943     for (Int_t iter=0; iter<5; iter++){
1944       for (Int_t iseed=0;iseed<registered;iseed++){
1945         Int_t index = sort[iseed];
1946         if (signedseed[index]) continue;
1947         Int_t labelsall[1000];
1948         Int_t nlabelsall=0;
1949         Int_t naccepted=0;;
1950         Int_t sLayer = seedlayer[index];
1951         Int_t ncl   = 0;
1952         Int_t nused = 0;
1953         Int_t nlayers =0;
1954         Int_t findable   = 0;
1955         for (Int_t jLayer=0;jLayer<6;jLayer++){
1956           if (TMath::Abs(seed[index][jLayer].fYref[0]/xcl[jLayer])<0.15)
1957             findable++;
1958           if (seed[index][jLayer].IsOK()){
1959             seed[index][jLayer].UpdateUsed();
1960             ncl   +=seed[index][jLayer].fN2;
1961             nused +=seed[index][jLayer].fNUsed;
1962             nlayers++;
1963             //cooking label
1964             for (Int_t itime=0;itime<25;itime++){
1965               if (seed[index][jLayer].fUsable[itime]){
1966                 naccepted++;
1967                 for (Int_t ilab=0;ilab<3;ilab++){
1968                   Int_t tindex = seed[index][jLayer].fClusters[itime]->GetLabel(ilab);
1969                   if (tindex>=0){
1970                     labelsall[nlabelsall] = tindex;
1971                     nlabelsall++;
1972                   }
1973                 }
1974               }
1975             }
1976           }
1977         }
1978         //
1979         if (nused>30) continue;
1980         //
1981         if (iter==0){
1982           if (nlayers<6) continue;
1983           if (TMath::Log(0.000000001+seedquality2[index])<-5.) continue;   // gold
1984         }
1985         //
1986         if (iter==1){
1987           if (nlayers<findable) continue;
1988           if (TMath::Log(0.000000001+seedquality2[index])<-4.) continue;  //
1989         }
1990         //
1991         //
1992         if (iter==2){
1993           if (nlayers==findable || nlayers==6) continue;
1994           if (TMath::Log(0.000000001+seedquality2[index])<-6.) continue;
1995         }
1996         //
1997         if (iter==3){
1998           if (TMath::Log(0.000000001+seedquality2[index])<-5.) continue;
1999         }
2000         //
2001         if (iter==4){
2002           if (TMath::Log(0.000000001+seedquality2[index])-nused/(nlayers-3.)<-15.) continue;
2003         }
2004         //
2005         signedseed[index] = kTRUE;
2006         //
2007         Int_t labels[1000], outlab[1000];
2008         Int_t nlab=0;
2009         for (Int_t iLayer=0;iLayer<6;iLayer++){
2010           if (seed[index][iLayer].IsOK()){
2011             if (seed[index][iLayer].fLabels[0]>=0) {
2012               labels[nlab] = seed[index][iLayer].fLabels[0];
2013               nlab++;
2014             }
2015             if (seed[index][iLayer].fLabels[1]>=0) {
2016               labels[nlab] = seed[index][iLayer].fLabels[1];
2017               nlab++;
2018             }
2019           }
2020         }
2021         Freq(nlab,labels,outlab,kFALSE);
2022         Int_t label  = outlab[0];
2023         Int_t frequency  = outlab[1];
2024         Freq(nlabelsall,labelsall,outlab,kFALSE);
2025         Int_t label1 = outlab[0];
2026         Int_t label2 = outlab[2];
2027         Float_t fakeratio = (naccepted-outlab[1])/Float_t(naccepted);
2028         Float_t ratio = Float_t(nused)/Float_t(ncl);
2029         if (ratio<0.25){
2030           for (Int_t jLayer=0;jLayer<6;jLayer++){
2031             if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.2 )
2032               seed[index][jLayer].UseClusters();  //sign gold
2033           }
2034         }
2035         //
2036         Int_t eventNr = esd->GetEventNumber();
2037         TTreeSRedirector& cstream = *fDebugStreamer;
2038         //
2039         // register seed
2040         //
2041         AliTRDtrack * track = RegisterSeed(seed[index],seedparams[index]);
2042         AliTRDtrack dummy;
2043         if (!track) track=&dummy;
2044         else{
2045           AliESDtrack esdtrack;
2046           esdtrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
2047           esdtrack.SetLabel(label);
2048           esd->AddTrack(&esdtrack);
2049           TTreeSRedirector& cstream = *fDebugStreamer;
2050           if (AliTRDReconstructor::StreamLevel()>0)
2051             cstream<<"Tracks"<<
2052               "EventNr="<<eventNr<<
2053               "ESD.="<<&esdtrack<<
2054               "trd.="<<track<<
2055               "trdback.="<<track<<
2056               "\n";
2057         }
2058         if (AliTRDReconstructor::StreamLevel()>0)
2059           cstream<<"Seeds2"<<
2060           "Iter="<<iter<<
2061           "Track.="<<track<<
2062           "Like="<<seedquality[index]<<
2063           "LikeF="<<seedquality2[index]<<
2064           "S0.="<<&seed[index][0]<<
2065           "S1.="<<&seed[index][1]<<
2066           "S2.="<<&seed[index][2]<<
2067           "S3.="<<&seed[index][3]<<
2068           "S4.="<<&seed[index][4]<<
2069           "S5.="<<&seed[index][5]<<
2070           "Label="<<label<<
2071           "Label1="<<label1<<
2072           "Label2="<<label2<<
2073           "FakeRatio="<<fakeratio<<
2074           "Freq="<<frequency<<
2075           "Ncl="<<ncl<<
2076           "Nlayers="<<nlayers<<
2077           "Findable="<<findable<<
2078           "NUsed="<<nused<<
2079           "sLayer="<<sLayer<<
2080           "EventNr="<<eventNr<<
2081           "\n";
2082       }
2083     }
2084   }        // end of loop over sectors
2085   delete [] pseed;
2086 }
2087
2088 //_____________________________________________________________________________
2089 Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *ClusterTree) const
2090 {
2091   //
2092   // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
2093   // from the file. The names of the cluster tree and branches
2094   // should match the ones used in AliTRDclusterizer::WriteClusters()
2095   //
2096   Int_t nsize = Int_t(ClusterTree->GetTotBytes()/(sizeof(AliTRDcluster)));
2097   TObjArray *clusterArray = new TObjArray(nsize+1000);
2098
2099   TBranch *branch=ClusterTree->GetBranch("TRDcluster");
2100   if (!branch) {
2101     Error("ReadClusters","Can't get the branch !");
2102     return 1;
2103   }
2104   branch->SetAddress(&clusterArray);
2105
2106   Int_t nEntries = (Int_t) ClusterTree->GetEntries();
2107   //  printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
2108
2109   // Loop through all entries in the tree
2110   Int_t nbytes = 0;
2111   AliTRDcluster *c = 0;
2112   //  printf("\n");
2113   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2114
2115     // Import the tree
2116     nbytes += ClusterTree->GetEvent(iEntry);
2117
2118     // Get the number of points in the detector
2119     Int_t nCluster = clusterArray->GetEntriesFast();
2120 //    printf("\r Read %d clusters from entry %d", nCluster, iEntry);
2121
2122     // Loop through all TRD digits
2123     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
2124       c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
2125       AliTRDcluster *co = c;
2126       array->AddLast(co);
2127       //      delete clusterArray->RemoveAt(iCluster);
2128       clusterArray->RemoveAt(iCluster);
2129     }
2130   }
2131 //   cout<<"Allocated"<<nsize<<"\tLoaded"<<array->GetEntriesFast()<<"\n";
2132
2133   delete clusterArray;
2134
2135   return 0;
2136 }
2137
2138 //__________________________________________________________________
2139 Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint& p) const
2140 {
2141   //
2142   // Get track space point with index i
2143   // Origin: C.Cheshkov
2144   //
2145
2146   AliTRDcluster *cl = (AliTRDcluster*)fClusters->UncheckedAt(index);
2147   Int_t  idet = cl->GetDetector();
2148   Int_t  isector = fGeom->GetSector(idet);
2149   Int_t  ichamber= fGeom->GetChamber(idet);
2150   Int_t  iplan   = fGeom->GetPlane(idet);
2151   Double_t local[3];
2152   local[0]=GetX(isector,iplan,cl->GetLocalTimeBin());
2153   local[1]=cl->GetY();
2154   local[2]=cl->GetZ();
2155   Double_t global[3];
2156   fGeom->RotateBack(idet,local,global);
2157   p.SetXYZ(global[0],global[1],global[2]);
2158   AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
2159   switch (iplan) {
2160   case 0:
2161     iLayer = AliAlignObj::kTRD1;
2162     break;
2163   case 1:
2164     iLayer = AliAlignObj::kTRD2;
2165     break;
2166   case 2:
2167     iLayer = AliAlignObj::kTRD3;
2168     break;
2169   case 3:
2170     iLayer = AliAlignObj::kTRD4;
2171     break;
2172   case 4:
2173     iLayer = AliAlignObj::kTRD5;
2174     break;
2175   case 5:
2176     iLayer = AliAlignObj::kTRD6;
2177     break;
2178   };
2179   Int_t modId = isector*fGeom->Ncham()+ichamber;
2180   UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
2181   p.SetVolumeID(volid);
2182
2183   return kTRUE;
2184
2185 }
2186
2187 //__________________________________________________________________
2188 void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const
2189 {
2190   //
2191   // This cooks a label. Mmmmh, smells good...
2192   //
2193
2194   Int_t label=123456789, index, i, j;
2195   Int_t ncl=pt->GetNumberOfClusters();
2196   const Int_t kRange = fTrSec[0]->GetOuterTimeBin()+1;
2197
2198   Bool_t labelAdded;
2199
2200   //  Int_t s[kRange][2];
2201   Int_t **s = new Int_t* [kRange];
2202   for (i=0; i<kRange; i++) {
2203     s[i] = new Int_t[2];
2204   }
2205   for (i=0; i<kRange; i++) {
2206     s[i][0]=-1;
2207     s[i][1]=0;
2208   }
2209
2210   Int_t t0,t1,t2;
2211   for (i=0; i<ncl; i++) {
2212     index=pt->GetClusterIndex(i);
2213     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2214     t0=c->GetLabel(0);
2215     t1=c->GetLabel(1);
2216     t2=c->GetLabel(2);
2217   }
2218
2219   for (i=0; i<ncl; i++) {
2220     index=pt->GetClusterIndex(i);
2221     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2222     for (Int_t k=0; k<3; k++) {
2223       label=c->GetLabel(k);
2224       labelAdded=kFALSE; j=0;
2225       if (label >= 0) {
2226         while ( (!labelAdded) && ( j < kRange ) ) {
2227           if (s[j][0]==label || s[j][1]==0) {
2228             s[j][0]=label;
2229             s[j][1]=s[j][1]+1;
2230             labelAdded=kTRUE;
2231           }
2232           j++;
2233         }
2234       }
2235     }
2236   }
2237
2238   Int_t max=0;
2239   label = -123456789;
2240
2241   for (i=0; i<kRange; i++) {
2242     if (s[i][1]>max) {
2243       max=s[i][1]; label=s[i][0];
2244     }
2245   }
2246
2247   for (i=0; i<kRange; i++) {
2248     delete []s[i];
2249   }
2250
2251   delete []s;
2252
2253   if ((1.- Float_t(max)/ncl) > wrong) label=-label;
2254
2255   pt->SetLabel(label);
2256
2257 }
2258
2259
2260 //__________________________________________________________________
2261 void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const
2262 {
2263   //
2264   // Use clusters, but don't abuse them!
2265   //
2266   const Float_t kmaxchi2 =18;
2267   const Float_t kmincl   =10;
2268   AliTRDtrack * track  = (AliTRDtrack*)t;
2269   //
2270   Int_t ncl=t->GetNumberOfClusters();
2271   for (Int_t i=from; i<ncl; i++) {
2272     Int_t index = t->GetClusterIndex(i);
2273     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
2274     //
2275     Int_t iplane = fGeom->GetPlane(c->GetDetector());
2276     if (track->fTracklets[iplane].GetChi2()>kmaxchi2) continue;
2277     if (track->fTracklets[iplane].GetN()<kmincl) continue;
2278     if (!(c->IsUsed())) c->Use();
2279   }
2280 }
2281
2282
2283 //_____________________________________________________________________
2284 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
2285 {
2286   // Parametrised "expected" error of the cluster reconstruction in Y
2287
2288   Double_t s = 0.08 * 0.08;
2289   return s;
2290 }
2291
2292 //_____________________________________________________________________
2293 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
2294 {
2295   // Parametrised "expected" error of the cluster reconstruction in Z
2296
2297   Double_t s = 9 * 9 /12.;
2298   return s;
2299 }
2300
2301 //_____________________________________________________________________
2302 Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
2303 {
2304   //
2305   // Returns radial position which corresponds to time bin <localTB>
2306   // in tracking sector <sector> and plane <plane>
2307   //
2308
2309   Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
2310   Int_t pl = fTrSec[sector]->GetLayerNumber(index);
2311   return fTrSec[sector]->GetLayer(pl)->GetX();
2312
2313 }
2314
2315
2316 //_______________________________________________________
2317 AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x,
2318                                                               Double_t dx, Double_t rho, Double_t radLength, Int_t tbIndex, Int_t plane)
2319 {
2320   //
2321   // AliTRDpropagationLayer constructor
2322   //
2323
2324   fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = radLength;
2325   fClusters = NULL; fIndex = NULL; fTimeBinIndex = tbIndex;
2326   fPlane = plane;
2327
2328   for(Int_t i=0; i < (Int_t) kZones; i++) {
2329     fZc[i]=0; fZmax[i] = 0;
2330   }
2331
2332   fYmax = 0;
2333
2334   if(fTimeBinIndex >= 0) {
2335     fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
2336     fIndex = new UInt_t[kMaxClusterPerTimeBin];
2337   }
2338
2339   for (Int_t i=0;i<5;i++) fIsHole[i] = kFALSE;
2340   fHole = kFALSE;
2341   fHoleZc = 0;
2342   fHoleZmax = 0;
2343   fHoleYc = 0;
2344   fHoleYmax = 0;
2345   fHoleRho = 0;
2346   fHoleX0 = 0;
2347
2348 }
2349
2350 //_______________________________________________________
2351 void AliTRDtracker::AliTRDpropagationLayer::SetHole(
2352           Double_t Zmax, Double_t Ymax, Double_t rho,
2353           Double_t radLength, Double_t Yc, Double_t Zc)
2354 {
2355   //
2356   // Sets hole in the layer
2357   //
2358   fHole = kTRUE;
2359   fHoleZc = Zc;
2360   fHoleZmax = Zmax;
2361   fHoleYc = Yc;
2362   fHoleYmax = Ymax;
2363   fHoleRho = rho;
2364   fHoleX0 = radLength;
2365 }
2366
2367
2368 //_______________________________________________________
2369 AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs)
2370 {
2371   //
2372   // AliTRDtrackingSector Constructor
2373   //
2374   AliTRDpadPlane *padPlane = 0;
2375
2376   fGeom = geo;
2377   fGeomSector = gs;
2378   fN = 0;
2379   //
2380   // get holes description from geometry
2381   Bool_t holes[AliTRDgeometry::kNcham];
2382   //printf("sector\t%d\t",gs);
2383   for (Int_t icham=0; icham<AliTRDgeometry::kNcham;icham++){
2384     holes[icham] = fGeom->IsHole(0,icham,gs);
2385     //printf("%d",holes[icham]);
2386   }
2387   //printf("\n");
2388
2389   for(UInt_t i=0; i < kMaxTimeBinIndex; i++) fTimeBinIndex[i] = -1;
2390
2391
2392   AliTRDpropagationLayer* ppl;
2393
2394   Double_t x, dx, rho, radLength;
2395   //  Int_t    steps;
2396
2397   // add layers for each of the planes
2398   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
2399   //Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2400
2401   Int_t    tbIndex;
2402   const Int_t  kNchambers = AliTRDgeometry::Ncham();
2403   Double_t  ymax = 0;
2404   Double_t ymaxsensitive=0;
2405   Double_t *zc = new Double_t[kNchambers];
2406   Double_t *zmax = new Double_t[kNchambers];
2407   Double_t *zmaxsensitive = new Double_t[kNchambers];
2408
2409   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
2410   if (!commonParam)
2411   {
2412     printf("<AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector> ");
2413     printf("Could not get common params\n");
2414     return;
2415   }
2416
2417   for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
2418
2419     ymax          = fGeom->GetChamberWidth(plane)/2.;
2420     // Modidified for new pad plane class, 22.04.05 (C.B.)
2421     padPlane = commonParam->GetPadPlane(plane,0);
2422     ymaxsensitive = (padPlane->GetColSize(1)*padPlane->GetNcols()-4)/2.;
2423     for(Int_t ch = 0; ch < kNchambers; ch++) {
2424       zmax[ch] = fGeom->GetChamberLength(plane,ch)/2;
2425       //
2426       // Modidified for new pad plane class, 22.04.05 (C.B.)
2427       Float_t pad = padPlane->GetRowSize(1);
2428       Float_t row0 = commonParam->GetRow0(plane,ch,0);
2429       Int_t nPads = commonParam->GetRowMax(plane,ch,0);
2430       zmaxsensitive[ch] = Float_t(nPads)*pad/2.;
2431       zc[ch] = -(pad * nPads)/2 + row0;
2432     }
2433
2434     dx  = AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
2435         / AliTRDcalibDB::Instance()->GetSamplingFrequency();
2436     rho = 0.00295 * 0.85; radLength = 11.0;
2437
2438     Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(plane);
2439     //Double_t xbottom = x0 - dxDrift;
2440     //Double_t xtop = x0 + dxAmp;
2441     //
2442     Int_t nTimeBins =  AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
2443     for (Int_t iTime = 0; iTime<nTimeBins; iTime++){
2444       Double_t xlayer  = iTime*dx - dxAmp;
2445       //if (xlayer<0) xlayer=dxAmp/2.;
2446       x = x0 - xlayer;
2447       //
2448       tbIndex = CookTimeBinIndex(plane, iTime);
2449       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex, plane);
2450       ppl->SetYmax(ymax,ymaxsensitive);
2451       ppl->SetZ(zc, zmax, zmaxsensitive);
2452       ppl->SetHoles(holes);
2453       InsertLayer(ppl);
2454     }
2455   }
2456
2457   MapTimeBinLayers();
2458   delete [] zc;
2459   delete [] zmax;
2460   delete [] zmaxsensitive;
2461
2462 }
2463
2464 //______________________________________________________
2465
2466 Int_t  AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t localTB) const
2467 {
2468   //
2469   // depending on the digitization parameters calculates "global"
2470   // time bin index for timebin <localTB> in plane <plane>
2471   //
2472   //
2473   Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
2474   Int_t gtb = (plane+1) * tbPerPlane - localTB -1;
2475   if (localTB<0) return -1;
2476   if (gtb<0) return -1;
2477   return gtb;
2478 }
2479
2480 //______________________________________________________
2481
2482 void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers()
2483 {
2484   //
2485   // For all sensitive time bins sets corresponding layer index
2486   // in the array fTimeBins
2487   //
2488
2489   Int_t index;
2490
2491   for(Int_t i = 0; i < fN; i++) {
2492     index = fLayers[i]->GetTimeBinIndex();
2493
2494     //    printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX());
2495
2496     if(index < 0) continue;
2497     if(index >= (Int_t) kMaxTimeBinIndex) {
2498       printf("*** AliTRDtracker::MapTimeBinLayers: \n");
2499       printf("    index %d exceeds allowed maximum of %d!\n",
2500              index, kMaxTimeBinIndex-1);
2501       continue;
2502     }
2503     fTimeBinIndex[index] = i;
2504   }
2505 }
2506
2507
2508 //______________________________________________________
2509
2510
2511 Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const
2512 {
2513   //
2514   // Returns the number of time bin which in radial position is closest to <x>
2515   //
2516
2517   if(x >= fLayers[fN-1]->GetX()) return fN-1;
2518   if(x <= fLayers[0]->GetX()) return 0;
2519
2520   Int_t b=0, e=fN-1, m=(b+e)/2;
2521   for (; b<e; m=(b+e)/2) {
2522     if (x > fLayers[m]->GetX()) b=m+1;
2523     else e=m;
2524   }
2525   if(TMath::Abs(x - fLayers[m]->GetX()) >
2526      TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
2527   else return m;
2528
2529 }
2530
2531 //______________________________________________________
2532
2533 Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const
2534 {
2535   //
2536   // Returns number of the innermost SENSITIVE propagation layer
2537   //
2538
2539   return GetLayerNumber(0);
2540 }
2541
2542 //______________________________________________________
2543
2544 Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const
2545 {
2546   //
2547   // Returns number of the outermost SENSITIVE time bin
2548   //
2549
2550   return GetLayerNumber(GetNumberOfTimeBins() - 1);
2551 }
2552
2553 //______________________________________________________
2554
2555 Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const
2556 {
2557   //
2558   // Returns number of SENSITIVE time bins
2559   //
2560
2561   Int_t tb, layer;
2562   for(tb = kMaxTimeBinIndex-1; tb >=0; tb--) {
2563     layer = GetLayerNumber(tb);
2564     if(layer>=0) break;
2565   }
2566   return tb+1;
2567 }
2568
2569 //______________________________________________________
2570
2571 void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
2572 {
2573   //
2574   // Insert layer <pl> in fLayers array.
2575   // Layers are sorted according to X coordinate.
2576
2577   if ( fN == ((Int_t) kMaxLayersPerSector)) {
2578     printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
2579     return;
2580   }
2581   if (fN==0) {fLayers[fN++] = pl; return;}
2582   Int_t i=Find(pl->GetX());
2583
2584   memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
2585   fLayers[i]=pl; fN++;
2586
2587 }
2588
2589 //______________________________________________________
2590
2591 Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const
2592 {
2593   //
2594   // Returns index of the propagation layer nearest to X
2595   //
2596
2597   if (x <= fLayers[0]->GetX()) return 0;
2598   if (x > fLayers[fN-1]->GetX()) return fN;
2599   Int_t b=0, e=fN-1, m=(b+e)/2;
2600   for (; b<e; m=(b+e)/2) {
2601     if (x > fLayers[m]->GetX()) b=m+1;
2602     else e=m;
2603   }
2604   return m;
2605 }
2606
2607
2608
2609
2610
2611 //______________________________________________________
2612 void AliTRDtracker::AliTRDpropagationLayer::SetZ(Double_t* center, Double_t *w, Double_t *wsensitive )
2613 {
2614   //
2615   // set centers and the width of sectors
2616   for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
2617     fZc[icham] = center[icham];
2618     fZmax[icham] = w[icham];
2619     fZmaxSensitive[icham] = wsensitive[icham];
2620     //   printf("chamber\t%d\tzc\t%f\tzmax\t%f\tzsens\t%f\n",icham,fZc[icham],fZmax[icham],fZmaxSensitive[icham]);
2621   }
2622 }
2623 //______________________________________________________
2624
2625 void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
2626 {
2627   //
2628   // set centers and the width of sectors
2629   fHole = kFALSE;
2630   for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
2631     fIsHole[icham] = holes[icham];
2632     if (holes[icham]) fHole = kTRUE;
2633   }
2634 }
2635
2636
2637
2638
2639
2640 //______________________________________________________
2641
2642 void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c,
2643                                                           UInt_t index) {
2644
2645 // Insert cluster in cluster array.
2646 // Clusters are sorted according to Y coordinate.
2647
2648   if(fTimeBinIndex < 0) {
2649     printf("*** attempt to insert cluster into non-sensitive time bin!\n");
2650     return;
2651   }
2652
2653   if (fN== (Int_t) kMaxClusterPerTimeBin) {
2654     printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n");
2655     return;
2656   }
2657   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
2658   Int_t i=Find(c->GetY());
2659   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
2660   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
2661   fIndex[i]=index; fClusters[i]=c; fN++;
2662 }
2663
2664 //______________________________________________________
2665
2666 Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const {
2667
2668 // Returns index of the cluster nearest in Y
2669
2670   if (fN<=0) return 0;
2671   if (y <= fClusters[0]->GetY()) return 0;
2672   if (y > fClusters[fN-1]->GetY()) return fN;
2673   Int_t b=0, e=fN-1, m=(b+e)/2;
2674   for (; b<e; m=(b+e)/2) {
2675     if (y > fClusters[m]->GetY()) b=m+1;
2676     else e=m;
2677   }
2678   return m;
2679 }
2680
2681 Int_t AliTRDtracker::AliTRDpropagationLayer::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad, Float_t maxroadz) const
2682 {
2683   //
2684   // Returns index of the cluster nearest to the given y,z
2685   //
2686   Int_t index = -1;
2687   Int_t maxn = fN;
2688   Float_t mindist = maxroad;
2689   //
2690   for (Int_t i=Find(y-maxroad); i<maxn; i++) {
2691     AliTRDcluster* c=(AliTRDcluster*)(fClusters[i]);
2692     Float_t ycl = c->GetY();
2693     //
2694     if (ycl > y+maxroad) break;
2695     if (TMath::Abs(c->GetZ()-z) > maxroadz) continue;
2696     if (TMath::Abs(ycl-y)<mindist){
2697       mindist = TMath::Abs(ycl-y);
2698       index = fIndex[i];
2699     }
2700   }
2701   return index;
2702 }
2703
2704
2705 //---------------------------------------------------------
2706
2707 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
2708 //
2709 //  Returns correction factor for tilted pads geometry
2710 //
2711   Int_t det = c->GetDetector();
2712   Int_t plane = fGeom->GetPlane(det);
2713   AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(plane,0);
2714   Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
2715
2716   if(fNoTilt) h01 = 0;
2717   return h01;
2718 }
2719
2720
2721 void AliTRDtracker::CookdEdxTimBin(AliTRDtrack& TRDtrack)
2722 {
2723   // *** ADDED TO GET MORE INFORMATION FOR TRD PID  ---- PS
2724   // This is setting fdEdxPlane and fTimBinPlane
2725   // Sums up the charge in each plane for track TRDtrack and also get the
2726   // Time bin for Max. Cluster
2727   // Prashant Shukla (shukla@physi.uni-heidelberg.de)
2728
2729   Double_t  clscharge[kNPlane], maxclscharge[kNPlane];
2730   Int_t  nCluster[kNPlane], timebin[kNPlane];
2731
2732   //Initialization of cluster charge per plane.
2733   for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
2734     clscharge[iPlane] = 0.0;
2735     nCluster[iPlane] = 0;
2736     timebin[iPlane] = -1;
2737     maxclscharge[iPlane] = 0.0;
2738   }
2739
2740   // Loop through all clusters associated to track TRDtrack
2741   Int_t nClus = TRDtrack.GetNumberOfClusters();  // from Kalmantrack
2742   for (Int_t iClus = 0; iClus < nClus; iClus++) {
2743     Double_t charge = TRDtrack.GetClusterdQdl(iClus);
2744     Int_t index = TRDtrack.GetClusterIndex(iClus);
2745     AliTRDcluster *pTRDcluster = (AliTRDcluster *) GetCluster(index);
2746     if (!pTRDcluster) continue;
2747     Int_t tb = pTRDcluster->GetLocalTimeBin();
2748     if (!tb) continue;
2749     Int_t detector = pTRDcluster->GetDetector();
2750     Int_t iPlane   = fGeom->GetPlane(detector);
2751     clscharge[iPlane] = clscharge[iPlane]+charge;
2752     if(charge > maxclscharge[iPlane]) {
2753       maxclscharge[iPlane] = charge;
2754       timebin[iPlane] = tb;
2755     }
2756     nCluster[iPlane]++;
2757   } // end of loop over cluster
2758
2759   // Setting the fdEdxPlane and fTimBinPlane variabales
2760   Double_t totalCharge = 0;
2761   for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
2762     // Quality control of TRD track.
2763     if (nCluster[iPlane]<= 5) {
2764       clscharge[iPlane]=0.0;
2765       timebin[iPlane]=-1;
2766     }
2767     if (nCluster[iPlane]) clscharge[iPlane] /= nCluster[iPlane];
2768     TRDtrack.SetPIDsignals(clscharge[iPlane], iPlane);
2769     TRDtrack.SetPIDTimBin(timebin[iPlane], iPlane);
2770     totalCharge= totalCharge+clscharge[iPlane];
2771   }
2772   //  Int_t i;
2773   //  Int_t nc=TRDtrack.GetNumberOfClusters();
2774   //  Float_t dedx=0;
2775   //  for (i=0; i<nc; i++) dedx += TRDtrack.GetClusterdQdl(i);
2776   //  dedx /= nc;
2777   //  for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
2778   //    TRDtrack.SetPIDsignals(dedx, iPlane);
2779   //    TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
2780   //  }
2781
2782 } // end of function
2783
2784
2785 Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters,AliTRDtracklet&tracklet)
2786 {
2787   //
2788   //
2789   //  try to find nearest clusters to the track in timebins from t0 to t1
2790   //
2791   //
2792   //
2793   // correction coeficients   - depends on TRD parameters  - to be changed according it
2794   //
2795
2796   Double_t x[100],yt[100],zt[100];
2797   Double_t xmean=0;   //reference x
2798   Double_t dz[10][100],dy[10][100];
2799   Float_t zmean[100], nmean[100];
2800   Int_t    clfound=0;
2801   Int_t    indexes[10][100];    // indexes of the clusters in the road
2802   AliTRDcluster *cl[10][100];   // pointers to the clusters in the road
2803   Int_t    best[10][100];       // index of best matching cluster
2804   //
2805   //
2806
2807   for (Int_t it=0;it<=t1-t0; it++){
2808     x[it]=0;
2809     yt[it]=0;
2810     zt[it]=0;
2811     clusters[it+t0]=-2;
2812     zmean[it]=0;
2813     nmean[it]=0;
2814     //
2815     for (Int_t ih=0;ih<10;ih++){
2816       indexes[ih][it]=-2;              //reset indexes1
2817       cl[ih][it]=0;
2818       dz[ih][it]=-100;
2819       dy[ih][it]=-100;
2820       best[ih][it]=0;
2821     }
2822   }
2823   //
2824   Double_t x0 = track->GetX();
2825   Double_t sigmaz = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
2826   Int_t nall=0;
2827   Int_t nfound=0;
2828   Double_t h01 =0;
2829   Int_t plane =-1;
2830   Int_t detector =-1;
2831   Float_t padlength=0;
2832   AliTRDtrack track2(*track);
2833   Float_t snpy = track->GetSnp();
2834   Float_t tany = TMath::Sqrt(snpy*snpy/(1.-snpy*snpy));
2835   if (snpy<0) tany*=-1;
2836   //
2837   Double_t sy2=ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
2838   Double_t sz2=ExpectedSigmaZ2(x0,track->GetTgl());
2839   Double_t road = 15.*sqrt(track->GetSigmaY2() + sy2);
2840   if (road>6.) road=6.;
2841
2842   //
2843   for (Int_t it=0;it<t1-t0;it++){
2844     Double_t maxChi2[2]={fgkMaxChi2,fgkMaxChi2};
2845     AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(it+t0));
2846     if (timeBin==0) continue;  // no indexes1
2847     Int_t maxn = timeBin;
2848     x[it] = timeBin.GetX();
2849     track2.PropagateTo(x[it]);
2850     yt[it] = track2.GetY();
2851     zt[it] = track2.GetZ();
2852
2853     Double_t  y=yt[it],z=zt[it];
2854     Double_t chi2 =1000000;
2855     nall++;
2856     //
2857     // find 2 nearest cluster at given time bin
2858     //
2859     //
2860     for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
2861       AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
2862       h01 = GetTiltFactor(c);
2863       if (plane<0){
2864         Int_t det = c->GetDetector();
2865         plane = fGeom->GetPlane(det);
2866         padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
2867       }
2868       //      if (c->GetLocalTimeBin()==0) continue;
2869       if (c->GetY() > y+road) break;
2870       if((c->GetZ()-z)*(c->GetZ()-z) > 12. * sz2) continue;
2871
2872       Double_t dist = TMath::Abs(c->GetZ()-z);
2873       if (dist> (0.5*padlength+6.*sigmaz)) continue;   // 6 sigma boundary cut
2874       Double_t cost = 0;
2875       //
2876       if (dist> (0.5*padlength-sigmaz)){   //  sigma boundary cost function
2877         cost =  (dist-0.5*padlength)/(2.*sigmaz);
2878         if (cost>-1) cost= (cost+1.)*(cost+1.);
2879         else cost=0;
2880       }
2881       //      Int_t label = TMath::Abs(track->GetLabel());
2882       //      if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
2883       chi2=track2.GetPredictedChi2(c,h01)+cost;
2884       //
2885       clfound++;
2886       if (chi2 > maxChi2[1]) continue;
2887       detector = c->GetDetector();
2888
2889       for (Int_t ih=2;ih<9; ih++){  //store the clusters in the road
2890         if (cl[ih][it]==0){
2891           cl[ih][it] = c;
2892           indexes[ih][it] =timeBin.GetIndex(i);   // index - 9 - reserved for outliers
2893           break;
2894         }
2895       }
2896       //
2897       if (chi2 <maxChi2[0]){
2898         maxChi2[1]     = maxChi2[0];
2899         maxChi2[0]     = chi2;
2900         indexes[1][it] = indexes[0][it];
2901         cl[1][it]      = cl[0][it];
2902         indexes[0][it] = timeBin.GetIndex(i);
2903         cl[0][it]      = c;
2904         continue;
2905       }
2906       maxChi2[1]=chi2;
2907       cl[1][it] = c;
2908       indexes[1][it] =timeBin.GetIndex(i);
2909     }
2910     if (cl[0][it]){
2911       nfound++;
2912       xmean += x[it];
2913     }
2914   }
2915   //
2916   if (nfound<4) return 0;
2917   xmean /=Float_t(nfound);     // middle x
2918   track2.PropagateTo(xmean);   // propagate track to the center
2919   //
2920   // choose one of the variants
2921   //
2922   Int_t changes[10];
2923   Float_t sumz      = 0;
2924   Float_t sum       = 0;
2925   Double_t sumdy    = 0;
2926   Double_t sumdy2   = 0;
2927   Double_t sumx     = 0;
2928   Double_t sumxy    = 0;
2929   Double_t sumx2    = 0;
2930   Double_t mpads    = 0;
2931   //
2932   Int_t   ngood[10];
2933   Int_t   nbad[10];
2934   //
2935   Double_t meanz[10];
2936   Double_t moffset[10];    // mean offset
2937   Double_t mean[10];       // mean value
2938   Double_t angle[10];      // angle
2939   //
2940   Double_t smoffset[10];   // sigma of mean offset
2941   Double_t smean[10];      // sigma of mean value
2942   Double_t sangle[10];     // sigma of angle
2943   Double_t smeanangle[10]; // correlation
2944   //
2945   Double_t sigmas[10];
2946   Double_t tchi2s[10];      // chi2s for tracklet
2947   //
2948   // calculate zmean
2949   //
2950   for (Int_t it=0;it<t1-t0;it++){
2951     if (!cl[0][it]) continue;
2952     for (Int_t dt=-3;dt<=3;dt++){
2953       if (it+dt<0) continue;
2954       if (it+dt>t1-t0) continue;
2955       if (!cl[0][it+dt]) continue;
2956       zmean[it]+=cl[0][it+dt]->GetZ();
2957       nmean[it]+=1.;
2958     }
2959     zmean[it]/=nmean[it];
2960   }
2961   //
2962   for (Int_t it=0; it<t1-t0;it++){
2963     best[0][it]=0;
2964     for (Int_t ih=0;ih<10;ih++){
2965       dz[ih][it]=-100;
2966       dy[ih][it]=-100;
2967       if (!cl[ih][it]) continue;
2968       Double_t  xcluster = cl[ih][it]->GetX();
2969       Double_t ytrack,ztrack;
2970       track2.GetProlongation(xcluster, ytrack, ztrack );
2971       dz[ih][it]  = cl[ih][it]->GetZ()- ztrack;                               // calculate distance from track  in z
2972       dy[ih][it]  = cl[ih][it]->GetY()+ dz[ih][it]*h01  -ytrack;     //                                in y
2973     }
2974     // minimize changes
2975     if (!cl[0][it]) continue;
2976     if (TMath::Abs(cl[0][it]->GetZ()-zmean[it])> padlength*0.8 &&cl[1][it])
2977       if (TMath::Abs(cl[1][it]->GetZ()-zmean[it])< padlength*0.5){
2978         best[0][it]=1;
2979       }
2980   }
2981   //
2982   // iterative choosing of "best path"
2983   //
2984   //
2985   Int_t label = TMath::Abs(track->GetLabel());
2986   Int_t bestiter=0;
2987   //
2988   for (Int_t iter=0;iter<9;iter++){
2989     //
2990     changes[iter]= 0;
2991     sumz      = 0; sum=0; sumdy=0;sumdy2=0;sumx=0;sumx2=0;sumxy=0;mpads=0; ngood[iter]=0; nbad[iter]=0;
2992     // linear fit
2993     for (Int_t it=0;it<t1-t0;it++){
2994       if (!cl[best[iter][it]][it]) continue;
2995       //calculates pad-row changes
2996       Double_t zbefore= cl[best[iter][it]][it]->GetZ();
2997       Double_t zafter = cl[best[iter][it]][it]->GetZ();
2998       for (Int_t itd = it-1; itd>=0;itd--) {
2999         if (cl[best[iter][itd]][itd]) {
3000           zbefore= cl[best[iter][itd]][itd]->GetZ();
3001           break;
3002         }
3003       }
3004       for (Int_t itd = it+1; itd<t1-t0;itd++) {
3005         if (cl[best[iter][itd]][itd]) {
3006           zafter= cl[best[iter][itd]][itd]->GetZ();
3007           break;
3008         }
3009       }
3010       if (TMath::Abs(cl[best[iter][it]][it]->GetZ()-zbefore)>0.1&&TMath::Abs(cl[best[iter][it]][it]->GetZ()-zafter)>0.1) changes[iter]++;
3011       //
3012       Double_t dx = x[it]-xmean;  // distance to reference x
3013       sumz += cl[best[iter][it]][it]->GetZ();
3014       sum++;
3015       sumdy += dy[best[iter][it]][it];
3016       sumdy2+= dy[best[iter][it]][it]*dy[best[iter][it]][it];
3017       sumx  += dx;
3018       sumx2 += dx*dx;
3019       sumxy  += dx*dy[best[iter][it]][it];
3020       mpads += cl[best[iter][it]][it]->GetNPads();
3021       if (cl[best[iter][it]][it]->GetLabel(0)==label || cl[best[iter][it]][it]->GetLabel(1)==label||cl[best[iter][it]][it]->GetLabel(2)==label){
3022         ngood[iter]++;
3023       }
3024       else{
3025         nbad[iter]++;
3026       }
3027     }
3028     //
3029     // calculates line parameters
3030     //
3031     Double_t det  = sum*sumx2-sumx*sumx;
3032     angle[iter]   = (sum*sumxy-sumx*sumdy)/det;
3033     mean[iter]    = (sumx2*sumdy-sumx*sumxy)/det;
3034     meanz[iter]   = sumz/sum;
3035     moffset[iter] = sumdy/sum;
3036     mpads        /= sum;                         // mean number of pads
3037     //
3038     //
3039     Double_t  sigma2 = 0;   // normalized residuals - for line fit
3040     Double_t  sigma1 = 0;   // normalized residuals - constant fit
3041     //
3042     for (Int_t it=0;it<t1-t0;it++){
3043       if (!cl[best[iter][it]][it]) continue;
3044       Double_t dx = x[it]-xmean;
3045       Double_t ytr = mean[iter]+angle[iter]*dx;
3046       sigma2 += (dy[best[iter][it]][it]-ytr)*(dy[best[iter][it]][it]-ytr);
3047       sigma1 +=  (dy[best[iter][it]][it]-moffset[iter])*(dy[best[iter][it]][it]-moffset[iter]);
3048       sum++;
3049     }
3050     sigma2      /=(sum-2);                    // normalized residuals
3051     sigma1      /=(sum-1);                    // normalized residuals
3052     //
3053     smean[iter]       = sigma2*(sumx2/det);   // estimated error2 of mean
3054     sangle[iter]      = sigma2*(sum/det);     // estimated error2 of angle
3055     smeanangle[iter]  = sigma2*(-sumx/det);   // correlation
3056     //
3057     //
3058     sigmas[iter]  = TMath::Sqrt(sigma1);      //
3059     smoffset[iter]= (sigma1/sum)+0.01*0.01;             // sigma of mean offset + unisochronity sigma
3060     //
3061     // iterative choosing of "better path"
3062     //
3063     for (Int_t it=0;it<t1-t0;it++){
3064       if (!cl[best[iter][it]][it]) continue;
3065       //
3066       Double_t sigmatr2 = smoffset[iter]+0.5*tany*tany;             //add unisochronity + angular effect contribution
3067       Double_t sweight  = 1./sigmatr2+1./track->GetSigmaY2();
3068       Double_t weighty  = (moffset[iter]/sigmatr2)/sweight;         // weighted mean
3069       Double_t sigmacl  = TMath::Sqrt(sigma1*sigma1+track->GetSigmaY2());   //
3070       Double_t mindist=100000;
3071       Int_t ihbest=0;
3072       for (Int_t ih=0;ih<10;ih++){
3073         if (!cl[ih][it]) break;
3074         Double_t dist2 = (dy[ih][it]-weighty)/sigmacl;
3075         dist2*=dist2;    //chi2 distance
3076         if (dist2<mindist){
3077           mindist = dist2;
3078           ihbest =ih;
3079         }
3080       }
3081       best[iter+1][it]=ihbest;
3082     }
3083     //
3084     //  update best hypothesy if better chi2 according tracklet position and angle
3085     //
3086     Double_t sy2 = smean[iter]  + track->GetSigmaY2();
3087     Double_t sa2 = sangle[iter] + track->fCee;
3088     Double_t say = track->fCey;
3089     //    Double_t chi20 = mean[bestiter]*mean[bestiter]/sy2+angle[bestiter]*angle[bestiter]/sa2;
3090     // Double_t chi21 = mean[iter]*mean[iter]/sy2+angle[iter]*angle[iter]/sa2;
3091
3092     Double_t detchi    = sy2*sa2-say*say;
3093     Double_t invers[3] = {sa2/detchi, sy2/detchi, -say/detchi};   //inverse value of covariance matrix
3094
3095     Double_t chi20 = mean[bestiter]*mean[bestiter]*invers[0]+angle[bestiter]*angle[bestiter]*invers[1]+
3096       2.*mean[bestiter]*angle[bestiter]*invers[2];
3097     Double_t chi21 = mean[iter]*mean[iter]*invers[0]+angle[iter]*angle[iter]*invers[1]+
3098       2*mean[iter]*angle[iter]*invers[2];
3099     tchi2s[iter] =chi21;
3100     //
3101     if (changes[iter]<=changes[bestiter] && chi21<chi20) {
3102       bestiter =iter;
3103     }
3104   }
3105   //
3106   //set clusters
3107   //
3108   Double_t sigma2 = sigmas[0];   // choose as sigma  from 0 iteration
3109   Short_t maxpos    = -1;
3110   Float_t maxcharge =  0;
3111   Short_t maxpos4    = -1;
3112   Float_t maxcharge4 =  0;
3113   Short_t maxpos5    = -1;
3114   Float_t maxcharge5 =  0;
3115
3116   //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
3117   //if (tchi2s[bestiter]>25.) sigma2=1000.;  // dont'accept
3118
3119   Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(AliTRDcalibDB::Instance()->GetVdrift(0,0,0));
3120   Double_t expectederr = sigma2*sigma2+0.01*0.01;
3121   if (mpads>3.5) expectederr  +=   (mpads-3.5)*0.04;
3122   if (changes[bestiter]>1) expectederr+=   changes[bestiter]*0.01;
3123   expectederr+=(0.03*(tany-exB)*(tany-exB))*15;
3124   //  if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
3125   //expectederr+=10000;
3126   for (Int_t it=0;it<t1-t0;it++){
3127     if (!cl[best[bestiter][it]][it]) continue;
3128     cl[best[bestiter][it]][it]->SetSigmaY2(expectederr);  // set cluster error
3129     if (!cl[best[bestiter][it]][it]->IsUsed()){
3130       cl[best[bestiter][it]][it]->SetY( cl[best[bestiter][it]][it]->GetY());
3131       //      cl[best[bestiter][it]][it]->Use();
3132     }
3133     //
3134     //  time bins with maximal charge
3135     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge){
3136       maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3137       maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3138     }
3139
3140     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge4){
3141       if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=4){
3142         maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3143         maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3144       }
3145     }
3146     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge5){
3147       if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=5){
3148         maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3149         maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3150       }
3151     }
3152     //
3153     //  time bins with maximal charge
3154     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge){
3155       maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3156       maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3157     }
3158
3159     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge4){
3160       if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=4){
3161         maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3162         maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3163       }
3164     }
3165     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge5){
3166       if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=5){
3167         maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
3168         maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
3169       }
3170     }
3171     clusters[it+t0] = indexes[best[bestiter][it]][it];
3172     //if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>4 && cl[best[bestiter][it]][it]->GetLocalTimeBin()<18) clusters[it+t0] = indexes[best[bestiter][it]][it];    //Test
3173   }
3174   //
3175   // set tracklet parameters
3176   //
3177   Double_t trackleterr2 = smoffset[bestiter]+0.01*0.01;
3178   if (mpads>3.5) trackleterr2  +=   (mpads-3.5)*0.04;
3179   trackleterr2+=   changes[bestiter]*0.01;
3180   trackleterr2*=   TMath::Max(14.-nfound,1.);
3181   trackleterr2+=   0.2*(tany-exB)*(tany-exB);
3182   //
3183   tracklet.Set(xmean, track2.GetY()+moffset[bestiter], meanz[bestiter], track2.GetAlpha(), trackleterr2);  //set tracklet parameters
3184   tracklet.SetTilt(h01);
3185   tracklet.SetP0(mean[bestiter]);
3186   tracklet.SetP1(angle[bestiter]);
3187   tracklet.SetN(nfound);
3188   tracklet.SetNCross(changes[bestiter]);
3189   tracklet.SetPlane(plane);
3190   tracklet.SetSigma2(expectederr);
3191   tracklet.SetChi2(tchi2s[bestiter]);
3192   tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
3193   track->fTracklets[plane] = tracklet;
3194   track->fNWrong+=nbad[0];
3195   //
3196   // Debuging part
3197   //
3198   TClonesArray array0("AliTRDcluster");
3199   TClonesArray array1("AliTRDcluster");
3200   array0.ExpandCreateFast(t1-t0+1);
3201   array1.ExpandCreateFast(t1-t0+1);
3202   TTreeSRedirector& cstream = *fDebugStreamer;
3203   AliTRDcluster dummy;
3204   Double_t dy0[100];
3205   Double_t dyb[100];
3206
3207   for (Int_t it=0;it<t1-t0;it++){
3208     dy0[it] = dy[0][it];
3209     dyb[it] = dy[best[bestiter][it]][it];
3210     if(cl[0][it]) {
3211       new(array0[it]) AliTRDcluster(*cl[0][it]);
3212     }
3213     else{
3214       new(array0[it]) AliTRDcluster(dummy);
3215     }
3216     if(cl[best[bestiter][it]][it]) {
3217       new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
3218     }
3219     else{
3220       new(array1[it]) AliTRDcluster(dummy);
3221     }
3222   }
3223   TGraph graph0(t1-t0,x,dy0);
3224   TGraph graph1(t1-t0,x,dyb);
3225   TGraph graphy(t1-t0,x,yt);
3226   TGraph graphz(t1-t0,x,zt);
3227   //
3228   //
3229   if (AliTRDReconstructor::StreamLevel()>0)
3230   cstream<<"tracklet"<<
3231     "track.="<<track<<                                       // track parameters
3232     "tany="<<tany<<                                          // tangent of the local track angle
3233     "xmean="<<xmean<<                                        // xmean - reference x of tracklet
3234     "tilt="<<h01<<                                           // tilt angle
3235     "nall="<<nall<<                                          // number of foundable clusters
3236     "nfound="<<nfound<<                                      // number of found clusters
3237     "clfound="<<clfound<<                                    // total number of found clusters in road
3238     "mpads="<<mpads<<                                        // mean number of pads per cluster
3239     "plane="<<plane<<                                        // plane number
3240     "detector="<<detector<<                                  // detector number
3241     "road="<<road<<                                          // the width of the used road
3242     "graph0.="<<&graph0<<                                    // x - y = dy for closest cluster
3243     "graph1.="<<&graph1<<                                    // x - y = dy for second closest cluster
3244     "graphy.="<<&graphy<<                                    // y position of the track
3245     "graphz.="<<&graphz<<                                    // z position of the track
3246     //    "fCl.="<<&array0<<                                       // closest cluster
3247     //"fCl2.="<<&array1<<                                      // second closest cluster
3248     "maxpos="<<maxpos<<                                      // maximal charge postion
3249     "maxcharge="<<maxcharge<<                                // maximal charge
3250     "maxpos4="<<maxpos4<<                                    // maximal charge postion - after bin 4
3251     "maxcharge4="<<maxcharge4<<                              // maximal charge         - after bin 4
3252     "maxpos5="<<maxpos5<<                                    // maximal charge postion - after bin 5
3253     "maxcharge5="<<maxcharge5<<                              // maximal charge         - after bin 5
3254     //
3255     "bestiter="<<bestiter<<                                  // best iteration number
3256     "tracklet.="<<&tracklet<<                                // corrspond to the best iteration
3257     "tchi20="<<tchi2s[0]<<                                   // chi2 of cluster in the 0 iteration
3258     "tchi2b="<<tchi2s[bestiter]<<                            // chi2 of cluster in the best  iteration
3259     "sigmas0="<<sigmas[0]<<                                  // residuals sigma
3260     "sigmasb="<<sigmas[bestiter]<<                           // residulas sigma
3261     //
3262     "ngood0="<<ngood[0]<<                                    // number of good clusters in 0 iteration
3263     "nbad0="<<nbad[0]<<                                      // number of bad clusters in 0 iteration
3264     "ngoodb="<<ngood[bestiter]<<                             //                        in  best iteration
3265     "nbadb="<<nbad[bestiter]<<                               //                        in  best iteration
3266     //
3267     "changes0="<<changes[0]<<                                // changes of pardrows in iteration number 0
3268     "changesb="<<changes[bestiter]<<                         // changes of pardrows in best iteration
3269     //
3270     "moffset0="<<moffset[0]<<                                // offset fixing angle in iter=0
3271     "smoffset0="<<smoffset[0]<<                              // sigma of offset fixing angle in iter=0
3272     "moffsetb="<<moffset[bestiter]<<                         // offset fixing angle in iter=best
3273     "smoffsetb="<<smoffset[bestiter]<<                       // sigma of offset fixing angle in iter=best
3274     //
3275     "mean0="<<mean[0]<<                                      // mean dy in iter=0;
3276     "smean0="<<smean[0]<<                                    // sigma of mean dy in iter=0
3277     "meanb="<<mean[bestiter]<<                               // mean dy in iter=best
3278     "smeanb="<<smean[bestiter]<<                             // sigma of mean dy in iter=best
3279     //
3280     "angle0="<<angle[0]<<                                    // angle deviation in the iteration number 0
3281     "sangle0="<<sangle[0]<<                                  // sigma of angular deviation in iteration number 0
3282     "angleb="<<angle[bestiter]<<                             // angle deviation in the best iteration
3283     "sangleb="<<sangle[bestiter]<<                           // sigma of angle deviation in the best iteration
3284     //
3285     "expectederr="<<expectederr<<                            // expected error of cluster position
3286     "\n";
3287   //
3288   //
3289   return nfound;
3290 }
3291
3292
3293 Int_t  AliTRDtracker::Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down)
3294 {
3295   //
3296   //  Sort eleements according occurancy
3297   //  The size of output array has is 2*n
3298   //
3299   Int_t * sindexS = new Int_t[n];     // temp array for sorting
3300   Int_t * sindexF = new Int_t[2*n];
3301   for (Int_t i=0;i<n;i++) sindexF[i]=0;
3302   //
3303   TMath::Sort(n,inlist, sindexS, down);
3304   Int_t last      = inlist[sindexS[0]];
3305   Int_t val       = last;
3306   sindexF[0]      = 1;
3307   sindexF[0+n]    = last;
3308   Int_t countPos  = 0;
3309   //
3310   //  find frequency
3311   for(Int_t i=1;i<n; i++){
3312     val = inlist[sindexS[i]];
3313     if (last == val)   sindexF[countPos]++;
3314     else{
3315       countPos++;
3316       sindexF[countPos+n] = val;
3317       sindexF[countPos]++;
3318       last =val;
3319     }
3320   }
3321   if (last==val) countPos++;
3322   // sort according frequency
3323   TMath::Sort(countPos, sindexF, sindexS, kTRUE);
3324   for (Int_t i=0;i<countPos;i++){
3325     outlist[2*i  ] = sindexF[sindexS[i]+n];
3326     outlist[2*i+1] = sindexF[sindexS[i]];
3327   }
3328   delete [] sindexS;
3329   delete [] sindexF;
3330
3331   return countPos;
3332 }
3333
3334 AliTRDtrack * AliTRDtracker::RegisterSeed(AliTRDseed * seeds, Double_t * params)
3335 {
3336   //
3337   //
3338   //
3339   Double_t alpha=AliTRDgeometry::GetAlpha();
3340   Double_t shift=AliTRDgeometry::GetAlpha()/2.;
3341   Double_t c[15];
3342   c[0] = 0.2;
3343   c[1] = 0  ; c[2] = 2;
3344   c[3] = 0  ; c[4] = 0; c[5] = 0.02;
3345   c[6] = 0  ; c[7] = 0; c[8] = 0;      c[9] = 0.1;
3346   c[10] = 0  ; c[11] = 0; c[12] = 0;   c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
3347   //
3348   Int_t index =0;
3349   AliTRDcluster *cl =0;
3350   for (Int_t ilayer=0;ilayer<6;ilayer++){
3351     if (seeds[ilayer].IsOK()){
3352       for (Int_t itime=22;itime>0;itime--){
3353         if (seeds[ilayer].fIndexes[itime]>0){
3354           index = seeds[ilayer].fIndexes[itime];
3355           cl = seeds[ilayer].fClusters[itime];
3356           break;
3357         }
3358       }
3359     }
3360     if (index>0) break;
3361   }
3362   if (cl==0) return 0;
3363   AliTRDtrack * track  = new AliTRDtrack(cl,index,&params[1],c, params[0],params[6]*alpha+shift);
3364   track->PropagateTo(params[0]-5.);
3365   track->ResetCovariance(1);
3366   //
3367   Int_t rc=FollowBackProlongation(*track);
3368   if (rc<30) {
3369     delete track;
3370     track =0;
3371   }else{
3372     track->CookdEdx();
3373     CookdEdxTimBin(*track);
3374     CookLabel(track, 0.9);
3375   }
3376   return track;
3377 }
3378
3379
3380
3381
3382
3383
3384 AliTRDseed::AliTRDseed()
3385 {
3386   //
3387   //
3388   fTilt =0;         // tilting angle
3389   fPadLength = 0;   // pad length
3390   fX0 = 0;           // x0 position
3391   for (Int_t i=0;i<25;i++){
3392     fX[i]=0;        // !x position
3393     fY[i]=0;        // !y position
3394     fZ[i]=0;        // !z position
3395     fIndexes[i]=0;  // !indexes
3396     fClusters[i]=0; // !clusters
3397   }
3398   for (Int_t i=0;i<2;i++){
3399     fYref[i]=0;      // reference y
3400     fZref[i]=0;      // reference z
3401     fYfit[i]=0;      // y fit position +derivation
3402     fYfitR[i]=0;      // y fit position +derivation
3403     fZfit[i]=0;      // z fit position
3404     fZfitR[i]=0;      // z fit position
3405     fLabels[i]=0;    // labels
3406   }
3407   fSigmaY  = 0;
3408   fSigmaY2 = 0;
3409   fMeanz=0;         // mean vaue of z
3410   fZProb=0;         // max probbable z
3411   fMPads=0;
3412   //
3413   fN=0;            // number of associated clusters
3414   fN2=0;            // number of not crossed
3415   fNUsed=0;        // number of used clusters
3416   fNChange=0;      // change z counter
3417 }
3418
3419 void AliTRDseed::Reset(){
3420   //
3421   // reset seed
3422   //
3423   for (Int_t i=0;i<25;i++){
3424     fX[i]=0;        // !x position
3425     fY[i]=0;        // !y position
3426     fZ[i]=0;        // !z position
3427     fIndexes[i]=0;  // !indexes
3428     fClusters[i]=0; // !clusters
3429     fUsable[i]  = kFALSE;
3430   }
3431   for (Int_t i=0;i<2;i++){
3432     fYref[i]=0;      // reference y
3433     fZref[i]=0;      // reference z
3434     fYfit[i]=0;      // y fit position +derivation
3435     fYfitR[i]=0;      // y fit position +derivation
3436     fZfit[i]=0;      // z fit position
3437     fZfitR[i]=0;      // z fit position
3438     fLabels[i]=-1;    // labels
3439   }
3440   fSigmaY =0;         //"robust" sigma in y
3441   fSigmaY2=0;         //"robust" sigma in y
3442   fMeanz =0;         // mean vaue of z
3443   fZProb =0;         // max probbable z
3444   fMPads =0;
3445   //
3446   fN=0;            // number of associated clusters
3447   fN2=0;            // number of not crossed
3448   fNUsed=0;        // number of used clusters
3449   fNChange=0;      // change z counter
3450 }
3451
3452 void AliTRDseed::CookLabels(){
3453   //
3454   // cook 2 labels for seed
3455   //
3456   Int_t labels[200];
3457   Int_t out[200];
3458   Int_t nlab =0;
3459   for (Int_t i=0;i<25;i++){
3460     if (!fClusters[i]) continue;
3461     for (Int_t ilab=0;ilab<3;ilab++){
3462       if (fClusters[i]->GetLabel(ilab)>=0){
3463         labels[nlab] = fClusters[i]->GetLabel(ilab);
3464         nlab++;
3465       }
3466     }
3467   }
3468   Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
3469   fLabels[0] = out[0];
3470   if (nlab2>1 && out[3]>1) fLabels[1] =out[2];
3471 }
3472
3473 void   AliTRDseed::UseClusters()
3474 {
3475   //
3476   // use clusters
3477   //
3478    for (Int_t i=0;i<25;i++){
3479      if (!fClusters[i]) continue;
3480      if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
3481    }
3482 }
3483
3484
3485 void        AliTRDseed::Update(){
3486   //
3487   //
3488   //
3489   const Float_t kRatio = 0.8;
3490   const Int_t   kClmin        = 6;
3491   const Float_t kmaxtan  = 2;
3492   if (TMath::Abs(fYref[1])>kmaxtan) return;             // too much inclined track
3493   //
3494   Float_t  sigmaexp = 0.05+TMath::Abs(fYref[1]*0.25);   // expected r.m.s in y direction
3495   Float_t  ycrosscor = fPadLength*fTilt*0.5;             // y correction for crossing
3496   fNChange =0;
3497   //
3498   Double_t sumw, sumwx,sumwx2;
3499   Double_t sumwy, sumwxy, sumwz,sumwxz;
3500   Int_t    zints[25];        // histograming of the z coordinate - get 1 and second max probable coodinates in z
3501   Int_t    zouts[50];        //
3502   Float_t  allowedz[25];     // allowed z for given time bin
3503   Float_t  yres[25];         // residuals from reference
3504   Float_t  anglecor = fTilt*fZref[1];  //correction to the angle
3505   //
3506   //
3507   fN=0; fN2 =0;
3508   for (Int_t i=0;i<25;i++){
3509     yres[i] =10000;
3510     if (!fClusters[i]) continue;
3511     yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i];   // residual y
3512     zints[fN] = Int_t(fZ[i]);
3513     fN++;
3514   }
3515   if (fN<kClmin) return;
3516   Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
3517   fZProb   = zouts[0];
3518   if (nz<=1) zouts[3]=0;
3519   if (zouts[1]+zouts[3]<kClmin) return;
3520   //
3521   if (TMath::Abs(zouts[0]-zouts[2])>12.) zouts[3]=0;   // z distance bigger than pad - length
3522   //
3523   Int_t  breaktime = -1;
3524   Bool_t mbefore   = kFALSE;
3525   Int_t  cumul[25][2];
3526   Int_t  counts[2]={0,0};
3527   //
3528   if (zouts[3]>=3){
3529     //
3530     // find the break time allowing one chage on pad-rows with maximal numebr of accepted clusters
3531     //
3532     fNChange=1;
3533     for (Int_t i=0;i<25;i++){
3534       cumul[i][0] = counts[0];
3535       cumul[i][1] = counts[1];
3536       if (TMath::Abs(fZ[i]-zouts[0])<2) counts[0]++;
3537       if (TMath::Abs(fZ[i]-zouts[2])<2) counts[1]++;
3538     }
3539     Int_t  maxcount  = 0;
3540     for (Int_t i=0;i<24;i++) {
3541       Int_t after  = cumul[24][0]-cumul[i][0];
3542       Int_t before = cumul[i][1];
3543       if (after+before>maxcount) {
3544         maxcount=after+before;
3545         breaktime=i;
3546         mbefore=kFALSE;
3547       }
3548       after  = cumul[24][1]-cumul[i][1];
3549       before = cumul[i][0];
3550       if (after+before>maxcount) {
3551         maxcount=after+before;
3552         breaktime=i;
3553         mbefore=kTRUE;
3554       }
3555     }
3556     breaktime-=1;
3557   }
3558   for (Int_t i=0;i<25;i++){
3559     if (i>breaktime)  allowedz[i] =   mbefore  ? zouts[2]:zouts[0];
3560     if (i<=breaktime) allowedz[i] = (!mbefore) ? zouts[2]:zouts[0];
3561   }
3562   if ( (allowedz[0]>allowedz[24] && fZref[1]<0) || (allowedz[0]<allowedz[24] &&  fZref[1]>0)){
3563     //
3564     // tracklet z-direction not in correspondance with track z direction
3565     //
3566     fNChange =0;
3567     for (Int_t i=0;i<25;i++){
3568       allowedz[i] =  zouts[0];  //only longest taken
3569     }
3570   }
3571   //
3572   if (fNChange>0){
3573     //
3574     // cross pad -row tracklet  - take the step change into account
3575     //
3576     for (Int_t i=0;i<25;i++){
3577       if (!fClusters[i]) continue;
3578       if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
3579       yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i];   // residual y
3580       if (TMath::Abs(fZ[i]-fZProb)>2){
3581         if (fZ[i]>fZProb) yres[i]+=fTilt*fPadLength;
3582         if (fZ[i]<fZProb) yres[i]-=fTilt*fPadLength;
3583       }
3584     }
3585   }
3586   //
3587   Double_t yres2[25];
3588   Double_t mean,sigma;
3589   for (Int_t i=0;i<25;i++){
3590     if (!fClusters[i]) continue;
3591     if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
3592     yres2[fN2] =  yres[i];
3593     fN2++;
3594   }
3595   if (fN2<kClmin){
3596     fN2 = 0;
3597     return;
3598   }
3599   EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*kRatio-2));
3600   if (sigma<sigmaexp*0.8) sigma=sigmaexp;
3601   fSigmaY = sigma;
3602   //
3603   //
3604   // reset sums
3605   sumw=0; sumwx=0; sumwx2=0;
3606   sumwy=0; sumwxy=0; sumwz=0;sumwxz=0;
3607   fN2 =0;
3608   fMeanz =0;
3609   fMPads =0;
3610   //
3611   for (Int_t i=0;i<25;i++){
3612     fUsable[i]=kFALSE;
3613     if (!fClusters[i]) continue;
3614     if (TMath::Abs(fZ[i]-allowedz[i])>2)  continue;
3615     if (TMath::Abs(yres[i]-mean)>4.*sigma) continue;
3616     fUsable[i] = kTRUE;
3617     fN2++;
3618     fMPads+=fClusters[i]->GetNPads();
3619     Float_t weight =1;
3620     if (fClusters[i]->GetNPads()>4) weight=0.5;
3621     if (fClusters[i]->GetNPads()>5) weight=0.2;
3622     //
3623     Double_t x = fX[i];
3624     sumw+=weight; sumwx+=x*weight; sumwx2+=x*x*weight;
3625     sumwy+=weight*yres[i];  sumwxy+=weight*(yres[i])*x;
3626     sumwz+=weight*fZ[i];    sumwxz+=weight*fZ[i]*x;
3627   }
3628   if (fN2<kClmin){
3629     fN2 = 0;
3630     return;
3631   }
3632   fMeanz       = sumwz/sumw;
3633   Float_t correction =0;
3634   if (fNChange>0){
3635     // tracklet on boundary
3636     if (fMeanz<fZProb) correction =   ycrosscor;
3637     if (fMeanz>fZProb) correction =  -ycrosscor;
3638   }
3639   Double_t det = sumw*sumwx2-sumwx*sumwx;
3640   fYfitR[0]    = (sumwx2*sumwy-sumwx*sumwxy)/det;
3641   fYfitR[1]    = (sumw*sumwxy-sumwx*sumwy)/det;
3642   //
3643   fSigmaY2     =0;
3644   for (Int_t i=0;i<25;i++){
3645     if (!fUsable[i]) continue;
3646     Float_t delta = yres[i]-fYfitR[0]-fYfitR[1]*fX[i];
3647     fSigmaY2+=delta*delta;
3648   }
3649   fSigmaY2 = TMath::Sqrt(fSigmaY2/Float_t(fN2-2));
3650   //
3651   fZfitR[0]    = (sumwx2*sumwz-sumwx*sumwxz)/det;
3652   fZfitR[1]    = (sumw*sumwxz-sumwx*sumwz)/det;
3653   fZfit[0]     = (sumwx2*sumwz-sumwx*sumwxz)/det;
3654   fZfit[1]     = (sumw*sumwxz-sumwx*sumwz)/det;
3655   fYfitR[0]   += fYref[0]+correction;
3656   fYfitR[1]   += fYref[1];
3657   fYfit[0]     = fYfitR[0];
3658   fYfit[1]     = fYfitR[1];
3659   //
3660   //
3661   UpdateUsed();
3662 }
3663
3664
3665
3666
3667
3668
3669 void AliTRDseed::UpdateUsed(){
3670   //
3671   fNUsed =0;
3672   for (Int_t i=0;i<25;i++){
3673      if (!fClusters[i]) continue;
3674      if ((fClusters[i]->IsUsed())) fNUsed++;
3675   }
3676 }
3677
3678
3679 void AliTRDseed::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh)
3680 {
3681   //
3682   // robust estimator in 1D case MI version
3683   //
3684   //for the univariate case
3685   //estimates of location and scatter are returned in mean and sigma parameters
3686   //the algorithm works on the same principle as in multivariate case -
3687   //it finds a subset of size hh with smallest sigma, and then returns mean and
3688   //sigma of this subset
3689
3690   if (hh==0)
3691     hh=(nvectors+2)/2;
3692   Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
3693   Int_t *index=new Int_t[nvectors];
3694   TMath::Sort(nvectors, data, index, kFALSE);
3695   //
3696   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
3697   Double_t factor = faclts[nquant-1];
3698   //
3699   //
3700   Double_t sumx  =0;
3701   Double_t sumx2 =0;
3702   Int_t    bestindex = -1;
3703   Double_t bestmean  = 0;
3704   Double_t bestsigma = data[index[nvectors-1]]-data[index[0]];   // maximal possible sigma
3705   for (Int_t i=0; i<hh; i++){
3706     sumx  += data[index[i]];
3707     sumx2 += data[index[i]]*data[index[i]];
3708   }
3709   //
3710   Double_t norm = 1./Double_t(hh);
3711   Double_t norm2 = 1./Double_t(hh-1);
3712   for (Int_t i=hh; i<nvectors; i++){
3713     Double_t cmean  = sumx*norm;
3714     Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
3715     if (csigma<bestsigma){
3716       bestmean  = cmean;
3717       bestsigma = csigma;
3718       bestindex = i-hh;
3719     }
3720     //
3721     //
3722     sumx  += data[index[i]]-data[index[i-hh]];
3723     sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
3724   }
3725
3726   Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
3727   mean  = bestmean;
3728   sigma = bstd;
3729   delete [] index;
3730 }
3731
3732
3733 Float_t   AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror){
3734   //
3735   //
3736   //
3737   TLinearFitter fitterT2(4,"hyp4");  // fitting with tilting pads - kz not fixed
3738   fitterT2.StoreData(kTRUE);
3739   Float_t xref2 = (cseed[2].fX0+cseed[3].fX0)*0.5; // reference x0 for z
3740   //
3741   Int_t npointsT =0;
3742   fitterT2.ClearPoints();
3743   for (Int_t iLayer=0; iLayer<6;iLayer++){
3744     if (!cseed[iLayer].IsOK()) continue;
3745     Double_t tilt = cseed[iLayer].fTilt;
3746
3747     for (Int_t itime=0;itime<25;itime++){
3748       if (!cseed[iLayer].fUsable[itime]) continue;
3749       Double_t x   = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2;  // x relative to the midle chamber
3750       Double_t y   = cseed[iLayer].fY[itime];
3751       Double_t z   = cseed[iLayer].fZ[itime];
3752       // tilted rieman
3753       //
3754       Double_t uvt[6];
3755       Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0;      // global x
3756       Double_t t = 1./(x2*x2+y*y);
3757       uvt[1]  = t;    // t
3758       uvt[0]  = 2.*x2*uvt[1];      // u
3759       uvt[2]  = 2.0*tilt*uvt[1];
3760       uvt[3]  = 2.0*tilt*x*uvt[1];
3761       uvt[4]  = 2.0*(y+tilt*z)*uvt[1];
3762       //
3763       Double_t error = 2*uvt[1];
3764       if (terror) error*=cseed[iLayer].fSigmaY;
3765       else {error *=0.2;} //default error
3766       fitterT2.AddPoint(uvt,uvt[4],error);
3767       npointsT++;
3768     }
3769   }
3770   fitterT2.Eval();
3771   Double_t rpolz0 = fitterT2.GetParameter(3);
3772   Double_t rpolz1 = fitterT2.GetParameter(4);
3773   //
3774   // linear fitter  - not possible to make boundaries
3775   // non accept non possible z and dzdx combination
3776   //
3777   Bool_t   acceptablez =kTRUE;
3778   for (Int_t iLayer=0; iLayer<6;iLayer++){
3779     if (cseed[iLayer].IsOK()){
3780       Double_t zT2 =  rpolz0+rpolz1*(cseed[iLayer].fX0 - xref2);
3781       if (TMath::Abs(cseed[iLayer].fZProb-zT2)>cseed[iLayer].fPadLength*0.5+1)
3782         acceptablez = kFALSE;
3783     }
3784   }
3785   if (!acceptablez){
3786     Double_t zmf  = cseed[2].fZref[0]+cseed[2].fZref[1]*(xref2-cseed[2].fX0);
3787     Double_t dzmf = (cseed[2].fZref[1]+ cseed[3].fZref[1])*0.5;
3788     fitterT2.FixParameter(3,zmf);
3789     fitterT2.FixParameter(4,dzmf);
3790     fitterT2.Eval();
3791     fitterT2.ReleaseParameter(3);
3792     fitterT2.ReleaseParameter(4);
3793     rpolz0 = fitterT2.GetParameter(3);
3794     rpolz1 = fitterT2.GetParameter(4);
3795   }
3796   //
3797   Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
3798   Double_t params[3];
3799   params[0]     =  fitterT2.GetParameter(0);
3800   params[1]     =  fitterT2.GetParameter(1);
3801   params[2]     =  fitterT2.GetParameter(2);
3802   Double_t curvature     =  1+params[1]*params[1]-params[2]*params[0];
3803   for (Int_t iLayer = 0; iLayer<6;iLayer++){
3804     Double_t  x = cseed[iLayer].fX0;
3805     Double_t  y=0,dy=0, z=0, dz=0;
3806     // y
3807     Double_t res2 = (x*params[0]+params[1]);
3808     res2*=res2;
3809     res2 = 1.-params[2]*params[0]+params[1]*params[1]-res2;
3810     if (res2>=0){
3811       res2 = TMath::Sqrt(res2);
3812       y    = (1-res2)/params[0];
3813     }
3814     //dy
3815     Double_t x0 = -params[1]/params[0];
3816     if (-params[2]*params[0]+params[1]*params[1]+1>0){
3817       Double_t rm1  = params[0]/TMath::Sqrt(-params[2]*params[0]+params[1]*params[1]+1);
3818       if ( 1./(rm1*rm1)-(x-x0)*(x-x0)>0){
3819         Double_t res = (x-x0)/TMath::Sqrt(1./(rm1*rm1)-(x-x0)*(x-x0));
3820         if (params[0]<0) res*=-1.;
3821         dy = res;
3822       }
3823     }
3824     z  = rpolz0+rpolz1*(x-xref2);
3825     dz = rpolz1;
3826     cseed[iLayer].fYref[0] = y;
3827     cseed[iLayer].fYref[1] = dy;
3828     cseed[iLayer].fZref[0] = z;
3829     cseed[iLayer].fZref[1] = dz;
3830     cseed[iLayer].fC  = curvature;
3831     //
3832   }
3833   return chi2TR;
3834 }