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