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