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")
66 fReconstructor = new AliTRDReconstructor();
67 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
68 fGeo = new AliTRDgeometry();
71 //________________________________________________________
72 AliTRDtrackingResolution::~AliTRDtrackingResolution()
75 delete fReconstructor;
76 if(gGeoManager) delete gGeoManager;
80 //________________________________________________________
81 void AliTRDtrackingResolution::CreateOutputObjects()
84 OpenFile(0, "RECREATE");
86 fContainer = Histos();
88 // cluster to tracklet residuals [2]
89 fContainer->AddAt(new TH2I("fYClRes", "Clusters Residuals", 21, -21., 21., 100, -.5, .5), kClusterYResidual);
90 // tracklet to Riemann fit residuals [2]
91 fContainer->AddAt(new TH2I("fYTrkltRRes", "Tracklet Riemann Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanYResidual);
92 fContainer->AddAt(new TH2I("fAngleTrkltRRes", "Tracklet Riemann Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletRiemanAngleResidual);
93 fContainer->AddAt(new TH2I("fYTrkltKRes", "Tracklet Kalman Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanYResidual);
94 fContainer->AddAt(new TH2I("fAngleTrkltKRes", "Tracklet Kalman Angular Residuals", 21, -21., 21., 100, -.5, .5), kTrackletKalmanAngleResidual);
98 // tracklet resolution [0]
99 fContainer->AddAt(new TH2I("fY", "Tracklet Resolution", 21, -21., 21., 100, -.5, .5), kTrackletYResolution);
100 // tracklet angular resolution [1]
101 fContainer->AddAt(new TH2I("fPhi", "Tracklet Angular Resolution", 21, -21., 21., 100, -10., 10.), kTrackletAngleResolution);
103 // Riemann track resolution [y, z, angular]
104 fContainer->AddAt(new TH2I("fYRT", "Track Riemann Y Resolution", 21, -21., 21., 100, -.5, .5), kTrackRYResolution);
105 fContainer->AddAt(new TH2I("fZRT", "Track Riemann Z Resolution", 21, -21., 21., 100, -.5, .5), kTrackRZResolution);
106 fContainer->AddAt(new TH2I("fPhiRT", "Track Riemann Angular Resolution", 21, -21., 21., 100, -10., 10.), kTrackRAngleResolution);
108 // Kalman track resolution [y, z, angular]
109 fContainer->AddAt(new TH2I("fYKT", "", 21, -21., 21., 100, -.5, .5), kTrackKYResolution);
110 fContainer->AddAt(new TH2I("fZKT", "", 21, -21., 21., 100, -.5, .5), kTrackKZResolution);
111 fContainer->AddAt(new TH2I("fPhiKT", "", 21, -21., 21., 100, -10., 10.), kTrackKAngleResolution);
114 // CREATE GRAPHS for DISPLAY
116 // define iterator over graphs
117 Int_t jgraph = (Int_t)kGraphStart;
118 TH2I *h2 = (TH2I *)(fContainer->At(kClusterYResidual));
119 // clusters tracklet residuals (mean-phi)
120 TH1 *h = new TH1I("h", "", 100, -40., 40.);
121 h->GetXaxis()->SetTitle("#Phi [deg]");
122 h->GetYaxis()->SetTitle("Clusters Residuals : #sigma/#mu [mm]");
123 h->GetYaxis()->SetRangeUser(-.05, 1.);
124 fContainer->AddAt(h, jgraph++);
126 TGraphErrors *g = new TGraphErrors(h2->GetNbinsX());
127 g->SetLineColor(kGreen);
128 g->SetMarkerStyle(22);
129 g->SetMarkerColor(kGreen);
130 g->SetNameTitle("clm", "Residuals Clusters-Tracklet Mean");
131 fContainer->AddAt(g, jgraph++);
133 // clusters tracklet residuals (sigma-phi)
134 g = new TGraphErrors(h2->GetNbinsX());
135 g->SetLineColor(kRed);
136 g->SetMarkerStyle(23);
137 g->SetMarkerColor(kRed);
138 g->SetNameTitle("cls", "Residuals Clusters-Tracklet Sigma");
139 fContainer->AddAt(g, jgraph++);
142 // tracklet y resolution
143 h2 = (TH2I*)fContainer->At(kTrackletYResolution);
144 h = new TH1I("h", "", 100, -40., 40.);
145 h->GetXaxis()->SetTitle("#Phi [deg]");
146 h->GetYaxis()->SetTitle("Tracklet Resolution : #sigma/#mu [mm]");
147 h->GetYaxis()->SetRangeUser(-.05, 1.);
148 fContainer->AddAt(h, jgraph++);
150 g = new TGraphErrors(h2->GetNbinsX());
151 g->SetLineColor(kGreen);
152 g->SetMarkerStyle(22);
153 g->SetMarkerColor(kGreen);
154 g->SetNameTitle("trkltym", "Resolution Tracklet Y Mean");
155 fContainer->AddAt(g, jgraph++);
156 g = new TGraphErrors(h2->GetNbinsX());
157 g->SetLineColor(kRed);
158 g->SetMarkerStyle(22);
159 g->SetMarkerColor(kRed);
160 g->SetNameTitle("trkltys", "Resolution Tracklet Y Sigma");
161 fContainer->AddAt(g, jgraph++);
163 // tracklet phi resolution
164 h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
165 h = new TH1I("h", "", 100, -40., 40.);
166 h->GetXaxis()->SetTitle("#Phi [deg]");
167 h->GetYaxis()->SetTitle("Tracklet Angular Resolution : #sigma/#mu [deg]");
168 h->GetYaxis()->SetRangeUser(-.05, .2);
169 fContainer->AddAt(h, jgraph++);
171 g = new TGraphErrors(h2->GetNbinsX());
172 g->SetLineColor(kGreen);
173 g->SetMarkerStyle(22);
174 g->SetMarkerColor(kGreen);
175 g->SetNameTitle("trkltam", "Resolution Tracklet Y Mean");
176 fContainer->AddAt(g, jgraph++);
177 g = new TGraphErrors(h2->GetNbinsX());
178 g->SetLineColor(kRed);
179 g->SetMarkerStyle(22);
180 g->SetMarkerColor(kRed);
181 g->SetNameTitle("trkltas", "Angle Resolution Tracklet Sigma");
182 fContainer->AddAt(g, jgraph++);
186 //________________________________________________________
187 void AliTRDtrackingResolution::Exec(Option_t *)
189 // spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire})
190 // angular Resolution: res = Tracklet angle - TrackRef Angle
192 Int_t nTrackInfos = fTracks->GetEntriesFast();
194 printf("Event[%d] TrackInfos[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTrackInfos);
198 Double_t p, dy, dphi, dymc, dzmc, dphimc;
199 Float_t fP[kNLayers], fX[kNLayers], fY[kNLayers], fZ[kNLayers], fPhi[kNLayers], fTheta[kNLayers]; // phi/theta angle per layer
200 Bool_t fMCMap[kNLayers], fLayerMap[kNLayers]; // layer map
202 AliTRDpadPlane *pp = 0x0;
203 AliTrackPoint tr[kNLayers], tk[kNLayers];
204 AliExternalTrackParam *fOp = 0x0;
205 AliTRDtrackV1 *fTrack = 0x0;
206 AliTRDtrackInfo *fInfo = 0x0;
207 for(Int_t iTI = 0; iTI < nTrackInfos; iTI++){
208 // check if ESD and MC-Information are available
209 if(!(fInfo = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(iTI)))) continue;
210 if(!(fTrack = fInfo->GetTRDtrack())) continue;
211 if(!(fOp = fInfo->GetOuterParam())) continue;
212 pdg = fInfo->GetPDG();
214 if(fDebugLevel>=3) printf("\tDoing track[%d] NTrackRefs[%d]\n", iTI, fInfo->GetNTrackRefs());
218 memset(fP, 0, kNLayers*sizeof(Float_t));
219 memset(fX, 0, kNLayers*sizeof(Float_t));
220 memset(fY, 0, kNLayers*sizeof(Float_t));
221 memset(fZ, 0, kNLayers*sizeof(Float_t));
222 memset(fPhi, 0, kNLayers*sizeof(Float_t));
223 memset(fTheta, 0, kNLayers*sizeof(Float_t));
224 memset(fLayerMap, 0, kNLayers*sizeof(Bool_t));
225 memset(fMCMap, 0, kNLayers*sizeof(Bool_t));
226 AliTRDseedV1 *fTracklet = 0x0;
227 for(Int_t iplane = 0; iplane < kNLayers; iplane++){
228 if(!(fTracklet = fTrack->GetTracklet(iplane))) continue;
229 if(!fTracklet->IsOK()) continue;
232 fLayerMap[iplane] = kTRUE;
233 tr[npts].SetXYZ(fTracklet->GetX0(), 0., 0.);
234 tk[npts].SetXYZ(fTracklet->GetX0(), fTracklet->GetYfit(0), fTracklet->GetZfit(0));
237 if(fDebugLevel>=4) printf("\t\tLy[%d] X0[%6.3f] Ncl[%d]\n", iplane, fTracklet->GetX0(), fTracklet->GetN());
239 // define reference values
241 fX[iplane] = fTracklet->GetX0();
242 fY[iplane] = fTracklet->GetYref(0);
243 fZ[iplane] = fTracklet->GetZref(0);
244 fPhi[iplane] = TMath::ATan(fTracklet->GetYref(1));
245 fTheta[iplane] = TMath::ATan(fTracklet->GetZref(1));
248 // RESOLUTION (compare to MC)
250 if(fInfo->GetNTrackRefs() >= 2){
251 Double_t pmc, ymc, zmc, phiMC, thetaMC;
252 if(Resolution(fTracklet, fInfo, pmc, ymc, zmc, phiMC, thetaMC)){
253 fMCMap[iplane] = kTRUE;
257 fPhi[iplane] = phiMC;
258 fTheta[iplane] = thetaMC;
262 Float_t phi = fPhi[iplane]*TMath::RadToDeg();
263 //Float_t theta = fTheta[iplane]*TMath::RadToDeg();
265 // Do clusters residuals
266 if(!fTracklet->Fit(kFALSE)) continue;
267 AliTRDcluster *c = 0x0;
268 for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
269 if(!(c = fTracklet->GetClusters(ic))) continue;
271 dy = fTracklet->GetYat(c->GetX()) - c->GetY();
272 ((TH2I*)fContainer->At(kClusterYResidual))->Fill(phi, dy);
275 Float_t q = c->GetQ();
276 // Get z-position with respect to anode wire
277 AliTRDSimParam *simParam = AliTRDSimParam::Instance();
278 Int_t det = c->GetDetector();
279 Float_t x = c->GetX();
280 Float_t z = fZ[iplane]-(fX[iplane]-x)*TMath::Tan(fTheta[iplane]);
281 Int_t stack = fGeo->GetStack(det);
282 pp = fGeo->GetPadPlane(iplane, stack);
283 Float_t row0 = pp->GetRow0();
284 Float_t d = row0 - z + simParam->GetAnodeWireOffset();
285 d -= ((Int_t)(2 * d)) / 2.0;
286 //if (d > 0.25) d = 0.5 - d;
288 (*fDebugStream) << "ResidualClusters"
292 << "phi=" << fPhi[iplane]
293 << "tht=" << fTheta[iplane]
306 // this protection we might drop TODO
307 if(fTrack->GetNumberOfTracklets() < 6) continue;
309 AliTRDtrackerV1::FitRiemanTilt(fTrack, 0x0, kTRUE, npts, tr);
311 for(Int_t ip=0; ip<kNLayers; ip++){
312 if(!fLayerMap[ip]) continue;
313 fTracklet = fTrack->GetTracklet(ip);
314 // recalculate fit based on the new tilt correction
317 dy = fTracklet->GetYfit(0) - tr[iref].GetY();
318 ((TH2I*)fContainer->At(kTrackletRiemanYResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dy);
320 dphi = fTracklet->GetYfit(1)- fTracklet->GetYref(1);
321 ((TH2I*)fContainer->At(kTrackletRiemanAngleResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dphi);
324 dymc = fY[ip] - tr[iref].GetY();
325 ((TH2I*)fContainer->At(kTrackRYResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dymc);
327 dzmc = fZ[ip] - tr[iref].GetZ();
328 ((TH2I*)fContainer->At(kTrackRZResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dzmc);
330 dphimc = fPhi[ip] - fTracklet->GetYfit(1);
331 ((TH2I*)fContainer->At(kTrackRAngleResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dphimc);
337 (*fDebugStream) << "RiemannTrack"
339 << "mc=" << fMCMap[ip]
341 << "phi=" << fPhi[ip]
342 << "tht=" << fTheta[ip]
347 << "dphimc="<< dphimc
352 // if(!gGeoManager) TGeoManager::Import("geometry.root");
353 // AliTRDtrackerV1::FitKalman(fTrack, 0x0, kFALSE, nc, tr);
354 // for(Int_t ip=0; ip<nc; ip++){
355 // dy = cl[ip].GetY() - tr[ip].GetY();
356 // ((TH2I*)fContainer->At(kTrackletKalmanYResidual))->Fill(phi*TMath::RadToDeg(), dy);
357 // dz = cl[ip].GetZ() - tr[ip].GetZ();
358 // if(fDebugLevel>=1){
359 // (*fDebugStream) << "KalmanTrack"
362 // /* << "phi=" << phi
363 // << "theta=" << theta
364 // << "dphi=" << dphi*/
371 PostData(0, fContainer);
374 //________________________________________________________
375 void AliTRDtrackingResolution::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t *opt)
377 //sprintf(opt, "pl");
380 first = (Int_t)kGraphStart; last = first+3;
383 first = (Int_t)kGraphStart+3; last = first+3;
386 first = (Int_t)kGraphStart+6; last = first+3;
389 first = (Int_t)kGraphStart; last = first;
395 //________________________________________________________
396 Bool_t AliTRDtrackingResolution::Resolution(AliTRDseedV1 *tracklet, AliTRDtrackInfo *fInfo, Double_t &p, Double_t &ymc, Double_t &zmc, Double_t &phi, Double_t &theta)
399 AliTrackReference *fTrackRefs[2] = {0x0, 0x0}, *tempTrackRef = 0x0;
401 // check for 2 track ref where the radial position has a distance less than 3.7mm
403 for(Int_t itr = 0; itr < fInfo->GetNTrackRefs(); itr++){
404 if(!(tempTrackRef = fInfo->GetTrackRef(itr))) continue;
405 if(fDebugLevel>=5) printf("TrackRef %d: x = %f\n", itr, tempTrackRef->LocalX());
406 if(TMath::Abs(tracklet->GetX0() - tempTrackRef->LocalX()) > 3.7) continue;
407 fTrackRefs[nFound++] = tempTrackRef;
408 if(nFound == 2) break;
411 if(fDebugLevel>=4) printf("\t\tFound track crossing [%d] refX[%6.3f]\n", nFound, nFound>0 ? fTrackRefs[0]->LocalX() : 0.);
414 // We found 2 track refs for the tracklet, get y and z at the anode wire by a linear approximation
418 Double_t dx = fTrackRefs[1]->LocalX() - fTrackRefs[0]->LocalX();
420 if(fDebugLevel>=4) printf("\t\ttrack ref in the wrong order refX0[%6.3f] refX1[%6.3f]\n", fTrackRefs[0]->LocalX(), fTrackRefs[1]->LocalX());
423 Double_t dydx = (fTrackRefs[1]->LocalY() - fTrackRefs[0]->LocalY()) / dx;
424 Double_t dzdx = (fTrackRefs[1]->Z() - fTrackRefs[0]->Z()) / dx;
425 Double_t dx0 = fTrackRefs[1]->LocalX() - tracklet->GetX0();
426 ymc = fTrackRefs[1]->LocalY() - dydx*dx0;
427 zmc = fTrackRefs[1]->Z() - dzdx*dx0;
429 // recalculate tracklet based on the MC info
430 tracklet->SetZref(0, zmc);
431 tracklet->SetZref(1, -dzdx); // TODO
433 Double_t dy = tracklet->GetYfit(0) - ymc;
434 Double_t dz = tracklet->GetZfit(0) - zmc;
436 //res_y *= 100; // in mm
437 p = fTrackRefs[0]->P();
439 phi = TMath::ATan(dydx);
440 theta = TMath::ATan(dzdx);
441 Double_t dphi = TMath::ATan(tracklet->GetYfit(1)) - phi;
442 if(fDebugLevel>=4) printf("\t\tdx[%6.4f] dy[%6.4f] dz[%6.4f] dphi[%6.4f] \n", dx, dy, dz, dphi);
445 if(TMath::Abs(dx-3.7)<1.E-3){
446 ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(phi*TMath::RadToDeg(), dy);
447 ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(phi*TMath::RadToDeg(), dphi*TMath::RadToDeg());
451 Int_t iplane = tracklet->GetPlane();
452 Int_t pdg = fInfo->GetPDG();
453 (*fDebugStream) << "ResolutionTrklt"
468 AliTRDpadPlane *pp = 0x0;
469 AliTRDcluster *c = 0x0;
470 for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
471 if(!(c = tracklet->GetClusters(ic))) continue;
473 pp = fGeo->GetPadPlane(iplane, fGeo->GetStack(c->GetDetector()));
474 z0 = pp->GetRow0() + AliTRDSimParam::Instance()->GetAnodeWireOffset();
476 dx = tracklet->GetX0() - c->GetX();
477 Float_t yt = ymc - dx*dydx;
478 Float_t zt = zmc - dx*dzdx;
479 Float_t yc = c->GetY() - TMath::Tan(tracklet->GetTilt()) * (c->GetZ() - zt);
482 d -= ((Int_t)(2 * d)) / 2.0;
483 (*fDebugStream) << "ResolutionClstr"
498 //________________________________________________________
499 Bool_t AliTRDtrackingResolution::PostProcess()
501 //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
504 Printf("ERROR: list not available");
510 TGraphErrors *gm = 0x0, *gs = 0x0;
511 TF1 f("f1", "gaus", -.5, .5);
512 // define iterator over graphs
513 Int_t jgraph = (Int_t)kGraphStart;
515 //PROCESS RESIDUAL DISTRIBUTIONS
517 // Clusters residuals
519 TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
520 h2 = (TH2I *)(fContainer->At(kClusterYResidual));
521 jgraph++; //skip the frame histo
522 gm = (TGraphErrors*)fContainer->At(jgraph++);
523 gs = (TGraphErrors*)fContainer->At(jgraph++);
524 for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
525 Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
526 Double_t dphi = h2->GetXaxis()->GetBinWidth(ibin)/2;
527 h = h2->ProjectionY("py", ibin, ibin);
529 gm->SetPoint(ibin - 1, phi, 10.*fc.GetParameter(1));
530 gm->SetPointError(ibin - 1, dphi, 10.*fc.GetParError(1));
531 gs->SetPoint(ibin - 1, phi, 10.*fc.GetParameter(2));
532 gs->SetPointError(ibin - 1, dphi, 10.*fc.GetParError(2));
537 //PROCESS RESOLUTION DISTRIBUTIONS
539 // tracklet y resolution
540 h2 = (TH2I*)fContainer->At(kTrackletYResolution);
541 jgraph++; //skip the frame histo
542 gm = (TGraphErrors*)fContainer->At(jgraph++);
543 gs = (TGraphErrors*)fContainer->At(jgraph++);
544 for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
545 Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
546 f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
547 h = h2->ProjectionY("py", iphi, iphi);
549 Int_t jphi = iphi -1;
550 gm->SetPoint(jphi, phi, 10.*f.GetParameter(1));
551 gm->SetPointError(jphi, 0., 10.*f.GetParError(1));
552 gs->SetPoint(jphi, phi, 10.*f.GetParameter(2));
553 gs->SetPointError(jphi, 0., 10.*f.GetParError(2));
557 // tracklet phi resolution
558 h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
559 jgraph++; //skip the frame histo
560 gm = (TGraphErrors*)fContainer->At(jgraph++);
561 gs = (TGraphErrors*)fContainer->At(jgraph++);
562 for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
563 Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
564 f.SetParameter(1, 0.);f.SetParameter(2, 2.e-2);
565 h = h2->ProjectionY("py", iphi, iphi);
566 h->Fit(&f, "QN", "", -.5, .5);
567 Int_t jphi = iphi -1;
568 gm->SetPoint(jphi, phi, f.GetParameter(1));
569 gm->SetPointError(jphi, 0., f.GetParError(1));
570 gs->SetPoint(jphi, phi, f.GetParameter(2));
571 gs->SetPointError(jphi, 0., f.GetParError(2));
580 //________________________________________________________
581 void AliTRDtrackingResolution::Terminate(Option_t *)
588 if(HasPostProcess()) PostProcess();
591 //________________________________________________________
592 void AliTRDtrackingResolution::Fit(TH1 *h, TF1 *f)
594 // Helper function to avoid duplication of code
595 // Make first guesses on the fit parameters
597 // find the intial parameters of the fit !! (thanks George)
598 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
600 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
601 f->SetParLimits(0, 0., 3.*sum);
602 f->SetParameter(0, .9*sum);
604 f->SetParLimits(1, -.1, .1);
605 f->SetParameter(1, 0.);
607 f->SetParLimits(2, 0., 1.e-1);
608 f->SetParameter(2, 2.e-2);
610 f->SetParLimits(3, 0., sum);
611 f->SetParameter(3, .1*sum);
613 f->SetParLimits(4, -.3, .3);
614 f->SetParameter(4, 0.);
616 f->SetParLimits(5, 0., 1.e2);
617 f->SetParameter(5, 2.e-1);
619 h->Fit(f, "QN", "", -0.5, 0.5);
622 //________________________________________________________
623 TObjArray* AliTRDtrackingResolution::Histos()
625 if(!fContainer) fContainer = new TObjArray(25);
630 //________________________________________________________
631 void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
634 fReconstructor->SetRecoParam(r);