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