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