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