TRD/AliTRDtrackerDebug.cxx

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