HBTAN/AliHBTTrackPoints.cxx

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