]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDtrackerDebug.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fix coverity
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackerDebug.cxx
CommitLineData
eb38ed55 1/**************************************************************************
e3cf3d02 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**************************************************************************/
eb38ed55 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
eb38ed55 28#include "TFile.h"
29#include "TTree.h"
30#include "TTreeStream.h"
31#include "TLinearFitter.h"
bb56afff 32#include "TGraph.h"
33#include "TCanvas.h"
34#include "TMath.h"
eb38ed55 35
36#include "AliLog.h"
4ecadb52 37#include "AliRieman.h"
4d6aee34 38
eb38ed55 39#include "AliTRDgeometry.h"
40#include "AliTRDtrackV1.h"
41#include "AliTRDseedV1.h"
eb38ed55 42#include "AliTRDcluster.h"
bb56afff 43#include "AliTRDgeometry.h"
eb38ed55 44
4d6aee34 45#include "AliTRDtrackerDebug.h"
46
eb38ed55 47ClassImp(AliTRDtrackerDebug)
48
bb56afff 49Int_t AliTRDtrackerDebug::fgEventNumber = 0;
50Int_t AliTRDtrackerDebug::fgTrackNumber = 0;
51Int_t AliTRDtrackerDebug::fgCandidateNumber = 0;
52
eb38ed55 53//____________________________________________________
54AliTRDtrackerDebug::AliTRDtrackerDebug() : AliTRDtrackerV1()
4d6aee34 55 ,fOutputStreamer(NULL)
56 ,fTree(NULL)
57 ,fTracklet(NULL)
58 ,fTrack(NULL)
e3cf3d02 59 ,fNClusters(0)
60 ,fAlpha(0.)
eb38ed55 61{
62 //
e3cf3d02 63 // Default constructor
64 //
65 fOutputStreamer = new TTreeSRedirector("TRD.Debug.root");
eb38ed55 66}
67
68//____________________________________________________
69AliTRDtrackerDebug::~AliTRDtrackerDebug()
70{
e3cf3d02 71 // destructor
72
73 delete fOutputStreamer;
eb38ed55 74}
75
76
eb38ed55 77//____________________________________________________
78Bool_t AliTRDtrackerDebug::Init()
79{
80// steer linking data for various debug streams
e3cf3d02 81 fTrack = new AliTRDtrackV1();
82 fTree->SetBranchAddress("ncl", &fNClusters);
83 fTree->SetBranchAddress("track.", &fTrack);
84 return kTRUE;
eb38ed55 85}
86
87//____________________________________________________
88Bool_t AliTRDtrackerDebug::Open(const char *method)
89{
e3cf3d02 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;
eb38ed55 102}
103
104//____________________________________________________
105Int_t AliTRDtrackerDebug::Process()
106{
107// steer debug process threads
e3cf3d02 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
4d6aee34 116 const AliTRDseedV1 *tracklet = NULL;
e3cf3d02 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;
eb38ed55 127}
128
129
130//____________________________________________________
131void AliTRDtrackerDebug::ResidualsClustersTrack(const AliTRDseedV1 *tracklet)
132{
133// Calculate averange distances from clusters to the TRD track
e3cf3d02 134
135 Double_t x[3];
4d6aee34 136 AliTRDcluster *c = NULL;
e3cf3d02 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
41f024fb 142 if(!PropagateToX(*fTrack, xc, 2.)) continue;
e3cf3d02 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 }
eb38ed55 160}
161
162//____________________________________________________
163void AliTRDtrackerDebug::ResidualsClustersTracklet(const AliTRDseedV1 *tracklet) const
164{
165// Calculates distances from clusters to tracklets
e3cf3d02 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
4d6aee34 173 AliTRDcluster *c = NULL;
e3cf3d02 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 :
eb38ed55 180 //ResidualsClustersTracklet->Draw("TMath::Abs(10.*dy):TMath::ATan(ys)*TMath::RadToDeg()>>h(20, -40, 40)", "", "prof");
e3cf3d02 181 TTreeSRedirector &cstreamer = *fOutputStreamer;
182 cstreamer << "ResidualsClustersTracklet"
183 << "c.=" << c
184 << "ys=" << ys
185 << "dy=" << dy
186 << "\n";
187 }
eb38ed55 188}
189
190//____________________________________________________
191void AliTRDtrackerDebug::ResidualsClustersParametrisation(const AliTRDseedV1 *tracklet) const
192{
193// Calculates distances from clusters to Rieman fit.
e3cf3d02 194
195 // store cluster positions
196 Double_t x0 = tracklet->GetX0();
4d6aee34 197 AliTRDcluster *c = NULL;
e3cf3d02 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 }
eb38ed55 227}
228
229
230//____________________________________________________
231void AliTRDtrackerDebug::ResidualsTrackletsTrack() const
232{
233// Calculates distances from tracklets to the TRD track.
e3cf3d02 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
4d6aee34 239 //AliTRDseedV1 tracklet[6] = {NULL, NULL, NULL, NULL, NULL, NULL};
240 AliTRDseedV1 tracklet[6];
241 const AliTRDseedV1 *ctracklet = NULL;
e3cf3d02 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
4d6aee34 265 AliTRDtrackerV1::FitRiemanTilt(NULL, &tracklet[0], kTRUE);
e3cf3d02 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
edda7a9d 275 AliTRDseedV1 *ptr(NULL);
276 if(!(ptr = fTrack->GetTracklet(ip))) continue;
277 tracklet[ip] = (*ptr);
e3cf3d02 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);
f29f13a6 290 ts.Fit(kTRUE);
e3cf3d02 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 }
eb38ed55 313}
314
bb56afff 315//____________________________________________________
316void 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//
e3cf3d02 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++)
4d6aee34 340 tracklets[iPlane] = NULL;
e3cf3d02 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 }
bb56afff 369}
370//____________________________________________________
371void 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//
e3cf3d02 381 TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
382 fTree = (TTree *)(debfile->Get("MakeSeeds2"));
383 if(!fTree) return;
384 Int_t nEntries = fTree->GetEntries();
4d6aee34 385 TLinearFitter *tiltedRiemanFitter = NULL;
e3cf3d02 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 }
bb56afff 410}
411
412//____________________________________________________
413Float_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//
e3cf3d02 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;
bb56afff 424}
425
426//____________________________________________________
427Float_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//
e3cf3d02 434 Float_t curvature = 1.0 + b*b - c*a;
435 if (curvature > 0.0)
436 curvature = a / TMath::Sqrt(curvature);
437 return curvature;
bb56afff 438}
439
440//____________________________________________________
441Float_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//
e3cf3d02 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;
bb56afff 454}
455
456//____________________________________________________
457void AliTRDtrackerDebug::AnalyseMinMax()
458{
4d6aee34 459 // Development function related to the old tracking code
e3cf3d02 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 }
4d6aee34 470 AliTRDseedV1 *cseed[4] = {NULL, NULL, NULL, NULL};
471 AliTRDcluster *c[4] = {NULL, NULL, NULL, NULL};
e3cf3d02 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 }
bb56afff 496}
497
498//____________________________________________________
c202e52e 499TCanvas* AliTRDtrackerDebug::PlotSeedingConfiguration(const Char_t *direction, Int_t event, Int_t candidate){
bb56afff 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//
e3cf3d02 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));
4d6aee34 514 return NULL;
e3cf3d02 515 }
516
517 TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
518 if(!debfile){
519 AliError("File TRD.TrackerDebug.root not found!");
4d6aee34 520 return NULL;
e3cf3d02 521 }
522 fTree = (TTree *)(debfile->Get("MakeSeeds0"));
523 if(!fTree){
524 AliError("Tree MakeSeeds0 not found in File TRD.TrackerDebug.root.");
4d6aee34 525 return NULL;
e3cf3d02 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
4d6aee34 535 AliTRDcluster *c[4] = {NULL, NULL, NULL, NULL};
536 AliRieman *rim = NULL;
e3cf3d02 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;
4d6aee34 595 return NULL;
e3cf3d02 596 }
bb56afff 597}
598
599//____________________________________________________
c202e52e 600TCanvas* AliTRDtrackerDebug::PlotFullTrackFit(Int_t event, Int_t candidate, Int_t iteration, const Char_t *direction){
bb56afff 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//
e3cf3d02 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));
4d6aee34 617 return NULL;
e3cf3d02 618 }
619
620 TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
621 if(!debfile){
622 AliError("File TRD.TrackerDebug.root not found.");
4d6aee34 623 return NULL;
e3cf3d02 624 }
4d6aee34 625 TString *treename = NULL;
e3cf3d02 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;
4d6aee34 634 return NULL;
e3cf3d02 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
4d6aee34 646 TLinearFitter *fitter = NULL;
647 AliTRDseedV1 *tracklet[6] = {NULL, NULL, NULL, NULL, NULL, NULL};
e3cf3d02 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++;
41f024fb 673 AliTRDcluster *cl(NULL);
e3cf3d02 674 for(Int_t itb = 0; itb < 30; itb++){
675 if(!tracklet[iLayer]->IsUsable(itb)) continue;
41f024fb 676 if(!(cl = tracklet[iLayer]->GetClusters(itb))) continue;
677
e3cf3d02 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;
4d6aee34 776 return NULL;
e3cf3d02 777 }
bb56afff 778}
779