ANALYSIS/AliTrackPoints.cxx

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