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 /*
  17 $Log$
  18 Revision 1.26  2003/04/10 10:36:54  hristov
  19 Code for unified TPC/TRD tracking (S.Radomski)
  20
  21 Revision 1.25  2003/03/19 17:14:11  hristov
  22 Load/UnloadClusters added to the base class and the derived classes changed correspondingly. Possibility to give 2 input files for ITS and TPC tracks in PropagateBack. TRD tracker uses fEventN from the base class (T.Kuhr)
  23
  24 Revision 1.24  2003/02/19 09:02:28  hristov
  25 Track time measurement (S.Radomski)
  26
  27 Revision 1.23  2003/02/10 14:06:10  cblume
  28 Add tracking without tilted pads as option
  29
  30 Revision 1.22  2003/01/30 15:19:58  cblume
  31 New set of  parameters
  32
  33 Revision 1.21  2003/01/27 16:34:49  cblume
  34 Update of tracking by Sergei and Chuncheng
  35
  36 Revision 1.20  2002/11/07 15:52:09  cblume
  37 Update of tracking code for tilted pads
  38
  39 Revision 1.19  2002/10/22 15:53:08  alibrary
  40 Introducing Riostream.h
  41
  42 Revision 1.18  2002/10/14 14:57:44  hristov
  43 Merging the VirtualMC branch to the main development branch (HEAD)
  44
  45 Revision 1.14.6.2  2002/07/24 10:09:31  alibrary
  46 Updating VirtualMC
  47
  48 Revision 1.17  2002/06/13 12:09:58  hristov
  49 Minor corrections
  50
  51 Revision 1.16  2002/06/12 09:54:36  cblume
  52 Update of tracking code provided by Sergei
  53
  54 Revision 1.14  2001/11/14 10:50:46  cblume
  55 Changes in digits IO. Add merging of summable digits
  56
  57 Revision 1.13  2001/05/30 12:17:47  hristov
  58 Loop variables declared once
  59
  60 Revision 1.12  2001/05/28 17:07:58  hristov
  61 Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
  62
  63 Revision 1.8  2000/12/20 13:00:44  cblume
  64 Modifications for the HP-compiler
  65
  66 Revision 1.7  2000/12/08 16:07:02  cblume
  67 Update of the tracking by Sergei
  68
  69 Revision 1.6  2000/11/30 17:38:08  cblume
  70 Changes to get in line with new STEER and EVGEN
  71
  72 Revision 1.5  2000/11/14 14:40:27  cblume
  73 Correction for the Sun compiler (kTRUE and kFALSE)
  74
  75 Revision 1.4  2000/11/10 14:57:52  cblume
  76 Changes in the geometry constants for the DEC compiler
  77
  78 Revision 1.3  2000/10/15 23:40:01  cblume
  79 Remove AliTRDconst
  80
  81 Revision 1.2  2000/10/06 16:49:46  cblume
  82 Made Getters const
  83
  84 Revision 1.1.2.2  2000/10/04 16:34:58  cblume
  85 Replace include files by forward declarations
  86
  87 Revision 1.1.2.1  2000/09/22 14:47:52  cblume
  88 Add the tracking code
  89
  90 */
  91
  92 #include <Riostream.h>
  93
  94 #include <TFile.h>
  95 #include <TBranch.h>
  96 #include <TTree.h>
  97 #include <TObjArray.h>
  98
  99 #include "AliTRDgeometry.h"
 100 #include "AliTRDparameter.h"
 101 #include "AliTRDgeometryDetail.h"
 102 #include "AliTRDcluster.h"
 103 #include "AliTRDtrack.h"
 104 #include "AliTRDPartID.h"
 105 #include "../TPC/AliTPCtrack.h"
 106
 107 #include "AliTRDtracker.h"
 108
 109 ClassImp(AliTRDtracker)
 110
 111   const  Float_t     AliTRDtracker::fSeedDepth          = 0.5;
 112   const  Float_t     AliTRDtracker::fSeedStep           = 0.10;
 113   const  Float_t     AliTRDtracker::fSeedGap            = 0.25;
 114
 115   const  Float_t     AliTRDtracker::fMaxSeedDeltaZ12    = 40.;
 116   const  Float_t     AliTRDtracker::fMaxSeedDeltaZ      = 25.;
 117   const  Float_t     AliTRDtracker::fMaxSeedC           = 0.0052;
 118   const  Float_t     AliTRDtracker::fMaxSeedTan         = 1.2;
 119   const  Float_t     AliTRDtracker::fMaxSeedVertexZ     = 150.;
 120
 121   const  Double_t    AliTRDtracker::fSeedErrorSY        = 0.2;
 122   const  Double_t    AliTRDtracker::fSeedErrorSY3       = 2.5;
 123   const  Double_t    AliTRDtracker::fSeedErrorSZ        = 0.1;
 124
 125   const  Float_t     AliTRDtracker::fMinClustersInSeed  = 0.7;
 126
 127   const  Float_t     AliTRDtracker::fMinClustersInTrack = 0.5;
 128   const  Float_t     AliTRDtracker::fMinFractionOfFoundClusters = 0.8;
 129
 130   const  Float_t     AliTRDtracker::fSkipDepth          = 0.3;
 131   const  Float_t     AliTRDtracker::fLabelFraction      = 0.8;
 132   const  Float_t     AliTRDtracker::fWideRoad           = 20.;
 133
 134   const  Double_t    AliTRDtracker::fMaxChi2            = 12.;
 135
 136 const Int_t AliTRDtracker::kFirstPlane = 5;
 137 const Int_t AliTRDtracker::kLastPlane = 17;
 138
 139
 140 //____________________________________________________________________
 141 AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
 142 {
 143   //
 144   //  Main constructor
 145   //
 146
 147   Float_t fTzero = 0;
 148
 149   fAddTRDseeds = kFALSE;
 150   fGeom = NULL;
 151   fNoTilt = kFALSE;
 152
 153   TDirectory *savedir=gDirectory;
 154   TFile *in=(TFile*)geomfile;
 155   if (!in->IsOpen()) {
 156     printf("AliTRDtracker::AliTRDtracker(): geometry file is not open!\n");
 157     printf("    DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
 158   }
 159   else {
 160     in->cd();
 161 //    in->ls();
 162     fGeom = (AliTRDgeometry*) in->Get("TRDgeometry");
 163     fPar  = (AliTRDparameter*) in->Get("TRDparameter");
 164 //    fGeom->Dump();
 165   }
 166
 167   if(fGeom) {
 168     //    fTzero = geo->GetT0();
 169     printf("Found geometry version %d on file \n", fGeom->IsVersion());
 170   }
 171   else {
 172     printf("AliTRDtracker::AliTRDtracker(): cann't find TRD geometry!\n");
 173     printf("    DETAIL TRD geometry and DEFAULT TRD parameter will be used\n");
 174     fGeom = new AliTRDgeometryDetail();
 175     fPar = new AliTRDparameter();
 176   }
 177
 178   savedir->cd();
 179
 180
 181   //  fGeom->SetT0(fTzero);
 182
 183   fNclusters = 0;
 184   fClusters  = new TObjArray(2000);
 185   fNseeds    = 0;
 186   fSeeds     = new TObjArray(2000);
 187   fNtracks   = 0;
 188   fTracks    = new TObjArray(1000);
 189
 190   for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) {
 191     Int_t tr_s = CookSectorIndex(geom_s);
 192     fTrSec[tr_s] = new AliTRDtrackingSector(fGeom, geom_s, fPar);
 193   }
 194
 195   Float_t tilt_angle = TMath::Abs(fPar->GetTiltingAngle());
 196   if(tilt_angle < 0.1) {
 197     fNoTilt = kTRUE;
 198   }
 199
 200   fSY2corr = 0.2;
 201   fSZ2corr = 120.;
 202
 203   if(fNoTilt && (tilt_angle > 0.1)) fSY2corr = fSY2corr + tilt_angle * 0.05;
 204
 205
 206   // calculate max gap on track
 207
 208   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
 209   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
 210
 211   Double_t dx = (Double_t) fPar->GetTimeBinSize();
 212   Int_t tbAmp = fPar->GetTimeBefore();
 213   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
 214   if(kTRUE) maxAmp = 0;  // intentional until we change the parameter class
 215   Int_t tbDrift = fPar->GetTimeMax();
 216   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
 217
 218   tbDrift = TMath::Min(tbDrift,maxDrift);
 219   tbAmp = TMath::Min(tbAmp,maxAmp);
 220
 221   fTimeBinsPerPlane = tbAmp + tbDrift;
 222   fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fSkipDepth);
 223
 224   fVocal = kFALSE;
 225
 226
 227   // Barrel Tracks [SR, 03.04.2003]
 228
 229   fBarrelFile = 0;
 230   fBarrelTree = 0;
 231   fBarrelArray = 0;
 232   fBarrelTrack = 0;
 233
 234   savedir->cd();
 235 }
 236
 237 //___________________________________________________________________
 238 AliTRDtracker::~AliTRDtracker()
 239 {
 240   delete fClusters;
 241   delete fTracks;
 242   delete fSeeds;
 243   delete fGeom;
 244   delete fPar;
 245
 246   for(Int_t geom_s = 0; geom_s < kTRACKING_SECTORS; geom_s++) {
 247     delete fTrSec[geom_s];
 248   }
 249 }
 250
 251 //_____________________________________________________________________
 252
 253 void AliTRDtracker::SetBarrelTree(const char *mode) {
 254   //
 255   //
 256   //
 257
 258   if (!IsStoringBarrel()) return;
 259
 260   TDirectory *sav = gDirectory;
 261   if (!fBarrelFile) fBarrelFile = new TFile("AliBarrelTracks.root", "UPDATE");
 262
 263   char buff[40];
 264   sprintf(buff,  "BarrelTRD_%d_%s", GetEventNumber(), mode);
 265
 266   fBarrelFile->cd();
 267   fBarrelTree = new TTree(buff, "Barrel TPC tracks");
 268
 269   Int_t nRefs = kLastPlane - kFirstPlane + 1;
 270
 271   if (!fBarrelArray) fBarrelArray = new TClonesArray("AliBarrelTrack", nRefs);
 272   for(Int_t i=0; i<nRefs; i++) new((*fBarrelArray)[i]) AliBarrelTrack();
 273
 274   fBarrelTree->Branch("tracks", &fBarrelArray);
 275   sav->cd();
 276 }
 277
 278 //_____________________________________________________________________
 279
 280 void AliTRDtracker::StoreBarrelTrack(AliTRDtrack *ps, Int_t refPlane, Int_t isIn) {
 281   //
 282   //
 283   //
 284
 285   if (!IsStoringBarrel()) return;
 286
 287   static Int_t nClusters;
 288   static Int_t nWrong;
 289   static Double_t chi2;
 290   static Int_t index;
 291   static Bool_t wasLast = kTRUE;
 292
 293   Int_t newClusters, newWrong;
 294   Double_t newChi2;
 295
 296   if (wasLast) {
 297
 298     fBarrelArray->Clear();
 299     nClusters = nWrong = 0;
 300     chi2 = 0.0;
 301     index = 0;
 302     wasLast = kFALSE;
 303   }
 304
 305   fBarrelTrack = (AliBarrelTrack*)(*fBarrelArray)[index++];
 306   ps->GetBarrelTrack(fBarrelTrack);
 307
 308   newClusters = ps->GetNumberOfClusters() - nClusters;
 309   newWrong = ps->GetNWrong() - nWrong;
 310   newChi2 = ps->GetChi2() - chi2;
 311
 312   nClusters =  ps->GetNumberOfClusters();
 313   nWrong = ps->GetNWrong();
 314   chi2 = ps->GetChi2();
 315
 316   if (refPlane != kLastPlane) {
 317     fBarrelTrack->SetNClusters(newClusters, newChi2);
 318     fBarrelTrack->SetNWrongClusters(newWrong);
 319   } else {
 320     wasLast = kTRUE;
 321   }
 322
 323   fBarrelTrack->SetRefPlane(refPlane, isIn);
 324 }
 325
 326 //_____________________________________________________________________
 327
 328 Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) {
 329   //
 330   // Rotates the track when necessary
 331   //
 332
 333   Double_t alpha = AliTRDgeometry::GetAlpha();
 334   Double_t y = track->GetY();
 335   Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
 336
 337   Int_t ns = AliTRDgeometry::kNsect;
 338   //Int_t s=Int_t(track->GetAlpha()/alpha)%ns;
 339
 340   if (y > ymax) {
 341     //s = (s+1) % ns;
 342     if (!track->Rotate(alpha)) return kFALSE;
 343   } else if (y <-ymax) {
 344     //s = (s-1+ns) % ns;
 345     if (!track->Rotate(-alpha)) return kFALSE;
 346   }
 347
 348   return kTRUE;
 349 }
 350
 351 //_____________________________________________________________________
 352 inline Double_t f1trd(Double_t x1,Double_t y1,
 353                       Double_t x2,Double_t y2,
 354                       Double_t x3,Double_t y3)
 355 {
 356   //
 357   // Initial approximation of the track curvature
 358   //
 359   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 360   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 361                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 362   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 363                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 364
 365   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 366
 367   return -xr*yr/sqrt(xr*xr+yr*yr);
 368 }
 369
 370 //_____________________________________________________________________
 371 inline Double_t f2trd(Double_t x1,Double_t y1,
 372                       Double_t x2,Double_t y2,
 373                       Double_t x3,Double_t y3)
 374 {
 375   //
 376   // Initial approximation of the track curvature times X coordinate
 377   // of the center of curvature
 378   //
 379
 380   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 381   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 382                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 383   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 384                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 385
 386   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 387
 388   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 389 }
 390
 391 //_____________________________________________________________________
 392 inline Double_t f3trd(Double_t x1,Double_t y1,
 393                       Double_t x2,Double_t y2,
 394                       Double_t z1,Double_t z2)
 395 {
 396   //
 397   // Initial approximation of the tangent of the track dip angle
 398   //
 399
 400   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 401 }
 402
 403 //___________________________________________________________________
 404 Int_t AliTRDtracker::Clusters2Tracks(const TFile *inp, TFile *out)
 405 {
 406   //
 407   // Finds tracks within the TRD. File <inp> is expected to contain seeds
 408   // at the outer part of the TRD. If <inp> is NULL, the seeds
 409   // are found within the TRD if fAddTRDseeds is TRUE.
 410   // The tracks are propagated to the innermost time bin
 411   // of the TRD and stored in file <out>.
 412   //
 413
 414   LoadEvent();
 415
 416   TDirectory *savedir=gDirectory;
 417
 418   char   tname[100];
 419
 420   if (!out->IsOpen()) {
 421     cerr<<"AliTRDtracker::Clusters2Tracks(): output file is not open !\n";
 422     return 1;
 423   }
 424
 425   sprintf(tname,"seedTRDtoTPC_%d",GetEventNumber());
 426   TTree tpc_tree(tname,"Tree with seeds from TRD at outer TPC pad row");
 427   AliTPCtrack *iotrack=0;
 428   tpc_tree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
 429
 430   sprintf(tname,"TreeT%d_TRD",GetEventNumber());
 431   TTree trd_tree(tname,"TRD tracks at inner TRD time bin");
 432   AliTRDtrack *iotrack_trd=0;
 433   trd_tree.Branch("tracks","AliTRDtrack",&iotrack_trd,32000,0);
 434
 435   Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 436   Float_t foundMin = fMinClustersInTrack * timeBins;
 437
 438   if (inp) {
 439      TFile *in=(TFile*)inp;
 440      if (!in->IsOpen()) {
 441         cerr<<"AliTRDtracker::Clusters2Tracks(): file with seeds is not open !\n";
 442         cerr<<" ... going for seeds finding inside the TRD\n";
 443      }
 444      else {
 445        in->cd();
 446        sprintf(tname,"TRDb_%d",GetEventNumber());
 447        TTree *seedTree=(TTree*)in->Get(tname);
 448        if (!seedTree) {
 449          cerr<<"AliTRDtracker::Clusters2Tracks(): ";
 450          cerr<<"can't get a tree with track seeds !\n";
 451          return 3;
 452        }
 453        AliTRDtrack *seed=new AliTRDtrack;
 454        seedTree->SetBranchAddress("tracks",&seed);
 455
 456        Int_t n=(Int_t)seedTree->GetEntries();
 457        for (Int_t i=0; i<n; i++) {
 458          seedTree->GetEvent(i);
 459          seed->ResetCovariance();
 460          AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
 461          fSeeds->AddLast(tr);
 462          fNseeds++;
 463        }
 464        delete seed;
 465        delete seedTree;
 466      }
 467   }
 468
 469   out->cd();
 470
 471
 472   // find tracks from loaded seeds
 473
 474   Int_t nseed=fSeeds->GetEntriesFast();
 475   Int_t i, found = 0;
 476   Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
 477
 478   for (i=0; i<nseed; i++) {
 479     AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
 480     FollowProlongation(t, innerTB);
 481     if (t.GetNumberOfClusters() >= foundMin) {
 482       UseClusters(&t);
 483       CookLabel(pt, 1-fLabelFraction);
 484       //      t.CookdEdx();
 485     }
 486     iotrack_trd = pt;
 487     trd_tree.Fill();
 488     found++;
 489 //    cout<<found<<'\r';
 490
 491     if(PropagateToTPC(t)) {
 492       AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
 493       iotrack = tpc;
 494       tpc_tree.Fill();
 495       delete tpc;
 496     }
 497     delete fSeeds->RemoveAt(i);
 498     fNseeds--;
 499   }
 500
 501   cout<<"Number of loaded seeds: "<<nseed<<endl;
 502   cout<<"Number of found tracks from loaded seeds: "<<found<<endl;
 503
 504   // after tracks from loaded seeds are found and the corresponding
 505   // clusters are used, look for additional seeds from TRD
 506
 507   if(fAddTRDseeds) {
 508     // Find tracks for the seeds in the TRD
 509     Int_t timeBins = fTrSec[0]->GetNumberOfTimeBins();
 510
 511     Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep);
 512     Int_t gap = (Int_t) (timeBins * fSeedGap);
 513     Int_t step = (Int_t) (timeBins * fSeedStep);
 514
 515     // make a first turn with tight cut on initial curvature
 516     for(Int_t turn = 1; turn <= 2; turn++) {
 517       if(turn == 2) {
 518         nSteps = (Int_t) (fSeedDepth / (3*fSeedStep));
 519         step = (Int_t) (timeBins * (3*fSeedStep));
 520       }
 521       for(Int_t i=0; i<nSteps; i++) {
 522         Int_t outer=timeBins-1-i*step;
 523         Int_t inner=outer-gap;
 524
 525         nseed=fSeeds->GetEntriesFast();
 526
 527         MakeSeeds(inner, outer, turn);
 528
 529         nseed=fSeeds->GetEntriesFast();
 530         printf("\n turn %d, step %d: number of seeds for TRD inward %d\n",
 531                turn, i, nseed);
 532
 533         for (Int_t i=0; i<nseed; i++) {
 534           AliTRDtrack *pt=(AliTRDtrack*)fSeeds->UncheckedAt(i), &t=*pt;
 535           FollowProlongation(t,innerTB);
 536           if (t.GetNumberOfClusters() >= foundMin) {
 537             UseClusters(&t);
 538             CookLabel(pt, 1-fLabelFraction);
 539             t.CookdEdx();
 540             found++;
 541 //            cout<<found<<'\r';
 542             iotrack_trd = pt;
 543             trd_tree.Fill();
 544             if(PropagateToTPC(t)) {
 545               AliTPCtrack *tpc = new AliTPCtrack(*pt,pt->GetAlpha());
 546               iotrack = tpc;
 547               tpc_tree.Fill();
 548               delete tpc;
 549             }
 550           }
 551           delete fSeeds->RemoveAt(i);
 552           fNseeds--;
 553         }
 554       }
 555     }
 556   }
 557   tpc_tree.Write();
 558   trd_tree.Write();
 559
 560   cout<<"Total number of found tracks: "<<found<<endl;
 561
 562   UnloadEvent();
 563
 564   savedir->cd();
 565
 566   return 0;
 567 }
 568
 569
 570
 571 //_____________________________________________________________________________
 572 Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) {
 573   //
 574   // Reads seeds from file <inp>. The seeds are AliTPCtrack's found and
 575   // backpropagated by the TPC tracker. Each seed is first propagated
 576   // to the TRD, and then its prolongation is searched in the TRD.
 577   // If sufficiently long continuation of the track is found in the TRD
 578   // the track is updated, otherwise it's stored as originaly defined
 579   // by the TPC tracker.
 580   //
 581
 582   LoadEvent();
 583
 584   TDirectory *savedir=gDirectory;
 585
 586   TFile *in=(TFile*)inp;
 587
 588   if (!in->IsOpen()) {
 589      cerr<<"AliTRDtracker::PropagateBack(): ";
 590      cerr<<"file with back propagated TPC tracks is not open !\n";
 591      return 1;
 592   }
 593
 594   if (!out->IsOpen()) {
 595      cerr<<"AliTRDtracker::PropagateBack(): ";
 596      cerr<<"file for back propagated TRD tracks is not open !\n";
 597      return 2;
 598   }
 599
 600   in->cd();
 601   char   tname[100];
 602   sprintf(tname,"seedsTPCtoTRD_%d",GetEventNumber());
 603   TTree *seedTree=(TTree*)in->Get(tname);
 604   if (!seedTree) {
 605      cerr<<"AliTRDtracker::PropagateBack(): ";
 606      cerr<<"can't get a tree with seeds from TPC !\n";
 607      cerr<<"check if your version of TPC tracker creates tree "<<tname<<"\n";
 608      return 3;
 609   }
 610
 611   AliTPCtrack *seed=new AliTPCtrack;
 612   seedTree->SetBranchAddress("tracks",&seed);
 613
 614   Int_t n=(Int_t)seedTree->GetEntries();
 615   for (Int_t i=0; i<n; i++) {
 616      seedTree->GetEvent(i);
 617      Int_t lbl = seed->GetLabel();
 618      AliTRDtrack *tr = new AliTRDtrack(*seed,seed->GetAlpha());
 619      tr->SetSeedLabel(lbl);
 620      fSeeds->AddLast(tr);
 621      fNseeds++;
 622   }
 623
 624   delete seed;
 625   delete seedTree;
 626
 627   out->cd();
 628
 629   AliTPCtrack *otrack=0;
 630
 631   sprintf(tname,"seedsTRDtoTOF1_%d",GetEventNumber());
 632   TTree tofTree1(tname,"Tracks back propagated through TPC and TRD");
 633   tofTree1.Branch("tracks","AliTPCtrack",&otrack,32000,0);
 634
 635   sprintf(tname,"seedsTRDtoTOF2_%d",GetEventNumber());
 636   TTree tofTree2(tname,"Tracks back propagated through TPC and TRD");
 637   tofTree2.Branch("tracks","AliTPCtrack",&otrack,32000,0);
 638
 639   sprintf(tname,"seedsTRDtoPHOS_%d",GetEventNumber());
 640   TTree phosTree(tname,"Tracks back propagated through TPC and TRD");
 641   phosTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
 642
 643   sprintf(tname,"seedsTRDtoRICH_%d",GetEventNumber());
 644   TTree richTree(tname,"Tracks back propagated through TPC and TRD");
 645   richTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
 646
 647   sprintf(tname,"TRDb_%d",GetEventNumber());
 648   TTree trdTree(tname,"Back propagated TRD tracks at outer TRD time bin");
 649   AliTRDtrack *otrack_trd=0;
 650   trdTree.Branch("tracks","AliTRDtrack",&otrack_trd,32000,0);
 651
 652   if (IsStoringBarrel()) SetBarrelTree("back");
 653   out->cd();
 654
 655   Int_t found=0;
 656   Int_t nseed=fSeeds->GetEntriesFast();
 657
 658   //  Float_t foundMin = fMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan();
 659   Float_t foundMin = 40;
 660
 661   Int_t outermost_tb  = fTrSec[0]->GetOuterTimeBin();
 662
 663   for (Int_t i=0; i<nseed; i++) {
 664
 665     AliTRDtrack *ps=(AliTRDtrack*)fSeeds->UncheckedAt(i), &s=*ps;
 666     Int_t expectedClr = FollowBackProlongation(s);
 667
 668     if (IsStoringBarrel()) {
 669       StoreBarrelTrack(ps, kLastPlane, kTrackBack);
 670       fBarrelTree->Fill();
 671     }
 672
 673     Int_t foundClr = s.GetNumberOfClusters();
 674     Int_t last_tb = fTrSec[0]->GetLayerNumber(s.GetX());
 675
 676     //    printf("seed %d: found %d out of %d expected clusters, Min is %f\n",
 677     //     i, foundClr, expectedClr, foundMin);
 678
 679     if (foundClr >= foundMin) {
 680       if(foundClr >= 2) {
 681         s.CookdEdx();
 682         CookLabel(ps, 1-fLabelFraction);
 683         UseClusters(ps);
 684       }
 685
 686       // Propagate to outer reference plane [SR, GSI, 18.02.2003]
 687       ps->PropagateTo(364.8);
 688       otrack_trd=ps;
 689       trdTree.Fill();
 690       found++;
 691 //      cout<<found<<'\r';
 692     }
 693
 694     if(((expectedClr < 10) && (last_tb == outermost_tb)) ||
 695        ((expectedClr >= 10) &&
 696         (((Float_t) foundClr) / ((Float_t) expectedClr) >=
 697          fMinFractionOfFoundClusters) && (last_tb == outermost_tb))) {
 698
 699       Double_t x_tof = 375.5;
 700
 701       if(PropagateToOuterPlane(s,x_tof)) {
 702         AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
 703         otrack = pt;
 704         tofTree1.Fill();
 705         delete pt;
 706
 707         x_tof = 381.5;
 708
 709         if(PropagateToOuterPlane(s,x_tof)) {
 710           AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
 711           otrack = pt;
 712           tofTree2.Fill();
 713           delete pt;
 714
 715           Double_t x_phos = 460.;
 716
 717           if(PropagateToOuterPlane(s,x_phos)) {
 718             AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
 719             otrack = pt;
 720             phosTree.Fill();
 721             delete pt;
 722
 723             Double_t x_rich = 490+1.267;
 724
 725             if(PropagateToOuterPlane(s,x_rich)) {
 726               AliTPCtrack *pt = new AliTPCtrack(*ps,ps->GetAlpha());
 727               otrack = pt;
 728               richTree.Fill();
 729               delete pt;
 730             }
 731           }
 732         }
 733       }
 734     }
 735   }
 736
 737
 738   out->cd();
 739   tofTree1.Write();
 740   tofTree2.Write();
 741   phosTree.Write();
 742   richTree.Write();
 743   trdTree.Write();
 744
 745   if (IsStoringBarrel()) { // [SR, 03.04.2003]
 746     fBarrelFile->cd();
 747     fBarrelTree->Write();
 748     fBarrelFile->Flush();
 749   }
 750
 751   savedir->cd();
 752   cerr<<"Number of seeds: "<<nseed<<endl;
 753   cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
 754
 755   UnloadEvent();
 756
 757   return 0;
 758
 759 }
 760
 761 //_____________________________________________________________________________
 762 Int_t AliTRDtracker::PropagateBack(AliESD* event) {
 763   //
 764   // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
 765   // backpropagated by the TPC tracker. Each seed is first propagated
 766   // to the TRD, and then its prolongation is searched in the TRD.
 767   // If sufficiently long continuation of the track is found in the TRD
 768   // the track is updated, otherwise it's stored as originaly defined
 769   // by the TPC tracker.
 770   //
 771
 772   Int_t found=0;
 773   Float_t foundMin = 40;
 774
 775   Int_t n = event->GetNumberOfTracks();
 776   for (Int_t i=0; i<n; i++) {
 777     AliESDtrack* seed=event->GetTrack(i);
 778     ULong_t status=seed->GetStatus();
 779     if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
 780     if ( (status & AliESDtrack::kTRDout) != 0 ) continue;
 781
 782     Int_t lbl = seed->GetLabel();
 783     AliTRDtrack *track = new AliTRDtrack(*seed);
 784     track->SetSeedLabel(lbl);
 785     fNseeds++;
 786
 787     Int_t expectedClr = FollowBackProlongation(*track);
 788
 789     Int_t foundClr = track->GetNumberOfClusters();
 790     if (foundClr >= foundMin) {
 791       if(foundClr >= 2) {
 792         track->CookdEdx();
 793 //      CookLabel(track, 1-fLabelFraction);
 794         UseClusters(track);
 795       }
 796
 797       // Propagate to outer reference plane [SR, GSI, 18.02.2003]
 798 //      track->PropagateTo(364.8);  why?
 799
 800       seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
 801       found++;
 802     }
 803
 804   }
 805
 806   cerr<<"Number of seeds: "<<fNseeds<<endl;
 807   cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
 808
 809   return 0;
 810
 811 }
 812
 813
 814 //---------------------------------------------------------------------------
 815 Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
 816 {
 817   // Starting from current position on track=t this function tries
 818   // to extrapolate the track up to timeBin=0 and to confirm prolongation
 819   // if a close cluster is found. Returns the number of clusters
 820   // expected to be found in sensitive layers
 821
 822   Float_t  wIndex, wTB, wChi2;
 823   Float_t  wYrt, wYclosest, wYcorrect, wYwindow;
 824   Float_t  wZrt, wZclosest, wZcorrect, wZwindow;
 825   Float_t  wPx, wPy, wPz, wC;
 826   Double_t Px, Py, Pz;
 827   Float_t  wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
 828
 829   Int_t trackIndex = t.GetLabel();
 830
 831   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
 832
 833   Int_t try_again=fMaxGap;
 834
 835   Double_t alpha=t.GetAlpha();
 836   TVector2::Phi_0_2pi(alpha);
 837
 838   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
 839   Double_t rad_length, rho, x, dx, y, ymax, z;
 840
 841   Int_t expectedNumberOfClusters = 0;
 842   Bool_t lookForCluster;
 843
 844   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
 845
 846
 847   for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
 848
 849     y = t.GetY(); z = t.GetZ();
 850
 851     // first propagate to the inner surface of the current time bin
 852     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
 853     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
 854     if(!t.PropagateTo(x,rad_length,rho)) break;
 855     y = t.GetY();
 856     ymax = x*TMath::Tan(0.5*alpha);
 857     if (y > ymax) {
 858       s = (s+1) % ns;
 859       if (!t.Rotate(alpha)) break;
 860       if(!t.PropagateTo(x,rad_length,rho)) break;
 861     } else if (y <-ymax) {
 862       s = (s-1+ns) % ns;
 863       if (!t.Rotate(-alpha)) break;
 864       if(!t.PropagateTo(x,rad_length,rho)) break;
 865     }
 866
 867     y = t.GetY(); z = t.GetZ();
 868
 869     // now propagate to the middle plane of the next time bin
 870     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
 871     x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
 872     if(!t.PropagateTo(x,rad_length,rho)) break;
 873     y = t.GetY();
 874     ymax = x*TMath::Tan(0.5*alpha);
 875     if (y > ymax) {
 876       s = (s+1) % ns;
 877       if (!t.Rotate(alpha)) break;
 878       if(!t.PropagateTo(x,rad_length,rho)) break;
 879     } else if (y <-ymax) {
 880       s = (s-1+ns) % ns;
 881       if (!t.Rotate(-alpha)) break;
 882       if(!t.PropagateTo(x,rad_length,rho)) break;
 883     }
 884
 885
 886     if(lookForCluster) {
 887
 888       expectedNumberOfClusters++;
 889       wIndex = (Float_t) t.GetLabel();
 890       wTB = nr;
 891
 892       AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr-1));
 893
 894       Double_t sy2=ExpectedSigmaY2(x,t.GetTgl(),t.GetPt());
 895       Double_t sz2=ExpectedSigmaZ2(x,t.GetTgl());
 896
 897       Double_t road;
 898       if((t.GetSigmaY2() + sy2) > 0) road=10.*sqrt(t.GetSigmaY2() + sy2);
 899       else return expectedNumberOfClusters;
 900
 901       wYrt = (Float_t) y;
 902       wZrt = (Float_t) z;
 903       wYwindow = (Float_t) road;
 904       t.GetPxPyPz(Px,Py,Pz);
 905       wPx = (Float_t) Px;
 906       wPy = (Float_t) Py;
 907       wPz = (Float_t) Pz;
 908       wC  = (Float_t) t.GetC();
 909       wSigmaC2 = (Float_t) t.GetSigmaC2();
 910       wSigmaTgl2    = (Float_t) t.GetSigmaTgl2();
 911       wSigmaY2 = (Float_t) t.GetSigmaY2();
 912       wSigmaZ2 = (Float_t) t.GetSigmaZ2();
 913       wChi2 = -1;
 914
 915       if (road>fWideRoad) {
 916         if (t.GetNumberOfClusters()>4)
 917           cerr<<t.GetNumberOfClusters()
 918               <<"FindProlongation warning: Too broad road !\n";
 919         return 0;
 920       }
 921
 922       AliTRDcluster *cl=0;
 923       UInt_t index=0;
 924
 925       Double_t max_chi2=fMaxChi2;
 926
 927       wYclosest = 12345678;
 928       wYcorrect = 12345678;
 929       wZclosest = 12345678;
 930       wZcorrect = 12345678;
 931       wZwindow  = TMath::Sqrt(2.25 * 12 * sz2);
 932
 933       // Find the closest correct cluster for debugging purposes
 934       if (time_bin) {
 935         Float_t minDY = 1000000;
 936         for (Int_t i=0; i<time_bin; i++) {
 937           AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
 938           if((c->GetLabel(0) != trackIndex) &&
 939              (c->GetLabel(1) != trackIndex) &&
 940              (c->GetLabel(2) != trackIndex)) continue;
 941           if(TMath::Abs(c->GetY() - y) > minDY) continue;
 942           minDY = TMath::Abs(c->GetY() - y);
 943           wYcorrect = c->GetY();
 944           wZcorrect = c->GetZ();
 945
 946           Double_t h01 = GetTiltFactor(c);
 947           wChi2 = t.GetPredictedChi2(c, h01);
 948         }
 949       }
 950
 951       // Now go for the real cluster search
 952
 953       if (time_bin) {
 954
 955         for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
 956           AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
 957           if (c->GetY() > y+road) break;
 958           if (c->IsUsed() > 0) continue;
 959           if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
 960
 961           Double_t h01 = GetTiltFactor(c);
 962           Double_t chi2=t.GetPredictedChi2(c,h01);
 963
 964           if (chi2 > max_chi2) continue;
 965           max_chi2=chi2;
 966           cl=c;
 967           index=time_bin.GetIndex(i);
 968         }
 969
 970         if(!cl) {
 971
 972           for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
 973             AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
 974
 975             if (c->GetY() > y+road) break;
 976             if (c->IsUsed() > 0) continue;
 977             if((c->GetZ()-z)*(c->GetZ()-z) > 12 * sz2) continue;
 978
 979             Double_t h01 = GetTiltFactor(c);
 980             Double_t chi2=t.GetPredictedChi2(c, h01);
 981
 982             if (chi2 > max_chi2) continue;
 983             max_chi2=chi2;
 984             cl=c;
 985             index=time_bin.GetIndex(i);
 986           }
 987         }
 988
 989
 990         if (cl) {
 991           wYclosest = cl->GetY();
 992           wZclosest = cl->GetZ();
 993           Double_t h01 = GetTiltFactor(cl);
 994
 995           t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
 996           if(!t.Update(cl,max_chi2,index,h01)) {
 997             if(!try_again--) return 0;
 998           }
 999           else try_again=fMaxGap;
1000         }
1001         else {
1002           if (try_again==0) break;
1003           try_again--;
1004         }
1005
1006         /*
1007         if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
1008
1009           printf(" %f", wIndex);       //1
1010           printf(" %f", wTB);          //2
1011           printf(" %f", wYrt);         //3
1012           printf(" %f", wYclosest);    //4
1013           printf(" %f", wYcorrect);    //5
1014           printf(" %f", wYwindow);     //6
1015           printf(" %f", wZrt);         //7
1016           printf(" %f", wZclosest);    //8
1017           printf(" %f", wZcorrect);    //9
1018           printf(" %f", wZwindow);     //10
1019           printf(" %f", wPx);          //11
1020           printf(" %f", wPy);          //12
1021           printf(" %f", wPz);          //13
1022           printf(" %f", wSigmaC2*1000000);  //14
1023           printf(" %f", wSigmaTgl2*1000);   //15
1024           printf(" %f", wSigmaY2);     //16
1025           //      printf(" %f", wSigmaZ2);     //17
1026           printf(" %f", wChi2);     //17
1027           printf(" %f", wC);           //18
1028           printf("\n");
1029         }
1030         */
1031       }
1032     }
1033   }
1034   return expectedNumberOfClusters;
1035
1036
1037 }
1038
1039 //___________________________________________________________________
1040
1041 Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
1042 {
1043   // Starting from current radial position of track <t> this function
1044   // extrapolates the track up to outer timebin and in the sensitive
1045   // layers confirms prolongation if a close cluster is found.
1046   // Returns the number of clusters expected to be found in sensitive layers
1047
1048   Float_t  wIndex, wTB, wChi2;
1049   Float_t  wYrt, wYclosest, wYcorrect, wYwindow;
1050   Float_t  wZrt, wZclosest, wZcorrect, wZwindow;
1051   Float_t  wPx, wPy, wPz, wC;
1052   Double_t Px, Py, Pz;
1053   Float_t  wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
1054
1055   Int_t trackIndex = t.GetLabel();
1056   Int_t try_again=fMaxGap;
1057
1058   Double_t alpha=t.GetAlpha();
1059   TVector2::Phi_0_2pi(alpha);
1060
1061   Int_t s;
1062
1063   Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
1064   Double_t rad_length, rho, x, dx, y, ymax = 0, z;
1065   Bool_t lookForCluster;
1066
1067   Int_t expectedNumberOfClusters = 0;
1068   x = t.GetX();
1069
1070   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1071
1072   Int_t nRefPlane = kFirstPlane;
1073   Bool_t isNewLayer = kFALSE;
1074
1075   Double_t chi2;
1076   Double_t minDY;
1077
1078   for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB+1; nr++) {
1079
1080     y = t.GetY();
1081     z = t.GetZ();
1082
1083     // first propagate to the outer surface of the current time bin
1084
1085     s = t.GetSector();
1086     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
1087     x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2;
1088     y = t.GetY();
1089     z = t.GetZ();
1090
1091     if(!t.PropagateTo(x,rad_length,rho)) break;
1092     if (!AdjustSector(&t)) break;
1093     s = t.GetSector();
1094     if (!t.PropagateTo(x,rad_length,rho)) break;
1095
1096     y = t.GetY();
1097     z = t.GetZ();
1098
1099     // Barrel Tracks [SR, 04.04.2003]
1100
1101     s = t.GetSector();
1102     if (fTrSec[s]->GetLayer(nr)->IsSensitive() !=
1103         fTrSec[s]->GetLayer(nr+1)->IsSensitive() ) {
1104
1105 //      if (IsStoringBarrel()) StoreBarrelTrack(&t, nRefPlane++, kTrackBack);
1106     }
1107
1108     if (fTrSec[s]->GetLayer(nr-1)->IsSensitive() &&
1109           ! fTrSec[s]->GetLayer(nr)->IsSensitive()) {
1110       isNewLayer = kTRUE;
1111     } else {isNewLayer = kFALSE;}
1112
1113     y = t.GetY();
1114     z = t.GetZ();
1115
1116     // now propagate to the middle plane of the next time bin
1117     fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
1118
1119     x = fTrSec[s]->GetLayer(nr+1)->GetX();
1120       if(!t.PropagateTo(x,rad_length,rho)) break;
1121     if (!AdjustSector(&t)) break;
1122     s = t.GetSector();
1123       if(!t.PropagateTo(x,rad_length,rho)) break;
1124
1125     y = t.GetY();
1126     z = t.GetZ();
1127
1128     if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax);
1129     //    printf("label %d, pl %d, lookForCluster %d \n",
1130     //     trackIndex, nr+1, lookForCluster);
1131
1132     if(lookForCluster) {
1133       expectedNumberOfClusters++;
1134
1135       wIndex = (Float_t) t.GetLabel();
1136       wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex();
1137
1138       AliTRDpropagationLayer& time_bin=*(fTrSec[s]->GetLayer(nr+1));
1139       Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
1140       Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
1141       if((t.GetSigmaY2() + sy2) < 0) break;
1142       Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2);
1143       Double_t y=t.GetY(), z=t.GetZ();
1144
1145       wYrt = (Float_t) y;
1146       wZrt = (Float_t) z;
1147       wYwindow = (Float_t) road;
1148       t.GetPxPyPz(Px,Py,Pz);
1149       wPx = (Float_t) Px;
1150       wPy = (Float_t) Py;
1151       wPz = (Float_t) Pz;
1152       wC  = (Float_t) t.GetC();
1153       wSigmaC2 = (Float_t) t.GetSigmaC2();
1154       wSigmaTgl2    = (Float_t) t.GetSigmaTgl2();
1155       wSigmaY2 = (Float_t) t.GetSigmaY2();
1156       wSigmaZ2 = (Float_t) t.GetSigmaZ2();
1157       wChi2 = -1;
1158
1159       if (road>fWideRoad) {
1160         if (t.GetNumberOfClusters()>4)
1161           cerr<<t.GetNumberOfClusters()
1162               <<"FindProlongation warning: Too broad road !\n";
1163         return 0;
1164       }
1165
1166       AliTRDcluster *cl=0;
1167       UInt_t index=0;
1168
1169       Double_t max_chi2=fMaxChi2;
1170
1171       if (isNewLayer) {
1172         road = 3 * road;
1173         //sz2 = 3 * sz2;
1174         max_chi2 = 10 * fMaxChi2;
1175       }
1176
1177       if (nRefPlane == kFirstPlane) max_chi2 = 20 * fMaxChi2;
1178       if (nRefPlane == kFirstPlane+2) max_chi2 = 15 * fMaxChi2;
1179       if (t.GetNRotate() > 0) max_chi2 = 3 * max_chi2;
1180
1181
1182       wYclosest = 12345678;
1183       wYcorrect = 12345678;
1184       wZclosest = 12345678;
1185       wZcorrect = 12345678;
1186       wZwindow  = TMath::Sqrt(2.25 * 12 * sz2);
1187
1188       // Find the closest correct cluster for debugging purposes
1189       if (time_bin) {
1190         minDY = 1000000;
1191         for (Int_t i=0; i<time_bin; i++) {
1192           AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
1193           if((c->GetLabel(0) != trackIndex) &&
1194              (c->GetLabel(1) != trackIndex) &&
1195              (c->GetLabel(2) != trackIndex)) continue;
1196           if(TMath::Abs(c->GetY() - y) > minDY) continue;
1197           //minDY = TMath::Abs(c->GetY() - y);
1198           minDY = c->GetY() - y;
1199           wYcorrect = c->GetY();
1200           wZcorrect = c->GetZ();
1201
1202           Double_t h01 = GetTiltFactor(c);
1203           wChi2 = t.GetPredictedChi2(c, h01);
1204         }
1205       }
1206
1207       // Now go for the real cluster search
1208
1209       if (time_bin) {
1210
1211         for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
1212           AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
1213           if (c->GetY() > y+road) break;
1214           if (c->IsUsed() > 0) continue;
1215           if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
1216
1217           Double_t h01 = GetTiltFactor(c);
1218           chi2=t.GetPredictedChi2(c,h01);
1219
1220           if (chi2 > max_chi2) continue;
1221           max_chi2=chi2;
1222           cl=c;
1223           index=time_bin.GetIndex(i);
1224
1225           //check is correct
1226           if((c->GetLabel(0) != trackIndex) &&
1227              (c->GetLabel(1) != trackIndex) &&
1228              (c->GetLabel(2) != trackIndex)) t.AddNWrong();
1229         }
1230
1231         if(!cl) {
1232
1233           for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
1234             AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
1235
1236             if (c->GetY() > y+road) break;
1237             if (c->IsUsed() > 0) continue;
1238             if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
1239
1240             Double_t h01 = GetTiltFactor(c);
1241             chi2=t.GetPredictedChi2(c,h01);
1242
1243             if (chi2 > max_chi2) continue;
1244             max_chi2=chi2;
1245             cl=c;
1246             index=time_bin.GetIndex(i);
1247           }
1248         }
1249
1250         if (cl) {
1251           wYclosest = cl->GetY();
1252           wZclosest = cl->GetZ();
1253
1254           t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters());
1255           Double_t h01 = GetTiltFactor(cl);
1256           if(!t.Update(cl,max_chi2,index,h01)) {
1257             if(!try_again--) return 0;
1258           }
1259           else try_again=fMaxGap;
1260         }
1261         else {
1262           if (try_again==0) break;
1263           try_again--;
1264
1265           //if (minDY < 1000000 && isNewLayer)
1266             //cout << "\t" << nRefPlane << "\t" << "\t" << t.GetNRotate() <<  "\t" <<
1267             //  road << "\t" << minDY << "\t" << chi2 << "\t" << wChi2 << "\t" << max_chi2 << endl;
1268
1269         }
1270
1271         isNewLayer = kFALSE;
1272
1273         /*
1274         if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
1275
1276           printf(" %f", wIndex);       //1
1277           printf(" %f", wTB);          //2
1278           printf(" %f", wYrt);         //3
1279           printf(" %f", wYclosest);    //4
1280           printf(" %f", wYcorrect);    //5
1281           printf(" %f", wYwindow);     //6
1282           printf(" %f", wZrt);         //7
1283           printf(" %f", wZclosest);    //8
1284           printf(" %f", wZcorrect);    //9
1285           printf(" %f", wZwindow);     //10
1286           printf(" %f", wPx);          //11
1287           printf(" %f", wPy);          //12
1288           printf(" %f", wPz);          //13
1289           printf(" %f", wSigmaC2*1000000);  //14
1290           printf(" %f", wSigmaTgl2*1000);   //15
1291           printf(" %f", wSigmaY2);     //16
1292           //      printf(" %f", wSigmaZ2);     //17
1293           printf(" %f", wChi2);     //17
1294           printf(" %f", wC);           //18
1295           printf("\n");
1296         }
1297         */
1298       }
1299     }
1300   }
1301   return expectedNumberOfClusters;
1302
1303
1304 }
1305
1306 //___________________________________________________________________
1307
1308 Int_t AliTRDtracker::PropagateToOuterPlane(AliTRDtrack& t, Double_t xToGo)
1309 {
1310   // Starting from current radial position of track <t> this function
1311   // extrapolates the track up to radial position <xToGo>.
1312   // Returns 1 if track reaches the plane, and 0 otherwise
1313
1314   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1315
1316   Double_t alpha=t.GetAlpha();
1317
1318   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1319   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1320
1321   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1322
1323   Bool_t lookForCluster;
1324   Double_t rad_length, rho, x, dx, y, ymax, z;
1325
1326   x = t.GetX();
1327
1328   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1329
1330   Int_t plToGo = fTrSec[0]->GetLayerNumber(xToGo);
1331
1332   for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr<plToGo; nr++) {
1333
1334     y = t.GetY(); z = t.GetZ();
1335
1336     // first propagate to the outer surface of the current time bin
1337     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
1338     x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
1339     if(!t.PropagateTo(x,rad_length,rho)) return 0;
1340     y = t.GetY();
1341     ymax = x*TMath::Tan(0.5*alpha);
1342     if (y > ymax) {
1343       s = (s+1) % ns;
1344       if (!t.Rotate(alpha)) return 0;
1345     } else if (y <-ymax) {
1346       s = (s-1+ns) % ns;
1347       if (!t.Rotate(-alpha)) return 0;
1348     }
1349     if(!t.PropagateTo(x,rad_length,rho)) return 0;
1350
1351     y = t.GetY(); z = t.GetZ();
1352
1353     // now propagate to the middle plane of the next time bin
1354     fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
1355     x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
1356     if(!t.PropagateTo(x,rad_length,rho)) return 0;
1357     y = t.GetY();
1358     ymax = x*TMath::Tan(0.5*alpha);
1359     if (y > ymax) {
1360       s = (s+1) % ns;
1361       if (!t.Rotate(alpha)) return 0;
1362     } else if (y <-ymax) {
1363       s = (s-1+ns) % ns;
1364       if (!t.Rotate(-alpha)) return 0;
1365     }
1366     if(!t.PropagateTo(x,rad_length,rho)) return 0;
1367   }
1368   return 1;
1369 }
1370
1371 //___________________________________________________________________
1372
1373 Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t)
1374 {
1375   // Starting from current radial position of track <t> this function
1376   // extrapolates the track up to radial position of the outermost
1377   // padrow of the TPC.
1378   // Returns 1 if track reaches the TPC, and 0 otherwise
1379
1380   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);
1381
1382   Double_t alpha=t.GetAlpha();
1383   TVector2::Phi_0_2pi(alpha);
1384
1385   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
1386
1387   Bool_t lookForCluster;
1388   Double_t rad_length, rho, x, dx, y, ymax, z;
1389
1390   x = t.GetX();
1391
1392   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
1393   Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055);
1394
1395   for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) {
1396
1397     y = t.GetY();
1398     z = t.GetZ();
1399
1400     // first propagate to the outer surface of the current time bin
1401     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
1402     x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2;
1403
1404     if(!t.PropagateTo(x,rad_length,rho)) return 0;
1405     AdjustSector(&t);
1406     if(!t.PropagateTo(x,rad_length,rho)) return 0;
1407
1408     y = t.GetY();
1409     z = t.GetZ();
1410
1411     // now propagate to the middle plane of the next time bin
1412     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
1413     x = fTrSec[s]->GetLayer(nr-1)->GetX();
1414
1415     if(!t.PropagateTo(x,rad_length,rho)) return 0;
1416     AdjustSector(&t);
1417     if(!t.PropagateTo(x,rad_length,rho)) return 0;
1418   }
1419   return 1;
1420 }
1421
1422 //_____________________________________________________________________________
1423 void AliTRDtracker::LoadEvent()
1424 {
1425   // Fills clusters into TRD tracking_sectors
1426   // Note that the numbering scheme for the TRD tracking_sectors
1427   // differs from that of TRD sectors
1428
1429   ReadClusters(fClusters);
1430   Int_t ncl=fClusters->GetEntriesFast();
1431   cout<<"\n LoadSectors: sorting "<<ncl<<" clusters"<<endl;
1432
1433   UInt_t index;
1434   while (ncl--) {
1435 //    printf("\r %d left  ",ncl);
1436     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
1437     Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin();
1438     Int_t sector=fGeom->GetSector(detector);
1439     Int_t plane=fGeom->GetPlane(detector);
1440
1441     Int_t tracking_sector = CookSectorIndex(sector);
1442
1443     Int_t gtb = fTrSec[tracking_sector]->CookTimeBinIndex(plane,local_time_bin);
1444     if(gtb < 0) continue;
1445     Int_t layer = fTrSec[tracking_sector]->GetLayerNumber(gtb);
1446
1447     index=ncl;
1448     fTrSec[tracking_sector]->GetLayer(layer)->InsertCluster(c,index);
1449   }
1450   printf("\r\n");
1451
1452 }
1453
1454 //_____________________________________________________________________________
1455 void AliTRDtracker::UnloadEvent()
1456 {
1457   //
1458   // Clears the arrays of clusters and tracks. Resets sectors and timebins
1459   //
1460
1461   Int_t i, nentr;
1462
1463   nentr = fClusters->GetEntriesFast();
1464   for (i = 0; i < nentr; i++) delete fClusters->RemoveAt(i);
1465
1466   nentr = fSeeds->GetEntriesFast();
1467   for (i = 0; i < nentr; i++) delete fSeeds->RemoveAt(i);
1468
1469   nentr = fTracks->GetEntriesFast();
1470   for (i = 0; i < nentr; i++) delete fTracks->RemoveAt(i);
1471
1472   Int_t nsec = AliTRDgeometry::kNsect;
1473
1474   for (i = 0; i < nsec; i++) {
1475     for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
1476       fTrSec[i]->GetLayer(pl)->Clear();
1477     }
1478   }
1479
1480 }
1481
1482 //__________________________________________________________________________
1483 void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn)
1484 {
1485   // Creates track seeds using clusters in timeBins=i1,i2
1486
1487   if(turn > 2) {
1488     cerr<<"MakeSeeds: turn "<<turn<<" exceeds the limit of 2"<<endl;
1489     return;
1490   }
1491
1492   Double_t x[5], c[15];
1493   Int_t max_sec=AliTRDgeometry::kNsect;
1494
1495   Double_t alpha=AliTRDgeometry::GetAlpha();
1496   Double_t shift=AliTRDgeometry::GetAlpha()/2.;
1497   Double_t cs=cos(alpha), sn=sin(alpha);
1498   Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);
1499
1500
1501   Int_t i2 = fTrSec[0]->GetLayerNumber(inner);
1502   Int_t i1 = fTrSec[0]->GetLayerNumber(outer);
1503
1504   Double_t x1 =fTrSec[0]->GetX(i1);
1505   Double_t xx2=fTrSec[0]->GetX(i2);
1506
1507   for (Int_t ns=0; ns<max_sec; ns++) {
1508
1509     Int_t nl2 = *(fTrSec[(ns-2+max_sec)%max_sec]->GetLayer(i2));
1510     Int_t nl=(*fTrSec[(ns-1+max_sec)%max_sec]->GetLayer(i2));
1511     Int_t nm=(*fTrSec[ns]->GetLayer(i2));
1512     Int_t nu=(*fTrSec[(ns+1)%max_sec]->GetLayer(i2));
1513     Int_t nu2=(*fTrSec[(ns+2)%max_sec]->GetLayer(i2));
1514
1515     AliTRDpropagationLayer& r1=*(fTrSec[ns]->GetLayer(i1));
1516
1517     for (Int_t is=0; is < r1; is++) {
1518       Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
1519
1520       for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
1521
1522         const AliTRDcluster *cl;
1523         Double_t x2,   y2,   z2;
1524         Double_t x3=0., y3=0.;
1525
1526         if (js<nl2) {
1527           if(turn != 2) continue;
1528           AliTRDpropagationLayer& r2=*(fTrSec[(ns-2+max_sec)%max_sec]->GetLayer(i2));
1529           cl=r2[js];
1530           y2=cl->GetY(); z2=cl->GetZ();
1531
1532           x2= xx2*cs2+y2*sn2;
1533           y2=-xx2*sn2+y2*cs2;
1534         }
1535         else if (js<nl2+nl) {
1536           if(turn != 1) continue;
1537           AliTRDpropagationLayer& r2=*(fTrSec[(ns-1+max_sec)%max_sec]->GetLayer(i2));
1538           cl=r2[js-nl2];
1539           y2=cl->GetY(); z2=cl->GetZ();
1540
1541           x2= xx2*cs+y2*sn;
1542           y2=-xx2*sn+y2*cs;
1543         }
1544         else if (js<nl2+nl+nm) {
1545           if(turn != 1) continue;
1546           AliTRDpropagationLayer& r2=*(fTrSec[ns]->GetLayer(i2));
1547           cl=r2[js-nl2-nl];
1548           x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
1549         }
1550         else if (js<nl2+nl+nm+nu) {
1551           if(turn != 1) continue;
1552           AliTRDpropagationLayer& r2=*(fTrSec[(ns+1)%max_sec]->GetLayer(i2));
1553           cl=r2[js-nl2-nl-nm];
1554           y2=cl->GetY(); z2=cl->GetZ();
1555
1556           x2=xx2*cs-y2*sn;
1557           y2=xx2*sn+y2*cs;
1558         }
1559         else {
1560           if(turn != 2) continue;
1561           AliTRDpropagationLayer& r2=*(fTrSec[(ns+2)%max_sec]->GetLayer(i2));
1562           cl=r2[js-nl2-nl-nm-nu];
1563           y2=cl->GetY(); z2=cl->GetZ();
1564
1565           x2=xx2*cs2-y2*sn2;
1566           y2=xx2*sn2+y2*cs2;
1567         }
1568
1569         if(TMath::Abs(z1-z2) > fMaxSeedDeltaZ12) continue;
1570
1571         Double_t zz=z1 - z1/x1*(x1-x2);
1572
1573         if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue;
1574
1575         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1576         if (d==0.) {cerr<<"TRD MakeSeeds: Straight seed !\n"; continue;}
1577
1578         x[0]=y1;
1579         x[1]=z1;
1580         x[4]=f1trd(x1,y1,x2,y2,x3,y3);
1581
1582         if (TMath::Abs(x[4]) > fMaxSeedC) continue;
1583
1584         x[2]=f2trd(x1,y1,x2,y2,x3,y3);
1585
1586         if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1587
1588         x[3]=f3trd(x1,y1,x2,y2,z1,z2);
1589
1590         if (TMath::Abs(x[3]) > fMaxSeedTan) continue;
1591
1592         Double_t a=asin(x[2]);
1593         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
1594
1595         if (TMath::Abs(zv)>fMaxSeedVertexZ) continue;
1596
1597         Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
1598         Double_t sy2=cl->GetSigmaY2(),     sz2=cl->GetSigmaZ2();
1599         Double_t sy3=fSeedErrorSY3, sy=fSeedErrorSY, sz=fSeedErrorSZ;
1600
1601         // Tilt changes
1602         Double_t h01 = GetTiltFactor(r1[is]);
1603         Double_t xu_factor = 100.;
1604         if(fNoTilt) {
1605           h01 = 0;
1606           xu_factor = 1;
1607         }
1608
1609         sy1=sy1+sz1*h01*h01;
1610         Double_t syz=sz1*(-h01);
1611         // end of tilt changes
1612
1613         Double_t f40=(f1trd(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1614         Double_t f42=(f1trd(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1615         Double_t f43=(f1trd(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1616         Double_t f20=(f2trd(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1617         Double_t f22=(f2trd(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1618         Double_t f23=(f2trd(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1619         Double_t f30=(f3trd(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1620         Double_t f31=(f3trd(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1621         Double_t f32=(f3trd(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1622         Double_t f34=(f3trd(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1623
1624
1625         c[0]=sy1;
1626         //        c[1]=0.;       c[2]=sz1;
1627         c[1]=syz;       c[2]=sz1*xu_factor;
1628         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1629         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
1630                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1631         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1632         c[13]=f30*sy1*f40+f32*sy2*f42;
1633         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1634
1635         UInt_t index=r1.GetIndex(is);
1636
1637         AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift);
1638
1639         Int_t rc=FollowProlongation(*track, i2);
1640
1641         if ((rc < 1) ||
1642             (track->GetNumberOfClusters() <
1643              (outer-inner)*fMinClustersInSeed)) delete track;
1644         else {
1645           fSeeds->AddLast(track); fNseeds++;
1646 //          cerr<<"\r found seed "<<fNseeds;
1647         }
1648       }
1649     }
1650   }
1651 }
1652
1653 //_____________________________________________________________________________
1654 void AliTRDtracker::ReadClusters(TObjArray *array, const TFile *inp)
1655 {
1656   //
1657   // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0)
1658   // from the file. The names of the cluster tree and branches
1659   // should match the ones used in AliTRDclusterizer::WriteClusters()
1660   //
1661
1662   TDirectory *savedir=gDirectory;
1663
1664   if (inp) {
1665      TFile *in=(TFile*)inp;
1666      if (!in->IsOpen()) {
1667         cerr<<"AliTRDtracker::ReadClusters(): input file is not open !\n";
1668         return;
1669      }
1670      else{
1671        in->cd();
1672      }
1673   }
1674
1675   Char_t treeName[12];
1676   sprintf(treeName,"TreeR%d_TRD",GetEventNumber());
1677   TTree *ClusterTree = (TTree*) gDirectory->Get(treeName);
1678
1679   TObjArray *ClusterArray = new TObjArray(400);
1680
1681   ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray);
1682
1683   Int_t nEntries = (Int_t) ClusterTree->GetEntries();
1684   printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
1685
1686   // Loop through all entries in the tree
1687   Int_t nbytes;
1688   AliTRDcluster *c = 0;
1689   printf("\n");
1690
1691   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
1692
1693     // Import the tree
1694     nbytes += ClusterTree->GetEvent(iEntry);
1695
1696     // Get the number of points in the detector
1697     Int_t nCluster = ClusterArray->GetEntriesFast();
1698 //    printf("\r Read %d clusters from entry %d", nCluster, iEntry);
1699
1700     // Loop through all TRD digits
1701     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
1702       c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
1703       AliTRDcluster *co = new AliTRDcluster(*c);
1704       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
1705       Int_t ltb = co->GetLocalTimeBin();
1706       if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
1707       else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
1708       array->AddLast(co);
1709       delete ClusterArray->RemoveAt(iCluster);
1710     }
1711   }
1712
1713   delete ClusterArray;
1714   savedir->cd();
1715
1716 }
1717
1718 //______________________________________________________________________
1719 void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename)
1720 {
1721   //
1722   // Reads AliTRDclusters from file <filename>. The names of the cluster
1723   // tree and branches should match the ones used in
1724   // AliTRDclusterizer::WriteClusters()
1725   // if <array> == 0, clusters are added into AliTRDtracker fCluster array
1726   //
1727
1728   TDirectory *savedir=gDirectory;
1729
1730   TFile *file = TFile::Open(filename);
1731   if (!file->IsOpen()) {
1732     cerr<<"Can't open file with TRD clusters"<<endl;
1733     return;
1734   }
1735
1736   Char_t treeName[12];
1737   sprintf(treeName,"TreeR%d_TRD",GetEventNumber());
1738   TTree *ClusterTree = (TTree*) gDirectory->Get(treeName);
1739
1740   if (!ClusterTree) {
1741      cerr<<"AliTRDtracker::ReadClusters(): ";
1742      cerr<<"can't get a tree with clusters !\n";
1743      return;
1744   }
1745
1746   TObjArray *ClusterArray = new TObjArray(400);
1747
1748   ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray);
1749
1750   Int_t nEntries = (Int_t) ClusterTree->GetEntries();
1751   cout<<"found "<<nEntries<<" in ClusterTree"<<endl;
1752
1753   // Loop through all entries in the tree
1754   Int_t nbytes;
1755   AliTRDcluster *c = 0;
1756
1757   printf("\n");
1758
1759   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
1760
1761     // Import the tree
1762     nbytes += ClusterTree->GetEvent(iEntry);
1763
1764     // Get the number of points in the detector
1765     Int_t nCluster = ClusterArray->GetEntriesFast();
1766     printf("\n Read %d clusters from entry %d", nCluster, iEntry);
1767
1768     // Loop through all TRD digits
1769     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) {
1770       c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
1771       AliTRDcluster *co = new AliTRDcluster(*c);
1772       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
1773       Int_t ltb = co->GetLocalTimeBin();
1774       if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2());
1775       else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr);
1776       array->AddLast(co);
1777       delete ClusterArray->RemoveAt(iCluster);
1778     }
1779   }
1780
1781   file->Close();
1782   delete ClusterArray;
1783   savedir->cd();
1784
1785 }
1786
1787
1788 //__________________________________________________________________
1789 void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const {
1790
1791   Int_t label=123456789, index, i, j;
1792   Int_t ncl=pt->GetNumberOfClusters();
1793   const Int_t range = fTrSec[0]->GetOuterTimeBin()+1;
1794
1795   Bool_t label_added;
1796
1797   //  Int_t s[range][2];
1798   Int_t **s = new Int_t* [range];
1799   for (i=0; i<range; i++) {
1800     s[i] = new Int_t[2];
1801   }
1802   for (i=0; i<range; i++) {
1803     s[i][0]=-1;
1804     s[i][1]=0;
1805   }
1806
1807   Int_t t0,t1,t2;
1808   for (i=0; i<ncl; i++) {
1809     index=pt->GetClusterIndex(i);
1810     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
1811     t0=c->GetLabel(0);
1812     t1=c->GetLabel(1);
1813     t2=c->GetLabel(2);
1814   }
1815
1816   for (i=0; i<ncl; i++) {
1817     index=pt->GetClusterIndex(i);
1818     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
1819     for (Int_t k=0; k<3; k++) {
1820       label=c->GetLabel(k);
1821       label_added=kFALSE; j=0;
1822       if (label >= 0) {
1823         while ( (!label_added) && ( j < range ) ) {
1824           if (s[j][0]==label || s[j][1]==0) {
1825             s[j][0]=label;
1826             s[j][1]=s[j][1]+1;
1827             label_added=kTRUE;
1828           }
1829           j++;
1830         }
1831       }
1832     }
1833   }
1834
1835   Int_t max=0;
1836   label = -123456789;
1837
1838   for (i=0; i<range; i++) {
1839     if (s[i][1]>max) {
1840       max=s[i][1]; label=s[i][0];
1841     }
1842   }
1843
1844   for (i=0; i<range; i++) {
1845     delete []s[i];
1846   }
1847
1848   delete []s;
1849
1850   if ((1.- Float_t(max)/ncl) > wrong) label=-label;
1851
1852   pt->SetLabel(label);
1853
1854 }
1855
1856
1857 //__________________________________________________________________
1858 void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const {
1859   Int_t ncl=t->GetNumberOfClusters();
1860   for (Int_t i=from; i<ncl; i++) {
1861     Int_t index = t->GetClusterIndex(i);
1862     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(index);
1863     c->Use();
1864   }
1865 }
1866
1867
1868 //_____________________________________________________________________
1869 Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt)
1870 {
1871   // Parametrised "expected" error of the cluster reconstruction in Y
1872
1873   Double_t s = 0.08 * 0.08;
1874   return s;
1875 }
1876
1877 //_____________________________________________________________________
1878 Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl)
1879 {
1880   // Parametrised "expected" error of the cluster reconstruction in Z
1881
1882   Double_t s = 9 * 9 /12.;
1883   return s;
1884 }
1885
1886
1887 //_____________________________________________________________________
1888 Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t local_tb) const
1889 {
1890   //
1891   // Returns radial position which corresponds to time bin <local_tb>
1892   // in tracking sector <sector> and plane <plane>
1893   //
1894
1895   Int_t index = fTrSec[sector]->CookTimeBinIndex(plane, local_tb);
1896   Int_t pl = fTrSec[sector]->GetLayerNumber(index);
1897   return fTrSec[sector]->GetLayer(pl)->GetX();
1898
1899 }
1900
1901
1902 //_______________________________________________________
1903 AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x,
1904                Double_t dx, Double_t rho, Double_t rad_length, Int_t tb_index)
1905 {
1906   //
1907   // AliTRDpropagationLayer constructor
1908   //
1909
1910   fN = 0; fX = x; fdX = dx; fRho = rho; fX0 = rad_length;
1911   fClusters = NULL; fIndex = NULL; fTimeBinIndex = tb_index;
1912
1913
1914   for(Int_t i=0; i < (Int_t) kZONES; i++) {
1915     fZc[i]=0; fZmax[i] = 0;
1916   }
1917
1918   fYmax = 0;
1919
1920   if(fTimeBinIndex >= 0) {
1921     fClusters = new AliTRDcluster*[kMAX_CLUSTER_PER_TIME_BIN];
1922     fIndex = new UInt_t[kMAX_CLUSTER_PER_TIME_BIN];
1923   }
1924
1925   fHole = kFALSE;
1926   fHoleZc = 0;
1927   fHoleZmax = 0;
1928   fHoleYc = 0;
1929   fHoleYmax = 0;
1930   fHoleRho = 0;
1931   fHoleX0 = 0;
1932
1933 }
1934
1935 //_______________________________________________________
1936 void AliTRDtracker::AliTRDpropagationLayer::SetHole(
1937           Double_t Zmax, Double_t Ymax, Double_t rho,
1938           Double_t rad_length, Double_t Yc, Double_t Zc)
1939 {
1940   //
1941   // Sets hole in the layer
1942   //
1943
1944   fHole = kTRUE;
1945   fHoleZc = Zc;
1946   fHoleZmax = Zmax;
1947   fHoleYc = Yc;
1948   fHoleYmax = Ymax;
1949   fHoleRho = rho;
1950   fHoleX0 = rad_length;
1951 }
1952
1953
1954 //_______________________________________________________
1955 AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, AliTRDparameter* par)
1956 {
1957   //
1958   // AliTRDtrackingSector Constructor
1959   //
1960
1961   fGeom = geo;
1962   fPar = par;
1963   fGeomSector = gs;
1964   fTzeroShift = 0.13;
1965   fN = 0;
1966
1967   for(UInt_t i=0; i < kMAX_TIME_BIN_INDEX; i++) fTimeBinIndex[i] = -1;
1968
1969
1970   AliTRDpropagationLayer* ppl;
1971
1972   Double_t x, xin, xout, dx, rho, rad_length;
1973   Int_t    steps;
1974
1975   // set time bins in the gas of the TPC
1976
1977   xin = 246.055; xout = 254.055; steps = 20; dx = (xout-xin)/steps;
1978   rho = 0.9e-3;  rad_length = 28.94;
1979
1980   for(Int_t i=0; i<steps; i++) {
1981     x = xin + i*dx + dx/2;
1982     ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
1983     InsertLayer(ppl);
1984   }
1985
1986   // set time bins in the outer field cage vessel
1987
1988   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; rad_length = 44.77; // Tedlar
1989   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1990   InsertLayer(ppl);
1991
1992   dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; rad_length = 44.86; // prepreg
1993   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
1994   InsertLayer(ppl);
1995
1996   dx = 2.; xin = xout; xout = xin + dx; rho = 1.45*0.02; rad_length = 41.28; // Nomex
1997   steps = 5; dx = (xout - xin)/steps;
1998   for(Int_t i=0; i<steps; i++) {
1999     x = xin + i*dx + dx/2;
2000     ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
2001     InsertLayer(ppl);
2002   }
2003
2004   dx = 0.02; xin = xout; xout = xin + dx; rho = 1.45; rad_length = 44.86; // prepreg
2005   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
2006   InsertLayer(ppl);
2007
2008   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; rad_length = 44.77; // Tedlar
2009   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
2010   InsertLayer(ppl);
2011
2012
2013   // set time bins in CO2
2014
2015   xin = xout; xout = 275.0;
2016   steps = 50; dx = (xout - xin)/steps;
2017   rho = 1.977e-3;  rad_length = 36.2;
2018
2019   for(Int_t i=0; i<steps; i++) {
2020     x = xin + i*dx + dx/2;
2021     ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
2022     InsertLayer(ppl);
2023   }
2024
2025   // set time bins in the outer containment vessel
2026
2027   dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; rad_length = 24.01; // Al
2028   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
2029   InsertLayer(ppl);
2030
2031   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; rad_length = 44.77; // Tedlar
2032   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
2033   InsertLayer(ppl);
2034
2035   dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; rad_length = 44.86; // prepreg
2036   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
2037   InsertLayer(ppl);
2038
2039   dx = 3.; xin = xout; xout = xin + dx; rho = 1.45*0.02; rad_length = 41.28; // Nomex
2040   steps = 10; dx = (xout - xin)/steps;
2041   for(Int_t i=0; i<steps; i++) {
2042     x = xin + i*dx + dx/2;
2043     ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
2044     InsertLayer(ppl);
2045   }
2046
2047   dx = 0.06; xin = xout; xout = xin + dx; rho = 1.45; rad_length = 44.86; // prepreg
2048   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
2049   InsertLayer(ppl);
2050
2051   dx = 50e-4; xin = xout; xout = xin + dx; rho = 1.71; rad_length = 44.77; // Tedlar
2052   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
2053   InsertLayer(ppl);
2054
2055   dx = 50e-4; xin = xout; xout = xin + dx; rho = 2.7; rad_length = 24.01; // Al
2056   ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
2057   InsertLayer(ppl);
2058
2059   Double_t xtrd = (Double_t) fGeom->Rmin();
2060
2061   // add layers between TPC and TRD (Air temporarily)
2062   xin = xout; xout = xtrd;
2063   steps = 50; dx = (xout - xin)/steps;
2064   rho = 1.2e-3;  rad_length = 36.66;
2065
2066   for(Int_t i=0; i<steps; i++) {
2067     x = xin + i*dx + dx/2;
2068     ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
2069     InsertLayer(ppl);
2070   }
2071
2072
2073   Double_t alpha=AliTRDgeometry::GetAlpha();
2074
2075   // add layers for each of the planes
2076
2077   Double_t dxRo = (Double_t) fGeom->CroHght();    // Rohacell
2078   Double_t dxSpace = (Double_t) fGeom->Cspace();  // Spacing between planes
2079   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
2080   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2081   Double_t dxRad = (Double_t) fGeom->CraHght();   // Radiator
2082   Double_t dxTEC = dxRad + dxDrift + dxAmp + dxRo;
2083   Double_t dxPlane = dxTEC + dxSpace;
2084
2085   Int_t tb, tb_index;
2086   const Int_t  nChambers = AliTRDgeometry::Ncham();
2087   Double_t  Ymax = 0, holeYmax = 0;
2088   Double_t *  Zc  = new Double_t[nChambers];
2089   Double_t * Zmax = new Double_t[nChambers];
2090   Double_t  holeZmax = 1000.;   // the whole sector is missing
2091
2092   for(Int_t plane = 0; plane < AliTRDgeometry::Nplan(); plane++) {
2093
2094     // Radiator
2095     xin = xtrd + plane * dxPlane; xout = xin + dxRad;
2096     steps = 12; dx = (xout - xin)/steps; rho = 0.074; rad_length = 40.6;
2097     for(Int_t i=0; i<steps; i++) {
2098       x = xin + i*dx + dx/2;
2099       ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
2100       if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2101         holeYmax = x*TMath::Tan(0.5*alpha);
2102         ppl->SetHole(holeYmax, holeZmax);
2103       }
2104       if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2105         holeYmax = x*TMath::Tan(0.5*alpha);
2106         ppl->SetHole(holeYmax, holeZmax);
2107       }
2108       InsertLayer(ppl);
2109     }
2110
2111     Ymax = fGeom->GetChamberWidth(plane)/2;
2112     for(Int_t ch = 0; ch < nChambers; ch++) {
2113       Zmax[ch] = fGeom->GetChamberLength(plane,ch)/2;
2114       Float_t pad = fPar->GetRowPadSize(plane,ch,0);
2115       Float_t row0 = fPar->GetRow0(plane,ch,0);
2116       Int_t nPads = fPar->GetRowMax(plane,ch,0);
2117       Zc[ch] = (pad * nPads)/2 + row0 - pad/2;
2118     }
2119
2120     dx = fPar->GetTimeBinSize();
2121     rho = 0.00295 * 0.85; rad_length = 11.0;
2122
2123     Double_t x0 = (Double_t) fPar->GetTime0(plane);
2124     Double_t xbottom = x0 - dxDrift;
2125     Double_t xtop = x0 + dxAmp;
2126
2127     // Amplification region
2128
2129     steps = (Int_t) (dxAmp/dx);
2130
2131     for(tb = 0; tb < steps; tb++) {
2132       x = x0 + tb * dx + dx/2;
2133       tb_index = CookTimeBinIndex(plane, -tb-1);
2134       ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,tb_index);
2135       ppl->SetYmax(Ymax);
2136       for(Int_t ch = 0; ch < nChambers; ch++) {
2137         ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
2138       }
2139       if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2140         holeYmax = x*TMath::Tan(0.5*alpha);
2141         ppl->SetHole(holeYmax, holeZmax);
2142       }
2143       if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2144         holeYmax = x*TMath::Tan(0.5*alpha);
2145         ppl->SetHole(holeYmax, holeZmax);
2146       }
2147       InsertLayer(ppl);
2148     }
2149     tb_index = CookTimeBinIndex(plane, -steps);
2150     x = (x + dx/2 + xtop)/2;
2151     dx = 2*(xtop-x);
2152     ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,tb_index);
2153     ppl->SetYmax(Ymax);
2154     for(Int_t ch = 0; ch < nChambers; ch++) {
2155       ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
2156     }
2157     if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2158       holeYmax = x*TMath::Tan(0.5*alpha);
2159       ppl->SetHole(holeYmax, holeZmax);
2160     }
2161     if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2162       holeYmax = x*TMath::Tan(0.5*alpha);
2163       ppl->SetHole(holeYmax, holeZmax);
2164     }
2165     InsertLayer(ppl);
2166
2167     // Drift region
2168     dx = fPar->GetTimeBinSize();
2169     steps = (Int_t) (dxDrift/dx);
2170
2171     for(tb = 0; tb < steps; tb++) {
2172       x = x0 - tb * dx - dx/2;
2173       tb_index = CookTimeBinIndex(plane, tb);
2174
2175       ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,tb_index);
2176       ppl->SetYmax(Ymax);
2177       for(Int_t ch = 0; ch < nChambers; ch++) {
2178         ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
2179       }
2180       if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2181         holeYmax = x*TMath::Tan(0.5*alpha);
2182         ppl->SetHole(holeYmax, holeZmax);
2183       }
2184       if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2185         holeYmax = x*TMath::Tan(0.5*alpha);
2186         ppl->SetHole(holeYmax, holeZmax);
2187       }
2188       InsertLayer(ppl);
2189     }
2190     tb_index = CookTimeBinIndex(plane, steps);
2191     x = (x - dx/2 + xbottom)/2;
2192     dx = 2*(x-xbottom);
2193     ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,tb_index);
2194     ppl->SetYmax(Ymax);
2195     for(Int_t ch = 0; ch < nChambers; ch++) {
2196       ppl->SetZmax(ch, Zc[ch], Zmax[ch]);
2197     }
2198     if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2199       holeYmax = x*TMath::Tan(0.5*alpha);
2200       ppl->SetHole(holeYmax, holeZmax);
2201     }
2202     if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2203       holeYmax = x*TMath::Tan(0.5*alpha);
2204       ppl->SetHole(holeYmax, holeZmax);
2205     }
2206     InsertLayer(ppl);
2207
2208     // Pad Plane
2209     xin = xtop; dx = 0.025; xout = xin + dx; rho = 1.7; rad_length = 33.0;
2210     ppl = new AliTRDpropagationLayer(xin+dx/2,dx,rho,rad_length,-1);
2211     if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2212       holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha);
2213       ppl->SetHole(holeYmax, holeZmax);
2214     }
2215     if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2216       holeYmax = (xin+dx/2)*TMath::Tan(0.5*alpha);
2217       ppl->SetHole(holeYmax, holeZmax);
2218     }
2219     InsertLayer(ppl);
2220
2221     // Rohacell
2222     xin = xout; xout = xtrd + (plane + 1) * dxPlane - dxSpace;
2223     steps = 5; dx = (xout - xin)/steps; rho = 0.074; rad_length = 40.6;
2224     for(Int_t i=0; i<steps; i++) {
2225       x = xin + i*dx + dx/2;
2226       ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
2227       if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2228         holeYmax = x*TMath::Tan(0.5*alpha);
2229         ppl->SetHole(holeYmax, holeZmax);
2230       }
2231       if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2232         holeYmax = x*TMath::Tan(0.5*alpha);
2233         ppl->SetHole(holeYmax, holeZmax);
2234       }
2235       InsertLayer(ppl);
2236     }
2237
2238     // Space between the chambers, air
2239     xin = xout; xout = xtrd + (plane + 1) * dxPlane;
2240     steps = 5; dx = (xout - xin)/steps; rho = 1.29e-3; rad_length = 36.66;
2241     for(Int_t i=0; i<steps; i++) {
2242       x = xin + i*dx + dx/2;
2243       ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
2244       if((fGeom->GetPHOShole()) && (fGeomSector >= 2) && (fGeomSector <= 6)) {
2245         holeYmax = x*TMath::Tan(0.5*alpha);
2246         ppl->SetHole(holeYmax, holeZmax);
2247       }
2248       if((fGeom->GetRICHhole()) && (fGeomSector >= 12) && (fGeomSector <= 14)) {
2249         holeYmax = x*TMath::Tan(0.5*alpha);
2250         ppl->SetHole(holeYmax, holeZmax);
2251       }
2252       InsertLayer(ppl);
2253     }
2254   }
2255
2256   // Space between the TRD and RICH
2257   Double_t xRICH = 500.;
2258   xin = xout; xout = xRICH;
2259   steps = 200; dx = (xout - xin)/steps; rho = 1.29e-3; rad_length = 36.66;
2260   for(Int_t i=0; i<steps; i++) {
2261     x = xin + i*dx + dx/2;
2262     ppl = new AliTRDpropagationLayer(x,dx,rho,rad_length,-1);
2263     InsertLayer(ppl);
2264   }
2265
2266   MapTimeBinLayers();
2267   delete [] Zc;
2268   delete [] Zmax;
2269
2270 }
2271
2272 //______________________________________________________
2273
2274 Int_t  AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t local_tb) const
2275 {
2276   //
2277   // depending on the digitization parameters calculates "global"
2278   // time bin index for timebin <local_tb> in plane <plane>
2279   //
2280
2281   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
2282   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
2283   Double_t dx = (Double_t) fPar->GetTimeBinSize();
2284
2285   Int_t tbAmp = fPar->GetTimeBefore();
2286   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
2287   if(kTRUE) maxAmp = 0;   // intentional until we change parameter class
2288   Int_t tbDrift = fPar->GetTimeMax();
2289   Int_t maxDrift = (Int_t) ((dxDrift+0.000001)/dx);
2290
2291   Int_t tb_per_plane = TMath::Min(tbAmp,maxAmp) + TMath::Min(tbDrift,maxDrift);
2292
2293   Int_t gtb = (plane+1) * tb_per_plane - local_tb - 1 - TMath::Min(tbAmp,maxAmp);
2294
2295   if((local_tb < 0) &&
2296      (TMath::Abs(local_tb) > TMath::Min(tbAmp,maxAmp))) return -1;
2297   if(local_tb >= TMath::Min(tbDrift,maxDrift)) return -1;
2298
2299   return gtb;
2300
2301
2302 }
2303
2304 //______________________________________________________
2305
2306 void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers()
2307 {
2308   //
2309   // For all sensitive time bins sets corresponding layer index
2310   // in the array fTimeBins
2311   //
2312
2313   Int_t index;
2314
2315   for(Int_t i = 0; i < fN; i++) {
2316     index = fLayers[i]->GetTimeBinIndex();
2317
2318     //    printf("gtb %d -> pl %d -> x %f \n", index, i, fLayers[i]->GetX());
2319
2320     if(index < 0) continue;
2321     if(index >= (Int_t) kMAX_TIME_BIN_INDEX) {
2322       printf("*** AliTRDtracker::MapTimeBinLayers: \n");
2323       printf("    index %d exceeds allowed maximum of %d!\n",
2324              index, kMAX_TIME_BIN_INDEX-1);
2325       continue;
2326     }
2327     fTimeBinIndex[index] = i;
2328   }
2329
2330   Double_t x1, dx1, x2, dx2, gap;
2331
2332   for(Int_t i = 0; i < fN-1; i++) {
2333     x1 = fLayers[i]->GetX();
2334     dx1 = fLayers[i]->GetdX();
2335     x2 = fLayers[i+1]->GetX();
2336     dx2 = fLayers[i+1]->GetdX();
2337     gap = (x2 - dx2/2) - (x1 + dx1/2);
2338     if(gap < -0.01) {
2339       printf("*** warning: layers %d and %d are overlayed:\n",i,i+1);
2340       printf("             %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2);
2341     }
2342     if(gap > 0.01) {
2343       printf("*** warning: layers %d and %d have a large gap:\n",i,i+1);
2344       printf("             (%f - %f) - (%f + %f) = %f\n",
2345              x2, dx2/2, x1, dx1, gap);
2346     }
2347   }
2348 }
2349
2350
2351 //______________________________________________________
2352
2353
2354 Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const
2355 {
2356   //
2357   // Returns the number of time bin which in radial position is closest to <x>
2358   //
2359
2360   if(x >= fLayers[fN-1]->GetX()) return fN-1;
2361   if(x <= fLayers[0]->GetX()) return 0;
2362
2363   Int_t b=0, e=fN-1, m=(b+e)/2;
2364   for (; b<e; m=(b+e)/2) {
2365     if (x > fLayers[m]->GetX()) b=m+1;
2366     else e=m;
2367   }
2368   if(TMath::Abs(x - fLayers[m]->GetX()) >
2369      TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
2370   else return m;
2371
2372 }
2373
2374 //______________________________________________________
2375
2376 Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const
2377 {
2378   //
2379   // Returns number of the innermost SENSITIVE propagation layer
2380   //
2381
2382   return GetLayerNumber(0);
2383 }
2384
2385 //______________________________________________________
2386
2387 Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const
2388 {
2389   //
2390   // Returns number of the outermost SENSITIVE time bin
2391   //
2392
2393   return GetLayerNumber(GetNumberOfTimeBins() - 1);
2394 }
2395
2396 //______________________________________________________
2397
2398 Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const
2399 {
2400   //
2401   // Returns number of SENSITIVE time bins
2402   //
2403
2404   Int_t tb, layer;
2405   for(tb = kMAX_TIME_BIN_INDEX-1; tb >=0; tb--) {
2406     layer = GetLayerNumber(tb);
2407     if(layer>=0) break;
2408   }
2409   return tb+1;
2410 }
2411
2412 //______________________________________________________
2413
2414 void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
2415 {
2416   //
2417   // Insert layer <pl> in fLayers array.
2418   // Layers are sorted according to X coordinate.
2419
2420   if ( fN == ((Int_t) kMAX_LAYERS_PER_SECTOR)) {
2421     printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
2422     return;
2423   }
2424   if (fN==0) {fLayers[fN++] = pl; return;}
2425   Int_t i=Find(pl->GetX());
2426
2427   memmove(fLayers+i+1 ,fLayers+i,(fN-i)*sizeof(AliTRDpropagationLayer*));
2428   fLayers[i]=pl; fN++;
2429
2430 }
2431
2432 //______________________________________________________
2433
2434 Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const
2435 {
2436   //
2437   // Returns index of the propagation layer nearest to X
2438   //
2439
2440   if (x <= fLayers[0]->GetX()) return 0;
2441   if (x > fLayers[fN-1]->GetX()) return fN;
2442   Int_t b=0, e=fN-1, m=(b+e)/2;
2443   for (; b<e; m=(b+e)/2) {
2444     if (x > fLayers[m]->GetX()) b=m+1;
2445     else e=m;
2446   }
2447   return m;
2448 }
2449
2450 //______________________________________________________
2451
2452 void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
2453         Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &rad_length,
2454         Bool_t &lookForCluster) const
2455 {
2456   //
2457   // Returns radial step <dx>, density <rho>, rad. length <rad_length>,
2458   // and sensitivity <lookForCluster> in point <y,z>
2459   //
2460
2461   dx  = fdX;
2462   rho = fRho;
2463   rad_length  = fX0;
2464   lookForCluster = kFALSE;
2465
2466   // check dead regions
2467   if(fTimeBinIndex >= 0) {
2468     for(Int_t ch = 0; ch < (Int_t) kZONES; ch++) {
2469       if(TMath::Abs(z - fZc[ch]) < fZmax[ch])
2470         lookForCluster = kTRUE;
2471       //  else { rho = 1.7; rad_length = 33.0; } // G10
2472     }
2473     if(TMath::Abs(y) > fYmax) lookForCluster = kFALSE;
2474     if(!lookForCluster) {
2475       //      rho = 1.7; rad_length = 33.0; // G10
2476     }
2477   }
2478
2479   // check hole
2480   if(fHole && (TMath::Abs(y - fHoleYc) < fHoleYmax) &&
2481               (TMath::Abs(z - fHoleZc) < fHoleZmax)) {
2482     lookForCluster = kFALSE;
2483     rho = fHoleRho;
2484     rad_length  = fHoleX0;
2485   }
2486
2487   return;
2488 }
2489
2490 //______________________________________________________
2491
2492 void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c,
2493                                                           UInt_t index) {
2494
2495 // Insert cluster in cluster array.
2496 // Clusters are sorted according to Y coordinate.
2497
2498   if(fTimeBinIndex < 0) {
2499     printf("*** attempt to insert cluster into non-sensitive time bin!\n");
2500     return;
2501   }
2502
2503   if (fN== (Int_t) kMAX_CLUSTER_PER_TIME_BIN) {
2504     printf("AliTRDpropagationLayer::InsertCluster(): Too many clusters !\n");
2505     return;
2506   }
2507   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
2508   Int_t i=Find(c->GetY());
2509   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
2510   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
2511   fIndex[i]=index; fClusters[i]=c; fN++;
2512 }
2513
2514 //______________________________________________________
2515
2516 Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Double_t y) const {
2517
2518 // Returns index of the cluster nearest in Y
2519
2520   if (y <= fClusters[0]->GetY()) return 0;
2521   if (y > fClusters[fN-1]->GetY()) return fN;
2522   Int_t b=0, e=fN-1, m=(b+e)/2;
2523   for (; b<e; m=(b+e)/2) {
2524     if (y > fClusters[m]->GetY()) b=m+1;
2525     else e=m;
2526   }
2527   return m;
2528 }
2529
2530 //---------------------------------------------------------
2531
2532 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
2533 //
2534 //  Returns correction factor for tilted pads geometry
2535 //
2536
2537   Double_t h01 = sin(TMath::Pi() / 180.0 * fPar->GetTiltingAngle());
2538   Int_t det = c->GetDetector();
2539   Int_t plane = fGeom->GetPlane(det);
2540
2541   if((plane == 1) || (plane == 3) || (plane == 5)) h01=-h01;
2542
2543   if(fNoTilt) h01 = 0;
2544
2545   return h01;
2546 }
2547
2548
2549
2550
2551
2552