Version number ++
[u/mrichter/AliRoot.git] / 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 "AliESD.h"
45 #include "AliRunLoader.h"
46 #include "AliTPCtrack.h"
47 #include "TTree.h"
48 #include "TBranch.h"
49 #include "TH2D.h"
50
51 ClassImp(AliTrackPoints)
52
53 Int_t AliTrackPoints::fgDebug = 0;
54
55 AliTrackPoints::AliTrackPoints():
56  fN(0),
57  fX(0x0),
58  fY(0x0),
59  fZ(0x0)
60 {
61   //constructor
62 }
63
64 /***************************************************************/
65
66 AliTrackPoints::AliTrackPoints(AliTrackPoints::ETypes type, AliESDtrack* track, Float_t mf):
67  fN(0),
68  fX(0x0),
69  fY(0x0),
70  fZ(0x0)
71 {
72   //constructor 
73   //tupe -  what kind of track points should be calculated
74   //mf - magnetic field in [kG] = [T]*10.0
75   switch (type)
76    {
77      case kITS:
78        fN = 6;
79        fX = new Float_t[fN];
80        fY = new Float_t[fN];
81        fZ = new Float_t[fN];
82        MakeITSPoints(track);
83        break;
84
85      case kITSInnerFromVertexOuterFromTPC:
86        fN = 6;
87        fX = new Float_t[fN];
88        fY = new Float_t[fN];
89        fZ = new Float_t[fN];
90        MakeITSPointsInnerFromVertexOuterFromTPC(track,mf);
91        break;
92
93      default:
94        Info("AliTrackPoints","Not recognized type");
95    }
96    
97 }
98
99 /***************************************************************/
100
101 AliTrackPoints::AliTrackPoints(Int_t n, AliESDtrack* track, Float_t mf, Float_t dr, Float_t r0):
102  fN(n),
103  fX(new Float_t[fN]),
104  fY(new Float_t[fN]),
105  fZ(new Float_t[fN])
106 {
107   //constructor
108   //mf - magnetic field in kG - needed to calculate curvature out of Pt
109   //r0 - starting radius
110   //dr - calculate points every dr cm, default every 30cm
111   if (track == 0x0)
112    {
113      Error("AliTrackPoints","ESD track is null");
114      fN = 0;
115      delete [] fX;
116      delete [] fY;
117      delete [] fZ;
118      fX = fY = fZ = 0x0;
119      return;
120    }
121
122   if ( ((track->GetStatus() & AliESDtrack::kTPCrefit) == kFALSE)&& 
123        ((track->GetStatus() & AliESDtrack::kTPCin) == kFALSE)  )
124    {
125      //could happend: its stand alone tracking
126      AliDebug(3,"This ESD track does not contain TPC information");
127        
128      fN = 0;
129      delete [] fX;
130      delete [] fY;
131      delete [] fZ;
132      fX = fY = fZ = 0x0;
133      
134      return;
135    }
136   
137   Double_t alpha,x;
138   Double_t par[5];
139   track->GetInnerExternalParameters(alpha,x,par);//get properties of the track
140   if (par[4] == 0)
141    {
142      Error("AliTrackPoints","This ESD track seem not to contain TPC information (curv is 0)");
143      return;
144    }
145
146   if (mf == 0.0)
147    {
148      Error("AliTrackPoints","Zero Magnetic field passed as parameter.");
149      return;
150    }
151    
152   //Double_t alpha = track->GetInnerAlpha();
153   Double_t cc = 1000./0.299792458/mf;//conversion constant
154   Double_t c=par[4]/cc;
155   
156   MakePoints(dr,r0,x,par,c,alpha);
157   
158 }
159
160 /***************************************************************/
161
162 AliTrackPoints::AliTrackPoints(Int_t n, AliTPCtrack* track, Float_t dr, Float_t r0):
163  fN(n),
164  fX(new Float_t[fN]),
165  fY(new Float_t[fN]),
166  fZ(new Float_t[fN])
167 {
168   //constructor
169   //r0 starting radius
170   //dr - calculate points every dr cm, default every 30cm
171   if (track == 0x0)
172    {
173      Error("AliTrackPoints","TPC track is null");
174      fN = 0;
175      delete [] fX;
176      delete [] fY;
177      delete [] fZ;
178      fX = fY = fZ = 0x0;
179      return;
180    }
181   track->PropagateTo(r0);
182
183  //* This formation is now fixed in the following way:                         *
184  //*      external param0:   local Y-coordinate of a track (cm)                *
185  //*      external param1:   local Z-coordinate of a track (cm)                *
186  //*      external param2:   local sine of the track momentum azimuth angle    *
187  //*      external param3:   tangent of the track momentum dip angle           *
188  //*      external param4:   1/pt (1/(GeV/c))                                  *
189     
190   Double_t x = 0;
191   Double_t par[5];
192   track->GetExternalParameters(x,par);     //get properties of the track
193   
194   Double_t alpha = track->GetAlpha();
195   Double_t c=track->GetC();
196   MakePoints(dr,r0,x,par,c,alpha);
197 }
198
199 /***************************************************************/
200
201 AliTrackPoints::~AliTrackPoints()
202 {
203   //destructor
204   delete [] fX;
205   delete [] fY;
206   delete [] fZ;
207 }
208
209 /***************************************************************/
210
211 void AliTrackPoints::MakePoints( Float_t dr, Float_t r0, Double_t x, Double_t* par, Double_t c, Double_t alpha)
212 {
213   //Calculates points starting at radius r0
214   //spacing every dr (in radial direction)
215   // according to track parameters
216   // x - position in sector local reference frame. x i parallel to R and sector is symmetric with respect to x
217   // par - track external parameters; array with 5 elements;  look at AliTPCtrack.h or AliESDtrack.h for their meaning
218   // c - track curvature
219   // alpha - sector's rotation angle (phi) == angle needed for local to global transformation
220   
221   Double_t y = par[0];
222   Double_t z0 = par[1];
223       
224   Double_t phi0local = TMath::ATan2(y,x);
225   Double_t phi0global = phi0local + alpha;
226   
227   if (phi0local<0) phi0local+=2*TMath::Pi();
228   if (phi0local>=2.*TMath::Pi()) phi0local-=2*TMath::Pi();
229   
230   if (phi0global<0) phi0global+=2*TMath::Pi();
231   if (phi0global>=2.*TMath::Pi()) phi0global-=2*TMath::Pi();
232   
233   Double_t r = TMath::Hypot(x,y);
234   
235
236   AliDebug(9,Form("Radius0 %f, Real Radius %f",r0,r)); 
237   
238   AliDebug(5,Form("Phi Global at first padraw %f, Phi locat %f",phi0global,phi0local));
239     
240   Double_t eta = x*c - par[2] ;//par[2] = fX*C - eta; eta==fP2 ; C==fP4
241   
242   //this calculattions are assuming new (current) model 
243   Double_t tmp = par[2];
244   tmp = 1. - tmp*tmp;
245   tmp = c*y + TMath::Sqrt(tmp);
246   Double_t dca=(TMath::Hypot(eta,tmp) - 1. )/TMath::Abs(c);
247
248   //Here is old model Cold=Cnew/2.
249   Double_t dcasq = dca*dca;
250   Double_t c2 = c/2.;
251   Double_t cst1 = (1.+c2*dca)*dca;//first constant
252   Double_t cst2 = 1. + 2.*c2*dca;//second constant
253       
254   Double_t factorPhi0 = TMath::ASin((c2*r + cst1/r)/cst2);
255   Double_t factorZ0 = TMath::ASin(c2*TMath::Sqrt((r*r-dcasq)/cst2))*par[3]/c2;
256  
257   for(Int_t i = 0; i<fN; i++)
258    {
259      Double_t rc = r0 + i*dr;
260      Double_t ftmp = (c2*rc + cst1/rc)/cst2;
261      if (ftmp > 1.0) 
262       {
263         AliDebug(1,Form("ASin argument > 1 %f:",ftmp));
264         ftmp=1.0;
265       }
266      else if (ftmp < -1.0) 
267       {
268         AliDebug(1,Form("ASin argument < -1 %f:",ftmp));
269         ftmp=-1.0;
270       }
271       
272      Double_t factorPhi = TMath::ASin( ftmp );//factor phi od rc
273      Double_t phi = phi0global + factorPhi - factorPhi0;
274      
275      ftmp = (rc*rc-dcasq)/cst2;
276      if (ftmp < 0.0) 
277       {
278         AliDebug(1,Form("Sqrt argument < 0: %f",ftmp));
279         ftmp=0.0;
280       }
281      
282      ftmp = c2*TMath::Sqrt(ftmp);
283      if (ftmp > 1.0) 
284       {
285         AliDebug(1,Form("ASin argument > 1: %f",ftmp));
286         ftmp=1.0;
287       }
288      else if (ftmp < -1.0) 
289       {
290         AliDebug(2,Form("ASin argument < -1: %f",ftmp));
291         ftmp=-1.0;
292       }
293      Double_t factorZ = TMath::ASin(ftmp)*par[3]/c2;
294      fZ[i] = z0 + factorZ - factorZ0;
295      fX[i] = rc*TMath::Cos(phi);
296      fY[i] = rc*TMath::Sin(phi);
297      
298      AliDebug(3,Form("AliTrackPoints","X %f Y %f Z %f R asked %f R obtained %f",
299              fX[i],fY[i],fZ[i],rc,TMath::Hypot(fX[i],fY[i])));
300    }
301 }
302
303 /***************************************************************/
304
305 void AliTrackPoints::MakeITSPoints(AliESDtrack* track)
306 {
307 //Calculates points in ITS
308 // z=R*Pz/Pt
309  AliITStrackV2 itstrack(*track,kTRUE);
310  Double_t x,y,z;
311  static const Double_t kR[6] = {4.0, 7.0, 14.9, 23.8, 39.1, 43.6};
312  for (Int_t i = 0; i < 6; i++)
313   {
314     itstrack.GetGlobalXYZat(kR[i],x,y,z);
315     fX[i] = x;
316     fY[i] = y;
317     fZ[i] = z;
318 //    Info("MakeITSPoints","X %f Y %f Z %f R asked %f R obtained %f",
319 //             fX[i],fY[i],fZ[i],r[i],TMath::Hypot(fX[i],fY[i]));
320   }   
321  
322 }
323
324 /***************************************************************/
325
326 void AliTrackPoints::MakeITSPointsInnerFromVertexOuterFromTPC(AliESDtrack* track, Float_t mf)
327 {
328 //makes trackpoints for ITS  
329 //for 3 inner layers calculates out of the vector at vertex
330 //for 3 outer ---------------//------------------ at inner TPC  
331
332  static const Double_t kR[6] = {4.0, 7.0, 14.9, 23.8, 39.1, 43.6};
333  AliITStrackV2 itstrack(*track,kTRUE);
334  Double_t x,y,z;
335  for (Int_t i = 0; i < 3; i++)
336   {
337     itstrack.GetGlobalXYZat(kR[i],x,y,z);
338     fX[i] = x;
339     fY[i] = y;
340     fZ[i] = z;
341     AliDebug(3,Form("X %f Y %f Z %f R asked %f R obtained %f",
342                 fX[i],fY[i],fZ[i],kR[i],TMath::Hypot(fX[i],fY[i])));
343   }   
344  
345  for (Int_t i = 3; i < 6; i++)
346   {
347     Float_t ax,ay,az;
348     AliTrackPoints tmptp(1,track,mf,0,kR[i]);
349     tmptp.PositionAt(0,ax,ay,az);
350     fX[i] = ax;
351     fY[i] = ay;
352     fZ[i] = az;
353     AliDebug(3,Form("X %f Y %f Z %f R asked %f R obtained %f",
354                 fX[i],fY[i],fZ[i],kR[i],TMath::Hypot(fX[i],fY[i])));
355   }
356  
357 }
358     
359 /***************************************************************/
360
361 void AliTrackPoints::PositionAt(Int_t n, Float_t &x,Float_t &y,Float_t &z)
362 {
363   //returns position at point n
364   if ((n<0) || (n>=fN))
365    {
366      Error("PositionAt","Point %d out of range",n);
367      return;
368    }
369   
370   x = fX[n];
371   y = fY[n];
372   z = fZ[n];
373   AliDebug(2,Form("n %d; X %f; Y %f; Z %f",n,x,y,z));
374
375 }
376
377 /***************************************************************/
378
379 void AliTrackPoints::Move(Float_t x, Float_t y, Float_t z)
380 {
381 //Moves all points about vector
382  for (Int_t i = 0; i<fN; i++)
383    {
384      fX[i]+=x;
385      fY[i]+=y;
386      fZ[i]+=z;
387    }   
388 }
389
390 /***************************************************************/
391
392 Double_t AliTrackPoints::AvarageDistance(const AliTrackPoints& tr)
393 {
394   //returns the aritmethic avarage distance between two tracks
395 //  Info("AvarageDistance","Entered");
396   if ( (fN <= 0) || (tr.fN <=0) )
397    {
398      AliDebug(1,"One of tracks is empty");
399      return -1;
400    }
401
402   if (fN != tr.fN)
403    {
404      Warning("AvarageDistance","Number of points is not equal");
405      return -1;
406    }
407    
408   Double_t sum = 0;
409   for (Int_t i = 0; i<fN; i++)
410    {
411      AliDebug(10,Form("radii: %f %f",TMath::Hypot(fX[i],fY[i]),TMath::Hypot(tr.fX[i],tr.fY[i])));
412 //       Float_t r1sq = fX[i]*fX[i]+fY[i]*fY[i];
413 //       Float_t r2sq = tr.fX[i]*tr.fX[i]+tr.fY[i]*tr.fY[i];
414
415       
416      Double_t dx = fX[i]-tr.fX[i];
417      Double_t dy = fY[i]-tr.fY[i];
418      Double_t dz = fZ[i]-tr.fZ[i];
419      sum+=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
420      
421      AliDebug(2,Form("Diff: x ,y z: %f , %f, %f",dx,dy,dz));
422      AliDebug(2,Form("xxyyzz %f %f %f %f %f %f",
423             fX[i],tr.fX[i],fY[i],tr.fY[i],fZ[i],tr.fZ[i]));
424    }
425    
426   Double_t retval = sum/((Double_t)fN);
427   AliDebug(1,Form("Avarage distance is %f.",retval));
428
429   return retval;
430 }
431
432 /***************************************************************/
433
434 void AliTrackPoints::Print(Option_t* /*option*/) const
435 {
436   // Prints the coordinates of the track points
437   Info("Print","There is %d points",fN);
438   for(Int_t i = 0; i < fN; i++)
439    {
440      Info("Print","%d: %f %f %f",i,fX[i],fY[i],fZ[i]);
441    }
442
443 }
444
445 /***************************************************************/
446
447 void AliTrackPoints::Testesd(Int_t entr,const char* fname )
448 {
449   // This is a former macro converted to member function.
450   // It is used to test the functionality of the class using ESD tracks
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 nTrackPoints = 170;
476   AliTrackPoints* tp = new AliTrackPoints(nTrackPoints,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<nTrackPoints;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 void AliTrackPoints::Testtpc(Int_t entr)
574 {
575   // This is a former macro converted to member function.
576   // It is used to test the functionality of the class using TPC tracks
577   delete gAlice;
578   gAlice = 0x0;
579   AliRunLoader* rl = AliRunLoader::Open();
580   AliLoader* l = rl->GetLoader("TPCLoader");
581   rl->LoadgAlice();
582   AliKalmanTrack::SetFieldMap(rl->GetAliRun()->Field());
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 nTrackPoints = 160;
589   AliTrackPoints* tp = new AliTrackPoints(nTrackPoints,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<nTrackPoints;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