fix mem leak and compiler warning
[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 "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