Using SetMaxDCAToVertexXY instead of the removed SetDCAToVertex
[u/mrichter/AliRoot.git] / TRD / AliTRDtrackerDebug.cxx
CommitLineData
eb38ed55 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"
bb56afff 34#include "TGraph.h"
35#include "TCanvas.h"
36#include "TMath.h"
eb38ed55 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"
bb56afff 44#include "AliTRDgeometry.h"
eb38ed55 45
46ClassImp(AliTRDtrackerDebug)
47
bb56afff 48Int_t AliTRDtrackerDebug::fgEventNumber = 0;
49Int_t AliTRDtrackerDebug::fgTrackNumber = 0;
50Int_t AliTRDtrackerDebug::fgCandidateNumber = 0;
51
eb38ed55 52//____________________________________________________
53AliTRDtrackerDebug::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//____________________________________________________
68AliTRDtrackerDebug::~AliTRDtrackerDebug()
69{
70 // destructor
71
72 delete fOutputStreamer;
73}
74
75
76//____________________________________________________
c202e52e 77void AliTRDtrackerDebug::Draw(Option_t *)
eb38ed55 78{
79// steer draw function
80}
81
82
83//____________________________________________________
84Bool_t AliTRDtrackerDebug::Init()
85{
86// steer linking data for various debug streams
eb38ed55 87 fTrack = new AliTRDtrackV1();
88 fTree->SetBranchAddress("ncl", &fNClusters);
89 fTree->SetBranchAddress("track.", &fTrack);
90 return kTRUE;
91}
92
93//____________________________________________________
94Bool_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//____________________________________________________
111Int_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//____________________________________________________
137void 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//____________________________________________________
169void 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//____________________________________________________
197void 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//____________________________________________________
237void 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
bb56afff 319//____________________________________________________
320void 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//____________________________________________________
375void 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//____________________________________________________
417Float_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//____________________________________________________
431Float_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//____________________________________________________
445Float_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//____________________________________________________
461void 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//____________________________________________________
c202e52e 503TCanvas* AliTRDtrackerDebug::PlotSeedingConfiguration(const Char_t *direction, Int_t event, Int_t candidate){
bb56afff 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//____________________________________________________
c202e52e 604TCanvas* AliTRDtrackerDebug::PlotFullTrackFit(Int_t event, Int_t candidate, Int_t iteration, const Char_t *direction){
bb56afff 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