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