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