fix coverity
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackerDebug.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: AliTRDtrackerDebug.cxx 23810 2008-02-08 09:00:27Z hristov $ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //  Tracker debug streamer                                                   //
21 //                                                                           //
22 //  Authors:                                                                 //
23 //    Alex Bercuci <A.Bercuci@gsi.de>                                        //
24 //    Markus Fasel <M.Fasel@gsi.de>                                          //
25 //                                                                           //
26 ///////////////////////////////////////////////////////////////////////////////
27
28 #include "TFile.h"
29 #include "TTree.h"
30 #include "TTreeStream.h"
31 #include "TLinearFitter.h"
32 #include "TGraph.h"
33 #include "TCanvas.h"
34 #include "TMath.h"
35
36 #include "AliLog.h"
37 #include "AliRieman.h"
38
39 #include "AliTRDgeometry.h"
40 #include "AliTRDtrackV1.h"
41 #include "AliTRDseedV1.h"
42 #include "AliTRDcluster.h"
43 #include "AliTRDgeometry.h"
44
45 #include "AliTRDtrackerDebug.h"
46
47 ClassImp(AliTRDtrackerDebug)
48
49 Int_t AliTRDtrackerDebug::fgEventNumber = 0;
50 Int_t AliTRDtrackerDebug::fgTrackNumber = 0;
51 Int_t AliTRDtrackerDebug::fgCandidateNumber = 0;
52
53 //____________________________________________________
54 AliTRDtrackerDebug::AliTRDtrackerDebug() : AliTRDtrackerV1()
55   ,fOutputStreamer(NULL)
56   ,fTree(NULL)
57   ,fTracklet(NULL)
58   ,fTrack(NULL)
59   ,fNClusters(0)
60   ,fAlpha(0.)
61 {
62         //
63   // Default constructor
64   //
65   fOutputStreamer = new TTreeSRedirector("TRD.Debug.root");
66 }
67
68 //____________________________________________________
69 AliTRDtrackerDebug::~AliTRDtrackerDebug()
70 {
71   // destructor
72   
73   delete fOutputStreamer;
74 }
75
76
77 //____________________________________________________
78 void AliTRDtrackerDebug::Draw(Option_t *) 
79 {
80 // steer draw function
81 }
82
83
84 //____________________________________________________
85 Bool_t AliTRDtrackerDebug::Init()
86 {
87 // steer linking data for various debug streams 
88   fTrack = new AliTRDtrackV1();
89   fTree->SetBranchAddress("ncl", &fNClusters);
90   fTree->SetBranchAddress("track.", &fTrack);
91   return kTRUE;
92 }
93
94 //____________________________________________________
95 Bool_t AliTRDtrackerDebug::Open(const char *method)
96 {
97   // Connect to the tracker debug file
98   
99   TDirectory *savedir = gDirectory; 
100   TFile::Open("TRD.TrackerDebugger.root");
101   fTree = (TTree*)gFile->Get(method);
102   if(!fTree){
103     AliInfo(Form("Can not find debug stream for the %s method.\n", method));
104     savedir->cd();
105     return kFALSE;
106   }
107   savedir->cd();
108   return kTRUE;
109 }
110
111 //____________________________________________________
112 Int_t AliTRDtrackerDebug::Process()
113 {
114 // steer debug process threads
115   
116   for(int it = 0; it<fTree->GetEntries(); it++){
117     if(!fTree->GetEntry(it)) continue;
118     if(!fNClusters) continue;
119     fAlpha = fTrack->GetAlpha();
120     //printf("Processing track %d [%d] ...\n", it, fNClusters);
121     ResidualsTrackletsTrack();
122
123     const AliTRDseedV1 *tracklet = NULL;
124     for(int ip = 5; ip>=0; ip--){
125       if(!(tracklet = fTrack->GetTracklet(ip))) continue;
126       if(!tracklet->GetN()) continue;
127       
128       ResidualsClustersTrack(tracklet);
129       ResidualsClustersTracklet(tracklet);
130       ResidualsClustersParametrisation(tracklet);
131     }
132   }
133   return kTRUE;
134 }
135
136
137 //____________________________________________________
138 void AliTRDtrackerDebug::ResidualsClustersTrack(const AliTRDseedV1 *tracklet)
139 {
140 // Calculate averange distances from clusters to the TRD track  
141   
142   Double_t x[3]; 
143   AliTRDcluster *c = NULL;
144   for(int ic=0; ic<35/*AliTRDseed:knTimebins*/; ic++){
145     if(!(c = tracklet->GetClusters(ic))) continue;
146     Double_t xc = c->GetX(), yc = c->GetY(), zc = c->GetZ();
147
148     // propagate track to cluster 
149     if(!PropagateToX(*fTrack, xc, 2.)) continue; 
150     fTrack->GetXYZ(x);
151     
152     // transform to local tracking coordinates
153     //Double_t xg =  x[0] * TMath::Cos(fAlpha) + x[1] * TMath::Sin(fAlpha); 
154     Double_t yg = -x[0] * TMath::Sin(fAlpha) + x[1] * TMath::Cos(fAlpha);
155
156     // apply tilt pad correction
157     yc+= (zc - x[2]) * tracklet->GetTilt();
158     
159     Double_t dy = yc-yg;
160
161     TTreeSRedirector &cstreamer = *fOutputStreamer;
162     cstreamer << "ResidualsClustersTrack"
163       << "c.="   << c
164       << "dy="   << dy
165       << "\n";
166   }
167 }
168
169 //____________________________________________________
170 void AliTRDtrackerDebug::ResidualsClustersTracklet(const AliTRDseedV1 *tracklet) const
171 {
172 // Calculates distances from clusters to tracklets
173   
174   Double_t x0 = tracklet->GetX0(), 
175           y0 = tracklet->GetYfit(0), 
176           ys = tracklet->GetYfit(1);
177           //z0 = tracklet->GetZfit(0), 
178           //zs = tracklet->GetZfit(1);
179   
180   AliTRDcluster *c = NULL;
181   for(int ic=0; ic<35/*AliTRDseed:knTimebins*/; ic++){
182     if(!(c = tracklet->GetClusters(ic))) continue;
183     Double_t xc = c->GetX(), yc = c->GetY()/*, zc = c->GetZ()*/;
184     Double_t dy = yc- (y0-(x0-xc)*ys);
185
186     //To draw  use : 
187     //ResidualsClustersTracklet->Draw("TMath::Abs(10.*dy):TMath::ATan(ys)*TMath::RadToDeg()>>h(20, -40, 40)", "", "prof");
188     TTreeSRedirector &cstreamer = *fOutputStreamer;
189     cstreamer << "ResidualsClustersTracklet"
190       << "c.="   << c
191       << "ys="   << ys
192       << "dy="   << dy
193       << "\n";
194   }
195 }
196
197 //____________________________________________________
198 void AliTRDtrackerDebug::ResidualsClustersParametrisation(const AliTRDseedV1 *tracklet) const
199 {
200 // Calculates distances from clusters to Rieman fit.
201   
202   // store cluster positions
203   Double_t x0 = tracklet->GetX0();
204   AliTRDcluster *c = NULL;
205   
206   Double_t x[2]; Int_t ncl, mcl, jc;
207   TLinearFitter fitter(3, "hyp2");
208   for(int ic=0; ic<35/*AliTRDseed:knTimebins*/; ic++){
209     if(!(c = tracklet->GetClusters(ic))) continue;
210     Double_t xc = c->GetX(), yc = c->GetY()/*, zc = c->GetZ()*/;
211     
212     jc = ic; ncl = 0; mcl=0; fitter.ClearPoints();
213     while(ncl<6){
214       // update index
215       mcl++;
216       jc = ic + ((mcl&1)?-1:1)*(mcl>>1);
217
218       if(jc<0 || jc>=35) continue;
219       if(!(c = tracklet->GetClusters(jc))) continue;
220
221       x[0] = c->GetX()-x0;
222       x[1] = x[0]*x[0];
223       fitter.AddPoint(x, c->GetY(), c->GetSigmaY2());
224       ncl++;
225     }
226     fitter.Eval();
227     Double_t dy = yc - fitter.GetParameter(0) -fitter.GetParameter(1) * (xc-x0) - fitter.GetParameter(2)* (xc-x0)*(xc-x0); 
228   
229     TTreeSRedirector &cstreamer = *fOutputStreamer;
230     cstreamer << "ResidualsClustersParametrisation"
231       << "dy="   << dy
232       << "\n";
233   }
234 }
235
236
237 //____________________________________________________
238 void AliTRDtrackerDebug::ResidualsTrackletsTrack() const
239 {
240 // Calculates distances from tracklets to the TRD track.
241   
242   if(fTrack->GetNumberOfTracklets() < 6) return;
243
244   // build a working copy of the tracklets attached to the track 
245   // and initialize working variables fX, fY and fZ
246   //AliTRDseedV1 tracklet[6] = {NULL, NULL, NULL, NULL, NULL, NULL};
247   AliTRDseedV1 tracklet[6];
248   const AliTRDseedV1 *ctracklet = NULL;
249   for(int ip = 0; ip<6; ip++){
250     if(!(ctracklet = fTrack->GetTracklet(ip))) continue;
251     tracklet[ip] = (*ctracklet); 
252 //              Double_t x0 = tracklet[ip].GetX0();
253 //              for(int ic=0; ic<AliTRDseedV1:knTimebins; ic++){
254 //                      if(!(c = tracklet[ip].GetClusters(ic))) continue;
255 //                      Double_t xc = c->GetX(), yc = c->GetY(), zc = c->GetZ();
256 //                      tracklet[ip].SetX(ic, xc-x0);
257 //                      tracklet[ip].SetY(ic, yc);
258 //                      tracklet[ip].SetZ(ic, zc);
259 //              }
260   }
261   
262   // Do a Rieman fit (with tilt correction) for all tracklets 
263   // except the one which is tested. 
264   // (Based on AliTRDseed::IsOK() return false)
265   for(int ip=0; ip<6; ip++){
266     // reset tracklet to be tested
267     Double_t x0 = tracklet[ip].GetX0();
268     new(&tracklet[ip]) AliTRDseedV1();
269     tracklet[ip].SetX0(x0);
270
271     // fit Rieman with tilt correction
272     AliTRDtrackerV1::FitRiemanTilt(NULL, &tracklet[0], kTRUE);
273
274     // make a copy of the fit result
275     Double_t 
276       y0   = tracklet[ip].GetYref(0),
277       dydx = tracklet[ip].GetYref(1),
278       z0   = tracklet[ip].GetZref(0),
279       dzdx = tracklet[ip].GetZref(1);
280
281     // restore tracklet
282     tracklet[ip] = (*fTrack->GetTracklet(ip)); 
283 //              for(int ic=0; ic<AliTRDseedV1:knTimebins; ic++){
284 //                      if(!(c = tracklet[ip].GetClusters(ic))) continue;
285 //                      Double_t xc = c->GetX(), yc = c->GetY(), zc = c->GetZ();
286 //                      tracklet[ip].SetX(ic, xc-x0);
287 //                      tracklet[ip].SetY(ic, yc);
288 //                      tracklet[ip].SetZ(ic, zc);
289 //              }               
290     
291     // fit clusters
292     AliTRDseedV1 ts(tracklet[ip]);
293     ts.SetYref(0, y0); ts.SetYref(1, dydx);
294     ts.SetZref(0, z0); ts.SetZref(1, dzdx);
295     ts.Fit(kTRUE);
296
297     // save results for plotting
298     Int_t plane   = tracklet[ip].GetPlane();
299     Double_t dy   = tracklet[ip].GetYfit(0) - ts.GetYfit(0);
300     Double_t tgy  = tracklet[ip].GetYfit(1);
301     Double_t dtgy = (tracklet[ip].GetYfit(1) - ts.GetYfit(1))/(1. + tracklet[ip].GetYfit(1) * ts.GetYfit(1));
302     Double_t dz   = tracklet[ip].GetZfit(0) - ts.GetZfit(0);
303     Double_t tgz  = tracklet[ip].GetZfit(1);
304     Double_t dtgz = (tracklet[ip].GetZfit(1) - ts.GetZfit(1))/(1. + tracklet[ip].GetZfit(1) * ts.GetZfit(1));
305     TTreeSRedirector &cstreamer = *fOutputStreamer;
306     cstreamer << "ResidualsTrackletsTrack"
307       << "ref.="   << &tracklet[ip]
308       << "fit.="   << &ts
309       << "plane="  << plane
310       << "dy="     << dy
311       << "tgy="    << tgy
312       << "dtgy="   << dtgy
313       << "dz="     << dz
314       << "tgz="    << tgz
315       << "dtgz="   << dtgz
316       << "\n";
317   }
318 }
319
320 //____________________________________________________
321 void AliTRDtrackerDebug::AnalyseFindable(Char_t *treename){
322 //
323 // Calculates the number of findable tracklets defined as the number of tracklets
324 // per track candidate where the tan phi_tracklet is below 0.15 (maximum inclination
325 // in y-direction.
326 // 
327 // Parameters:  -the treename (this method can be used for all trees which store the
328 //                               tracklets
329 // Output:              -void
330 //
331 // A new TTree containing the number of findable tracklets and the number of clusters
332 // attached to the full track is stored to disk
333 //
334   // Link the File
335   TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
336   fTree = (TTree *)(debfile->Get(treename));
337   if(!fTree){
338     AliError(Form("Tree %s not found in file TRDdebug.root. Abborting!", treename));
339     debfile->Close();
340     return;
341   }
342   
343   AliTRDseedV1 *tracklets[kNPlanes];
344   for(Int_t iPlane = 0; iPlane < AliTRDtrackerV1::kNPlanes; iPlane++)
345     tracklets[iPlane] = NULL;
346   for(Int_t iPlane = 0; iPlane < kNPlanes; iPlane++)
347     fTree->SetBranchAddress(Form("S%d.", iPlane), &tracklets[iPlane]);
348   fTree->SetBranchAddress("EventNumber", &fgEventNumber);
349   fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
350   
351   Int_t findable = 0, nClusters = 0;
352   Int_t nEntries = fTree->GetEntriesFast();
353   for(Int_t iEntry = 0; iEntry < nEntries; iEntry++){
354     printf("Entry %d\n", iEntry);
355     fTree->GetEntry(iEntry);
356     findable = 0;
357     nClusters = 0;
358     // Calculate Findable
359     for(Int_t iPlane = 0; iPlane < kNPlanes; iPlane++){
360       if (TMath::Abs(tracklets[iPlane]->GetYref(0) / tracklets[iPlane]->GetX0()) < 0.15) findable++;
361       if (!tracklets[iPlane]->IsOK()) continue;
362       nClusters += tracklets[iPlane]->GetN2();
363     }
364     
365     // Fill Histogramms
366     TTreeSRedirector &cstreamer = *fOutputStreamer;
367     cstreamer << "AnalyseFindable"
368       << "EventNumber="         << fgEventNumber
369       << "CandidateNumber="     << fgCandidateNumber
370       << "Findable="                    << findable
371       << "NClusters="                   << nClusters
372       << "\n";
373   }
374 }
375 //____________________________________________________
376 void AliTRDtrackerDebug::AnalyseTiltedRiemanFit(){
377 //
378 // Creating a Data Set for the method FitTiltedRieman containing usefull variables
379 // Each variable can be addressed to tracks later. Data can be processed later.
380 //
381 // Parameters: -
382 // Output:     -
383 //
384 // TODO: Match everything with Init and Process
385 //
386   TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
387   fTree = (TTree *)(debfile->Get("MakeSeeds2"));
388   if(!fTree) return;
389   Int_t nEntries = fTree->GetEntries();
390   TLinearFitter *tiltedRiemanFitter = NULL;
391   fTree->SetBranchAddress("FitterT.", &tiltedRiemanFitter);
392   fTree->SetBranchAddress("EventNumber", &fgEventNumber);
393   fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
394   for(Int_t entry = 0; entry < nEntries; entry++){
395     fTree->GetEntry(entry);
396     Double_t a = tiltedRiemanFitter->GetParameter(0);
397     Double_t b = tiltedRiemanFitter->GetParameter(1);
398     Double_t c = tiltedRiemanFitter->GetParameter(2);
399     Double_t offset = tiltedRiemanFitter->GetParameter(3);
400     Double_t slope  = tiltedRiemanFitter->GetParameter(4);
401     Float_t radius = GetTrackRadius(a, b, c);
402     Float_t curvature = GetTrackCurvature(a, b, c);
403     Float_t dca = GetDCA(a, b, c);
404     TTreeSRedirector &cstreamer = *fOutputStreamer;
405     cstreamer << "AnalyseTiltedRiemanFit"
406     << "EventNumber="           << fgEventNumber
407     << "CandidateNumber=" << fgCandidateNumber
408     << "Radius="                                        << radius
409     << "Curvature="                             << curvature
410     << "DCA="                                                   << dca
411     << "Offset="                                        << offset
412     << "Slope="                                         << slope
413     << "\n";
414   }
415 }
416
417 //____________________________________________________
418 Float_t AliTRDtrackerDebug::GetTrackRadius(Float_t a, Float_t b, Float_t c) const {
419 //
420 // Calculates the track radius using the parameters given by the tilted Rieman fit 
421 //
422 // Parameters: The three parameters from the Rieman fit
423 // Output:     The track radius
424 //
425   Float_t radius = 0;
426   if(1.0 + b*b - c*a > 0.0)
427     radius = TMath::Sqrt(1.0 + b*b - c*a )/a;
428   return radius;
429 }
430
431 //____________________________________________________
432 Float_t AliTRDtrackerDebug::GetTrackCurvature(Float_t a, Float_t b, Float_t c) const {
433 //
434 // Calculates the track curvature using the parameters given by the linear fitter 
435 //
436 // Parameters:  the three parameters from the tilted Rieman fitter
437 // Output:              the full track curvature
438 //
439   Float_t curvature =  1.0 + b*b - c*a;
440   if (curvature > 0.0) 
441     curvature  =  a / TMath::Sqrt(curvature);
442   return curvature;
443 }
444
445 //____________________________________________________
446 Float_t AliTRDtrackerDebug::GetDCA(Float_t a, Float_t b, Float_t c) const {
447 //
448 // Calculates the Distance to Clostest Approach for the Vertex using the paramters
449 // given by the tilted Rieman fit 
450 //
451 // Parameters: the three parameters from the tilted Rieman fitter
452 // Output:     the Distance to Closest Approach
453 //
454   Float_t dca  =  0.0;
455   if (1.0 + b*b - c*a > 0.0) {
456     dca = -c / (TMath::Sqrt(1.0 + b*b - c*a) + TMath::Sqrt(1.0 + b*b));
457   }
458   return dca;
459 }
460
461 //____________________________________________________
462 void AliTRDtrackerDebug::AnalyseMinMax()
463 {
464   // Development function related to the old tracking code
465   TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
466   if(!debfile){
467     AliError("File TRD.TrackerDebug.root not found!");
468     return; 
469   }
470   fTree = (TTree *)(debfile->Get("MakeSeeds0"));
471   if(!fTree){
472     AliError("Tree MakeSeeds0 not found in File TRD.TrackerDebug.root.");
473     return;
474   }
475   AliTRDseedV1 *cseed[4] = {NULL, NULL, NULL, NULL};
476   AliTRDcluster *c[4] = {NULL, NULL, NULL, NULL};
477   for(Int_t il = 0; il < 4; il++){
478     fTree->SetBranchAddress(Form("Seed%d.", il),        &cseed[il]);
479     fTree->SetBranchAddress(Form("c%d.",il), &c[il]);
480   }
481   fTree->SetBranchAddress("CandidateNumber",    &fgCandidateNumber);
482   fTree->SetBranchAddress("EventNumber",        &fgEventNumber);
483   Int_t entries = fTree->GetEntries();
484   for(Int_t ientry = 0; ientry < entries; ientry++){
485     fTree->GetEntry(ientry);
486     Float_t minmax[2] = { -100.0,  100.0 };
487     for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
488       Float_t max = c[iLayer]->GetZ() + cseed[iLayer]->GetPadLength() * 0.5 + 1.0 - cseed[iLayer]->GetZref(0);
489       if (max < minmax[1]) minmax[1] = max;
490       Float_t min = c[iLayer]->GetZ()-cseed[iLayer]->GetPadLength() * 0.5 - 1.0 - cseed[iLayer]->GetZref(0);
491       if (min > minmax[0]) minmax[0] = min;
492     }
493     TTreeSRedirector &cstreamer = *fOutputStreamer;
494     cstreamer << "AnalyseMinMaxLayer"
495     << "EventNumber="                           << fgEventNumber
496     << "CandidateNumber="               << fgCandidateNumber
497     << "Min="                                                           << minmax[0]
498     << "Max="                                                           << minmax[1]
499     << "\n";
500   }
501 }
502
503 //____________________________________________________
504 TCanvas* AliTRDtrackerDebug::PlotSeedingConfiguration(const Char_t *direction, Int_t event, Int_t candidate){
505 //
506 // Plots the four seeding clusters, the helix fit and the reference Points for
507 // a given combination consisting of eventnr. and candidatenr.
508 //
509 // Parameters:  -direction (y or z)
510 //                              -Event Nr
511 //              -Candidate that has to be plotted
512 //
513   const Float_t kxmin = 280;
514   const Float_t kxmax = 380;
515   const Float_t kxdelta = (kxmax - kxmin)/1000;
516   
517   if((strcmp(direction, "y") != 0) && (strcmp(direction, "z") != 0)){
518     AliError(Form("Direction %s does not exist. Abborting!", direction));
519     return NULL;
520   }
521
522   TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
523   if(!debfile){
524     AliError("File TRD.TrackerDebug.root not found!");
525     return NULL; 
526   }
527   fTree = (TTree *)(debfile->Get("MakeSeeds0"));
528   if(!fTree){
529     AliError("Tree MakeSeeds0 not found in File TRD.TrackerDebug.root.");
530     return NULL;
531   }
532   
533   TGraph *seedcl = new TGraph(4);
534   TGraph *seedRef = new TGraph(4);
535   TGraph *riemanFit = new TGraph(1000);
536   seedcl->SetMarkerStyle(20);
537   seedcl->SetMarkerColor(kRed);
538   seedRef->SetMarkerStyle(2);
539
540   AliTRDcluster *c[4] = {NULL, NULL, NULL, NULL};
541   AliRieman *rim = NULL;
542   Bool_t found = kFALSE;
543   for(Int_t il = 0; il < 4; il++) fTree->SetBranchAddress(Form("c%d.",il), &c[il]);
544   fTree->SetBranchAddress("EventNumber", &fgEventNumber);
545   fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
546   fTree->SetBranchAddress("RiemanFitter.", &rim);
547   Int_t entries = fTree->GetEntries();
548   for(Int_t entry = 0; entry < entries; entry++){
549     fTree->GetEntry(entry);
550     if(fgEventNumber < event) continue;
551     if(fgEventNumber > event) break;
552     // EventNumber matching: Do the same for the candidate number
553     if(fgCandidateNumber < candidate) continue;
554     if(fgCandidateNumber > candidate) break;
555     found = kTRUE;
556     Int_t nPoints = 0;
557     for(Int_t il = 0; il < 4; il++){
558       Float_t cluster = 0.0, reference = 0.0;
559       if(!strcmp(direction, "y")){
560         cluster = c[il]->GetY();
561         reference = rim->GetYat(c[il]->GetX());
562       }
563       else{
564         cluster = c[il]->GetZ();
565         reference = rim->GetZat(c[il]->GetX());
566       }
567       seedcl->SetPoint(nPoints, cluster, c[il]->GetX());
568       seedRef->SetPoint(nPoints, reference , c[il]->GetX());
569       nPoints++;
570     }
571     // evaluate the fitter Function numerically
572     nPoints = 0;
573     for(Int_t ipt = 0; ipt < 1000; ipt++){
574       Float_t x = kxmin + ipt * kxdelta;
575       Float_t point = 0.0;
576       if(!strcmp(direction, "y"))
577         point = rim->GetYat(x);
578       else
579         point = rim->GetZat(x);
580       riemanFit->SetPoint(nPoints++, point, x);
581     }
582     // We reached the End: break
583     break;
584   }
585   if(found){
586     seedcl->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
587     seedRef->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
588     riemanFit->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
589     TCanvas *c1 = new TCanvas();
590     seedcl->Draw("ap");
591     seedRef->Draw("psame");
592     riemanFit->Draw("lpsame");
593     return c1;
594   }
595   else{
596     AliError(Form("Combination consisting of event %d and candidate %d not found", event, candidate));
597     delete seedcl;
598     delete seedRef;
599     delete riemanFit;
600     return NULL;
601   }
602 }
603
604 //____________________________________________________
605 TCanvas* AliTRDtrackerDebug::PlotFullTrackFit(Int_t event, Int_t candidate, Int_t iteration, const Char_t *direction){
606 //
607 // Plots the tracklets (clusters and reference in y direction) and the fitted function for several iterations
608 // in the function ImproveSeedQuality (default is before ImproveSeedQuality)
609 // 
610 // Parameters: -Event Number
611 //             -Candidate Number
612 //             -Iteration Number in ImproveSeedQuality (default: -1 = before ImproveSeedQuality)
613 //                         -direction (default: y)
614 // Output:     -TCanvas (containing the Picture);
615 //
616   const Float_t kxmin = 280;
617   const Float_t kxmax = 380;
618   const Float_t kxdelta = (kxmax - kxmin)/1000;
619   
620   if(strcmp(direction, "y") && strcmp(direction, "z")){
621     AliError(Form("Direction %s does not exist. Abborting!", direction));
622     return NULL;
623   }
624
625   TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
626   if(!debfile){
627     AliError("File TRD.TrackerDebug.root not found.");
628     return NULL;
629   }
630   TString *treename = NULL;
631   if(iteration > -1)
632     treename = new TString("ImproveSeedQuality");
633   else
634     treename = new TString("MakeSeeds1");
635   fTree = (TTree *)(debfile->Get(treename->Data()));
636   if(!fTree){
637     AliError(Form("Tree %s not found in File TRD.TrackerDebug.root.", treename->Data()));
638     delete treename;
639     return NULL;
640   }
641   delete treename;
642
643   TGraph *fitfun = new TGraph(1000);
644   // Prepare containers
645   Float_t x0[AliTRDtrackerV1::kNPlanes],
646       refP[AliTRDtrackerV1::kNPlanes],
647       clx[AliTRDtrackerV1::kNPlanes * AliTRDtrackerV1::kNTimeBins],
648       clp[AliTRDtrackerV1::kNPlanes * AliTRDtrackerV1::kNTimeBins];
649   Int_t nLayers = 0, ncls = 0;
650   
651   TLinearFitter *fitter = NULL;
652   AliTRDseedV1 *tracklet[6] = {NULL, NULL, NULL, NULL, NULL, NULL};
653   for(Int_t iLayer = 0; iLayer < 6; iLayer++)
654     fTree->SetBranchAddress(Form("S%d.", iLayer), &tracklet[iLayer]);
655   fTree->SetBranchAddress("FitterT.", &fitter);
656   fTree->SetBranchAddress("EventNumber", &fgEventNumber);
657   fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
658   
659   Int_t nEntries = fTree->GetEntriesFast();
660   Bool_t found = kFALSE;
661   for(Int_t entry = 0; entry < nEntries; entry++){
662     fTree->GetEntry(entry);
663     if(fgEventNumber < event) continue;
664     if(fgEventNumber > event) break;
665     // EventNumber matching: Do the same for the candidate number
666     if(fgCandidateNumber < candidate) continue;
667     if(fgCandidateNumber > candidate) break;
668     found = kTRUE;
669     
670     for(Int_t iLayer = 0; iLayer < 6; iLayer++){
671       if(!tracklet[iLayer]->IsOK()) continue;
672       x0[nLayers] = tracklet[iLayer]->GetX0();
673       if(!strcmp(direction, "y"))
674         refP[nLayers] = tracklet[iLayer]->GetYref(0);
675       else
676         refP[nLayers] = tracklet[iLayer]->GetZref(0);
677       nLayers++;
678       AliTRDcluster *cl(NULL);
679       for(Int_t itb = 0; itb < 30; itb++){
680         if(!tracklet[iLayer]->IsUsable(itb)) continue;
681         if(!(cl = tracklet[iLayer]->GetClusters(itb))) continue;
682         
683         if(!strcmp(direction, "y"))
684           clp[ncls] = cl->GetY();
685         else
686           clp[ncls] = cl->GetZ();
687         clx[ncls] = cl->GetX();
688         ncls++;
689       }
690     }
691     // Add function derived by the tilted Rieman fit (Defined by the curvature)
692     Int_t nPoints = 0;
693     if(!strcmp(direction, "y")){
694       Double_t a = fitter->GetParameter(0);
695       Double_t b = fitter->GetParameter(1);
696       Double_t c = fitter->GetParameter(2);
697       Double_t curvature =  1.0 + b*b - c*a;
698       if (curvature > 0.0) {
699         curvature  =  a / TMath::Sqrt(curvature);
700       }
701       // Numerical evaluation of the function:
702       for(Int_t ipt = 0; ipt < 1000; ipt++){
703         Float_t x = kxmin + ipt * kxdelta;
704         Double_t res = (x * a + b);                                                             // = (x - x0)/y0
705         res *= res;
706         res  = 1.0 - c * a + b * b - res;                                       // = (R^2 - (x - x0)^2)/y0^2
707         Double_t y = 0.;
708         if (res >= 0) {
709           res = TMath::Sqrt(res);
710           y    = (1.0 - res) / a;
711         }
712         fitfun->SetPoint(nPoints++, y, x);
713       }
714     }
715     else{
716       Double_t offset   = fitter->GetParameter(3);
717       Double_t slope    = fitter->GetParameter(4);       
718       // calculate the reference x (defined as medium between layer 2 and layer 3)
719       // same procedure as in the tracker code
720       Float_t medx = 0, xref = 0;
721       Int_t startIndex = 5, nDistances = 0;
722       for(Int_t iLayer = 5; iLayer > 0; iLayer--){
723         if(tracklet[iLayer]->IsOK() && tracklet[iLayer - 1]->IsOK()){
724           medx += tracklet[iLayer]->GetX0() - tracklet[iLayer - 1]->GetX0();
725           startIndex = iLayer - 1;
726           nDistances++;
727         }
728       }
729       if(nDistances){
730         medx /= nDistances;
731       }
732       else{
733         Float_t xpos[2];        memset(xpos, 0, sizeof(Float_t) * 2);
734         Int_t ien = 0, idiff = 0;
735         for(Int_t iLayer = 5; iLayer > 0; iLayer--){
736           if(tracklet[iLayer]->IsOK()){
737             xpos[ien++] = tracklet[iLayer]->GetX0();
738             startIndex = iLayer;
739           }
740           if(ien)
741             idiff++;
742           if(ien >=2)
743             break;
744         }
745         medx = (xpos[0] - xpos[1])/idiff;
746       }
747       xref = tracklet[startIndex]->GetX0() + medx * (2.5 - startIndex) - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
748
749       for(Int_t ipt = 0; ipt < 1000; ipt++){
750         Float_t x = kxmin + ipt * kxdelta;
751         Float_t z = offset + slope * (x - xref);
752         fitfun->SetPoint(nPoints++, z, x);
753       }
754     }
755     break;
756   }
757   if(found){
758     TGraph *trGraph             = new TGraph(ncls);
759     TGraph *refPoints   = new TGraph(nLayers);
760     trGraph->SetMarkerStyle(20);
761     trGraph->SetMarkerColor(kRed);
762     refPoints->SetMarkerStyle(21);
763     refPoints->SetMarkerColor(kBlue);
764     // fill the graphs
765     for(Int_t iLayer = 0; iLayer < nLayers; iLayer++)
766       refPoints->SetPoint(iLayer, refP[iLayer], x0[iLayer]);
767     for(Int_t icls = 0; icls < ncls; icls++)
768       trGraph->SetPoint(icls, clp[icls], clx[icls]);
769     TCanvas *c1 = new TCanvas();
770     trGraph->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
771     refPoints->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
772     fitfun->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
773     trGraph->Draw("ap");
774     refPoints->Draw("psame");
775     fitfun->Draw("lpsame");
776     return c1;
777   }
778   else{
779     AliError(Form("Combination consisting of event %d and candidate %d not found", event, candidate));
780     delete fitfun;
781     return NULL;
782   }
783 }
784