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 //                                                                           //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <Riostream.h>
  25 #include <TFile.h>
  26 #include <TBranch.h>
  27 #include <TTree.h>
  28 #include <TObjArray.h>
  29
  30 #include "AliTRDgeometry.h"
  31 #include "AliTRDparameter.h"
  32 #include "AliTRDpadPlane.h"
  33 #include "AliTRDgeometryDetail.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
  49 //
  50
  51 ClassImp(AliTRDtracker)
  52 ClassImp(AliTRDseed)
  53
  54   const  Float_t     AliTRDtracker::fgkSeedDepth          = 0.5;
  55   const  Float_t     AliTRDtracker::fgkSeedStep           = 0.10;
  56   const  Float_t     AliTRDtracker::fgkSeedGap            = 0.25;
  57
  58  const  Float_t     AliTRDtracker::fgkMaxSeedDeltaZ12    = 40.;
  59   const  Float_t     AliTRDtracker::fgkMaxSeedDeltaZ      = 25.;
  60   const  Float_t     AliTRDtracker::fgkMaxSeedC           = 0.0052;
  61   const  Float_t     AliTRDtracker::fgkMaxSeedTan         = 1.2;
  62  const  Float_t     AliTRDtracker::fgkMaxSeedVertexZ     = 150.;
  63
  64   const  Double_t    AliTRDtracker::fgkSeedErrorSY        = 0.2;
  65   const  Double_t    AliTRDtracker::fgkSeedErrorSY3       = 2.5;
  66   const  Double_t    AliTRDtracker::fgkSeedErrorSZ        = 0.1;
  67
  68   const  Float_t     AliTRDtracker::fgkMinClustersInSeed  = 0.7;
  69
  70   const  Float_t     AliTRDtracker::fgkMinClustersInTrack = 0.5;
  71   const  Float_t     AliTRDtracker::fgkMinFractionOfFoundClusters = 0.8;
  72
  73   const  Float_t     AliTRDtracker::fgkSkipDepth          = 0.3;
  74   const  Float_t     AliTRDtracker::fgkLabelFraction      = 0.8;
  75   const  Float_t     AliTRDtracker::fgkWideRoad           = 20.;
  76
  77   const  Double_t    AliTRDtracker::fgkMaxChi2            = 12.;
  78
  79 //   const  Double_t    AliTRDtracker::fgkOffset             = -0.012;
  80 //   const  Double_t    AliTRDtracker::fgkOffsetX            = 0.35;
  81 //   const  Double_t    AliTRDtracker::fgkCoef               = 0.00;
  82 //   const  Double_t    AliTRDtracker::fgkMean               = 8.;
  83 //   const  Double_t    AliTRDtracker::fgkDriftCorrection    = 1.07;
  84 //   const  Double_t    AliTRDtracker::fgkExB                = 0.072;
  85
  86   const  Double_t    AliTRDtracker::fgkOffset             = -0.019;
  87   const  Double_t    AliTRDtracker::fgkOffsetX            = 0.26;       // "time offset"
  88 //  const  Double_t    AliTRDtracker::fgkCoef               = 0.0096;   // angular shift
  89   const  Double_t    AliTRDtracker::fgkCoef               = 0.0106;   // angular shift
  90   const  Double_t    AliTRDtracker::fgkMean               = 0.;
  91   const  Double_t    AliTRDtracker::fgkDriftCorrection    = 1.055;   // drift coefficient correction
  92   const  Double_t    AliTRDtracker::fgkExB                = 0.072;  // ExB angle - for error parameterization
  93
  94
  95 //     poscorrection =  fgkCoef*(GetLocalTimeBin() - fgkMean)+fgkOffset;
  96
  97 const Int_t AliTRDtracker::fgkFirstPlane = 5;
  98 const Int_t AliTRDtracker::fgkLastPlane = 17;
  99
 100 //____________________________________________________________________
 101 AliTRDtracker::AliTRDtracker():AliTracker(),
 102                                fGeom(0),
 103                                fPar(0),
 104                                fNclusters(0),
 105                                fClusters(0),
 106                                fNseeds(0),
 107                                fSeeds(0),
 108                                fNtracks(0),
 109                                fTracks(0),
 110                                fSY2corr(0),
 111                                fSZ2corr(0),
 112                                fTimeBinsPerPlane(0),
 113                                fMaxGap(0),
 114                                fVocal(kFALSE),
 115                                fAddTRDseeds(kFALSE),
 116                                fNoTilt(kFALSE)
 117 {
 118   // Default constructor
 119
 120   for(Int_t i=0;i<kTrackingSectors;i++) fTrSec[i]=0;
 121   for(Int_t j=0;j<5;j++)
 122     for(Int_t k=0;k<18;k++) fHoles[j][k]=kFALSE;
 123   fDebugStreamer = 0;
 124 }
 125 //____________________________________________________________________
 126 AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
 127 {
 128   //
 129   //  Main constructor
 130   //
 131
 132   //Float_t fTzero = 0;
 133
 134   fAddTRDseeds = kFALSE;
 135   fGeom = NULL;
 136   fNoTilt = kFALSE;
 137
 138   TDirectory *savedir=gDirectory;
 139   TFile *in=(TFile*)geomfile;
 140   if (!in->IsOpen()) {
 141     printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n");
 142     printf("    DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
 143   }
 144   else {
 145     in->cd();
 146 //    in->ls();
 147     fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
 148     fPar  = (AliTRDparameter*) in->Get("TRDparameter");
 149 //    fGeom->Dump();
 150   }
 151
 152   if(fGeom) {
 153     //    fTzero = geo->GetT0();
 154     //    printf("Found geometry version %d on file \n", fGeom->IsVersion());
 155   }
 156   else {
 157     printf("AliTRDtracker::AliTRDtracker(): can't find TRD geometry!\n");
 158     //printf("The DETAIL TRD geometry will be used\n");
 159     //fGeom = new AliTRDgeometryDetail();
 160     fGeom = new AliTRDgeometryDetail();
 161     fGeom->SetPHOShole();
 162     fGeom->SetRICHhole();
 163   }
 164
 165   if (!fPar) {
 166     printf("AliTRDtracker::AliTRDtracker(): can't find TRD parameter!\n");
 167     printf("The DEFAULT TRD parameter will be used\n");
 168     fPar = new AliTRDparameter("Pica","Vyjebana");
 169   }
 170   fPar = new AliTRDparameter("Pica","Vyjebana");
 171   fPar->Init();
 172
 173   savedir->cd();
 174
 175   //  fGeom->SetT0(fTzero);
 176
 177   fNclusters = 0;
 178   fClusters  = new TObjArray(2000);
 179   fNseeds    = 0;
 180   fSeeds     = new TObjArray(2000);
 181   fNtracks   = 0;
 182   fTracks    = new TObjArray(1000);
 183
 184   for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
 185     Int_t trS = CookSectorIndex(geomS);
 186     fTrSec[trS] = new AliTRDtrackingSector(fGeom, geomS, fPar);
 187     for (Int_t icham=0;icham<AliTRDgeometry::kNcham; icham++){
 188       fHoles[icham][trS]=fGeom->IsHole(0,icham,geomS);
 189     }
 190   }
 191   AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
 192   Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
 193   //  Float_t tiltAngle = TMath::Abs(fPar->GetTiltingAngle());
 194   if(tiltAngle < 0.1) {
 195     fNoTilt = kTRUE;
 196   }
 197
 198   fSY2corr = 0.2;
 199   fSZ2corr = 120.;
 200
 201   if(fNoTilt && (tiltAngle > 0.1)) fSY2corr = fSY2corr + tiltAngle * 0.05;
 202
 203
 204   // calculate max gap on track
 205
 206   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
 207   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
 208
 209   Double_t dx = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
 210                          / AliTRDcalibDB::Instance()->GetSamplingFrequency();
 211
 212   Int_t tbAmp = fPar->GetTimeBefore();
 213   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
 214   if(kTRUE) maxAmp = 0;  // intentional until we change the parameter class
 215   Int_t tbDrift = fPar->GetTimeMax();
 216   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx)+4; // MI change  - take also last time bins
 217
 218   tbDrift = TMath::Min(tbDrift,maxDrift);
 219   tbAmp = TMath::Min(tbAmp,maxAmp);
 220
 221   fTimeBinsPerPlane = tbAmp + tbDrift;
 222   fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fgkSkipDepth);
 223
 224   fVocal = kFALSE;
 225
 226   fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
 227
 228   savedir->cd();
 229 }
 230
 231 //___________________________________________________________________
 232 AliTRDtracker::~AliTRDtracker()
 233 {
 234   //
 235   // Destructor of AliTRDtracker
 236   //
 237
 238   if (fClusters) {
 239     fClusters->Delete();
 240     delete fClusters;
 241   }
 242   if (fTracks) {
 243     fTracks->Delete();
 244     delete fTracks;
 245   }
 246   if (fSeeds) {
 247     fSeeds->Delete();
 248     delete fSeeds;
 249   }
 250   delete fGeom;
 251   delete fPar;
 252
 253   for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
 254     delete fTrSec[geomS];
 255   }
 256   if (fDebugStreamer) {
 257     //fDebugStreamer->Close();
 258     delete fDebugStreamer;
 259   }
 260 }
 261
 262 //_____________________________________________________________________
 263
 264 Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) {
 265   //
 266   // Rotates the track when necessary
 267   //
 268
 269   Double_t alpha = AliTRDgeometry::GetAlpha();
 270   Double_t y = track->GetY();
 271   Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
 272
 273   //Int_t ns = AliTRDgeometry::kNsect;
 274   //Int_t s=Int_t(track->GetAlpha()/alpha)%ns;
 275
 276   if (y > ymax) {
 277     //s = (s+1) % ns;
 278     if (!track->Rotate(alpha)) return kFALSE;
 279   } else if (y <-ymax) {
 280     //s = (s-1+ns) % ns;
 281     if (!track->Rotate(-alpha)) return kFALSE;
 282   }
 283
 284   return kTRUE;
 285 }
 286
 287 //_____________________________________________________________________
 288 inline Double_t f1trd(Double_t x1,Double_t y1,
 289                       Double_t x2,Double_t y2,
 290                       Double_t x3,Double_t y3)
 291 {
 292   //
 293   // Initial approximation of the track curvature
 294   //
 295   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 296   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 297                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 298   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 299                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 300
 301   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 302
 303   return -xr*yr/sqrt(xr*xr+yr*yr);
 304 }
 305
 306 //_____________________________________________________________________
 307 inline Double_t f2trd(Double_t x1,Double_t y1,
 308                       Double_t x2,Double_t y2,
 309                       Double_t x3,Double_t y3)
 310 {
 311   //
 312   // Initial approximation of the track curvature times X coordinate
 313   // of the center of curvature
 314   //
 315
 316   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 317   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 318                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 319   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 320                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 321
 322   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 323
 324   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 325 }
 326
 327 //_____________________________________________________________________
 328 inline Double_t f3trd(Double_t x1,Double_t y1,
 329                       Double_t x2,Double_t y2,
 330                       Double_t z1,Double_t z2)
 331 {
 332   //
 333   // Initial approximation of the tangent of the track dip angle
 334   //
 335
 336   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 337 }
 338
 339
 340 AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin, UInt_t &index){
 341   //
 342   //try to find cluster in the backup list
 343   //
 344   AliTRDcluster * cl =0;
 345   UInt_t *indexes = track->GetBackupIndexes();
 346   for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
 347     if (indexes[i]==0) break;
 348     AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
 349     if (!cli) break;
 350     if (cli->GetLocalTimeBin()!=timebin) continue;
 351     Int_t iplane = fGeom->GetPlane(cli->GetDetector());
 352     if (iplane==plane) {
 353       cl = cli;
 354       index = indexes[i];
 355       break;
 356     }
 357   }
 358   return cl;
 359 }
 360
 361
 362 Int_t  AliTRDtracker::GetLastPlane(AliTRDtrack * track){
 363   //
 364   //return last updated plane
 365   Int_t lastplane=0;
 366   UInt_t *indexes = track->GetBackupIndexes();
 367   for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
 368     AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
 369     if (!cli) break;
 370     Int_t iplane = fGeom->GetPlane(cli->GetDetector());
 371     if (iplane>lastplane) {
 372       lastplane = iplane;
 373     }
 374   }
 375   return lastplane;
 376 }
 377 //___________________________________________________________________
 378 Int_t AliTRDtracker::Clusters2Tracks(AliESD* event)
 379 {
 380   //
 381   // Finds tracks within the TRD. The ESD event is expected to contain seeds
 382   // at the outer part of the TRD. The seeds
 383   // are found within the TRD if fAddTRDseeds is TRUE.
 384   // The tracks are propagated to the innermost time bin
 385   // of the TRD and the ESD event is updated
 386   //
 387
 388   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 389   Float_t foundMin = fgkMinClustersInTrack * timeBins;
 390   Int_t nseed = 0;
 391   Int_t found = 0;
 392   Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
 393
 394   Int_t n = event->GetNumberOfTracks();
 395   for (Int_t i=0; i<n; i++) {
 396     AliESDtrack* seed=event->GetTrack(i);
 397     ULong_t status=seed->GetStatus();
 398     if ( (status & AliESDtrack::kTRDout ) == 0 ) continue;
 399     if ( (status & AliESDtrack::kTRDin) != 0 ) continue;
 400     nseed++;
 401
 402     AliTRDtrack* seed2 = new AliTRDtrack(*seed);
 403     //seed2->ResetCovariance();
 404     AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
 405     AliTRDtrack &t=*pt;
 406     FollowProlongation(t, innerTB);
 407     if (t.GetNumberOfClusters() >= foundMin) {
 408       UseClusters(&t);
 409       CookLabel(pt, 1-fgkLabelFraction);
 410       //      t.CookdEdx();
 411     }
 412     found++;
 413 //    cout<<found<<'\r';
 414
 415     if(PropagateToTPC(t)) {
 416       seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
 417     }
 418     delete seed2;
 419     delete pt;
 420   }
 421
 422   cout<<"Number of loaded seeds: "<<nseed<<endl;
 423   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
 424
 425   // after tracks from loaded seeds are found and the corresponding
 426   // clusters are used, look for additional seeds from TRD
 427
 428   if(fAddTRDseeds) {
 429     // Find tracks for the seeds in the TRD
 430     Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 431
 432     Int_t nSteps = (Int_t) (fgkSeedDepth / fgkSeedStep);
 433     Int_t gap = (Int_t) (timeBins * fgkSeedGap);
 434     Int_t step = (Int_t) (timeBins * fgkSeedStep);
 435
 436     // make a first turn with tight cut on initial curvature
 437     for(Int_t turn = 1; turn <= 2; turn++) {
 438       if(turn == 2) {
 439         nSteps = (Int_t) (fgkSeedDepth / (3*fgkSeedStep));
 440         step = (Int_t) (timeBins * (3*fgkSeedStep));
 441       }
 442       for(Int_t i=0; i<nSteps; i++) {
 443         Int_t outer=timeBins-1-i*step;
 444         Int_t inner=outer-gap;
 445
 446         nseed=fSeeds->GetEntriesFast();
 447
 448         MakeSeeds(inner, outer, turn);
 449
 450         nseed=fSeeds->GetEntriesFast();
 451         //        printf("\n turn %d, step %d: number of seeds for TRD inward %d\n",
 452         //               turn, i, nseed);
 453
 454         for (Int_t i=0; i<nseed; i++) {
 455           AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
 456           FollowProlongation(t,innerTB);
 457           if (t.GetNumberOfClusters() >= foundMin) {
 458             UseClusters(&t);
 459             CookLabel(pt, 1-fgkLabelFraction);
 460             t.CookdEdx();
 461             found++;
 462 //            cout<<found<<'\r';
 463             if(PropagateToTPC(t)) {
 464               AliESDtrack track;
 465               track.UpdateTrackParams(pt,AliESDtrack::kTRDin);
 466               event->AddTrack(&track);
 467               //              track.SetTRDtrack(new AliTRDtrack(*pt));
 468             }
 469           }
 470           delete fSeeds->RemoveAt(i);
 471           fNseeds--;
 472         }
 473       }
 474     }
 475   }
 476
 477   cout<<"Total number of found tracks: "<<found<<endl;
 478
 479   return 0;
 480 }
 481
 482
 483
 484 //_____________________________________________________________________________
 485 Int_t AliTRDtracker::PropagateBack(AliESD* event) {
 486   //
 487   // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
 488   // backpropagated by the TPC tracker. Each seed is first propagated
 489   // to the TRD, and then its prolongation is searched in the TRD.
 490   // If sufficiently long continuation of the track is found in the TRD
 491   // the track is updated, otherwise it's stored as originaly defined
 492   // by the TPC tracker.
 493   //
 494
 495   Int_t found=0;
 496   Float_t foundMin = 20;
 497   Int_t n = event->GetNumberOfTracks();
 498   //
 499   //Sort tracks
 500   Float_t *quality =new Float_t[n];
 501   Int_t *index   =new Int_t[n];
 502   for (Int_t i=0; i<n; i++) {
 503     AliESDtrack* seed=event->GetTrack(i);
 504     Double_t covariance[15];
 505     seed->GetExternalCovariance(covariance);
 506     quality[i] = covariance[0]+covariance[2];
 507   }
 508   TMath::Sort(n,quality,index,kFALSE);
 509   //
 510   for (Int_t i=0; i<n; i++) {
 511     //    AliESDtrack* seed=event->GetTrack(i);
 512     AliESDtrack* seed=event->GetTrack(index[i]);
 513
 514     ULong_t status=seed->GetStatus();
 515     if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
 516     if ( (status & AliESDtrack::kTRDout) != 0 ) continue;
 517
 518     Int_t lbl = seed->GetLabel();
 519     AliTRDtrack *track = new AliTRDtrack(*seed);
 520     track->SetSeedLabel(lbl);
 521     seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); //make backup
 522     fNseeds++;
 523     Float_t p4     = track->GetC();
 524     //
 525     Int_t expectedClr = FollowBackProlongationG(*track);
 526     /*
 527       // only debug purpose
 528     if (track->GetNumberOfClusters()<expectedClr/3){
 529       AliTRDtrack *track1 = new AliTRDtrack(*seed);
 530       track1->SetSeedLabel(lbl);
 531       FollowBackProlongation(*track1);
 532       AliTRDtrack *track2= new AliTRDtrack(*seed);
 533       track->SetSeedLabel(lbl);
 534       FollowBackProlongation(*track2);
 535       delete track1;
 536       delete track2;
 537     }
 538     */
 539     if (TMath::Abs(track->GetC()-p4)/TMath::Abs(p4)<0.2 || TMath::Abs(track->GetPt())>0.8 ) {
 540       //
 541       //make backup for back propagation
 542       //
 543       Int_t foundClr = track->GetNumberOfClusters();
 544       if (foundClr >= foundMin) {
 545         track->CookdEdx();
 546         CookdEdxTimBin(*track);
 547         CookLabel(track, 1-fgkLabelFraction);
 548         if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
 549         if(track->GetChi2()/track->GetNumberOfClusters()<4) {   // sign only gold tracks
 550           if (seed->GetKinkIndex(0)==0&&TMath::Abs(track->GetPt())<1.5 ) UseClusters(track);
 551         }
 552         Bool_t isGold = kFALSE;
 553
 554         if (track->GetChi2()/track->GetNumberOfClusters()<5) {  //full gold track
 555           // seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
 556            if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
 557           isGold = kTRUE;
 558         }
 559         if (!isGold && track->GetNCross()==0&&track->GetChi2()/track->GetNumberOfClusters()<7){ //almost gold track
 560           //      seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
 561           if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
 562           isGold = kTRUE;
 563         }
 564         if (!isGold && track->GetBackupTrack()){
 565           if (track->GetBackupTrack()->GetNumberOfClusters()>foundMin&&
 566               (track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1))<7){
 567             seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
 568             isGold = kTRUE;
 569           }
 570         }
 571         if (track->StatusForTOF()>0 &&track->fNCross==0 && Float_t(track->fN)/Float_t(track->fNExpected)>0.4){
 572           seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
 573         }
 574       }
 575     }
 576     // Debug part of tracking
 577     TTreeSRedirector& cstream = *fDebugStreamer;
 578     Int_t eventNr = event->GetEventNumber();
 579     if (track->GetBackupTrack()){
 580       cstream<<"Tracks"<<
 581         "EventNr="<<eventNr<<
 582         "ESD.="<<seed<<
 583         "trd.="<<track<<
 584         "trdback.="<<track->GetBackupTrack()<<
 585         "\n";
 586     }else{
 587       cstream<<"Tracks"<<
 588         "EventNr="<<eventNr<<
 589         "ESD.="<<seed<<
 590         "trd.="<<track<<
 591         "trdback.="<<track<<
 592         "\n";
 593     }
 594     //
 595     //Propagation to the TOF (I.Belikov)
 596     if (track->GetStop()==kFALSE){
 597
 598       Double_t xtof=371.;
 599       Double_t c2=track->GetC()*xtof - track->GetEta();
 600       if (TMath::Abs(c2)>=0.99) {
 601         delete track;
 602         continue;
 603       }
 604       Double_t xTOF0 = 365. ;
 605       PropagateToOuterPlane(*track,xTOF0);
 606       //
 607       //energy losses taken to the account - check one more time
 608       c2=track->GetC()*xtof - track->GetEta();
 609       if (TMath::Abs(c2)>=0.99) {
 610         delete track;
 611         continue;
 612       }
 613
 614       //
 615       Double_t ymax=xtof*TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
 616       Double_t y=track->GetYat(xtof);
 617       if (y > ymax) {
 618         if (!track->Rotate(AliTRDgeometry::GetAlpha())) {
 619           delete track;
 620           continue;
 621         }
 622       } else if (y <-ymax) {
 623         if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
 624           delete track;
 625           continue;
 626         }
 627       }
 628
 629       if (track->PropagateTo(xtof)) {
 630         seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
 631         for (Int_t i=0;i<kNPlane;i++) {
 632            seed->SetTRDsignals(track->GetPIDsignals(i),i);
 633            seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
 634         }
 635         //      seed->SetTRDtrack(new AliTRDtrack(*track));
 636         if (track->GetNumberOfClusters()>foundMin) found++;
 637       }
 638     }else{
 639       if (track->GetNumberOfClusters()>15&&track->GetNumberOfClusters()>0.5*expectedClr){
 640         seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
 641         //seed->SetStatus(AliESDtrack::kTRDStop);
 642         for (Int_t i=0;i<kNPlane;i++) {
 643            seed->SetTRDsignals(track->GetPIDsignals(i),i);
 644            seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
 645         }
 646         //seed->SetTRDtrack(new AliTRDtrack(*track));
 647         found++;
 648       }
 649     }
 650     seed->SetTRDQuality(track->StatusForTOF());
 651     seed->SetTRDBudget(track->fBudget[0]);
 652
 653     delete track;
 654     //
 655     //End of propagation to the TOF
 656     //if (foundClr>foundMin)
 657     //  seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
 658
 659
 660   }
 661
 662   cerr<<"Number of seeds: "<<fNseeds<<endl;
 663   cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
 664
 665   //  MakeSeedsMI(3,5,event); //new seeding
 666
 667
 668   fSeeds->Clear(); fNseeds=0;
 669   delete [] index;
 670   delete [] quality;
 671
 672   return 0;
 673
 674 }
 675
 676 //_____________________________________________________________________________
 677 Int_t AliTRDtracker::RefitInward(AliESD* event)
 678 {
 679   //
 680   // Refits tracks within the TRD. The ESD event is expected to contain seeds
 681   // at the outer part of the TRD.
 682   // The tracks are propagated to the innermost time bin
 683   // of the TRD and the ESD event is updated
 684   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
 685   //
 686
 687   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 688   Float_t foundMin = fgkMinClustersInTrack * timeBins;
 689   Int_t nseed = 0;
 690   Int_t found = 0;
 691   Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
 692   AliTRDtrack seed2;
 693
 694   Int_t n = event->GetNumberOfTracks();
 695   for (Int_t i=0; i<n; i++) {
 696     AliESDtrack* seed=event->GetTrack(i);
 697     new(&seed2) AliTRDtrack(*seed);
 698     if (seed2.GetX()<270){
 699       seed->UpdateTrackParams(&seed2, AliESDtrack::kTRDbackup); // backup TPC track - only update
 700       continue;
 701     }
 702
 703     ULong_t status=seed->GetStatus();
 704     if ( (status & AliESDtrack::kTRDout ) == 0 ) {
 705       continue;
 706     }
 707     if ( (status & AliESDtrack::kTRDin) != 0 ) {
 708       continue;
 709     }
 710     nseed++;
 711 //     if (1/seed2.Get1Pt()>1.5&& seed2.GetX()>260.) {
 712 //       Double_t oldx = seed2.GetX();
 713 //       seed2.PropagateTo(500.);
 714 //       seed2.ResetCovariance(1.);
 715 //       seed2.PropagateTo(oldx);
 716 //     }
 717 //     else{
 718 //       seed2.ResetCovariance(5.);
 719 //     }
 720
 721     AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
 722     UInt_t * indexes2 = seed2.GetIndexes();
 723     for (Int_t i=0;i<kNPlane;i++) {
 724       pt->SetPIDsignals(seed2.GetPIDsignals(i),i);
 725       pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i);
 726     }
 727
 728     UInt_t * indexes3 = pt->GetBackupIndexes();
 729     for (Int_t i=0;i<200;i++) {
 730       if (indexes2[i]==0) break;
 731       indexes3[i] = indexes2[i];
 732     }
 733     //AliTRDtrack *pt = seed2;
 734     AliTRDtrack &t=*pt;
 735     FollowProlongationG(t, innerTB);
 736     if (t.GetNumberOfClusters() >= foundMin) {
 737       //      UseClusters(&t);
 738       //CookLabel(pt, 1-fgkLabelFraction);
 739       t.CookdEdx();
 740       CookdEdxTimBin(t);
 741     }
 742     found++;
 743 //    cout<<found<<'\r';
 744
 745     if(PropagateToTPC(t)) {
 746       seed->UpdateTrackParams(pt, AliESDtrack::kTRDrefit);
 747       for (Int_t i=0;i<kNPlane;i++) {
 748         seed->SetTRDsignals(pt->GetPIDsignals(i),i);
 749         seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
 750       }
 751     }else{
 752       //if not prolongation to TPC - propagate without update
 753       AliTRDtrack* seed2 = new AliTRDtrack(*seed);
 754       seed2->ResetCovariance(5.);
 755       AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha());
 756       delete seed2;
 757       if (PropagateToTPC(*pt2)) {
 758         //pt2->CookdEdx(0.,1.);
 759         pt2->CookdEdx( ); // Modification by PS
 760         CookdEdxTimBin(*pt2);
 761         seed->UpdateTrackParams(pt2, AliESDtrack::kTRDrefit);
 762         for (Int_t i=0;i<kNPlane;i++) {
 763           seed->SetTRDsignals(pt2->GetPIDsignals(i),i);
 764           seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
 765         }
 766       }
 767       delete pt2;
 768     }
 769     delete pt;
 770   }
 771
 772   cout<<"Number of loaded seeds: "<<nseed<<endl;
 773   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
 774
 775   return 0;
 776
 777 }
 778
 779
 780 //---------------------------------------------------------------------------
 781 Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
 782 {
 783   // Starting from current position on track=t this function tries
 784   // to extrapolate the track up to timeBin=0 and to confirm prolongation
 785   // if a close cluster is found. Returns the number of clusters
 786   // expected to be found in sensitive layers
 787
 788   Float_t  wIndex, wTB, wChi2;
 789   Float_t  wYrt, wYclosest, wYcorrect, wYwindow;
 790   Float_t  wZrt, wZclosest, wZcorrect, wZwindow;
 791   Float_t  wPx, wPy, wPz, wC;
 792   Double_t px, py, pz;
 793   Float_t  wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
 794   Int_t lastplane = GetLastPlane(&t);
 795   Double_t dxsample = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
 796                          / AliTRDcalibDB::Instance()->GetSamplingFrequency();
 797   Int_t trackIndex = t.GetLabel();
 798
 799   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
 800
 801   Int_t tryAgain=fMaxGap;
 802
 803   Double_t alpha=t.GetAlpha();
 804   alpha = TVector2::Phi_0_2pi(alpha);
 805
 806   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
 807   Double_t radLength, rho, x, dx, y, ymax, z;
 808
 809   Int_t expectedNumberOfClusters = 0;
 810   Bool_t lookForCluster;
 811
 812   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
 813
 814
 815   for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
 816
 817     y = t.GetY(); z = t.GetZ();
 818
 819     // first propagate to the inner surface of the current time bin
 820     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
 821     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
 822     if(!t.PropagateTo(x,radLength,rho)) break;
 823     y = t.GetY();
 824     ymax = x*TMath::Tan(0.5*alpha);
 825     if (y > ymax) {
 826       s = (s+1) % ns;
 827       if (!t.Rotate(alpha)) break;
 828       if(!t.PropagateTo(x,radLength,rho)) break;
 829     } else if (y <-ymax) {
 830       s = (s-1+ns) % ns;
 831       if (!t.Rotate(-alpha)) break;
 832       if(!t.PropagateTo(x,radLength,rho)) break;
 833     }
 834
 835     y = t.GetY(); z = t.GetZ();
 836
 837     // now propagate to the middle plane of the next time bin
 838     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
 839     x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
 840     if(!t.PropagateTo(x,radLength,rho)) break;
 841     y = t.GetY();
 842     ymax = x*TMath::Tan(0.5*alpha);
 843     if (y > ymax) {
 844       s = (s+1) % ns;
 845       if (!t.Rotate(alpha)) break;
 846       if(!t.PropagateTo(x,radLength,rho)) break;
 847     } else if (y <-ymax) {
 848       s = (s-1+ns) % ns;
 849       if (!t.Rotate(-alpha)) break;
 850       if(!t.PropagateTo(x,radLength,rho)) break;
 851     }
 852
 853
 854     if(lookForCluster) {
 855
 856       expectedNumberOfClusters++;
 857       wIndex = (Float_t) t.GetLabel();
 858       wTB = nr;
 859
 860       AliTRDpropagationLayer& timeBin=*(fTrSec[s]->GetLayer(nr-1));
 861
 862       Double_t sy2=ExpectedSigmaY2(x,t.GetTgl(),t.GetPt());
 863       Double_t sz2=ExpectedSigmaZ2(x,t.GetTgl());
 864
 865       Double_t road;
 866       if((t.GetSigmaY2() + sy2) > 0) road=10.*sqrt(t.GetSigmaY2() + sy2);
 867       else return expectedNumberOfClusters;
 868
 869       wYrt = (Float_t) y;
 870       wZrt = (Float_t) z;
 871       wYwindow = (Float_t) road;
 872       t.GetPxPyPz(px,py,pz);
 873       wPx = (Float_t) px;
 874       wPy = (Float_t) py;
 875       wPz = (Float_t) pz;
 876       wC  = (Float_t) t.GetC();
 877       wSigmaC2 = (Float_t) t.GetSigmaC2();
 878       wSigmaTgl2    = (Float_t) t.GetSigmaTgl2();
 879       wSigmaY2 = (Float_t) t.GetSigmaY2();
 880       wSigmaZ2 = (Float_t) t.GetSigmaZ2();
 881       wChi2 = -1;
 882
 883
 884       AliTRDcluster *cl=0;
 885       UInt_t index=0;
 886
 887       Double_t maxChi2=fgkMaxChi2;
 888
 889       wYclosest = 12345678;
 890       wYcorrect = 12345678;
 891       wZclosest = 12345678;
 892       wZcorrect = 12345678;
 893       wZwindow  = TMath::Sqrt(2.25 * 12 * sz2);
 894
 895       // Find the closest correct cluster for debugging purposes
 896       if (timeBin&&fVocal) {
 897         Float_t minDY = 1000000;
 898         for (Int_t i=0; i<timeBin; i++) {
 899           AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
 900           if((c->GetLabel(0) != trackIndex) &&
 901              (c->GetLabel(1) != trackIndex) &&
 902              (c->GetLabel(2) != trackIndex)) continue;
 903           if(TMath::Abs(c->GetY() - y) > minDY) continue;
 904           minDY = TMath::Abs(c->GetY() - y);
 905           wYcorrect = c->GetY();
 906           wZcorrect = c->GetZ();
 907
 908           Double_t h01 = GetTiltFactor(c);
 909           wChi2 = t.GetPredictedChi2(c, h01);
 910         }
 911       }
 912
 913       // Now go for the real cluster search
 914
 915       if (timeBin) {
 916         //
 917         //find cluster in history
 918         cl =0;
 919
 920         AliTRDcluster * cl0 = timeBin[0];
 921         if (!cl0) {
 922           continue;
 923         }
 924         Int_t plane = fGeom->GetPlane(cl0->GetDetector());
 925         if (plane>lastplane) continue;
 926         Int_t timebin = cl0->GetLocalTimeBin();
 927         AliTRDcluster * cl2= GetCluster(&t,plane, timebin,index);
 928         if (cl2) {
 929           cl =cl2;
 930           Double_t h01 = GetTiltFactor(cl);
 931           maxChi2=t.GetPredictedChi2(cl,h01);
 932         }
 933         if ((!cl) && road>fgkWideRoad) {
 934           //if (t.GetNumberOfClusters()>4)
 935           //  cerr<<t.GetNumberOfClusters()
 936           //    <<"FindProlongation warning: Too broad road !\n";
 937           continue;
 938         }
 939
 940         if (cl) {
 941
 942           wYclosest = cl->GetY();
 943           wZclosest = cl->GetZ();
 944           Double_t h01 = GetTiltFactor(cl);
 945
 946           //if (cl->GetNPads()<5)
 947             t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
 948           Int_t det = cl->GetDetector();
 949           Int_t plane = fGeom->GetPlane(det);
 950
 951           if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
 952             //if(!t.Update(cl,maxChi2,index,h01)) {
 953             //if(!tryAgain--) return 0;
 954           }
 955           else tryAgain=fMaxGap;
 956         }
 957         else {
 958           //if (tryAgain==0) break;
 959           tryAgain--;
 960         }
 961       }
 962     }
 963   }
 964   return expectedNumberOfClusters;
 965
 966
 967 }
 968
 969
 970
 971
 972 //---------------------------------------------------------------------------
 973 Int_t AliTRDtracker::FollowProlongationG(AliTRDtrack& t, Int_t rf)
 974 {
 975   // Starting from current position on track=t this function tries
 976   // to extrapolate the track up to timeBin=0 and to confirm prolongation
 977   // if a close cluster is found. Returns the number of clusters
 978   // expected to be found in sensitive layers
 979   // GeoManager used to estimate mean density
 980   Int_t sector;
 981   Int_t lastplane = GetLastPlane(&t);
 982   Int_t tryAgain=fMaxGap;
 983   Double_t alpha=t.GetAlpha();
 984   alpha = TVector2::Phi_0_2pi(alpha);
 985   Double_t radLength = 0.0;
 986   Double_t rho = 0.0;
 987   Double_t x, dx;
 988   //, y, ymax, z;
 989   Int_t expectedNumberOfClusters = 0;
 990   Double_t dxsample = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
 991                          / AliTRDcalibDB::Instance()->GetSamplingFrequency();
 992   //
 993   //
 994   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
 995   Double_t tanmax = TMath::Tan(0.5*alpha);
 996
 997   for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
 998     //
 999     // propagate track in  non active layers
1000     //
1001     if (!(fTrSec[0]->GetLayer(nr)->IsSensitive())){
1002       Double_t xyz0[3],xyz1[3],param[7],x,y,z;
1003       t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);   //starting global position
1004       while (nr >rf && (!(fTrSec[0]->GetLayer(nr)->IsSensitive()))){
1005         x = fTrSec[0]->GetLayer(nr)->GetX();
1006         nr--;
1007         if (!t.GetProlongation(x,y,z)) break;
1008         if (TMath::Abs(y)>x*tanmax){
1009           nr--;
1010           break;
1011         }
1012       }
1013       nr++;
1014       x = fTrSec[0]->GetLayer(nr)->GetX();
1015       if (!t.GetProlongation(x,y,z)) break;
1016       //
1017       // minimal mean and maximal budget scan
1018       Float_t minbudget  =10000;
1019       Float_t meanbudget =0;
1020       Float_t maxbudget  =-1;
1021 //      Float_t normbudget =0;
1022 //       for (Int_t idy=-1;idy<=1;idy++)
1023 //      for (Int_t idz=-1;idz<=1;idz++){
1024       for (Int_t idy=0;idy<1;idy++)
1025         for (Int_t idz=0;idz<1;idz++){
1026           Double_t y2 = y+idy*TMath::Min(TMath::Sqrt(t.GetSigmaY2()),1.);
1027           Double_t z2 = z+idz*TMath::Min(TMath::Sqrt(t.GetSigmaY2()),1.);
1028           xyz1[0] = x*TMath::Cos(t.GetAlpha())-y2*TMath::Sin(t.GetAlpha());
1029           xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y2*TMath::Cos(t.GetAlpha());
1030           xyz1[2] = z2;
1031           AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1032           Float_t budget = param[0]*param[4];
1033           meanbudget+=budget;
1034           if (budget<minbudget) minbudget=budget;
1035           if (budget>maxbudget) maxbudget=budget;
1036         }
1037       t.fBudget[0]+=minbudget;
1038       t.fBudget[1]+=meanbudget/9.;
1039       t.fBudget[2]+=minbudget;
1040       //
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       AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1046
1047       t.PropagateTo(x,param[1],param[0]);
1048       t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);   //end  global position
1049       AdjustSector(&t);
1050       continue;
1051     }
1052      //
1053     //
1054     // stop tracking for highly inclined tracks
1055     if (!AdjustSector(&t)) break;
1056     if (TMath::Abs(t.GetSnp())>0.95) break;
1057     //
1058     // propagate and update track in active layers
1059     //
1060     Int_t nr0 = nr;  //first active layer
1061     if (nr >rf  && (fTrSec[0]->GetLayer(nr)->IsSensitive())){
1062       Double_t xyz0[3],xyz1[3],param[7],x,y,z;
1063       t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);   //starting global position
1064       while (nr >rf && ((fTrSec[0]->GetLayer(nr)->IsSensitive()))){
1065         x = fTrSec[0]->GetLayer(nr)->GetX();
1066         nr--;
1067         if (!t.GetProlongation(x,y,z)) break;
1068         if (TMath::Abs(y)>x*tanmax){
1069           nr--;
1070           break;
1071         }
1072       }
1073       //      nr++;
1074       x = fTrSec[0]->GetLayer(nr)->GetX();
1075       if (!t.GetProlongation(x,y,z)) break;
1076       xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
1077       xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
1078       xyz1[2] = z;
1079       // end global position
1080       AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1081       rho = param[0];
1082       radLength = param[1];   // get mean propagation parameters
1083     }
1084     //
1085     // propagate and update
1086     if (nr0-nr<5){
1087       // short tracklet - do not update - edge effect
1088       x = fTrSec[0]->GetLayer(nr)->GetX();
1089       t.PropagateTo(x,radLength,rho);
1090       AdjustSector(&t);
1091       continue;
1092     }
1093     sector = t.GetSector();
1094     //
1095     //
1096     for (Int_t ilayer=nr0;ilayer>=nr;ilayer--) {
1097       expectedNumberOfClusters++;
1098       t.fNExpected++;
1099       if (t.fX>345) t.fNExpectedLast++;
1100       AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(ilayer));
1101       AliTRDcluster *cl=0;
1102       UInt_t index=0;
1103       Double_t maxChi2=fgkMaxChi2;
1104       dx = (fTrSec[sector]->GetLayer(ilayer+1))->GetX()-timeBin.GetX();
1105       x = timeBin.GetX();
1106       t.PropagateTo(x,radLength,rho);
1107       // Now go for the real cluster search
1108       if (timeBin) {
1109         AliTRDcluster * cl0 = timeBin[0];
1110         if (!cl0) continue;         // no clusters in given time bin
1111         Int_t plane = fGeom->GetPlane(cl0->GetDetector());
1112         if (plane>lastplane) continue;
1113         Int_t timebin = cl0->GetLocalTimeBin();
1114         AliTRDcluster * cl2= GetCluster(&t,plane, timebin,index);
1115         //
1116         if (cl2) {
1117           cl =cl2;
1118           Double_t h01 = GetTiltFactor(cl);
1119           maxChi2=t.GetPredictedChi2(cl,h01);
1120         }
1121
1122         if (cl) {
1123           //      if (cl->GetNPads()<5)
1124             t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1125           Double_t h01 = GetTiltFactor(cl);
1126           Int_t det = cl->GetDetector();
1127           Int_t plane = fGeom->GetPlane(det);
1128           if (t.fX>345){
1129             t.fNLast++;
1130             t.fChi2Last+=maxChi2;
1131           }
1132           if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1133             if(!t.Update(cl,maxChi2,index,h01)) {
1134               //if(!tryAgain--) return 0;
1135             }
1136           }
1137           else tryAgain=fMaxGap;
1138           //
1139         }
1140       }
1141     }
1142   }
1143   return expectedNumberOfClusters;
1144
1145
1146 }
1147
1148 //___________________________________________________________________
1149
1150 Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
1151 {
1152   // Starting from current radial position of track <t> this function
1153   // extrapolates the track up to outer timebin and in the sensitive
1154   // layers confirms prolongation if a close cluster is found.
1155   // Returns the number of clusters expected to be found in sensitive layers
1156
1157   Int_t tryAgain=fMaxGap;
1158
1159   Double_t alpha=t.GetAlpha();
1160   TVector2::Phi_0_2pi(alpha);
1161
1162   Int_t s;
1163
1164   Int_t clusters[1000];
1165   for (Int_t i=0;i<1000;i++) clusters[i]=-1;
1166
1167   Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
1168   //Double_t radLength, rho, x, dx, y, ymax = 0, z;
1169   Double_t radLength, rho, x, dx, y, z;
1170   Bool_t lookForCluster;
1171   Double_t dxsample = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
1172                          / AliTRDcalibDB::Instance()->GetSamplingFrequency();
1173
1174   Int_t expectedNumberOfClusters = 0;
1175   x = t.GetX();
1176
1177   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1178
1179   //  Int_t zone =0;
1180   Int_t nr;
1181   Float_t ratio0=0;
1182   AliTRDtracklet tracklet;
1183   //
1184   for (nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB+1; nr++) {
1185     y = t.GetY();
1186     z = t.GetZ();
1187     // first propagate to the outer surface of the current time bin
1188
1189     s = t.GetSector();
1190     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1191     x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2;
1192     y = t.GetY();
1193     z = t.GetZ();
1194
1195     if(!t.PropagateTo(x,radLength,rho)) break;
1196     //
1197     // MI -fix untill correct material desription will be implemented
1198     //
1199     //Int_t nrotate = t.GetNRotate();
1200     if (!AdjustSector(&t)) break;
1201     //
1202     //
1203     y = t.GetY();
1204     z = t.GetZ();
1205     s = t.GetSector();
1206
1207     // now propagate to the middle plane of the next time bin
1208     fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1209 //     if (nrotate!=t.GetNRotate()){
1210 //       rho = 1000*2.7; radLength = 24.01;  //TEMPORARY - aluminium in between z - will be detected using GeoModeler in future versions
1211 //     }
1212     x = fTrSec[s]->GetLayer(nr+1)->GetX();
1213     if(!t.PropagateTo(x,radLength,rho)) break;
1214     if (!AdjustSector(&t)) break;
1215     s = t.GetSector();
1216     //    if(!t.PropagateTo(x,radLength,rho)) break;
1217
1218     if (TMath::Abs(t.GetSnp())>0.95) break;
1219
1220     y = t.GetY();
1221     z = t.GetZ();
1222
1223     if(lookForCluster) {
1224       if (clusters[nr]==-1) {
1225         Float_t  ncl   = FindClusters(s,nr,nr+30,&t,clusters,tracklet);
1226         ratio0 = ncl/Float_t(fTimeBinsPerPlane);
1227         Float_t  ratio1 = Float_t(t.fN+1)/Float_t(t.fNExpected+1.);
1228         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){
1229           t.MakeBackupTrack();                            // make backup of the track until is gold
1230         }
1231 //      if (ncl>4){
1232 //        t.PropagateTo(tracklet.GetX());
1233 //        t.UpdateMI(tracklet);
1234 //        nr = fTrSec[0]->GetLayerNumber(t.GetX())+1;
1235 //        continue;
1236 //      }
1237       }
1238
1239       expectedNumberOfClusters++;
1240       t.fNExpected++;
1241       if (t.fX>345) t.fNExpectedLast++;
1242
1243       AliTRDpropagationLayer& timeBin=*(fTrSec[s]->GetLayer(nr+1));
1244       Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
1245       if((t.GetSigmaY2() + sy2) < 0) {
1246         printf("problem\n");
1247         break;
1248       }
1249       Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2);
1250
1251       if (road>fgkWideRoad) {
1252         return 0;
1253       }
1254
1255       AliTRDcluster *cl=0;
1256       UInt_t index=0;
1257       Double_t maxChi2=fgkMaxChi2;
1258
1259       // Now go for the real cluster search
1260
1261       if (timeBin) {
1262
1263         if (clusters[nr+1]>0) {
1264           index = clusters[nr+1];
1265           cl    = (AliTRDcluster*)GetCluster(index);
1266           Double_t h01 = GetTiltFactor(cl);
1267           maxChi2=t.GetPredictedChi2(cl,h01);
1268         }
1269
1270         if (cl) {
1271           //      if (cl->GetNPads()<5)
1272             t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1273           Double_t h01 = GetTiltFactor(cl);
1274           Int_t det = cl->GetDetector();
1275           Int_t plane = fGeom->GetPlane(det);
1276           if (t.fX>345){
1277             t.fNLast++;
1278             t.fChi2Last+=maxChi2;
1279           }
1280           if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1281             if(!t.Update(cl,maxChi2,index,h01)) {
1282               //if(!tryAgain--) return 0;
1283             }
1284           }
1285           else tryAgain=fMaxGap;
1286           //
1287
1288           if (cl->GetLocalTimeBin()==1&&t.fN>20 && float(t.fChi2)/float(t.fN)<5){
1289             Float_t  ratio1 = Float_t(t.fN)/Float_t(t.fNExpected);
1290             if (tracklet.GetChi2()<18&&ratio0>0.8&&ratio1>0.6 &&ratio0+ratio1>1.5 && t.GetNCross()==0 && TMath::Abs(t.GetSnp())<0.85){
1291               t.MakeBackupTrack();                            // make backup of the track until is gold
1292             }
1293           }
1294
1295         }
1296         else {
1297           // if (tryAgain==0) break;
1298           //tryAgain--;
1299         }
1300       }
1301     }
1302   }
1303   if (nr<outerTB)
1304     t.SetStop(kTRUE);
1305   else
1306     t.SetStop(kFALSE);
1307   return expectedNumberOfClusters;
1308 }
1309
1310
1311
1312 //___________________________________________________________________
1313 Int_t AliTRDtracker::FollowBackProlongationG(AliTRDtrack& t)
1314 {
1315
1316   // Starting from current radial position of track <t> this function
1317   // extrapolates the track up to outer timebin and in the sensitive
1318   // layers confirms prolongation if a close cluster is found.
1319   // Returns the number of clusters expected to be found in sensitive layers
1320   // Use GEO manager for material Description
1321   Int_t tryAgain=fMaxGap;
1322
1323   Double_t alpha=t.GetAlpha();
1324   TVector2::Phi_0_2pi(alpha);
1325   Int_t sector;
1326   Int_t clusters[1000];
1327   for (Int_t i=0;i<1000;i++) clusters[i]=-1;
1328   Double_t dxsample = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
1329     / AliTRDcalibDB::Instance()->GetSamplingFrequency();
1330   Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
1331   Double_t radLength = 0.0;
1332   Double_t rho = 0.0;
1333   Double_t x, dx; //y, z;
1334
1335   Int_t expectedNumberOfClusters = 0;
1336   x = t.GetX();
1337
1338   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1339   Double_t tanmax = TMath::Tan(0.5*alpha);
1340   Int_t nr;
1341   Float_t ratio0=0;
1342   AliTRDtracklet tracklet;
1343   //
1344   //
1345   //
1346   for (nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB+1; nr++) {
1347     //
1348     // propagate track in  non active layers
1349     //
1350     if (!(fTrSec[0]->GetLayer(nr)->IsSensitive())){
1351       Double_t xyz0[3],xyz1[3],param[7],x,y,z;
1352       t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);   //starting global position
1353       while (nr <outerTB && (!(fTrSec[0]->GetLayer(nr)->IsSensitive()))){
1354         x = fTrSec[0]->GetLayer(nr)->GetX();
1355         nr++;
1356         if (!t.GetProlongation(x,y,z)) break;
1357         if (TMath::Abs(y)>x*tanmax){
1358           nr++;
1359           break;
1360         }
1361       }
1362       nr--;
1363       x = fTrSec[0]->GetLayer(nr)->GetX();
1364       if (!t.GetProlongation(x,y,z)) break;
1365       // minimal mean and maximal budget scan
1366       Float_t minbudget  =10000;
1367       Float_t meanbudget =0;
1368       Float_t maxbudget  =-1;
1369 //      Float_t normbudget =0;
1370 //       for (Int_t idy=-1;idy<=1;idy++)
1371 //      for (Int_t idz=-1;idz<=1;idz++){
1372       for (Int_t idy=0;idy<1;idy++)
1373         for (Int_t idz=0;idz<1;idz++){
1374           Double_t y2 = y+idy*TMath::Min(TMath::Sqrt(t.GetSigmaY2()),1.);
1375           Double_t z2 = z+idz*TMath::Min(TMath::Sqrt(t.GetSigmaZ2()),1.);
1376
1377           xyz1[0] = x*TMath::Cos(t.GetAlpha())-y2*TMath::Sin(t.GetAlpha());
1378           xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y2*TMath::Cos(t.GetAlpha());
1379           xyz1[2] = z2;
1380           AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1381           Float_t budget = param[0]*param[4];
1382           meanbudget+=budget;
1383           if (budget<minbudget) minbudget=budget;
1384           if (budget>maxbudget) maxbudget=budget;
1385         }
1386       t.fBudget[0]+=minbudget;
1387       t.fBudget[1]+=meanbudget/9.;
1388       t.fBudget[2]+=minbudget;
1389
1390       xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
1391       xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
1392       xyz1[2] = z;
1393       AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1394       t.PropagateTo(x,param[1],param[0]);
1395       t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);   //end  global position
1396       AdjustSector(&t);
1397       continue;
1398     }
1399     //
1400     //
1401     // stop tracking for highly inclined tracks
1402     if (!AdjustSector(&t)) break;
1403     if (TMath::Abs(t.GetSnp())>0.95) break;
1404     //
1405     // propagate and update track in active layers
1406     //
1407     Int_t nr0 = nr;  //first active layer
1408     if (nr <outerTB && (fTrSec[0]->GetLayer(nr)->IsSensitive())){
1409       Double_t xyz0[3],xyz1[3],param[7],x,y,z;
1410       t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);   //starting global position
1411       while (nr <outerTB && ((fTrSec[0]->GetLayer(nr)->IsSensitive()))){
1412         x = fTrSec[0]->GetLayer(nr)->GetX();
1413         nr++;
1414         if (!t.GetProlongation(x,y,z)) break;
1415         if (TMath::Abs(y)>(x*tanmax)){
1416           nr++;
1417           break;
1418         }
1419       }
1420       x = fTrSec[0]->GetLayer(nr)->GetX();
1421       if (!t.GetProlongation(x,y,z)) break;
1422        // minimal mean and maximal budget scan
1423       Float_t minbudget  =10000;
1424       Float_t meanbudget =0;
1425       Float_t maxbudget  =-1;
1426       //      Float_t normbudget =0;
1427       //     for (Int_t idy=-1;idy<=1;idy++)
1428       //        for (Int_t idz=-1;idz<=1;idz++){
1429       for (Int_t idy=0;idy<1;idy++)
1430         for (Int_t idz=0;idz<1;idz++){
1431           Double_t y2 = y+idy*TMath::Min(TMath::Sqrt(t.GetSigmaY2()),1.);
1432           Double_t z2 = z+idz*TMath::Min(TMath::Sqrt(t.GetSigmaZ2()),1.);
1433
1434           xyz1[0] = x*TMath::Cos(t.GetAlpha())-y2*TMath::Sin(t.GetAlpha());
1435           xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y2*TMath::Cos(t.GetAlpha());
1436           xyz1[2] = z2;
1437           AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1438           Float_t budget = param[0]*param[4];
1439           meanbudget+=budget;
1440           if (budget<minbudget) minbudget=budget;
1441           if (budget>maxbudget) maxbudget=budget;
1442         }
1443       t.fBudget[0]+=minbudget;
1444       t.fBudget[1]+=meanbudget/9.;
1445       t.fBudget[2]+=minbudget;
1446       //
1447       xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
1448       xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
1449       xyz1[2] = z;
1450       // end global position
1451       AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1452       rho = param[0];
1453       radLength = param[1];   // get mean propagation parameters
1454     }
1455     //
1456     //
1457     if (nr-nr0<5){
1458       // short tracklet - do not update - edge effect
1459       x = fTrSec[0]->GetLayer(nr+1)->GetX();
1460       t.PropagateTo(x,radLength,rho);
1461       AdjustSector(&t);
1462       continue;
1463     }
1464     //
1465     //
1466     sector = t.GetSector();
1467     Float_t  ncl   = FindClusters(sector,nr0,nr,&t,clusters,tracklet);
1468     if (tracklet.GetN()-2*tracklet.GetNCross()<10) continue;
1469     ratio0 = ncl/Float_t(fTimeBinsPerPlane);
1470     Float_t  ratio1 = Float_t(t.fN+1)/Float_t(t.fNExpected+1.);
1471     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){
1472       t.MakeBackupTrack();                            // make backup of the track until is gold
1473     }
1474     //
1475     //
1476     for (Int_t ilayer=nr0;ilayer<=nr;ilayer++) {
1477       expectedNumberOfClusters++;
1478       t.fNExpected++;
1479       if (t.fX>345) t.fNExpectedLast++;
1480       AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(ilayer));
1481       AliTRDcluster *cl=0;
1482       UInt_t index=0;
1483       Double_t maxChi2=fgkMaxChi2;
1484       dx = (fTrSec[sector]->GetLayer(ilayer-1))->GetX()-timeBin.GetX();
1485       x = timeBin.GetX();
1486       t.PropagateTo(x,radLength,rho);
1487       // Now go for the real cluster search
1488       if (timeBin) {
1489         if (clusters[ilayer]>0) {
1490           index = clusters[ilayer];
1491           cl    = (AliTRDcluster*)GetCluster(index);
1492           Double_t h01 = GetTiltFactor(cl);
1493           maxChi2=t.GetPredictedChi2(cl,h01);
1494         }
1495
1496         if (cl) {
1497           //      if (cl->GetNPads()<5)
1498             t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1499           Double_t h01 = GetTiltFactor(cl);
1500           Int_t det = cl->GetDetector();
1501           Int_t plane = fGeom->GetPlane(det);
1502           if (t.fX>345){
1503             t.fNLast++;
1504             t.fChi2Last+=maxChi2;
1505           }
1506           if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
1507             if(!t.Update(cl,maxChi2,index,h01)) {
1508               //if(!tryAgain--) return 0;
1509             }
1510           }
1511           else tryAgain=fMaxGap;
1512           //
1513
1514           if (cl->GetLocalTimeBin()==1&&t.fN>20 && float(t.fChi2)/float(t.fN)<5){
1515             Float_t  ratio1 = Float_t(t.fN)/Float_t(t.fNExpected);
1516             if (tracklet.GetChi2()<18&&ratio0>0.8&&ratio1>0.6 &&ratio0+ratio1>1.5 && t.GetNCross()==0 && TMath::Abs(t.GetSnp())<0.85){
1517               t.MakeBackupTrack();                            // make backup of the track until is gold
1518             }
1519           }
1520           // reset material budget if 2 consecutive gold
1521           if (plane>0)
1522             if (t.fTracklets[plane].GetN()+t.fTracklets[plane-1].GetN()>20){
1523               t.fBudget[2] = 0;
1524             }
1525         }
1526       }
1527     }
1528   }
1529   //
1530   if (nr<outerTB)
1531     t.SetStop(kTRUE);
1532   else
1533     t.SetStop(kFALSE);
1534   return expectedNumberOfClusters;
1535 }
1536
1537 //---------------------------------------------------------------------------
1538 Int_t AliTRDtracker::Refit(AliTRDtrack& t, Int_t rf)
1539 {
1540   // Starting from current position on track=t this function tries
1541   // to extrapolate the track up to timeBin=0 and to reuse already
1542   // assigned clusters. Returns the number of clusters
1543   // expected to be found in sensitive layers
1544   // get indices of assigned clusters for each layer
1545   // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
1546   Double_t dxsample = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
1547     / AliTRDcalibDB::Instance()->GetSamplingFrequency();
1548   Int_t iCluster[90];
1549   for (Int_t i = 0; i < 90; i++) iCluster[i] = 0;
1550   for (Int_t i = 0; i < t.GetNumberOfClusters(); i++) {
1551     Int_t index = t.GetClusterIndex(i);
1552     AliTRDcluster *cl=(AliTRDcluster*) GetCluster(index);
1553     if (!cl) continue;
1554     Int_t detector=cl->GetDetector();
1555     Int_t localTimeBin=cl->GetLocalTimeBin();
1556     Int_t sector=fGeom->GetSector(detector);
1557     Int_t plane=fGeom->GetPlane(detector);
1558
1559     Int_t trackingSector = CookSectorIndex(sector);
1560
1561     Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
1562     if(gtb < 0) continue;
1563     Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1564     iCluster[layer] = index;
1565   }
1566   t.ResetClusters();
1567
1568   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1569
1570   Double_t alpha=t.GetAlpha();
1571   alpha = TVector2::Phi_0_2pi(alpha);
1572
1573   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1574   Double_t radLength, rho, x, dx, y, ymax, z;
1575
1576   Int_t expectedNumberOfClusters = 0;
1577   Bool_t lookForCluster;
1578
1579   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1580
1581
1582   for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
1583
1584     y = t.GetY(); z = t.GetZ();
1585
1586     // first propagate to the inner surface of the current time bin
1587     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1588     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
1589     if(!t.PropagateTo(x,radLength,rho)) break;
1590     y = t.GetY();
1591     ymax = x*TMath::Tan(0.5*alpha);
1592     if (y > ymax) {
1593       s = (s+1) % ns;
1594       if (!t.Rotate(alpha)) break;
1595       if(!t.PropagateTo(x,radLength,rho)) break;
1596     } else if (y <-ymax) {
1597       s = (s-1+ns) % ns;
1598       if (!t.Rotate(-alpha)) break;
1599       if(!t.PropagateTo(x,radLength,rho)) break;
1600     }
1601
1602     y = t.GetY(); z = t.GetZ();
1603
1604     // now propagate to the middle plane of the next time bin
1605     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1606     x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
1607     if(!t.PropagateTo(x,radLength,rho)) break;
1608     y = t.GetY();
1609     ymax = x*TMath::Tan(0.5*alpha);
1610     if (y > ymax) {
1611       s = (s+1) % ns;
1612       if (!t.Rotate(alpha)) break;
1613       if(!t.PropagateTo(x,radLength,rho)) break;
1614     } else if (y <-ymax) {
1615       s = (s-1+ns) % ns;
1616       if (!t.Rotate(-alpha)) break;
1617       if(!t.PropagateTo(x,radLength,rho)) break;
1618     }
1619
1620     if(lookForCluster) expectedNumberOfClusters++;
1621
1622     // use assigned cluster
1623     if (!iCluster[nr-1]) continue;
1624     AliTRDcluster *cl=(AliTRDcluster*)GetCluster(iCluster[nr-1]);
1625     Double_t h01 = GetTiltFactor(cl);
1626     Double_t chi2=t.GetPredictedChi2(cl, h01);
1627     //if (cl->GetNPads()<5)
1628       t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
1629
1630       //t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
1631     t.Update(cl,chi2,iCluster[nr-1],h01);
1632   }
1633
1634   return expectedNumberOfClusters;
1635 }
1636
1637 //___________________________________________________________________
1638
1639 Int_t AliTRDtracker::PropagateToOuterPlane(AliTRDtrack& t, Double_t xToGo)
1640 {
1641   // Starting from current radial position of track <t> this function
1642   // extrapolates the track up to radial position <xToGo>.
1643   // Returns 1 if track reaches the plane, and 0 otherwise
1644
1645   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1646
1647   Double_t alpha=t.GetAlpha();
1648
1649   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1650   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1651
1652   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1653
1654   Bool_t lookForCluster;
1655   Double_t radLength, rho, x, dx, y, ymax, z;
1656
1657   x = t.GetX();
1658
1659   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1660
1661   Int_t plToGo = fTrSec[0]->GetLayerNumber(xToGo);
1662
1663   for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr<plToGo; nr++) {
1664
1665     y = t.GetY(); z = t.GetZ();
1666
1667     // first propagate to the outer surface of the current time bin
1668     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1669     x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
1670     if(!t.PropagateTo(x,radLength,rho)) return 0;
1671     y = t.GetY();
1672     ymax = x*TMath::Tan(0.5*alpha);
1673     if (y > ymax) {
1674       s = (s+1) % ns;
1675       if (!t.Rotate(alpha)) return 0;
1676     } else if (y <-ymax) {
1677       s = (s-1+ns) % ns;
1678       if (!t.Rotate(-alpha)) return 0;
1679     }
1680     if(!t.PropagateTo(x,radLength,rho)) return 0;
1681
1682     y = t.GetY(); z = t.GetZ();
1683
1684     // now propagate to the middle plane of the next time bin
1685     fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1686     x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
1687     if(!t.PropagateTo(x,radLength,rho)) return 0;
1688     y = t.GetY();
1689     ymax = x*TMath::Tan(0.5*alpha);
1690     if (y > ymax) {
1691       s = (s+1) % ns;
1692       if (!t.Rotate(alpha)) return 0;
1693     } else if (y <-ymax) {
1694       s = (s-1+ns) % ns;
1695       if (!t.Rotate(-alpha)) return 0;
1696     }
1697     if(!t.PropagateTo(x,radLength,rho)) return 0;
1698   }
1699   return 1;
1700 }
1701
1702 //___________________________________________________________________
1703
1704 Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t)
1705 {
1706   // Starting from current radial position of track <t> this function
1707   // extrapolates the track up to radial position of the outermost
1708   // padrow of the TPC.
1709   // Returns 1 if track reaches the TPC, and 0 otherwise
1710
1711   //Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1712
1713   Double_t alpha=t.GetAlpha();
1714   alpha = TVector2::Phi_0_2pi(alpha);
1715
1716   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1717
1718   Bool_t lookForCluster;
1719   Double_t radLength, rho, x, dx, y, /*ymax,*/ z;
1720
1721   x = t.GetX();
1722
1723   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1724   Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055);
1725
1726   for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) {
1727
1728     y = t.GetY();
1729     z = t.GetZ();
1730
1731     // first propagate to the outer surface of the current time bin
1732     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1733     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2;
1734
1735     if(!t.PropagateTo(x,radLength,rho)) return 0;
1736     AdjustSector(&t);
1737     if(!t.PropagateTo(x,radLength,rho)) return 0;
1738
1739     y = t.GetY();
1740     z = t.GetZ();
1741
1742     // now propagate to the middle plane of the next time bin
1743     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
1744     x = fTrSec[s]->GetLayer(nr-1)->GetX();
1745
1746     if(!t.PropagateTo(x,radLength,rho)) return 0;
1747     AdjustSector(&t);
1748     if(!t.PropagateTo(x,radLength,rho)) return 0;
1749   }
1750   return 1;
1751 }
1752
1753 //_____________________________________________________________________________
1754 Int_t AliTRDtracker::LoadClusters(TTree *cTree)
1755 {
1756   // Fills clusters into TRD tracking_sectors
1757   // Note that the numbering scheme for the TRD tracking_sectors
1758   // differs from that of TRD sectors
1759   cout<<"\n Read Sectors  clusters"<<endl;
1760   if (ReadClusters(fClusters,cTree)) {
1761      Error("LoadClusters","Problem with reading the clusters !");
1762      return 1;
1763   }
1764   Int_t ncl=fClusters->GetEntriesFast();
1765   fNclusters=ncl;
1766   cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1767
1768   UInt_t index;
1769   for (Int_t ichamber=0;ichamber<5;ichamber++)
1770     for (Int_t isector=0;isector<18;isector++){
1771       fHoles[ichamber][isector]=kTRUE;
1772     }
1773
1774
1775   while (ncl--) {
1776 //    printf("\r %d left  ",ncl);
1777     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
1778     Int_t detector=c->GetDetector();
1779     Int_t localTimeBin=c->GetLocalTimeBin();
1780     Int_t sector=fGeom->GetSector(detector);
1781     Int_t plane=fGeom->GetPlane(detector);
1782
1783     Int_t trackingSector = CookSectorIndex(sector);
1784     if (c->GetLabel(0)>0){
1785       Int_t chamber = fGeom->GetChamber(detector);
1786       fHoles[chamber][trackingSector]=kFALSE;
1787     }
1788
1789     Int_t gtb = fTrSec[trackingSector]->CookTimeBinIndex(plane,localTimeBin);
1790     if(gtb < 0) continue;
1791     Int_t layer = fTrSec[trackingSector]->GetLayerNumber(gtb);
1792
1793     index=ncl;
1794     //
1795     // apply pos correction
1796     Float_t poscor =  fgkCoef*(c->GetLocalTimeBin() - fgkMean)+fgkOffset;
1797     c->SetY(c->GetY()-poscor);
1798     fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
1799   }
1800   //  printf("\r\n");
1801   //
1802   //
1803   /*
1804   for (Int_t isector=0;isector<18;isector++){
1805     for (Int_t ichamber=0;ichamber<5;ichamber++)
1806       if (fHoles[ichamber][isector]!=fGeom->IsHole(0,ichamber,17-isector))
1807         printf("Problem \t%d\t%d\t%d\t%d\n",isector,ichamber,fHoles[ichamber][isector],
1808              fGeom->IsHole(0,ichamber,17-isector));
1809   }
1810   */
1811   return 0;
1812 }
1813
1814 //_____________________________________________________________________________
1815 void AliTRDtracker::UnloadClusters()
1816 {
1817   //
1818   // Clears the arrays of clusters and tracks. Resets sectors and timebins
1819   //
1820
1821   Int_t i, nentr;
1822
1823   nentr = fClusters->GetEntriesFast();
1824   for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i);
1825   fNclusters = 0;
1826
1827   nentr = fSeeds->GetEntriesFast();
1828   for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i);
1829
1830   nentr = fTracks->GetEntriesFast();
1831   for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i);
1832
1833   Int_t nsec = AliTRDgeometry::kNsect;
1834
1835   for (i = 0; i < nsec; i++) {
1836     for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1837       fTrSec[i]->GetLayer(pl)->Clear();
1838     }
1839   }
1840
1841 }
1842
1843 //__________________________________________________________________________
1844 void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn)
1845 {
1846   // Creates track seeds using clusters in timeBins=i1,i2
1847
1848   if(turn > 2) {
1849     cerr<<"MakeSeeds: turn "<<turn<<" exceeds the limit of 2"<<endl;
1850     return;
1851   }
1852
1853   Double_t x[5], c[15];
1854   Int_t maxSec=AliTRDgeometry::kNsect;
1855   Double_t alpha=AliTRDgeometry::GetAlpha();
1856   Double_t shift=AliTRDgeometry::GetAlpha()/2.;
1857   Double_t cs=cos(alpha), sn=sin(alpha);
1858   Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
1859   Int_t i2 = fTrSec[0]->GetLayerNumber(inner);
1860   Int_t i1 = fTrSec[0]->GetLayerNumber(outer);
1861   Double_t x1 =fTrSec[0]->GetX(i1);
1862   Double_t xx2=fTrSec[0]->GetX(i2);
1863
1864   for (Int_t ns=0; ns<maxSec; ns++) {
1865
1866     Int_t nl2 = *(fTrSec[(ns-2+maxSec)%maxSec]->GetLayer(i2));
1867     Int_t nl=(*fTrSec[(ns-1+maxSec)%maxSec]->GetLayer(i2));
1868     Int_t nm=(*fTrSec[ns]->GetLayer(i2));
1869     Int_t nu=(*fTrSec[(ns+1)%maxSec]->GetLayer(i2));
1870     Int_t nu2=(*fTrSec[(ns+2)%maxSec]->GetLayer(i2));
1871
1872     AliTRDpropagationLayer& r1=*(fTrSec[ns]->GetLayer(i1));
1873
1874     for (Int_t is=0; is < r1; is++) {
1875       Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
1876
1877       for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
1878
1879         const AliTRDcluster *cl;
1880         Double_t x2,   y2,   z2;
1881         Double_t x3=0., y3=0.;
1882
1883         if (js<nl2) {
1884           if(turn != 2) continue;
1885           AliTRDpropagationLayer& r2=*(fTrSec[(ns-2+maxSec)%maxSec]->GetLayer(i2));
1886           cl=r2[js];
1887           y2=cl->GetY(); z2=cl->GetZ();
1888
1889           x2= xx2*cs2+y2*sn2;
1890           y2=-xx2*sn2+y2*cs2;
1891         }
1892         else if (js<nl2+nl) {
1893           if(turn != 1) continue;
1894           AliTRDpropagationLayer& r2=*(fTrSec[(ns-1+maxSec)%maxSec]->GetLayer(i2));
1895           cl=r2[js-nl2];
1896           y2=cl->GetY(); z2=cl->GetZ();
1897
1898           x2= xx2*cs+y2*sn;
1899           y2=-xx2*sn+y2*cs;
1900         }
1901         else if (js<nl2+nl+nm) {
1902           if(turn != 1) continue;
1903           AliTRDpropagationLayer& r2=*(fTrSec[ns]->GetLayer(i2));
1904           cl=r2[js-nl2-nl];
1905           x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
1906         }
1907         else if (js<nl2+nl+nm+nu) {
1908           if(turn != 1) continue;
1909           AliTRDpropagationLayer& r2=*(fTrSec[(ns+1)%maxSec]->GetLayer(i2));
1910           cl=r2[js-nl2-nl-nm];
1911           y2=cl->GetY(); z2=cl->GetZ();
1912
1913           x2=xx2*cs-y2*sn;
1914           y2=xx2*sn+y2*cs;
1915         }
1916         else {
1917           if(turn != 2) continue;
1918           AliTRDpropagationLayer& r2=*(fTrSec[(ns+2)%maxSec]->GetLayer(i2));
1919           cl=r2[js-nl2-nl-nm-nu];
1920           y2=cl->GetY(); z2=cl->GetZ();
1921
1922           x2=xx2*cs2-y2*sn2;
1923           y2=xx2*sn2+y2*cs2;
1924         }
1925
1926         if(TMath::Abs(z1-z2) > fgkMaxSeedDeltaZ12) continue;
1927
1928         Double_t zz=z1 - z1/x1*(x1-x2);
1929
1930         if (TMath::Abs(zz-z2)>fgkMaxSeedDeltaZ) continue;
1931
1932         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1933         if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
1934
1935         x[0]=y1;
1936         x[1]=z1;
1937         x[4]=f1trd(x1,y1,x2,y2,x3,y3);
1938
1939         if (TMath::Abs(x[4]) > fgkMaxSeedC) continue;
1940
1941         x[2]=f2trd(x1,y1,x2,y2,x3,y3);
1942
1943         if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1944
1945         x[3]=f3trd(x1,y1,x2,y2,z1,z2);
1946
1947         if (TMath::Abs(x[3]) > fgkMaxSeedTan) continue;
1948
1949         Double_t a=asin(x[2]);
1950         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
1951
1952         if (TMath::Abs(zv)>fgkMaxSeedVertexZ) continue;
1953
1954         Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
1955         Double_t sy2=cl->GetSigmaY2(),     sz2=cl->GetSigmaZ2();
1956         Double_t sy3=fgkSeedErrorSY3, sy=fgkSeedErrorSY, sz=fgkSeedErrorSZ;
1957
1958         // Tilt changes
1959         Double_t h01 = GetTiltFactor(r1[is]);
1960         Double_t xuFactor = 100.;
1961         if(fNoTilt) {
1962           h01 = 0;
1963           xuFactor = 1;
1964         }
1965
1966         sy1=sy1+sz1*h01*h01;
1967         Double_t syz=sz1*(-h01);
1968         // end of tilt changes
1969
1970         Double_t f40=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1971         Double_t f42=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1972         Double_t f43=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1973         Double_t f20=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1974         Double_t f22=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1975         Double_t f23=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1976         Double_t f30=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1977         Double_t f31=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1978         Double_t f32=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1979         Double_t f34=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1980
1981
1982         c[0]=sy1;
1983         //        c[1]=0.;       c[2]=sz1;
1984         c[1]=syz;       c[2]=sz1*xuFactor;
1985         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1986         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
1987                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1988         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1989         c[13]=f30*sy1*f40+f32*sy2*f42;
1990         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1991
1992         UInt_t index=r1.GetIndex(is);
1993
1994         AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
1995
1996         Int_t rc=FollowProlongation(*track, i2);
1997
1998         if ((rc < 1) ||
1999             (track->GetNumberOfClusters() <
2000              (outer-inner)*fgkMinClustersInSeed)) delete track;
2001         else {
2002           fSeeds->AddLast(track); fNseeds++;
2003 //          cerr<<"\r found seed "<<fNseeds;
2004         }
2005       }
2006     }
2007   }
2008 }
2009 //__________________________________________________________________________
2010 void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD * esd)
2011 {
2012   //
2013   // Creates  seeds using clusters between  position inner plane  and outer plane
2014   //
2015   const Double_t maxtheta = 1;
2016   const Double_t maxphi   = 2.0;
2017   //
2018   const Double_t kRoad0y  =  6;     // road for middle cluster
2019   const Double_t kRoad0z  =  8.5;   // road for middle cluster
2020   //
2021   const Double_t kRoad1y  =  2;    // road in y for seeded cluster
2022   const Double_t kRoad1z  =  20;    // road in z for seeded cluster
2023   //
2024   const Double_t kRoad2y  =  3;    // road in y for extrapolated cluster
2025   const Double_t kRoad2z  =  20;   // road in z for extrapolated cluster
2026   const Int_t    maxseed  = 3000;
2027   Int_t maxSec=AliTRDgeometry::kNsect;
2028
2029   //
2030   // linear fitters in planes
2031   TLinearFitter fitterTC(2,"hyp2");  // fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
2032   TLinearFitter fitterT2(4,"hyp4");  // fitting with tilting pads - kz not fixed
2033   fitterTC.StoreData(kTRUE);
2034   fitterT2.StoreData(kTRUE);
2035   AliRieman rieman(1000);   // rieman fitter
2036   AliRieman rieman2(1000);   // rieman fitter
2037   //
2038   // find the maximal and minimal layer for the planes
2039   //
2040   Int_t layers[6][2];
2041   AliTRDpropagationLayer* reflayers[6];
2042   for (Int_t i=0;i<6;i++){layers[i][0]=10000; layers[i][1]=0;}
2043   for (Int_t ns=0;ns<maxSec;ns++){
2044     for (Int_t ilayer=0;ilayer<fTrSec[ns]->GetNumberOfLayers();ilayer++){
2045       AliTRDpropagationLayer& layer=*(fTrSec[ns]->GetLayer(ilayer));
2046       if (layer==0) continue;
2047       Int_t det   = layer[0]->GetDetector();
2048       Int_t plane = fGeom->GetPlane(det);
2049       if (ilayer<layers[plane][0]) layers[plane][0] = ilayer;
2050       if (ilayer>layers[plane][1]) layers[plane][1] = ilayer;
2051     }
2052   }
2053   //
2054   AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(0,0);
2055   Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
2056   Double_t hL[6];         // tilting angle
2057   Double_t xcl[6];        // x - position of reference cluster
2058   Double_t ycl[6];        // y - position of reference cluster
2059   Double_t zcl[6];        // z - position of reference cluster
2060   AliTRDcluster *cl[6]={0,0,0,0,0,0};    // seeding clusters
2061   Float_t padlength[6]={10,10,10,10,10,10};   //current pad-length
2062   Double_t chi2R =0, chi2Z=0;
2063   Double_t chi2RF =0, chi2ZF=0;
2064   //
2065   Int_t nclusters;     // total number of clusters
2066   for (Int_t i=0;i<6;i++) {hL[i]=h01; if (i%2==1) hL[i]*=-1.;}
2067   //
2068   //
2069   //         registered seed
2070   AliTRDseed *pseed = new AliTRDseed[maxseed*6];
2071   AliTRDseed *seed[maxseed];
2072   for (Int_t iseed=0;iseed<maxseed;iseed++) seed[iseed]= &pseed[iseed*6];
2073   AliTRDseed *cseed = seed[0];
2074   //
2075   Double_t   seedquality[maxseed];
2076   Double_t   seedquality2[maxseed];
2077   Double_t   seedparams[maxseed][7];
2078   Int_t      seedlayer[maxseed];
2079   Int_t      registered =0;
2080   Int_t      sort[maxseed];
2081   //
2082   // seeding part
2083   //
2084   for (Int_t ns = 0; ns<maxSec; ns++){         //loop over sectors
2085   //for (Int_t ns = 0; ns<5; ns++){         //loop over sectors
2086     registered = 0;   // reset registerd seed counter
2087     cseed      = seed[registered];
2088     Float_t iter=0;
2089     for (Int_t sLayer=2; sLayer>=0;sLayer--){
2090       //for (Int_t dseed=5;dseed<15; dseed+=3){  //loop over central seeding time bins
2091       iter+=1.;
2092       Int_t dseed = 5+Int_t(iter)*3;
2093       // Initialize seeding layers
2094       for (Int_t ilayer=0;ilayer<6;ilayer++){
2095         reflayers[ilayer] = fTrSec[ns]->GetLayer(layers[ilayer][1]-dseed);
2096         xcl[ilayer]       = reflayers[ilayer]->GetX();
2097       }
2098       //
2099       Double_t xref                 = (xcl[sLayer+1] + xcl[sLayer+2])*0.5;
2100       AliTRDpropagationLayer& layer0=*reflayers[sLayer+0];
2101       AliTRDpropagationLayer& layer1=*reflayers[sLayer+1];
2102       AliTRDpropagationLayer& layer2=*reflayers[sLayer+2];
2103       AliTRDpropagationLayer& layer3=*reflayers[sLayer+3];
2104       //
2105       Int_t maxn3  = layer3;
2106       for (Int_t icl3=0;icl3<maxn3;icl3++){
2107         AliTRDcluster *cl3 = layer3[icl3];
2108         if (!cl3) continue;
2109         padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2()*12.);
2110         ycl[sLayer+3] = cl3->GetY();
2111         zcl[sLayer+3] = cl3->GetZ();
2112         Float_t yymin0 = ycl[sLayer+3] - 1- maxphi *(xcl[sLayer+3]-xcl[sLayer+0]);
2113         Float_t yymax0 = ycl[sLayer+3] + 1+ maxphi *(xcl[sLayer+3]-xcl[sLayer+0]);
2114         Int_t   maxn0 = layer0;  //
2115         for (Int_t icl0=layer0.Find(yymin0);icl0<maxn0;icl0++){
2116           AliTRDcluster *cl0 = layer0[icl0];
2117           if (!cl0) continue;
2118           if (cl3->IsUsed()&&cl0->IsUsed()) continue;
2119           ycl[sLayer+0] = cl0->GetY();
2120           zcl[sLayer+0] = cl0->GetZ();
2121           if ( ycl[sLayer+0]>yymax0) break;
2122           Double_t tanphi   = (ycl[sLayer+3]-ycl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
2123           if (TMath::Abs(tanphi)>maxphi) continue;
2124           Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
2125           if (TMath::Abs(tantheta)>maxtheta) continue;
2126           padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2()*12.);
2127           //
2128           // expected position in 1 layer
2129           Double_t y1exp = ycl[sLayer+0]+(tanphi)  *(xcl[sLayer+1]-xcl[sLayer+0]);
2130           Double_t z1exp = zcl[sLayer+0]+(tantheta)*(xcl[sLayer+1]-xcl[sLayer+0]);
2131           Float_t yymin1 = y1exp - kRoad0y-tanphi;
2132           Float_t yymax1 = y1exp + kRoad0y+tanphi;
2133           Int_t   maxn1  = layer1;  //
2134           //
2135           for (Int_t icl1=layer1.Find(yymin1);icl1<maxn1;icl1++){
2136             AliTRDcluster *cl1 = layer1[icl1];
2137             if (!cl1) continue;
2138             Int_t nusedCl = 0;
2139             if (cl3->IsUsed()) nusedCl++;
2140             if (cl0->IsUsed()) nusedCl++;
2141             if (cl1->IsUsed()) nusedCl++;
2142             if (nusedCl>1) continue;
2143             ycl[sLayer+1] = cl1->GetY();
2144             zcl[sLayer+1] = cl1->GetZ();
2145             if ( ycl[sLayer+1]>yymax1) break;
2146             if (TMath::Abs(ycl[sLayer+1]-y1exp)>kRoad0y+tanphi) continue;
2147             if (TMath::Abs(zcl[sLayer+1]-z1exp)>kRoad0z)        continue;
2148             padlength[sLayer+1] = TMath::Sqrt(cl1->GetSigmaZ2()*12.);
2149             //
2150             Double_t y2exp  = ycl[sLayer+0]+(tanphi)  *(xcl[sLayer+2]-xcl[sLayer+0])+(ycl[sLayer+1]-y1exp);
2151             Double_t z2exp  = zcl[sLayer+0]+(tantheta)*(xcl[sLayer+2]-xcl[sLayer+0]);
2152             Int_t    index2 = layer2.FindNearestCluster(y2exp,z2exp,kRoad1y,  kRoad1z);
2153             if (index2<=0) continue;
2154             AliTRDcluster *cl2 = (AliTRDcluster*)GetCluster(index2);
2155             padlength[sLayer+2] = TMath::Sqrt(cl2->GetSigmaZ2()*12.);
2156             ycl[sLayer+2] = cl2->GetY();
2157             zcl[sLayer+2] = cl2->GetZ();
2158             if (TMath::Abs(cl2->GetZ()-z2exp)>kRoad0z)        continue;
2159             //
2160             rieman.Reset();
2161             rieman.AddPoint(xcl[sLayer+0],ycl[sLayer+0],zcl[sLayer+0],1,10);
2162             rieman.AddPoint(xcl[sLayer+1],ycl[sLayer+1],zcl[sLayer+1],1,10);
2163             rieman.AddPoint(xcl[sLayer+3],ycl[sLayer+3],zcl[sLayer+3],1,10);
2164             rieman.AddPoint(xcl[sLayer+2],ycl[sLayer+2],zcl[sLayer+2],1,10);
2165             rieman.Update();
2166             //
2167             // reset fitter
2168             for (Int_t iLayer=0;iLayer<6;iLayer++){
2169               cseed[iLayer].Reset();
2170             }
2171             chi2Z =0.; chi2R=0.;
2172             for (Int_t iLayer=0;iLayer<4;iLayer++){
2173               cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
2174               chi2Z += (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer])*
2175                 (cseed[sLayer+iLayer].fZref[0]- zcl[sLayer+iLayer]);
2176               cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
2177               cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
2178               chi2R += (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer])*
2179                 (cseed[sLayer+iLayer].fYref[0]- ycl[sLayer+iLayer]);
2180               cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
2181             }
2182             if (TMath::Sqrt(chi2R)>1./iter) continue;
2183             if (TMath::Sqrt(chi2Z)>7./iter) continue;
2184             //
2185             //
2186             //
2187             Float_t minmax[2]={-100,100};
2188             for (Int_t iLayer=0;iLayer<4;iLayer++){
2189               Float_t max = zcl[sLayer+iLayer]+padlength[sLayer+iLayer]*0.5+1 -cseed[sLayer+iLayer].fZref[0];
2190               if (max<minmax[1]) minmax[1]=max;
2191               Float_t min = zcl[sLayer+iLayer]-padlength[sLayer+iLayer]*0.5-1 -cseed[sLayer+iLayer].fZref[0];
2192               if (min>minmax[0]) minmax[0]=min;
2193             }
2194             Bool_t isFake = kFALSE;
2195             if (cl0->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
2196             if (cl1->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
2197             if (cl2->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
2198             if ((!isFake) || (icl3%10)==0 ){  //debugging print
2199               TTreeSRedirector& cstream = *fDebugStreamer;
2200               cstream<<"Seeds0"<<
2201                 "isFake="<<isFake<<
2202                 "Cl0.="<<cl0<<
2203                 "Cl1.="<<cl1<<
2204                 "Cl2.="<<cl2<<
2205                 "Cl3.="<<cl3<<
2206                 "Xref="<<xref<<
2207                 "X0="<<xcl[sLayer+0]<<
2208                 "X1="<<xcl[sLayer+1]<<
2209                 "X2="<<xcl[sLayer+2]<<
2210                 "X3="<<xcl[sLayer+3]<<
2211                 "Y2exp="<<y2exp<<
2212                 "Z2exp="<<z2exp<<
2213                 "Chi2R="<<chi2R<<
2214                 "Chi2Z="<<chi2Z<<
2215                 "Seed0.="<<&cseed[sLayer+0]<<
2216                 "Seed1.="<<&cseed[sLayer+1]<<
2217                 "Seed2.="<<&cseed[sLayer+2]<<
2218                 "Seed3.="<<&cseed[sLayer+3]<<
2219                 "Zmin="<<minmax[0]<<
2220                 "Zmax="<<minmax[1]<<
2221                 "\n";
2222             }
2223
2224             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2225             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2226             //<<<<<<<<<<<<<<<<<<    FIT SEEDING PART                  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2227             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2228             cl[sLayer+0] = cl0;
2229             cl[sLayer+1] = cl1;
2230             cl[sLayer+2] = cl2;
2231             cl[sLayer+3] = cl3;
2232             Bool_t isOK=kTRUE;
2233             for (Int_t jLayer=0;jLayer<4;jLayer++){
2234               cseed[sLayer+jLayer].fTilt = hL[sLayer+jLayer];
2235               cseed[sLayer+jLayer].fPadLength = padlength[sLayer+jLayer];
2236               cseed[sLayer+jLayer].fX0   = xcl[sLayer+jLayer];
2237               for (Int_t iter=0; iter<2; iter++){
2238                 //
2239                 // in iteration 0 we try only one pad-row
2240                 // if quality not sufficient we try 2 pad-rows - about 5% of tracks cross 2 pad-rows
2241                 //
2242                 AliTRDseed tseed = cseed[sLayer+jLayer];
2243                 Float_t    roadz  = padlength[sLayer+jLayer]*0.5;
2244                 if (iter>0) roadz = padlength[sLayer+jLayer];
2245                 //
2246                 Float_t quality =10000;
2247                 for (Int_t iTime=2;iTime<20;iTime++){
2248                   AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[sLayer+jLayer][1]-iTime));
2249                   Double_t dxlayer= layer.GetX()-xcl[sLayer+jLayer];
2250                   Double_t zexp   = cl[sLayer+jLayer]->GetZ() ;
2251                   if (iter>0){
2252                     // try 2 pad-rows in second iteration
2253                     zexp  = tseed.fZref[0]+ tseed.fZref[1]*dxlayer;
2254                     if (zexp>cl[sLayer+jLayer]->GetZ()) zexp = cl[sLayer+jLayer]->GetZ()+padlength[sLayer+jLayer]*0.5;
2255                     if (zexp<cl[sLayer+jLayer]->GetZ()) zexp = cl[sLayer+jLayer]->GetZ()-padlength[sLayer+jLayer]*0.5;
2256                   }
2257                   //
2258                   Double_t yexp  =  tseed.fYref[0]+
2259                     tseed.fYref[1]*dxlayer;
2260                   Int_t    index = layer.FindNearestCluster(yexp,zexp,kRoad1y, roadz);
2261                   if (index<=0) continue;
2262                   AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
2263                   //
2264                   tseed.fIndexes[iTime]  = index;
2265                   tseed.fClusters[iTime] = cl;   // register cluster
2266                   tseed.fX[iTime] = dxlayer;     // register cluster
2267                   tseed.fY[iTime] = cl->GetY();  // register cluster
2268                   tseed.fZ[iTime] = cl->GetZ();  // register cluster
2269                 }
2270                 tseed.Update();
2271                 //count the number of clusters and distortions into quality
2272                 Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
2273                 Float_t tquality   = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
2274                   TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
2275                   2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
2276                 if (iter==0 && tseed.isOK()) {
2277                   cseed[sLayer+jLayer] = tseed;
2278                   quality = tquality;
2279                   if (tquality<5) break;
2280                 }
2281                 if (tseed.isOK() && tquality<quality)
2282                   cseed[sLayer+jLayer] = tseed;
2283               }
2284               if (!cseed[sLayer+jLayer].isOK()){
2285                 isOK = kFALSE;
2286                 break;
2287               }
2288               cseed[sLayer+jLayer].CookLabels();
2289               cseed[sLayer+jLayer].UpdateUsed();
2290               nusedCl+= cseed[sLayer+jLayer].fNUsed;
2291               if (nusedCl>25){
2292                 isOK = kFALSE;
2293                 break;
2294               }
2295             }
2296             //
2297             if (!isOK) continue;
2298             nclusters=0;
2299             for (Int_t iLayer=0;iLayer<4;iLayer++){
2300               if (cseed[sLayer+iLayer].isOK()){
2301                 nclusters+=cseed[sLayer+iLayer].fN2;
2302               }
2303             }
2304             //
2305             // iteration 0
2306             rieman.Reset();
2307             for (Int_t iLayer=0;iLayer<4;iLayer++){
2308               rieman.AddPoint(xcl[sLayer+iLayer],cseed[sLayer+iLayer].fYfitR[0],
2309                               cseed[sLayer+iLayer].fZProb,1,10);
2310             }
2311             rieman.Update();
2312             //
2313             //
2314             chi2R =0; chi2Z=0;
2315             for (Int_t iLayer=0;iLayer<4;iLayer++){
2316               cseed[sLayer+iLayer].fYref[0] = rieman.GetYat(xcl[sLayer+iLayer]);
2317               chi2R += (cseed[sLayer+iLayer].fYref[0]-cseed[sLayer+iLayer].fYfitR[0])*
2318                 (cseed[sLayer+iLayer].fYref[0]-cseed[sLayer+iLayer].fYfitR[0]);
2319               cseed[sLayer+iLayer].fYref[1] = rieman.GetDYat(xcl[sLayer+iLayer]);
2320               cseed[sLayer+iLayer].fZref[0] = rieman.GetZat(xcl[sLayer+iLayer]);
2321               chi2Z += (cseed[sLayer+iLayer].fZref[0]- cseed[sLayer+iLayer].fMeanz)*
2322                 (cseed[sLayer+iLayer].fZref[0]- cseed[sLayer+iLayer].fMeanz);
2323               cseed[sLayer+iLayer].fZref[1] = rieman.GetDZat(xcl[sLayer+iLayer]);
2324             }
2325             Double_t curv = rieman.GetC();
2326             //
2327             // likelihoods
2328             //
2329             Double_t sumda =
2330               TMath::Abs(cseed[sLayer+0].fYfitR[1]- cseed[sLayer+0].fYref[1])+
2331               TMath::Abs(cseed[sLayer+1].fYfitR[1]- cseed[sLayer+1].fYref[1])+
2332               TMath::Abs(cseed[sLayer+2].fYfitR[1]- cseed[sLayer+2].fYref[1])+
2333               TMath::Abs(cseed[sLayer+3].fYfitR[1]- cseed[sLayer+3].fYref[1]);
2334             Double_t likea = TMath::Exp(-sumda*10.6);
2335             Double_t likechi2 = 0.0000000001;
2336             if (chi2R<0.5) likechi2+=TMath::Exp(-TMath::Sqrt(chi2R)*7.73);
2337             Double_t likechi2z = TMath::Exp(-chi2Z*0.088)/TMath::Exp(-chi2Z*0.019);
2338             Double_t likeN    = TMath::Exp(-(72-nclusters)*0.19);
2339             Double_t like     = likea*likechi2*likechi2z*likeN;
2340             //
2341             Double_t likePrimY = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fYref[1]-130*curv)*1.9);
2342             Double_t likePrimZ = TMath::Exp(-TMath::Abs(cseed[sLayer+0].fZref[1]-
2343                                                         cseed[sLayer+0].fZref[0]/xcl[sLayer+0])*5.9);
2344             Double_t likePrim  = TMath::Max(likePrimY*likePrimZ,0.0005);
2345
2346             seedquality[registered]  = like;
2347             seedlayer[registered]    = sLayer;
2348             if (TMath::Log(0.000000000000001+like)<-15) continue;
2349             AliTRDseed seedb[6];
2350             for (Int_t iLayer=0;iLayer<6;iLayer++){
2351               seedb[iLayer] = cseed[iLayer];
2352             }
2353             //
2354             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2355             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2356             //<<<<<<<<<<<<<<<   FULL TRACK FIT PART         <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2357             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
2358             //
2359             Int_t nlayers            = 0;
2360             Int_t nusedf             = 0;
2361             Int_t findable           = 0;
2362             //
2363             // add new layers  - avoid long extrapolation
2364             //
2365             Int_t tLayer[2]={0,0};
2366             if (sLayer==2) {tLayer[0]=1; tLayer[1]=0;}
2367             if (sLayer==1) {tLayer[0]=5; tLayer[1]=0;}
2368             if (sLayer==0) {tLayer[0]=4; tLayer[1]=5;}
2369             //
2370             for (Int_t iLayer=0;iLayer<2;iLayer++){
2371               Int_t jLayer = tLayer[iLayer];      // set tracking layer
2372               cseed[jLayer].Reset();
2373               cseed[jLayer].fTilt    = hL[jLayer];
2374               cseed[jLayer].fPadLength = padlength[jLayer];
2375               cseed[jLayer].fX0      = xcl[jLayer];
2376               // get pad length and rough cluster
2377               Int_t    indexdummy = reflayers[jLayer]->FindNearestCluster(cseed[jLayer].fYref[0],
2378                                                                           cseed[jLayer].fZref[0],kRoad2y,kRoad2z);
2379               if (indexdummy<=0) continue;
2380               AliTRDcluster *cldummy = (AliTRDcluster*)GetCluster(indexdummy);
2381               padlength[jLayer]      = TMath::Sqrt(cldummy->GetSigmaZ2()*12.);
2382             }
2383             AliTRDseed::FitRiemanTilt(cseed, kTRUE);
2384             //
2385             for (Int_t iLayer=0;iLayer<2;iLayer++){
2386               Int_t jLayer = tLayer[iLayer];      // set tracking layer
2387               if ( (jLayer==0) && !(cseed[1].isOK())) continue;  // break not allowed
2388               if ( (jLayer==5) && !(cseed[4].isOK())) continue;  // break not allowed
2389               Float_t  zexp  = cseed[jLayer].fZref[0];
2390               Double_t zroad =  padlength[jLayer]*0.5+1.;
2391               //
2392               //
2393               for (Int_t iter=0;iter<2;iter++){
2394                 AliTRDseed tseed = cseed[jLayer];
2395                 Float_t quality = 10000;
2396                 for (Int_t iTime=2;iTime<20;iTime++){
2397                   AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[jLayer][1]-iTime));
2398                   Double_t dxlayer     = layer.GetX()-xcl[jLayer];
2399                   Double_t yexp        = tseed.fYref[0]+tseed.fYref[1]*dxlayer;
2400                   Float_t  yroad       = kRoad1y;
2401                   Int_t    index = layer.FindNearestCluster(yexp,zexp, yroad, zroad);
2402                   if (index<=0) continue;
2403                   AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
2404                   //
2405                   tseed.fIndexes[iTime]  = index;
2406                   tseed.fClusters[iTime] = cl;   // register cluster
2407                   tseed.fX[iTime] = dxlayer;     // register cluster
2408                   tseed.fY[iTime] = cl->GetY();  // register cluster
2409                   tseed.fZ[iTime] = cl->GetZ();  // register cluster
2410                 }
2411                 tseed.Update();
2412                 if (tseed.isOK()){
2413                   Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
2414                   Float_t tquality   = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
2415                     TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
2416                     2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
2417                   //
2418                   if (tquality<quality){
2419                     cseed[jLayer]=tseed;
2420                     quality = tquality;
2421                   }
2422                 }
2423                 zroad*=2.;
2424               }
2425               if ( cseed[jLayer].isOK()){
2426                 cseed[jLayer].CookLabels();
2427                 cseed[jLayer].UpdateUsed();
2428                 nusedf+= cseed[jLayer].fNUsed;
2429                 AliTRDseed::FitRiemanTilt(cseed, kTRUE);
2430               }
2431             }
2432             //
2433             //
2434             // make copy
2435             AliTRDseed bseed[6];
2436             for (Int_t jLayer=0;jLayer<6;jLayer++){
2437               bseed[jLayer] = cseed[jLayer];
2438             }
2439             Float_t lastquality = 10000;
2440             Float_t lastchi2    = 10000;
2441             Float_t chi2        = 1000;
2442
2443             //
2444             for (Int_t iter =0; iter<4;iter++){
2445               //
2446               // sort tracklets according "quality", try to "improve" 4 worst
2447               //
2448               Float_t sumquality = 0;
2449               Float_t squality[6];
2450               Int_t   sortindexes[6];
2451               for (Int_t jLayer=0;jLayer<6;jLayer++){
2452                 if (bseed[jLayer].isOK()){
2453                   AliTRDseed &tseed = bseed[jLayer];
2454                   Double_t zcor  =  tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
2455                   Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
2456                   Float_t tquality  = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
2457                     TMath::Abs(tseed.fYfit[0]-(tseed.fYref[0]-zcor))/0.2+
2458                     2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
2459                   squality[jLayer] = tquality;
2460                 }
2461                 else  squality[jLayer]=-1;
2462                 sumquality +=squality[jLayer];
2463               }
2464
2465               if (sumquality>=lastquality ||  chi2>lastchi2) break;
2466               lastquality = sumquality;
2467               lastchi2    = chi2;
2468               if (iter>0){
2469                 for (Int_t jLayer=0;jLayer<6;jLayer++){
2470                   cseed[jLayer] = bseed[jLayer];
2471                 }
2472               }
2473               TMath::Sort(6,squality,sortindexes,kFALSE);
2474               //
2475               //
2476               for (Int_t jLayer=5;jLayer>1;jLayer--){
2477                 Int_t bLayer = sortindexes[jLayer];
2478                 AliTRDseed tseed = bseed[bLayer];
2479                 for (Int_t iTime=2;iTime<20;iTime++){
2480                   AliTRDpropagationLayer& layer = *(fTrSec[ns]->GetLayer(layers[bLayer][1]-iTime));
2481                   Double_t dxlayer= layer.GetX()-xcl[bLayer];
2482                   //
2483                   Double_t zexp  =  tseed.fZref[0];
2484                   Double_t zcor  =  tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
2485                   //
2486                   Float_t  roadz = padlength[bLayer]+1;
2487                   if (TMath::Abs(tseed.fZProb-zexp)> padlength[bLayer]*0.5) {roadz = padlength[bLayer]*0.5;}
2488                   if (tseed.fZfit[1]*tseed.fZref[1]<0) {roadz = padlength[bLayer]*0.5;}
2489                   if (TMath::Abs(tseed.fZProb-zexp)<0.1*padlength[bLayer]) {
2490                     zexp = tseed.fZProb;
2491                     roadz = padlength[bLayer]*0.5;
2492                   }
2493                   //
2494                   Double_t yexp  =  tseed.fYref[0]+
2495                     tseed.fYref[1]*dxlayer-zcor;
2496                   Int_t    index = layer.FindNearestCluster(yexp,zexp,kRoad1y, roadz);
2497                   if (index<=0) continue;
2498                   AliTRDcluster *cl = (AliTRDcluster*)GetCluster(index);
2499                   //
2500                   tseed.fIndexes[iTime]  = index;
2501                   tseed.fClusters[iTime] = cl;   // register cluster
2502                   tseed.fX[iTime] = dxlayer;     // register cluster
2503                   tseed.fY[iTime] = cl->GetY();  // register cluster
2504                   tseed.fZ[iTime] = cl->GetZ();  // register cluster
2505                 }
2506                 tseed.Update();
2507                 if (tseed.isOK()) {
2508                   Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
2509                   Double_t zcor  =  tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
2510                   //
2511                   Float_t tquality   = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
2512                     TMath::Abs(tseed.fYfit[0]-(tseed.fYref[0]-zcor))/0.2+
2513                     2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
2514                   //
2515                   if (tquality<squality[bLayer])
2516                     bseed[bLayer] = tseed;
2517                 }
2518               }
2519               chi2 = AliTRDseed::FitRiemanTilt(bseed, kTRUE);
2520             }
2521             //
2522             //
2523             //
2524             nclusters  = 0;
2525             nlayers    = 0;
2526             findable   = 0;
2527             for (Int_t iLayer=0;iLayer<6;iLayer++) {
2528               if (TMath::Abs(cseed[iLayer].fYref[0]/cseed[iLayer].fX0)<0.15)
2529                 findable++;
2530               if (cseed[iLayer].isOK()){
2531                 nclusters+=cseed[iLayer].fN2;
2532                 nlayers++;
2533               }
2534             }
2535             if (nlayers<3) continue;
2536             rieman.Reset();
2537             for (Int_t iLayer=0;iLayer<6;iLayer++){
2538               if (cseed[iLayer].isOK()) rieman.AddPoint(xcl[iLayer],cseed[iLayer].fYfitR[0],
2539                                                                    cseed[iLayer].fZProb,1,10);
2540             }
2541             rieman.Update();
2542             //
2543             chi2RF =0;
2544             chi2ZF =0;
2545             for (Int_t iLayer=0;iLayer<6;iLayer++){
2546               if (cseed[iLayer].isOK()){
2547                 cseed[iLayer].fYref[0] = rieman.GetYat(xcl[iLayer]);
2548                 chi2RF += (cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0])*
2549                   (cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0]);
2550                 cseed[iLayer].fYref[1] = rieman.GetDYat(xcl[iLayer]);
2551                 cseed[iLayer].fZref[0] = rieman.GetZat(xcl[iLayer]);
2552                 chi2ZF += (cseed[iLayer].fZref[0]- cseed[iLayer].fMeanz)*
2553                   (cseed[iLayer].fZref[0]- cseed[iLayer].fMeanz);
2554                 cseed[iLayer].fZref[1] = rieman.GetDZat(xcl[iLayer]);
2555               }
2556             }
2557             chi2RF/=TMath::Max((nlayers-3.),1.);
2558             chi2ZF/=TMath::Max((nlayers-3.),1.);
2559             curv = rieman.GetC();
2560
2561             //
2562
2563             Double_t xref2    = (xcl[2]+xcl[3])*0.5;  // middle of the chamber
2564             Double_t dzmf     = rieman.GetDZat(xref2);
2565             Double_t zmf      = rieman.GetZat(xref2);
2566             //
2567             // fit hyperplane
2568             //
2569             Int_t npointsT =0;
2570             fitterTC.ClearPoints();
2571             fitterT2.ClearPoints();
2572             rieman2.Reset();
2573             for (Int_t iLayer=0; iLayer<6;iLayer++){
2574               if (!cseed[iLayer].isOK()) continue;
2575               for (Int_t itime=0;itime<25;itime++){
2576                 if (!cseed[iLayer].fUsable[itime]) continue;
2577                 Double_t x   = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2;  // x relative to the midle chamber
2578                 Double_t y   = cseed[iLayer].fY[itime];
2579                 Double_t z   = cseed[iLayer].fZ[itime];
2580                 // ExB correction to the correction
2581                 // tilted rieman
2582                 //
2583                 Double_t uvt[6];
2584                 Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0;      // global x
2585                 //
2586                 Double_t t = 1./(x2*x2+y*y);
2587                 uvt[1]  = t;    // t
2588                 uvt[0]  = 2.*x2*uvt[1];      // u
2589                 //
2590                 uvt[2]  = 2.0*hL[iLayer]*uvt[1];
2591                 uvt[3]  = 2.0*hL[iLayer]*x*uvt[1];
2592                 uvt[4]  = 2.0*(y+hL[iLayer]*z)*uvt[1];
2593                 //
2594                 Double_t error = 2*0.2*uvt[1];
2595                 fitterT2.AddPoint(uvt,uvt[4],error);
2596                 //
2597                 // constrained rieman
2598                 //
2599                 z =cseed[iLayer].fZ[itime];
2600                 uvt[0]  = 2.*x2*t;           // u
2601                 uvt[1]  = 2*hL[iLayer]*x2*uvt[1];
2602                 uvt[2]  = 2*(y+hL[iLayer]*(z-GetZ()))*t;
2603                 fitterTC.AddPoint(uvt,uvt[2],error);
2604                 //
2605                 rieman2.AddPoint(x2,y,z,1,10);
2606                 npointsT++;
2607               }
2608             }
2609             rieman2.Update();
2610             fitterTC.Eval();
2611             fitterT2.Eval();
2612             Double_t rpolz0 = fitterT2.GetParameter(3);
2613             Double_t rpolz1 = fitterT2.GetParameter(4);
2614             //
2615             // linear fitter  - not possible to make boundaries
2616             // non accept non possible z and dzdx combination
2617             //
2618             Bool_t   acceptablez =kTRUE;
2619             for (Int_t iLayer=0; iLayer<6;iLayer++){
2620               if (cseed[iLayer].isOK()){
2621                 Double_t zT2 =  rpolz0+rpolz1*(xcl[iLayer] - xref2);
2622                 if (TMath::Abs(cseed[iLayer].fZProb-zT2)>padlength[iLayer]*0.5+1)
2623                   acceptablez = kFALSE;
2624               }
2625             }
2626             if (!acceptablez){
2627               fitterT2.FixParameter(3,zmf);
2628               fitterT2.FixParameter(4,dzmf);
2629               fitterT2.Eval();
2630               fitterT2.ReleaseParameter(3);
2631               fitterT2.ReleaseParameter(4);
2632               rpolz0 = fitterT2.GetParameter(3);
2633               rpolz1 = fitterT2.GetParameter(4);
2634             }
2635             //
2636             Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
2637             Double_t chi2TC = fitterTC.GetChisquare()/Float_t(npointsT);
2638             //
2639             Double_t polz1c = fitterTC.GetParameter(2);
2640             Double_t polz0c = polz1c*xref2;
2641             //
2642             Double_t aC     =  fitterTC.GetParameter(0);
2643             Double_t bC     =  fitterTC.GetParameter(1);
2644             Double_t CC     =  aC/TMath::Sqrt(bC*bC+1.);     // curvature
2645             //
2646             Double_t aR     =  fitterT2.GetParameter(0);
2647             Double_t bR     =  fitterT2.GetParameter(1);
2648             Double_t dR     =  fitterT2.GetParameter(2);
2649             Double_t CR     =  1+bR*bR-dR*aR;
2650             Double_t dca    =  0.;
2651             if (CR>0){
2652               dca = -dR/(TMath::Sqrt(1+bR*bR-dR*aR)+TMath::Sqrt(1+bR*bR));
2653               CR  = aR/TMath::Sqrt(CR);
2654             }
2655             //
2656             Double_t chi2ZT2=0, chi2ZTC=0;
2657             for (Int_t iLayer=0; iLayer<6;iLayer++){
2658               if (cseed[iLayer].isOK()){
2659                 Double_t zT2 =  rpolz0+rpolz1*(xcl[iLayer] - xref2);
2660                 Double_t zTC =  polz0c+polz1c*(xcl[iLayer] - xref2);
2661                 chi2ZT2 += TMath::Abs(cseed[iLayer].fMeanz-zT2);
2662                 chi2ZTC += TMath::Abs(cseed[iLayer].fMeanz-zTC);
2663               }
2664             }
2665             chi2ZT2/=TMath::Max((nlayers-3.),1.);
2666             chi2ZTC/=TMath::Max((nlayers-3.),1.);
2667             //
2668             //
2669             //
2670             AliTRDseed::FitRiemanTilt(cseed, kTRUE);
2671             Float_t sumdaf = 0;
2672             for (Int_t iLayer=0;iLayer<6;iLayer++){
2673               if (cseed[iLayer].isOK())
2674                 sumdaf += TMath::Abs((cseed[iLayer].fYfit[1]-cseed[iLayer].fYref[1])/cseed[iLayer].fSigmaY2);
2675             }
2676             sumdaf /= Float_t (nlayers-2.);
2677             //
2678             // likelihoods for full track
2679             //
2680             Double_t likezf      = TMath::Exp(-chi2ZF*0.14);
2681             Double_t likechi2C   = TMath::Exp(-chi2TC*0.677);
2682             Double_t likechi2TR  = TMath::Exp(-chi2TR*0.78);
2683             Double_t likeaf      = TMath::Exp(-sumdaf*3.23);
2684             seedquality2[registered] = likezf*likechi2TR*likeaf;
2685 //          Bool_t isGold = kFALSE;
2686 //
2687 //          if (nlayers == 6        && TMath::Log(0.000000001+seedquality2[index])<-5.) isGold =kTRUE;   // gold
2688 //          if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE;   // gold
2689 //          if (isGold &&nusedf<10){
2690 //            for (Int_t jLayer=0;jLayer<6;jLayer++){
2691 //              if ( seed[index][jLayer].isOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
2692 //                seed[index][jLayer].UseClusters();  //sign gold
2693 //            }
2694 //          }
2695             //
2696             //
2697             //
2698             Int_t index0=0;
2699             if (!cseed[0].isOK()){
2700               index0 = 1;
2701               if (!cseed[1].isOK()) index0 = 2;
2702             }
2703             seedparams[registered][0] = cseed[index0].fX0;
2704             seedparams[registered][1] = cseed[index0].fYref[0];
2705             seedparams[registered][2] = cseed[index0].fZref[0];
2706             seedparams[registered][5] = CR;
2707             seedparams[registered][3] = cseed[index0].fX0*CR - TMath::Sin(TMath::ATan(cseed[0].fYref[1]));
2708             seedparams[registered][4] = cseed[index0].fZref[1]/
2709               TMath::Sqrt(1+cseed[index0].fYref[1]*cseed[index0].fYref[1]);
2710             seedparams[registered][6] = ns;
2711             //
2712             //
2713             Int_t labels[12], outlab[24];
2714             Int_t nlab=0;
2715             for (Int_t iLayer=0;iLayer<6;iLayer++){
2716               if (!cseed[iLayer].isOK()) continue;
2717               if (cseed[iLayer].fLabels[0]>=0) {
2718                 labels[nlab] = cseed[iLayer].fLabels[0];
2719                 nlab++;
2720               }
2721               if (cseed[iLayer].fLabels[1]>=0) {
2722                 labels[nlab] = cseed[iLayer].fLabels[1];
2723                 nlab++;
2724               }
2725             }
2726             Freq(nlab,labels,outlab,kFALSE);
2727             Int_t label = outlab[0];
2728             Int_t frequency  = outlab[1];
2729             for (Int_t iLayer=0;iLayer<6;iLayer++){
2730               cseed[iLayer].fFreq  = frequency;
2731               cseed[iLayer].fC     = CR;
2732               cseed[iLayer].fCC     = CC;
2733               cseed[iLayer].fChi2  = chi2TR;
2734               cseed[iLayer].fChi2Z = chi2ZF;
2735             }
2736             //
2737             if (1||(!isFake)){  //debugging print
2738               Float_t zvertex = GetZ();
2739               TTreeSRedirector& cstream = *fDebugStreamer;
2740               cstream<<"Seeds1"<<
2741                 "isFake="<<isFake<<
2742                 "Vertex="<<zvertex<<
2743                 "Rieman2.="<<&rieman2<<
2744                 "Rieman.="<<&rieman<<
2745                 "Xref="<<xref<<
2746                 "X0="<<xcl[0]<<
2747                 "X1="<<xcl[1]<<
2748                 "X2="<<xcl[2]<<
2749                 "X3="<<xcl[3]<<
2750                 "X4="<<xcl[4]<<
2751                 "X5="<<xcl[5]<<
2752                 "Chi2R="<<chi2R<<
2753                 "Chi2Z="<<chi2Z<<
2754                 "Chi2RF="<<chi2RF<<                          //chi2 of trackletes on full track
2755                 "Chi2ZF="<<chi2ZF<<                          //chi2 z on tracklets on full track
2756                 "Chi2ZT2="<<chi2ZT2<<                        //chi2 z on tracklets on full track  - rieman tilt
2757                 "Chi2ZTC="<<chi2ZTC<<                        //chi2 z on tracklets on full track  - rieman tilt const
2758                 //
2759                 "Chi2TR="<<chi2TR<<                           //chi2 without vertex constrain
2760                 "Chi2TC="<<chi2TC<<                           //chi2 with    vertex constrain
2761                 "C="<<curv<<                                  // non constrained - no tilt correction
2762                 "DR="<<dR<<                                   // DR parameter          - tilt correction
2763                 "DCA="<<dca<<                                 // DCA                   - tilt correction
2764                 "CR="<<CR<<                                   // non constrained curvature - tilt correction
2765                 "CC="<<CC<<                                   // constrained curvature
2766                 "Polz0="<<polz0c<<
2767                 "Polz1="<<polz1c<<
2768                 "RPolz0="<<rpolz0<<
2769                 "RPolz1="<<rpolz1<<
2770                 "Ncl="<<nclusters<<
2771                 "Nlayers="<<nlayers<<
2772                 "NUsedS="<<nusedCl<<
2773                 "NUsed="<<nusedf<<
2774                 "Findable="<<findable<<
2775                 "Like="<<like<<
2776                 "LikePrim="<<likePrim<<
2777                 "Likechi2C="<<likechi2C<<
2778                 "Likechi2TR="<<likechi2TR<<
2779                 "Likezf="<<likezf<<
2780                 "LikeF="<<seedquality2[registered]<<
2781                 "S0.="<<&cseed[0]<<
2782                 "S1.="<<&cseed[1]<<
2783                 "S2.="<<&cseed[2]<<
2784                 "S3.="<<&cseed[3]<<
2785                 "S4.="<<&cseed[4]<<
2786                 "S5.="<<&cseed[5]<<
2787                 "SB0.="<<&seedb[0]<<
2788                 "SB1.="<<&seedb[1]<<
2789                 "SB2.="<<&seedb[2]<<
2790                 "SB3.="<<&seedb[3]<<
2791                 "SB4.="<<&seedb[4]<<
2792                 "SB5.="<<&seedb[5]<<
2793                 "Label="<<label<<
2794                 "Freq="<<frequency<<
2795                 "sLayer="<<sLayer<<
2796                 "\n";
2797             }
2798             if (registered<maxseed-1) {
2799               registered++;
2800               cseed = seed[registered];
2801             }
2802           }// end of loop over layer 1
2803         }  // end of loop over layer 0
2804       }    // end of loop over layer 3
2805     }      // end of loop over seeding time bins
2806     //
2807     // choos best
2808     //
2809     TMath::Sort(registered,seedquality2,sort,kTRUE);
2810     Bool_t signedseed[maxseed];
2811     for (Int_t i=0;i<registered;i++){
2812       signedseed[i]= kFALSE;
2813     }
2814     for (Int_t iter=0; iter<5; iter++){
2815       for (Int_t iseed=0;iseed<registered;iseed++){
2816         Int_t index = sort[iseed];
2817         if (signedseed[index]) continue;
2818         Int_t labelsall[1000];
2819         Int_t nlabelsall=0;
2820         Int_t naccepted=0;;
2821         Int_t sLayer = seedlayer[index];
2822         Int_t ncl   = 0;
2823         Int_t nused = 0;
2824         Int_t nlayers =0;
2825         Int_t findable   = 0;
2826         for (Int_t jLayer=0;jLayer<6;jLayer++){
2827           if (TMath::Abs(seed[index][jLayer].fYref[0]/xcl[jLayer])<0.15)
2828             findable++;
2829           if (seed[index][jLayer].isOK()){
2830             seed[index][jLayer].UpdateUsed();
2831             ncl   +=seed[index][jLayer].fN2;
2832             nused +=seed[index][jLayer].fNUsed;
2833             nlayers++;
2834             //cooking label
2835             for (Int_t itime=0;itime<25;itime++){
2836               if (seed[index][jLayer].fUsable[itime]){
2837                 naccepted++;
2838                 for (Int_t ilab=0;ilab<3;ilab++){
2839                   Int_t tindex = seed[index][jLayer].fClusters[itime]->GetLabel(ilab);
2840                   if (tindex>=0){
2841                     labelsall[nlabelsall] = tindex;
2842                     nlabelsall++;
2843                   }
2844                 }
2845               }
2846             }
2847           }
2848         }
2849         //
2850         if (nused>30) continue;
2851         //
2852         if (iter==0){
2853           if (nlayers<6) continue;
2854           if (TMath::Log(0.000000001+seedquality2[index])<-5.) continue;   // gold
2855         }
2856         //
2857         if (iter==1){
2858           if (nlayers<findable) continue;
2859           if (TMath::Log(0.000000001+seedquality2[index])<-4.) continue;  //
2860         }
2861         //
2862         //
2863         if (iter==2){
2864           if (nlayers==findable || nlayers==6) continue;
2865           if (TMath::Log(0.000000001+seedquality2[index])<-6.) continue;
2866         }
2867         //
2868         if (iter==3){
2869           if (TMath::Log(0.000000001+seedquality2[index])<-5.) continue;
2870         }
2871         //
2872         if (iter==4){
2873           if (TMath::Log(0.000000001+seedquality2[index])-nused/(nlayers-3.)<-15.) continue;
2874         }
2875         //
2876         signedseed[index] = kTRUE;
2877         //
2878         Int_t labels[1000], outlab[1000];
2879         Int_t nlab=0;
2880         for (Int_t iLayer=0;iLayer<6;iLayer++){
2881           if (seed[index][iLayer].isOK()){
2882             if (seed[index][iLayer].fLabels[0]>=0) {
2883               labels[nlab] = seed[index][iLayer].fLabels[0];
2884               nlab++;
2885             }
2886             if (seed[index][iLayer].fLabels[1]>=0) {
2887               labels[nlab] = seed[index][iLayer].fLabels[1];
2888               nlab++;
2889             }
2890           }
2891         }
2892         Freq(nlab,labels,outlab,kFALSE);
2893         Int_t label  = outlab[0];
2894         Int_t frequency  = outlab[1];
2895         Freq(nlabelsall,labelsall,outlab,kFALSE);
2896         Int_t label1 = outlab[0];
2897         Int_t label2 = outlab[2];
2898         Float_t fakeratio = (naccepted-outlab[1])/Float_t(naccepted);
2899         Float_t ratio = Float_t(nused)/Float_t(ncl);
2900         if (ratio<0.25){
2901           for (Int_t jLayer=0;jLayer<6;jLayer++){
2902             if ( seed[index][jLayer].isOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.2 )
2903               seed[index][jLayer].UseClusters();  //sign gold
2904           }
2905         }
2906         //
2907         Int_t eventNr = esd->GetEventNumber();
2908         TTreeSRedirector& cstream = *fDebugStreamer;
2909         //
2910         // register seed
2911         //
2912         AliTRDtrack * track = RegisterSeed(seed[index],seedparams[index]);
2913         AliTRDtrack dummy;
2914         if (!track) track=&dummy;
2915         else{
2916           AliESDtrack esdtrack;
2917           esdtrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
2918           esdtrack.SetLabel(label);
2919           esd->AddTrack(&esdtrack);
2920           TTreeSRedirector& cstream = *fDebugStreamer;
2921           cstream<<"Tracks"<<
2922             "EventNr="<<eventNr<<
2923             "ESD.="<<&esdtrack<<
2924             "trd.="<<track<<
2925             "trdback.="<<track<<
2926             "\n";
2927         }
2928
2929         cstream<<"Seeds2"<<
2930           "Iter="<<iter<<
2931           "Track.="<<track<<
2932           "Like="<<seedquality[index]<<
2933           "LikeF="<<seedquality2[index]<<
2934           "S0.="<<&seed[index][0]<<
2935           "S1.="<<&seed[index][1]<<
2936           "S2.="<<&seed[index][2]<<
2937           "S3.="<<&seed[index][3]<<
2938           "S4.="<<&seed[index][4]<<
2939           "S5.="<<&seed[index][5]<<
2940           "Label="<<label<<
2941           "Label1="<<label1<<
2942           "Label2="<<label2<<
2943           "FakeRatio="<<fakeratio<<
2944           "Freq="<<frequency<<
2945           "Ncl="<<ncl<<
2946           "Nlayers="<<nlayers<<
2947           "Findable="<<findable<<
2948           "NUsed="<<nused<<
2949           "sLayer="<<sLayer<<
2950           "EventNr="<<eventNr<<
2951           "\n";
2952       }
2953     }
2954   }        // end of loop over sectors
2955   delete [] pseed;
2956 }
2957
2958 //_____________________________________________________________________________
2959 Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *ClusterTree) const
2960 {
2961   //
2962   // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
2963   // from the file. The names of the cluster tree and branches
2964   // should match the ones used in AliTRDclusterizer::WriteClusters()
2965   //
2966   Int_t nsize = Int_t(ClusterTree->GetTotBytes()/(sizeof(AliTRDcluster)));
2967   TObjArray *clusterArray = new TObjArray(nsize+1000);
2968
2969   TBranch *branch=ClusterTree->GetBranch("TRDcluster");
2970   if (!branch) {
2971     Error("ReadClusters","Can't get the branch !");
2972     return 1;
2973   }
2974   branch->SetAddress(&clusterArray);
2975
2976   Int_t nEntries = (Int_t) ClusterTree->GetEntries();
2977   //  printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
2978
2979   // Loop through all entries in the tree
2980   Int_t nbytes = 0;
2981   AliTRDcluster *c = 0;
2982   //  printf("\n");
2983   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
2984
2985     // Import the tree
2986     nbytes += ClusterTree->GetEvent(iEntry);
2987
2988     // Get the number of points in the detector
2989     Int_t nCluster = clusterArray->GetEntriesFast();
2990 //    printf("\r Read %d clusters from entry %d", nCluster, iEntry);
2991
2992     // Loop through all TRD digits
2993     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
2994       c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
2995 //       if (c->GetNPads()>3&&(iCluster%3>0)) {
2996 //      delete clusterArray->RemoveAt(iCluster);
2997 //      continue;
2998 //       }
2999       //      AliTRDcluster *co = new AliTRDcluster(*c);  //remove unnecesary coping - + clusters are together in memory
3000       AliTRDcluster *co = c;
3001       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
3002       Int_t ltb = co->GetLocalTimeBin();
3003       if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
3004       else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
3005       array->AddLast(co);
3006       //      delete clusterArray->RemoveAt(iCluster);
3007       clusterArray->RemoveAt(iCluster);
3008     }
3009   }
3010 //   cout<<"Allocated"<<nsize<<"\tLoaded"<<array->GetEntriesFast()<<"\n";
3011
3012   delete clusterArray;
3013
3014   return 0;
3015 }
3016
3017 //__________________________________________________________________
3018 Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint& p) const
3019 {
3020   //
3021   // Get track space point with index i
3022   // Origin: C.Cheshkov
3023   //
3024
3025   AliTRDcluster *cl = (AliTRDcluster*)fClusters->UncheckedAt(index);
3026   Int_t  idet = cl->GetDetector();
3027   Int_t  isector = fGeom->GetSector(idet);
3028   Int_t  ichamber= fGeom->GetChamber(idet);
3029   Int_t  iplan   = fGeom->GetPlane(idet);
3030   Double_t local[3];
3031   local[0]=GetX(isector,iplan,cl->GetLocalTimeBin());
3032   local[1]=cl->GetY();
3033   local[2]=cl->GetZ();
3034   Double_t global[3];
3035   fGeom->RotateBack(idet,local,global);
3036   p.SetXYZ(global[0],global[1],global[2]);
3037   AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
3038   switch (iplan) {
3039   case 0:
3040     iLayer = AliAlignObj::kTRD1;
3041     break;
3042   case 1:
3043     iLayer = AliAlignObj::kTRD2;
3044     break;
3045   case 2:
3046     iLayer = AliAlignObj::kTRD3;
3047     break;
3048   case 3:
3049     iLayer = AliAlignObj::kTRD4;
3050     break;
3051   case 4:
3052     iLayer = AliAlignObj::kTRD5;
3053     break;
3054   case 5:
3055     iLayer = AliAlignObj::kTRD6;
3056     break;
3057   };
3058   Int_t modId = isector*fGeom->Ncham()+ichamber;
3059   UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
3060   p.SetVolumeID(volid);
3061
3062   return kTRUE;
3063
3064 }
3065
3066 //__________________________________________________________________
3067 void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const
3068 {
3069   //
3070   // This cooks a label. Mmmmh, smells good...
3071   //
3072
3073   Int_t label=123456789, index, i, j;
3074   Int_t ncl=pt->GetNumberOfClusters();
3075   const Int_t kRange = fTrSec[0]->GetOuterTimeBin()+1;
3076
3077   Bool_t labelAdded;
3078
3079   //  Int_t s[kRange][2];
3080   Int_t **s = new Int_t* [kRange];
3081   for (i=0; i<kRange; i++) {
3082     s[i] = new Int_t[2];
3083   }
3084   for (i=0; i<kRange; i++) {
3085     s[i][0]=-1;
3086     s[i][1]=0;
3087   }
3088
3089   Int_t t0,t1,t2;
3090   for (i=0; i<ncl; i++) {
3091     index=pt->GetClusterIndex(i);
3092     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
3093     t0=c->GetLabel(0);
3094     t1=c->GetLabel(1);
3095     t2=c->GetLabel(2);
3096   }
3097
3098   for (i=0; i<ncl; i++) {
3099     index=pt->GetClusterIndex(i);
3100     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
3101     for (Int_t k=0; k<3; k++) {
3102       label=c->GetLabel(k);
3103       labelAdded=kFALSE; j=0;
3104       if (label >= 0) {
3105         while ( (!labelAdded) && ( j < kRange ) ) {
3106           if (s[j][0]==label || s[j][1]==0) {
3107             s[j][0]=label;
3108             s[j][1]=s[j][1]+1;
3109             labelAdded=kTRUE;
3110           }
3111           j++;
3112         }
3113       }
3114     }
3115   }
3116
3117   Int_t max=0;
3118   label = -123456789;
3119
3120   for (i=0; i<kRange; i++) {
3121     if (s[i][1]>max) {
3122       max=s[i][1]; label=s[i][0];
3123     }
3124   }
3125
3126   for (i=0; i<kRange; i++) {
3127     delete []s[i];
3128   }
3129
3130   delete []s;
3131
3132   if ((1.- Float_t(max)/ncl) > wrong) label=-label;
3133
3134   pt->SetLabel(label);
3135
3136 }
3137
3138
3139 //__________________________________________________________________
3140 void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const
3141 {
3142   //
3143   // Use clusters, but don't abuse them!
3144   //
3145   const Float_t kmaxchi2 =18;
3146   const Float_t kmincl   =10;
3147   AliTRDtrack * track  = (AliTRDtrack*)t;
3148   //
3149   Int_t ncl=t->GetNumberOfClusters();
3150   for (Int_t i=from; i<ncl; i++) {
3151     Int_t index = t->GetClusterIndex(i);
3152     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
3153     //
3154     Int_t iplane = fGeom->GetPlane(c->GetDetector());
3155     if (track->fTracklets[iplane].GetChi2()>kmaxchi2) continue;
3156     if (track->fTracklets[iplane].GetN()<kmincl) continue;
3157     if (!(c->IsUsed())) c->Use();
3158   }
3159 }
3160
3161
3162 //_____________________________________________________________________
3163 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
3164 {
3165   // Parametrised "expected" error of the cluster reconstruction in Y
3166
3167   Double_t s = 0.08 * 0.08;
3168   return s;
3169 }
3170
3171 //_____________________________________________________________________
3172 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
3173 {
3174   // Parametrised "expected" error of the cluster reconstruction in Z
3175
3176   Double_t s = 9 * 9 /12.;
3177   return s;
3178 }
3179
3180 //_____________________________________________________________________
3181 Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
3182 {
3183   //
3184   // Returns radial position which corresponds to time bin <localTB>
3185   // in tracking sector <sector> and plane <plane>
3186   //
3187
3188   Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, localTB);
3189   Int_t pl = fTrSec[sector]->GetLayerNumber(index);
3190   return fTrSec[sector]->GetLayer(pl)->GetX();
3191
3192 }
3193
3194
3195 //_______________________________________________________
3196 AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x,
3197                Double_t dx, Double_t rho, Double_t radLength, Int_t tbIndex)
3198 {
3199   //
3200   // AliTRDpropagationLayer constructor
3201   //
3202
3203   fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = radLength;
3204   fClusters = NULL; fIndex = NULL; fTimeBinIndex = tbIndex;
3205
3206
3207   for(Int_t i=0; i < (Int_t) kZones; i++) {
3208     fZc[i]=0; fZmax[i] = 0;
3209   }
3210
3211   fYmax = 0;
3212
3213   if(fTimeBinIndex >= 0) {
3214     fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
3215     fIndex = new UInt_t[kMaxClusterPerTimeBin];
3216   }
3217
3218   for (Int_t i=0;i<5;i++) fIsHole[i] = kFALSE;
3219   fHole = kFALSE;
3220   fHoleZc = 0;
3221   fHoleZmax = 0;
3222   fHoleYc = 0;
3223   fHoleYmax = 0;
3224   fHoleRho = 0;
3225   fHoleX0 = 0;
3226
3227 }
3228
3229 //_______________________________________________________
3230 void AliTRDtracker::AliTRDpropagationLayer::SetHole(
3231           Double_t Zmax, Double_t Ymax, Double_t rho,
3232           Double_t radLength, Double_t Yc, Double_t Zc)
3233 {
3234   //
3235   // Sets hole in the layer
3236   //
3237   fHole = kTRUE;
3238   fHoleZc = Zc;
3239   fHoleZmax = Zmax;
3240   fHoleYc = Yc;
3241   fHoleYmax = Ymax;
3242   fHoleRho = rho;
3243   fHoleX0 = radLength;
3244 }
3245
3246
3247 //_______________________________________________________
3248 AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, AliTRDparameter* par)
3249 {
3250   //
3251   // AliTRDtrackingSector Constructor
3252   //
3253   AliTRDpadPlane *padPlane = 0;
3254
3255   fGeom = geo;
3256   fPar = par;
3257   fGeomSector = gs;
3258   fTzeroShift = 0.13;
3259   fN = 0;
3260   //
3261   // get holes description from geometry
3262   Bool_t holes[AliTRDgeometry::kNcham];
3263   //printf("sector\t%d\t",gs);
3264   for (Int_t icham=0; icham<AliTRDgeometry::kNcham;icham++){
3265     holes[icham] = fGeom->IsHole(0,icham,gs);
3266     //printf("%d",holes[icham]);
3267   }
3268   //printf("\n");
3269
3270   for(UInt_t i=0; i < kMaxTimeBinIndex; i++) fTimeBinIndex[i] = -1;
3271
3272
3273   AliTRDpropagationLayer* ppl;
3274
3275   Double_t x, xin, xout, dx, rho, radLength;
3276   Int_t    steps;
3277
3278   // set time bins in the gas of the TPC
3279
3280   xin = 246.055; xout = 254.055; steps = 20; dx = (xout-xin)/steps;
3281   rho = 0.9e-3;  radLength = 28.94;
3282
3283   for(Int_t i=0; i<steps; i++) {
3284     x = xin + i*dx + dx/2;
3285     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
3286     InsertLayer(ppl);
3287   }
3288
3289   // set time bins in the outer field cage vessel
3290
3291   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
3292   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3293   InsertLayer(ppl);
3294
3295   dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
3296   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3297   InsertLayer(ppl);
3298
3299   dx = 2.; xin = xout; xout = xin + dx; rho = 1.45*0.02; radLength = 41.28; // Nomex
3300   steps = 5; dx = (xout - xin)/steps;
3301   for(Int_t i=0; i<steps; i++) {
3302     x = xin + i*dx + dx/2;
3303     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
3304     InsertLayer(ppl);
3305   }
3306
3307   dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
3308   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3309   InsertLayer(ppl);
3310
3311   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
3312   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3313   InsertLayer(ppl);
3314
3315
3316   // set time bins in CO2
3317
3318   xin = xout; xout = 275.0;
3319   steps = 50; dx = (xout - xin)/steps;
3320   rho = 1.977e-3;  radLength = 36.2;
3321
3322   for(Int_t i=0; i<steps; i++) {
3323     x = xin + i*dx + dx/2;
3324     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
3325     InsertLayer(ppl);
3326   }
3327
3328   // set time bins in the outer containment vessel
3329
3330   dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; radLength = 24.01; // Al
3331   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3332   InsertLayer(ppl);
3333
3334   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
3335   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3336   InsertLayer(ppl);
3337
3338   dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
3339   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3340   InsertLayer(ppl);
3341
3342   dx = 3.; xin = xout; xout = xin + dx; rho = 1.45*0.02; radLength = 41.28; // Nomex
3343   steps = 10; dx = (xout - xin)/steps;
3344   for(Int_t i=0; i<steps; i++) {
3345     x = xin + i*dx + dx/2;
3346     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
3347     InsertLayer(ppl);
3348   }
3349
3350   dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; radLength = 44.86; // prepreg
3351   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3352   InsertLayer(ppl);
3353
3354   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; radLength = 44.77; // Tedlar
3355   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3356   InsertLayer(ppl);
3357
3358   dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; radLength = 24.01; // Al
3359   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3360   InsertLayer(ppl);
3361
3362   Double_t xtrd = (Double_t) fGeom->Rmin();
3363
3364   // add layers between TPC and TRD (Air temporarily)
3365   xin = xout; xout = xtrd;
3366   steps = 50; dx = (xout - xin)/steps;
3367   rho = 1.2e-3;  radLength = 36.66;
3368
3369   for(Int_t i=0; i<steps; i++) {
3370     x = xin + i*dx + dx/2;
3371     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
3372     InsertLayer(ppl);
3373   }
3374
3375
3376   //  Double_t alpha=AliTRDgeometry::GetAlpha();
3377
3378   // add layers for each of the planes
3379
3380   Double_t dxRo = (Double_t) fGeom->CroHght();    // Rohacell
3381   Double_t dxSpace = (Double_t) fGeom->Cspace();  // Spacing between planes
3382   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
3383   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
3384   Double_t dxRad = (Double_t) fGeom->CraHght();   // Radiator
3385   Double_t dxTEC = dxRad + dxDrift + dxAmp + dxRo;
3386   Double_t dxPlane = dxTEC + dxSpace;
3387
3388   Int_t tb, tbIndex;
3389   const Int_t  kNchambers = AliTRDgeometry::Ncham();
3390   Double_t  ymax = 0;
3391   //, holeYmax = 0;
3392   Double_t ymaxsensitive=0;
3393   Double_t *zc = new Double_t[kNchambers];
3394   Double_t *zmax = new Double_t[kNchambers];
3395   Double_t *zmaxsensitive = new Double_t[kNchambers];
3396   //  Double_t  holeZmax = 1000.;   // the whole sector is missing
3397
3398   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
3399   if (!commonParam)
3400   {
3401     printf("<AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector> ");
3402     printf("Could not get common params\n");
3403     return;
3404   }
3405
3406   for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
3407     //
3408     // Radiator
3409     xin = xtrd + plane * dxPlane; xout = xin + dxRad;
3410     steps = 12; dx = (xout - xin)/steps; rho = 0.074; radLength = 40.6;
3411     for(Int_t i=0; i<steps; i++) {
3412       x = xin + i*dx + dx/2;
3413       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
3414       InsertLayer(ppl);
3415     }
3416
3417     ymax          = fGeom->GetChamberWidth(plane)/2.;
3418     // Modidified for new pad plane class, 22.04.05 (C.B.)
3419     // ymaxsensitive = (fPar->GetColPadSize(plane)*fPar->GetColMax(plane)-4)/2.;
3420     padPlane = commonParam->GetPadPlane(plane,0);
3421     ymaxsensitive = (padPlane->GetColSize(1)*padPlane->GetNcols()-4)/2.;
3422
3423     //    ymaxsensitive = (fPar->GetColPadSize(plane)*fPar->GetColMax(plane)-4)/2.;
3424
3425     for(Int_t ch = 0; ch < kNchambers; ch++) {
3426       zmax[ch] = fGeom->GetChamberLength(plane,ch)/2;
3427       //
3428       // Modidified for new pad plane class, 22.04.05 (C.B.)
3429       //Float_t pad = fPar->GetRowPadSize(plane,ch,0);
3430       Float_t pad = padPlane->GetRowSize(1);
3431       //Float_t pad = fPar->GetRowPadSize(plane,ch,0);
3432       Float_t row0 = commonParam->GetRow0(plane,ch,0);
3433       Int_t nPads = commonParam->GetRowMax(plane,ch,0);
3434       zmaxsensitive[ch] = Float_t(nPads)*pad/2.;
3435       //      zc[ch] = (pad * nPads)/2 + row0 - pad/2;
3436       //      zc[ch] = (pad * nPads)/2 + row0;
3437       zc[ch] = -(pad * nPads)/2 + row0;
3438       //zc[ch] = row0+zmax[ch]-AliTRDgeometry::RpadW();
3439
3440     }
3441
3442     dx  = fgkDriftCorrection*fPar->GetDriftVelocity()
3443         / AliTRDcalibDB::Instance()->GetSamplingFrequency();
3444     rho = 0.00295 * 0.85; radLength = 11.0;
3445
3446     Double_t x0 = (Double_t) AliTRDgeometry::GetTime0(plane);
3447     Double_t xbottom = x0 - dxDrift;
3448     Double_t xtop = x0 + dxAmp;
3449     //
3450     // Amplification region
3451     steps = (Int_t) (dxAmp/dx);
3452
3453     for(tb = 0; tb < steps; tb++) {
3454       x = x0 + tb * dx + dx/2+ fgkOffsetX;
3455       tbIndex = CookTimeBinIndex(plane, -tb-1);
3456       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
3457       ppl->SetYmax(ymax,ymaxsensitive);
3458       ppl->SetZ(zc, zmax, zmaxsensitive);
3459       ppl->SetHoles(holes);
3460       InsertLayer(ppl);
3461     }
3462     tbIndex = CookTimeBinIndex(plane, -steps);
3463     x = (x + dx/2 + xtop)/2;
3464     dx = 2*(xtop-x);
3465     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
3466     ppl->SetYmax(ymax,ymaxsensitive);
3467     ppl->SetZ(zc, zmax,zmaxsensitive);
3468     ppl->SetHoles(holes);
3469     InsertLayer(ppl);
3470
3471     // Drift region
3472
3473     dx = fgkDriftCorrection*fPar->GetDriftVelocity()
3474        / AliTRDcalibDB::Instance()->GetSamplingFrequency();
3475     steps = (Int_t) (dxDrift/dx)+3;
3476
3477     for(tb = 0; tb < steps; tb++) {
3478       x = x0 - tb * dx - dx/2 + fgkOffsetX;                      //temporary fix - fix it the parameters
3479       tbIndex = CookTimeBinIndex(plane, tb);
3480
3481       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
3482       ppl->SetYmax(ymax,ymaxsensitive);
3483       ppl->SetZ(zc, zmax, zmaxsensitive);
3484       ppl->SetHoles(holes);
3485       InsertLayer(ppl);
3486     }
3487     tbIndex = CookTimeBinIndex(plane, steps);
3488     x = (x - dx/2 + xbottom)/2;
3489     dx = 2*(x-xbottom);
3490     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
3491     ppl->SetYmax(ymax,ymaxsensitive);
3492     ppl->SetZ(zc, zmax, zmaxsensitive);
3493     ppl->SetHoles(holes);
3494     InsertLayer(ppl);
3495
3496     // Pad Plane
3497     xin = xtop; dx = 0.025; xout = xin + dx; rho = 1.7; radLength = 33.0;
3498     ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,radLength,-1);
3499     ppl->SetYmax(ymax,ymaxsensitive);
3500     ppl->SetZ(zc, zmax,zmax);
3501     ppl->SetHoles(holes);
3502     InsertLayer(ppl);
3503
3504     // Rohacell
3505     xin = xout; xout = xtrd + (plane + 1) * dxPlane - dxSpace;
3506     steps = 5; dx = (xout - xin)/steps; rho = 0.074; radLength = 40.6;
3507     for(Int_t i=0; i<steps; i++) {
3508       x = xin + i*dx + dx/2;
3509       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
3510       ppl->SetYmax(ymax,ymaxsensitive);
3511       ppl->SetZ(zc, zmax,zmax);
3512       ppl->SetHoles(holes);
3513       InsertLayer(ppl);
3514     }
3515
3516     // Space between the chambers, air
3517     xin = xout; xout = xtrd + (plane + 1) * dxPlane;
3518     steps = 5; dx = (xout - xin)/steps; rho = 1.29e-3; radLength = 36.66;
3519     for(Int_t i=0; i<steps; i++) {
3520       x = xin + i*dx + dx/2;
3521       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
3522       InsertLayer(ppl);
3523     }
3524   }
3525
3526   // Space between the TRD and RICH
3527   Double_t xRICH = 500.;
3528   xin = xout; xout = xRICH;
3529   steps = 200; dx = (xout - xin)/steps; rho = 1.29e-3; radLength = 36.66;
3530   for(Int_t i=0; i<steps; i++) {
3531     x = xin + i*dx + dx/2;
3532     ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,-1);
3533     InsertLayer(ppl);
3534   }
3535
3536   MapTimeBinLayers();
3537   delete [] zc;
3538   delete [] zmax;
3539   delete [] zmaxsensitive;
3540
3541 }
3542
3543 //______________________________________________________
3544
3545 Int_t  AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t localTB) const
3546 {
3547   //
3548   // depending on the digitization parameters calculates "global"
3549   // time bin index for timebin <localTB> in plane <plane>
3550   //
3551
3552   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
3553   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
3554
3555   Double_t dx = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
3556                          / AliTRDcalibDB::Instance()->GetSamplingFrequency();
3557
3558   Int_t tbAmp = fPar->GetTimeBefore();
3559   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
3560   if(kTRUE) maxAmp = 0;   // intentional until we change parameter class
3561   Int_t tbDrift = fPar->GetTimeMax();
3562   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx)+4; // MI change  - take also last time bins
3563
3564   Int_t tbPerPlane = TMath::Min(tbAmp,maxAmp) + TMath::Min(tbDrift,maxDrift);
3565
3566   Int_t gtb = (plane+1) * tbPerPlane - localTB - 1 - TMath::Min(tbAmp,maxAmp);
3567
3568   if((localTB < 0) &&
3569      (TMath::Abs(localTB) > TMath::Min(tbAmp,maxAmp))) return -1;
3570   if(localTB >= TMath::Min(tbDrift,maxDrift)) return -1;
3571
3572   return gtb;
3573
3574
3575 }
3576
3577 //______________________________________________________
3578
3579 void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers()
3580 {
3581   //
3582   // For all sensitive time bins sets corresponding layer index
3583   // in the array fTimeBins
3584   //
3585
3586   Int_t index;
3587
3588   for(Int_t i = 0; i < fN; i++) {
3589     index = fLayers[i]->GetTimeBinIndex();
3590
3591     //    printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX());
3592
3593     if(index < 0) continue;
3594     if(index >= (Int_t) kMaxTimeBinIndex) {
3595       printf("*** AliTRDtracker::MapTimeBinLayers: \n");
3596       printf("    index %d exceeds allowed maximum of %d!\n",
3597              index, kMaxTimeBinIndex-1);
3598       continue;
3599     }
3600     fTimeBinIndex[index] = i;
3601   }
3602
3603   Double_t x1, dx1, x2, dx2, gap;
3604
3605   for(Int_t i = 0; i < fN-1; i++) {
3606     x1 = fLayers[i]->GetX();
3607     dx1 = fLayers[i]->GetdX();
3608     x2 = fLayers[i+1]->GetX();
3609     dx2 = fLayers[i+1]->GetdX();
3610     gap = (x2 - dx2/2) - (x1 + dx1/2);
3611 //     if(gap < -0.01) {
3612 //       printf("*** warning: layers %d and %d are overlayed:\n",i,i+1);
3613 //       printf("             %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2);
3614 //     }
3615 //     if(gap > 0.01) {
3616 //       printf("*** warning: layers %d and %d have a large gap:\n",i,i+1);
3617 //       printf("             (%f - %f) - (%f + %f) = %f\n",
3618 //              x2, dx2/2, x1, dx1, gap);
3619 //     }
3620   }
3621 }
3622
3623
3624 //______________________________________________________
3625
3626
3627 Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const
3628 {
3629   //
3630   // Returns the number of time bin which in radial position is closest to <x>
3631   //
3632
3633   if(x >= fLayers[fN-1]->GetX()) return fN-1;
3634   if(x <= fLayers[0]->GetX()) return 0;
3635
3636   Int_t b=0, e=fN-1, m=(b+e)/2;
3637   for (; b<e; m=(b+e)/2) {
3638     if (x > fLayers[m]->GetX()) b=m+1;
3639     else e=m;
3640   }
3641   if(TMath::Abs(x - fLayers[m]->GetX()) >
3642      TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
3643   else return m;
3644
3645 }
3646
3647 //______________________________________________________
3648
3649 Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const
3650 {
3651   //
3652   // Returns number of the innermost SENSITIVE propagation layer
3653   //
3654
3655   return GetLayerNumber(0);
3656 }
3657
3658 //______________________________________________________
3659
3660 Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const
3661 {
3662   //
3663   // Returns number of the outermost SENSITIVE time bin
3664   //
3665
3666   return GetLayerNumber(GetNumberOfTimeBins() - 1);
3667 }
3668
3669 //______________________________________________________
3670
3671 Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const
3672 {
3673   //
3674   // Returns number of SENSITIVE time bins
3675   //
3676
3677   Int_t tb, layer;
3678   for(tb = kMaxTimeBinIndex-1; tb >=0; tb--) {
3679     layer = GetLayerNumber(tb);
3680     if(layer>=0) break;
3681   }
3682   return tb+1;
3683 }
3684
3685 //______________________________________________________
3686
3687 void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
3688 {
3689   //
3690   // Insert layer <pl> in fLayers array.
3691   // Layers are sorted according to X coordinate.
3692
3693   if ( fN == ((Int_t) kMaxLayersPerSector)) {
3694     printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
3695     return;
3696   }
3697   if (fN==0) {fLayers[fN++] = pl; return;}
3698   Int_t i=Find(pl->GetX());
3699
3700   memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
3701   fLayers[i]=pl; fN++;
3702
3703 }
3704
3705 //______________________________________________________
3706
3707 Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const
3708 {
3709   //
3710   // Returns index of the propagation layer nearest to X
3711   //
3712
3713   if (x <= fLayers[0]->GetX()) return 0;
3714   if (x > fLayers[fN-1]->GetX()) return fN;
3715   Int_t b=0, e=fN-1, m=(b+e)/2;
3716   for (; b<e; m=(b+e)/2) {
3717     if (x > fLayers[m]->GetX()) b=m+1;
3718     else e=m;
3719   }
3720   return m;
3721 }
3722
3723
3724
3725
3726
3727 //______________________________________________________
3728 void AliTRDtracker::AliTRDpropagationLayer::SetZ(Double_t* center, Double_t *w, Double_t *wsensitive )
3729 {
3730   //
3731   // set centers and the width of sectors
3732   for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
3733     fZc[icham] = center[icham];
3734     fZmax[icham] = w[icham];
3735     fZmaxSensitive[icham] = wsensitive[icham];
3736     //   printf("chamber\t%d\tzc\t%f\tzmax\t%f\tzsens\t%f\n",icham,fZc[icham],fZmax[icham],fZmaxSensitive[icham]);
3737   }
3738 }
3739 //______________________________________________________
3740
3741 void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
3742 {
3743   //
3744   // set centers and the width of sectors
3745   fHole = kFALSE;
3746   for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
3747     fIsHole[icham] = holes[icham];
3748     if (holes[icham]) fHole = kTRUE;
3749   }
3750 }
3751
3752
3753
3754 Bool_t AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
3755         Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &radLength,
3756         Bool_t &lookForCluster) const
3757 {
3758   //
3759   // Returns radial step <dx>, density <rho>, rad. length <radLength>,
3760   // and sensitivity <lookForCluster> in point <y,z>
3761   //
3762
3763   Double_t alpha =  AliTRDgeometry::GetAlpha();
3764   Double_t ymax  =  fX*TMath::Tan(0.5*alpha);
3765
3766
3767   dx  = fdX;
3768   rho = fRho;
3769   radLength  = fX0;
3770   lookForCluster = kFALSE;
3771   Bool_t cross =kFALSE;
3772   //
3773   //
3774   if ( (ymax-TMath::Abs(y))<3.){   //cross material
3775     rho*=40.;
3776     radLength*=40.;
3777     cross=kTRUE;
3778   }
3779   //
3780   // check dead regions in sensitive volume
3781     //
3782   Int_t zone=-1;
3783   for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
3784     if (TMath::Abs(z - fZc[ch]) > fZmax[ch]) continue;  //not in given zone
3785     //
3786     if  (TMath::Abs(z - fZc[ch]) < fZmaxSensitive[ch]){
3787       if (fTimeBinIndex>=0) lookForCluster = !(fIsHole[zone]);
3788       if(TMath::Abs(y) > fYmaxSensitive){
3789         lookForCluster = kFALSE;
3790       }
3791       if (fIsHole[zone]) {
3792         //if hole
3793         rho = 1.29e-3;
3794         radLength = 36.66;
3795       }
3796     }else{
3797       cross = kTRUE; rho = 2.7; radLength = 24.01;  //aluminium in between
3798     }
3799   }
3800   //
3801   if (fTimeBinIndex>=0) return cross;
3802   //
3803   //
3804   // check hole
3805   if (fHole==kFALSE) return cross;
3806   //
3807   for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
3808     if  (TMath::Abs(z - fZc[ch]) < fZmax[ch]){
3809       if (fIsHole[ch]) {
3810         //if hole
3811         rho = 1.29e-3;
3812         radLength = 36.66;
3813       }
3814     }
3815   }
3816   return cross;
3817 }
3818
3819 Int_t  AliTRDtracker::AliTRDpropagationLayer::GetZone( Double_t z) const
3820 {
3821   //
3822   //
3823   if (fTimeBinIndex < 0) return -20;  //unknown
3824   Int_t zone=-10;   // dead zone
3825   for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
3826     if(TMath::Abs(z - fZc[ch]) < fZmax[ch])
3827       zone = ch;
3828   }
3829   return zone;
3830 }
3831
3832
3833 //______________________________________________________
3834
3835 void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c,
3836                                                           UInt_t index) {
3837
3838 // Insert cluster in cluster array.
3839 // Clusters are sorted according to Y coordinate.
3840
3841   if(fTimeBinIndex < 0) {
3842     printf("*** attempt to insert cluster into non-sensitive time bin!\n");
3843     return;
3844   }
3845
3846   if (fN== (Int_t) kMaxClusterPerTimeBin) {
3847     printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n");
3848     return;
3849   }
3850   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
3851   Int_t i=Find(c->GetY());
3852   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
3853   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
3854   fIndex[i]=index; fClusters[i]=c; fN++;
3855 }
3856
3857 //______________________________________________________
3858
3859 Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const {
3860
3861 // Returns index of the cluster nearest in Y
3862
3863   if (fN<=0) return 0;
3864   if (y <= fClusters[0]->GetY()) return 0;
3865   if (y > fClusters[fN-1]->GetY()) return fN;
3866   Int_t b=0, e=fN-1, m=(b+e)/2;
3867   for (; b<e; m=(b+e)/2) {
3868     if (y > fClusters[m]->GetY()) b=m+1;
3869     else e=m;
3870   }
3871   return m;
3872 }
3873
3874 Int_t AliTRDtracker::AliTRDpropagationLayer::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad, Float_t maxroadz) const
3875 {
3876   //
3877   // Returns index of the cluster nearest to the given y,z
3878   //
3879   Int_t index = -1;
3880   Int_t maxn = fN;
3881   Float_t mindist = maxroad;
3882   //
3883   for (Int_t i=Find(y-maxroad); i<maxn; i++) {
3884     AliTRDcluster* c=(AliTRDcluster*)(fClusters[i]);
3885     Float_t ycl = c->GetY();
3886     //
3887     if (ycl > y+maxroad) break;
3888     if (TMath::Abs(c->GetZ()-z) > maxroadz) continue;
3889     if (TMath::Abs(ycl-y)<mindist){
3890       mindist = TMath::Abs(ycl-y);
3891       index = fIndex[i];
3892     }
3893   }
3894   return index;
3895 }
3896
3897
3898 //---------------------------------------------------------
3899
3900 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
3901 //
3902 //  Returns correction factor for tilted pads geometry
3903 //
3904   Int_t det = c->GetDetector();
3905   Int_t plane = fGeom->GetPlane(det);
3906   AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(plane,0);
3907   Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
3908
3909   if(fNoTilt) h01 = 0;
3910   return h01;
3911 }
3912
3913
3914 void AliTRDtracker::CookdEdxTimBin(AliTRDtrack& TRDtrack)
3915 {
3916   // *** ADDED TO GET MORE INFORMATION FOR TRD PID  ---- PS
3917   // This is setting fdEdxPlane and fTimBinPlane
3918   // Sums up the charge in each plane for track TRDtrack and also get the
3919   // Time bin for Max. Cluster
3920   // Prashant Shukla (shukla@physi.uni-heidelberg.de)
3921
3922   //  const Int_t kNPlane = AliTRDgeometry::Nplan();
3923   //  const Int_t kNPlane = 6;
3924   Double_t  clscharge[kNPlane], maxclscharge[kNPlane];
3925   Int_t  nCluster[kNPlane], timebin[kNPlane];
3926
3927   //Initialization of cluster charge per plane.
3928   for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
3929     clscharge[iPlane] = 0.0;
3930     nCluster[iPlane] = 0;
3931     timebin[iPlane] = -1;
3932     maxclscharge[iPlane] = 0.0;
3933   }
3934
3935   // Loop through all clusters associated to track TRDtrack
3936   Int_t nClus = TRDtrack.GetNumberOfClusters();  // from Kalmantrack
3937   for (Int_t iClus = 0; iClus < nClus; iClus++) {
3938     Double_t charge = TRDtrack.GetClusterdQdl(iClus);
3939     Int_t index = TRDtrack.GetClusterIndex(iClus);
3940     AliTRDcluster *TRDcluster = (AliTRDcluster *) GetCluster(index);
3941     if (!TRDcluster) continue;
3942     Int_t tb = TRDcluster->GetLocalTimeBin();
3943     if (!tb) continue;
3944     Int_t detector = TRDcluster->GetDetector();
3945     Int_t iPlane   = fGeom->GetPlane(detector);
3946     clscharge[iPlane] = clscharge[iPlane]+charge;
3947     if(charge > maxclscharge[iPlane]) {
3948       maxclscharge[iPlane] = charge;
3949       timebin[iPlane] = tb;
3950     }
3951     nCluster[iPlane]++;
3952   } // end of loop over cluster
3953
3954   // Setting the fdEdxPlane and fTimBinPlane variabales
3955   Double_t Total_ch = 0;
3956   for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
3957     // Quality control of TRD track.
3958     if (nCluster[iPlane]<= 5) {
3959       clscharge[iPlane]=0.0;
3960       timebin[iPlane]=-1;
3961     }
3962     if (nCluster[iPlane]) clscharge[iPlane] /= nCluster[iPlane];
3963     TRDtrack.SetPIDsignals(clscharge[iPlane], iPlane);
3964     TRDtrack.SetPIDTimBin(timebin[iPlane], iPlane);
3965     Total_ch= Total_ch+clscharge[iPlane];
3966   }
3967   //  Int_t i;
3968   //  Int_t nc=TRDtrack.GetNumberOfClusters();
3969   //  Float_t dedx=0;
3970   //  for (i=0; i<nc; i++) dedx += TRDtrack.GetClusterdQdl(i);
3971   //  dedx /= nc;
3972   //  for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
3973   //    TRDtrack.SetPIDsignals(dedx, iPlane);
3974   //    TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
3975   //  }
3976
3977 } // end of function
3978
3979
3980 Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters,AliTRDtracklet&tracklet)
3981 {
3982   //
3983   //
3984   //  try to find nearest clusters to the track in timebins from t0 to t1
3985   //
3986   //
3987   //
3988   // correction coeficients   - depends on TRD parameters  - to be changed according it
3989   //
3990
3991   Double_t x[100],yt[100],zt[100];
3992   Double_t xmean=0;   //reference x
3993   Double_t dz[10][100],dy[10][100];
3994   Float_t zmean[100], nmean[100];
3995   Int_t    clfound=0;
3996   Int_t    indexes[10][100];    // indexes of the clusters in the road
3997   AliTRDcluster *cl[10][100];   // pointers to the clusters in the road
3998   Int_t    best[10][100];       // index of best matching cluster
3999   //
4000   //
4001
4002   for (Int_t it=0;it<=t1-t0; it++){
4003     x[it]=0;
4004     yt[it]=0;
4005     zt[it]=0;
4006     clusters[it+t0]=-2;
4007     zmean[it]=0;
4008     nmean[it]=0;
4009     //
4010     for (Int_t ih=0;ih<10;ih++){
4011       indexes[ih][it]=-2;              //reset indexes1
4012       cl[ih][it]=0;
4013       dz[ih][it]=-100;
4014       dy[ih][it]=-100;
4015       best[ih][it]=0;
4016     }
4017   }
4018   //
4019   Double_t x0 = track->GetX();
4020   Double_t sigmaz = TMath::Sqrt(TMath::Abs(track->GetSigmaZ2()));
4021   Int_t nall=0;
4022   Int_t nfound=0;
4023   Double_t h01 =0;
4024   Int_t plane =-1;
4025   Float_t padlength=0;
4026   AliTRDtrack track2(*track);
4027   Float_t snpy = track->GetSnp();
4028   Float_t tany = TMath::Sqrt(snpy*snpy/(1.-snpy*snpy));
4029   if (snpy<0) tany*=-1;
4030   //
4031   Double_t sy2=ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
4032   Double_t sz2=ExpectedSigmaZ2(x0,track->GetTgl());
4033   Double_t road = 15.*sqrt(track->GetSigmaY2() + sy2);
4034   if (road>6.) road=6.;
4035
4036   //
4037   for (Int_t it=0;it<t1-t0;it++){
4038     Double_t maxChi2[2]={fgkMaxChi2,fgkMaxChi2};
4039     AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(it+t0));
4040     if (timeBin==0) continue;  // no indexes1
4041     Int_t maxn = timeBin;
4042     x[it] = timeBin.GetX();
4043     track2.PropagateTo(x[it]);
4044     yt[it] = track2.GetY();
4045     zt[it] = track2.GetZ();
4046
4047     Double_t  y=yt[it],z=zt[it];
4048     Double_t chi2 =1000000;
4049     nall++;
4050     //
4051     // find 2 nearest cluster at given time bin
4052     //
4053     //
4054     for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
4055       AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
4056       h01 = GetTiltFactor(c);
4057       if (plane<0){
4058         Int_t det = c->GetDetector();
4059         plane = fGeom->GetPlane(det);
4060         padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
4061       }
4062       //      if (c->GetLocalTimeBin()==0) continue;
4063       if (c->GetY() > y+road) break;
4064       if((c->GetZ()-z)*(c->GetZ()-z) > 12. * sz2) continue;
4065
4066       Double_t dist = TMath::Abs(c->GetZ()-z);
4067       if (dist> (0.5*padlength+6.*sigmaz)) continue;   // 6 sigma boundary cut
4068       Double_t cost = 0;
4069       //
4070       if (dist> (0.5*padlength-sigmaz)){   //  sigma boundary cost function
4071         cost =  (dist-0.5*padlength)/(2.*sigmaz);
4072         if (cost>-1) cost= (cost+1.)*(cost+1.);
4073         else cost=0;
4074       }
4075       //      Int_t label = TMath::Abs(track->GetLabel());
4076       //      if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
4077       chi2=track2.GetPredictedChi2(c,h01)+cost;
4078       //
4079       clfound++;
4080       if (chi2 > maxChi2[1]) continue;
4081
4082       for (Int_t ih=2;ih<9; ih++){  //store the clusters in the road
4083         if (cl[ih][it]==0){
4084           cl[ih][it] = c;
4085           indexes[ih][it] =timeBin.GetIndex(i);   // index - 9 - reserved for outliers
4086           break;
4087         }
4088       }
4089       //
4090       if (chi2 <maxChi2[0]){
4091         maxChi2[1]     = maxChi2[0];
4092         maxChi2[0]     = chi2;
4093         indexes[1][it] = indexes[0][it];
4094         cl[1][it]      = cl[0][it];
4095         indexes[0][it] = timeBin.GetIndex(i);
4096         cl[0][it]      = c;
4097         continue;
4098       }
4099       maxChi2[1]=chi2;
4100       cl[1][it] = c;
4101       indexes[1][it] =timeBin.GetIndex(i);
4102     }
4103     if (cl[0][it]){
4104       nfound++;
4105       xmean += x[it];
4106     }
4107   }
4108   //
4109   if (nfound<4) return 0;
4110   xmean /=Float_t(nfound);     // middle x
4111   track2.PropagateTo(xmean);   // propagate track to the center
4112   //
4113   // choose one of the variants
4114   //
4115   Int_t changes[10];
4116   Float_t sumz      = 0;
4117   Float_t sum       = 0;
4118   Double_t sumdy    = 0;
4119   Double_t sumdy2   = 0;
4120   Double_t sumx     = 0;
4121   Double_t sumxy    = 0;
4122   Double_t sumx2    = 0;
4123   Double_t mpads    = 0;
4124   //
4125   Int_t   ngood[10];
4126   Int_t   nbad[10];
4127   //
4128   Double_t meanz[10];
4129   Double_t moffset[10];    // mean offset
4130   Double_t mean[10];       // mean value
4131   Double_t angle[10];      // angle
4132   //
4133   Double_t smoffset[10];   // sigma of mean offset
4134   Double_t smean[10];      // sigma of mean value
4135   Double_t sangle[10];     // sigma of angle
4136   Double_t smeanangle[10]; // correlation
4137   //
4138   Double_t sigmas[10];
4139   Double_t tchi2s[10];      // chi2s for tracklet
4140   //
4141   // calculate zmean
4142   //
4143   for (Int_t it=0;it<t1-t0;it++){
4144     if (!cl[0][it]) continue;
4145     for (Int_t dt=-3;dt<=3;dt++){
4146       if (it+dt<0) continue;
4147       if (it+dt>t1-t0) continue;
4148       if (!cl[0][it+dt]) continue;
4149       zmean[it]+=cl[0][it+dt]->GetZ();
4150       nmean[it]+=1.;
4151     }
4152     zmean[it]/=nmean[it];
4153   }
4154   //
4155   for (Int_t it=0; it<t1-t0;it++){
4156     best[0][it]=0;
4157     for (Int_t ih=0;ih<10;ih++){
4158       dz[ih][it]=-100;
4159       dy[ih][it]=-100;
4160       if (!cl[ih][it]) continue;
4161       //Float_t poscor =  fgkCoef*(cl[ih][it]->GetLocalTimeBin() - fgkMean)+fgkOffset;
4162       Float_t poscor =  0;   // applied during loading of clusters
4163       if (cl[ih][it]->IsUsed()) poscor=0;  // correction already applied
4164       dz[ih][it]  = cl[ih][it]->GetZ()- zt[it];                               // calculate distance from track  in z
4165       dy[ih][it]  = cl[ih][it]->GetY()+ dz[ih][it]*h01 - poscor  -yt[it];     //                                in y
4166     }
4167     // minimize changes
4168     if (!cl[0][it]) continue;
4169     if (TMath::Abs(cl[0][it]->GetZ()-zmean[it])> padlength*0.8 &&cl[1][it])
4170       if (TMath::Abs(cl[1][it]->GetZ()-zmean[it])< padlength*0.5){
4171         best[0][it]=1;
4172       }
4173   }
4174   //
4175   // iterative choosing of "best path"
4176   //
4177   //
4178   Int_t label = TMath::Abs(track->GetLabel());
4179   Int_t bestiter=0;
4180   //
4181   for (Int_t iter=0;iter<9;iter++){
4182     //
4183     changes[iter]= 0;
4184     sumz      = 0; sum=0; sumdy=0;sumdy2=0;sumx=0;sumx2=0;sumxy=0;mpads=0; ngood[iter]=0; nbad[iter]=0;
4185     // linear fit
4186     for (Int_t it=0;it<t1-t0;it++){
4187       if (!cl[best[iter][it]][it]) continue;
4188       //calculates pad-row changes
4189       Double_t zbefore= cl[best[iter][it]][it]->GetZ();
4190       Double_t zafter = cl[best[iter][it]][it]->GetZ();
4191       for (Int_t itd = it-1; itd>=0;itd--) {
4192         if (cl[best[iter][itd]][itd]) {
4193           zbefore= cl[best[iter][itd]][itd]->GetZ();
4194           break;
4195         }
4196       }
4197       for (Int_t itd = it+1; itd<t1-t0;itd++) {
4198         if (cl[best[iter][itd]][itd]) {
4199           zafter= cl[best[iter][itd]][itd]->GetZ();
4200           break;
4201         }
4202       }
4203       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]++;
4204       //
4205       Double_t dx = x[it]-xmean;  // distance to reference x
4206       sumz += cl[best[iter][it]][it]->GetZ();
4207       sum++;
4208       sumdy += dy[best[iter][it]][it];
4209       sumdy2+= dy[best[iter][it]][it]*dy[best[iter][it]][it];
4210       sumx  += dx;
4211       sumx2 += dx*dx;
4212       sumxy  += dx*dy[best[iter][it]][it];
4213       mpads += cl[best[iter][it]][it]->GetNPads();
4214       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){
4215         ngood[iter]++;
4216       }
4217       else{
4218         nbad[iter]++;
4219       }
4220     }
4221     //
4222     // calculates line parameters
4223     //
4224     Double_t det  = sum*sumx2-sumx*sumx;
4225     angle[iter]   = (sum*sumxy-sumx*sumdy)/det;
4226     mean[iter]    = (sumx2*sumdy-sumx*sumxy)/det;
4227     meanz[iter]   = sumz/sum;
4228     moffset[iter] = sumdy/sum;
4229     mpads        /= sum;                         // mean number of pads
4230     //
4231     //
4232     Double_t  sigma2 = 0;   // normalized residuals - for line fit
4233     Double_t  sigma1 = 0;   // normalized residuals - constant fit
4234     //
4235     for (Int_t it=0;it<t1-t0;it++){
4236       if (!cl[best[iter][it]][it]) continue;
4237       Double_t dx = x[it]-xmean;
4238       Double_t ytr = mean[iter]+angle[iter]*dx;
4239       sigma2 += (dy[best[iter][it]][it]-ytr)*(dy[best[iter][it]][it]-ytr);
4240       sigma1 +=  (dy[best[iter][it]][it]-moffset[iter])*(dy[best[iter][it]][it]-moffset[iter]);
4241       sum++;
4242     }
4243     sigma2      /=(sum-2);                    // normalized residuals
4244     sigma1      /=(sum-1);                    // normalized residuals
4245     //
4246     smean[iter]       = sigma2*(sumx2/det);   // estimated error2 of mean
4247     sangle[iter]      = sigma2*(sum/det);     // estimated error2 of angle
4248     smeanangle[iter]  = sigma2*(-sumx/det);   // correlation
4249     //
4250     //
4251     sigmas[iter]  = TMath::Sqrt(sigma1);      //
4252     smoffset[iter]= (sigma1/sum)+0.01*0.01;             // sigma of mean offset + unisochronity sigma
4253     //
4254     // iterative choosing of "better path"
4255     //
4256     for (Int_t it=0;it<t1-t0;it++){
4257       if (!cl[best[iter][it]][it]) continue;
4258       //
4259       Double_t sigmatr2 = smoffset[iter]+0.5*tany*tany;             //add unisochronity + angular effect contribution
4260       Double_t sweight  = 1./sigmatr2+1./track->GetSigmaY2();
4261       Double_t weighty  = (moffset[iter]/sigmatr2)/sweight;         // weighted mean
4262       Double_t sigmacl  = TMath::Sqrt(sigma1*sigma1+track->GetSigmaY2());   //
4263       Double_t mindist=100000;
4264       Int_t ihbest=0;
4265       for (Int_t ih=0;ih<10;ih++){
4266         if (!cl[ih][it]) break;
4267         Double_t dist2 = (dy[ih][it]-weighty)/sigmacl;
4268         dist2*=dist2;    //chi2 distance
4269         if (dist2<mindist){
4270           mindist = dist2;
4271           ihbest =ih;
4272         }
4273       }
4274       best[iter+1][it]=ihbest;
4275     }
4276     //
4277     //  update best hypothesy if better chi2 according tracklet position and angle
4278     //
4279     Double_t sy2 = smean[iter]  + track->GetSigmaY2();
4280     Double_t sa2 = sangle[iter] + track->fCee;
4281     Double_t say = track->fCey;
4282     //    Double_t chi20 = mean[bestiter]*mean[bestiter]/sy2+angle[bestiter]*angle[bestiter]/sa2;
4283     // Double_t chi21 = mean[iter]*mean[iter]/sy2+angle[iter]*angle[iter]/sa2;
4284
4285     Double_t detchi    = sy2*sa2-say*say;
4286     Double_t invers[3] = {sa2/detchi, sy2/detchi, -say/detchi};   //inverse value of covariance matrix
4287
4288     Double_t chi20 = mean[bestiter]*mean[bestiter]*invers[0]+angle[bestiter]*angle[bestiter]*invers[1]+
4289       2.*mean[bestiter]*angle[bestiter]*invers[2];
4290     Double_t chi21 = mean[iter]*mean[iter]*invers[0]+angle[iter]*angle[iter]*invers[1]+
4291       2*mean[iter]*angle[iter]*invers[2];
4292     tchi2s[iter] =chi21;
4293     //
4294     if (changes[iter]<=changes[bestiter] && chi21<chi20) {
4295       bestiter =iter;
4296     }
4297   }
4298   //
4299   //set clusters
4300   //
4301   Double_t sigma2 = sigmas[0];   // choose as sigma  from 0 iteration
4302   Short_t maxpos    = -1;
4303   Float_t maxcharge =  0;
4304   Short_t maxpos4    = -1;
4305   Float_t maxcharge4 =  0;
4306   Short_t maxpos5    = -1;
4307   Float_t maxcharge5 =  0;
4308
4309   //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
4310   //if (tchi2s[bestiter]>25.) sigma2=1000.;  // dont'accept
4311
4312   Double_t expectederr = sigma2*sigma2+0.01*0.01;
4313   if (mpads>3.5) expectederr  +=   (mpads-3.5)*0.04;
4314   if (changes[bestiter]>1) expectederr+=   changes[bestiter]*0.01;
4315   expectederr+=(0.03*(tany-fgkExB)*(tany-fgkExB))*15;
4316   //  if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
4317   //expectederr+=10000;
4318   for (Int_t it=0;it<t1-t0;it++){
4319     if (!cl[best[bestiter][it]][it]) continue;
4320     //    Float_t poscor =  fgkCoef*(cl[best[bestiter][it]][it]->GetLocalTimeBin() - fgkMean)+fgkOffset;
4321     Float_t poscor =  0;           //applied during loading of cluster
4322     cl[best[bestiter][it]][it]->SetSigmaY2(expectederr);  // set cluster error
4323     if (!cl[best[bestiter][it]][it]->IsUsed()){
4324       cl[best[bestiter][it]][it]->SetY( cl[best[bestiter][it]][it]->GetY()-poscor);  // ExB corrction correction
4325       //      cl[best[bestiter][it]][it]->Use();
4326     }
4327     //
4328     //  time bins with maximal charge
4329     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge){
4330       maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4331       maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
4332     }
4333
4334     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge4){
4335       if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=4){
4336         maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4337         maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
4338       }
4339     }
4340     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge5){
4341       if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=5){
4342         maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4343         maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
4344       }
4345     }
4346     //
4347     //  time bins with maximal charge
4348     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge){
4349       maxcharge = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4350       maxpos = cl[best[bestiter][it]][it]->GetLocalTimeBin();
4351     }
4352
4353     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge4){
4354       if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=4){
4355         maxcharge4 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4356         maxpos4 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
4357       }
4358     }
4359     if (TMath::Abs(cl[best[bestiter][it]][it]->GetQ())> maxcharge5){
4360       if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>=5){
4361         maxcharge5 = TMath::Abs(cl[best[bestiter][it]][it]->GetQ());
4362         maxpos5 = cl[best[bestiter][it]][it]->GetLocalTimeBin();
4363       }
4364     }
4365     clusters[it+t0] = indexes[best[bestiter][it]][it];
4366     //if (cl[best[bestiter][it]][it]->GetLocalTimeBin()>4 && cl[best[bestiter][it]][it]->GetLocalTimeBin()<18) clusters[it+t0] = indexes[best[bestiter][it]][it];    //Test
4367   }
4368   //
4369   // set tracklet parameters
4370   //
4371   Double_t trackleterr2 = smoffset[bestiter]+0.01*0.01;
4372   if (mpads>3.5) trackleterr2  +=   (mpads-3.5)*0.04;
4373   trackleterr2+=   changes[bestiter]*0.01;
4374   trackleterr2*=   TMath::Max(14.-nfound,1.);
4375   trackleterr2+=   0.2*(tany-fgkExB)*(tany-fgkExB);
4376   //
4377   tracklet.Set(xmean, track2.GetY()+moffset[bestiter], meanz[bestiter], track2.GetAlpha(), trackleterr2);  //set tracklet parameters
4378   tracklet.SetTilt(h01);
4379   tracklet.SetP0(mean[bestiter]);
4380   tracklet.SetP1(angle[bestiter]);
4381   tracklet.SetN(nfound);
4382   tracklet.SetNCross(changes[bestiter]);
4383   tracklet.SetPlane(plane);
4384   tracklet.SetSigma2(expectederr);
4385   tracklet.SetChi2(tchi2s[bestiter]);
4386   tracklet.SetMaxPos(maxpos,maxpos4,maxpos5);
4387   track->fTracklets[plane] = tracklet;
4388   track->fNWrong+=nbad[0];
4389   //
4390   // Debuging part
4391   //
4392   TClonesArray array0("AliTRDcluster");
4393   TClonesArray array1("AliTRDcluster");
4394   array0.ExpandCreateFast(t1-t0+1);
4395   array1.ExpandCreateFast(t1-t0+1);
4396   TTreeSRedirector& cstream = *fDebugStreamer;
4397   AliTRDcluster dummy;
4398   Double_t dy0[100];
4399   Double_t dyb[100];
4400
4401   for (Int_t it=0;it<t1-t0;it++){
4402     dy0[it] = dy[0][it];
4403     dyb[it] = dy[best[bestiter][it]][it];
4404     if(cl[0][it]) {
4405       new(array0[it]) AliTRDcluster(*cl[0][it]);
4406     }
4407     else{
4408       new(array0[it]) AliTRDcluster(dummy);
4409     }
4410     if(cl[best[bestiter][it]][it]) {
4411       new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
4412     }
4413     else{
4414       new(array1[it]) AliTRDcluster(dummy);
4415     }
4416   }
4417   TGraph graph0(t1-t0,x,dy0);
4418   TGraph graph1(t1-t0,x,dyb);
4419   TGraph graphy(t1-t0,x,yt);
4420   TGraph graphz(t1-t0,x,zt);
4421   //
4422   //
4423   cstream<<"tracklet"<<
4424     "track.="<<track<<                                       // track parameters
4425     "tany="<<tany<<                                          // tangent of the local track angle
4426     "xmean="<<xmean<<                                        // xmean - reference x of tracklet
4427     "tilt="<<h01<<                                           // tilt angle
4428     "nall="<<nall<<                                          // number of foundable clusters
4429     "nfound="<<nfound<<                                      // number of found clusters
4430     "clfound="<<clfound<<                                    // total number of found clusters in road
4431     "mpads="<<mpads<<                                        // mean number of pads per cluster
4432     "plane="<<plane<<                                        // plane number
4433     "road="<<road<<                                          // the width of the used road
4434     "graph0.="<<&graph0<<                                    // x - y = dy for closest cluster
4435     "graph1.="<<&graph1<<                                    // x - y = dy for second closest cluster
4436     "graphy.="<<&graphy<<                                    // y position of the track
4437     "graphz.="<<&graphz<<                                    // z position of the track
4438     //    "fCl.="<<&array0<<                                       // closest cluster
4439     // "fCl2.="<<&array1<<                                      // second closest cluster
4440     "maxpos="<<maxpos<<                                      // maximal charge postion
4441     "maxcharge="<<maxcharge<<                                // maximal charge
4442     "maxpos4="<<maxpos4<<                                    // maximal charge postion - after bin 4
4443     "maxcharge4="<<maxcharge4<<                              // maximal charge         - after bin 4
4444     "maxpos5="<<maxpos5<<                                    // maximal charge postion - after bin 5
4445     "maxcharge5="<<maxcharge5<<                              // maximal charge         - after bin 5
4446     //
4447     "bestiter="<<bestiter<<                                  // best iteration number
4448     "tracklet.="<<&tracklet<<                                // corrspond to the best iteration
4449     "tchi20="<<tchi2s[0]<<                                   // chi2 of cluster in the 0 iteration
4450     "tchi2b="<<tchi2s[bestiter]<<                            // chi2 of cluster in the best  iteration
4451     "sigmas0="<<sigmas[0]<<                                  // residuals sigma
4452     "sigmasb="<<sigmas[bestiter]<<                           // residulas sigma
4453     //
4454     "ngood0="<<ngood[0]<<                                    // number of good clusters in 0 iteration
4455     "nbad0="<<nbad[0]<<                                      // number of bad clusters in 0 iteration
4456     "ngoodb="<<ngood[bestiter]<<                             //                        in  best iteration
4457     "nbadb="<<nbad[bestiter]<<                               //                        in  best iteration
4458     //
4459     "changes0="<<changes[0]<<                                // changes of pardrows in iteration number 0
4460     "changesb="<<changes[bestiter]<<                         // changes of pardrows in best iteration
4461     //
4462     "moffset0="<<moffset[0]<<                                // offset fixing angle in iter=0
4463     "smoffset0="<<smoffset[0]<<                              // sigma of offset fixing angle in iter=0
4464     "moffsetb="<<moffset[bestiter]<<                         // offset fixing angle in iter=best
4465     "smoffsetb="<<smoffset[bestiter]<<                       // sigma of offset fixing angle in iter=best
4466     //
4467     "mean0="<<mean[0]<<                                      // mean dy in iter=0;
4468     "smean0="<<smean[0]<<                                    // sigma of mean dy in iter=0
4469     "meanb="<<mean[bestiter]<<                               // mean dy in iter=best
4470     "smeanb="<<smean[bestiter]<<                             // sigma of mean dy in iter=best
4471     //
4472     "angle0="<<angle[0]<<                                    // angle deviation in the iteration number 0
4473     "sangle0="<<sangle[0]<<                                  // sigma of angular deviation in iteration number 0
4474     "angleb="<<angle[bestiter]<<                             // angle deviation in the best iteration
4475     "sangleb="<<sangle[bestiter]<<                           // sigma of angle deviation in the best iteration
4476     //
4477     "expectederr="<<expectederr<<                            // expected error of cluster position
4478     "\n";
4479   //
4480   //
4481   return nfound;
4482 }
4483
4484
4485 Int_t  AliTRDtracker::Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down)
4486 {
4487   //
4488   //  Sort eleements according occurancy
4489   //  The size of output array has is 2*n
4490   //
4491   Int_t * sindexS = new Int_t[n];     // temp array for sorting
4492   Int_t * sindexF = new Int_t[2*n];
4493   for (Int_t i=0;i<n;i++) sindexF[i]=0;
4494   //
4495   TMath::Sort(n,inlist, sindexS, down);
4496   Int_t last      = inlist[sindexS[0]];
4497   Int_t val       = last;
4498   sindexF[0]      = 1;
4499   sindexF[0+n]    = last;
4500   Int_t countPos  = 0;
4501   //
4502   //  find frequency
4503   for(Int_t i=1;i<n; i++){
4504     val = inlist[sindexS[i]];
4505     if (last == val)   sindexF[countPos]++;
4506     else{
4507       countPos++;
4508       sindexF[countPos+n] = val;
4509       sindexF[countPos]++;
4510       last =val;
4511     }
4512   }
4513   if (last==val) countPos++;
4514   // sort according frequency
4515   TMath::Sort(countPos, sindexF, sindexS, kTRUE);
4516   for (Int_t i=0;i<countPos;i++){
4517     outlist[2*i  ] = sindexF[sindexS[i]+n];
4518     outlist[2*i+1] = sindexF[sindexS[i]];
4519   }
4520   delete [] sindexS;
4521   delete [] sindexF;
4522
4523   return countPos;
4524 }
4525
4526 AliTRDtrack * AliTRDtracker::RegisterSeed(AliTRDseed * seeds, Double_t * params)
4527 {
4528   //
4529   //
4530   //
4531   Double_t alpha=AliTRDgeometry::GetAlpha();
4532   Double_t shift=AliTRDgeometry::GetAlpha()/2.;
4533   Double_t c[15];
4534   c[0] = 0.2;
4535   c[1] = 0  ; c[2] = 2;
4536   c[3] = 0  ; c[4] = 0; c[5] = 0.02;
4537   c[6] = 0  ; c[7] = 0; c[8] = 0;      c[9] = 0.1;
4538   c[10] = 0  ; c[11] = 0; c[12] = 0;   c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
4539   //
4540   Int_t index =0;
4541   AliTRDcluster *cl =0;
4542   for (Int_t ilayer=0;ilayer<6;ilayer++){
4543     if (seeds[ilayer].isOK()){
4544       for (Int_t itime=22;itime>0;itime--){
4545         if (seeds[ilayer].fIndexes[itime]>0){
4546           index = seeds[ilayer].fIndexes[itime];
4547           cl = seeds[ilayer].fClusters[itime];
4548           break;
4549         }
4550       }
4551     }
4552     if (index>0) break;
4553   }
4554   if (cl==0) return 0;
4555   AliTRDtrack * track  = new AliTRDtrack(cl,index,&params[1],c, params[0],params[6]*alpha+shift);
4556   track->PropagateTo(params[0]-5.);
4557   track->ResetCovariance(1);
4558   //
4559   Int_t rc=FollowBackProlongationG(*track);
4560   if (rc<30) {
4561     delete track;
4562     track =0;
4563   }else{
4564     track->CookdEdx();
4565     CookdEdxTimBin(*track);
4566     CookLabel(track, 0.9);
4567   }
4568   return track;
4569 }
4570
4571
4572
4573
4574
4575
4576 AliTRDseed::AliTRDseed()
4577 {
4578   //
4579   //
4580   fTilt =0;         // tilting angle
4581   fPadLength = 0;   // pad length
4582   fX0 = 0;           // x0 position
4583   for (Int_t i=0;i<25;i++){
4584     fX[i]=0;        // !x position
4585     fY[i]=0;        // !y position
4586     fZ[i]=0;        // !z position
4587     fIndexes[i]=0;  // !indexes
4588     fClusters[i]=0; // !clusters
4589   }
4590   for (Int_t i=0;i<2;i++){
4591     fYref[i]=0;      // reference y
4592     fZref[i]=0;      // reference z
4593     fYfit[i]=0;      // y fit position +derivation
4594     fYfitR[i]=0;      // y fit position +derivation
4595     fZfit[i]=0;      // z fit position
4596     fZfitR[i]=0;      // z fit position
4597     fLabels[i]=0;    // labels
4598   }
4599   fSigmaY  = 0;
4600   fSigmaY2 = 0;
4601   fMeanz=0;         // mean vaue of z
4602   fZProb=0;         // max probbable z
4603   fMPads=0;
4604   //
4605   fN=0;            // number of associated clusters
4606   fN2=0;            // number of not crossed
4607   fNUsed=0;        // number of used clusters
4608   fNChange=0;      // change z counter
4609 }
4610
4611 void AliTRDseed::Reset(){
4612   //
4613   // reset seed
4614   //
4615   for (Int_t i=0;i<25;i++){
4616     fX[i]=0;        // !x position
4617     fY[i]=0;        // !y position
4618     fZ[i]=0;        // !z position
4619     fIndexes[i]=0;  // !indexes
4620     fClusters[i]=0; // !clusters
4621     fUsable[i]  = kFALSE;
4622   }
4623   for (Int_t i=0;i<2;i++){
4624     fYref[i]=0;      // reference y
4625     fZref[i]=0;      // reference z
4626     fYfit[i]=0;      // y fit position +derivation
4627     fYfitR[i]=0;      // y fit position +derivation
4628     fZfit[i]=0;      // z fit position
4629     fZfitR[i]=0;      // z fit position
4630     fLabels[i]=-1;    // labels
4631   }
4632   fSigmaY =0;         //"robust" sigma in y
4633   fSigmaY2=0;         //"robust" sigma in y
4634   fMeanz =0;         // mean vaue of z
4635   fZProb =0;         // max probbable z
4636   fMPads =0;
4637   //
4638   fN=0;            // number of associated clusters
4639   fN2=0;            // number of not crossed
4640   fNUsed=0;        // number of used clusters
4641   fNChange=0;      // change z counter
4642 }
4643
4644 void AliTRDseed::CookLabels(){
4645   //
4646   // cook 2 labels for seed
4647   //
4648   Int_t labels[200];
4649   Int_t out[200];
4650   Int_t nlab =0;
4651   for (Int_t i=0;i<25;i++){
4652     if (!fClusters[i]) continue;
4653     for (Int_t ilab=0;ilab<3;ilab++){
4654       if (fClusters[i]->GetLabel(ilab)>=0){
4655         labels[nlab] = fClusters[i]->GetLabel(ilab);
4656         nlab++;
4657       }
4658     }
4659   }
4660   Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
4661   fLabels[0] = out[0];
4662   if (nlab2>1 && out[3]>1) fLabels[1] =out[2];
4663 }
4664
4665 void   AliTRDseed::UseClusters()
4666 {
4667   //
4668   // use clusters
4669   //
4670    for (Int_t i=0;i<25;i++){
4671      if (!fClusters[i]) continue;
4672      if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
4673    }
4674 }
4675
4676
4677 void        AliTRDseed::Update(){
4678   //
4679   //
4680   //
4681   const Float_t ratio = 0.8;
4682   const Int_t   kClmin        = 6;
4683   const Float_t kmaxtan  = 2;
4684   if (TMath::Abs(fYref[1])>kmaxtan) return;             // too much inclined track
4685   //
4686   Float_t  sigmaexp = 0.05+TMath::Abs(fYref[1]*0.25);   // expected r.m.s in y direction
4687   Float_t  ycrosscor = fPadLength*fTilt*0.5;             // y correction for crossing
4688   fNChange =0;
4689   //
4690   Double_t sumw, sumwx,sumwx2;
4691   Double_t sumwy, sumwxy, sumwz,sumwxz;
4692   Int_t    zints[25];        // histograming of the z coordinate - get 1 and second max probable coodinates in z
4693   Int_t    zouts[50];        //
4694   Float_t  allowedz[25];     // allowed z for given time bin
4695   Float_t  yres[25];         // residuals from reference
4696   Float_t  anglecor = fTilt*fZref[1];  //correction to the angle
4697   //
4698   //
4699   fN=0; fN2 =0;
4700   for (Int_t i=0;i<25;i++){
4701     yres[i] =10000;
4702     if (!fClusters[i]) continue;
4703     yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i];   // residual y
4704     zints[fN] = Int_t(fZ[i]);
4705     fN++;
4706   }
4707   if (fN<kClmin) return;
4708   Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
4709   fZProb   = zouts[0];
4710   if (nz<=1) zouts[3]=0;
4711   if (zouts[1]+zouts[3]<kClmin) return;
4712   //
4713   if (TMath::Abs(zouts[0]-zouts[2])>12.) zouts[3]=0;   // z distance bigger than pad - length
4714   //
4715   Int_t  breaktime = -1;
4716   Bool_t mbefore   = kFALSE;
4717   Int_t  cumul[25][2];
4718   Int_t  counts[2]={0,0};
4719   //
4720   if (zouts[3]>=3){
4721     //
4722     // find the break time allowing one chage on pad-rows with maximal numebr of accepted clusters
4723     //
4724     fNChange=1;
4725     for (Int_t i=0;i<25;i++){
4726       cumul[i][0] = counts[0];
4727       cumul[i][1] = counts[1];
4728       if (TMath::Abs(fZ[i]-zouts[0])<2) counts[0]++;
4729       if (TMath::Abs(fZ[i]-zouts[2])<2) counts[1]++;
4730     }
4731     Int_t  maxcount  = 0;
4732     for (Int_t i=0;i<24;i++) {
4733       Int_t after  = cumul[24][0]-cumul[i][0];
4734       Int_t before = cumul[i][1];
4735       if (after+before>maxcount) {
4736         maxcount=after+before;
4737         breaktime=i;
4738         mbefore=kFALSE;
4739       }
4740       after  = cumul[24][1]-cumul[i][1];
4741       before = cumul[i][0];
4742       if (after+before>maxcount) {
4743         maxcount=after+before;
4744         breaktime=i;
4745         mbefore=kTRUE;
4746       }
4747     }
4748     breaktime-=1;
4749   }
4750   for (Int_t i=0;i<25;i++){
4751     if (i>breaktime)  allowedz[i] =   mbefore  ? zouts[2]:zouts[0];
4752     if (i<=breaktime) allowedz[i] = (!mbefore) ? zouts[2]:zouts[0];
4753   }
4754   if ( (allowedz[0]>allowedz[24] && fZref[1]<0) || (allowedz[0]<allowedz[24] &&  fZref[1]>0)){
4755     //
4756     // tracklet z-direction not in correspondance with track z direction
4757     //
4758     fNChange =0;
4759     for (Int_t i=0;i<25;i++){
4760       allowedz[i] =  zouts[0];  //only longest taken
4761     }
4762   }
4763   //
4764   if (fNChange>0){
4765     //
4766     // cross pad -row tracklet  - take the step change into account
4767     //
4768     for (Int_t i=0;i<25;i++){
4769       if (!fClusters[i]) continue;
4770       if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
4771       yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i];   // residual y
4772       if (TMath::Abs(fZ[i]-fZProb)>2){
4773         if (fZ[i]>fZProb) yres[i]+=fTilt*fPadLength;
4774         if (fZ[i]<fZProb) yres[i]-=fTilt*fPadLength;
4775       }
4776     }
4777   }
4778   //
4779   Double_t yres2[25];
4780   Double_t mean,sigma;
4781   for (Int_t i=0;i<25;i++){
4782     if (!fClusters[i]) continue;
4783     if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
4784     yres2[fN2] =  yres[i];
4785     fN2++;
4786   }
4787   if (fN2<kClmin){
4788     fN2 = 0;
4789     return;
4790   }
4791   EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*ratio-2));
4792   if (sigma<sigmaexp*0.8) sigma=sigmaexp;
4793   fSigmaY = sigma;
4794   //
4795   //
4796   // reset sums
4797   sumw=0; sumwx=0; sumwx2=0;
4798   sumwy=0; sumwxy=0; sumwz=0;sumwxz=0;
4799   fN2 =0;
4800   fMeanz =0;
4801   fMPads =0;
4802   //
4803   for (Int_t i=0;i<25;i++){
4804     fUsable[i]=kFALSE;
4805     if (!fClusters[i]) continue;
4806     if (TMath::Abs(fZ[i]-allowedz[i])>2)  continue;
4807     if (TMath::Abs(yres[i]-mean)>4.*sigma) continue;
4808     fUsable[i] = kTRUE;
4809     fN2++;
4810     fMPads+=fClusters[i]->GetNPads();
4811     Float_t weight =1;
4812     if (fClusters[i]->GetNPads()>4) weight=0.5;
4813     if (fClusters[i]->GetNPads()>5) weight=0.2;
4814     //
4815     Double_t x = fX[i];
4816     sumw+=weight; sumwx+=x*weight; sumwx2+=x*x*weight;
4817     sumwy+=weight*yres[i];  sumwxy+=weight*(yres[i])*x;
4818     sumwz+=weight*fZ[i];    sumwxz+=weight*fZ[i]*x;
4819   }
4820   if (fN2<kClmin){
4821     fN2 = 0;
4822     return;
4823   }
4824   fMeanz       = sumwz/sumw;
4825   Float_t correction =0;
4826   if (fNChange>0){
4827     // tracklet on boundary
4828     if (fMeanz<fZProb) correction =   ycrosscor;
4829     if (fMeanz>fZProb) correction =  -ycrosscor;
4830   }
4831   Double_t det = sumw*sumwx2-sumwx*sumwx;
4832   fYfitR[0]    = (sumwx2*sumwy-sumwx*sumwxy)/det;
4833   fYfitR[1]    = (sumw*sumwxy-sumwx*sumwy)/det;
4834   //
4835   fSigmaY2     =0;
4836   for (Int_t i=0;i<25;i++){
4837     if (!fUsable[i]) continue;
4838     Float_t delta = yres[i]-fYfitR[0]-fYfitR[1]*fX[i];
4839     fSigmaY2+=delta*delta;
4840   }
4841   fSigmaY2 = TMath::Sqrt(fSigmaY2/Float_t(fN2-2));
4842   //
4843   fZfitR[0]    = (sumwx2*sumwz-sumwx*sumwxz)/det;
4844   fZfitR[1]    = (sumw*sumwxz-sumwx*sumwz)/det;
4845   fZfit[0]     = (sumwx2*sumwz-sumwx*sumwxz)/det;
4846   fZfit[1]     = (sumw*sumwxz-sumwx*sumwz)/det;
4847   fYfitR[0]   += fYref[0]+correction;
4848   fYfitR[1]   += fYref[1];
4849   fYfit[0]     = fYfitR[0];
4850   fYfit[1]     = fYfitR[1];
4851   //
4852   //
4853   UpdateUsed();
4854 }
4855
4856
4857
4858
4859
4860
4861 void AliTRDseed::UpdateUsed(){
4862   //
4863   fNUsed =0;
4864   for (Int_t i=0;i<25;i++){
4865      if (!fClusters[i]) continue;
4866      if ((fClusters[i]->IsUsed())) fNUsed++;
4867   }
4868 }
4869
4870
4871 void AliTRDseed::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh)
4872 {
4873   //
4874   // robust estimator in 1D case MI version
4875   //
4876   //for the univariate case
4877   //estimates of location and scatter are returned in mean and sigma parameters
4878   //the algorithm works on the same principle as in multivariate case -
4879   //it finds a subset of size hh with smallest sigma, and then returns mean and
4880   //sigma of this subset
4881
4882   if (hh==0)
4883     hh=(nvectors+2)/2;
4884   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};
4885   Int_t *index=new Int_t[nvectors];
4886   TMath::Sort(nvectors, data, index, kFALSE);
4887   //
4888   Int_t    nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
4889   Double_t factor = faclts[nquant-1];
4890   //
4891   //
4892   Double_t sumx  =0;
4893   Double_t sumx2 =0;
4894   Int_t    bestindex = -1;
4895   Double_t bestmean  = 0;
4896   Double_t bestsigma = data[index[nvectors-1]]-data[index[0]];   // maximal possible sigma
4897   for (Int_t i=0; i<hh; i++){
4898     sumx  += data[index[i]];
4899     sumx2 += data[index[i]]*data[index[i]];
4900   }
4901   //
4902   Double_t norm = 1./Double_t(hh);
4903   Double_t norm2 = 1./Double_t(hh-1);
4904   for (Int_t i=hh; i<nvectors; i++){
4905     Double_t cmean  = sumx*norm;
4906     Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
4907     if (csigma<bestsigma){
4908       bestmean  = cmean;
4909       bestsigma = csigma;
4910       bestindex = i-hh;
4911     }
4912     //
4913     //
4914     sumx  += data[index[i]]-data[index[i-hh]];
4915     sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
4916   }
4917
4918   Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
4919   mean  = bestmean;
4920   sigma = bstd;
4921   delete [] index;
4922 }
4923
4924
4925 Float_t   AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror){
4926   //
4927   //
4928   //
4929   TLinearFitter fitterT2(4,"hyp4");  // fitting with tilting pads - kz not fixed
4930   fitterT2.StoreData(kTRUE);
4931   Float_t xref2 = (cseed[2].fX0+cseed[3].fX0)*0.5; // reference x0 for z
4932   //
4933   Int_t npointsT =0;
4934   fitterT2.ClearPoints();
4935   for (Int_t iLayer=0; iLayer<6;iLayer++){
4936     if (!cseed[iLayer].isOK()) continue;
4937     Double_t tilt = cseed[iLayer].fTilt;
4938
4939     for (Int_t itime=0;itime<25;itime++){
4940       if (!cseed[iLayer].fUsable[itime]) continue;
4941       Double_t x   = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2;  // x relative to the midle chamber
4942       Double_t y   = cseed[iLayer].fY[itime];
4943       Double_t z   = cseed[iLayer].fZ[itime];
4944       // tilted rieman
4945       //
4946       Double_t uvt[6];
4947       Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0;      // global x
4948       Double_t t = 1./(x2*x2+y*y);
4949       uvt[1]  = t;    // t
4950       uvt[0]  = 2.*x2*uvt[1];      // u
4951       uvt[2]  = 2.0*tilt*uvt[1];
4952       uvt[3]  = 2.0*tilt*x*uvt[1];
4953       uvt[4]  = 2.0*(y+tilt*z)*uvt[1];
4954       //
4955       Double_t error = 2*uvt[1];
4956       if (terror) error*=cseed[iLayer].fSigmaY;
4957       else {error *=0.2;} //default error
4958       fitterT2.AddPoint(uvt,uvt[4],error);
4959       npointsT++;
4960     }
4961   }
4962   fitterT2.Eval();
4963   Double_t rpolz0 = fitterT2.GetParameter(3);
4964   Double_t rpolz1 = fitterT2.GetParameter(4);
4965   //
4966   // linear fitter  - not possible to make boundaries
4967   // non accept non possible z and dzdx combination
4968   //
4969   Bool_t   acceptablez =kTRUE;
4970   for (Int_t iLayer=0; iLayer<6;iLayer++){
4971     if (cseed[iLayer].isOK()){
4972       Double_t zT2 =  rpolz0+rpolz1*(cseed[iLayer].fX0 - xref2);
4973       if (TMath::Abs(cseed[iLayer].fZProb-zT2)>cseed[iLayer].fPadLength*0.5+1)
4974         acceptablez = kFALSE;
4975     }
4976   }
4977   if (!acceptablez){
4978     Double_t zmf  = cseed[2].fZref[0]+cseed[2].fZref[1]*(xref2-cseed[2].fX0);
4979     Double_t dzmf = (cseed[2].fZref[1]+ cseed[3].fZref[1])*0.5;
4980     fitterT2.FixParameter(3,zmf);
4981     fitterT2.FixParameter(4,dzmf);
4982     fitterT2.Eval();
4983     fitterT2.ReleaseParameter(3);
4984     fitterT2.ReleaseParameter(4);
4985     rpolz0 = fitterT2.GetParameter(3);
4986     rpolz1 = fitterT2.GetParameter(4);
4987   }
4988   //
4989   Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
4990   Double_t params[3];
4991   params[0]     =  fitterT2.GetParameter(0);
4992   params[1]     =  fitterT2.GetParameter(1);
4993   params[2]     =  fitterT2.GetParameter(2);
4994   Double_t CR     =  1+params[1]*params[1]-params[2]*params[0];
4995   for (Int_t iLayer = 0; iLayer<6;iLayer++){
4996     Double_t  x = cseed[iLayer].fX0;
4997     Double_t  y=0,dy=0, z=0, dz=0;
4998     // y
4999     Double_t res2 = (x*params[0]+params[1]);
5000     res2*=res2;
5001     res2 = 1.-params[2]*params[0]+params[1]*params[1]-res2;
5002     if (res2>=0){
5003       res2 = TMath::Sqrt(res2);
5004       y    = (1-res2)/params[0];
5005     }
5006     //dy
5007     Double_t x0 = -params[1]/params[0];
5008     if (-params[2]*params[0]+params[1]*params[1]+1>0){
5009       Double_t Rm1  = params[0]/TMath::Sqrt(-params[2]*params[0]+params[1]*params[1]+1);
5010       if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)>0){
5011         Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0));
5012         if (params[0]<0) res*=-1.;
5013         dy = res;
5014       }
5015     }
5016     z  = rpolz0+rpolz1*(x-xref2);
5017     dz = rpolz1;
5018     cseed[iLayer].fYref[0] = y;
5019     cseed[iLayer].fYref[1] = dy;
5020     cseed[iLayer].fZref[0] = z;
5021     cseed[iLayer].fZref[1] = dz;
5022     cseed[iLayer].fC  = CR;
5023     //
5024   }
5025   return chi2TR;
5026 }