Major update for the TRD tracking code
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDresolution.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercialf purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15
16 /* $Id: AliTRDresolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 //  TRD tracking resolution                                               //
21 //
22 // The class performs resolution and residual studies 
23 // of the TRD tracks for the following quantities :
24 //   - spatial position (y, [z])
25 //   - angular (phi) tracklet
26 //   - momentum at the track level
27 // 
28 // The class has to be used for regular detector performance checks using the official macros:
29 //   - $ALICE_ROOT/TRD/qaRec/run.C
30 //   - $ALICE_ROOT/TRD/qaRec/makeResults.C
31 // 
32 // For stand alone usage please refer to the following example: 
33 // {  
34 //   gSystem->Load("libANALYSIS.so");
35 //   gSystem->Load("libTRDqaRec.so");
36 //   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);
43 // }  
44 //
45 //  Authors:                                                              //
46 //    Alexandru Bercuci <A.Bercuci@gsi.de>                                //
47 //    Markus Fasel <M.Fasel@gsi.de>                                       //
48 //                                                                        //
49 ////////////////////////////////////////////////////////////////////////////
50
51 #include <cstring>
52
53 #include <TROOT.h>
54 #include <TSystem.h>
55 #include <TPDGCode.h>
56 #include <TObjArray.h>
57 #include <TH2.h>
58 #include <TH1.h>
59 #include <TF1.h>
60 #include <TCanvas.h>
61 #include <TBox.h>
62 #include <TProfile.h>
63 #include <TGraphErrors.h>
64 #include <TMath.h>
65 #include "TTreeStream.h"
66 #include "TGeoManager.h"
67
68 #include "AliAnalysisManager.h"
69 #include "AliTrackReference.h"
70 #include "AliTrackPointArray.h"
71 #include "AliCDBManager.h"
72
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"
82
83 #include "AliTRDtrackInfo/AliTRDclusterInfo.h"
84 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
85 #include "AliTRDresolution.h"
86
87 ClassImp(AliTRDresolution)
88
89 //________________________________________________________
90 AliTRDresolution::AliTRDresolution()
91   :AliTRDrecoTask("Resolution", "Spatial and Momentum Resolution")
92   ,fStatus(0)
93   ,fReconstructor(0x0)
94   ,fGeo(0x0)
95   ,fGraphS(0x0)
96   ,fGraphM(0x0)
97   ,fCl(0x0)
98   ,fTrklt(0x0)
99   ,fMCcl(0x0)
100   ,fMCtrklt(0x0)
101 {
102   fReconstructor = new AliTRDReconstructor();
103   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
104   fGeo = new AliTRDgeometry();
105
106   InitFunctorList();
107
108   DefineOutput(1, TObjArray::Class()); // cluster2track
109   DefineOutput(2, TObjArray::Class()); // tracklet2track
110   DefineOutput(3, TObjArray::Class()); // cluster2mc
111   DefineOutput(4, TObjArray::Class()); // tracklet2mc
112 }
113
114 //________________________________________________________
115 AliTRDresolution::~AliTRDresolution()
116 {
117   if(fGraphS){fGraphS->Delete(); delete fGraphS;}
118   if(fGraphM){fGraphM->Delete(); delete fGraphM;}
119   delete fGeo;
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;}
126 }
127
128
129 //________________________________________________________
130 void AliTRDresolution::CreateOutputObjects()
131 {
132   // spatial resolution
133   OpenFile(0, "RECREATE");
134
135   fContainer = Histos();
136
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);
145 }
146
147 //________________________________________________________
148 void AliTRDresolution::Exec(Option_t *opt)
149 {
150   fCl->Delete();
151   fTrklt->Delete();
152   fMCcl->Delete();
153   fMCtrklt->Delete();
154
155   AliTRDrecoTask::Exec(opt);
156
157   PostData(1, fCl);
158   PostData(2, fTrklt);
159   PostData(3, fMCcl);
160   PostData(4, fMCtrklt);
161 }
162
163 //________________________________________________________
164 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
165 {
166   if(track) fTrack = track;
167   if(!fTrack){
168     AliWarning("No Track defined.");
169     return 0x0;
170   }
171   TH1 *h = 0x0;
172   if(!(h = ((TH2I*)fContainer->At(kCluster)))){
173     AliWarning("No output histogram defined.");
174     return 0x0;
175   }
176
177   Double_t cov[3];
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();
184
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
202       dy = yt - yc;
203
204       h->Fill(dydx, dy/TMath::Sqrt(cov[0] + c->GetSigmaY2()));
205   
206       if(fDebugLevel>=1){
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;
215
216 /*        AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
217         fCl->Add(clInfo);
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
225           << "\n";*/
226       }
227     }
228   }
229   return h;
230 }
231
232
233 //________________________________________________________
234 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
235 {
236   if(track) fTrack = track;
237   if(!fTrack){
238     AliWarning("No Track defined.");
239     return 0x0;
240   }
241   TH1 *h = 0x0;
242   if(!(h = ((TH2I*)fContainer->At(kTrackletY)))){
243     AliWarning("No output histogram defined.");
244     return 0x0;
245   }
246
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]));
261   }
262   return h;
263 }
264
265 //________________________________________________________
266 TH1* AliTRDresolution::PlotTrackletPhi(const AliTRDtrackV1 *track)
267 {
268   if(track) fTrack = track;
269   if(!fTrack){
270     AliWarning("No Track defined.");
271     return 0x0;
272   }
273   TH1 *h = 0x0;
274   if(!(h = ((TH2I*)fContainer->At(kTrackletPhi)))){
275     AliWarning("No output histogram defined.");
276     return 0x0;
277   }
278
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)));
284   }
285   return h;
286 }
287
288
289 //________________________________________________________
290 TH1* AliTRDresolution::PlotResolution(const AliTRDtrackV1 *track)
291 {
292   if(!HasMCdata()){ 
293     AliWarning("No MC defined. Results will not be available.");
294     return 0x0;
295   }
296   if(track) fTrack = track;
297   if(!fTrack){
298     AliWarning("No Track defined.");
299     return 0x0;
300   }
301   TH1 *h = 0x0;
302   UChar_t s;
303   Int_t pdg = fMC->GetPDG(), det=-1;
304   Int_t label = fMC->GetLabel();
305   Double_t x, y, z;
306   Float_t p, pt, x0, y0, z0, dx, dy, dz, dydx, dzdx;
307   Double_t covR[3];
308
309   if(fDebugLevel>=1){
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"
313       << "pdg="  << pdg
314       << "dx0="  << DX[0]
315       << "dx1="  << DX[1]
316       << "dx2="  << DX[2]
317       << "dy0="  << DY[0]
318       << "dy1="  << DY[1]
319       << "dy2="  << DY[2]
320       << "dz0="  << DZ[0]
321       << "dz1="  << DZ[1]
322       << "dz2="  << DZ[2]
323       << "dpt0=" << DPt[0]
324       << "dpt1=" << DPt[1]
325       << "dpt2=" << DPt[2]
326       << "\n";
327   }
328
329   AliTRDseedV1 *fTracklet = 0x0;  
330   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
331     if(!(fTracklet = fTrack->GetTracklet(ily)) ||
332        !fTracklet->IsOK()) continue;
333
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     dx  = x0 - x;
339     if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, pt, s)) continue;
340     if(fDebugLevel>=1){
341       (*fDebugStream) << "MC"
342         << "det="     << det
343         << "pdg="     << pdg
344         << "pt="      << pt
345         << "x0="      << x0
346         << "y0="      << y0
347         << "z0="      << z0
348         << "dydx="    << dydx
349         << "dzdx="    << dzdx
350         << "\n";
351     }
352     // MC track position at reference radial position
353     Float_t yt = y0 - dx*dydx;
354     Float_t zt = z0 - dx*dzdx;
355     p = pt*(1.+dzdx*dzdx); // pt -> p
356
357     // add Kalman residuals for y, z and pt
358     Float_t yr = fTracklet->GetYref(0) - dx*fTracklet->GetYref(1);
359     dy = yt - yr;
360     Float_t zr = fTracklet->GetZref(0) - dx*fTracklet->GetZref(1);
361     dz = zt - zr;
362     Float_t tgl = fTracklet->GetTgl();
363     Float_t dpt = pt - fTracklet->GetMomentum()/(1.+tgl*tgl);
364     fTracklet->GetCovRef(covR);
365
366     ((TH2I*)fContainer->At(kMCtrackY))->Fill(dydx, dy);
367     ((TH2I*)fContainer->At(kMCtrackYPull))->Fill(dydx, dy/TMath::Sqrt(covR[0]));
368     if(ily==0){
369       ((TH2I*)fContainer->At(kMCtrackZIn))->Fill(dzdx, dz);
370       ((TH2I*)fContainer->At(kMCtrackZInPull))->Fill(dzdx, dz/TMath::Sqrt(covR[2]));
371     } else if(ily==AliTRDgeometry::kNlayer-1) {
372       ((TH2I*)fContainer->At(kMCtrackZOut))->Fill(dzdx, dz);
373       ((TH2I*)fContainer->At(kMCtrackZOutPull))->Fill(dzdx, dz/TMath::Sqrt(covR[2]));
374     }
375     if(pdg!=kElectron && pdg!=kPositron) ((TH2I*)fContainer->At(kMCtrackPt))->Fill(1./pt, dpt);
376     // Fill Debug stream for Kalman track
377     if(fDebugLevel>=1){
378       Float_t dydxr = fTracklet->GetYref(1);
379       (*fDebugStream) << "MCtrack"
380         << "x="       << x
381         << "dpt="     << dpt
382         << "dy="      << dy
383         << "dz="      << dz
384         << "dydxr="   << dydxr
385         << "dzdxr="   << tgl
386         << "s2y="     << covR[0]
387         << "s2z="     << covR[2]
388         << "\n";
389     }
390
391     // recalculate tracklet based on the MC info
392     AliTRDseedV1 tt(*fTracklet);
393     tt.SetZref(0, z0);
394     tt.SetZref(1, dzdx); 
395     tt.Fit(kTRUE);
396     x= tt.GetX(); // the true one 
397     dx = x0 - x;
398     yt = y0 - dx*dydx;
399     Bool_t rc = tt.IsRowCross(); 
400     
401     // add tracklet residuals for y and dydx
402     Float_t yf = tt.GetY();
403     dy = yt - yf; dz = 100.;
404     Float_t dphi   = (tt.GetYfit(1) - dydx);
405     dphi /= 1.- tt.GetYfit(1)*dydx;
406     Double_t s2y = tt.GetS2Y(), s2z = tt.GetS2Z();
407     if(!rc){
408       ((TH2I*)fContainer->At(kMCtrackletY))->Fill(dydx, dy);
409       ((TH2I*)fContainer->At(kMCtrackletYPull))->Fill(dydx, dy/TMath::Sqrt(s2y));
410       ((TH2I*)fContainer->At(kMCtrackletPhi))->Fill(dydx, dphi*TMath::RadToDeg());
411     } else {
412       // add tracklet residuals for z
413       dz = tt.GetZ() - (z0 - dx*dzdx) ;
414       ((TH2I*)fContainer->At(kMCtrackletZ))->Fill(dzdx, dz);
415       ((TH2I*)fContainer->At(kMCtrackletZPull))->Fill(dzdx, dz/TMath::Sqrt(s2z));
416     }
417   
418     // Fill Debug stream for tracklet
419     if(fDebugLevel>=1){
420       (*fDebugStream) << "MCtracklet"
421         << "rc="    << rc
422         << "x="     << x
423         << "dy="    << dy
424         << "dz="    << dz
425         << "dphi="  << dphi
426         << "s2y="   << s2y
427         << "s2z="   << s2z
428         << "\n";
429     }
430
431     Int_t istk = AliTRDgeometry::GetStack(det); 
432     AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
433     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
434     Float_t tilt = fTracklet->GetTilt();
435
436     AliTRDcluster *c = 0x0;
437     fTracklet->ResetClusterIter(kFALSE);
438     while((c = fTracklet->PrevCluster())){
439       Float_t  q = TMath::Abs(c->GetQ());
440       //AliTRDseedV1::GetClusterXY(c,x,y);
441       x = c->GetX(); y = c->GetY(); z = c->GetZ();
442       Float_t xc = x;
443       Float_t yc = y;
444       Float_t zc = z;
445       dx = x0 - xc; 
446       yt = y0 - dx*dydx;
447       zt = z0 - dx*dzdx;
448       dy = yt - (yc - tilt*(zc-zt));
449
450       // Fill Histograms
451       if(q>20. && q<250.) ((TH2I*)fContainer->At(kMCcluster))->Fill(dydx, dy);
452
453       // Fill calibration container
454       Float_t d = zr0 - zt;
455       d -= ((Int_t)(2 * d)) / 2.0;
456       if (d > 0.25) d  = 0.5 - d;
457       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
458       fMCcl->Add(clInfo);
459       clInfo->SetCluster(c);
460       clInfo->SetMC(pdg, label);
461       clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
462       clInfo->SetResolution(dy);
463       clInfo->SetAnisochronity(d);
464       clInfo->SetDriftLength(dx-.5*AliTRDgeometry::CamHght());
465       clInfo->SetTilt(tilt);
466
467       // Fill Debug Tree
468       if(fDebugLevel>=2){
469         //clInfo->Print();
470         (*fDebugStream) << "MCcluster"
471           <<"clInfo.=" << clInfo
472           << "\n";
473       }
474     }
475   }
476   return h;
477 }
478
479
480 //________________________________________________________
481 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
482 {
483   Float_t y[2] = {0., 0.};
484   TBox *b = 0x0;
485   TAxis *ax = 0x0;
486   TGraphErrors *g = 0x0;
487   switch(ifig){
488   case kCluster:
489     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
490     g->Draw("apl");
491     ax = g->GetHistogram()->GetYaxis();
492     y[0] = -0.5; y[1] = 2.5;
493     ax->SetRangeUser(y[0], y[1]);
494     ax->SetTitle("Cluster-Track Pools #sigma/#mu [mm]");
495     ax = g->GetHistogram()->GetXaxis();
496     ax->SetTitle("tg(#phi)");
497     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
498     g->Draw("pl");
499     b = new TBox(-.15, y[0], .15, y[1]);
500     b->SetFillStyle(3002);b->SetFillColor(kGreen);
501     b->SetLineColor(0); b->Draw();
502     return kTRUE;
503   case kTrackletY:
504     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
505     g->Draw("apl");
506     ax = g->GetHistogram()->GetYaxis();
507     ax->SetRangeUser(-.5, 3.);
508     ax->SetTitle("Tracklet-Track Y-Pools #sigma/#mu [mm]");
509     ax = g->GetHistogram()->GetXaxis();
510     ax->SetTitle("tg(#phi)");
511     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
512     g->Draw("pl");
513     b = new TBox(-.15, y[0], .15, y[1]);
514     b->SetFillStyle(3002);b->SetFillColor(kGreen);
515     b->SetLineColor(0); b->Draw();
516     return kTRUE;
517   case kTrackletPhi:
518     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
519     g->Draw("apl");
520     ax = g->GetHistogram()->GetYaxis();
521     ax->SetRangeUser(-.5, 2.);
522     ax->SetTitle("Tracklet-Track Phi-Residuals #sigma/#mu [rad]");
523     ax = g->GetHistogram()->GetXaxis();
524     ax->SetTitle("tg(#phi)");
525     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
526     g->Draw("pl");
527     b = new TBox(-.15, y[0], .15, y[1]);
528     b->SetFillStyle(3002);b->SetFillColor(kGreen);
529     b->SetLineColor(0); b->Draw();
530     return kTRUE;
531   case kMCcluster:
532     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
533     ax = g->GetHistogram()->GetYaxis();
534     y[0] = -.05; y[1] = 0.6;
535     ax->SetRangeUser(y[0], y[1]);
536     ax->SetTitle("Y_{cluster} #sigma/#mu [mm]");
537     ax = g->GetHistogram()->GetXaxis();
538     ax->SetRangeUser(-.3, .3);
539     ax->SetTitle("tg(#phi)");
540     g->Draw("apl");
541     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
542     g->Draw("pl");
543     b = new TBox(-.15, y[0], .15, y[1]);
544     b->SetFillStyle(3002);b->SetFillColor(kBlue);
545     b->SetLineColor(0); b->Draw();
546     return kTRUE;
547   case kMCtrackletY:
548     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
549     ax = g->GetHistogram()->GetYaxis();
550     y[0] = -.05; y[1] = 0.25;
551     ax->SetRangeUser(y[0], y[1]);
552     ax->SetTitle("Y_{tracklet} #sigma/#mu [mm]");
553     ax = g->GetHistogram()->GetXaxis();
554     ax->SetRangeUser(-.2, .2);
555     ax->SetTitle("tg(#phi)");
556     g->Draw("apl");
557     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
558     g->Draw("pl");
559     b = new TBox(-.15, y[0], .15, y[1]);
560     b->SetFillStyle(3002);b->SetFillColor(kBlue);
561     b->SetLineColor(0); b->Draw();
562     return kTRUE;
563   case kMCtrackletZ:
564     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
565     ax = g->GetHistogram()->GetYaxis();
566     ax->SetRangeUser(-.5, 1.);
567     ax->SetTitle("Z_{tracklet} #sigma/#mu [mm]");
568     ax = g->GetHistogram()->GetXaxis();
569     ax->SetTitle("tg(#theta)");
570     g->Draw("apl");
571     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
572     g->Draw("pl");
573     return kTRUE;
574   case kMCtrackletPhi:
575     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
576     ax = g->GetHistogram()->GetYaxis();
577     y[0] = -.05; y[1] = .2;
578     ax->SetRangeUser(y[0], y[1]);
579     ax->SetTitle("#Phi_{tracklet} #sigma/#mu [deg]");
580     ax = g->GetHistogram()->GetXaxis();
581     ax->SetTitle("tg(#phi)");
582     g->Draw("apl");
583     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
584     g->Draw("pl");
585     return kTRUE;
586   case kMCtrackY:
587     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
588     ax = g->GetHistogram()->GetYaxis();
589     y[0] = -.05; y[1] = 0.25;
590     ax->SetRangeUser(y[0], y[1]);
591     ax->SetTitle("Y_{track} #sigma/#mu [mm]");
592     ax = g->GetHistogram()->GetXaxis();
593     ax->SetRangeUser(-.2, .2);
594     ax->SetTitle("tg(#phi)");
595     g->Draw("apl");
596     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
597     g->Draw("pl");
598     b = new TBox(-.15, y[0], .15, y[1]);
599     b->SetFillStyle(3002);b->SetFillColor(kBlue);
600     b->SetLineColor(0); b->Draw();
601     return kTRUE;
602   case kMCtrackZIn:
603     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
604     ax = g->GetHistogram()->GetYaxis();
605     ax->SetRangeUser(-.5, 2.);
606     ax->SetTitle("Z_{track} #sigma/#mu [mm]");
607     ax = g->GetHistogram()->GetXaxis();
608     ax->SetTitle("tg(#theta)");
609     g->Draw("apl");
610     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
611     g->Draw("pl");
612     return kTRUE;
613   case kMCtrackPt:
614     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
615     ax = g->GetHistogram()->GetYaxis();
616     ax->SetRangeUser(-.5, 2.);
617     ax->SetTitle("#epsilon_{P_{t}}^{track} / #mu [%]");
618     ax = g->GetHistogram()->GetXaxis();
619     ax->SetTitle("1/p_{t}");
620     g->Draw("apl");
621     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
622     g->Draw("pl");
623     return kTRUE;
624   }
625   AliInfo(Form("Reference plot [%d] missing result", ifig));
626   return kFALSE;
627 }
628
629
630 //________________________________________________________
631 Bool_t AliTRDresolution::PostProcess()
632 {
633   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
634   if (!fContainer) {
635     Printf("ERROR: list not available");
636     return kFALSE;
637   }
638   fNRefFigures = fContainer->GetEntriesFast();
639   TGraphErrors *gm = 0x0, *gs = 0x0;
640   if(!fGraphS){ 
641     fGraphS = new TObjArray(fNRefFigures);
642     fGraphS->SetOwner();
643     for(Int_t ig=0; ig<fNRefFigures; ig++){
644       gs = new TGraphErrors();
645       gs->SetLineColor(kRed);
646       gs->SetMarkerStyle(23);
647       gs->SetMarkerColor(kRed);
648       gs->SetNameTitle(Form("s_%d", ig), "");
649       fGraphS->AddAt(gs, ig);
650     }
651   }
652   if(!fGraphM){ 
653     fGraphM = new TObjArray(fNRefFigures);
654     fGraphM->SetOwner();
655     for(Int_t ig=0; ig<fNRefFigures; ig++){
656       gm = new TGraphErrors();
657       gm->SetLineColor(kBlue);
658       gm->SetMarkerStyle(7);
659       gm->SetMarkerColor(kBlue);
660       gm->SetNameTitle(Form("m_%d", ig), "");
661       fGraphM->AddAt(gm, ig);
662     }
663   }
664
665   TH2I *h2 = 0x0;
666   TH1D *h = 0x0;
667
668   // define models
669   TF1 f("f1", "gaus", -.5, .5);  
670
671   TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
672
673   TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
674
675   TCanvas *c = 0x0;
676   if(IsVisual()) c = new TCanvas("c", Form("%s Visual", GetName()), 500, 500);
677   char opt[5];
678   sprintf(opt, "%c%c", IsVerbose() ? ' ' : 'Q', IsVisual() ? ' ': 'N');
679
680
681   //PROCESS RESIDUAL DISTRIBUTIONS
682
683   // Clusters residuals
684   h2 = (TH2I *)(fContainer->At(kCluster));
685   gm = (TGraphErrors*)fGraphM->At(kCluster);
686   gs = (TGraphErrors*)fGraphS->At(kCluster);
687   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
688     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
689     h = h2->ProjectionY("py", ibin, ibin);
690     if(h->GetEntries()<100) continue;
691     AdjustF1(h, &f);
692
693     if(IsVisual()){c->cd(); c->SetLogy();}
694     h->Fit(&f, opt, "", -0.5, 0.5);
695     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
696     
697     Int_t ip = gm->GetN();
698     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
699     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
700     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
701     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
702   }
703
704   // Tracklet y residuals
705   h2 = (TH2I *)(fContainer->At(kTrackletY));
706   gm = (TGraphErrors*)fGraphM->At(kTrackletY);
707   gs = (TGraphErrors*)fGraphS->At(kTrackletY);
708   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
709     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
710     h = h2->ProjectionY("py", ibin, ibin);
711     if(h->GetEntries()<100) continue;
712     AdjustF1(h, &f);
713
714     if(IsVisual()){c->cd(); c->SetLogy();}
715     h->Fit(&f, opt, "", -0.5, 0.5);
716     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
717     
718     Int_t ip = gm->GetN();
719     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
720     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
721     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
722     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
723   }
724
725   // Tracklet phi residuals
726   h2 = (TH2I *)(fContainer->At(kTrackletPhi));
727   gm = (TGraphErrors*)fGraphM->At(kTrackletPhi);
728   gs = (TGraphErrors*)fGraphS->At(kTrackletPhi);
729   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
730     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
731     h = h2->ProjectionY("py", ibin, ibin);
732     if(h->GetEntries()<100) continue;
733     AdjustF1(h, &f);
734
735     if(IsVisual()){c->cd(); c->SetLogy();}
736     h->Fit(&f, opt, "", -0.5, 0.5);
737     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
738     
739     Int_t ip = gm->GetN();
740     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
741     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
742     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
743     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
744   }
745
746   if(!HasMCdata()){
747     if(c) delete c;
748     return kTRUE;
749   }
750
751   //PROCESS MC RESIDUAL DISTRIBUTIONS
752
753   // cluster y resolution
754   h2 = (TH2I*)fContainer->At(kMCcluster);
755   gm = (TGraphErrors*)fGraphM->At(kMCcluster);
756   gs = (TGraphErrors*)fGraphS->At(kMCcluster);
757   for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
758     h = h2->ProjectionY("py", iphi, iphi);
759     if(h->GetEntries()<100) continue;
760     AdjustF1(h, &f);
761
762     if(IsVisual()){c->cd(); c->SetLogy();}
763     h->Fit(&f, opt, "", -0.5, 0.5);
764     if(IsVerbose()){
765       printf("phi[%d] mean[%e] sigma[%e]\n\n", iphi, 10.*f.GetParameter(1), 10.*f.GetParameter(2));
766     }
767     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
768
769     Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
770     Int_t ip = gm->GetN();
771     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
772     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
773     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
774     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
775   }
776
777   // tracklet y resolution
778   h2 = (TH2I*)fContainer->At(kMCtrackletY);
779   gm = (TGraphErrors*)fGraphM->At(kMCtrackletY);
780   gs = (TGraphErrors*)fGraphS->At(kMCtrackletY);
781   for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
782     h = h2->ProjectionY("py", iphi, iphi);
783     if(h->GetEntries()<100) continue;
784     AdjustF1(h, &f);
785
786     if(IsVisual()){c->cd(); c->SetLogy();}
787     h->Fit(&f, opt, "", -0.5, 0.5);
788     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
789
790     Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
791     Int_t ip = gm->GetN();
792     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
793     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
794     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
795     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
796   }
797
798   // tracklet z resolution
799   h2 = (TH2I*)fContainer->At(kMCtrackletZ);
800   gm = (TGraphErrors*)fGraphM->At(kMCtrackletZ);
801   gs = (TGraphErrors*)fGraphS->At(kMCtrackletZ);
802   for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
803     h = h2->ProjectionY("py", iphi, iphi);
804     if(h->GetEntries()<100) continue;
805     AdjustF1(h, &fb);
806
807     if(IsVisual()){c->cd(); c->SetLogy();}
808     h->Fit(&fb, opt, "", -0.5, 0.5);
809     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
810
811     Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
812     Int_t ip = gm->GetN();
813     gm->SetPoint(ip, phi, 10.*fb.GetParameter(1));
814     gm->SetPointError(ip, 0., 10.*fb.GetParError(1));
815     gs->SetPoint(ip, phi, 10.*fb.GetParameter(2));
816     gs->SetPointError(ip, 0., 10.*fb.GetParError(2));
817   }
818
819   //tracklet phi resolution
820   h2 = (TH2I*)fContainer->At(kMCtrackletPhi);
821   gm = (TGraphErrors*)fGraphM->At(kMCtrackletPhi);
822   gs = (TGraphErrors*)fGraphS->At(kMCtrackletPhi);
823   for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
824     h = h2->ProjectionY("py", iphi, iphi);
825     if(h->GetEntries()<100) continue;
826
827     if(IsVisual()){c->cd(); c->SetLogy();}
828     h->Fit(&f, opt, "", -0.5, 0.5);
829     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
830
831     Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
832     Int_t ip = gm->GetN();
833     gm->SetPoint(ip, phi, f.GetParameter(1));
834     gm->SetPointError(ip, 0., f.GetParError(1));
835     gs->SetPoint(ip, phi, f.GetParameter(2));
836     gs->SetPointError(ip, 0., f.GetParError(2));
837   }
838
839   // track y resolution
840   h2 = (TH2I*)fContainer->At(kMCtrackY);
841   gm = (TGraphErrors*)fGraphM->At(kMCtrackY);
842   gs = (TGraphErrors*)fGraphS->At(kMCtrackY);
843   for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
844     h = h2->ProjectionY("py", iphi, iphi);
845     if(h->GetEntries()<100) continue;
846     AdjustF1(h, &f);
847
848     if(IsVisual()){c->cd(); c->SetLogy();}
849     h->Fit(&f, opt, "", -0.5, 0.5);
850     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
851
852     Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
853     Int_t ip = gm->GetN();
854     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
855     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
856     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
857     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
858   }
859
860   // track z resolution
861   h2 = (TH2I*)fContainer->At(kMCtrackZIn);
862   gm = (TGraphErrors*)fGraphM->At(kMCtrackZIn);
863   gs = (TGraphErrors*)fGraphS->At(kMCtrackZIn);
864   for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
865     h = h2->ProjectionY("pz", iphi, iphi);
866     if(h->GetEntries()<70) continue;
867     AdjustF1(h, &f);
868
869     if(IsVisual()){c->cd(); c->SetLogy();}
870     h->Fit(&f, opt, "", -0.5, 0.5);
871     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
872
873     Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
874     Int_t ip = gm->GetN();
875     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
876     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
877     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
878     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
879   }
880
881   // track Pt resolution
882   h2 = (TH2I*)fContainer->At(kMCtrackPt);
883   TAxis *ax = h2->GetXaxis();
884   gm = (TGraphErrors*)fGraphM->At(kMCtrackPt);
885   gs = (TGraphErrors*)fGraphS->At(kMCtrackPt);
886   TF1 fg("fg", "gaus", -1.5, 1.5);
887   TF1 fl("fl", "landau", -4., 15.);
888   TF1 fgl("fgl", "gaus(0)+landau(3)", -5., 20.);
889   for(Int_t ip=1; ip<=ax->GetNbins(); ip++){
890     h = h2->ProjectionY("ppt", ip, ip);
891     if(h->GetEntries()<70) continue;
892
893     h->Fit(&fg, "QN", "", -1.5, 1.5);
894     fgl.SetParameter(0, fg.GetParameter(0));
895     fgl.SetParameter(1, fg.GetParameter(1));
896     fgl.SetParameter(2, fg.GetParameter(2));
897     h->Fit(&fl, "QN", "", -4., 15.);
898     fgl.SetParameter(3, fl.GetParameter(0));
899     fgl.SetParameter(4, fl.GetParameter(1));
900     fgl.SetParameter(5, fl.GetParameter(2));
901
902     if(IsVisual()){c->cd(); c->SetLogy();}
903     h->Fit(&fgl, opt, "", -5., 20.);
904     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
905
906     Float_t invpt = ax->GetBinCenter(ip);
907     Int_t ip = gm->GetN();
908     gm->SetPoint(ip, invpt, fgl.GetParameter(1));
909     gm->SetPointError(ip, 0., fgl.GetParError(1));
910     gs->SetPoint(ip, invpt, fgl.GetParameter(2)*invpt);
911     gs->SetPointError(ip, 0., fgl.GetParError(2));
912     // fgl.GetParameter(4) // Landau MPV
913     // fgl.GetParameter(5) // Landau Sigma
914   }
915
916
917   if(c) delete c;
918   return kTRUE;
919 }
920
921
922 //________________________________________________________
923 void AliTRDresolution::Terminate(Option_t *)
924 {
925   if(fDebugStream){ 
926     delete fDebugStream;
927     fDebugStream = 0x0;
928     fDebugLevel = 0;
929   }
930   if(HasPostProcess()) PostProcess();
931 }
932
933 //________________________________________________________
934 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
935 {
936 // Helper function to avoid duplication of code
937 // Make first guesses on the fit parameters
938
939   // find the intial parameters of the fit !! (thanks George)
940   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
941   Double_t sum = 0.;
942   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
943   f->SetParLimits(0, 0., 3.*sum);
944   f->SetParameter(0, .9*sum);
945
946   f->SetParLimits(1, -.2, .2);
947   f->SetParameter(1, -0.1);
948
949   f->SetParLimits(2, 0., 4.e-1);
950   f->SetParameter(2, 2.e-2);
951   if(f->GetNpar() <= 4) return;
952
953   f->SetParLimits(3, 0., sum);
954   f->SetParameter(3, .1*sum);
955
956   f->SetParLimits(4, -.3, .3);
957   f->SetParameter(4, 0.);
958
959   f->SetParLimits(5, 0., 1.e2);
960   f->SetParameter(5, 2.e-1);
961 }
962
963 //________________________________________________________
964 TObjArray* AliTRDresolution::Histos()
965 {
966   if(fContainer) return fContainer;
967
968   fContainer  = new TObjArray(kNhistos);
969   //fContainer->SetOwner(kTRUE);
970
971   TH1 *h = 0x0;
972   // cluster to tracklet residuals [2]
973   if(!(h = (TH2I*)gROOT->FindObject("hCl"))){
974     h = new TH2I("hCl", "Clusters-Track Residuals", 21, -.33, .33, 100, -.5, .5);
975     h->GetXaxis()->SetTitle("tg(#phi)");
976     h->GetYaxis()->SetTitle("#Delta y [cm]");
977     h->GetZaxis()->SetTitle("entries");
978   } else h->Reset();
979   fContainer->AddAt(h, kCluster);
980
981   // tracklet to track residuals [2]
982   if(!(h = (TH2I*)gROOT->FindObject("hTrkltY"))){
983     h = new TH2I("hTrkltY", "Tracklets-Track Residuals (Y)", 21, -.33, .33, 100, -.5, .5);
984     h->GetXaxis()->SetTitle("tg(#phi)");
985     h->GetYaxis()->SetTitle("#Delta y [cm]");
986     h->GetZaxis()->SetTitle("entries");
987   } else h->Reset();
988   fContainer->AddAt(h, kTrackletY);
989
990   // tracklet to track residuals angular [2]
991   if(!(h = (TH2I*)gROOT->FindObject("hTrkltPhi"))){
992     h = new TH2I("hTrkltPhi", "Tracklets-Track Residuals (#Phi)", 21, -.33, .33, 100, -.5, .5);
993     h->GetXaxis()->SetTitle("tg(#phi)");
994     h->GetYaxis()->SetTitle("#Delta phi [#circ]");
995     h->GetZaxis()->SetTitle("entries");
996   } else h->Reset();
997   fContainer->AddAt(h, kTrackletPhi);
998
999
1000   // Resolution histos
1001   if(!HasMCdata()) return fContainer;
1002
1003   // cluster y resolution [0]
1004   if(!(h = (TH2I*)gROOT->FindObject("hMCcl"))){
1005     h = new TH2I("hMCcl", "Cluster Resolution", 31, -.48, .48, 100, -.5, .5);
1006     h->GetXaxis()->SetTitle("tg(#phi)");
1007     h->GetYaxis()->SetTitle("#Delta y [cm]");
1008     h->GetZaxis()->SetTitle("entries");
1009   } else h->Reset();
1010   fContainer->AddAt(h, kMCcluster);
1011
1012   // tracklet y resolution [0]
1013   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltY"))){
1014     h = new TH2I("hMCtrkltY", "Tracklet Resolution (Y)", 31, -.48, .48, 100, -.5, .5);
1015     h->GetXaxis()->SetTitle("tg(#phi)");
1016     h->GetYaxis()->SetTitle("#Delta y [cm]");
1017     h->GetZaxis()->SetTitle("entries");
1018   } else h->Reset();
1019   fContainer->AddAt(h, kMCtrackletY);
1020
1021   // tracklet y resolution [0]
1022   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltYPull"))){
1023     h = new TH2I("hMCtrkltYPull", "Tracklet Pulls (Y)", 31, -.48, .48, 100, -3.5, 3.5);
1024     h->GetXaxis()->SetTitle("tg(#phi)");
1025     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1026     h->GetZaxis()->SetTitle("entries");
1027   } else h->Reset();
1028   fContainer->AddAt(h, kMCtrackletYPull);
1029
1030   // tracklet y resolution [0]
1031   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltZ"))){
1032     h = new TH2I("hMCtrkltZ", "Tracklet Resolution (Z)", 31, -.48, .48, 100, -.5, .5);
1033     h->GetXaxis()->SetTitle("tg(#theta)");
1034     h->GetYaxis()->SetTitle("#Delta z [cm]");
1035     h->GetZaxis()->SetTitle("entries");
1036   } else h->Reset();
1037   fContainer->AddAt(h, kMCtrackletZ);
1038
1039   // tracklet y resolution [0]
1040   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltZPull"))){
1041     h = new TH2I("hMCtrkltZ", "Tracklet Pulls (Z)", 31, -.48, .48, 100, -3.5, 3.5);
1042     h->GetXaxis()->SetTitle("tg(#theta)");
1043     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1044     h->GetZaxis()->SetTitle("entries");
1045   } else h->Reset();
1046   fContainer->AddAt(h, kMCtrackletZPull);
1047
1048   // tracklet angular resolution [1]
1049   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkltPhi"))){
1050     h = new TH2I("hMCtrkltPhi", "Tracklet Resolution (#Phi)", 31, -.48, .48, 100, -10., 10.);
1051     h->GetXaxis()->SetTitle("tg(#phi)");
1052     h->GetYaxis()->SetTitle("#Delta #phi [deg]");
1053     h->GetZaxis()->SetTitle("entries");
1054   } else h->Reset();
1055   fContainer->AddAt(h, kMCtrackletPhi);
1056
1057   // Kalman track y resolution
1058   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkY"))){
1059     h = new TH2I("hMCtrkY", "Kalman Track Resolution (Y)", 31, -.48, .48, 100, -.5, .5);
1060     h->GetXaxis()->SetTitle("tg(#phi)");
1061     h->GetYaxis()->SetTitle("#Delta y [cm]");
1062     h->GetZaxis()->SetTitle("entries");
1063   } else h->Reset();
1064   fContainer->AddAt(h, kMCtrackY);
1065
1066   // Kalman track y resolution
1067   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkYPull"))){
1068     h = new TH2I("hMCtrkYPull", "Kalman Track Pulls (Y)", 31, -.48, .48, 100, -3.5, 3.5);
1069     h->GetXaxis()->SetTitle("tg(#phi)");
1070     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1071     h->GetZaxis()->SetTitle("entries");
1072   } else h->Reset();
1073   fContainer->AddAt(h, kMCtrackYPull);
1074
1075   // Kalman track Z resolution
1076   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZIn"))){
1077     h = new TH2I("hMCtrkZIn", "Kalman Track Resolution (Zin)", 20, -1., 1., 100, -1.5, 1.5);
1078     h->GetXaxis()->SetTitle("tg(#theta)");
1079     h->GetYaxis()->SetTitle("#Delta z [cm]");
1080     h->GetZaxis()->SetTitle("entries");
1081   } else h->Reset();
1082   fContainer->AddAt(h, kMCtrackZIn);
1083
1084   // Kalman track Z resolution
1085   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZOut"))){
1086     h = new TH2I("hMCtrkZOut", "Kalman Track Resolution (Zout)", 20, -1., 1., 100, -1.5, 1.5);
1087     h->GetXaxis()->SetTitle("tg(#theta)");
1088     h->GetYaxis()->SetTitle("#Delta z [cm]");
1089     h->GetZaxis()->SetTitle("entries");
1090   } else h->Reset();
1091   fContainer->AddAt(h, kMCtrackZOut);
1092
1093   // Kalman track Z resolution
1094   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZInPull"))){
1095     h = new TH2I("hMCtrkZInPull", "Kalman Track Pulls (Zin)", 20, -1., 1., 100, -4.5, 4.5);
1096     h->GetXaxis()->SetTitle("tg(#theta)");
1097     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1098     h->GetZaxis()->SetTitle("entries");
1099   } else h->Reset();
1100   fContainer->AddAt(h, kMCtrackZInPull);
1101
1102   // Kalman track Z resolution
1103   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkZOutPull"))){
1104     h = new TH2I("hMCtrkZOutPull", "Kalman Track Pulls (Zout)", 20, -1., 1., 100, -4.5, 4.5);
1105     h->GetXaxis()->SetTitle("tg(#theta)");
1106     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1107     h->GetZaxis()->SetTitle("entries");
1108   } else h->Reset();
1109   fContainer->AddAt(h, kMCtrackZOutPull);
1110
1111   // Kalman track Pt resolution
1112   if(!(h = (TH2I*)gROOT->FindObject("hMCtrkPt"))){
1113     h = new TH2I("hMCtrkPt", "Kalman Track Resolution (Pt)", 100, 0., 2., 150, -5., 20.);
1114     h->GetXaxis()->SetTitle("1/p_{t}");
1115     h->GetYaxis()->SetTitle("#Delta p_{t} [GeV/c]");
1116     h->GetZaxis()->SetTitle("entries");
1117   } else h->Reset();
1118   fContainer->AddAt(h, kMCtrackPt);
1119
1120   return fContainer;
1121 }
1122
1123
1124 //________________________________________________________
1125 void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
1126 {
1127
1128   fReconstructor->SetRecoParam(r);
1129 }