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 // TRD tracking resolution //
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
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
32 // For stand alone usage please refer to the following example:
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);
46 // Alexandru Bercuci <A.Bercuci@gsi.de> //
47 // Markus Fasel <M.Fasel@gsi.de> //
49 ////////////////////////////////////////////////////////////////////////////
54 #include <TObjArray.h>
60 #include <TGraphErrors.h>
62 #include "TTreeStream.h"
63 #include "TGeoManager.h"
65 #include "AliAnalysisManager.h"
66 #include "AliTrackReference.h"
67 #include "AliTrackPointArray.h"
68 #include "AliCDBManager.h"
70 #include "AliTRDSimParam.h"
71 #include "AliTRDgeometry.h"
72 #include "AliTRDpadPlane.h"
73 #include "AliTRDcluster.h"
74 #include "AliTRDseedV1.h"
75 #include "AliTRDtrackV1.h"
76 #include "AliTRDtrackerV1.h"
77 #include "AliTRDReconstructor.h"
78 #include "AliTRDrecoParam.h"
80 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
81 #include "AliTRDtrackingResolution.h"
83 ClassImp(AliTRDtrackingResolution)
85 //________________________________________________________
86 AliTRDtrackingResolution::AliTRDtrackingResolution()
87 :AliTRDrecoTask("Resolution", "Tracking Resolution")
94 fReconstructor = new AliTRDReconstructor();
95 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
96 fGeo = new AliTRDgeometry();
100 //________________________________________________________
101 AliTRDtrackingResolution::~AliTRDtrackingResolution()
103 if(fGraphS){fGraphS->Delete(); delete fGraphS;}
104 if(fGraphM){fGraphM->Delete(); delete fGraphM;}
106 delete fReconstructor;
107 if(gGeoManager) delete gGeoManager;
111 //________________________________________________________
112 void AliTRDtrackingResolution::CreateOutputObjects()
114 // spatial resolution
115 OpenFile(0, "RECREATE");
117 fContainer = Histos();
120 // //________________________________________________________
121 // void AliTRDtrackingResolution::Exec(Option_t *)
123 // // spatial Resolution: res = pos_{Tracklet}(x = x_{Anode wire}) - pos_{TrackRef}(x = x_{Anode wire})
124 // // angular Resolution: res = Tracklet angle - TrackRef Angle
126 // Int_t nTrackInfos = fTracks->GetEntriesFast();
127 // if(fDebugLevel>=2 && nTrackInfos){
128 // printf("Event[%d] TrackInfos[%d]\n", (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTrackInfos);
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
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.;
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;
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));
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 // theta += fTheta[iplane]; nly++;
199 // if(fTracklet->GetNChange()) ncrs++;
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;
207 // dy = fTracklet->GetYat(c->GetX()) - c->GetY();
208 // ((TH2I*)fContainer->At(kClusterYResidual))->Fill(phi, dy);
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;
224 // (*fDebugStream) << "ResidualClusters"
225 // << "ly=" << iplane
226 // << "stk=" << stack
228 // << "phi=" << fPhi[iplane]
229 // << "tht=" << fTheta[iplane]
240 // if(nly) theta /= nly;
241 // if(fDebugLevel>=1){
242 // (*fDebugStream) << "TrackStatistics"
244 // << "ncrs=" << ncrs
245 // << "tht=" << theta
250 // // // this protection we might drop TODO
251 // // if(fTrack->GetNumberOfTracklets() < 6) continue;
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();
261 // // dy = fTracklet->GetYfit(0) - tr[iref].GetY();
262 // // ((TH2I*)fContainer->At(kTrackletRiemanYResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dy);
264 // // dphi = fTracklet->GetYfit(1)- fTracklet->GetYref(1);
265 // // ((TH2I*)fContainer->At(kTrackletRiemanAngleResidual))->Fill(fPhi[ip]*TMath::RadToDeg(), dphi);
267 // // if(HasMCdata()){
268 // // dymc = fY[ip] - tr[iref].GetY();
269 // // ((TH2I*)fContainer->At(kTrackRYResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dymc);
271 // // dzmc = fZ[ip] - tr[iref].GetZ();
272 // // ((TH2I*)fContainer->At(kTrackRZResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dzmc);
274 // // dphimc = fPhi[ip] - fTracklet->GetYfit(1);
275 // // ((TH2I*)fContainer->At(kTrackRAngleResolution))->Fill(fPhi[ip]*TMath::RadToDeg(), dphimc);
280 // // if(fDebugLevel>=1){
281 // // (*fDebugStream) << "RiemannTrack"
283 // // << "mc=" << fMCMap[ip]
284 // // << "p=" << fP[ip]
285 // // << "phi=" << fPhi[ip]
286 // // << "tht=" << fTheta[ip]
288 // // << "dphi=" << dphi
289 // // << "dymc=" << dymc
290 // // << "dzmc=" << dzmc
291 // // << "dphimc="<< dphimc
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"
306 // // /* << "phi=" << phi
307 // // << "theta=" << theta
308 // // << "dphi=" << dphi*/
315 // PostData(0, fContainer);
318 //________________________________________________________
319 TH1* AliTRDtrackingResolution::PlotClusterResiduals(const AliTRDtrackV1 *track)
321 if(track) fTrack = track;
323 AliWarning("No Track defined.");
327 if(!(h = ((TH2I*)fContainer->At(kClusterResidual)))){
328 AliWarning("No output histogram defined.");
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;
338 x0 = fTracklet->GetX0();
340 // retrive the track angle with the chamber
341 if(fMC) fMC->GetDirections(x0, y0, z0, dydx, dzdx);
343 y0 = fTracklet->GetYref(0);
344 z0 = fTracklet->GetZref(0);
345 dydx = fTracklet->GetYref(1);
346 dzdx = fTracklet->GetZref(1);
349 AliTRDseedV1 trklt(*fTracklet);
350 if(!trklt.Fit(kFALSE)) continue;
352 AliTRDcluster *c = 0x0;
353 fTracklet->ResetClusterIter(kFALSE);
354 while((c = fTracklet->PrevCluster())){
355 dy = trklt.GetYat(c->GetX()) - c->GetY();
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;
372 (*fDebugStream) << "ClusterResidual"
390 //________________________________________________________
391 TH1* AliTRDtrackingResolution::PlotResolution(const AliTRDtrackV1 *track)
394 AliWarning("No MC defined. Results will not be available.");
397 if(track) fTrack = track;
399 AliWarning("No Track defined.");
403 if(!(h = ((TH2I*)fContainer->At(kClusterResolution)))){
404 AliWarning("No output histogram defined.");
407 if(!(h = ((TH2I*)fContainer->At(kTrackletYResolution)))){
408 AliWarning("No output histogram defined.");
411 if(!(h = ((TH2I*)fContainer->At(kTrackletZResolution)))){
412 AliWarning("No output histogram defined.");
415 if(!(h = ((TH2I*)fContainer->At(kTrackletAngleResolution)))){
416 AliWarning("No output histogram defined.");
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);
432 // recalculate tracklet based on the MC info
433 AliTRDseedV1 tt(*fTracklet);
436 if(!tt.Fit()) continue;
437 dy = tt.GetYfit(0) - y0;
438 Float_t dphi = TMath::ATan(tt.GetYfit(1)) - TMath::ATan(dydx);
440 Bool_t cross = fTracklet->GetNChange();
442 Double_t *xyz = tt.GetCrossXYZ();
443 dz = xyz[2] - (z0 - (x0 - xyz[0])*dzdx) ;
444 ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(dzdx, dz);
448 ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(dydx, dy);
449 ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(dydx, dphi*TMath::RadToDeg());
453 Float_t p = fMC->GetTrackRef() ? fMC->GetTrackRef()->P() : -1.;
454 (*fDebugStream) << "TrkltResolution"
469 Int_t istk = AliTRDgeometry::GetStack(det);
470 AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
471 Float_t zr0 = pp->GetRow0() + AliTRDSimParam::Instance()->GetAnodeWireOffset();
472 Float_t tilt = fTracklet->GetTilt();
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));
487 if(q>100.) ((TH2I*)fContainer->At(kClusterResolution))->Fill(dydx, dy);
491 Float_t d = zr0 - zt;
492 d -= ((Int_t)(2 * d)) / 2.0;
493 if (d > 0.25) d = 0.5 - d;
495 (*fDebugStream) << "ClusterResolution"
512 //________________________________________________________
513 void AliTRDtrackingResolution::GetRefFigure(Int_t ifig)
516 TGraphErrors *g = 0x0;
518 case kClusterResidual:
519 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
521 ax = g->GetHistogram()->GetYaxis();
522 ax->SetRangeUser(-.5, 1.);
523 ax->SetTitle("Clusters Y Residuals #sigma/#mu [mm]");
524 ax = g->GetHistogram()->GetXaxis();
525 ax->SetTitle("tg(#phi)");
526 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
529 case kClusterResolution:
530 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
531 ax = g->GetHistogram()->GetYaxis();
532 ax->SetRangeUser(-.5, 1.);
533 ax->SetTitle("Cluster Y Resolution #sigma/#mu [mm]");
534 ax = g->GetHistogram()->GetXaxis();
535 ax->SetTitle("tg(#phi)");
537 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
540 case kTrackletYResolution:
541 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
542 ax = g->GetHistogram()->GetYaxis();
543 ax->SetRangeUser(-.5, 1.);
544 ax->SetTitle("Tracklet Y Resolution #sigma/#mu [mm]");
545 ax = g->GetHistogram()->GetXaxis();
546 ax->SetTitle("#tg(phi)");
548 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
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)");
559 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
562 case kTrackletAngleResolution:
563 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
564 ax = g->GetHistogram()->GetYaxis();
565 ax->SetRangeUser(-.05, .2);
566 ax->SetTitle("Tracklet Angular Resolution #sigma/#mu [deg]");
567 ax = g->GetHistogram()->GetXaxis();
568 ax->SetTitle("tg(#phi)");
570 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
574 AliInfo(Form("Reference plot [%d] not implemented yet", ifig));
577 AliInfo(Form("Reference plot [%d] missing result", ifig));
581 //________________________________________________________
582 Bool_t AliTRDtrackingResolution::PostProcess()
584 //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
586 Printf("ERROR: list not available");
589 fNRefFigures = fContainer->GetEntriesFast();
591 fGraphS = new TObjArray(fNRefFigures);
595 fGraphM = new TObjArray(fNRefFigures);
601 TGraphErrors *gm = 0x0, *gs = 0x0;
604 TF1 f("f1", "gaus", -.5, .5);
606 TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
608 TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
611 if(IsVisual()) c = new TCanvas("c", Form("%s Visual", GetName()), 500, 500);
613 sprintf(opt, "%c%c", IsVerbose() ? ' ' : 'Q', IsVisual() ? ' ': 'N');
616 //PROCESS RESIDUAL DISTRIBUTIONS
618 // Clusters residuals
619 h2 = (TH2I *)(fContainer->At(kClusterResidual));
620 gm = new TGraphErrors(h2->GetNbinsX());
621 gm->SetLineColor(kBlue);
622 gm->SetMarkerStyle(7);
623 gm->SetMarkerColor(kBlue);
624 gm->SetNameTitle("clm", "");
625 fGraphM->AddAt(gm, kClusterResidual);
626 gs = new TGraphErrors(h2->GetNbinsX());
627 gs->SetLineColor(kRed);
628 gs->SetMarkerStyle(23);
629 gs->SetMarkerColor(kRed);
630 gs->SetNameTitle("cls", "");
631 fGraphS->AddAt(gs, kClusterResidual);
632 for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
633 Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
634 h = h2->ProjectionY("py", ibin, ibin);
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);}
641 gm->SetPoint(ibin - 1, TMath::Tan(phi*TMath::DegToRad()), 10.*fc.GetParameter(1));
642 gm->SetPointError(ibin - 1, 0., 10.*fc.GetParError(1));
643 gs->SetPoint(ibin - 1, TMath::Tan(phi*TMath::DegToRad()), 10.*fc.GetParameter(2));
644 gs->SetPointError(ibin - 1, 0., 10.*fc.GetParError(2));
648 //PROCESS RESOLUTION DISTRIBUTIONS
651 // cluster y resolution
652 h2 = (TH2I*)fContainer->At(kClusterResolution);
653 gm = new TGraphErrors(h2->GetNbinsX());
654 gm->SetLineColor(kBlue);
655 gm->SetMarkerStyle(7);
656 gm->SetMarkerColor(kBlue);
657 gm->SetNameTitle("clym", "");
658 fGraphM->AddAt(gm, kClusterResolution);
659 gs = new TGraphErrors(h2->GetNbinsX());
660 gs->SetLineColor(kRed);
661 gs->SetMarkerStyle(23);
662 gs->SetMarkerColor(kRed);
663 gs->SetNameTitle("clys", "");
664 fGraphS->AddAt(gs, kClusterResolution);
665 for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
666 h = h2->ProjectionY("py", iphi, iphi);
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);}
673 Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
674 Int_t jphi = iphi -1;
675 gm->SetPoint(jphi, TMath::Tan(phi*TMath::DegToRad()), 10.*fb.GetParameter(1));
676 gm->SetPointError(jphi, 0., 10.*fb.GetParError(1));
677 gs->SetPoint(jphi, TMath::Tan(phi*TMath::DegToRad()), 10.*fb.GetParameter(2));
678 gs->SetPointError(jphi, 0., 10.*fb.GetParError(2));
681 // tracklet y resolution
682 h2 = (TH2I*)fContainer->At(kTrackletYResolution);
683 gm = new TGraphErrors(h2->GetNbinsX());
684 gm->SetLineColor(kBlue);
685 gm->SetMarkerStyle(7);
686 gm->SetMarkerColor(kBlue);
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);
695 for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
696 h = h2->ProjectionY("py", iphi, iphi);
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);}
703 Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
704 Int_t jphi = iphi -1;
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));
711 // tracklet phi resolution
712 h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
713 gm = new TGraphErrors(h2->GetNbinsX());
714 gm->SetLineColor(kBlue);
715 gm->SetMarkerStyle(7);
716 gm->SetMarkerColor(kBlue);
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);
725 for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
726 h = h2->ProjectionY("py", iphi, iphi);
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);}
732 Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
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));
746 //________________________________________________________
747 void AliTRDtrackingResolution::Terminate(Option_t *)
754 if(HasPostProcess()) PostProcess();
757 //________________________________________________________
758 void AliTRDtrackingResolution::AdjustF1(TH1 *h, TF1 *f)
760 // Helper function to avoid duplication of code
761 // Make first guesses on the fit parameters
763 // find the intial parameters of the fit !! (thanks George)
764 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
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);
770 f->SetParLimits(1, -.2, .2);
771 f->SetParameter(1, 0.);
773 f->SetParLimits(2, 0., 4.e-1);
774 f->SetParameter(2, 2.e-2);
775 if(f->GetNpar() <= 4) return;
777 f->SetParLimits(3, 0., sum);
778 f->SetParameter(3, .1*sum);
780 f->SetParLimits(4, -.3, .3);
781 f->SetParameter(4, 0.);
783 f->SetParLimits(5, 0., 1.e2);
784 f->SetParameter(5, 2.e-1);
787 //________________________________________________________
788 TObjArray* AliTRDtrackingResolution::Histos()
790 if(fContainer) return fContainer;
792 fContainer = new TObjArray(5);
794 // cluster to tracklet residuals [2]
795 fContainer->AddAt(new TH2I("fYClRes", "Clusters Residuals", 21, -.33, .33, 100, -.5, .5), kClusterResidual);
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);
804 // cluster y resolution [0]
805 fContainer->AddAt(new TH2I("fCY", "Cluster Resolution", 31, -.48, .48, 100, -.5, .5), kClusterResolution);
806 // tracklet y resolution [0]
807 fContainer->AddAt(new TH2I("fY", "Tracklet Resolution", 31, -.48, .48, 100, -.5, .5), kTrackletYResolution);
808 // tracklet y resolution [0]
809 fContainer->AddAt(new TH2I("fY", "Tracklet Resolution", 31, -.48, .48, 100, -.5, .5), kTrackletZResolution);
810 // tracklet angular resolution [1]
811 fContainer->AddAt(new TH2I("fPhi", "Tracklet Angular Resolution", 31, -.48, .48, 100, -10., 10.), kTrackletAngleResolution);
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);
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);
827 //________________________________________________________
828 void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
831 fReconstructor->SetRecoParam(r);