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