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