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