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: AliTRDresolution.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 // AliTRDresolution *res = new AliTRDresolution();
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 ////////////////////////////////////////////////////////////////////////////
56 #include <TObjArray.h>
63 #include <TGraphErrors.h>
65 #include "TTreeStream.h"
66 #include "TGeoManager.h"
68 #include "AliAnalysisManager.h"
69 #include "AliTrackReference.h"
70 #include "AliTrackPointArray.h"
71 #include "AliCDBManager.h"
73 #include "AliTRDSimParam.h"
74 #include "AliTRDgeometry.h"
75 #include "AliTRDpadPlane.h"
76 #include "AliTRDcluster.h"
77 #include "AliTRDseedV1.h"
78 #include "AliTRDtrackV1.h"
79 #include "AliTRDtrackerV1.h"
80 #include "AliTRDReconstructor.h"
81 #include "AliTRDrecoParam.h"
83 #include "AliTRDtrackInfo/AliTRDclusterInfo.h"
84 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
85 #include "AliTRDresolution.h"
87 ClassImp(AliTRDresolution)
89 //________________________________________________________
90 AliTRDresolution::AliTRDresolution()
91 :AliTRDrecoTask("Resolution", "Spatial and Momentum Resolution")
102 fReconstructor = new AliTRDReconstructor();
103 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
104 fGeo = new AliTRDgeometry();
108 DefineOutput(1, TObjArray::Class()); // cluster2track
109 DefineOutput(2, TObjArray::Class()); // tracklet2track
110 DefineOutput(3, TObjArray::Class()); // cluster2mc
111 DefineOutput(4, TObjArray::Class()); // tracklet2mc
114 //________________________________________________________
115 AliTRDresolution::~AliTRDresolution()
117 if(fGraphS){fGraphS->Delete(); delete fGraphS;}
118 if(fGraphM){fGraphM->Delete(); delete fGraphM;}
120 delete fReconstructor;
121 if(gGeoManager) delete gGeoManager;
122 if(fCl){fCl->Delete(); delete fCl;}
123 if(fTrklt){fTrklt->Delete(); delete fTrklt;}
124 if(fMCcl){fMCcl->Delete(); delete fMCcl;}
125 if(fMCtrklt){fMCtrklt->Delete(); delete fMCtrklt;}
129 //________________________________________________________
130 void AliTRDresolution::CreateOutputObjects()
132 // spatial resolution
133 OpenFile(0, "RECREATE");
135 fContainer = Histos();
137 fCl = new TObjArray();
138 fCl->SetOwner(kTRUE);
139 fTrklt = new TObjArray();
140 fTrklt->SetOwner(kTRUE);
141 fMCcl = new TObjArray();
142 fMCcl->SetOwner(kTRUE);
143 fMCtrklt = new TObjArray();
144 fMCtrklt->SetOwner(kTRUE);
147 //________________________________________________________
148 void AliTRDresolution::Exec(Option_t *opt)
155 AliTRDrecoTask::Exec(opt);
160 PostData(4, fMCtrklt);
163 //________________________________________________________
164 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
166 if(track) fTrack = track;
168 AliWarning("No Track defined.");
172 if(!(h = ((TH2I*)fContainer->At(kCluster)))){
173 AliWarning("No output histogram defined.");
178 Float_t x0, y0, z0, dy, dydx, dzdx;
179 AliTRDseedV1 *fTracklet = 0x0;
180 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
181 if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
182 if(!fTracklet->IsOK()) continue;
183 x0 = fTracklet->GetX0();
185 // retrive the track angle with the chamber
186 y0 = fTracklet->GetYref(0);
187 z0 = fTracklet->GetZref(0);
188 dydx = fTracklet->GetYref(1);
189 dzdx = fTracklet->GetZref(1);
190 fTracklet->GetCovRef(cov);
191 Float_t tilt = fTracklet->GetTilt();
192 AliTRDcluster *c = 0x0;
193 fTracklet->ResetClusterIter(kFALSE);
194 while((c = fTracklet->PrevCluster())){
195 Float_t xc = c->GetX();
196 Float_t yc = c->GetY();
197 Float_t zc = c->GetZ();
198 Float_t dx = x0 - xc;
199 Float_t yt = y0 - dx*dydx;
200 Float_t zt = z0 - dx*dzdx;
201 yc -= tilt*(zc-zt); // tilt correction
204 h->Fill(dydx, dy/TMath::Sqrt(cov[0] + c->GetSigmaY2()));
207 // Get z-position with respect to anode wire
208 //AliTRDSimParam *simParam = AliTRDSimParam::Instance();
209 Int_t istk = fGeo->GetStack(c->GetDetector());
210 AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
211 Float_t row0 = pp->GetRow0();
212 Float_t d = row0 - zt + pp->GetAnodeWireOffset();
213 d -= ((Int_t)(2 * d)) / 2.0;
214 if (d > 0.25) d = 0.5 - d;
216 /* AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
218 clInfo->SetCluster(c);
219 clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
220 clInfo->SetResolution(dy);
221 clInfo->SetAnisochronity(d);
222 clInfo->SetDriftLength(dx);
223 (*fDebugStream) << "ClusterResiduals"
224 <<"clInfo.=" << clInfo
233 //________________________________________________________
234 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
236 if(track) fTrack = track;
238 AliWarning("No Track defined.");
242 if(!(h = ((TH2I*)fContainer->At(kTrackletY)))){
243 AliWarning("No output histogram defined.");
247 Double_t cov[3], covR[3];
248 Float_t x, y0, dx, dy, dydx;
249 AliTRDseedV1 *fTracklet = 0x0;
250 for(Int_t il=AliTRDgeometry::kNlayer; il--;){
251 if(!(fTracklet = fTrack->GetTracklet(il))) continue;
252 if(!fTracklet->IsOK()) continue;
253 y0 = fTracklet->GetYref(0);
254 dydx = fTracklet->GetYref(1);
255 x = fTracklet->GetX();
256 dx = fTracklet->GetX0() - x;
257 dy = y0-dx*dydx - fTracklet->GetY();
258 fTracklet->GetCovAt(x, cov);
259 fTracklet->GetCovRef(covR);
260 h->Fill(dydx, dy/*/TMath::Sqrt(cov[0] + covR[0])*/);
265 //________________________________________________________
266 TH1* AliTRDresolution::PlotTrackletPhi(const AliTRDtrackV1 *track)
268 if(track) fTrack = track;
270 AliWarning("No Track defined.");
274 if(!(h = ((TH2I*)fContainer->At(kTrackletPhi)))){
275 AliWarning("No output histogram defined.");
279 AliTRDseedV1 *tracklet = 0x0;
280 for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
281 if(!(tracklet = fTrack->GetTracklet(il))) continue;
282 if(!tracklet->IsOK()) continue;
283 h->Fill(tracklet->GetYref(1), TMath::ATan(tracklet->GetYref(1))-TMath::ATan(tracklet->GetYfit(1)));
289 //________________________________________________________
290 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
293 AliWarning("No MC defined. Results will not be available.");
296 if(track) fTrack = track;
298 AliWarning("No Track defined.");
303 Int_t pdg = fMC->GetPDG(), det=-1;
304 Int_t label = fMC->GetLabel();
306 Float_t p, pt, x0, y0, z0, dx, dy, dz, dydx, dzdx;
307 Double_t covR[3], cov[3];
310 Double_t DX[12], DY[12], DZ[12], DPt[12], COV[12][15];
311 fMC->PropagateKalman(DX, DY, DZ, DPt, COV);
312 (*fDebugStream) << "MCkalman"
329 AliTRDseedV1 *fTracklet = 0x0;
330 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
331 if(!(fTracklet = fTrack->GetTracklet(ily)))/* ||
332 !fTracklet->IsOK())*/ continue;
334 det = fTracklet->GetDetector();
335 x0 = fTracklet->GetX0();
336 //radial shift with respect to the MC reference (radial position of the pad plane)
337 x= fTracklet->GetX();
338 if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, pt, s)) continue;
339 // MC track position at reference radial position
342 (*fDebugStream) << "MC"
354 Float_t yt = y0 - dx*dydx;
355 Float_t zt = z0 - dx*dzdx;
356 p = pt*(1.+dzdx*dzdx); // pt -> p
358 // add Kalman residuals for y, z and pt
359 dx = fTracklet->GetX0() - x;
360 Float_t yr = fTracklet->GetYref(0) - dx*fTracklet->GetYref(1);
362 Float_t zr = fTracklet->GetZref(0) - dx*fTracklet->GetZref(1);
364 Float_t tgl = fTracklet->GetTgl();
365 Float_t dpt = pt - fTracklet->GetMomentum()/(1.+tgl*tgl);
366 fTracklet->GetCovRef(covR);
368 ((TH2I*)fContainer->At(kMCtrackY))->Fill(dydx, dy);
369 ((TH2I*)fContainer->At(kMCtrackYPull))->Fill(dydx, dy/TMath::Sqrt(covR[0]));
371 ((TH2I*)fContainer->At(kMCtrackZIn))->Fill(dzdx, dz);
372 ((TH2I*)fContainer->At(kMCtrackZInPull))->Fill(dzdx, dz/TMath::Sqrt(covR[2]));
373 } else if(ily==AliTRDgeometry::kNlayer-1) {
374 ((TH2I*)fContainer->At(kMCtrackZOut))->Fill(dzdx, dz);
375 ((TH2I*)fContainer->At(kMCtrackZOutPull))->Fill(dzdx, dz/TMath::Sqrt(covR[2]));
377 if(pdg!=kElectron && pdg!=kPositron) ((TH2I*)fContainer->At(kMCtrackPt))->Fill(1./pt, dpt);
378 // Fill Debug stream for Kalman track
380 Float_t dydxr = fTracklet->GetYref(1);
381 (*fDebugStream) << "MCtrack"
393 // recalculate tracklet based on the MC info
394 AliTRDseedV1 tt(*fTracklet);
398 x= tt.GetX(); // the true one
401 Bool_t rc = tt.IsRowCross();
403 // add tracklet residuals for y and dydx
404 Float_t yf = tt.GetY();
405 dy = yt - yf; dz = 100.;
406 Float_t dphi = (tt.GetYfit(1) - dydx);
407 dphi /= 1.- tt.GetYfit(1)*dydx;
410 ((TH2I*)fContainer->At(kMCtrackletY))->Fill(dydx, dy);
411 if(cov[0]>0.) ((TH2I*)fContainer->At(kMCtrackletYPull))->Fill(dydx, dy/TMath::Sqrt(cov[0]));
412 ((TH2I*)fContainer->At(kMCtrackletPhi))->Fill(dydx, dphi*TMath::RadToDeg());
414 // add tracklet residuals for z
415 dz = tt.GetZ() - (z0 - dx*dzdx) ;
416 ((TH2I*)fContainer->At(kMCtrackletZ))->Fill(dzdx, dz);
417 if(cov[2]>0.) ((TH2I*)fContainer->At(kMCtrackletZPull))->Fill(dzdx, dz/TMath::Sqrt(cov[2]));
420 // Fill Debug stream for tracklet
422 (*fDebugStream) << "MCtracklet"
433 Int_t istk = AliTRDgeometry::GetStack(det);
434 AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
435 Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
436 Float_t tilt = fTracklet->GetTilt();
438 AliTRDcluster *c = 0x0;
439 fTracklet->ResetClusterIter(kFALSE);
440 while((c = fTracklet->PrevCluster())){
441 Float_t q = TMath::Abs(c->GetQ());
442 AliTRDseedV1::GetClusterXY(c,x,y);
443 //x = c->GetX(); y = c->GetY();
446 Float_t zc = c->GetZ();
450 dy = yt - (yc - tilt*(zc-zt));
453 if(q>20. && q<250.) ((TH2I*)fContainer->At(kMCcluster))->Fill(dydx, dy);
455 // Fill calibration container
456 Float_t d = zr0 - zt;
457 d -= ((Int_t)(2 * d)) / 2.0;
458 if (d > 0.25) d = 0.5 - d;
459 AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
461 clInfo->SetCluster(c);
462 clInfo->SetMC(pdg, label);
463 clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
464 clInfo->SetResolution(dy);
465 clInfo->SetAnisochronity(d);
466 clInfo->SetDriftLength(dx-.5*AliTRDgeometry::CamHght());
467 clInfo->SetTilt(tilt);
472 (*fDebugStream) << "MCcluster"
473 <<"clInfo.=" << clInfo
482 //________________________________________________________
483 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
485 Float_t y[2] = {0., 0.};
488 TGraphErrors *g = 0x0;
491 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
493 ax = g->GetHistogram()->GetYaxis();
494 y[0] = -0.5; y[1] = 2.5;
495 ax->SetRangeUser(y[0], y[1]);
496 ax->SetTitle("Cluster-Track Pulls #sigma/#mu [mm]");
497 ax = g->GetHistogram()->GetXaxis();
498 ax->SetTitle("tg(#phi)");
499 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
501 b = new TBox(-.15, y[0], .15, y[1]);
502 b->SetFillStyle(3002);b->SetFillColor(kGreen);
503 b->SetLineColor(0); b->Draw();
506 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
508 ax = g->GetHistogram()->GetYaxis();
509 ax->SetRangeUser(-.5, 3.);
510 ax->SetTitle("Tracklet-Track Y-Pulls #sigma/#mu [mm]");
511 ax = g->GetHistogram()->GetXaxis();
512 ax->SetTitle("tg(#phi)");
513 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
515 b = new TBox(-.15, y[0], .15, y[1]);
516 b->SetFillStyle(3002);b->SetFillColor(kGreen);
517 b->SetLineColor(0); b->Draw();
520 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
522 ax = g->GetHistogram()->GetYaxis();
523 ax->SetRangeUser(-.5, 2.);
524 ax->SetTitle("Tracklet-Track Phi-Residuals #sigma/#mu [rad]");
525 ax = g->GetHistogram()->GetXaxis();
526 ax->SetTitle("tg(#phi)");
527 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
529 b = new TBox(-.15, y[0], .15, y[1]);
530 b->SetFillStyle(3002);b->SetFillColor(kGreen);
531 b->SetLineColor(0); b->Draw();
534 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
535 ax = g->GetHistogram()->GetYaxis();
536 y[0] = -.05; y[1] = 0.6;
537 ax->SetRangeUser(y[0], y[1]);
538 ax->SetTitle("Y_{cluster} #sigma/#mu [mm]");
539 ax = g->GetHistogram()->GetXaxis();
540 ax->SetRangeUser(-.3, .3);
541 ax->SetTitle("tg(#phi)");
543 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
545 b = new TBox(-.15, y[0], .15, y[1]);
546 b->SetFillStyle(3002);b->SetFillColor(kBlue);
547 b->SetLineColor(0); b->Draw();
550 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
551 ax = g->GetHistogram()->GetYaxis();
552 y[0] = -.05; y[1] = 0.25;
553 ax->SetRangeUser(y[0], y[1]);
554 ax->SetTitle("Y_{tracklet} #sigma/#mu [mm]");
555 ax = g->GetHistogram()->GetXaxis();
556 ax->SetRangeUser(-.2, .2);
557 ax->SetTitle("tg(#phi)");
559 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
561 b = new TBox(-.15, y[0], .15, y[1]);
562 b->SetFillStyle(3002);b->SetFillColor(kBlue);
563 b->SetLineColor(0); b->Draw();
566 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
567 ax = g->GetHistogram()->GetYaxis();
568 ax->SetRangeUser(-.5, 1.);
569 ax->SetTitle("Z_{tracklet} #sigma/#mu [mm]");
570 ax = g->GetHistogram()->GetXaxis();
571 ax->SetTitle("tg(#theta)");
573 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
577 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
578 ax = g->GetHistogram()->GetYaxis();
579 y[0] = -.05; y[1] = .2;
580 ax->SetRangeUser(y[0], y[1]);
581 ax->SetTitle("#Phi_{tracklet} #sigma/#mu [deg]");
582 ax = g->GetHistogram()->GetXaxis();
583 ax->SetTitle("tg(#phi)");
585 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
589 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
590 ax = g->GetHistogram()->GetYaxis();
591 y[0] = -.05; y[1] = 0.25;
592 ax->SetRangeUser(y[0], y[1]);
593 ax->SetTitle("Y_{track} #sigma/#mu [mm]");
594 ax = g->GetHistogram()->GetXaxis();
595 ax->SetRangeUser(-.2, .2);
596 ax->SetTitle("tg(#phi)");
598 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
600 b = new TBox(-.15, y[0], .15, y[1]);
601 b->SetFillStyle(3002);b->SetFillColor(kBlue);
602 b->SetLineColor(0); b->Draw();
606 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
607 ax = g->GetHistogram()->GetYaxis();
608 ax->SetRangeUser(-.5, 2.);
609 ax->SetTitle(Form("Z_{track}^{%s} #sigma/#mu [mm]", ifig==kMCtrackZIn ? "in" : "out"));
610 ax = g->GetHistogram()->GetXaxis();
611 ax->SetTitle("tg(#theta)");
613 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
617 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
618 ax = g->GetHistogram()->GetYaxis();
619 ax->SetRangeUser(-.5, 2.);
620 ax->SetTitle("#epsilon_{P_{t}}^{track} / #mu [%]");
621 ax = g->GetHistogram()->GetXaxis();
622 ax->SetTitle("1/p_{t}");
624 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
627 case kMCtrackletYPull:
628 case kMCtrackletZPull:
630 case kMCtrackZInPull:
631 case kMCtrackZOutPull:
632 if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
633 ax = g->GetHistogram()->GetYaxis();
634 ax->SetRangeUser(-.5, 2.);
635 ax->SetTitle("MC Pulls");
637 if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
641 AliInfo(Form("Reference plot [%d] missing result", ifig));
646 //________________________________________________________
647 Bool_t AliTRDresolution::PostProcess()
649 //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
651 Printf("ERROR: list not available");
654 fNRefFigures = fContainer->GetEntriesFast();
655 TGraphErrors *gm= 0x0, *gs= 0x0;
657 fGraphS = new TObjArray(fNRefFigures);
659 for(Int_t ig=0; ig<fNRefFigures; ig++){
660 gs = new TGraphErrors();
661 gs->SetLineColor(kRed);
662 gs->SetMarkerStyle(23);
663 gs->SetMarkerColor(kRed);
664 gs->SetNameTitle(Form("s_%d", ig), "");
665 fGraphS->AddAt(gs, ig);
669 fGraphM = new TObjArray(fNRefFigures);
671 for(Int_t ig=0; ig<fNRefFigures; ig++){
672 gm = new TGraphErrors();
673 gm->SetLineColor(kBlue);
674 gm->SetMarkerStyle(7);
675 gm->SetMarkerColor(kBlue);
676 gm->SetNameTitle(Form("m_%d", ig), "");
677 fGraphM->AddAt(gm, ig);
683 TF1 f("f1", "gaus", -.5, .5);
684 // gauss on a constant background
685 TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
686 // gauss on a gauss background
687 TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
690 //PROCESS EXPERIMENTAL DISTRIBUTIONS
692 // Clusters residuals
693 Process(kCluster, &f);
695 // Tracklet y residuals
696 Process(kTrackletY, &f);
698 // Tracklet phi residuals
699 Process(kTrackletPhi, &f);
701 if(!HasMCdata()) return kTRUE;
704 //PROCESS MC RESIDUAL DISTRIBUTIONS
706 // cluster y resolution
707 Process(kMCcluster, &f);
709 // tracklet resolution
710 Process(kMCtrackletY, &f); // y
711 Process(kMCtrackletZ, &f); // z
712 Process(kMCtrackletPhi, &f); // phi
715 Process(kMCtrackletYPull, &f); // y
716 Process(kMCtrackletZPull, &f); // z
719 Process(kMCtrackY, &f); // y
720 Process(kMCtrackZIn, &f); // z towards TPC
721 Process(kMCtrackZOut, &f); // z towards TOF
724 Process(kMCtrackYPull, &f); // y
725 Process(kMCtrackZInPull, &f); // z towards TPC
726 Process(kMCtrackZOutPull, &f); // z towards TOF
728 // track Pt resolution
729 TH2I *h2 = (TH2I*)fContainer->At(kMCtrackPt);
730 TAxis *ax = h2->GetXaxis();
731 gm = (TGraphErrors*)fGraphM->At(kMCtrackPt);
732 gs = (TGraphErrors*)fGraphS->At(kMCtrackPt);
733 TF1 fg("fg", "gaus", -1.5, 1.5);
734 TF1 fl("fl", "landau", -4., 15.);
735 TF1 fgl("fgl", "gaus(0)+landau(3)", -5., 20.);
736 for(Int_t ip=1; ip<=ax->GetNbins(); ip++){
737 TH1D *h = h2->ProjectionY("ppt", ip, ip);
738 if(h->GetEntries()<70) continue;
740 h->Fit(&fg, "QN", "", -1.5, 1.5);
741 fgl.SetParameter(0, fg.GetParameter(0));
742 fgl.SetParameter(1, fg.GetParameter(1));
743 fgl.SetParameter(2, fg.GetParameter(2));
744 h->Fit(&fl, "QN", "", -4., 15.);
745 fgl.SetParameter(3, fl.GetParameter(0));
746 fgl.SetParameter(4, fl.GetParameter(1));
747 fgl.SetParameter(5, fl.GetParameter(2));
749 h->Fit(&fgl, "NQ", "", -5., 20.);
751 Float_t invpt = ax->GetBinCenter(ip);
752 Int_t ip = gm->GetN();
753 gm->SetPoint(ip, invpt, fgl.GetParameter(1));
754 gm->SetPointError(ip, 0., fgl.GetParError(1));
755 gs->SetPoint(ip, invpt, fgl.GetParameter(2)*invpt);
756 gs->SetPointError(ip, 0., fgl.GetParError(2));
757 // fgl.GetParameter(4) // Landau MPV
758 // fgl.GetParameter(5) // Landau Sigma
766 //________________________________________________________
767 void AliTRDresolution::Terminate(Option_t *)
774 if(HasPostProcess()) PostProcess();
777 //________________________________________________________
778 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
780 // Helper function to avoid duplication of code
781 // Make first guesses on the fit parameters
783 // find the intial parameters of the fit !! (thanks George)
784 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
786 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
787 f->SetParLimits(0, 0., 3.*sum);
788 f->SetParameter(0, .9*sum);
790 f->SetParLimits(1, -.2, .2);
791 f->SetParameter(1, -0.1);
793 f->SetParLimits(2, 0., 4.e-1);
794 f->SetParameter(2, 2.e-2);
795 if(f->GetNpar() <= 4) return;
797 f->SetParLimits(3, 0., sum);
798 f->SetParameter(3, .1*sum);
800 f->SetParLimits(4, -.3, .3);
801 f->SetParameter(4, 0.);
803 f->SetParLimits(5, 0., 1.e2);
804 f->SetParameter(5, 2.e-1);
807 //________________________________________________________
808 TObjArray* AliTRDresolution::Histos()
810 if(fContainer) return fContainer;
812 fContainer = new TObjArray(kNhistos);
813 //fContainer->SetOwner(kTRUE);
816 // cluster to tracklet residuals [2]
817 if(!(h = (TH2I*)gROOT->FindObject("hCl"))){
818 h = new TH2I("hCl", "Clusters-Track Residuals", 21, -.33, .33, 100, -.5, .5);
819 h->GetXaxis()->SetTitle("tg(#phi)");
820 h->GetYaxis()->SetTitle("#Delta y [cm]");
821 h->GetZaxis()->SetTitle("entries");
823 fContainer->AddAt(h, kCluster);
825 // tracklet to track residuals [2]
826 if(!(h = (TH2I*)gROOT->FindObject("hTrkltY"))){
827 h = new TH2I("hTrkltY", "Tracklets-Track Residuals (Y)", 21, -.33, .33, 100, -.5, .5);
828 h->GetXaxis()->SetTitle("tg(#phi)");
829 h->GetYaxis()->SetTitle("#Delta y [cm]");
830 h->GetZaxis()->SetTitle("entries");
832 fContainer->AddAt(h, kTrackletY);
834 // tracklet to track residuals angular [2]
835 if(!(h = (TH2I*)gROOT->FindObject("hTrkltPhi"))){
836 h = new TH2I("hTrkltPhi", "Tracklets-Track Residuals (#Phi)", 21, -.33, .33, 100, -.5, .5);
837 h->GetXaxis()->SetTitle("tg(#phi)");
838 h->GetYaxis()->SetTitle("#Delta phi [#circ]");
839 h->GetZaxis()->SetTitle("entries");
841 fContainer->AddAt(h, kTrackletPhi);
845 if(!HasMCdata()) return fContainer;
847 // cluster y resolution [0]
848 if(!(h = (TH2I*)gROOT->FindObject("hMCcl"))){
849 h = new TH2I("hMCcl", "Cluster Resolution", 31, -.48, .48, 100, -.5, .5);
850 h->GetXaxis()->SetTitle("tg(#phi)");
851 h->GetYaxis()->SetTitle("#Delta y [cm]");
852 h->GetZaxis()->SetTitle("entries");
854 fContainer->AddAt(h, kMCcluster);
856 // tracklet y resolution [0]
857 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltY"))){
858 h = new TH2I("hMCtrkltY", "Tracklet Resolution (Y)", 31, -.48, .48, 100, -.2, .2);
859 h->GetXaxis()->SetTitle("tg(#phi)");
860 h->GetYaxis()->SetTitle("#Delta y [cm]");
861 h->GetZaxis()->SetTitle("entries");
863 fContainer->AddAt(h, kMCtrackletY);
865 // tracklet y resolution [0]
866 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltYPull"))){
867 h = new TH2I("hMCtrkltYPull", "Tracklet Pulls (Y)", 31, -.48, .48, 100, -3.5, 3.5);
868 h->GetXaxis()->SetTitle("tg(#phi)");
869 h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
870 h->GetZaxis()->SetTitle("entries");
872 fContainer->AddAt(h, kMCtrackletYPull);
874 // tracklet y resolution [0]
875 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltZ"))){
876 h = new TH2I("hMCtrkltZ", "Tracklet Resolution (Z)", 31, -.48, .48, 100, -.5, .5);
877 h->GetXaxis()->SetTitle("tg(#theta)");
878 h->GetYaxis()->SetTitle("#Delta z [cm]");
879 h->GetZaxis()->SetTitle("entries");
881 fContainer->AddAt(h, kMCtrackletZ);
883 // tracklet y resolution [0]
884 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltZPull"))){
885 h = new TH2I("hMCtrkltZPull", "Tracklet Pulls (Z)", 31, -.48, .48, 100, -3.5, 3.5);
886 h->GetXaxis()->SetTitle("tg(#theta)");
887 h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
888 h->GetZaxis()->SetTitle("entries");
890 fContainer->AddAt(h, kMCtrackletZPull);
892 // tracklet angular resolution [1]
893 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltPhi"))){
894 h = new TH2I("hMCtrkltPhi", "Tracklet Resolution (#Phi)", 31, -.48, .48, 100, -10., 10.);
895 h->GetXaxis()->SetTitle("tg(#phi)");
896 h->GetYaxis()->SetTitle("#Delta #phi [deg]");
897 h->GetZaxis()->SetTitle("entries");
899 fContainer->AddAt(h, kMCtrackletPhi);
901 // Kalman track y resolution
902 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkY"))){
903 h = new TH2I("hMCtrkY", "Kalman Track Resolution (Y)", 31, -.48, .48, 100, -.2, .2);
904 h->GetXaxis()->SetTitle("tg(#phi)");
905 h->GetYaxis()->SetTitle("#Delta y [cm]");
906 h->GetZaxis()->SetTitle("entries");
908 fContainer->AddAt(h, kMCtrackY);
910 // Kalman track y resolution
911 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkYPull"))){
912 h = new TH2I("hMCtrkYPull", "Kalman Track Pulls (Y)", 31, -.48, .48, 100, -3.5, 3.5);
913 h->GetXaxis()->SetTitle("tg(#phi)");
914 h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
915 h->GetZaxis()->SetTitle("entries");
917 fContainer->AddAt(h, kMCtrackYPull);
919 // Kalman track Z resolution
920 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZIn"))){
921 h = new TH2I("hMCtrkZIn", "Kalman Track Resolution (Zin)", 20, -1., 1., 100, -1., 1.);
922 h->GetXaxis()->SetTitle("tg(#theta)");
923 h->GetYaxis()->SetTitle("#Delta z [cm]");
924 h->GetZaxis()->SetTitle("entries");
926 fContainer->AddAt(h, kMCtrackZIn);
928 // Kalman track Z resolution
929 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZOut"))){
930 h = new TH2I("hMCtrkZOut", "Kalman Track Resolution (Zout)", 20, -1., 1., 100, -1., 1.);
931 h->GetXaxis()->SetTitle("tg(#theta)");
932 h->GetYaxis()->SetTitle("#Delta z [cm]");
933 h->GetZaxis()->SetTitle("entries");
935 fContainer->AddAt(h, kMCtrackZOut);
937 // Kalman track Z resolution
938 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZInPull"))){
939 h = new TH2I("hMCtrkZInPull", "Kalman Track Pulls (Zin)", 20, -1., 1., 100, -4.5, 4.5);
940 h->GetXaxis()->SetTitle("tg(#theta)");
941 h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
942 h->GetZaxis()->SetTitle("entries");
944 fContainer->AddAt(h, kMCtrackZInPull);
946 // Kalman track Z resolution
947 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZOutPull"))){
948 h = new TH2I("hMCtrkZOutPull", "Kalman Track Pulls (Zout)", 20, -1., 1., 100, -4.5, 4.5);
949 h->GetXaxis()->SetTitle("tg(#theta)");
950 h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
951 h->GetZaxis()->SetTitle("entries");
953 fContainer->AddAt(h, kMCtrackZOutPull);
955 // Kalman track Pt resolution
956 if(!(h = (TH2I*)gROOT->FindObject("hMCtrkPt"))){
957 h = new TH2I("hMCtrkPt", "Kalman Track Resolution (Pt)", 100, 0., 2., 150, -5., 20.);
958 h->GetXaxis()->SetTitle("1/p_{t}");
959 h->GetYaxis()->SetTitle("#Delta p_{t} [GeV/c]");
960 h->GetZaxis()->SetTitle("entries");
962 fContainer->AddAt(h, kMCtrackPt);
968 //________________________________________________________
969 Bool_t AliTRDresolution::Process(ETRDresolutionPlot ip, TF1 *f)
971 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
972 Bool_t kBUILD = kFALSE;
974 f = new TF1("f1", "gaus", -.5, .5);
979 if(!(h2 = (TH2I *)(fContainer->At(ip)))) return kFALSE;
980 TGraphErrors *gm = 0x0, *gs = 0x0;
981 if(!(gm=(TGraphErrors*)fGraphM->At(ip))) return kFALSE;
982 if(gm->GetN()) for(Int_t ip=gm->GetN(); ip--;) gm->RemovePoint(ip);
983 if(!(gs=(TGraphErrors*)fGraphS->At(ip))) return kFALSE;
984 if(gs->GetN()) for(Int_t ip=gs->GetN(); ip--;) gs->RemovePoint(ip);
985 for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
986 Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
987 TH1D *h = h2->ProjectionY("py", ibin, ibin);
988 if(h->GetEntries()<100) continue;
991 // if(IsVisual()){c->cd(); c->SetLogy();}
992 h->Fit(f, "Q", "", -0.5, 0.5);
993 /* if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}*/
995 Int_t ip = gm->GetN();
996 gm->SetPoint(ip, x, 10.*f->GetParameter(1));
997 gm->SetPointError(ip, 0., 10.*f->GetParError(1));
998 gs->SetPoint(ip, x, 10.*f->GetParameter(2));
999 gs->SetPointError(ip, 0., 10.*f->GetParError(2));
1002 if(kBUILD) delete f;
1006 //________________________________________________________
1007 void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
1010 fReconstructor->SetRecoParam(r);