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