TRD/AliTRDtrackerDebug.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliTRDtrackerDebug.cxx 23810 2008-02-08 09:00:27Z hristov $ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Tracker debug streamer                                                   //
  21 //                                                                           //
  22 //  Authors:                                                                 //
  23 //    Alex Bercuci <A.Bercuci@gsi.de>                                        //
  24 //    Markus Fasel <M.Fasel@gsi.de>                                          //
  25 //                                                                           //
  26 ///////////////////////////////////////////////////////////////////////////////
  27
  28 #include "AliTRDtrackerDebug.h"
  29
  30 #include "TFile.h"
  31 #include "TTree.h"
  32 #include "TTreeStream.h"
  33 #include "TLinearFitter.h"
  34 #include "TGraph.h"
  35 #include "TCanvas.h"
  36 #include "TMath.h"
  37
  38 #include "AliLog.h"
  39 #include "AliTRDgeometry.h"
  40 #include "AliTRDtrackV1.h"
  41 #include "AliTRDseedV1.h"
  42 #include "AliTRDseed.h"
  43 #include "AliTRDcluster.h"
  44 #include "AliTRDgeometry.h"
  45
  46 ClassImp(AliTRDtrackerDebug)
  47
  48 Int_t AliTRDtrackerDebug::fgEventNumber = 0;
  49 Int_t AliTRDtrackerDebug::fgTrackNumber = 0;
  50 Int_t AliTRDtrackerDebug::fgCandidateNumber = 0;
  51
  52 //____________________________________________________
  53 AliTRDtrackerDebug::AliTRDtrackerDebug() : AliTRDtrackerV1()
  54         ,fOutputStreamer(0x0)
  55         ,fTree(0x0)
  56         ,fTracklet(0x0)
  57         ,fTrack(0x0)
  58         ,fNClusters(0)
  59         ,fAlpha(0.)
  60 {
  61         //
  62         // Default constructor
  63         //
  64         fOutputStreamer = new TTreeSRedirector("TRD.Debug.root");
  65 }
  66
  67 //____________________________________________________
  68 AliTRDtrackerDebug::~AliTRDtrackerDebug()
  69 {
  70         // destructor
  71
  72         delete fOutputStreamer;
  73 }
  74
  75
  76 //____________________________________________________
  77 void AliTRDtrackerDebug::Draw(const Option_t *)
  78 {
  79 // steer draw function
  80 }
  81
  82
  83 //____________________________________________________
  84 Bool_t AliTRDtrackerDebug::Init()
  85 {
  86 // steer linking data for various debug streams
  87         fTrack = new AliTRDtrackV1();
  88         fTree->SetBranchAddress("ncl", &fNClusters);
  89         fTree->SetBranchAddress("track.", &fTrack);
  90         return kTRUE;
  91 }
  92
  93 //____________________________________________________
  94 Bool_t AliTRDtrackerDebug::Open(const char *method)
  95 {
  96         // Connect to the tracker debug file
  97
  98         TDirectory *savedir = gDirectory;
  99         TFile::Open("TRD.TrackerDebugger.root");
 100         fTree = (TTree*)gFile->Get(method);
 101         if(!fTree){
 102                 AliInfo(Form("Can not find debug stream for the %s method.\n", method));
 103                 savedir->cd();
 104                 return kFALSE;
 105         }
 106         savedir->cd();
 107         return kTRUE;
 108 }
 109
 110 //____________________________________________________
 111 Int_t AliTRDtrackerDebug::Process()
 112 {
 113 // steer debug process threads
 114
 115         for(int it = 0; it<fTree->GetEntries(); it++){
 116                 if(!fTree->GetEntry(it)) continue;
 117                 if(!fNClusters) continue;
 118                 fAlpha = fTrack->GetAlpha();
 119                 //printf("Processing track %d [%d] ...\n", it, fNClusters);
 120                 ResidualsTrackletsTrack();
 121
 122                 const AliTRDseedV1 *tracklet = 0x0;
 123                 for(int ip = 5; ip>=0; ip--){
 124                         if(!(tracklet = fTrack->GetTracklet(ip))) continue;
 125                         if(!tracklet->GetN()) continue;
 126
 127                         ResidualsClustersTrack(tracklet);
 128                         ResidualsClustersTracklet(tracklet);
 129                         ResidualsClustersParametrisation(tracklet);
 130                 }
 131         }
 132         return kTRUE;
 133 }
 134
 135
 136 //____________________________________________________
 137 void AliTRDtrackerDebug::ResidualsClustersTrack(const AliTRDseedV1 *tracklet)
 138 {
 139 // Calculate averange distances from clusters to the TRD track
 140
 141         Double_t x[3];
 142         AliTRDcluster *c = 0x0;
 143         for(int ic=0; ic<35/*AliTRDseed:knTimebins*/; ic++){
 144                 if(!(c = tracklet->GetClusters(ic))) continue;
 145                 Double_t xc = c->GetX(), yc = c->GetY(), zc = c->GetZ();
 146
 147                 // propagate track to cluster
 148                 PropagateToX(*fTrack, xc, 2.);
 149                 fTrack->GetXYZ(x);
 150
 151                 // transform to local tracking coordinates
 152                 //Double_t xg =  x[0] * TMath::Cos(fAlpha) + x[1] * TMath::Sin(fAlpha);
 153                 Double_t yg = -x[0] * TMath::Sin(fAlpha) + x[1] * TMath::Cos(fAlpha);
 154
 155                 // apply tilt pad correction
 156                 yc+= (zc - x[2]) * tracklet->GetTilt();
 157
 158                 Double_t dy = yc-yg;
 159
 160                 TTreeSRedirector &cstreamer = *fOutputStreamer;
 161                 cstreamer << "ResidualsClustersTrack"
 162                         << "c.="   << c
 163                         << "dy="   << dy
 164                         << "\n";
 165         }
 166 }
 167
 168 //____________________________________________________
 169 void AliTRDtrackerDebug::ResidualsClustersTracklet(const AliTRDseedV1 *tracklet) const
 170 {
 171 // Calculates distances from clusters to tracklets
 172
 173         Double_t x0 = tracklet->GetX0(),
 174                  y0 = tracklet->GetYfit(0),
 175                  ys = tracklet->GetYfit(1);
 176                  //z0 = tracklet->GetZfit(0),
 177                  //zs = tracklet->GetZfit(1);
 178
 179         AliTRDcluster *c = 0x0;
 180         for(int ic=0; ic<35/*AliTRDseed:knTimebins*/; ic++){
 181                 if(!(c = tracklet->GetClusters(ic))) continue;
 182                 Double_t xc = c->GetX(), yc = c->GetY()/*, zc = c->GetZ()*/;
 183                 Double_t dy = yc- (y0-(x0-xc)*ys);
 184
 185                 //To draw  use :
 186     //ResidualsClustersTracklet->Draw("TMath::Abs(10.*dy):TMath::ATan(ys)*TMath::RadToDeg()>>h(20, -40, 40)", "", "prof");
 187                 TTreeSRedirector &cstreamer = *fOutputStreamer;
 188                 cstreamer << "ResidualsClustersTracklet"
 189                         << "c.="   << c
 190                         << "ys="   << ys
 191                         << "dy="   << dy
 192                         << "\n";
 193         }
 194 }
 195
 196 //____________________________________________________
 197 void AliTRDtrackerDebug::ResidualsClustersParametrisation(const AliTRDseedV1 *tracklet) const
 198 {
 199 // Calculates distances from clusters to Rieman fit.
 200
 201         // store cluster positions
 202         Double_t x0 = tracklet->GetX0();
 203         AliTRDcluster *c = 0x0;
 204
 205         Double_t x[2]; Int_t ncl, mcl, jc;
 206         TLinearFitter fitter(3, "hyp2");
 207         for(int ic=0; ic<35/*AliTRDseed:knTimebins*/; ic++){
 208                 if(!(c = tracklet->GetClusters(ic))) continue;
 209                 Double_t xc = c->GetX(), yc = c->GetY()/*, zc = c->GetZ()*/;
 210
 211                 jc = ic; ncl = 0; mcl=0; fitter.ClearPoints();
 212                 while(ncl<6){
 213                         // update index
 214                         mcl++;
 215                         jc = ic + ((mcl&1)?-1:1)*(mcl>>1);
 216
 217                         if(jc<0 || jc>=35) continue;
 218                         if(!(c = tracklet->GetClusters(jc))) continue;
 219
 220                         x[0] = c->GetX()-x0;
 221                         x[1] = x[0]*x[0];
 222                         fitter.AddPoint(x, c->GetY(), c->GetSigmaY2());
 223                         ncl++;
 224                 }
 225                 fitter.Eval();
 226                 Double_t dy = yc - fitter.GetParameter(0) -fitter.GetParameter(1) * (xc-x0) - fitter.GetParameter(2)* (xc-x0)*(xc-x0);
 227
 228                 TTreeSRedirector &cstreamer = *fOutputStreamer;
 229                 cstreamer << "ResidualsClustersParametrisation"
 230                         << "dy="   << dy
 231                         << "\n";
 232         }
 233 }
 234
 235
 236 //____________________________________________________
 237 void AliTRDtrackerDebug::ResidualsTrackletsTrack() const
 238 {
 239 // Calculates distances from tracklets to the TRD track.
 240
 241         if(fTrack->GetNumberOfTracklets() < 6) return;
 242
 243         // build a working copy of the tracklets attached to the track
 244         // and initialize working variables fX, fY and fZ
 245         AliTRDcluster *c = 0x0;
 246         AliTRDseedV1 tracklet[6] = {0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
 247         const AliTRDseedV1 *ctracklet = 0x0;
 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<35/*AliTRDseed: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                 AliTRDseedV1::FitRiemanTilt(&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<35/*AliTRDseed: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.Update();
 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] = 0x0;
 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 = 0x0;
 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 //
 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] = {0x0, 0x0, 0x0, 0x0};
 475         AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0};
 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(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 0x0;
 519         }
 520
 521         TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
 522         if(!debfile){
 523                 AliError("File TRD.TrackerDebug.root not found!");
 524                 return 0x0;
 525         }
 526         fTree = (TTree *)(debfile->Get("MakeSeeds0"));
 527         if(!fTree){
 528                 AliError("Tree MakeSeeds0 not found in File TRD.TrackerDebug.root.");
 529                 return 0x0;
 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] = {0x0, 0x0, 0x0, 0x0};
 540         AliRieman *rim = 0x0;
 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 0x0;
 600         }
 601 }
 602
 603 //____________________________________________________
 604 TCanvas* AliTRDtrackerDebug::PlotFullTrackFit(Int_t event, Int_t candidate, Int_t iteration, 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 0x0;
 622         }
 623
 624         TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
 625         if(!debfile){
 626                 AliError("File TRD.TrackerDebug.root not found.");
 627                 return 0x0;
 628         }
 629         TString *treename = 0x0;
 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 0x0;
 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 = 0x0;
 651         AliTRDseedV1 *tracklet[6] = {0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
 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 0x0;
 779         }
 780 }
 781