TRD/AliTRDtrackingAnalysis.cxx

   1
   2 #include "AliTRDtrackingAnalysis.h"
   3
   4
   5 #include "TFile.h"
   6 #include "TTree.h"
   7 #include "TH1D.h"
   8 #include "TH2D.h"
   9 #include "TObjArray.h"
  10 #include "TCanvas.h"
  11 #include "TGeoMatrix.h"
  12 #include "TStyle.h"
  13 #include "TGraphErrors.h"
  14 #include "TF1.h"
  15 #include "TMath.h"
  16
  17 #include "AliRunLoader.h"
  18 #include "AliTRDgeometry.h"
  19 #include "AliRun.h"
  20 #include "AliESD.h"
  21 #include "AliESDtrack.h"
  22 #include "AliTrackReference.h"
  23
  24 #include "AliTRDcluster.h"
  25 #include "AliTRDCommonParam.h"
  26 #include "AliTRDpadPlane.h"
  27 #include "AliTRDcalibDB.h"
  28 #include "AliTracker.h"
  29 #include "AliTRDtracker.h"
  30 //#include "AliTRDtracklet.h"
  31
  32 #include "TGeoManager.h"
  33
  34 //////////////////////////////////////////////////////////////////////////////////////////
  35
  36 AliTRDtrackingAnalysis::AliTRDtrackingAnalysis():
  37   TObject(),
  38   fPath(0),
  39   fRefTPC(0),
  40   fRefTRD(0),
  41   fLoader(0),
  42   fEsdTree(0),
  43   fESD(0),
  44   fTracker(0),
  45   fDeltaPt(0),
  46   fDeltaZ(0),
  47   fDeltaX(0),
  48   fDeltaYPos(0),
  49   fDeltaYNeg(0),
  50   fNPoints(0),
  51   fNGood(0),
  52   fRefSpace(0),
  53   fGeo(0),
  54   fClY2(0),
  55   fClY3(0),
  56   fTgPhi(0),
  57   fGrResTgPhi(0),
  58   fGrMeanTgPhi(0),
  59   fTrklY(0),
  60   fTrklZ(0),
  61   fClZ(0),
  62   fClZZ(0),
  63   fClYY(0),
  64   fClYX(0),
  65   fNLabels(0),
  66   fBits(0),
  67   fRefDx(0),
  68   fClZXref(0),
  69   fClZXcl(0),
  70   fClPos(0)
  71 {
  72
  73   fDeltaX = new TH1D("deltaX", ";delta X (cm)", 100, -1, 1);
  74   fDeltaZ = new TH1D("deltaZ", ";delta Z (cm)", 100, -2, 2);
  75
  76   fDeltaYPos = new TH1D("deltaYpos", ";delta Y (mm)", 100, -1, 1);
  77   fDeltaYNeg = new TH1D("deltaYneg", ";delta Y (mm)", 100, -1, 1);
  78
  79   fNPoints = new TH1D("nPoints", ";np", 40, -0.5, 39.5);
  80   fNGood   = new TH1D("nGood", ";np", 40, -0.5, 39.5);
  81
  82   fDeltaPt = new TH1D("deltaPt", ";delta Pt/Pt (%)", 100, -10, 10);
  83   fRefSpace = new TH2D("refSpace", ";y;x", 120, -60, 60, 200, -4, 1);
  84
  85   fTrklY = new TH1D("trklY", ";delta Y (mm)", 100, -1, 1);
  86   fTrklZ = new TH1D("trklZ", ";delta Z (cm)", 100, -10, 10);
  87
  88
  89   // cluster studies
  90   fClY2 = new TH1D("clY2", ";delta Y (mm)", 100, -10, 10);
  91   fClY3 = new TH1D("clY3", ";delta Y (mm)", 100, -10, 10);
  92
  93   for(int i=0; i<12; i++) // bewere hidden constants in the code
  94     fClYTgPhi[i] = new TH1D(Form("clYtgPhi%d", i), ";delta Y (mm)", 100, -3, 3);
  95
  96   fTgPhi = new TH1D("tgPhi", ";Tg(#phi)", 100, -0.3, 0.3);
  97   fGrResTgPhi = new TGraphErrors();
  98   fGrMeanTgPhi = new TGraphErrors();
  99
 100   //fPullY2 = new TH1D("pullY2", ";pulls Y", 100, -5, 5);
 101   //fPullY3 = new TH1D("pullY3", ";pulls Y", 100, -5, 5);
 102
 103
 104   fClZ = new TH1D("clZ", ";delta Z (cm)", 200, -20, 20);
 105   fClZZ = new TH2D("clZZ", ";z ref;z cl", 600, -300, 300, 600, -300, 300);
 106
 107   fClYY = new TH2D("clYY", ";dY;dY", 100, -3, 3, 100, -3, 3);
 108   fClYX = new TH2D("clYX", ";Y;X", 250, -60, 60, 100, -4, 1);
 109
 110   fNLabels = new TH1D("clLabels", ";n labels", 10, -0.5, 9.5);
 111   fBits = new TH1D("bits", ";bits", 10, -0.5, 9.5);
 112
 113   fRefDx = new TH1D("refDX", ";delta X", 100, 0, 20);
 114   fClPos = new TH2D("clPos", ";z;y", 400, -400, 400, 120, -60, 60);
 115
 116   fClZXref = new TH2D("clZXref", ";z;x", 36, -54, 54, 300, 280, 380);
 117   fClZXcl =  new TH2D("clZXcl", ";z;x", 36, -54, 54, 300, 280, 380);
 118
 119   //fGeo = new AliTRDgeometry();
 120 }
 121
 122 //////////////////////////////////////////////////////////////////////////////////////////
 123
 124 void AliTRDtrackingAnalysis::DrawResolutionPt(int startEvent, int stopEvent) {
 125
 126   CheckFiles();
 127
 128   // loop over ESD events
 129   int nevents = fEsdTree->GetEntries();
 130
 131   for(int iEvent=startEvent; iEvent<nevents && iEvent < stopEvent; iEvent++) {
 132
 133     Info("Draw", "Event = %d", iEvent);
 134
 135     fEsdTree->GetEvent(iEvent);
 136     fLoader->GetEvent(iEvent);
 137     LoadRefs();
 138
 139     int nTracks = fESD->GetNumberOfTracks();
 140     for(int iTrack=0; iTrack<nTracks; iTrack++) {
 141
 142       //Info("Track", "Track = %d", iTrack);
 143       AliESDtrack *esdTrack = fESD->GetTrack(iTrack);
 144       if (!esdTrack->GetInnerParam()) continue;
 145       const AliExternalTrackParam *param = esdTrack->GetOuterParam();
 146       int status = esdTrack->GetStatus();
 147
 148       if (!(status & AliESDtrack::kTRDout)) continue;
 149       if (!(status & AliESDtrack::kTRDrefit)) continue;
 150       if (TMath::Abs(esdTrack->GetOuterParam()->GetPt()) < 1.0) continue;
 151
 152       int ch=0;
 153       while(param->GetX() > fGeo->GetTime0(ch)+2) ch++;
 154       fRefSpace->Fill(2.*ch+0.5, param->GetX() - fGeo->GetTime0(ch));
 155       //if (ch < 5) continue;
 156
 157       double lastX = 0;
 158       int label = abs(esdTrack->GetTRDLabel());
 159       int ntr = 0;
 160       int ngood = 0;
 161
 162       for(int iPoint=GetReference(label); iPoint<fRefTRD->GetEntries(); iPoint++) {
 163
 164         AliTrackReference *aRef = (AliTrackReference*)(*fRefTRD)[iPoint];
 165         if (aRef->GetTrack() != label) break;
 166         ntr++;
 167
 168         lastX = (aRef->LocalX() < lastX)? lastX : aRef->LocalX();
 169         double dx = aRef->LocalX() - param->GetX();
 170         if (TMath::Abs(dx) > 1.) continue;
 171         ngood++;
 172
 173         double bz=fESD->GetMagneticField();
 174         AliExternalTrackParam out(*param);
 175         out.PropagateTo(aRef->LocalX(),bz);
 176
 177         double dp = aRef->Pt() + out.GetPt();
 178         double dy = 10. * (aRef->LocalY() - out.GetY()); // in mm
 179
 180         fDeltaPt->Fill(100. * dp / aRef->Pt());
 181         //fDeltaPt->Fill(out.GetPt() / aRef->Pt());
 182         fDeltaX->Fill(dx);
 183
 184         if (esdTrack->GetSign() > 0) fDeltaYPos->Fill(dy);
 185         else fDeltaYNeg->Fill(dy);
 186
 187         fDeltaZ->Fill(aRef->Z() - out.GetZ());
 188       }
 189
 190       //if (ngood == 0) Info("X", "N = %d, X = %f, DX = %f", ntr, param->GetX(), param->GetX()-lastX);
 191
 192       fNPoints->Fill(ntr);
 193       fNGood->Fill(ngood);
 194     }
 195   }
 196
 197   new TCanvas();
 198   fDeltaPt->Draw();
 199
 200   new TCanvas();
 201   fDeltaX->Draw();
 202
 203   /*
 204   new TCanvas();
 205   fNPoints->Draw();
 206
 207   new TCanvas();
 208   fNGood->Draw();
 209   */
 210
 211   new TCanvas();
 212   fDeltaYPos->Draw();
 213
 214   new TCanvas();
 215   fDeltaYNeg->Draw();
 216
 217   new TCanvas();
 218   fDeltaZ->Draw();
 219
 220   new TCanvas();
 221   fRefSpace->Draw();
 222
 223 }
 224
 225 //////////////////////////////////////////////////////////////////////////////////////////
 226
 227 // void  AliTRDtrackingAnalysis::DrawTrackletResolution(int startEvent, int stopEvent) {
 228
 229 //   LoadRecPointsFile();
 230
 231 //   TFile *file = new TFile(Form("%s/TRD.Tracklets.root", fPath), "read");
 232 //   TTree *tree = (TTree*)file->Get("TRDtracklets");
 233
 234 //   AliTRDtracklet *tracklet = new AliTRDtracklet();
 235 //   tree->SetBranchAddress("tracklets", &tracklet);
 236
 237 //   for(int ev=startEvent; ev<stopEvent; ev++) {
 238
 239 //     gAlice->GetEvent(ev);
 240 //     LoadRefs();
 241
 242 //     int N = tree->GetEntries();
 243 //     for(int i=0; i<N; i++) {
 244
 245 //       tree->GetEntry(i);
 246
 247 //       Double_t yref, zref, tgphi;
 248 //       int stat = GetMCPosition(tracklet->GetLabel(), tracklet->GetX(), yref, zref, tgphi);
 249 //       if (stat < 0) continue;
 250
 251 //       int plane = tracklet->GetPlane();
 252 //       Double_t h01 = tracklet->GetTilt();
 253
 254 //       //printf("Tile = %f\tcorrection = %f um \n", h01, 1e4 * h01 * (tracklet->GetZ()-zref));
 255 //       //double dz = zref - tracklet->GetZ() - cls->GetZ();
 256
 257 //       fTrklY->Fill(10 * (tracklet->GetY() - yref));
 258 //       fTrklZ->Fill(tracklet->GetZ() - zref);
 259
 260 //       int ch=0;
 261 //       while(tracklet->GetX() > fGeo->GetTime0(ch)+2) ch++;
 262 //       fRefSpace->Fill(tracklet->GetY(), tracklet->GetX() - fGeo->GetTime0(ch));
 263 //     }
 264 //   }
 265
 266 //   new TCanvas();
 267 //   fTrklZ->Draw();
 268
 269 //   new TCanvas();
 270 //   fRefSpace->Draw();
 271
 272 //   gStyle->SetOptFit(1);
 273 //   new TCanvas();
 274 //   fTrklY->Draw();
 275 //   fTrklY->Fit("gaus");
 276
 277 //   //new TCanvas();
 278 //   //fClZ->Draw();
 279 // }
 280
 281 //////////////////////////////////////////////////////////////////////////////////////////
 282
 283 void  AliTRDtrackingAnalysis::DrawRecPointResolution(int startEvent, int stopEvent) {
 284
 285   LoadRecPointsFile();
 286   TObjArray *module = new TObjArray();
 287
 288   int nEvents = gAlice->GetEventsPerRun();
 289
 290   for(int ev=startEvent; ev<nEvents && ev < stopEvent; ev++) {
 291
 292     gAlice->GetEvent(ev);
 293     LoadRefs();
 294
 295     TTree *tree = fLoader->GetTreeR("TRD", 0);
 296     tree->SetBranchAddress("TRDcluster", &module);
 297
 298     Info("Res", "Refs Loaded");
 299
 300     int N = tree->GetEntries();
 301     for(int i=0; i<N; i++) {
 302
 303       tree->GetEntry(i);
 304       int m = module->GetEntries();
 305
 306       for(int j=0; j<m; j++) {
 307
 308         AliTRDcluster *cls = (AliTRDcluster*)module->At(j);
 309         if (cls->GetQ() < 10) continue;
 310         fTracker->Transform(cls);
 311         fClPos->Fill(cls->GetZ(), cls->GetY());
 312
 313         int plane = fGeo->GetPlane(cls->GetDetector());
 314
 315         int nl = 0;
 316         for(int k=0; k<3; k++) if (cls->GetLabel(k) > -1) nl++;
 317         fNLabels->Fill(nl);
 318
 319         Double_t yref, zref, tgphi;
 320         int stat = GetMCPosition(cls->GetLabel(0), cls->GetX(), yref, zref, tgphi);
 321         if (stat < 0) continue;
 322
 323         fClZXcl->Fill(cls->GetZ(), cls->GetX());
 324         fClZXref->Fill(zref, cls->GetX());
 325
 326         AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(plane,0);
 327         Double_t h01   = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
 328
 329         //double dz = zref - padPlane->GetRow0();
 330         double dz = zref - cls->GetZ();
 331         double dy = dz * h01;
 332         double yy = cls->GetY() - dy;
 333
 334         if (cls->From2pad()) fClY2->Fill(10 * (yy - yref));
 335         if (cls->From3pad()) fClY3->Fill(10 * (yy - yref));
 336
 337         int idx = GetPhiBin(tgphi);
 338         if (idx >= 0 && idx < 12) fClYTgPhi[idx]->Fill(10 * (yy - yref));
 339
 340         fClZZ->Fill(zref, cls->GetZ());
 341         fClZ->Fill(dz);
 342         fTgPhi->Fill(tgphi);
 343         fClYX->Fill(cls->GetY(), cls->GetX() - fGeo->GetTime0(plane));
 344       }
 345     }
 346   }
 347
 348   new TCanvas();
 349   fClYX->Draw();
 350
 351   //new TCanvas();
 352   //fNLabels->Draw();
 353   new TCanvas();
 354   fClZ->Draw();
 355
 356   new TCanvas();
 357   fTgPhi->Draw();
 358
 359   gStyle->SetOptFit(1);
 360
 361   new TCanvas();
 362   gPad->SetLogy();
 363   fClY2->Draw();
 364   fClY2->Fit("gaus", "", "", -2*fClY2->GetRMS(), 2*fClY2->GetRMS());
 365
 366   new TCanvas();
 367   gPad->SetLogy();
 368   fClY3->Draw();
 369   fClY3->Fit("gaus", "", "", -2*fClY3->GetRMS(), 2*fClY3->GetRMS());
 370
 371   //new TCanvas();
 372   //fRefSpace->Draw();
 373
 374   //new TCanvas();
 375   //fClZXcl->Draw();
 376
 377   //new TCanvas();
 378   //fClZXref->Draw();
 379
 380   /**/
 381   TCanvas *c = new TCanvas();
 382   c->Divide(4,3);
 383
 384   for(int i=0; i<12; i++) {
 385
 386     c->cd(i+1);
 387     fClYTgPhi[i]->Draw();
 388     if (fClYTgPhi[i]->GetSum() < 100) continue;
 389
 390     double mean = fClYTgPhi[i]->GetMean();
 391     double rms = fClYTgPhi[i]->GetRMS();
 392
 393     fClYTgPhi[i]->Fit("gaus", "", "", mean-2*rms, mean+2*rms);
 394     TF1 *f = fClYTgPhi[i]->GetFunction("gaus");
 395
 396     int n = fGrResTgPhi->GetN();
 397     fGrResTgPhi->SetPoint(n, GetPhi(i), f->GetParameter(2));
 398     fGrResTgPhi->SetPointError(n, 0, f->GetParError(2));
 399
 400     fGrMeanTgPhi->SetPoint(n, GetPhi(i), f->GetParameter(1));
 401     fGrMeanTgPhi->SetPointError(n, 0, f->GetParError(1));
 402   }
 403
 404   //gSystem->Sleep(1000);
 405   gStyle->SetOptStat(0);
 406
 407   c = new TCanvas();
 408   TH1D *dummy = new TH1D("dummy", "", 100, -0.3, 0.3);
 409
 410   dummy->SetTitle(";tg(#phi);resolution (mm)");
 411   dummy->SetMinimum(0);
 412   dummy->SetMaximum(1);
 413
 414   //c->cd();
 415   (dummy->Clone("dummy1"))->Draw();
 416
 417   fGrResTgPhi->Draw("PL");
 418   fGrResTgPhi->GetHistogram()->SetTitle(";tg(#phi);resolution (mm)");
 419   fGrResTgPhi->SetMarkerStyle(20);
 420
 421   c = new TCanvas();
 422   dummy->SetTitle(";tg(#phi);mean value (mm)");
 423   dummy->SetMinimum(-0.3);
 424   dummy->SetMaximum(0.3);
 425
 426   //c->cd();
 427   (dummy->Clone("dummy2"))->Draw();
 428
 429   fGrMeanTgPhi->Draw("PL");
 430   fGrMeanTgPhi->GetHistogram()->SetTitle(";tg(#phi);mean value (mm)");
 431   fGrMeanTgPhi->SetMarkerStyle(20);
 432   /**/
 433
 434
 435   //new TCanvas();
 436   //fClZZ->Draw("colz");
 437
 438   //new TCanvas();
 439   //fClPos->Draw("colz");
 440
 441   //new TCanvas();
 442   //fBits->Draw();
 443
 444   //new TCanvas();
 445   //fRefDx->Draw();
 446
 447 }
 448
 449 //////////////////////////////////////////////////////////////////////////////////////////
 450
 451 void AliTRDtrackingAnalysis::LoadRecPointsFile() {
 452
 453   char filename[256];
 454   sprintf(filename, "%s/galice.root", fPath);
 455
 456   fLoader = AliRunLoader::Open(filename);
 457   if (!fLoader) {
 458     Error("CheckFiles", "getting run loader from file %s/galice.root failed", filename);
 459     return;
 460   }
 461
 462   fLoader->LoadgAlice();
 463   gAlice = fLoader->GetAliRun();
 464
 465   if (!gAlice) {
 466     Error("CheckFiles", "no galice object found");
 467     return;
 468   }
 469
 470   fLoader->LoadKinematics();
 471   fLoader->LoadHeader();
 472   fLoader->LoadTrackRefs();
 473
 474   //TGeoManager::Import("/data/alice_u/radomski/condor/run_0/geometry.root");
 475   TGeoManager::Import(Form("%s/geometry.root", fPath));
 476
 477
 478   fLoader->CdGAFile();
 479   fGeo = (AliTRDgeometry*)gDirectory->Get("TRDgeometry");
 480   fTracker = new AliTRDtracker(gFile);
 481
 482   fLoader->LoadRecPoints("TRD");
 483
 484   AliTracker::SetFieldMap(gAlice->Field(), 1);
 485 }
 486
 487 //////////////////////////////////////////////////////////////////////////////////////////
 488
 489 void  AliTRDtrackingAnalysis::CheckFiles() {
 490
 491   // MC info
 492
 493   char filename[256];
 494   sprintf(filename, "%s/galice.root", fPath);
 495
 496   fLoader = AliRunLoader::Open(filename);
 497   if (!fLoader) {
 498     Error("CheckFiles", "getting run loader from file %s/galice.root failed", filename);
 499     return;
 500   }
 501
 502   fLoader->LoadgAlice();
 503   gAlice = fLoader->GetAliRun();
 504
 505   if (!gAlice) {
 506     Error("CheckFiles", "no galice object found");
 507     return;
 508   }
 509
 510   fLoader->LoadKinematics();
 511   fLoader->LoadHeader();
 512   fLoader->LoadTrackRefs();
 513
 514   fLoader->CdGAFile();
 515   fGeo = (AliTRDgeometry*)gDirectory->Get("TRDgeometry");
 516   //fGeo->ReadGeoMatrices();
 517
 518   // ESD
 519
 520   sprintf(filename,"%s/AliESDs.root", fPath);
 521   TFile *esdFile = new TFile(filename, "READ");
 522
 523   if (esdFile->IsZombie()) {
 524     Error("CheckFiles", "file not present: AliESDs.root");
 525     return;
 526   }
 527
 528   fEsdTree = (TTree*)esdFile->Get("esdTree");
 529   fESD = new AliESD();
 530
 531   fEsdTree->SetBranchAddress("ESD", &fESD);
 532 }
 533
 534 //////////////////////////////////////////////////////////////////////////////////////////
 535
 536 void  AliTRDtrackingAnalysis::LoadRefs() {
 537
 538   if (fRefTPC) delete fRefTPC;
 539   if (fRefTRD) delete fRefTRD;
 540
 541   fRefTPC = new TObjArray();
 542   fRefTRD = new TObjArray();
 543
 544   //fLoader->GetEvent(event);
 545   //AliStack* stack = gAlice->Stack();
 546   TTree *refTree = fLoader->TreeTR();
 547
 548   const int nBranch = 2;
 549   const char *brName[] = {"TPC", "TRD"};
 550   TClonesArray *clRefs = new TClonesArray("AliTrackReference");
 551
 552   for(int b=0; b<nBranch; b++) {
 553
 554     TBranch *branch = refTree->GetBranch(brName[b]);
 555     refTree->SetBranchAddress(brName[b],&clRefs);
 556
 557     int nEntries = branch->GetEntries();
 558     for(int iTrack = 0; iTrack < nEntries; iTrack++) {
 559
 560       refTree->GetEvent(iTrack);
 561       int nPoints =  clRefs->GetEntries();
 562       for(int iPoint=0; iPoint<nPoints; iPoint++) {
 563         AliTrackReference *ref = (AliTrackReference*)clRefs->At(iPoint);
 564         if (b == 0) fRefTPC->Add(new AliTrackReference(*ref));
 565         if (b == 1) fRefTRD->Add(new AliTrackReference(*ref));
 566       }
 567     }
 568   }
 569
 570   fRefTPC->Sort();
 571   fRefTRD->Sort();
 572
 573   for(int i=0; i<fRefTRD->GetEntries(); i++) {
 574     AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
 575     fLabels[i] = ref->GetTrack();
 576
 577     int p=0;
 578     while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
 579     fRefSpace->Fill(ref->LocalY(), ref->LocalX()-fGeo->GetTime0(p));
 580
 581     //for(int bit=0; bit<9; bit++) if (ref->TestBit(bit)) fBits->Fill(bit);
 582   }
 583
 584   delete clRefs;
 585   Info("LoadRefs", "TPC = %d\t TRD = %d", fRefTPC->GetEntries(), fRefTRD->GetEntries());
 586 }
 587
 588 //////////////////////////////////////////////////////////////////////////////////////////
 589
 590 Int_t AliTRDtrackingAnalysis::GetReference(Int_t label) {
 591
 592   int start = TMath::BinarySearch(fRefTRD->GetEntries(), fLabels, label);
 593
 594   while (start >= 0) {
 595     AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[start];
 596     if (ref->GetTrack() != label) return start+1;
 597     start--;
 598   }
 599
 600   return 0;
 601 }
 602
 603 //////////////////////////////////////////////////////////////////////////////////////////
 604
 605 int AliTRDtrackingAnalysis::GetMCPosition(Int_t label, Double_t x, Double_t &Y, Double_t &Z, Double_t &tgphi) {
 606
 607   double lowX = 100.;
 608   double highX = 100.;
 609   int idLow = -1;
 610   int idHigh = -1;
 611
 612   int nref= 0;
 613   int idx = GetReference(label);
 614   for(int i=idx; i<fRefTRD->GetEntries(); i++) {
 615
 616     AliTrackReference *ref = (AliTrackReference*)(*fRefTRD)[i];
 617     if (ref->GetTrack() != label) break;
 618     nref++;
 619
 620     //int p=0;
 621     //while(ref->LocalX() > fGeo->GetTime0(p)+2) p++;
 622     //if (p != layer) continue;
 623
 624     double dX = ref->LocalX()-x;
 625     if ( dX > 0 ) {
 626       if (dX < highX) {
 627         idHigh = i;
 628         highX = dX;
 629       }
 630     } else {
 631       dX = TMath::Abs(dX);
 632       if (dX < lowX) {
 633         idLow = i;
 634         lowX = dX;
 635       }
 636     }
 637   }
 638
 639   if (idLow == -1 || idHigh == -1) return -1;
 640
 641   AliTrackReference *refI = (AliTrackReference*)(*fRefTRD)[idLow];
 642   AliTrackReference *refO = (AliTrackReference*)(*fRefTRD)[idHigh];
 643
 644
 645   double dx = refO->LocalX() - refI->LocalX();
 646   double dy = refO->LocalY() - refI->LocalY();
 647   double dz = refO->Z() - refI->Z();
 648   double ddx = (x - refI->LocalX())/dx;
 649
 650   fRefDx->Fill(dx);
 651
 652   Y = refI->LocalY() + ddx * dy;
 653   Z = refI->Z() + ddx * dz;
 654
 655   tgphi = dy/dx;
 656
 657   return 0;
 658 }
 659
 660 //////////////////////////////////////////////////////////////////////////////////////////
 661 Int_t AliTRDtrackingAnalysis::GetPhiBin(Double_t phi)  {
 662   return (int)((phi+0.3)/0.05);
 663 }
 664
 665 //////////////////////////////////////////////////////////////////////////////////////////
 666
 667 Double_t AliTRDtrackingAnalysis::GetPhi(Int_t bin) {
 668   return bin * 0.05 - 0.3 + 0.025;
 669 }
 670 //////////////////////////////////////////////////////////////////////////////////////////