TRD/AliTRDtrackingAnalysis.cxx

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