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-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 **************************************************************************/
16 /* $Id: AliTRDtrackingResolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
18 ////////////////////////////////////////////////////////////////////////////
20 // Reconstruction QA //
23 // Markus Fasel <M.Fasel@gsi.de> //
25 ////////////////////////////////////////////////////////////////////////////
30 #include <TObjArray.h>
35 #include <TGraphErrors.h>
37 #include "TTreeStream.h"
38 #include "TGeoManager.h"
40 #include "AliAnalysisManager.h"
41 #include "AliTrackReference.h"
42 #include "AliTrackPointArray.h"
43 #include "AliCDBManager.h"
45 #include "AliTRDSimParam.h"
46 #include "AliTRDgeometry.h"
47 #include "AliTRDpadPlane.h"
48 #include "AliTRDcluster.h"
49 #include "AliTRDseedV1.h"
50 #include "AliTRDtrackV1.h"
51 #include "AliTRDtrackerV1.h"
52 #include "AliTRDReconstructor.h"
53 #include "AliTRDrecoParam.h"
55 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
56 #include "AliTRDtrackingResolution.h"
58 ClassImp(AliTRDtrackingResolution)
60 //________________________________________________________
61 AliTRDtrackingResolution::AliTRDtrackingResolution()
62 :AliTRDrecoTask("Resolution", "Tracking Resolution")
68 fReconstructor = new AliTRDReconstructor();
69 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
70 fGeo = new AliTRDgeometry();
73 //________________________________________________________
74 AliTRDtrackingResolution::~AliTRDtrackingResolution()
76 if(fGraphS){fGraphS->Delete(); delete fGraphS;}
77 if(fGraphM){fGraphM->Delete(); delete fGraphM;}
79 delete fReconstructor;
80 if(gGeoManager) delete gGeoManager;
84 //________________________________________________________
85 void AliTRDtrackingResolution::CreateOutputObjects()
88 OpenFile(0, "RECREATE");
90 fContainer = Histos();
92 // cluster to tracklet residuals [2]
93 fContainer->AddAt(new TH2I("fYClRes", "Clusters Residuals", 21, -21., 21., 100, -.5, .5), kClusterYResidual);
94 // // tracklet to Riemann fit residuals [2]
95 // fContainer->AddAt(new TH2I("fYTrkltRRes", "Tracklet Riemann Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanYResidual);
96 // fContainer->AddAt(new TH2I("fAngleTrkltRRes", "Tracklet Riemann Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanAngleResidual);
97 // fContainer->AddAt(new TH2I("fYTrkltKRes", "Tracklet Kalman Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanYResidual);
98 // fContainer->AddAt(new TH2I("fAngleTrkltKRes", "Tracklet Kalman Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanAngleResidual);
102 // cluster y resolution [0]
103 fContainer->AddAt(new TH2I("fCY", "Cluster Resolution", 31, -31., 31., 100, -.5, .5), kClusterYResolution);
104 // tracklet y resolution [0]
105 fContainer->AddAt(new TH2I("fY", "Tracklet Resolution", 31, -31., 31., 100, -.5, .5), kTrackletYResolution);
106 // tracklet angular resolution [1]
107 fContainer->AddAt(new TH2I("fPhi", "Tracklet Angular Resolution", 31, -31., 31., 100, -10., 10.), kTrackletAngleResolution);
109 // // Riemann track resolution [y, z, angular]
110 // fContainer->AddAt(new TH2I("fYRT", "Track Riemann Y Resolution", 21, -21., 21., 100, -.5, .5), kTrackRYResolution);
111 // fContainer->AddAt(new TH2I("fZRT", "Track Riemann Z Resolution", 21, -21., 21., 100, -.5, .5), kTrackRZResolution);
112 // fContainer->AddAt(new TH2I("fPhiRT", "Track Riemann Angular Resolution", 21, -21., 21., 100, -10., 10.), kTrackRAngleResolution);
114 // Kalman track resolution [y, z, angular]
115 // fContainer->AddAt(new TH2I("fYKT", "", 21, -21., 21., 100, -.5, .5), kTrackKYResolution);
116 // fContainer->AddAt(new TH2I("fZKT", "", 21, -21., 21., 100, -.5, .5), kTrackKZResolution);
117 // fContainer->AddAt(new TH2I("fPhiKT", "", 21, -21., 21., 100, -10., 10.), kTrackKAngleResolution);
121 //________________________________________________________
122 void AliTRDtrackingResolution::Exec(Option_t *)
124 // spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire})
125 // angular Resolution: res = Tracklet angle - TrackRef Angle
127 Int_t nTrackInfos = fTracks->GetEntriesFast();
128 if(fDebugLevel>=2 && nTrackInfos){
129 printf("Event[%d] TrackInfos[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTrackInfos);
133 Double_t p, dy/*, dphi, dymc, dzmc, dphimc*/;
134 Float_t fP[kNLayers], fX[kNLayers], fY[kNLayers], fZ[kNLayers], fPhi[kNLayers], fTheta[kNLayers]; // phi/theta angle per layer
135 Bool_t fMCMap[kNLayers], fLayerMap[kNLayers]; // layer map
137 AliTRDpadPlane *pp = 0x0;
138 AliTrackPoint tr[kNLayers], tk[kNLayers];
139 AliExternalTrackParam *fOp = 0x0;
140 AliTRDtrackV1 *fTrack = 0x0;
141 AliTRDtrackInfo *fInfo = 0x0;
142 for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){
143 // check if ESD and MC-Information are available
144 if(!(fInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iTI)))) continue;
145 if(!(fTrack = fInfo->GetTrack())) continue;
146 if(!(fOp = fInfo->GetOuterParam())) continue;
147 pdg = fInfo->GetPDG();
149 if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs());
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;
167 fLayerMap[iplane] = kTRUE;
168 tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.);
169 tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0));
172 if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN());
174 // define reference values
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));
183 // RESOLUTION (compare to MC)
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;
192 fPhi[iplane] = phiMC;
193 fTheta[iplane] = thetaMC;
197 Float_t phi = fPhi[iplane]*TMath::RadToDeg();
198 //Float_t theta = fTheta[iplane]*TMath::RadToDeg();
200 // Do clusters residuals
201 if(!fTracklet->Fit(kFALSE)) continue;
202 AliTRDcluster *c = 0x0;
203 for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
204 if(!(c = fTracklet->GetClusters(ic))) continue;
206 dy = fTracklet->GetYat(c->GetX()) - c->GetY();
207 ((TH2I*)fContainer->At(kClusterYResidual))->Fill(phi, dy);
210 Float_t q = c->GetQ();
211 // Get z-position with respect to anode wire
212 AliTRDSimParam *simParam = AliTRDSimParam::Instance();
213 Int_t det = c->GetDetector();
214 Float_t x = c->GetX();
215 Float_t z = fZ[iplane]-(fX[iplane]-x)*TMath::Tan(fTheta[iplane]);
216 Int_t stack = fGeo->GetStack(det);
217 pp = fGeo->GetPadPlane(iplane, stack);
218 Float_t row0 = pp->GetRow0();
219 Float_t d = row0 - z + simParam->GetAnodeWireOffset();
220 d -= ((Int_t)(2 * d)) / 2.0;
221 //if (d > 0.25) d = 0.5 - d;
223 (*fDebugStream) << "ResidualClusters"
227 << "phi=" << fPhi[iplane]
228 << "tht=" << fTheta[iplane]
241 // // this protection we might drop TODO
242 // if(fTrack->GetNumberOfTracklets() < 6) continue;
244 // AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr);
246 // for(Int_t ip=0; ip<kNLayers; ip++){
247 // if(!fLayerMap[ip]) continue;
248 // fTracklet = fTrack->GetTracklet(ip);
249 // // recalculate fit based on the new tilt correction
252 // dy = fTracklet->GetYfit(0) - tr[iref].GetY();
253 // ((TH2I*)fContainer->At(kTrackletRiemanYResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dy);
255 // dphi = fTracklet->GetYfit(1)- fTracklet->GetYref(1);
256 // ((TH2I*)fContainer->At(kTrackletRiemanAngleResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dphi);
259 // dymc = fY[ip] - tr[iref].GetY();
260 // ((TH2I*)fContainer->At(kTrackRYResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dymc);
262 // dzmc = fZ[ip] - tr[iref].GetZ();
263 // ((TH2I*)fContainer->At(kTrackRZResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dzmc);
265 // dphimc = fPhi[ip] - fTracklet->GetYfit(1);
266 // ((TH2I*)fContainer->At(kTrackRAngleResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dphimc);
271 // if(fDebugLevel>=1){
272 // (*fDebugStream) << "RiemannTrack"
274 // << "mc=" << fMCMap[ip]
276 // << "phi=" << fPhi[ip]
277 // << "tht=" << fTheta[ip]
279 // << "dphi=" << dphi
280 // << "dymc=" << dymc
281 // << "dzmc=" << dzmc
282 // << "dphimc="<< dphimc
287 // if(!gGeoManager) TGeoManager::Import("geometry.root");
288 // AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr);
289 // for(Int_t ip=0; ip<nc; ip++){
290 // dy = cl[ip].GetY() - tr[ip].GetY();
291 // ((TH2I*)fContainer->At(kTrackletKalmanYResidual))->Fill(phi*TMath::RadToDeg(), dy);
292 // dz = cl[ip].GetZ() - tr[ip].GetZ();
293 // if(fDebugLevel>=1){
294 // (*fDebugStream) << "KalmanTrack"
297 // /* << "phi=" << phi
298 // << "theta=" << theta
299 // << "dphi=" << dphi*/
306 PostData(0, fContainer);
309 //________________________________________________________
310 void AliTRDtrackingResolution::GetRefFigure(Int_t ifig)
313 TGraphErrors *g = 0x0;
315 case kClusterYResidual:
316 g = (TGraphErrors*)fGraphS->At(kClusterYResidual);
318 ax = g->GetHistogram()->GetYaxis();
319 ax->SetRangeUser(-.5, 1.);
320 ax->SetTitle("Clusters Y Residuals #sigma/#mu [mm]");
321 ax = g->GetHistogram()->GetXaxis();
322 ax->SetTitle("#phi [deg]");
323 ((TGraphErrors*)fGraphM->At(kClusterYResidual))->Draw("pl");
325 case kClusterYResolution:
326 g = (TGraphErrors*)fGraphS->At(kClusterYResolution);
327 ax = g->GetHistogram()->GetYaxis();
328 ax->SetRangeUser(-.5, 1.);
329 ax->SetTitle("Cluster Y Resolution #sigma/#mu [mm]");
330 ax = g->GetHistogram()->GetXaxis();
331 ax->SetTitle("#phi [deg]");
333 ((TGraphErrors*)fGraphM->At(kClusterYResolution))->Draw("pl");
335 case kTrackletYResolution:
336 g = (TGraphErrors*)fGraphS->At(kTrackletYResolution);
337 ax = g->GetHistogram()->GetYaxis();
338 ax->SetRangeUser(-.5, 1.);
339 ax->SetTitle("Tracklet Y Resolution #sigma/#mu [mm]");
340 ax = g->GetHistogram()->GetXaxis();
341 ax->SetTitle("#phi [deg]");
343 ((TGraphErrors*)fGraphM->At(kTrackletYResolution))->Draw("pl");
345 case kTrackletAngleResolution:
346 g = (TGraphErrors*)fGraphS->At(kTrackletAngleResolution);
347 ax = g->GetHistogram()->GetYaxis();
348 ax->SetRangeUser(-.1, 1.);
349 ax->SetTitle("Tracklet Angular Resolution #sigma/#mu [mrad]");
350 ax = g->GetHistogram()->GetXaxis();
351 ax->SetTitle("#phi [deg]");
353 ((TGraphErrors*)fGraphM->At(kTrackletAngleResolution))->Draw("pl");
356 AliInfo(Form("Reference plot [%d] not implemented yet", ifig));
362 //________________________________________________________
363 Bool_t AliTRDtrackingResolution::Resolution(AliTRDseedV1 *tracklet, AliTRDtrackInfo *fInfo, Double_t &p, Double_t &ymc, Double_t &zmc, Double_t &phi, Double_t &theta)
366 AliTrackReference *fTrackRefs[2] = {0x0, 0x0};
367 Float_t x0 = tracklet->GetX0();
368 Float_t tilt= tracklet->GetTilt();
369 Int_t cross = tracklet->GetNChange();
370 Int_t det = tracklet->GetDetector();
371 Int_t pdg = fInfo->GetPDG();
373 // check for 2 track ref where the radial position has a distance less than 3.7mm
375 for(Int_t itr = 0; itr < AliTRDtrackInfo::kNTrackRefs; itr++){
376 if(!(fTrackRefs[nFound] = fInfo->GetTrackRef(itr))) break;
378 if(fDebugLevel>=5) printf("\t\tref[%2d] x[%6.3f]\n", itr, fTrackRefs[nFound]->LocalX());
379 if(TMath::Abs(x0 - fTrackRefs[nFound]->LocalX()) > 3.7) continue;
381 if(nFound == 2) break;
384 if(fDebugLevel>=3) printf("\t\tMissing track ref x0[%6.3f] ly[%d] nref[%d]\n", x0, tracklet->GetPlane(), fInfo->GetMCinfo()->GetNTrackRefs());
387 // We found 2 track refs for the tracklet, get y and z at the anode wire by a linear approximation
391 Double_t dx = fTrackRefs[1]->LocalX() - fTrackRefs[0]->LocalX();
392 if(dx <= 0. || TMath::Abs(dx-3.7)>1.E-3){
393 if(fDebugLevel>=3) printf("\t\tTrack ref with wrong radial distances refX0[%6.3f] refX1[%6.3f]\n", fTrackRefs[0]->LocalX(), fTrackRefs[1]->LocalX());
397 Double_t dydx = (fTrackRefs[1]->LocalY() - fTrackRefs[0]->LocalY()) / dx;
398 Double_t dzdx = (fTrackRefs[1]->Z() - fTrackRefs[0]->Z()) / dx;
399 Double_t dx0 = fTrackRefs[1]->LocalX() - tracklet->GetX0();
400 ymc = fTrackRefs[1]->LocalY() - dydx*dx0;
401 zmc = fTrackRefs[1]->Z() - dzdx*dx0;
403 // recalculate tracklet based on the MC info
404 AliTRDseedV1 tt(*tracklet);
408 Double_t dy = tt.GetYfit(0) - ymc;
411 Double_t *xyz = tt.GetCrossXYZ();
412 dz = xyz[2] - (zmc - (x0 - xyz[0])*dzdx) ;
414 p = fTrackRefs[0]->P();
416 phi = TMath::ATan(dydx);
417 theta = TMath::ATan(dzdx);
418 Double_t dphi = TMath::ATan(tt.GetYfit(1)) - phi;
419 if(fDebugLevel>=4) printf("\t\tdx[%6.4f] dy[%6.4f] dz[%6.4f] dphi[%6.4f] \n", dx, dy, dz, dphi);
422 ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(phi*TMath::RadToDeg(), dy);
423 ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(phi*TMath::RadToDeg(), dphi*TMath::RadToDeg());
427 (*fDebugStream) << "ResolutionTrklt"
442 AliTRDpadPlane *pp = fGeo->GetPadPlane(AliTRDgeometry::GetLayer(det), AliTRDgeometry::GetStack(det));
443 Float_t z0 = pp->GetRow0() + AliTRDSimParam::Instance()->GetAnodeWireOffset();
445 AliTRDcluster *c = 0x0;
446 tracklet->ResetClusterIter(kFALSE);
447 while((c = tracklet->PrevCluster())){
448 Float_t q = TMath::Abs(c->GetQ());
449 Float_t xc = c->GetX();
450 Float_t yc = c->GetY();
451 Float_t zc = c->GetZ();
453 Float_t yt = ymc - dx*dydx;
454 Float_t zt = zmc - dx*dzdx;
455 dy = yt - (yc - tilt*(zc-zt));
458 if(q>100.) ((TH2I*)fContainer->At(kClusterYResolution))->Fill(phi*TMath::RadToDeg(), dy);
463 d -= ((Int_t)(2 * d)) / 2.0;
464 (*fDebugStream) << "ResolutionClstr"
480 //________________________________________________________
481 Bool_t AliTRDtrackingResolution::PostProcess()
483 //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
485 Printf("ERROR: list not available");
488 fNRefFigures = fContainer->GetEntriesFast();
490 fGraphS = new TObjArray(fNRefFigures);
494 fGraphM = new TObjArray(fNRefFigures);
500 TGraphErrors *gm = 0x0, *gs = 0x0;
503 TF1 f("f1", "gaus", -.5, .5);
505 TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
507 //PROCESS RESIDUAL DISTRIBUTIONS
509 // Clusters residuals
510 h2 = (TH2I *)(fContainer->At(kClusterYResidual));
511 gm = new TGraphErrors(h2->GetNbinsX());
512 gm->SetLineColor(kRed);
513 gm->SetMarkerStyle(23);
514 gm->SetMarkerColor(kRed);
515 gm->SetNameTitle("clm", "");
516 fGraphM->AddAt(gm, kClusterYResidual);
517 gs = new TGraphErrors(h2->GetNbinsX());
518 gs->SetLineColor(kRed);
519 gs->SetMarkerStyle(23);
520 gs->SetMarkerColor(kRed);
521 gs->SetNameTitle("cls", "");
522 fGraphS->AddAt(gs, kClusterYResidual);
523 for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
524 Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
525 Double_t dphi = h2->GetXaxis()->GetBinWidth(ibin)/2;
526 h = h2->ProjectionY("py", ibin, ibin);
528 gm->SetPoint(ibin - 1, phi, 10.*fc.GetParameter(1));
529 gm->SetPointError(ibin - 1, dphi, 10.*fc.GetParError(1));
530 gs->SetPoint(ibin - 1, phi, 10.*fc.GetParameter(2));
531 gs->SetPointError(ibin - 1, dphi, 10.*fc.GetParError(2));
535 //PROCESS RESOLUTION DISTRIBUTIONS
538 // cluster y resolution
539 h2 = (TH2I*)fContainer->At(kClusterYResolution);
540 gm = new TGraphErrors(h2->GetNbinsX());
541 gm->SetLineColor(kRed);
542 gm->SetMarkerStyle(23);
543 gm->SetMarkerColor(kRed);
544 gm->SetNameTitle("clym", "");
545 fGraphM->AddAt(gm, kClusterYResolution);
546 gs = new TGraphErrors(h2->GetNbinsX());
547 gs->SetLineColor(kRed);
548 gs->SetMarkerStyle(23);
549 gs->SetMarkerColor(kRed);
550 gs->SetNameTitle("clys", "");
551 fGraphS->AddAt(gs, kClusterYResolution);
552 for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
553 h = h2->ProjectionY("py", iphi, iphi);
555 Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
556 Int_t jphi = iphi -1;
557 gm->SetPoint(jphi, phi, 10.*f.GetParameter(1));
558 gm->SetPointError(jphi, 0., 10.*f.GetParError(1));
559 gs->SetPoint(jphi, phi, 10.*f.GetParameter(2));
560 gs->SetPointError(jphi, 0., 10.*f.GetParError(2));
563 // tracklet y resolution
564 h2 = (TH2I*)fContainer->At(kTrackletYResolution);
565 gm = new TGraphErrors(h2->GetNbinsX());
566 gm->SetLineColor(kRed);
567 gm->SetMarkerStyle(23);
568 gm->SetMarkerColor(kRed);
569 gm->SetNameTitle("trkltym", "");
570 fGraphM->AddAt(gm, kTrackletYResolution);
571 gs = new TGraphErrors(h2->GetNbinsX());
572 gs->SetLineColor(kRed);
573 gs->SetMarkerStyle(23);
574 gs->SetMarkerColor(kRed);
575 gs->SetNameTitle("trkltys", "");
576 fGraphS->AddAt(gs, kTrackletYResolution);
577 for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
578 h = h2->ProjectionY("py", iphi, iphi);
580 Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
581 Int_t jphi = iphi -1;
582 gm->SetPoint(jphi, phi, 10.*f.GetParameter(1));
583 gm->SetPointError(jphi, 0., 10.*f.GetParError(1));
584 gs->SetPoint(jphi, phi, 10.*f.GetParameter(2));
585 gs->SetPointError(jphi, 0., 10.*f.GetParError(2));
588 // tracklet phi resolution
589 h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
590 gm = new TGraphErrors(h2->GetNbinsX());
591 gm->SetLineColor(kRed);
592 gm->SetMarkerStyle(23);
593 gm->SetMarkerColor(kRed);
594 gm->SetNameTitle("trkltym", "");
595 fGraphM->AddAt(gm, kTrackletAngleResolution);
596 gs = new TGraphErrors(h2->GetNbinsX());
597 gs->SetLineColor(kRed);
598 gs->SetMarkerStyle(23);
599 gs->SetMarkerColor(kRed);
600 gs->SetNameTitle("trkltys", "");
601 fGraphS->AddAt(gs, kTrackletAngleResolution);
602 for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
603 h = h2->ProjectionY("py", iphi, iphi);
604 h->Fit(&f, "QN", "", -.5, .5);
605 Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
606 Int_t jphi = iphi -1;
607 gm->SetPoint(jphi, phi, f.GetParameter(1));
608 gm->SetPointError(jphi, 0., f.GetParError(1));
609 gs->SetPoint(jphi, phi, f.GetParameter(2));
610 gs->SetPointError(jphi, 0., f.GetParError(2));
618 //________________________________________________________
619 void AliTRDtrackingResolution::Terminate(Option_t *)
626 if(HasPostProcess()) PostProcess();
629 //________________________________________________________
630 void AliTRDtrackingResolution::Fit(TH1 *h, TF1 *f)
632 // Helper function to avoid duplication of code
633 // Make first guesses on the fit parameters
635 // find the intial parameters of the fit !! (thanks George)
636 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
638 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
639 f->SetParLimits(0, 0., 3.*sum);
640 f->SetParameter(0, .9*sum);
642 f->SetParLimits(1, -.1, .1);
643 f->SetParameter(1, 0.);
645 f->SetParLimits(2, 0., 1.e-1);
646 f->SetParameter(2, 2.e-2);
648 f->SetParLimits(3, 0., sum);
649 f->SetParameter(3, .1*sum);
651 f->SetParLimits(4, -.3, .3);
652 f->SetParameter(4, 0.);
654 f->SetParLimits(5, 0., 1.e2);
655 f->SetParameter(5, 2.e-1);
657 h->Fit(f, "QN", "", -0.5, 0.5);
660 //________________________________________________________
661 TObjArray* AliTRDtrackingResolution::Histos()
663 if(!fContainer) fContainer = new TObjArray(4);
668 //________________________________________________________
669 void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
672 fReconstructor->SetRecoParam(r);