- Move to the new attaching cluster to tracklet procedure in the stand
[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   AliTRDseedV1 tracklet[6] = {0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
246   const AliTRDseedV1 *ctracklet = 0x0;
247   for(int ip = 0; ip<6; ip++){
248     if(!(ctracklet = fTrack->GetTracklet(ip))) continue;
249     tracklet[ip] = (*ctracklet); 
250 //              Double_t x0 = tracklet[ip].GetX0();
251 //              for(int ic=0; ic<AliTRDseedV1:knTimebins; ic++){
252 //                      if(!(c = tracklet[ip].GetClusters(ic))) continue;
253 //                      Double_t xc = c->GetX(), yc = c->GetY(), zc = c->GetZ();
254 //                      tracklet[ip].SetX(ic, xc-x0);
255 //                      tracklet[ip].SetY(ic, yc);
256 //                      tracklet[ip].SetZ(ic, zc);
257 //              }
258   }
259   
260   // Do a Rieman fit (with tilt correction) for all tracklets 
261   // except the one which is tested. 
262   // (Based on AliTRDseed::IsOK() return false)
263   for(int ip=0; ip<6; ip++){
264     // reset tracklet to be tested
265     Double_t x0 = tracklet[ip].GetX0();
266     new(&tracklet[ip]) AliTRDseedV1();
267     tracklet[ip].SetX0(x0);
268
269     // fit Rieman with tilt correction
270     AliTRDtrackerV1::FitRiemanTilt(0x0, &tracklet[0], kTRUE);
271
272     // make a copy of the fit result
273     Double_t 
274       y0   = tracklet[ip].GetYref(0),
275       dydx = tracklet[ip].GetYref(1),
276       z0   = tracklet[ip].GetZref(0),
277       dzdx = tracklet[ip].GetZref(1);
278
279     // restore tracklet
280     tracklet[ip] = (*fTrack->GetTracklet(ip)); 
281 //              for(int ic=0; ic<AliTRDseedV1:knTimebins; ic++){
282 //                      if(!(c = tracklet[ip].GetClusters(ic))) continue;
283 //                      Double_t xc = c->GetX(), yc = c->GetY(), zc = c->GetZ();
284 //                      tracklet[ip].SetX(ic, xc-x0);
285 //                      tracklet[ip].SetY(ic, yc);
286 //                      tracklet[ip].SetZ(ic, zc);
287 //              }               
288     
289     // fit clusters
290     AliTRDseedV1 ts(tracklet[ip]);
291     ts.SetYref(0, y0); ts.SetYref(1, dydx);
292     ts.SetZref(0, z0); ts.SetZref(1, dzdx);
293     ts.Fit(kTRUE);
294
295     // save results for plotting
296     Int_t plane   = tracklet[ip].GetPlane();
297     Double_t dy   = tracklet[ip].GetYfit(0) - ts.GetYfit(0);
298     Double_t tgy  = tracklet[ip].GetYfit(1);
299     Double_t dtgy = (tracklet[ip].GetYfit(1) - ts.GetYfit(1))/(1. + tracklet[ip].GetYfit(1) * ts.GetYfit(1));
300     Double_t dz   = tracklet[ip].GetZfit(0) - ts.GetZfit(0);
301     Double_t tgz  = tracklet[ip].GetZfit(1);
302     Double_t dtgz = (tracklet[ip].GetZfit(1) - ts.GetZfit(1))/(1. + tracklet[ip].GetZfit(1) * ts.GetZfit(1));
303     TTreeSRedirector &cstreamer = *fOutputStreamer;
304     cstreamer << "ResidualsTrackletsTrack"
305       << "ref.="   << &tracklet[ip]
306       << "fit.="   << &ts
307       << "plane="  << plane
308       << "dy="     << dy
309       << "tgy="    << tgy
310       << "dtgy="   << dtgy
311       << "dz="     << dz
312       << "tgz="    << tgz
313       << "dtgz="   << dtgz
314       << "\n";
315   }
316 }
317
318 //____________________________________________________
319 void AliTRDtrackerDebug::AnalyseFindable(Char_t *treename){
320 //
321 // Calculates the number of findable tracklets defined as the number of tracklets
322 // per track candidate where the tan phi_tracklet is below 0.15 (maximum inclination
323 // in y-direction.
324 // 
325 // Parameters:  -the treename (this method can be used for all trees which store the
326 //                               tracklets
327 // Output:              -void
328 //
329 // A new TTree containing the number of findable tracklets and the number of clusters
330 // attached to the full track is stored to disk
331 //
332   // Link the File
333   TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
334   fTree = (TTree *)(debfile->Get(treename));
335   if(!fTree){
336     AliError(Form("Tree %s not found in file TRDdebug.root. Abborting!", treename));
337     debfile->Close();
338     return;
339   }
340   
341   AliTRDseedV1 *tracklets[kNPlanes];
342   for(Int_t iPlane = 0; iPlane < AliTRDtrackerV1::kNPlanes; iPlane++)
343     tracklets[iPlane] = 0x0;
344   for(Int_t iPlane = 0; iPlane < kNPlanes; iPlane++)
345     fTree->SetBranchAddress(Form("S%d.", iPlane), &tracklets[iPlane]);
346   fTree->SetBranchAddress("EventNumber", &fgEventNumber);
347   fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
348   
349   Int_t findable = 0, nClusters = 0;
350   Int_t nEntries = fTree->GetEntriesFast();
351   for(Int_t iEntry = 0; iEntry < nEntries; iEntry++){
352     printf("Entry %d\n", iEntry);
353     fTree->GetEntry(iEntry);
354     findable = 0;
355     nClusters = 0;
356     // Calculate Findable
357     for(Int_t iPlane = 0; iPlane < kNPlanes; iPlane++){
358       if (TMath::Abs(tracklets[iPlane]->GetYref(0) / tracklets[iPlane]->GetX0()) < 0.15) findable++;
359       if (!tracklets[iPlane]->IsOK()) continue;
360       nClusters += tracklets[iPlane]->GetN2();
361     }
362     
363     // Fill Histogramms
364     TTreeSRedirector &cstreamer = *fOutputStreamer;
365     cstreamer << "AnalyseFindable"
366       << "EventNumber="         << fgEventNumber
367       << "CandidateNumber="     << fgCandidateNumber
368       << "Findable="                    << findable
369       << "NClusters="                   << nClusters
370       << "\n";
371   }
372 }
373 //____________________________________________________
374 void AliTRDtrackerDebug::AnalyseTiltedRiemanFit(){
375 //
376 // Creating a Data Set for the method FitTiltedRieman containing usefull variables
377 // Each variable can be addressed to tracks later. Data can be processed later.
378 //
379 // Parameters: -
380 // Output:     -
381 //
382 // TODO: Match everything with Init and Process
383 //
384   TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
385   fTree = (TTree *)(debfile->Get("MakeSeeds2"));
386   if(!fTree) return;
387   Int_t nEntries = fTree->GetEntries();
388   TLinearFitter *tiltedRiemanFitter = 0x0;
389   fTree->SetBranchAddress("FitterT.", &tiltedRiemanFitter);
390   fTree->SetBranchAddress("EventNumber", &fgEventNumber);
391   fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
392   for(Int_t entry = 0; entry < nEntries; entry++){
393     fTree->GetEntry(entry);
394     Double_t a = tiltedRiemanFitter->GetParameter(0);
395     Double_t b = tiltedRiemanFitter->GetParameter(1);
396     Double_t c = tiltedRiemanFitter->GetParameter(2);
397     Double_t offset = tiltedRiemanFitter->GetParameter(3);
398     Double_t slope  = tiltedRiemanFitter->GetParameter(4);
399     Float_t radius = GetTrackRadius(a, b, c);
400     Float_t curvature = GetTrackCurvature(a, b, c);
401     Float_t dca = GetDCA(a, b, c);
402     TTreeSRedirector &cstreamer = *fOutputStreamer;
403     cstreamer << "AnalyseTiltedRiemanFit"
404     << "EventNumber="           << fgEventNumber
405     << "CandidateNumber=" << fgCandidateNumber
406     << "Radius="                                        << radius
407     << "Curvature="                             << curvature
408     << "DCA="                                                   << dca
409     << "Offset="                                        << offset
410     << "Slope="                                         << slope
411     << "\n";
412   }
413 }
414
415 //____________________________________________________
416 Float_t AliTRDtrackerDebug::GetTrackRadius(Float_t a, Float_t b, Float_t c) const {
417 //
418 // Calculates the track radius using the parameters given by the tilted Rieman fit 
419 //
420 // Parameters: The three parameters from the Rieman fit
421 // Output:     The track radius
422 //
423   Float_t radius = 0;
424   if(1.0 + b*b - c*a > 0.0)
425     radius = TMath::Sqrt(1.0 + b*b - c*a )/a;
426   return radius;
427 }
428
429 //____________________________________________________
430 Float_t AliTRDtrackerDebug::GetTrackCurvature(Float_t a, Float_t b, Float_t c) const {
431 //
432 // Calculates the track curvature using the parameters given by the linear fitter 
433 //
434 // Parameters:  the three parameters from the tilted Rieman fitter
435 // Output:              the full track curvature
436 //
437   Float_t curvature =  1.0 + b*b - c*a;
438   if (curvature > 0.0) 
439     curvature  =  a / TMath::Sqrt(curvature);
440   return curvature;
441 }
442
443 //____________________________________________________
444 Float_t AliTRDtrackerDebug::GetDCA(Float_t a, Float_t b, Float_t c) const {
445 //
446 // Calculates the Distance to Clostest Approach for the Vertex using the paramters
447 // given by the tilted Rieman fit 
448 //
449 // Parameters: the three parameters from the tilted Rieman fitter
450 // Output:     the Distance to Closest Approach
451 //
452   Float_t dca  =  0.0;
453   if (1.0 + b*b - c*a > 0.0) {
454     dca = -c / (TMath::Sqrt(1.0 + b*b - c*a) + TMath::Sqrt(1.0 + b*b));
455   }
456   return dca;
457 }
458
459 //____________________________________________________
460 void AliTRDtrackerDebug::AnalyseMinMax()
461 {
462 //
463   TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
464   if(!debfile){
465     AliError("File TRD.TrackerDebug.root not found!");
466     return; 
467   }
468   fTree = (TTree *)(debfile->Get("MakeSeeds0"));
469   if(!fTree){
470     AliError("Tree MakeSeeds0 not found in File TRD.TrackerDebug.root.");
471     return;
472   }
473   AliTRDseedV1 *cseed[4] = {0x0, 0x0, 0x0, 0x0};
474   AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0};
475   for(Int_t il = 0; il < 4; il++){
476     fTree->SetBranchAddress(Form("Seed%d.", il),        &cseed[il]);
477     fTree->SetBranchAddress(Form("c%d.",il), &c[il]);
478   }
479   fTree->SetBranchAddress("CandidateNumber",    &fgCandidateNumber);
480   fTree->SetBranchAddress("EventNumber",        &fgEventNumber);
481   Int_t entries = fTree->GetEntries();
482   for(Int_t ientry = 0; ientry < entries; ientry++){
483     fTree->GetEntry(ientry);
484     Float_t minmax[2] = { -100.0,  100.0 };
485     for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
486       Float_t max = c[iLayer]->GetZ() + cseed[iLayer]->GetPadLength() * 0.5 + 1.0 - cseed[iLayer]->GetZref(0);
487       if (max < minmax[1]) minmax[1] = max;
488       Float_t min = c[iLayer]->GetZ()-cseed[iLayer]->GetPadLength() * 0.5 - 1.0 - cseed[iLayer]->GetZref(0);
489       if (min > minmax[0]) minmax[0] = min;
490     }
491     TTreeSRedirector &cstreamer = *fOutputStreamer;
492     cstreamer << "AnalyseMinMaxLayer"
493     << "EventNumber="                           << fgEventNumber
494     << "CandidateNumber="               << fgCandidateNumber
495     << "Min="                                                           << minmax[0]
496     << "Max="                                                           << minmax[1]
497     << "\n";
498   }
499 }
500
501 //____________________________________________________
502 TCanvas* AliTRDtrackerDebug::PlotSeedingConfiguration(const Char_t *direction, Int_t event, Int_t candidate){
503 //
504 // Plots the four seeding clusters, the helix fit and the reference Points for
505 // a given combination consisting of eventnr. and candidatenr.
506 //
507 // Parameters:  -direction (y or z)
508 //                              -Event Nr
509 //              -Candidate that has to be plotted
510 //
511   const Float_t kxmin = 280;
512   const Float_t kxmax = 380;
513   const Float_t kxdelta = (kxmax - kxmin)/1000;
514   
515   if((strcmp(direction, "y") != 0) && (strcmp(direction, "z") != 0)){
516     AliError(Form("Direction %s does not exist. Abborting!", direction));
517     return 0x0;
518   }
519
520   TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
521   if(!debfile){
522     AliError("File TRD.TrackerDebug.root not found!");
523     return 0x0; 
524   }
525   fTree = (TTree *)(debfile->Get("MakeSeeds0"));
526   if(!fTree){
527     AliError("Tree MakeSeeds0 not found in File TRD.TrackerDebug.root.");
528     return 0x0;
529   }
530   
531   TGraph *seedcl = new TGraph(4);
532   TGraph *seedRef = new TGraph(4);
533   TGraph *riemanFit = new TGraph(1000);
534   seedcl->SetMarkerStyle(20);
535   seedcl->SetMarkerColor(kRed);
536   seedRef->SetMarkerStyle(2);
537
538   AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0};
539   AliRieman *rim = 0x0;
540   Bool_t found = kFALSE;
541   for(Int_t il = 0; il < 4; il++) fTree->SetBranchAddress(Form("c%d.",il), &c[il]);
542   fTree->SetBranchAddress("EventNumber", &fgEventNumber);
543   fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
544   fTree->SetBranchAddress("RiemanFitter.", &rim);
545   Int_t entries = fTree->GetEntries();
546   for(Int_t entry = 0; entry < entries; entry++){
547     fTree->GetEntry(entry);
548     if(fgEventNumber < event) continue;
549     if(fgEventNumber > event) break;
550     // EventNumber matching: Do the same for the candidate number
551     if(fgCandidateNumber < candidate) continue;
552     if(fgCandidateNumber > candidate) break;
553     found = kTRUE;
554     Int_t nPoints = 0;
555     for(Int_t il = 0; il < 4; il++){
556       Float_t cluster = 0.0, reference = 0.0;
557       if(!strcmp(direction, "y")){
558         cluster = c[il]->GetY();
559         reference = rim->GetYat(c[il]->GetX());
560       }
561       else{
562         cluster = c[il]->GetZ();
563         reference = rim->GetZat(c[il]->GetX());
564       }
565       seedcl->SetPoint(nPoints, cluster, c[il]->GetX());
566       seedRef->SetPoint(nPoints, reference , c[il]->GetX());
567       nPoints++;
568     }
569     // evaluate the fitter Function numerically
570     nPoints = 0;
571     for(Int_t ipt = 0; ipt < 1000; ipt++){
572       Float_t x = kxmin + ipt * kxdelta;
573       Float_t point = 0.0;
574       if(!strcmp(direction, "y"))
575         point = rim->GetYat(x);
576       else
577         point = rim->GetZat(x);
578       riemanFit->SetPoint(nPoints++, point, x);
579     }
580     // We reached the End: break
581     break;
582   }
583   if(found){
584     seedcl->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
585     seedRef->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
586     riemanFit->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
587     TCanvas *c1 = new TCanvas();
588     seedcl->Draw("ap");
589     seedRef->Draw("psame");
590     riemanFit->Draw("lpsame");
591     return c1;
592   }
593   else{
594     AliError(Form("Combination consisting of event %d and candidate %d not found", event, candidate));
595     delete seedcl;
596     delete seedRef;
597     delete riemanFit;
598     return 0x0;
599   }
600 }
601
602 //____________________________________________________
603 TCanvas* AliTRDtrackerDebug::PlotFullTrackFit(Int_t event, Int_t candidate, Int_t iteration, const Char_t *direction){
604 //
605 // Plots the tracklets (clusters and reference in y direction) and the fitted function for several iterations
606 // in the function ImproveSeedQuality (default is before ImproveSeedQuality)
607 // 
608 // Parameters: -Event Number
609 //             -Candidate Number
610 //             -Iteration Number in ImproveSeedQuality (default: -1 = before ImproveSeedQuality)
611 //                         -direction (default: y)
612 // Output:     -TCanvas (containing the Picture);
613 //
614   const Float_t kxmin = 280;
615   const Float_t kxmax = 380;
616   const Float_t kxdelta = (kxmax - kxmin)/1000;
617   
618   if(strcmp(direction, "y") && strcmp(direction, "z")){
619     AliError(Form("Direction %s does not exist. Abborting!", direction));
620     return 0x0;
621   }
622
623   TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
624   if(!debfile){
625     AliError("File TRD.TrackerDebug.root not found.");
626     return 0x0;
627   }
628   TString *treename = 0x0;
629   if(iteration > -1)
630     treename = new TString("ImproveSeedQuality");
631   else
632     treename = new TString("MakeSeeds1");
633   fTree = (TTree *)(debfile->Get(treename->Data()));
634   if(!fTree){
635     AliError(Form("Tree %s not found in File TRD.TrackerDebug.root.", treename->Data()));
636     delete treename;
637     return 0x0;
638   }
639   delete treename;
640
641   TGraph *fitfun = new TGraph(1000);
642   // Prepare containers
643   Float_t x0[AliTRDtrackerV1::kNPlanes],
644       refP[AliTRDtrackerV1::kNPlanes],
645       clx[AliTRDtrackerV1::kNPlanes * AliTRDtrackerV1::kNTimeBins],
646       clp[AliTRDtrackerV1::kNPlanes * AliTRDtrackerV1::kNTimeBins];
647   Int_t nLayers = 0, ncls = 0;
648   
649   TLinearFitter *fitter = 0x0;
650   AliTRDseedV1 *tracklet[6] = {0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
651   for(Int_t iLayer = 0; iLayer < 6; iLayer++)
652     fTree->SetBranchAddress(Form("S%d.", iLayer), &tracklet[iLayer]);
653   fTree->SetBranchAddress("FitterT.", &fitter);
654   fTree->SetBranchAddress("EventNumber", &fgEventNumber);
655   fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
656   
657   Int_t nEntries = fTree->GetEntriesFast();
658   Bool_t found = kFALSE;
659   for(Int_t entry = 0; entry < nEntries; entry++){
660     fTree->GetEntry(entry);
661     if(fgEventNumber < event) continue;
662     if(fgEventNumber > event) break;
663     // EventNumber matching: Do the same for the candidate number
664     if(fgCandidateNumber < candidate) continue;
665     if(fgCandidateNumber > candidate) break;
666     found = kTRUE;
667     
668     for(Int_t iLayer = 0; iLayer < 6; iLayer++){
669       if(!tracklet[iLayer]->IsOK()) continue;
670       x0[nLayers] = tracklet[iLayer]->GetX0();
671       if(!strcmp(direction, "y"))
672         refP[nLayers] = tracklet[iLayer]->GetYref(0);
673       else
674         refP[nLayers] = tracklet[iLayer]->GetZref(0);
675       nLayers++;
676       for(Int_t itb = 0; itb < 30; itb++){
677         if(!tracklet[iLayer]->IsUsable(itb)) continue;
678         AliTRDcluster *cl = tracklet[iLayer]->GetClusters(itb);
679         if(!strcmp(direction, "y"))
680           clp[ncls] = cl->GetY();
681         else
682           clp[ncls] = cl->GetZ();
683         clx[ncls] = cl->GetX();
684         ncls++;
685       }
686     }
687     // Add function derived by the tilted Rieman fit (Defined by the curvature)
688     Int_t nPoints = 0;
689     if(!strcmp(direction, "y")){
690       Double_t a = fitter->GetParameter(0);
691       Double_t b = fitter->GetParameter(1);
692       Double_t c = fitter->GetParameter(2);
693       Double_t curvature =  1.0 + b*b - c*a;
694       if (curvature > 0.0) {
695         curvature  =  a / TMath::Sqrt(curvature);
696       }
697       // Numerical evaluation of the function:
698       for(Int_t ipt = 0; ipt < 1000; ipt++){
699         Float_t x = kxmin + ipt * kxdelta;
700         Double_t res = (x * a + b);                                                             // = (x - x0)/y0
701         res *= res;
702         res  = 1.0 - c * a + b * b - res;                                       // = (R^2 - (x - x0)^2)/y0^2
703         Double_t y = 0.;
704         if (res >= 0) {
705           res = TMath::Sqrt(res);
706           y    = (1.0 - res) / a;
707         }
708         fitfun->SetPoint(nPoints++, y, x);
709       }
710     }
711     else{
712       Double_t offset   = fitter->GetParameter(3);
713       Double_t slope    = fitter->GetParameter(4);       
714       // calculate the reference x (defined as medium between layer 2 and layer 3)
715       // same procedure as in the tracker code
716       Float_t medx = 0, xref = 0;
717       Int_t startIndex = 5, nDistances = 0;
718       for(Int_t iLayer = 5; iLayer > 0; iLayer--){
719         if(tracklet[iLayer]->IsOK() && tracklet[iLayer - 1]->IsOK()){
720           medx += tracklet[iLayer]->GetX0() - tracklet[iLayer - 1]->GetX0();
721           startIndex = iLayer - 1;
722           nDistances++;
723         }
724       }
725       if(nDistances){
726         medx /= nDistances;
727       }
728       else{
729         Float_t xpos[2];        memset(xpos, 0, sizeof(Float_t) * 2);
730         Int_t ien = 0, idiff = 0;
731         for(Int_t iLayer = 5; iLayer > 0; iLayer--){
732           if(tracklet[iLayer]->IsOK()){
733             xpos[ien++] = tracklet[iLayer]->GetX0();
734             startIndex = iLayer;
735           }
736           if(ien)
737             idiff++;
738           if(ien >=2)
739             break;
740         }
741         medx = (xpos[0] - xpos[1])/idiff;
742       }
743       xref = tracklet[startIndex]->GetX0() + medx * (2.5 - startIndex) - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
744
745       for(Int_t ipt = 0; ipt < 1000; ipt++){
746         Float_t x = kxmin + ipt * kxdelta;
747         Float_t z = offset + slope * (x - xref);
748         fitfun->SetPoint(nPoints++, z, x);
749       }
750     }
751     break;
752   }
753   if(found){
754     TGraph *trGraph             = new TGraph(ncls);
755     TGraph *refPoints   = new TGraph(nLayers);
756     trGraph->SetMarkerStyle(20);
757     trGraph->SetMarkerColor(kRed);
758     refPoints->SetMarkerStyle(21);
759     refPoints->SetMarkerColor(kBlue);
760     // fill the graphs
761     for(Int_t iLayer = 0; iLayer < nLayers; iLayer++)
762       refPoints->SetPoint(iLayer, refP[iLayer], x0[iLayer]);
763     for(Int_t icls = 0; icls < ncls; icls++)
764       trGraph->SetPoint(icls, clp[icls], clx[icls]);
765     TCanvas *c1 = new TCanvas();
766     trGraph->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
767     refPoints->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
768     fitfun->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
769     trGraph->Draw("ap");
770     refPoints->Draw("psame");
771     fitfun->Draw("lpsame");
772     return c1;
773   }
774   else{
775     AliError(Form("Combination consisting of event %d and candidate %d not found", event, candidate));
776     delete fitfun;
777     return 0x0;
778   }
779 }
780