ANALYSIS/AliTrackPoints.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 #include "AliTrackPoints.h"
  19 #include "AliLog.h"
  20
  21 //_________________________________
  22 ////////////////////////////////////////////////////////////
  23 //                                                        //
  24 // class AliTrackPoints                                //
  25 //                                                        //
  26 // used by Anti-Merging cut                               //
  27 // contains set of poits the lay on track trajectory      //
  28 // according to reconstructed track parameters -          //
  29 // NOT CLUSTERS POSITIONS!!!                              //
  30 // Anti-Merging cut is applied only on tracks coming from //
  31 // different events (that are use to fill deniminators)   //
  32 //                                                        //
  33 ////////////////////////////////////////////////////////////
  34
  35 #include <TClonesArray.h>
  36 #include <TFile.h>
  37 #include <TMath.h>
  38
  39 #include "AliESDtrack.h"
  40 #include "AliTPCtrack.h"
  41 #include "AliTrackReference.h"
  42 #include "AliITStrackV2.h"
  43
  44 ClassImp(AliTrackPoints)
  45
  46 Int_t AliTrackPoints::fgDebug = 0;
  47
  48 AliTrackPoints::AliTrackPoints():
  49  fN(0),
  50  fX(0x0),
  51  fY(0x0),
  52  fZ(0x0)
  53 {
  54   //constructor
  55 }
  56 /***************************************************************/
  57
  58 AliTrackPoints::AliTrackPoints(AliTrackPoints::ETypes type, AliESDtrack* track, Float_t mf):
  59  fN(0),
  60  fX(0x0),
  61  fY(0x0),
  62  fZ(0x0)
  63 {
  64   //constructor
  65   //tupe -  what kind of track points should be calculated
  66   //mf - magnetic field in [kG] = [T]*10.0
  67   switch (type)
  68    {
  69      case kITS:
  70        fN = 6;
  71        fX = new Float_t[fN];
  72        fY = new Float_t[fN];
  73        fZ = new Float_t[fN];
  74        MakeITSPoints(track);
  75        break;
  76
  77      case kITSInnerFromVertexOuterFromTPC:
  78        fN = 6;
  79        fX = new Float_t[fN];
  80        fY = new Float_t[fN];
  81        fZ = new Float_t[fN];
  82        MakeITSPointsInnerFromVertexOuterFromTPC(track,mf);
  83        break;
  84
  85      default:
  86        Info("AliTrackPoints","Not recognized type");
  87    }
  88
  89 }
  90 /***************************************************************/
  91
  92 AliTrackPoints::AliTrackPoints(Int_t n, AliESDtrack* track, Float_t mf, Float_t dr, Float_t r0):
  93  fN(n),
  94  fX(new Float_t[fN]),
  95  fY(new Float_t[fN]),
  96  fZ(new Float_t[fN])
  97 {
  98   //constructor
  99   //mf - magnetic field in kG - needed to calculate curvature out of Pt
 100   //r0 - starting radius
 101   //dr - calculate points every dr cm, default every 30cm
 102   if (track == 0x0)
 103    {
 104      Error("AliTrackPoints","ESD track is null");
 105      fN = 0;
 106      delete [] fX;
 107      delete [] fY;
 108      delete [] fZ;
 109      fX = fY = fZ = 0x0;
 110      return;
 111    }
 112
 113   if ( ((track->GetStatus() & AliESDtrack::kTPCrefit) == kFALSE)&&
 114        ((track->GetStatus() & AliESDtrack::kTPCin) == kFALSE)  )
 115    {
 116      //could happend: its stand alone tracking
 117      AliDebug(3,"This ESD track does not contain TPC information");
 118
 119      fN = 0;
 120      delete [] fX;
 121      delete [] fY;
 122      delete [] fZ;
 123      fX = fY = fZ = 0x0;
 124
 125      return;
 126    }
 127
 128   Double_t x;
 129   Double_t par[5];
 130   track->GetInnerExternalParameters(x,par);     //get properties of the track
 131   if (par[4] == 0)
 132    {
 133      Error("AliTrackPoints","This ESD track seem not to contain TPC information (curv is 0)");
 134      return;
 135    }
 136
 137   if (mf == 0.0)
 138    {
 139      Error("AliTrackPoints","Zero Magnetic field passed as parameter.");
 140      return;
 141    }
 142
 143   Double_t alpha = track->GetInnerAlpha();
 144   Double_t cc = 1000./0.299792458/mf;//conversion constant
 145   Double_t c=par[4]/cc;
 146
 147   MakePoints(dr,r0,x,par,c,alpha);
 148
 149 }
 150 /***************************************************************/
 151
 152 AliTrackPoints::AliTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr, Float_t r0):
 153  fN(n),
 154  fX(new Float_t[fN]),
 155  fY(new Float_t[fN]),
 156  fZ(new Float_t[fN])
 157 {
 158   //constructor
 159   //r0 starting radius
 160   //dr - calculate points every dr cm, default every 30cm
 161   if (track == 0x0)
 162    {
 163      Error("AliTrackPoints","TPC track is null");
 164      fN = 0;
 165      delete [] fX;
 166      delete [] fY;
 167      delete [] fZ;
 168      fX = fY = fZ = 0x0;
 169      return;
 170    }
 171   track->PropagateTo(r0);
 172
 173  //* This formation is now fixed in the following way:                         *
 174  //*      external param0:   local Y-coordinate of a track (cm)                *
 175  //*      external param1:   local Z-coordinate of a track (cm)                *
 176  //*      external param2:   local sine of the track momentum azimuth angle    *
 177  //*      external param3:   tangent of the track momentum dip angle           *
 178  //*      external param4:   1/pt (1/(GeV/c))                                  *
 179
 180   Double_t x = 0;
 181   Double_t par[5];
 182   track->GetExternalParameters(x,par);     //get properties of the track
 183
 184   Double_t alpha = track->GetAlpha();
 185   Double_t c=track->GetC();
 186   MakePoints(dr,r0,x,par,c,alpha);
 187 }
 188 /***************************************************************/
 189
 190 AliTrackPoints::~AliTrackPoints()
 191 {
 192   //destructor
 193   delete [] fX;
 194   delete [] fY;
 195   delete [] fZ;
 196 }
 197 /***************************************************************/
 198
 199 void AliTrackPoints::MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t* par, Double_t c, Double_t alpha)
 200 {
 201   //Calculates points starting at radius r0
 202   //spacing every dr (in radial direction)
 203   // according to track parameters
 204   // x - position in sector local reference frame. x i parallel to R and sector is symmetric with respect to x
 205   // par - track external parameters; array with 5 elements;  look at AliTPCtrack.h or AliESDtrack.h for their meaning
 206   // c - track curvature
 207   // alpha - sector's rotation angle (phi) == angle needed for local to global transformation
 208
 209   Double_t y = par[0];
 210   Double_t z0 = par[1];
 211
 212   Double_t phi0local = TMath::ATan2(y,x);
 213   Double_t phi0global = phi0local + alpha;
 214
 215   if (phi0local<0) phi0local+=2*TMath::Pi();
 216   if (phi0local>=2.*TMath::Pi()) phi0local-=2*TMath::Pi();
 217
 218   if (phi0global<0) phi0global+=2*TMath::Pi();
 219   if (phi0global>=2.*TMath::Pi()) phi0global-=2*TMath::Pi();
 220
 221   Double_t r = TMath::Hypot(x,y);
 222
 223
 224   AliDebug(9,Form("Radius0 %f, Real Radius %f",r0,r));
 225
 226   AliDebug(5,Form("Phi Global at first padraw %f, Phi locat %f",phi0global,phi0local));
 227
 228   Double_t eta = x*c - par[2] ;//par[2] = fX*C - eta; eta==fP2 ; C==fP4
 229
 230   //this calculattions are assuming new (current) model
 231   Double_t tmp = par[2];
 232   tmp = 1. - tmp*tmp;
 233   tmp = c*y + TMath::Sqrt(tmp);
 234   Double_t dca=(TMath::Hypot(eta,tmp) - 1. )/TMath::Abs(c);
 235
 236   //Here is old model Cold=Cnew/2.
 237   Double_t dcasq = dca*dca;
 238   Double_t c2 = c/2.;
 239   Double_t cst1 = (1.+c2*dca)*dca;//first constant
 240   Double_t cst2 = 1. + 2.*c2*dca;//second constant
 241
 242   Double_t factorPhi0 = TMath::ASin((c2*r + cst1/r)/cst2);
 243   Double_t factorZ0 = TMath::ASin(c2*TMath::Sqrt((r*r-dcasq)/cst2))*par[3]/c2;
 244
 245   for(Int_t i = 0; i<fN; i++)
 246    {
 247      Double_t rc = r0 + i*dr;
 248      Double_t ftmp = (c2*rc + cst1/rc)/cst2;
 249      if (ftmp > 1.0)
 250       {
 251         AliDebug(1,Form("ASin argument > 1 %f:",ftmp));
 252         ftmp=1.0;
 253       }
 254      else if (ftmp < -1.0)
 255       {
 256         AliDebug(1,Form("ASin argument < -1 %f:",ftmp));
 257         ftmp=-1.0;
 258       }
 259
 260      Double_t factorPhi = TMath::ASin( ftmp );//factor phi od rc
 261      Double_t phi = phi0global + factorPhi - factorPhi0;
 262
 263      ftmp = (rc*rc-dcasq)/cst2;
 264      if (ftmp < 0.0)
 265       {
 266         AliDebug(1,Form("Sqrt argument < 0: %f",ftmp));
 267         ftmp=0.0;
 268       }
 269
 270      ftmp = c2*TMath::Sqrt(ftmp);
 271      if (ftmp > 1.0)
 272       {
 273         AliDebug(1,Form("ASin argument > 1: %f",ftmp));
 274         ftmp=1.0;
 275       }
 276      else if (ftmp < -1.0)
 277       {
 278         AliDebug(2,Form("ASin argument < -1: %f",ftmp));
 279         ftmp=-1.0;
 280       }
 281      Double_t factorZ = TMath::ASin(ftmp)*par[3]/c2;
 282      fZ[i] = z0 + factorZ - factorZ0;
 283      fX[i] = rc*TMath::Cos(phi);
 284      fY[i] = rc*TMath::Sin(phi);
 285
 286      AliDebug(3,Form("AliTrackPoints","X %f Y %f Z %f R asked %f R obtained %f",
 287              fX[i],fY[i],fZ[i],rc,TMath::Hypot(fX[i],fY[i])));
 288    }
 289 }
 290 /***************************************************************/
 291
 292 void AliTrackPoints::MakeITSPoints(AliESDtrack* track)
 293 {
 294 //Calculates points in ITS
 295 // z=R*Pz/Pt
 296  AliITStrackV2 itstrack(*track,kTRUE);
 297  Double_t x,y,z;
 298  static const Double_t r[6] = {4.0, 7.0, 14.9, 23.8, 39.1, 43.6};
 299  for (Int_t i = 0; i < 6; i++)
 300   {
 301     itstrack.GetGlobalXYZat(r[i],x,y,z);
 302     fX[i] = x;
 303     fY[i] = y;
 304     fZ[i] = z;
 305 //    Info("MakeITSPoints","X %f Y %f Z %f R asked %f R obtained %f",
 306 //             fX[i],fY[i],fZ[i],r[i],TMath::Hypot(fX[i],fY[i]));
 307   }
 308
 309 }
 310
 311 /***************************************************************/
 312 void AliTrackPoints::MakeITSPointsInnerFromVertexOuterFromTPC(AliESDtrack* track, Float_t mf)
 313 {
 314 //makes trackpoints for ITS
 315 //for 3 inner layers calculates out of the vector at vertex
 316 //for 3 outer ---------------//------------------ at inner TPC
 317
 318  static const Double_t r[6] = {4.0, 7.0, 14.9, 23.8, 39.1, 43.6};
 319  AliITStrackV2 itstrack(*track,kTRUE);
 320  Double_t x,y,z;
 321  for (Int_t i = 0; i < 3; i++)
 322   {
 323     itstrack.GetGlobalXYZat(r[i],x,y,z);
 324     fX[i] = x;
 325     fY[i] = y;
 326     fZ[i] = z;
 327     AliDebug(3,Form("X %f Y %f Z %f R asked %f R obtained %f",
 328                 fX[i],fY[i],fZ[i],r[i],TMath::Hypot(fX[i],fY[i])));
 329   }
 330
 331  for (Int_t i = 3; i < 6; i++)
 332   {
 333     Float_t ax,ay,az;
 334     AliTrackPoints tmptp(1,track,mf,0,r[i]);
 335     tmptp.PositionAt(0,ax,ay,az);
 336     fX[i] = ax;
 337     fY[i] = ay;
 338     fZ[i] = az;
 339     AliDebug(3,Form("X %f Y %f Z %f R asked %f R obtained %f",
 340                 fX[i],fY[i],fZ[i],r[i],TMath::Hypot(fX[i],fY[i])));
 341   }
 342
 343 }
 344
 345
 346 /***************************************************************/
 347
 348 void AliTrackPoints::PositionAt(Int_t n, Float_t &x,Float_t &y,Float_t &z)
 349 {
 350   //returns position at point n
 351   if ((n<0) || (n>=fN))
 352    {
 353      Error("PositionAt","Point %d out of range",n);
 354      return;
 355    }
 356
 357   x = fX[n];
 358   y = fY[n];
 359   z = fZ[n];
 360   AliDebug(2,Form("n %d; X %f; Y %f; Z %f",n,x,y,z));
 361
 362 }
 363 /***************************************************************/
 364
 365 void AliTrackPoints::Move(Float_t x, Float_t y, Float_t z)
 366 {
 367 //Moves all points about vector
 368  for (Int_t i = 0; i<fN; i++)
 369    {
 370      fX[i]+=x;
 371      fY[i]+=y;
 372      fZ[i]+=z;
 373    }
 374 }
 375 /***************************************************************/
 376
 377 Double_t AliTrackPoints::AvarageDistance(const AliTrackPoints& tr)
 378 {
 379   //returns the aritmethic avarage distance between two tracks
 380 //  Info("AvarageDistance","Entered");
 381   if ( (fN <= 0) || (tr.fN <=0) )
 382    {
 383      AliDebug(1,"One of tracks is empty");
 384      return -1;
 385    }
 386
 387   if (fN != tr.fN)
 388    {
 389      Warning("AvarageDistance","Number of points is not equal");
 390      return -1;
 391    }
 392
 393   Double_t sum = 0;
 394   for (Int_t i = 0; i<fN; i++)
 395    {
 396      AliDebug(10,Form("radii: %f %f",TMath::Hypot(fX[i],fY[i]),TMath::Hypot(tr.fX[i],tr.fY[i])));
 397 //       Float_t r1sq = fX[i]*fX[i]+fY[i]*fY[i];
 398 //       Float_t r2sq = tr.fX[i]*tr.fX[i]+tr.fY[i]*tr.fY[i];
 399
 400
 401      Double_t dx = fX[i]-tr.fX[i];
 402      Double_t dy = fY[i]-tr.fY[i];
 403      Double_t dz = fZ[i]-tr.fZ[i];
 404      sum+=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 405
 406      AliDebug(2,Form("Diff: x ,y z: %f , %f, %f",dx,dy,dz));
 407      AliDebug(2,Form("xxyyzz %f %f %f %f %f %f",
 408             fX[i],tr.fX[i],fY[i],tr.fY[i],fZ[i],tr.fZ[i]));
 409    }
 410
 411   Double_t retval = sum/((Double_t)fN);
 412   AliDebug(1,Form("Avarage distance is %f.",retval));
 413
 414   return retval;
 415 }
 416 /***************************************************************/
 417
 418 void AliTrackPoints::Print(Option_t* /*option*/) const
 419 {
 420
 421   Info("Print","There is %d points",fN);
 422   for(Int_t i = 0; i < fN; i++)
 423    {
 424      Info("Print","%d: %f %f %f",i,fX[i],fY[i],fZ[i]);
 425    }
 426
 427 }
 428
 429 /***************************************************************/
 430 /***************************************************************/
 431 /***************************************************************/
 432 /***************************************************************/
 433 /***************************************************************/
 434 /***************************************************************/
 435
 436 #include "AliRun.h"
 437 #include "AliESD.h"
 438 #include "AliRunLoader.h"
 439 #include "AliTPCtrack.h"
 440 #include "TTree.h"
 441 #include "TBranch.h"
 442 #include "TH2D.h"
 443 #include "TCanvas.h"
 444 #include "AliMagF.h"
 445
 446
 447
 448
 449 void AliTrackPoints::Testesd(Int_t entr,const char* fname )
 450 {
 451   delete gAlice;
 452   gAlice = 0x0;
 453   AliRunLoader* rl = AliRunLoader::Open();
 454   rl->LoadgAlice();
 455
 456   Float_t mf = rl->GetAliRun()->Field()->SolenoidField();
 457
 458
 459   TFile* fFile = TFile::Open(fname);
 460
 461   if (fFile == 0x0)
 462    {
 463      printf("testesd: There is no suche a ESD file\n");
 464      return;
 465    }
 466   AliESD* esd = dynamic_cast<AliESD*>(fFile->Get("0"));
 467   AliESDtrack *t = esd->GetTrack(entr);
 468   if (t == 0x0)
 469    {
 470      ::Error("testesd","Can not get track %d",entr);
 471      return;
 472    }
 473
 474
 475   Int_t N = 170;
 476   AliTrackPoints* tp = new AliTrackPoints(N,t,mf,1.);
 477
 478   Float_t xmin = -250;
 479   Float_t xmax =  250;
 480
 481   Float_t ymin = -250;
 482   Float_t ymax =  250;
 483
 484   Float_t zmin = -250;
 485   Float_t zmax =  250;
 486
 487   TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
 488   TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
 489   TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);
 490
 491   TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
 492   TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
 493   TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);
 494
 495   hxyt->SetDirectory(0x0);
 496   hxy->SetDirectory(0x0);
 497   hxyTR->SetDirectory(0x0);
 498
 499   hxzt->SetDirectory(0x0);
 500   hxz->SetDirectory(0x0);
 501   hxzTR->SetDirectory(0x0);
 502
 503   Float_t x,y,z;
 504
 505   for (Int_t i = 0;i<N;i++)
 506    {
 507       Double_t r = 84.1+i;
 508       tp->PositionAt(i,x,y,z);
 509       hxy->Fill(x,y);
 510       hxz->Fill(x,z);
 511       printf("Rdemanded %f\n",r);
 512       printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
 513
 514    }
 515
 516   rl->LoadTrackRefs();
 517   TTree* treeTR = rl->TreeTR();
 518   TBranch* b = treeTR->GetBranch("TPC");
 519
 520   TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
 521   AliTrackReference* tref;
 522   b->SetAddress(&trackrefs);
 523
 524   Int_t tlab = TMath::Abs(t->GetLabel());
 525
 526   Int_t netr = (Int_t)treeTR->GetEntries();
 527   printf("Found %d entries in TR tree\n",netr);
 528
 529   for (Int_t e = 0; e < netr; e++)
 530    {
 531      treeTR->GetEntry(e);
 532      tref = (AliTrackReference*)trackrefs->At(0);
 533      if (tref == 0x0) continue;
 534      if (tref->GetTrack() != tlab) continue;
 535
 536      printf("Found %d entries in TR array\n",trackrefs->GetEntries());
 537
 538      for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
 539       {
 540         tref = (AliTrackReference*)trackrefs->At(i);
 541         if (tref->GetTrack() != tlab) continue;
 542         x = tref->X();
 543         y = tref->Y();
 544         z = tref->Z();
 545         printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
 546
 547         hxzTR->Fill(x,z);
 548         hxyTR->Fill(x,y);
 549         for (Int_t j = 1; j < 10; j++)
 550          {
 551             hxyTR->Fill(x, y+j*0.1);
 552             hxyTR->Fill(x, y-j*0.1);
 553             hxyTR->Fill(x+j*0.1,y);
 554             hxyTR->Fill(x-j*0.1,y);
 555
 556             hxzTR->Fill(x,z-j*0.1);
 557             hxzTR->Fill(x,z+j*0.1);
 558             hxzTR->Fill(x-j*0.1,z);
 559             hxzTR->Fill(x+j*0.1,z);
 560          }
 561       }
 562       break;
 563     }
 564   hxy->Draw("");
 565 //  hxzt->Draw("same");
 566   hxyTR->Draw("same");
 567
 568   delete rl;
 569 }
 570
 571 /***************************************************************/
 572 /***************************************************************/
 573 /***************************************************************/
 574
 575 void AliTrackPoints::Testtpc(Int_t entr)
 576 {
 577   delete gAlice;
 578   gAlice = 0x0;
 579   AliRunLoader* rl = AliRunLoader::Open();
 580   AliLoader* l = rl->GetLoader("TPCLoader");
 581   rl->LoadgAlice();
 582   AliKalmanTrack::SetConvConst(100/0.299792458/0.2/rl->GetAliRun()->Field()->Factor());
 583   l->LoadTracks();
 584   AliTPCtrack* t = new AliTPCtrack();
 585   TBranch* b=l->TreeT()->GetBranch("tracks");
 586   b->SetAddress(&t);
 587   l->TreeT()->GetEntry(entr);
 588   Int_t N = 160;
 589   AliTrackPoints* tp = new AliTrackPoints(N,t,1.);
 590
 591   Float_t xmin = -250;
 592   Float_t xmax =  250;
 593
 594   Float_t ymin = -250;
 595   Float_t ymax =  250;
 596
 597   Float_t zmin = -250;
 598   Float_t zmax =  250;
 599
 600   TH2D* hxy = new TH2D("hxy","hxy",1000,xmin,xmax,1000,ymin,ymax);
 601   TH2D* hxyt = new TH2D("hxyt","hxyt",1000,xmin,xmax,1000,ymin,ymax);
 602   TH2D* hxyTR = new TH2D("hxyTR","hxyTR",1000,xmin,xmax,1000,ymin,ymax);
 603
 604   TH2D* hxz = new TH2D("hxz","hxz",1000,xmin,xmax,1000,zmin,zmax);
 605   TH2D* hxzt = new TH2D("hxzt","hxzt",1000,xmin,xmax,1000,zmin,zmax);
 606   TH2D* hxzTR = new TH2D("hxzTR","hxzTR",1000,xmin,xmax,1000,zmin,zmax);
 607
 608   hxyt->SetDirectory(0x0);
 609   hxy->SetDirectory(0x0);
 610   hxyTR->SetDirectory(0x0);
 611
 612   hxzt->SetDirectory(0x0);
 613   hxz->SetDirectory(0x0);
 614   hxzTR->SetDirectory(0x0);
 615
 616   Float_t x,y,z;
 617
 618   for (Int_t i = 0;i<N;i++)
 619    {
 620       Double_t r = 84.1+i;
 621       tp->PositionAt(i,x,y,z);
 622       hxy->Fill(x,y);
 623       hxz->Fill(x,z);
 624       printf("Rdemanded %f\n",r);
 625       printf("tpx %f tpy %f tpz %f Rt =%f\n", x,y,z,TMath::Hypot(x,y));
 626
 627       //BUT they are local!!!!
 628       t->PropagateTo(r);
 629 //      Double_t phi = t->Phi();
 630       Double_t rl = TMath::Hypot(t->GetX(),t->GetY());//real radius
 631
 632       Double_t alpha = t->GetAlpha();
 633       Double_t salpha = TMath::Sin(alpha);
 634       Double_t calpha = TMath::Cos(alpha);
 635       x = t->GetX()*calpha - t->GetY()*salpha;
 636       y = t->GetX()*salpha + t->GetY()*calpha;
 637       z = t->GetZ();
 638
 639       printf("tx %f ty %f tz %f Rt = %f R from XY %f\n",x,y,z,TMath::Hypot(x,y),rl);
 640
 641       printf("tpz - tz %f\n",z-t->GetZ());
 642       printf("\n");
 643       hxyt->Fill(x,y);
 644       hxzt->Fill(x,z);
 645
 646    }
 647
 648   rl->LoadTrackRefs();
 649   TTree* treeTR = rl->TreeTR();
 650   b = treeTR->GetBranch("TPC");
 651
 652   TClonesArray* trackrefs = new TClonesArray("AliTrackReference", 100);
 653   AliTrackReference* tref;
 654   b->SetAddress(&trackrefs);
 655
 656   Int_t tlab = TMath::Abs(t->GetLabel());
 657
 658   Int_t netr = (Int_t)treeTR->GetEntries();
 659   printf("Found %d entries in TR tree\n",netr);
 660
 661   for (Int_t e = 0; e < netr; e++)
 662    {
 663      treeTR->GetEntry(e);
 664      tref = (AliTrackReference*)trackrefs->At(0);
 665      if (tref == 0x0) continue;
 666      if (tref->GetTrack() != tlab) continue;
 667
 668      printf("Found %d entries in TR array\n",trackrefs->GetEntries());
 669
 670      for (Int_t i = 0; i < trackrefs->GetEntries(); i++)
 671       {
 672         tref = (AliTrackReference*)trackrefs->At(i);
 673         if (tref->GetTrack() != tlab) continue;
 674         x = tref->X();
 675         y = tref->Y();
 676         z = tref->Z();
 677         printf("Track Ref: x %f y %f z %f\n",tref->X(),tref->Y(),tref->Z());
 678
 679         hxzTR->Fill(x,z);
 680         hxyTR->Fill(x,y);
 681         for (Int_t j = 1; j < 10; j++)
 682          {
 683             hxyTR->Fill(x, y+j*0.1);
 684             hxyTR->Fill(x, y-j*0.1);
 685             hxyTR->Fill(x+j*0.1,y);
 686             hxyTR->Fill(x-j*0.1,y);
 687
 688             hxzTR->Fill(x,z-j*0.1);
 689             hxzTR->Fill(x,z+j*0.1);
 690             hxzTR->Fill(x-j*0.1,z);
 691             hxzTR->Fill(x+j*0.1,z);
 692          }
 693       }
 694       break;
 695     }
 696   hxz->Draw("");
 697 //  hxzt->Draw("same");
 698   hxzTR->Draw("same");
 699
 700   delete rl;
 701 }
 702