fix error - refit tracklets using tilt correction
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDtrackingResolution.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: AliTRDtrackingResolution.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 //   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);
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 <TObjArray.h>
56 #include <TH2.h>
57 #include <TH1.h>
58 #include <TF1.h>
59 #include <TCanvas.h>
60 #include <TBox.h>
61 #include <TProfile.h>
62 #include <TGraphErrors.h>
63 #include <TMath.h>
64 #include "TTreeStream.h"
65 #include "TGeoManager.h"
66
67 #include "AliAnalysisManager.h"
68 #include "AliTrackReference.h"
69 #include "AliTrackPointArray.h"
70 #include "AliCDBManager.h"
71
72 #include "AliTRDSimParam.h"
73 #include "AliTRDgeometry.h"
74 #include "AliTRDpadPlane.h"
75 #include "AliTRDcluster.h"
76 #include "AliTRDseedV1.h"
77 #include "AliTRDtrackV1.h"
78 #include "AliTRDtrackerV1.h"
79 #include "AliTRDReconstructor.h"
80 #include "AliTRDrecoParam.h"
81
82 #include "AliTRDtrackInfo/AliTRDclusterInfo.h"
83 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
84 #include "AliTRDtrackingResolution.h"
85
86 ClassImp(AliTRDtrackingResolution)
87
88 //________________________________________________________
89 AliTRDtrackingResolution::AliTRDtrackingResolution()
90   :AliTRDrecoTask("Resolution", "Tracking Resolution")
91   ,fStatus(0)
92   ,fReconstructor(0x0)
93   ,fGeo(0x0)
94   ,fGraphS(0x0)
95   ,fGraphM(0x0)
96   ,fClResiduals(0x0)
97   ,fTrkltResiduals(0x0)
98   ,fTrkltPhiResiduals(0x0)
99   ,fClResolution(0x0)
100   ,fTrkltResolution(0x0)
101 {
102   fReconstructor = new AliTRDReconstructor();
103   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
104   fGeo = new AliTRDgeometry();
105
106   InitFunctorList();
107
108   DefineOutput(1+kClusterResidual, TObjArray::Class());
109   DefineOutput(1+kTrackletYResidual, TObjArray::Class());
110   DefineOutput(1+kTrackletPhiResidual, TObjArray::Class());
111   DefineOutput(1+kClusterResolution, TObjArray::Class());
112   DefineOutput(1+kTrackletYResolution, TObjArray::Class());
113 }
114
115 //________________________________________________________
116 AliTRDtrackingResolution::~AliTRDtrackingResolution()
117 {
118   if(fGraphS){fGraphS->Delete(); delete fGraphS;}
119   if(fGraphM){fGraphM->Delete(); delete fGraphM;}
120   delete fGeo;
121   delete fReconstructor;
122   if(gGeoManager) delete gGeoManager;
123   if(fClResiduals){fClResiduals->Delete(); delete fClResiduals;}
124   if(fTrkltResiduals){fTrkltResiduals->Delete(); delete fTrkltResiduals;}
125   if(fTrkltPhiResiduals){fTrkltPhiResiduals->Delete(); delete fTrkltPhiResiduals;}
126   if(fClResolution){fClResolution->Delete(); delete fClResolution;}
127   if(fTrkltResolution){fTrkltResolution->Delete(); delete fTrkltResolution;}
128 }
129
130
131 //________________________________________________________
132 void AliTRDtrackingResolution::CreateOutputObjects()
133 {
134   // spatial resolution
135   OpenFile(0, "RECREATE");
136
137   fContainer = Histos();
138
139   fClResiduals = new TObjArray();
140   fClResiduals->SetOwner(kTRUE);
141   fTrkltResiduals = new TObjArray();
142   fTrkltResiduals->SetOwner(kTRUE);
143   fTrkltPhiResiduals = new TObjArray();
144   fTrkltPhiResiduals->SetOwner(kTRUE);
145   fClResolution = new TObjArray();
146   fClResolution->SetOwner(kTRUE);
147   fTrkltResolution = new TObjArray();
148   fTrkltResolution->SetOwner(kTRUE);
149 }
150
151 //________________________________________________________
152 void AliTRDtrackingResolution::Exec(Option_t *opt)
153 {
154   fClResiduals->Delete();
155   fTrkltResiduals->Delete();
156   fTrkltPhiResiduals->Delete();
157   fClResolution->Delete();
158   fTrkltResolution->Delete();
159
160   AliTRDrecoTask::Exec(opt);
161
162   PostData(1+kClusterResidual, fClResiduals);
163   PostData(1+kTrackletYResidual, fTrkltResiduals);
164   PostData(1+kTrackletPhiResidual, fTrkltPhiResiduals);
165   PostData(1+kClusterResolution, fClResolution);
166   PostData(1+kTrackletYResolution, fTrkltResolution);
167 }
168
169 //________________________________________________________
170 TH1* AliTRDtrackingResolution::PlotClusterResiduals(const AliTRDtrackV1 *track)
171 {
172   if(track) fTrack = track;
173   if(!fTrack){
174     AliWarning("No Track defined.");
175     return 0x0;
176   }
177   TH1 *h = 0x0;
178   if(!(h = ((TH2I*)fContainer->At(kClusterResidual)))){
179     AliWarning("No output histogram defined.");
180     return 0x0;
181   }
182
183   Float_t x0, y0, z0, dy, dydx, dzdx;
184   AliTRDseedV1 *fTracklet = 0x0;  
185   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
186     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
187     if(!fTracklet->IsOK()) continue;
188     x0 = fTracklet->GetX0();
189
190     // retrive the track angle with the chamber
191     y0   = fTracklet->GetYref(0);
192     z0   = fTracklet->GetZref(0);
193     dydx = fTracklet->GetYref(1);
194     dzdx = fTracklet->GetZref(1);
195     Float_t tilt = fTracklet->GetTilt();
196     AliTRDcluster *c = 0x0;
197     fTracklet->ResetClusterIter(kFALSE);
198     while((c = fTracklet->PrevCluster())){
199       Float_t xc = c->GetX();
200       Float_t yc = c->GetY();
201       Float_t zc = c->GetZ();
202       Float_t dx = x0 - xc; 
203       Float_t yt = y0 - dx*dydx;
204       Float_t zt = z0 - dx*dzdx; 
205       dy = yt - (yc - tilt*(zc-zt));
206
207       //dy = trklt.GetYat(c->GetX()) - c->GetY();
208       h->Fill(dydx, dy);
209   
210       if(fDebugLevel>=1){
211         // Get z-position with respect to anode wire
212         AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
213         Int_t istk = fGeo->GetStack(c->GetDetector());
214         AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
215         Float_t row0 = pp->GetRow0();
216         Float_t d  =  row0 - zt + simParam->GetAnodeWireOffset();
217         d -= ((Int_t)(2 * d)) / 2.0;
218         if (d > 0.25) d  = 0.5 - d;
219
220         AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
221         fClResiduals->Add(clInfo);
222         clInfo->SetCluster(c);
223         clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
224         clInfo->SetResolution(dy);
225         clInfo->SetAnisochronity(d);
226         clInfo->SetDriftLength(dx);
227         (*fDebugStream) << "ClusterResiduals"
228           <<"clInfo.=" << clInfo
229           << "\n";
230       }
231     }
232   }
233   return h;
234 }
235
236
237 //________________________________________________________
238 TH1* AliTRDtrackingResolution::PlotTrackletResiduals(const AliTRDtrackV1 *track)
239 {
240   if(track) fTrack = track;
241   if(!fTrack){
242     AliWarning("No Track defined.");
243     return 0x0;
244   }
245   TH1 *h = 0x0;
246   if(!(h = ((TH2I*)fContainer->At(kTrackletYResidual)))){
247     AliWarning("No output histogram defined.");
248     return 0x0;
249   }
250
251   // refit the track
252   AliRieman fRim(fTrack->GetNumberOfClusters());
253   Float_t x[AliTRDgeometry::kNlayer] = {-1., -1., -1., -1., -1., -1.}, y[AliTRDgeometry::kNlayer], dydx[AliTRDgeometry::kNlayer];
254   AliTRDseedV1 *tracklet = 0x0;
255   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
256     if(!(tracklet = fTrack->GetTracklet(il))) continue;
257     if(!tracklet->IsOK()) continue;
258     AliTRDcluster *c = 0x0;
259     tracklet->ResetClusterIter(kFALSE);
260     while((c = tracklet->PrevCluster())){
261       Float_t xc = c->GetX();
262       Float_t yc = c->GetY();
263       Float_t zc = c->GetZ();
264       Float_t zt = tracklet->GetZref(0) - (tracklet->GetX0()-xc)*tracklet->GetZref(1); 
265       yc -= tracklet->GetTilt()*(zc-zt);
266       fRim.AddPoint(xc, yc, zc, 1, 10);
267     }
268     tracklet->Fit(kTRUE);
269
270     x[il] = tracklet->GetX0();
271     y[il] = tracklet->GetYfit(0)-tracklet->GetYfit(1)*(tracklet->GetX0()-x[il]);
272     dydx[il] = tracklet->GetYref(1);
273   }
274   fRim.Update();
275
276   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
277     if(x[il] < 0.) continue;
278     Float_t dy = y[il]-fRim.GetYat(x[il])/*/sigma_track*/;
279     h->Fill(dydx[il], dy);
280
281     if(fDebugLevel>=1){
282       Double_t sigmay = fRim.GetErrY(x[il]);
283       Float_t yt = fRim.GetYat(x[il]);
284       (*fDebugStream) << "TrkltResiduals"
285         << "layer="  << il
286         << "x="      << x[il]
287         << "y="      << y[il]
288         << "yt="     << yt
289         << "dydx="   << dydx[il]
290         << "dy="     << dy
291         << "sigmay=" << sigmay
292         << "\n";
293     }
294   }
295   return h;
296 }
297
298 //________________________________________________________
299 TH1* AliTRDtrackingResolution::PlotTrackletPhiResiduals(const AliTRDtrackV1 *track)
300 {
301   if(track) fTrack = track;
302   if(!fTrack){
303     AliWarning("No Track defined.");
304     return 0x0;
305   }
306   TH1 *h = 0x0;
307   if(!(h = ((TH2I*)fContainer->At(kTrackletPhiResidual)))){
308     AliWarning("No output histogram defined.");
309     return 0x0;
310   }
311   // refit the track
312   AliRieman fRim(fTrack->GetNumberOfClusters());
313   Float_t x[AliTRDgeometry::kNlayer] = {-1., -1., -1., -1., -1., -1.}, y[AliTRDgeometry::kNlayer], dydx[AliTRDgeometry::kNlayer];
314   Float_t dydx_ref=0, dydx_fit=0, phiref=0, phifit=0, phidiff=0;
315   AliTRDseedV1 *tracklet = 0x0;
316   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
317     if(!(tracklet = fTrack->GetTracklet(il))) continue;
318     if(!tracklet->IsOK()) continue;
319     AliTRDcluster *c = 0x0;
320     tracklet->ResetClusterIter(kFALSE);
321     while((c = tracklet->PrevCluster())){
322       Float_t xc = c->GetX();
323       Float_t yc = c->GetY();
324       Float_t zc = c->GetZ();
325       Float_t zt = tracklet->GetZref(0) - (tracklet->GetX0()-xc)*tracklet->GetZref(1); 
326       yc -= tracklet->GetTilt()*(zc-zt);
327       fRim.AddPoint(xc, yc, zc, 1, 10);
328     }
329     tracklet->Fit(kTRUE);
330
331     x[il] = tracklet->GetX0();
332     y[il] = tracklet->GetYfit(0)-tracklet->GetYfit(1)*(tracklet->GetX0()-x[il]);
333     dydx[il] = tracklet->GetYref(1);
334   }
335   fRim.Update();
336
337   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
338     if(x[il] < 0.) continue;
339     if(!(tracklet = fTrack->GetTracklet(il))) continue;
340     if(!tracklet->IsOK()) continue;
341
342     dydx_ref = fRim.GetDYat(x[il]); 
343     dydx_fit = tracklet->GetYfit(1);
344
345     phiref = TMath::ATan(dydx_ref);//*TMath::RadToDeg();
346     phifit = TMath::ATan(dydx_fit);//*TMath::RadToDeg();
347     
348     phidiff = phiref-phifit; /*/sigma_phi*/;
349
350     h->Fill(dydx_ref, phidiff);
351
352
353     if(fDebugLevel>=1){
354       (*fDebugStream) << "TrkltPhiResiduals"
355         << "layer="  << il
356         << "dydx_fit="   << dydx_fit
357         << "dydx_ref="   << dydx_ref
358         << "phiref="     << phiref
359         << "phifit="     << phifit
360         << "phidiff="     << phidiff
361         << "\n";
362     }
363   }
364   return h;
365 }
366
367
368 //________________________________________________________
369 TH1* AliTRDtrackingResolution::PlotResolution(const AliTRDtrackV1 *track)
370 {
371   if(!HasMCdata()){ 
372     AliWarning("No MC defined. Results will not be available.");
373     return 0x0;
374   }
375   if(track) fTrack = track;
376   if(!fTrack){
377     AliWarning("No Track defined.");
378     return 0x0;
379   }
380   TH1 *h = 0x0;
381   if(!(h = ((TH2I*)fContainer->At(kClusterResolution)))){
382     AliWarning("No output histogram defined.");
383     return 0x0;
384   }
385   if(!(h = ((TH2I*)fContainer->At(kTrackletYResolution)))){
386     AliWarning("No output histogram defined.");
387     return 0x0;
388   }
389   if(!(h = ((TH2I*)fContainer->At(kTrackletZResolution)))){
390     AliWarning("No output histogram defined.");
391     return 0x0;
392   }
393   if(!(h = ((TH2I*)fContainer->At(kTrackletAngleResolution)))){
394     AliWarning("No output histogram defined.");
395     return 0x0;
396   }
397   //printf("called for %d tracklets ... \n", fTrack->GetNumberOfTracklets());
398   UChar_t s;
399   Int_t pdg = fMC->GetPDG(), det=-1;
400   Int_t label = fMC->GetLabel();
401   Float_t x0, y0, z0, dx, dy, dydx, dzdx;
402   AliTRDseedV1 *fTracklet = 0x0;  
403   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
404     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
405     if(!fTracklet->IsOK()) continue;
406     //printf("process layer[%d]\n", ily);
407     // retrive the track position and direction within the chamber
408     det = fTracklet->GetDetector();
409     x0  = fTracklet->GetX0();
410     if(!fMC->GetDirections(x0, y0, z0, dydx, dzdx, s)) continue; 
411
412     // recalculate tracklet based on the MC info
413     AliTRDseedV1 tt(*fTracklet);
414     tt.SetZref(0, z0);
415     tt.SetZref(1, dzdx); 
416     if(!tt.Fit(kTRUE)) continue;
417     //tt.Update();
418
419     dx = 0.;//x0 - tt.GetXref();
420     Float_t yt = y0 - dx*dydx;
421     Float_t yf = tt.GetYfit(0) - dx*tt.GetYfit(1);
422     dy = yf-yt;
423     Float_t dphi   = TMath::ATan(tt.GetYfit(1)) - TMath::ATan(dydx);
424     Float_t dz = 100.;
425     Bool_t  cross = fTracklet->GetNChange();
426     if(cross){
427       Double_t *xyz = tt.GetCrossXYZ();
428       dz = xyz[2] - (z0 - (x0 - xyz[0])*dzdx) ;
429       ((TH2I*)fContainer->At(kTrackletZResolution))->Fill(dzdx, dz);
430     }
431   
432     // Fill Histograms
433     ((TH2I*)fContainer->At(kTrackletYResolution))->Fill(dydx, dy);
434     ((TH2I*)fContainer->At(kTrackletAngleResolution))->Fill(dydx, dphi*TMath::RadToDeg());
435
436     // Fill Debug stream
437     if(fDebugLevel>=1){
438       Float_t p = fMC->GetTrackRef() ? fMC->GetTrackRef()->P() : -1.;
439       (*fDebugStream) << "TrkltResolution"
440         << "det="                 << det
441         << "pdg="     << pdg
442         << "p="       << p
443         << "ymc="     << y0
444         << "zmc="     << z0
445         << "dydx="    << dydx
446         << "dzdx="    << dzdx
447         << "cross="   << cross
448         << "dy="                  << dy
449         << "dz="                  << dz
450         << "dphi="              << dphi
451         << "\n";
452     }
453
454     Int_t istk = AliTRDgeometry::GetStack(det); 
455     AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
456     Float_t zr0 = pp->GetRow0() + AliTRDSimParam::Instance()->GetAnodeWireOffset();
457     Float_t tilt = fTracklet->GetTilt();
458
459     AliTRDcluster *c = 0x0;
460     fTracklet->ResetClusterIter(kFALSE);
461     while((c = fTracklet->PrevCluster())){
462       Float_t  q = TMath::Abs(c->GetQ());
463       Float_t xc = c->GetX();
464       Float_t yc = c->GetY();
465       Float_t zc = c->GetZ();
466       dx = x0 - xc; 
467       Float_t yt = y0 - dx*dydx;
468       Float_t zt = z0 - dx*dzdx; 
469       dy = yt - (yc - tilt*(zc-zt));
470
471       // Fill Histograms
472       if(q>100.) ((TH2I*)fContainer->At(kClusterResolution))->Fill(dydx, dy);
473       
474       // Fill Debug Tree
475       if(fDebugLevel>=1){
476         Float_t d = zr0 - zt;
477         d -= ((Int_t)(2 * d)) / 2.0;
478         if (d > 0.25) d  = 0.5 - d;
479
480         AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
481         fClResolution->Add(clInfo);
482         clInfo->SetCluster(c);
483         clInfo->SetMC(pdg, label);
484         clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
485         clInfo->SetResolution(dy);
486         clInfo->SetAnisochronity(d);
487         clInfo->SetDriftLength(x0-xc);
488         //clInfo->Print();
489         (*fDebugStream) << "ClusterResolution"
490           <<"clInfo.=" << clInfo
491           << "\n";
492       }
493     }
494   }
495   return h;
496 }
497
498
499 //________________________________________________________
500 void AliTRDtrackingResolution::GetRefFigure(Int_t ifig)
501 {
502   TBox *b = 0x0;
503   TAxis *ax = 0x0;
504   TGraphErrors *g = 0x0;
505   switch(ifig){
506   case kClusterResidual:
507     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
508     g->Draw("apl");
509     ax = g->GetHistogram()->GetYaxis();
510     ax->SetRangeUser(-.5, 2.5);
511     ax->SetTitle("Clusters Y Residuals #sigma/#mu [mm]");
512     ax = g->GetHistogram()->GetXaxis();
513     ax->SetTitle("tg(#phi)");
514     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
515     g->Draw("pl");
516     b = new TBox(-.15, -.5, .15, 1.);
517     b->SetFillStyle(3002);b->SetFillColor(kGreen);
518     b->SetLineColor(0); b->Draw();
519     return;
520   case kTrackletYResidual:
521     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
522     g->Draw("apl");
523     ax = g->GetHistogram()->GetYaxis();
524     ax->SetRangeUser(-.5, 3.);
525     ax->SetTitle("Tracklet Y Residuals #sigma/#mu [mm]");
526     ax = g->GetHistogram()->GetXaxis();
527     ax->SetTitle("tg(#phi)");
528     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
529     g->Draw("pl");
530     b = new TBox(-.15, -.5, .15, 1.);
531     b->SetFillStyle(3002);b->SetFillColor(kGreen);
532     b->SetLineColor(0); b->Draw();
533     return;
534   case kTrackletPhiResidual:
535     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
536     g->Draw("apl");
537     ax = g->GetHistogram()->GetYaxis();
538     ax->SetRangeUser(-.5, 2.);
539     ax->SetTitle("Tracklet Phi Residuals #sigma/#mu [rad]");
540     ax = g->GetHistogram()->GetXaxis();
541     ax->SetTitle("tg(#phi)");
542     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
543     g->Draw("pl");
544     b = new TBox(-.15, -.5, .15, 1.);
545     b->SetFillStyle(3002);b->SetFillColor(kGreen);
546     b->SetLineColor(0); b->Draw();
547     return;
548   case kClusterResolution:
549     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
550     ax = g->GetHistogram()->GetYaxis();
551     ax->SetRangeUser(-.5, 1.);
552     ax->SetTitle("Cluster Y Resolution #sigma/#mu [mm]");
553     ax = g->GetHistogram()->GetXaxis();
554     ax->SetTitle("tg(#phi)");
555     g->Draw("apl");
556     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
557     g->Draw("pl");
558     b = new TBox(-.15, -.5, .15, 1.);
559     b->SetFillStyle(3002);b->SetFillColor(kGreen);
560     b->SetLineColor(0); b->Draw();
561     return;
562   case kTrackletYResolution:
563     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
564     ax = g->GetHistogram()->GetYaxis();
565     ax->SetRangeUser(-.5, 1.);
566     ax->SetTitle("Tracklet Y Resolution #sigma/#mu [mm]");
567     ax = g->GetHistogram()->GetXaxis();
568     ax->SetTitle("tg(#phi)");
569     g->Draw("apl");
570     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
571     g->Draw("pl");
572     return;
573   case kTrackletZResolution:
574     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
575     ax = g->GetHistogram()->GetYaxis();
576     ax->SetRangeUser(-.5, 1.);
577     ax->SetTitle("Tracklet Z Resolution #sigma/#mu [mm]");
578     ax = g->GetHistogram()->GetXaxis();
579     ax->SetTitle("tg(#theta)");
580     g->Draw("apl");
581     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
582     g->Draw("pl");
583     return;
584   case kTrackletAngleResolution:
585     if(!(g = (TGraphErrors*)fGraphS->At(ifig))) break;
586     ax = g->GetHistogram()->GetYaxis();
587     ax->SetRangeUser(-.05, .2);
588     ax->SetTitle("Tracklet Angular Resolution #sigma/#mu [deg]");
589     ax = g->GetHistogram()->GetXaxis();
590     ax->SetTitle("tg(#phi)");
591     g->Draw("apl");
592     if(!(g = (TGraphErrors*)fGraphM->At(ifig))) break;
593     g->Draw("pl");
594     return;
595   default:
596     AliInfo(Form("Reference plot [%d] not implemented yet", ifig));
597     return;
598   }
599   AliInfo(Form("Reference plot [%d] missing result", ifig));
600 }
601
602
603 //________________________________________________________
604 Bool_t AliTRDtrackingResolution::PostProcess()
605 {
606   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
607   if (!fContainer) {
608     Printf("ERROR: list not available");
609     return kFALSE;
610   }
611   fNRefFigures = fContainer->GetEntriesFast();
612   if(!fGraphS){ 
613     fGraphS = new TObjArray(fNRefFigures);
614     fGraphS->SetOwner();
615   }
616   if(!fGraphM){ 
617     fGraphM = new TObjArray(fNRefFigures);
618     fGraphM->SetOwner();
619   }
620
621   TH2I *h2 = 0x0;
622   TH1D *h = 0x0;
623   TGraphErrors *gm = 0x0, *gs = 0x0;
624
625   // define models
626   TF1 f("f1", "gaus", -.5, .5);  
627
628   TF1 fb("fb", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]", -.5, .5);
629
630   TF1 fc("fc", "[0]*exp(-0.5*((x-[1])/[2])**2)+[3]*exp(-0.5*((x-[4])/[5])**2)", -.5, .5);
631
632   TCanvas *c = 0x0;
633   if(IsVisual()) c = new TCanvas("c", Form("%s Visual", GetName()), 500, 500);
634   char opt[5];
635   sprintf(opt, "%c%c", IsVerbose() ? ' ' : 'Q', IsVisual() ? ' ': 'N');
636
637
638   //PROCESS RESIDUAL DISTRIBUTIONS
639
640   // Clusters residuals
641   h2 = (TH2I *)(fContainer->At(kClusterResidual));
642   gm = new TGraphErrors();
643   gm->SetLineColor(kBlue);
644   gm->SetMarkerStyle(7);
645   gm->SetMarkerColor(kBlue);
646   gm->SetNameTitle("clm", "");
647   fGraphM->AddAt(gm, kClusterResidual);
648   gs = new TGraphErrors();
649   gs->SetLineColor(kRed);
650   gs->SetMarkerStyle(23);
651   gs->SetMarkerColor(kRed);
652   gs->SetNameTitle("cls", "");
653   fGraphS->AddAt(gs, kClusterResidual);
654   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
655     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
656     h = h2->ProjectionY("py", ibin, ibin);
657     if(h->GetEntries()<100) continue;
658     AdjustF1(h, &f);
659
660     if(IsVisual()){c->cd(); c->SetLogy();}
661     h->Fit(&f, opt, "", -0.5, 0.5);
662     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
663     
664     Int_t ip = gm->GetN();
665     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
666     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
667     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
668     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
669   }
670
671   // Tracklet y residuals
672   h2 = (TH2I *)(fContainer->At(kTrackletYResidual));
673   gm = new TGraphErrors();
674   gm->SetLineColor(kBlue);
675   gm->SetMarkerStyle(7);
676   gm->SetMarkerColor(kBlue);
677   gm->SetNameTitle("tktm", "");
678   fGraphM->AddAt(gm, kTrackletYResidual);
679   gs = new TGraphErrors();
680   gs->SetLineColor(kRed);
681   gs->SetMarkerStyle(23);
682   gs->SetMarkerColor(kRed);
683   gs->SetNameTitle("tkts", "");
684   fGraphS->AddAt(gs, kTrackletYResidual);
685   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
686     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
687     h = h2->ProjectionY("py", ibin, ibin);
688     if(h->GetEntries()<100) continue;
689     AdjustF1(h, &f);
690
691     if(IsVisual()){c->cd(); c->SetLogy();}
692     h->Fit(&f, opt, "", -0.5, 0.5);
693     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
694     
695     Int_t ip = gm->GetN();
696     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
697     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
698     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
699     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
700   }
701
702   // Tracklet phi residuals
703   h2 = (TH2I *)(fContainer->At(kTrackletPhiResidual));
704   gm = new TGraphErrors();
705   gm->SetLineColor(kBlue);
706   gm->SetMarkerStyle(7);
707   gm->SetMarkerColor(kBlue);
708   gm->SetNameTitle("tktphim", "");
709   fGraphM->AddAt(gm, kTrackletPhiResidual);
710   gs = new TGraphErrors();
711   gs->SetLineColor(kRed);
712   gs->SetMarkerStyle(23);
713   gs->SetMarkerColor(kRed);
714   gs->SetNameTitle("tktphis", "");
715   fGraphS->AddAt(gs, kTrackletPhiResidual);
716   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
717     Double_t phi = h2->GetXaxis()->GetBinCenter(ibin);
718     h = h2->ProjectionY("py", ibin, ibin);
719     if(h->GetEntries()<100) continue;
720     AdjustF1(h, &f);
721
722     if(IsVisual()){c->cd(); c->SetLogy();}
723     h->Fit(&f, opt, "", -0.5, 0.5);
724     if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
725     
726     Int_t ip = gm->GetN();
727     gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
728     gm->SetPointError(ip, 0., 10.*f.GetParError(1));
729     gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
730     gs->SetPointError(ip, 0., 10.*f.GetParError(2));
731   }
732
733
734   //PROCESS RESOLUTION DISTRIBUTIONS
735
736   if(HasMCdata()){
737     // cluster y resolution
738     h2 = (TH2I*)fContainer->At(kClusterResolution);
739     gm = new TGraphErrors();
740     gm->SetLineColor(kBlue);
741     gm->SetMarkerStyle(7);
742     gm->SetMarkerColor(kBlue);
743     gm->SetNameTitle("clym", "");
744     fGraphM->AddAt(gm, kClusterResolution);
745     gs = new TGraphErrors();
746     gs->SetLineColor(kRed);
747     gs->SetMarkerStyle(23);
748     gs->SetMarkerColor(kRed);
749     gs->SetNameTitle("clys", "");
750     fGraphS->AddAt(gs, kClusterResolution);
751     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
752       h = h2->ProjectionY("py", iphi, iphi);
753       if(h->GetEntries()<100) continue;
754       AdjustF1(h, &f);
755
756       if(IsVisual()){c->cd(); c->SetLogy();}
757       h->Fit(&f, opt, "", -0.5, 0.5);
758       if(IsVerbose()){
759         printf("phi[%d] mean[%e] sigma[%e]\n\n", iphi, 10.*f.GetParameter(1), 10.*f.GetParameter(2));
760       }
761       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
762
763       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
764       Int_t ip = gm->GetN();
765       gm->SetPoint(ip, phi, 10.*f.GetParameter(1));
766       gm->SetPointError(ip, 0., 10.*f.GetParError(1));
767       gs->SetPoint(ip, phi, 10.*f.GetParameter(2));
768       gs->SetPointError(ip, 0., 10.*f.GetParError(2));
769     }
770   
771     // tracklet y resolution
772     h2 = (TH2I*)fContainer->At(kTrackletYResolution);
773     gm = new TGraphErrors();
774     gm->SetLineColor(kBlue);
775     gm->SetMarkerStyle(7);
776     gm->SetMarkerColor(kBlue);
777     gm->SetNameTitle("trkltym", "");
778     fGraphM->AddAt(gm, kTrackletYResolution);
779     gs = new TGraphErrors();
780     gs->SetLineColor(kRed);
781     gs->SetMarkerStyle(23);
782     gs->SetMarkerColor(kRed);
783     gs->SetNameTitle("trkltys", "");
784     fGraphS->AddAt(gs, kTrackletYResolution);
785     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
786       h = h2->ProjectionY("py", iphi, iphi);
787       if(h->GetEntries()<100) continue;
788       AdjustF1(h, &fb);
789
790       if(IsVisual()){c->cd(); c->SetLogy();}
791       h->Fit(&fb, opt, "", -0.5, 0.5);
792       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
793
794       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
795       Int_t ip = gm->GetN();
796       gm->SetPoint(ip, phi, 10.*fb.GetParameter(1));
797       gm->SetPointError(ip, 0., 10.*fb.GetParError(1));
798       gs->SetPoint(ip, phi, 10.*fb.GetParameter(2));
799       gs->SetPointError(ip, 0., 10.*fb.GetParError(2));
800     }
801   
802     // tracklet z resolution
803     h2 = (TH2I*)fContainer->At(kTrackletZResolution);
804     gm = new TGraphErrors();
805     gm->SetLineColor(kBlue);
806     gm->SetMarkerStyle(7);
807     gm->SetMarkerColor(kBlue);
808     gm->SetNameTitle("trkltzm", "");
809     fGraphM->AddAt(gm, kTrackletZResolution);
810     gs = new TGraphErrors();
811     gs->SetLineColor(kRed);
812     gs->SetMarkerStyle(23);
813     gs->SetMarkerColor(kRed);
814     gs->SetNameTitle("trkltzs", "");
815     fGraphS->AddAt(gs, kTrackletZResolution);
816     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
817       h = h2->ProjectionY("py", iphi, iphi);
818       if(h->GetEntries()<100) continue;
819       AdjustF1(h, &fb);
820
821       if(IsVisual()){c->cd(); c->SetLogy();}
822       h->Fit(&fb, opt, "", -0.5, 0.5);
823       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
824
825       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
826       Int_t ip = gm->GetN();
827       gm->SetPoint(ip, phi, 10.*fb.GetParameter(1));
828       gm->SetPointError(ip, 0., 10.*fb.GetParError(1));
829       gs->SetPoint(ip, phi, 10.*fb.GetParameter(2));
830       gs->SetPointError(ip, 0., 10.*fb.GetParError(2));
831     }
832   
833     // tracklet phi resolution
834     h2 = (TH2I*)fContainer->At(kTrackletAngleResolution);
835     gm = new TGraphErrors();
836     gm->SetLineColor(kBlue);
837     gm->SetMarkerStyle(7);
838     gm->SetMarkerColor(kBlue);
839     gm->SetNameTitle("trkltym", "");
840     fGraphM->AddAt(gm, kTrackletAngleResolution);
841     gs = new TGraphErrors();
842     gs->SetLineColor(kRed);
843     gs->SetMarkerStyle(23);
844     gs->SetMarkerColor(kRed);
845     gs->SetNameTitle("trkltys", "");
846     fGraphS->AddAt(gs, kTrackletAngleResolution);
847     for(Int_t iphi=1; iphi<=h2->GetNbinsX(); iphi++){
848       h = h2->ProjectionY("py", iphi, iphi);
849       if(h->GetEntries()<100) continue;
850
851       if(IsVisual()){c->cd(); c->SetLogy();}
852       h->Fit(&f, opt, "", -0.5, 0.5);
853       if(IsVisual()){c->Modified(); c->Update(); gSystem->Sleep(500);}
854
855       Double_t phi = h2->GetXaxis()->GetBinCenter(iphi);
856       Int_t ip = gm->GetN();
857       gm->SetPoint(ip, phi, f.GetParameter(1));
858       gm->SetPointError(ip, 0., f.GetParError(1));
859       gs->SetPoint(ip, phi, f.GetParameter(2));
860       gs->SetPointError(ip, 0., f.GetParError(2));
861     }
862   }
863   if(c) delete c;
864
865   return kTRUE;
866 }
867
868
869 //________________________________________________________
870 void AliTRDtrackingResolution::Terminate(Option_t *)
871 {
872   if(fDebugStream){ 
873     delete fDebugStream;
874     fDebugStream = 0x0;
875     fDebugLevel = 0;
876   }
877   if(HasPostProcess()) PostProcess();
878 }
879
880 //________________________________________________________
881 void AliTRDtrackingResolution::AdjustF1(TH1 *h, TF1 *f)
882 {
883 // Helper function to avoid duplication of code
884 // Make first guesses on the fit parameters
885
886   // find the intial parameters of the fit !! (thanks George)
887   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
888   Double_t sum = 0.;
889   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
890   f->SetParLimits(0, 0., 3.*sum);
891   f->SetParameter(0, .9*sum);
892
893   f->SetParLimits(1, -.2, .2);
894   f->SetParameter(1, -0.1);
895
896   f->SetParLimits(2, 0., 4.e-1);
897   f->SetParameter(2, 2.e-2);
898   if(f->GetNpar() <= 4) return;
899
900   f->SetParLimits(3, 0., sum);
901   f->SetParameter(3, .1*sum);
902
903   f->SetParLimits(4, -.3, .3);
904   f->SetParameter(4, 0.);
905
906   f->SetParLimits(5, 0., 1.e2);
907   f->SetParameter(5, 2.e-1);
908 }
909
910 //________________________________________________________
911 TObjArray* AliTRDtrackingResolution::Histos()
912 {
913   if(fContainer) return fContainer;
914
915   fContainer  = new TObjArray(7);
916
917   TH1 *h = 0x0;
918   // cluster to tracklet residuals [2]
919   if(!(h = (TH2I*)gROOT->FindObject("fYCl"))){
920     h = new TH2I("fYCl", "Clusters Residuals", 21, -.33, .33, 100, -.5, .5);
921     h->GetXaxis()->SetTitle("tg(#phi)");
922     h->GetYaxis()->SetTitle("#Delta y [cm]");
923     h->GetZaxis()->SetTitle("entries");
924   } else h->Reset();
925   fContainer->AddAt(h, kClusterResidual);
926
927   // tracklet to track residuals [2]
928   if(!(h = (TH2I*)gROOT->FindObject("hTrkltYRez"))){ 
929     h = new TH2I("hTrkltYRez", "Tracklets", 21, -.33, .33, 100, -.5, .5);
930     h->GetXaxis()->SetTitle("tg(#phi)");
931     h->GetYaxis()->SetTitle("#Delta y [cm]");
932     h->GetZaxis()->SetTitle("entries");
933   } else h->Reset();
934   fContainer->AddAt(h, kTrackletYResidual);
935
936   // tracklet to track residuals angular [2]
937   if(!(h = (TH2I*)gROOT->FindObject("hTrkltPhiRez"))){ 
938     h = new TH2I("hTrkltPhiRez", "Tracklets", 21, -.33, .33, 100, -.5, .5);
939     h->GetXaxis()->SetTitle("tg(#phi)");
940     h->GetYaxis()->SetTitle("#Delta phi [#circ]");
941     h->GetZaxis()->SetTitle("entries");
942   } else h->Reset();
943   fContainer->AddAt(h, kTrackletPhiResidual);
944
945
946   // Resolution histos
947   if(HasMCdata()){
948     // cluster y resolution [0]
949     if(!(h = (TH2I*)gROOT->FindObject("fYClMC"))){
950       h = new TH2I("fYClMC", "Cluster Resolution", 31, -.48, .48, 100, -.5, .5);
951       h->GetXaxis()->SetTitle("tg(#phi)");
952       h->GetYaxis()->SetTitle("#Delta y [cm]");
953       h->GetZaxis()->SetTitle("entries");
954     } else h->Reset();
955     fContainer->AddAt(h, kClusterResolution);
956
957     // tracklet y resolution [0]
958     if(!(h = (TH2I*)gROOT->FindObject("fYTrkltMC"))){
959       h = new TH2I("fYTrkltMC", "Tracklet Resolution (Y)", 31, -.48, .48, 100, -.5, .5);
960       h->GetXaxis()->SetTitle("tg(#phi)");
961       h->GetYaxis()->SetTitle("#Delta y [cm]");
962       h->GetZaxis()->SetTitle("entries");
963     } else h->Reset();
964     fContainer->AddAt(h, kTrackletYResolution);
965
966     // tracklet y resolution [0]
967     if(!(h = (TH2I*)gROOT->FindObject("fZTrkltMC"))){
968       h = new TH2I("fZTrkltMC", "Tracklet Resolution (Z)", 31, -.48, .48, 100, -.5, .5);
969       h->GetXaxis()->SetTitle("tg(#theta)");
970       h->GetYaxis()->SetTitle("#Delta z [cm]");
971       h->GetZaxis()->SetTitle("entries");
972     } else h->Reset();
973     fContainer->AddAt(h, kTrackletZResolution);
974
975     // tracklet angular resolution [1]
976     if(!(h = (TH2I*)gROOT->FindObject("fPhiTrkltMC"))){
977       h = new TH2I("fPhiTrkltMC", "Tracklet Resolution (Angular)", 31, -.48, .48, 100, -10., 10.);
978       h->GetXaxis()->SetTitle("tg(#phi)");
979       h->GetYaxis()->SetTitle("#Delta #phi [deg]");
980       h->GetZaxis()->SetTitle("entries");
981     } else h->Reset();
982     fContainer->AddAt(h, kTrackletAngleResolution);
983   }
984   return fContainer;
985 }
986
987
988 //________________________________________________________
989 void AliTRDtrackingResolution::SetRecoParam(AliTRDrecoParam *r)
990 {
991
992   fReconstructor->SetRecoParam(r);
993 }