]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/qaRec/AliTRDtrackingResolution.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
- fix connection to the visualization framework
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDtrackingResolution.cxx
CommitLineData
77203477 1/**************************************************************************
d2381af5 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-commercialf 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**************************************************************************/
77203477 15
16/* $Id: AliTRDtrackingResolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
17
18////////////////////////////////////////////////////////////////////////////
19// //
017bd6af 20// TRD tracking resolution //
21//
22// The class performs resolution and residual studies
23// of the TRD tracks for the following quantities :
24// - spatial position (y, [z])
25// - angular (phi) tracklet
26// - momentum at the track level
27//
28// The class has to be used for regular detector performance checks using the official macros:
29// - $ALICE_ROOT/TRD/qaRec/run.C
30// - $ALICE_ROOT/TRD/qaRec/makeResults.C
31//
32// For stand alone usage please refer to the following example:
33// {
34// gSystem->Load("libANALYSIS.so");
35// gSystem->Load("libTRDqaRec.so");
36// AliTRDtrackingResolution *res = new AliTRDtrackingResolution();
37// //res->SetMCdata();
38// //res->SetVerbose();
39// //res->SetVisual();
40// res->Load("TRD.TaskResolution.root");
41// if(!res->PostProcess()) return;
42// res->GetRefFigure(0);
43// }
44//
77203477 45// Authors: //
017bd6af 46// Alexandru Bercuci <A.Bercuci@gsi.de> //
77203477 47// Markus Fasel <M.Fasel@gsi.de> //
48// //
49////////////////////////////////////////////////////////////////////////////
50
aaf47b30 51#include <cstring>
52
017bd6af 53#include <TSystem.h>
77203477 54#include <TObjArray.h>
7102d1b1 55#include <TH2.h>
56#include <TH1.h>
57#include <TF1.h>
017bd6af 58#include <TCanvas.h>
77203477 59#include <TProfile.h>
7102d1b1 60#include <TGraphErrors.h>
77203477 61#include <TMath.h>
aaf47b30 62#include "TTreeStream.h"
63#include "TGeoManager.h"
77203477 64
65#include "AliAnalysisManager.h"
77203477 66#include "AliTrackReference.h"
aaf47b30 67#include "AliTrackPointArray.h"
68#include "AliCDBManager.h"
69
9605ce80 70#include "AliTRDSimParam.h"
71#include "AliTRDgeometry.h"
72#include "AliTRDpadPlane.h"
aaf47b30 73#include "AliTRDcluster.h"
74#include "AliTRDseedV1.h"
75#include "AliTRDtrackV1.h"
76#include "AliTRDtrackerV1.h"
77#include "AliTRDReconstructor.h"
78#include "AliTRDrecoParam.h"
77203477 79
80#include "AliTRDtrackInfo/AliTRDtrackInfo.h"
81#include "AliTRDtrackingResolution.h"
82
77203477 83ClassImp(AliTRDtrackingResolution)
84
85//________________________________________________________
3d86166d 86AliTRDtrackingResolution::AliTRDtrackingResolution()
87 :AliTRDrecoTask("Resolution", "Tracking Resolution")
017bd6af 88 ,fStatus(0)
aaf47b30 89 ,fReconstructor(0x0)
9605ce80 90 ,fGeo(0x0)
b718144c 91 ,fGraphS(0x0)
92 ,fGraphM(0x0)
77203477 93{
aaf47b30 94 fReconstructor = new AliTRDReconstructor();
95 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
9605ce80 96 fGeo = new AliTRDgeometry();
de520d8f 97 InitFunctorList();
77203477 98}
99
ed383798 100//________________________________________________________
101AliTRDtrackingResolution::~AliTRDtrackingResolution()
102{
b718144c 103 if(fGraphS){fGraphS->Delete(); delete fGraphS;}
104 if(fGraphM){fGraphM->Delete(); delete fGraphM;}
9605ce80 105 delete fGeo;
ed383798 106 delete fReconstructor;
2c0cf367 107 if(gGeoManager) delete gGeoManager;
ed383798 108}
109
77203477 110
111//________________________________________________________
112void AliTRDtrackingResolution::CreateOutputObjects()
113{
114 // spatial resolution
77203477 115 OpenFile(0, "RECREATE");
39779ce6 116
3d86166d 117 fContainer = Histos();
77203477 118}
119
de520d8f 120// //________________________________________________________
121// void AliTRDtrackingResolution::Exec(Option_t *)
122// {
123// // spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire})
124// // angular Resolution: res = Tracklet angle - TrackRef Angle
b718144c 125//
de520d8f 126// Int_t nTrackInfos = fTracks->GetEntriesFast();
127// if(fDebugLevel>=2 && nTrackInfos){
128// printf("Event[%d] TrackInfos[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTrackInfos);
129// }
130// const Int_t kNLayers = AliTRDgeometry::kNlayer;
131// Int_t pdg, nly, ncrs;
132// Double_t p, dy, theta/*, dphi, dymc, dzmc, dphimc*/;
133// Float_t fP[kNLayers], fX[kNLayers], fY[kNLayers], fZ[kNLayers], fPhi[kNLayers], fTheta[kNLayers]; // phi/theta angle per layer
134// Bool_t fMCMap[kNLayers], fLayerMap[kNLayers]; // layer map
b718144c 135//
de520d8f 136// AliTRDpadPlane *pp = 0x0;
137// AliTrackPoint tr[kNLayers], tk[kNLayers];
138// AliExternalTrackParam *fOp = 0x0;
139// AliTRDtrackV1 *fTrack = 0x0;
140// AliTRDtrackInfo *fInfo = 0x0;
141// for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){
142// // check if ESD and MC-Information are available
143// if(!(fInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iTI)))) continue;
144// if(!(fTrack = fInfo->GetTrack())) continue;
145// if(!(fOp = fInfo->GetOuterParam())) continue;
146// pdg = fInfo->GetPDG();
147// nly = 0; ncrs = 0; theta = 0.;
b718144c 148//
de520d8f 149// if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs());
b718144c 150//
de520d8f 151// p = fOp->P();
152// Int_t npts = 0;
153// memset(fP, 0, kNLayers*sizeof(Float_t));
154// memset(fX, 0, kNLayers*sizeof(Float_t));
155// memset(fY, 0, kNLayers*sizeof(Float_t));
156// memset(fZ, 0, kNLayers*sizeof(Float_t));
157// memset(fPhi, 0, kNLayers*sizeof(Float_t));
158// memset(fTheta, 0, kNLayers*sizeof(Float_t));
159// memset(fLayerMap, 0, kNLayers*sizeof(Bool_t));
160// memset(fMCMap, 0, kNLayers*sizeof(Bool_t));
161// AliTRDseedV1 *fTracklet = 0x0;
162// for(Int_t iplane = 0; iplane < kNLayers; iplane++){
163// if(!(fTracklet = fTrack->GetTracklet(iplane))) continue;
164// if(!fTracklet->IsOK()) continue;
b718144c 165//
de520d8f 166// // Book point arrays
167// fLayerMap[iplane] = kTRUE;
168// tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.);
169// tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0));
170// npts++;
b718144c 171//
de520d8f 172// if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN());
173//
174// // define reference values
175// fP[iplane] = p;
176// fX[iplane] = fTracklet->GetX0();
177// fY[iplane] = fTracklet->GetYref(0);
178// fZ[iplane] = fTracklet->GetZref(0);
179// fPhi[iplane] = TMath::ATan(fTracklet->GetYref(1));
180// fTheta[iplane] = TMath::ATan(fTracklet->GetZref(1));
181//
182//
183// // RESOLUTION (compare to MC)
184// if(HasMCdata()){
185// if(fInfo->GetNTrackRefs() >= 2){
186// Double_t pmc, ymc, zmc, phiMC, thetaMC;
187// if(Resolution(fTracklet, fInfo, pmc, ymc, zmc, phiMC, thetaMC)){
188// fMCMap[iplane] = kTRUE;
189// fP[iplane] = pmc;
190// fY[iplane] = ymc;
191// fZ[iplane] = zmc;
192// fPhi[iplane] = phiMC;
193// fTheta[iplane] = thetaMC;
194// }
195// }
b718144c 196// }
de520d8f 197// Float_t phi = fPhi[iplane]*TMath::RadToDeg();
198// theta += fTheta[iplane]; nly++;
199// if(fTracklet->GetNChange()) ncrs++;
b718144c 200//
de520d8f 201// // Do clusters residuals
202// if(!fTracklet->Fit(kFALSE)) continue;
203// AliTRDcluster *c = 0x0;
204// for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
205// if(!(c = fTracklet->GetClusters(ic))) continue;
b718144c 206//
de520d8f 207// dy = fTracklet->GetYat(c->GetX()) - c->GetY();
208// ((TH2I*)fContainer->At(kClusterYResidual))->Fill(phi, dy);
209//
210// if(fDebugLevel>=2){
211// Float_t q = c->GetQ();
212// // Get z-position with respect to anode wire
213// AliTRDSimParam *simParam = AliTRDSimParam::Instance();
214// Int_t det = c->GetDetector();
215// Float_t x = c->GetX();
216// Float_t z = fZ[iplane]-(fX[iplane]-x)*TMath::Tan(fTheta[iplane]);
217// Int_t stack = fGeo->GetStack(det);
218// pp = fGeo->GetPadPlane(iplane, stack);
219// Float_t row0 = pp->GetRow0();
220// Float_t d = row0 - z + simParam->GetAnodeWireOffset();
221// d -= ((Int_t)(2 * d)) / 2.0;
222// if (d > 0.25) d = 0.5 - d;
223//
224// (*fDebugStream) << "ResidualClusters"
225// << "ly=" << iplane
226// << "stk=" << stack
227// << "pdg=" << pdg
228// << "phi=" << fPhi[iplane]
229// << "tht=" << fTheta[iplane]
230// << "q=" << q
231// << "x=" << x
232// << "z=" << z
233// << "d=" << d
234// << "dy=" << dy
235// << "\n";
236// }
b718144c 237// }
de520d8f 238// pp = 0x0;
239// }
240// if(nly) theta /= nly;
241// if(fDebugLevel>=1){
242// (*fDebugStream) << "TrackStatistics"
243// << "nly=" << nly
244// << "ncrs=" << ncrs
245// << "tht=" << theta
246// << "\n";
b718144c 247// }
de520d8f 248//
249//
250// // // this protection we might drop TODO
251// // if(fTrack->GetNumberOfTracklets() < 6) continue;
252// //
253// // AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr);
254// // Int_t iref = 0;
255// // for(Int_t ip=0; ip<kNLayers; ip++){
256// // if(!fLayerMap[ip]) continue;
257// // fTracklet = fTrack->GetTracklet(ip);
258// // // recalculate fit based on the new tilt correction
259// // fTracklet->Fit();
260// //
261// // dy = fTracklet->GetYfit(0) - tr[iref].GetY();
262// // ((TH2I*)fContainer->At(kTrackletRiemanYResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dy);
263// //
264// // dphi = fTracklet->GetYfit(1)- fTracklet->GetYref(1);
265// // ((TH2I*)fContainer->At(kTrackletRiemanAngleResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dphi);
266// //
267// // if(HasMCdata()){
268// // dymc = fY[ip] - tr[iref].GetY();
269// // ((TH2I*)fContainer->At(kTrackRYResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dymc);
270// //
271// // dzmc = fZ[ip] - tr[iref].GetZ();
272// // ((TH2I*)fContainer->At(kTrackRZResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dzmc);
273// //
274// // dphimc = fPhi[ip] - fTracklet->GetYfit(1);
275// // ((TH2I*)fContainer->At(kTrackRAngleResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dphimc);
276// // }
277// //
278// // iref++;
279// //
280// // if(fDebugLevel>=1){
281// // (*fDebugStream) << "RiemannTrack"
282// // << "ly=" << ip
283// // << "mc=" << fMCMap[ip]
284// // << "p=" << fP[ip]
285// // << "phi=" << fPhi[ip]
286// // << "tht=" << fTheta[ip]
287// // << "dy=" << dy
288// // << "dphi=" << dphi
289// // << "dymc=" << dymc
290// // << "dzmc=" << dzmc
291// // << "dphimc="<< dphimc
292// // << "\n";
293// // }
294// // }
295//
296// // if(!gGeoManager) TGeoManager::Import("geometry.root");
297// // AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr);
298// // for(Int_t ip=0; ip<nc; ip++){
299// // dy = cl[ip].GetY() - tr[ip].GetY();
300// // ((TH2I*)fContainer->At(kTrackletKalmanYResidual))->Fill(phi*TMath::RadToDeg(), dy);
301// // dz = cl[ip].GetZ() - tr[ip].GetZ();
302// // if(fDebugLevel>=1){
303// // (*fDebugStream) << "KalmanTrack"
304// // << "dy=" << dy
305// // << "dz=" << dz
306// // /* << "phi=" << phi
307// // << "theta=" << theta
308// // << "dphi=" << dphi*/
309// // << "\n";
310// // }
311// // }
312//
313//
314// }
315// PostData(0, fContainer);
316// }
aaf47b30 317
de520d8f 318//________________________________________________________
319TH1* AliTRDtrackingResolution::PlotClusterResiduals(const AliTRDtrackV1 *track)
320{
74b2e03d 321 if(track) fTrack = track;
322 if(!fTrack){
323 AliWarning("No Track defined.");
324 return 0x0;
de520d8f 325 }
326 TH1 *h = 0x0;
327 if(!(h = ((TH2I*)fContainer->At(kClusterResidual)))){
328 AliWarning("No output histogram defined.");
329 return 0x0;
330 }
331
332 Int_t pdg = fMC ? fMC->GetPDG() : 0;
333 Float_t x0, y0, z0, dy, dydx, dzdx;
334 AliTRDseedV1 *fTracklet = 0x0;
335 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
336 if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
337 if(!fTracklet->IsOK()) continue;
de520d8f 338 x0 = fTracklet->GetX0();
339
340 // retrive the track angle with the chamber
341 if(fMC) fMC->GetDirections(x0, y0, z0, dydx, dzdx);
342 else{
343 y0 = fTracklet->GetYref(0);
344 z0 = fTracklet->GetZref(0);
345 dydx = fTracklet->GetYref(1);
346 dzdx = fTracklet->GetZref(1);
347 }
348
107fde80 349 AliTRDseedV1 trklt(*fTracklet);
350 if(!trklt.Fit(kFALSE)) continue;
351
de520d8f 352 AliTRDcluster *c = 0x0;
353 fTracklet->ResetClusterIter(kFALSE);
354 while((c = fTracklet->PrevCluster())){
107fde80 355 dy = trklt.GetYat(c->GetX()) - c->GetY();
de520d8f 356 h->Fill(dydx, dy);
357
358 if(fDebugLevel>=1){
359 Float_t q = c->GetQ();
360 // Get z-position with respect to anode wire
361 AliTRDSimParam *simParam = AliTRDSimParam::Instance();
362 Int_t det = c->GetDetector();
363 Float_t x = c->GetX();
364 Float_t z = z0-(x0-x)*dzdx;
365 Int_t istk = fGeo->GetStack(det);
366 AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
367 Float_t row0 = pp->GetRow0();
368 Float_t d = row0 - z + simParam->GetAnodeWireOffset();
369 d -= ((Int_t)(2 * d)) / 2.0;
370 if (d > 0.25) d = 0.5 - d;
371
372 (*fDebugStream) << "ClusterResidual"
373 << "ly=" << ily
374 << "stk=" << istk
375 << "pdg=" << pdg
376 << "dydx=" << dydx
377 << "dzdx=" << dzdx
378 << "q=" << q
379 << "x=" << x
380 << "z=" << z
381 << "d=" << d
382 << "dy=" << dy
383 << "\n";
384 }
385 }
386 }
387 return h;
388}
389
390//________________________________________________________
391TH1* AliTRDtrackingResolution::PlotResolution(const AliTRDtrackV1 *track)
392{
393 if(!fMC){
74b2e03d 394 AliWarning("No MC defined. Results will not be available.");
de520d8f 395 return 0x0;
396 }
74b2e03d 397 if(track) fTrack = track;
398 if(!fTrack){
399 AliWarning("No Track defined.");
400 return 0x0;
de520d8f 401 }
402 TH1 *h = 0x0;
403 if(!(h = ((TH2I*)fContainer->At(kClusterResolution)))){
404 AliWarning("No output histogram defined.");
405 return 0x0;
406 }
407 if(!(h = ((TH2I*)fContainer->At(kTrackletYResolution)))){
408 AliWarning("No output histogram defined.");
409 return 0x0;
410 }
411 if(!(h = ((TH2I*)fContainer->At(kTrackletZResolution)))){
412 AliWarning("No output histogram defined.");
413 return 0x0;
414 }
415 if(!(h = ((TH2I*)fContainer->At(kTrackletAngleResolution)))){
416 AliWarning("No output histogram defined.");
417 return 0x0;
418 }
419 //printf("called for %d tracklets ... \n", fTrack->GetNumberOfTracklets());
420 Int_t pdg = fMC->GetPDG(), det=-1;
421 Float_t x0, y0, z0, dy, dydx, dzdx;
422 AliTRDseedV1 *fTracklet = 0x0;
423 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
424 if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
425 if(!fTracklet->IsOK()) continue;
426 //printf("process layer[%d]\n", ily);
427 // retrive the track position and direction within the chamber
428 det = fTracklet->GetDetector();
429 x0 = fTracklet->GetX0();
430 fMC->GetDirections(x0, y0, z0, dydx, dzdx);
431
432 // recalculate tracklet based on the MC info
433 AliTRDseedV1 tt(*fTracklet);
434 tt.SetZref(0, z0);
435 tt.SetZref(1, dzdx);
436 if(!tt.Fit()) continue;
437 dy = tt.GetYfit(0) - y0;
438 Float_t dphi = TMath::ATan(tt.GetYfit(1)) - TMath::ATan(dydx);
439 Float_t dz = 100.;
440 Bool_t cross = fTracklet->GetNChange();
441 if(cross){
442 Double_t *xyz = tt.GetCrossXYZ();
443 dz = xyz[2] - (z0 - (x0 - xyz[0])*dzdx) ;
444 ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(dzdx, dz);
445 }
446
447 // Fill Histograms
448 ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(dydx, dy);
449 ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(dydx, dphi*TMath::RadToDeg());
39779ce6 450
de520d8f 451 // Fill Debug stream
452 if(fDebugLevel>=1){
06b2847e 453 Float_t p = fMC->GetTrackRef() ? fMC->GetTrackRef()->P() : -1.;
de520d8f 454 (*fDebugStream) << "TrkltResolution"
455 << "det=" << det
456 << "pdg=" << pdg
457 << "p=" << p
458 << "ymc=" << y0
459 << "zmc=" << z0
460 << "dydx=" << dydx
461 << "dzdx=" << dzdx
462 << "cross=" << cross
463 << "dy=" << dy
464 << "dz=" << dz
465 << "dphi=" << dphi
466 << "\n";
467 }
39779ce6 468
de520d8f 469 Int_t istk = AliTRDgeometry::GetStack(det);
470 AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
cf194b94 471 Float_t zr0 = pp->GetRow0() + AliTRDSimParam::Instance()->GetAnodeWireOffset();
de520d8f 472 Float_t tilt = fTracklet->GetTilt();
473
474 AliTRDcluster *c = 0x0;
475 fTracklet->ResetClusterIter(kFALSE);
476 while((c = fTracklet->PrevCluster())){
477 Float_t q = TMath::Abs(c->GetQ());
478 Float_t xc = c->GetX();
479 Float_t yc = c->GetY();
480 Float_t zc = c->GetZ();
481 Float_t dx = x0 - xc;
482 Float_t yt = y0 - dx*dydx;
483 Float_t zt = z0 - dx*dzdx;
484 dy = yt - (yc - tilt*(zc-zt));
485
486 // Fill Histograms
487 if(q>100.) ((TH2I*)fContainer->At(kClusterResolution))->Fill(dydx, dy);
488
489 // Fill Debug Tree
490 if(fDebugLevel>=1){
cf194b94 491 Float_t d = zr0 - zt;
de520d8f 492 d -= ((Int_t)(2 * d)) / 2.0;
493 if (d > 0.25) d = 0.5 - d;
494
495 (*fDebugStream) << "ClusterResolution"
496 << "ly=" << ily
497 << "stk=" << istk
498 << "pdg=" << pdg
499 << "dydx=" << dydx
500 << "dzdx=" << dzdx
501 << "q=" << q
502 << "d=" << d
503 << "dy=" << dy
504 << "\n";
505 }
506 }
77203477 507 }
de520d8f 508 return h;
77203477 509}
510
de520d8f 511
d85cd79c 512//________________________________________________________
a391a274 513void AliTRDtrackingResolution::GetRefFigure(Int_t ifig)
d85cd79c 514{
a391a274 515 TAxis *ax = 0x0;
516 TGraphErrors *g = 0x0;
d85cd79c 517 switch(ifig){
de520d8f 518 case kClusterResidual:
519 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
a391a274 520 g->Draw("apl");
521 ax = g->GetHistogram()->GetYaxis();
b718144c 522 ax->SetRangeUser(-.5, 1.);
a391a274 523 ax->SetTitle("Clusters Y Residuals #sigma/#mu [mm]");
524 ax = g->GetHistogram()->GetXaxis();
017bd6af 525 ax->SetTitle("tg(#phi)");
de520d8f 526 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
017bd6af 527 g->Draw("pl");
528 return;
de520d8f 529 case kClusterResolution:
530 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
a391a274 531 ax = g->GetHistogram()->GetYaxis();
b718144c 532 ax->SetRangeUser(-.5, 1.);
533 ax->SetTitle("Cluster Y Resolution #sigma/#mu [mm]");
534 ax = g->GetHistogram()->GetXaxis();
017bd6af 535 ax->SetTitle("tg(#phi)");
b718144c 536 g->Draw("apl");
de520d8f 537 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
017bd6af 538 g->Draw("pl");
539 return;
b718144c 540 case kTrackletYResolution:
de520d8f 541 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
b718144c 542 ax = g->GetHistogram()->GetYaxis();
543 ax->SetRangeUser(-.5, 1.);
a391a274 544 ax->SetTitle("Tracklet Y Resolution #sigma/#mu [mm]");
545 ax = g->GetHistogram()->GetXaxis();
de520d8f 546 ax->SetTitle("#tg(phi)");
547 g->Draw("apl");
548 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
549 g->Draw("pl");
550 return;
551 case kTrackletZResolution:
552 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
553 ax = g->GetHistogram()->GetYaxis();
554 ax->SetRangeUser(-.5, 1.);
555 ax->SetTitle("Tracklet Z Resolution #sigma/#mu [mm]");
556 ax = g->GetHistogram()->GetXaxis();
557 ax->SetTitle("#tg(theta)");
a391a274 558 g->Draw("apl");
de520d8f 559 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
017bd6af 560 g->Draw("pl");
561 return;
b718144c 562 case kTrackletAngleResolution:
de520d8f 563 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
a391a274 564 ax = g->GetHistogram()->GetYaxis();
017bd6af 565 ax->SetRangeUser(-.05, .2);
566 ax->SetTitle("Tracklet Angular Resolution #sigma/#mu [deg]");
a391a274 567 ax = g->GetHistogram()->GetXaxis();
de520d8f 568 ax->SetTitle("tg(#phi)");
a391a274 569 g->Draw("apl");
de520d8f 570 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
017bd6af 571 g->Draw("pl");
572 return;
b718144c 573 default:
574 AliInfo(Form("Reference plot [%d] not implemented yet", ifig));
017bd6af 575 return;
d85cd79c 576 }
017bd6af 577 AliInfo(Form("Reference plot [%d] missing result", ifig));
d85cd79c 578}
579
39779ce6 580
77203477 581//________________________________________________________
d85cd79c 582Bool_t AliTRDtrackingResolution::PostProcess()
7102d1b1 583{
d85cd79c 584 //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
585 if (!fContainer) {
586 Printf("ERROR: list not available");
587 return kFALSE;
3d86166d 588 }
b718144c 589 fNRefFigures = fContainer->GetEntriesFast();
590 if(!fGraphS){
591 fGraphS = new TObjArray(fNRefFigures);
592 fGraphS->SetOwner();
593 }
594 if(!fGraphM){
595 fGraphM = new TObjArray(fNRefFigures);
596 fGraphM->SetOwner();
597 }
7102d1b1 598
d2381af5 599 TH2I *h2 = 0x0;
600 TH1D *h = 0x0;
d85cd79c 601 TGraphErrors *gm = 0x0, *gs = 0x0;
b718144c 602
603 // define models
d2381af5 604 TF1 f("f1", "gaus", -.5, .5);
b718144c 605
017bd6af 606 TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
607
b718144c 608 TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
874acced 609
017bd6af 610 TCanvas *c = 0x0;
611 if(IsVisual()) c = new TCanvas("c", Form("%s Visual", GetName()), 500, 500);
612 char opt[5];
613 sprintf(opt, "%c%c", IsVerbose() ? ' ' : 'Q', IsVisual() ? ' ': 'N');
614
615
874acced 616 //PROCESS RESIDUAL DISTRIBUTIONS
617
618 // Clusters residuals
de520d8f 619 h2 = (TH2I *)(fContainer->At(kClusterResidual));
b718144c 620 gm = new TGraphErrors(h2->GetNbinsX());
017bd6af 621 gm->SetLineColor(kBlue);
622 gm->SetMarkerStyle(7);
623 gm->SetMarkerColor(kBlue);
b718144c 624 gm->SetNameTitle("clm", "");
de520d8f 625 fGraphM->AddAt(gm, kClusterResidual);
b718144c 626 gs = new TGraphErrors(h2->GetNbinsX());
627 gs->SetLineColor(kRed);
628 gs->SetMarkerStyle(23);
629 gs->SetMarkerColor(kRed);
630 gs->SetNameTitle("cls", "");
de520d8f 631 fGraphS->AddAt(gs, kClusterResidual);
874acced 632 for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
633 Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
874acced 634 h = h2->ProjectionY("py", ibin, ibin);
017bd6af 635 AdjustF1(h, &fc);
636
637 if(IsVisual()){c->cd(); c->SetLogy();}
638 h->Fit(&fc, opt, "", -0.5, 0.5);
639 if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
640
641 gm->SetPoint(ibin - 1, TMath::Tan(phi*TMath::DegToRad()), 10.*fc.GetParameter(1));
b2238a96 642 gm->SetPointError(ibin - 1, 0., 10.*fc.GetParError(1));
017bd6af 643 gs->SetPoint(ibin - 1, TMath::Tan(phi*TMath::DegToRad()), 10.*fc.GetParameter(2));
b2238a96 644 gs->SetPointError(ibin - 1, 0., 10.*fc.GetParError(2));
874acced 645 }
d2381af5 646
874acced 647
648 //PROCESS RESOLUTION DISTRIBUTIONS
b718144c 649
874acced 650 if(HasMCdata()){
b718144c 651 // cluster y resolution
de520d8f 652 h2 = (TH2I*)fContainer->At(kClusterResolution);
b718144c 653 gm = new TGraphErrors(h2->GetNbinsX());
017bd6af 654 gm->SetLineColor(kBlue);
655 gm->SetMarkerStyle(7);
656 gm->SetMarkerColor(kBlue);
b718144c 657 gm->SetNameTitle("clym", "");
de520d8f 658 fGraphM->AddAt(gm, kClusterResolution);
b718144c 659 gs = new TGraphErrors(h2->GetNbinsX());
660 gs->SetLineColor(kRed);
661 gs->SetMarkerStyle(23);
662 gs->SetMarkerColor(kRed);
663 gs->SetNameTitle("clys", "");
de520d8f 664 fGraphS->AddAt(gs, kClusterResolution);
b718144c 665 for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
666 h = h2->ProjectionY("py", iphi, iphi);
017bd6af 667 AdjustF1(h, &fb);
668
669 if(IsVisual()){c->cd(); c->SetLogy();}
670 h->Fit(&fb, opt, "", -0.5, 0.5);
671 if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
672
b718144c 673 Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
674 Int_t jphi = iphi -1;
017bd6af 675 gm->SetPoint(jphi, TMath::Tan(phi*TMath::DegToRad()), 10.*fb.GetParameter(1));
b2238a96 676 gm->SetPointError(jphi, 0., 10.*fb.GetParError(1));
017bd6af 677 gs->SetPoint(jphi, TMath::Tan(phi*TMath::DegToRad()), 10.*fb.GetParameter(2));
b2238a96 678 gs->SetPointError(jphi, 0., 10.*fb.GetParError(2));
b718144c 679 }
680
874acced 681 // tracklet y resolution
3d86166d 682 h2 = (TH2I*)fContainer->At(kTrackletYResolution);
b718144c 683 gm = new TGraphErrors(h2->GetNbinsX());
017bd6af 684 gm->SetLineColor(kBlue);
685 gm->SetMarkerStyle(7);
686 gm->SetMarkerColor(kBlue);
b718144c 687 gm->SetNameTitle("trkltym", "");
688 fGraphM->AddAt(gm, kTrackletYResolution);
689 gs = new TGraphErrors(h2->GetNbinsX());
690 gs->SetLineColor(kRed);
691 gs->SetMarkerStyle(23);
692 gs->SetMarkerColor(kRed);
693 gs->SetNameTitle("trkltys", "");
694 fGraphS->AddAt(gs, kTrackletYResolution);
d2381af5 695 for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
d2381af5 696 h = h2->ProjectionY("py", iphi, iphi);
017bd6af 697 AdjustF1(h, &fb);
698
699 if(IsVisual()){c->cd(); c->SetLogy();}
700 h->Fit(&fb, opt, "", -0.5, 0.5);
701 if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
702
b718144c 703 Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
d2381af5 704 Int_t jphi = iphi -1;
017bd6af 705 gm->SetPoint(jphi, phi, 10.*fb.GetParameter(1));
706 gm->SetPointError(jphi, 0., 10.*fb.GetParError(1));
707 gs->SetPoint(jphi, phi, 10.*fb.GetParameter(2));
708 gs->SetPointError(jphi, 0., 10.*fb.GetParError(2));
d2381af5 709 }
d2381af5 710
874acced 711 // tracklet phi resolution
3d86166d 712 h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
b718144c 713 gm = new TGraphErrors(h2->GetNbinsX());
017bd6af 714 gm->SetLineColor(kBlue);
715 gm->SetMarkerStyle(7);
716 gm->SetMarkerColor(kBlue);
b718144c 717 gm->SetNameTitle("trkltym", "");
718 fGraphM->AddAt(gm, kTrackletAngleResolution);
719 gs = new TGraphErrors(h2->GetNbinsX());
720 gs->SetLineColor(kRed);
721 gs->SetMarkerStyle(23);
722 gs->SetMarkerColor(kRed);
723 gs->SetNameTitle("trkltys", "");
724 fGraphS->AddAt(gs, kTrackletAngleResolution);
d2381af5 725 for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
d2381af5 726 h = h2->ProjectionY("py", iphi, iphi);
017bd6af 727
728 if(IsVisual()){c->cd(); c->SetLogy();}
729 h->Fit(&f, opt, "", -0.5, 0.5);
730 if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
731
b718144c 732 Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
d2381af5 733 Int_t jphi = iphi -1;
734 gm->SetPoint(jphi, phi, f.GetParameter(1));
735 gm->SetPointError(jphi, 0., f.GetParError(1));
736 gs->SetPoint(jphi, phi, f.GetParameter(2));
737 gs->SetPointError(jphi, 0., f.GetParError(2));
738 }
d2381af5 739 }
017bd6af 740 if(c) delete c;
765bd0ab 741
d85cd79c 742 return kTRUE;
743}
744
745
746//________________________________________________________
747void AliTRDtrackingResolution::Terminate(Option_t *)
748{
749 if(fDebugStream){
750 delete fDebugStream;
751 fDebugStream = 0x0;
752 fDebugLevel = 0;
753 }
754 if(HasPostProcess()) PostProcess();
874acced 755}
d2381af5 756
3c3d9ff1 757//________________________________________________________
017bd6af 758void AliTRDtrackingResolution::AdjustF1(TH1 *h, TF1 *f)
3c3d9ff1 759{
760// Helper function to avoid duplication of code
761// Make first guesses on the fit parameters
762
763 // find the intial parameters of the fit !! (thanks George)
764 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
765 Double_t sum = 0.;
766 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
767 f->SetParLimits(0, 0., 3.*sum);
768 f->SetParameter(0, .9*sum);
769
017bd6af 770 f->SetParLimits(1, -.2, .2);
3c3d9ff1 771 f->SetParameter(1, 0.);
772
017bd6af 773 f->SetParLimits(2, 0., 4.e-1);
3c3d9ff1 774 f->SetParameter(2, 2.e-2);
017bd6af 775 if(f->GetNpar() <= 4) return;
3c3d9ff1 776
777 f->SetParLimits(3, 0., sum);
778 f->SetParameter(3, .1*sum);
779
780 f->SetParLimits(4, -.3, .3);
781 f->SetParameter(4, 0.);
782
783 f->SetParLimits(5, 0., 1.e2);
784 f->SetParameter(5, 2.e-1);
3c3d9ff1 785}
786
874acced 787//________________________________________________________
3d86166d 788TObjArray* AliTRDtrackingResolution::Histos()
874acced 789{
cf194b94 790 if(fContainer) return fContainer;
791
792 fContainer = new TObjArray(5);
793
794 // cluster to tracklet residuals [2]
107fde80 795 fContainer->AddAt(new TH2I("fYClRes", "Clusters Residuals", 21, -.33, .33, 100, -.5, .5), kClusterResidual);
cf194b94 796// // tracklet to Riemann fit residuals [2]
797// fContainer->AddAt(new TH2I("fYTrkltRRes", "Tracklet Riemann Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanYResidual);
798// fContainer->AddAt(new TH2I("fAngleTrkltRRes", "Tracklet Riemann Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanAngleResidual);
799// fContainer->AddAt(new TH2I("fYTrkltKRes", "Tracklet Kalman Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanYResidual);
800// fContainer->AddAt(new TH2I("fAngleTrkltKRes", "Tracklet Kalman Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanAngleResidual);
801
802 // Resolution histos
803 if(HasMCdata()){
804 // cluster y resolution [0]
107fde80 805 fContainer->AddAt(new TH2I("fCY", "Cluster Resolution", 31, -.48, .48, 100, -.5, .5), kClusterResolution);
cf194b94 806 // tracklet y resolution [0]
107fde80 807 fContainer->AddAt(new TH2I("fY", "Tracklet Resolution", 31, -.48, .48, 100, -.5, .5), kTrackletYResolution);
cf194b94 808 // tracklet y resolution [0]
107fde80 809 fContainer->AddAt(new TH2I("fY", "Tracklet Resolution", 31, -.48, .48, 100, -.5, .5), kTrackletZResolution);
cf194b94 810 // tracklet angular resolution [1]
107fde80 811 fContainer->AddAt(new TH2I("fPhi", "Tracklet Angular Resolution", 31, -.48, .48, 100, -10., 10.), kTrackletAngleResolution);
cf194b94 812
813// // Riemann track resolution [y, z, angular]
814// fContainer->AddAt(new TH2I("fYRT", "Track Riemann Y Resolution", 21, -21., 21., 100, -.5, .5), kTrackRYResolution);
815// fContainer->AddAt(new TH2I("fZRT", "Track Riemann Z Resolution", 21, -21., 21., 100, -.5, .5), kTrackRZResolution);
816// fContainer->AddAt(new TH2I("fPhiRT", "Track Riemann Angular Resolution", 21, -21., 21., 100, -10., 10.), kTrackRAngleResolution);
817//
818// Kalman track resolution [y, z, angular]
819// fContainer->AddAt(new TH2I("fYKT", "", 21, -21., 21., 100, -.5, .5), kTrackKYResolution);
820// fContainer->AddAt(new TH2I("fZKT", "", 21, -21., 21., 100, -.5, .5), kTrackKZResolution);
821// fContainer->AddAt(new TH2I("fPhiKT", "", 21, -21., 21., 100, -10., 10.), kTrackKAngleResolution);
822 }
3d86166d 823 return fContainer;
77203477 824}
825
aaf47b30 826
827//________________________________________________________
828void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
829{
3d86166d 830
aaf47b30 831 fReconstructor->SetRecoParam(r);
832}