1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliTRDtrackerDebug.cxx 23810 2008-02-08 09:00:27Z hristov $ */
18 ///////////////////////////////////////////////////////////////////////////////
20 // Tracker debug streamer //
23 // Alex Bercuci <A.Bercuci@gsi.de> //
24 // Markus Fasel <M.Fasel@gsi.de> //
26 ///////////////////////////////////////////////////////////////////////////////
28 #include "AliTRDtrackerDebug.h"
32 #include "TTreeStream.h"
33 #include "TLinearFitter.h"
39 #include "AliTRDgeometry.h"
40 #include "AliTRDtrackV1.h"
41 #include "AliTRDseedV1.h"
42 #include "AliTRDseed.h"
43 #include "AliTRDcluster.h"
44 #include "AliTRDgeometry.h"
46 ClassImp(AliTRDtrackerDebug)
48 Int_t AliTRDtrackerDebug::fgEventNumber = 0;
49 Int_t AliTRDtrackerDebug::fgTrackNumber = 0;
50 Int_t AliTRDtrackerDebug::fgCandidateNumber = 0;
52 //____________________________________________________
53 AliTRDtrackerDebug::AliTRDtrackerDebug() : AliTRDtrackerV1()
62 // Default constructor
64 fOutputStreamer = new TTreeSRedirector("TRD.Debug.root");
67 //____________________________________________________
68 AliTRDtrackerDebug::~AliTRDtrackerDebug()
72 delete fOutputStreamer;
76 //____________________________________________________
77 void AliTRDtrackerDebug::Draw(Option_t *)
79 // steer draw function
83 //____________________________________________________
84 Bool_t AliTRDtrackerDebug::Init()
86 // steer linking data for various debug streams
87 fTrack = new AliTRDtrackV1();
88 fTree->SetBranchAddress("ncl", &fNClusters);
89 fTree->SetBranchAddress("track.", &fTrack);
93 //____________________________________________________
94 Bool_t AliTRDtrackerDebug::Open(const char *method)
96 // Connect to the tracker debug file
98 TDirectory *savedir = gDirectory;
99 TFile::Open("TRD.TrackerDebugger.root");
100 fTree = (TTree*)gFile->Get(method);
102 AliInfo(Form("Can not find debug stream for the %s method.\n", method));
110 //____________________________________________________
111 Int_t AliTRDtrackerDebug::Process()
113 // steer debug process threads
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();
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;
127 ResidualsClustersTrack(tracklet);
128 ResidualsClustersTracklet(tracklet);
129 ResidualsClustersParametrisation(tracklet);
136 //____________________________________________________
137 void AliTRDtrackerDebug::ResidualsClustersTrack(const AliTRDseedV1 *tracklet)
139 // Calculate averange distances from clusters to the TRD track
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();
147 // propagate track to cluster
148 PropagateToX(*fTrack, xc, 2.);
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);
155 // apply tilt pad correction
156 yc+= (zc - x[2]) * tracklet->GetTilt();
160 TTreeSRedirector &cstreamer = *fOutputStreamer;
161 cstreamer << "ResidualsClustersTrack"
168 //____________________________________________________
169 void AliTRDtrackerDebug::ResidualsClustersTracklet(const AliTRDseedV1 *tracklet) const
171 // Calculates distances from clusters to tracklets
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);
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);
186 //ResidualsClustersTracklet->Draw("TMath::Abs(10.*dy):TMath::ATan(ys)*TMath::RadToDeg()>>h(20, -40, 40)", "", "prof");
187 TTreeSRedirector &cstreamer = *fOutputStreamer;
188 cstreamer << "ResidualsClustersTracklet"
196 //____________________________________________________
197 void AliTRDtrackerDebug::ResidualsClustersParametrisation(const AliTRDseedV1 *tracklet) const
199 // Calculates distances from clusters to Rieman fit.
201 // store cluster positions
202 Double_t x0 = tracklet->GetX0();
203 AliTRDcluster *c = 0x0;
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()*/;
211 jc = ic; ncl = 0; mcl=0; fitter.ClearPoints();
215 jc = ic + ((mcl&1)?-1:1)*(mcl>>1);
217 if(jc<0 || jc>=35) continue;
218 if(!(c = tracklet->GetClusters(jc))) continue;
222 fitter.AddPoint(x, c->GetY(), c->GetSigmaY2());
226 Double_t dy = yc - fitter.GetParameter(0) -fitter.GetParameter(1) * (xc-x0) - fitter.GetParameter(2)* (xc-x0)*(xc-x0);
228 TTreeSRedirector &cstreamer = *fOutputStreamer;
229 cstreamer << "ResidualsClustersParametrisation"
236 //____________________________________________________
237 void AliTRDtrackerDebug::ResidualsTrackletsTrack() const
239 // Calculates distances from tracklets to the TRD track.
241 if(fTrack->GetNumberOfTracklets() < 6) return;
243 // build a working copy of the tracklets attached to the track
244 // and initialize working variables fX, fY and fZ
245 AliTRDseedV1 tracklet[6] = {0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
246 const AliTRDseedV1 *ctracklet = 0x0;
247 for(int ip = 0; ip<6; ip++){
248 if(!(ctracklet = fTrack->GetTracklet(ip))) continue;
249 tracklet[ip] = (*ctracklet);
250 // Double_t x0 = tracklet[ip].GetX0();
251 // for(int ic=0; ic<AliTRDseedV1:knTimebins; ic++){
252 // if(!(c = tracklet[ip].GetClusters(ic))) continue;
253 // Double_t xc = c->GetX(), yc = c->GetY(), zc = c->GetZ();
254 // tracklet[ip].SetX(ic, xc-x0);
255 // tracklet[ip].SetY(ic, yc);
256 // tracklet[ip].SetZ(ic, zc);
260 // Do a Rieman fit (with tilt correction) for all tracklets
261 // except the one which is tested.
262 // (Based on AliTRDseed::IsOK() return false)
263 for(int ip=0; ip<6; ip++){
264 // reset tracklet to be tested
265 Double_t x0 = tracklet[ip].GetX0();
266 new(&tracklet[ip]) AliTRDseedV1();
267 tracklet[ip].SetX0(x0);
269 // fit Rieman with tilt correction
270 AliTRDtrackerV1::FitRiemanTilt(0x0, &tracklet[0], kTRUE);
272 // make a copy of the fit result
274 y0 = tracklet[ip].GetYref(0),
275 dydx = tracklet[ip].GetYref(1),
276 z0 = tracklet[ip].GetZref(0),
277 dzdx = tracklet[ip].GetZref(1);
280 tracklet[ip] = (*fTrack->GetTracklet(ip));
281 // for(int ic=0; ic<AliTRDseedV1:knTimebins; ic++){
282 // if(!(c = tracklet[ip].GetClusters(ic))) continue;
283 // Double_t xc = c->GetX(), yc = c->GetY(), zc = c->GetZ();
284 // tracklet[ip].SetX(ic, xc-x0);
285 // tracklet[ip].SetY(ic, yc);
286 // tracklet[ip].SetZ(ic, zc);
290 AliTRDseedV1 ts(tracklet[ip]);
291 ts.SetYref(0, y0); ts.SetYref(1, dydx);
292 ts.SetZref(0, z0); ts.SetZref(1, dzdx);
295 // save results for plotting
296 Int_t plane = tracklet[ip].GetPlane();
297 Double_t dy = tracklet[ip].GetYfit(0) - ts.GetYfit(0);
298 Double_t tgy = tracklet[ip].GetYfit(1);
299 Double_t dtgy = (tracklet[ip].GetYfit(1) - ts.GetYfit(1))/(1. + tracklet[ip].GetYfit(1) * ts.GetYfit(1));
300 Double_t dz = tracklet[ip].GetZfit(0) - ts.GetZfit(0);
301 Double_t tgz = tracklet[ip].GetZfit(1);
302 Double_t dtgz = (tracklet[ip].GetZfit(1) - ts.GetZfit(1))/(1. + tracklet[ip].GetZfit(1) * ts.GetZfit(1));
303 TTreeSRedirector &cstreamer = *fOutputStreamer;
304 cstreamer << "ResidualsTrackletsTrack"
305 << "ref.=" << &tracklet[ip]
318 //____________________________________________________
319 void AliTRDtrackerDebug::AnalyseFindable(Char_t *treename){
321 // Calculates the number of findable tracklets defined as the number of tracklets
322 // per track candidate where the tan phi_tracklet is below 0.15 (maximum inclination
325 // Parameters: -the treename (this method can be used for all trees which store the
329 // A new TTree containing the number of findable tracklets and the number of clusters
330 // attached to the full track is stored to disk
333 TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
334 fTree = (TTree *)(debfile->Get(treename));
336 AliError(Form("Tree %s not found in file TRDdebug.root. Abborting!", treename));
341 AliTRDseedV1 *tracklets[kNPlanes];
342 for(Int_t iPlane = 0; iPlane < AliTRDtrackerV1::kNPlanes; iPlane++)
343 tracklets[iPlane] = 0x0;
344 for(Int_t iPlane = 0; iPlane < kNPlanes; iPlane++)
345 fTree->SetBranchAddress(Form("S%d.", iPlane), &tracklets[iPlane]);
346 fTree->SetBranchAddress("EventNumber", &fgEventNumber);
347 fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
349 Int_t findable = 0, nClusters = 0;
350 Int_t nEntries = fTree->GetEntriesFast();
351 for(Int_t iEntry = 0; iEntry < nEntries; iEntry++){
352 printf("Entry %d\n", iEntry);
353 fTree->GetEntry(iEntry);
356 // Calculate Findable
357 for(Int_t iPlane = 0; iPlane < kNPlanes; iPlane++){
358 if (TMath::Abs(tracklets[iPlane]->GetYref(0) / tracklets[iPlane]->GetX0()) < 0.15) findable++;
359 if (!tracklets[iPlane]->IsOK()) continue;
360 nClusters += tracklets[iPlane]->GetN2();
364 TTreeSRedirector &cstreamer = *fOutputStreamer;
365 cstreamer << "AnalyseFindable"
366 << "EventNumber=" << fgEventNumber
367 << "CandidateNumber=" << fgCandidateNumber
368 << "Findable=" << findable
369 << "NClusters=" << nClusters
373 //____________________________________________________
374 void AliTRDtrackerDebug::AnalyseTiltedRiemanFit(){
376 // Creating a Data Set for the method FitTiltedRieman containing usefull variables
377 // Each variable can be addressed to tracks later. Data can be processed later.
382 // TODO: Match everything with Init and Process
384 TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
385 fTree = (TTree *)(debfile->Get("MakeSeeds2"));
387 Int_t nEntries = fTree->GetEntries();
388 TLinearFitter *tiltedRiemanFitter = 0x0;
389 fTree->SetBranchAddress("FitterT.", &tiltedRiemanFitter);
390 fTree->SetBranchAddress("EventNumber", &fgEventNumber);
391 fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
392 for(Int_t entry = 0; entry < nEntries; entry++){
393 fTree->GetEntry(entry);
394 Double_t a = tiltedRiemanFitter->GetParameter(0);
395 Double_t b = tiltedRiemanFitter->GetParameter(1);
396 Double_t c = tiltedRiemanFitter->GetParameter(2);
397 Double_t offset = tiltedRiemanFitter->GetParameter(3);
398 Double_t slope = tiltedRiemanFitter->GetParameter(4);
399 Float_t radius = GetTrackRadius(a, b, c);
400 Float_t curvature = GetTrackCurvature(a, b, c);
401 Float_t dca = GetDCA(a, b, c);
402 TTreeSRedirector &cstreamer = *fOutputStreamer;
403 cstreamer << "AnalyseTiltedRiemanFit"
404 << "EventNumber=" << fgEventNumber
405 << "CandidateNumber=" << fgCandidateNumber
406 << "Radius=" << radius
407 << "Curvature=" << curvature
409 << "Offset=" << offset
415 //____________________________________________________
416 Float_t AliTRDtrackerDebug::GetTrackRadius(Float_t a, Float_t b, Float_t c) const {
418 // Calculates the track radius using the parameters given by the tilted Rieman fit
420 // Parameters: The three parameters from the Rieman fit
421 // Output: The track radius
424 if(1.0 + b*b - c*a > 0.0)
425 radius = TMath::Sqrt(1.0 + b*b - c*a )/a;
429 //____________________________________________________
430 Float_t AliTRDtrackerDebug::GetTrackCurvature(Float_t a, Float_t b, Float_t c) const {
432 // Calculates the track curvature using the parameters given by the linear fitter
434 // Parameters: the three parameters from the tilted Rieman fitter
435 // Output: the full track curvature
437 Float_t curvature = 1.0 + b*b - c*a;
439 curvature = a / TMath::Sqrt(curvature);
443 //____________________________________________________
444 Float_t AliTRDtrackerDebug::GetDCA(Float_t a, Float_t b, Float_t c) const {
446 // Calculates the Distance to Clostest Approach for the Vertex using the paramters
447 // given by the tilted Rieman fit
449 // Parameters: the three parameters from the tilted Rieman fitter
450 // Output: the Distance to Closest Approach
453 if (1.0 + b*b - c*a > 0.0) {
454 dca = -c / (TMath::Sqrt(1.0 + b*b - c*a) + TMath::Sqrt(1.0 + b*b));
459 //____________________________________________________
460 void AliTRDtrackerDebug::AnalyseMinMax()
463 TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
465 AliError("File TRD.TrackerDebug.root not found!");
468 fTree = (TTree *)(debfile->Get("MakeSeeds0"));
470 AliError("Tree MakeSeeds0 not found in File TRD.TrackerDebug.root.");
473 AliTRDseedV1 *cseed[4] = {0x0, 0x0, 0x0, 0x0};
474 AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0};
475 for(Int_t il = 0; il < 4; il++){
476 fTree->SetBranchAddress(Form("Seed%d.", il), &cseed[il]);
477 fTree->SetBranchAddress(Form("c%d.",il), &c[il]);
479 fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
480 fTree->SetBranchAddress("EventNumber", &fgEventNumber);
481 Int_t entries = fTree->GetEntries();
482 for(Int_t ientry = 0; ientry < entries; ientry++){
483 fTree->GetEntry(ientry);
484 Float_t minmax[2] = { -100.0, 100.0 };
485 for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
486 Float_t max = c[iLayer]->GetZ() + cseed[iLayer]->GetPadLength() * 0.5 + 1.0 - cseed[iLayer]->GetZref(0);
487 if (max < minmax[1]) minmax[1] = max;
488 Float_t min = c[iLayer]->GetZ()-cseed[iLayer]->GetPadLength() * 0.5 - 1.0 - cseed[iLayer]->GetZref(0);
489 if (min > minmax[0]) minmax[0] = min;
491 TTreeSRedirector &cstreamer = *fOutputStreamer;
492 cstreamer << "AnalyseMinMaxLayer"
493 << "EventNumber=" << fgEventNumber
494 << "CandidateNumber=" << fgCandidateNumber
495 << "Min=" << minmax[0]
496 << "Max=" << minmax[1]
501 //____________________________________________________
502 TCanvas* AliTRDtrackerDebug::PlotSeedingConfiguration(const Char_t *direction, Int_t event, Int_t candidate){
504 // Plots the four seeding clusters, the helix fit and the reference Points for
505 // a given combination consisting of eventnr. and candidatenr.
507 // Parameters: -direction (y or z)
509 // -Candidate that has to be plotted
511 const Float_t kxmin = 280;
512 const Float_t kxmax = 380;
513 const Float_t kxdelta = (kxmax - kxmin)/1000;
515 if((strcmp(direction, "y") != 0) && (strcmp(direction, "z") != 0)){
516 AliError(Form("Direction %s does not exist. Abborting!", direction));
520 TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
522 AliError("File TRD.TrackerDebug.root not found!");
525 fTree = (TTree *)(debfile->Get("MakeSeeds0"));
527 AliError("Tree MakeSeeds0 not found in File TRD.TrackerDebug.root.");
531 TGraph *seedcl = new TGraph(4);
532 TGraph *seedRef = new TGraph(4);
533 TGraph *riemanFit = new TGraph(1000);
534 seedcl->SetMarkerStyle(20);
535 seedcl->SetMarkerColor(kRed);
536 seedRef->SetMarkerStyle(2);
538 AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0};
539 AliRieman *rim = 0x0;
540 Bool_t found = kFALSE;
541 for(Int_t il = 0; il < 4; il++) fTree->SetBranchAddress(Form("c%d.",il), &c[il]);
542 fTree->SetBranchAddress("EventNumber", &fgEventNumber);
543 fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
544 fTree->SetBranchAddress("RiemanFitter.", &rim);
545 Int_t entries = fTree->GetEntries();
546 for(Int_t entry = 0; entry < entries; entry++){
547 fTree->GetEntry(entry);
548 if(fgEventNumber < event) continue;
549 if(fgEventNumber > event) break;
550 // EventNumber matching: Do the same for the candidate number
551 if(fgCandidateNumber < candidate) continue;
552 if(fgCandidateNumber > candidate) break;
555 for(Int_t il = 0; il < 4; il++){
556 Float_t cluster = 0.0, reference = 0.0;
557 if(!strcmp(direction, "y")){
558 cluster = c[il]->GetY();
559 reference = rim->GetYat(c[il]->GetX());
562 cluster = c[il]->GetZ();
563 reference = rim->GetZat(c[il]->GetX());
565 seedcl->SetPoint(nPoints, cluster, c[il]->GetX());
566 seedRef->SetPoint(nPoints, reference , c[il]->GetX());
569 // evaluate the fitter Function numerically
571 for(Int_t ipt = 0; ipt < 1000; ipt++){
572 Float_t x = kxmin + ipt * kxdelta;
574 if(!strcmp(direction, "y"))
575 point = rim->GetYat(x);
577 point = rim->GetZat(x);
578 riemanFit->SetPoint(nPoints++, point, x);
580 // We reached the End: break
584 seedcl->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
585 seedRef->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
586 riemanFit->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
587 TCanvas *c1 = new TCanvas();
589 seedRef->Draw("psame");
590 riemanFit->Draw("lpsame");
594 AliError(Form("Combination consisting of event %d and candidate %d not found", event, candidate));
602 //____________________________________________________
603 TCanvas* AliTRDtrackerDebug::PlotFullTrackFit(Int_t event, Int_t candidate, Int_t iteration, const Char_t *direction){
605 // Plots the tracklets (clusters and reference in y direction) and the fitted function for several iterations
606 // in the function ImproveSeedQuality (default is before ImproveSeedQuality)
608 // Parameters: -Event Number
610 // -Iteration Number in ImproveSeedQuality (default: -1 = before ImproveSeedQuality)
611 // -direction (default: y)
612 // Output: -TCanvas (containing the Picture);
614 const Float_t kxmin = 280;
615 const Float_t kxmax = 380;
616 const Float_t kxdelta = (kxmax - kxmin)/1000;
618 if(strcmp(direction, "y") && strcmp(direction, "z")){
619 AliError(Form("Direction %s does not exist. Abborting!", direction));
623 TFile *debfile = TFile::Open("TRD.TrackerDebug.root");
625 AliError("File TRD.TrackerDebug.root not found.");
628 TString *treename = 0x0;
630 treename = new TString("ImproveSeedQuality");
632 treename = new TString("MakeSeeds1");
633 fTree = (TTree *)(debfile->Get(treename->Data()));
635 AliError(Form("Tree %s not found in File TRD.TrackerDebug.root.", treename->Data()));
641 TGraph *fitfun = new TGraph(1000);
642 // Prepare containers
643 Float_t x0[AliTRDtrackerV1::kNPlanes],
644 refP[AliTRDtrackerV1::kNPlanes],
645 clx[AliTRDtrackerV1::kNPlanes * AliTRDtrackerV1::kNTimeBins],
646 clp[AliTRDtrackerV1::kNPlanes * AliTRDtrackerV1::kNTimeBins];
647 Int_t nLayers = 0, ncls = 0;
649 TLinearFitter *fitter = 0x0;
650 AliTRDseedV1 *tracklet[6] = {0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
651 for(Int_t iLayer = 0; iLayer < 6; iLayer++)
652 fTree->SetBranchAddress(Form("S%d.", iLayer), &tracklet[iLayer]);
653 fTree->SetBranchAddress("FitterT.", &fitter);
654 fTree->SetBranchAddress("EventNumber", &fgEventNumber);
655 fTree->SetBranchAddress("CandidateNumber", &fgCandidateNumber);
657 Int_t nEntries = fTree->GetEntriesFast();
658 Bool_t found = kFALSE;
659 for(Int_t entry = 0; entry < nEntries; entry++){
660 fTree->GetEntry(entry);
661 if(fgEventNumber < event) continue;
662 if(fgEventNumber > event) break;
663 // EventNumber matching: Do the same for the candidate number
664 if(fgCandidateNumber < candidate) continue;
665 if(fgCandidateNumber > candidate) break;
668 for(Int_t iLayer = 0; iLayer < 6; iLayer++){
669 if(!tracklet[iLayer]->IsOK()) continue;
670 x0[nLayers] = tracklet[iLayer]->GetX0();
671 if(!strcmp(direction, "y"))
672 refP[nLayers] = tracklet[iLayer]->GetYref(0);
674 refP[nLayers] = tracklet[iLayer]->GetZref(0);
676 for(Int_t itb = 0; itb < 30; itb++){
677 if(!tracklet[iLayer]->IsUsable(itb)) continue;
678 AliTRDcluster *cl = tracklet[iLayer]->GetClusters(itb);
679 if(!strcmp(direction, "y"))
680 clp[ncls] = cl->GetY();
682 clp[ncls] = cl->GetZ();
683 clx[ncls] = cl->GetX();
687 // Add function derived by the tilted Rieman fit (Defined by the curvature)
689 if(!strcmp(direction, "y")){
690 Double_t a = fitter->GetParameter(0);
691 Double_t b = fitter->GetParameter(1);
692 Double_t c = fitter->GetParameter(2);
693 Double_t curvature = 1.0 + b*b - c*a;
694 if (curvature > 0.0) {
695 curvature = a / TMath::Sqrt(curvature);
697 // Numerical evaluation of the function:
698 for(Int_t ipt = 0; ipt < 1000; ipt++){
699 Float_t x = kxmin + ipt * kxdelta;
700 Double_t res = (x * a + b); // = (x - x0)/y0
702 res = 1.0 - c * a + b * b - res; // = (R^2 - (x - x0)^2)/y0^2
705 res = TMath::Sqrt(res);
708 fitfun->SetPoint(nPoints++, y, x);
712 Double_t offset = fitter->GetParameter(3);
713 Double_t slope = fitter->GetParameter(4);
714 // calculate the reference x (defined as medium between layer 2 and layer 3)
715 // same procedure as in the tracker code
716 Float_t medx = 0, xref = 0;
717 Int_t startIndex = 5, nDistances = 0;
718 for(Int_t iLayer = 5; iLayer > 0; iLayer--){
719 if(tracklet[iLayer]->IsOK() && tracklet[iLayer - 1]->IsOK()){
720 medx += tracklet[iLayer]->GetX0() - tracklet[iLayer - 1]->GetX0();
721 startIndex = iLayer - 1;
729 Float_t xpos[2]; memset(xpos, 0, sizeof(Float_t) * 2);
730 Int_t ien = 0, idiff = 0;
731 for(Int_t iLayer = 5; iLayer > 0; iLayer--){
732 if(tracklet[iLayer]->IsOK()){
733 xpos[ien++] = tracklet[iLayer]->GetX0();
741 medx = (xpos[0] - xpos[1])/idiff;
743 xref = tracklet[startIndex]->GetX0() + medx * (2.5 - startIndex) - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
745 for(Int_t ipt = 0; ipt < 1000; ipt++){
746 Float_t x = kxmin + ipt * kxdelta;
747 Float_t z = offset + slope * (x - xref);
748 fitfun->SetPoint(nPoints++, z, x);
754 TGraph *trGraph = new TGraph(ncls);
755 TGraph *refPoints = new TGraph(nLayers);
756 trGraph->SetMarkerStyle(20);
757 trGraph->SetMarkerColor(kRed);
758 refPoints->SetMarkerStyle(21);
759 refPoints->SetMarkerColor(kBlue);
761 for(Int_t iLayer = 0; iLayer < nLayers; iLayer++)
762 refPoints->SetPoint(iLayer, refP[iLayer], x0[iLayer]);
763 for(Int_t icls = 0; icls < ncls; icls++)
764 trGraph->SetPoint(icls, clp[icls], clx[icls]);
765 TCanvas *c1 = new TCanvas();
766 trGraph->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
767 refPoints->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
768 fitfun->SetTitle(Form("Event %d, Candidate %d\n", fgEventNumber, fgCandidateNumber));
770 refPoints->Draw("psame");
771 fitfun->Draw("lpsame");
775 AliError(Form("Combination consisting of event %d and candidate %d not found", event, candidate));