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