- call UserExec instead of Exec of base class AliTRDreco
[u/mrichter/AliRoot.git] / PWG1 / TRD / 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();
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 <TSystem.h>
52
53 #include <TROOT.h>
54 #include <TObjArray.h>
55 #include <TH3.h>
56 #include <TH2.h>
57 #include <TH1.h>
58 #include <TF1.h>
59 #include <TCanvas.h>
60 #include <TGaxis.h>
61 #include <TBox.h>
62 #include <TLegend.h>
63 #include <TGraphErrors.h>
64 #include <TGraphAsymmErrors.h>
65 #include <TMath.h>
66 #include <TMatrixT.h>
67 #include <TVectorT.h>
68 #include <TTreeStream.h>
69 #include <TGeoManager.h>
70 #include <TDatabasePDG.h>
71
72 #include "AliPID.h"
73 #include "AliESDtrack.h"
74
75 #include "AliTRDresolution.h"
76 #include "AliTRDgeometry.h"
77 #include "AliTRDpadPlane.h"
78 #include "AliTRDcluster.h"
79 #include "AliTRDseedV1.h"
80 #include "AliTRDtrackV1.h"
81 #include "AliTRDReconstructor.h"
82 #include "AliTRDrecoParam.h"
83 #include "AliTRDpidUtil.h"
84
85 #include "info/AliTRDclusterInfo.h"
86
87 ClassImp(AliTRDresolution)
88
89 UChar_t const AliTRDresolution::fgNhistos[kNviews] = {
90   2, 2, 5, 5,
91   2, 5, 12, 2, 14
92 };
93 UChar_t const AliTRDresolution::fgNproj[kNviews] = {
94   2, 2, 5, 5,
95   4, 7, 12, 2, 14
96 };
97 Char_t const * AliTRDresolution::fgPerformanceName[kNviews] = {
98      "Charge"
99     ,"Cluster2Track"
100     ,"Tracklet2Track"
101     ,"Tracklet2TPC" 
102     ,"Cluster2MC"
103     ,"Tracklet2MC"
104     ,"TPC2MC"
105     ,"TOF/HMPID2MC"
106     ,"TRD2MC"
107 };
108 UChar_t const AliTRDresolution::fgNcomp[kNprojs] = {
109   1, 1, //2, 
110   1, 1, //2, 
111   1, 1, 1, 1, 1, //5, 
112   1, 1, 1, 1, 1, 1, 1, //5,
113 // MC
114   1, 1, 1, 1,    //4, 
115   1, 1, 1, 1, 1, //5, 
116   1, 1, 1, 1, 1, 1, 1, 1, 11, 11, 11, 11, //12, 
117   1, 1, //2, 
118   1, 1, 1, 1, 1, 1, 1, 1, 66, 66, 66, 66, 66, 66 //14
119 };
120 Char_t const *AliTRDresolution::fgAxTitle[kNprojs][4] = {
121   // Charge
122   {"Impv", "x [cm]", "I_{mpv}", "x/x_{0}"}
123  ,{"dI/Impv", "x/x_{0}", "#delta I/I_{mpv}", "x[cm]"}
124   // Clusters to Kalman
125  ,{"Cluster2Track residuals", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
126  ,{"Cluster2Track  pulls", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
127   // TRD tracklet to Kalman fit
128  ,{"Tracklet2Track Y residuals", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
129  ,{"Tracklet2Track Y pulls", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
130  ,{"Tracklet2Track Z residuals", "tg(#theta)", "#mu_{z} [#mum]", "#sigma_{z} [#mum]"}
131  ,{"Tracklet2Track Z pulls", "tg(#theta)", "#mu_{z}", "#sigma_{z}"}
132  ,{"Tracklet2Track Phi residuals", "tg(#phi)", "#mu_{#phi} [mrad]", "#sigma_{#phi} [mrad]"}
133   // TPC track 2 first TRD tracklet
134  ,{"Tracklet2Track Y residuals @ TRDin", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
135  ,{"Tracklet2Track Y pulls @ TRDin", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
136  ,{"Tracklet2Track Z residuals @ TRDin", "tg(#theta)", "#mu_{z} [#mum]", "#sigma_{z} [#mum]"}
137  ,{"Tracklet2Track Z pulls @ TRDin", "tg(#theta)", "#mu_{z}", "#sigma_{z}"}
138  ,{"Tracklet2Track Phi residuals @ TRDin", "tg(#phi)", "#mu_{#phi} [mrad]", "#sigma_{#phi} [mrad]"}
139   // MC cluster
140  ,{"MC Cluster Y resolution (p_{t}<1 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
141  ,{"MC Cluster Y resolution (1<p_{t}<2 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
142  ,{"MC Cluster Y resolution (p_{t}>3 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
143  ,{"MC Cluster Y pulls", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
144   // MC tracklet
145  ,{"MC Tracklet Y resolution (p_{t}<1 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
146  ,{"MC Tracklet Y resolution (1<p_{t}<2 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
147  ,{"MC Tracklet Y resolution (p_{t}>3 GeV/c)", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y}[#mum]"}
148  ,{"MC Tracklet Y pulls", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
149  ,{"MC Tracklet Cross Z resolution", "tg(#theta)", "#mu_{z} [#mum]", "#sigma_{z} [#mum]"}
150  ,{"MC Tracklet Cross Z pulls", "tg(#theta)", "#mu_{z}", "#sigma_{z}"}
151  ,{"MC Tracklet Phi resolution", "tg(#phi)", "#mu_{#phi} [mrad]", "#sigma_{#phi} [mrad]"}
152   // MC track TPC
153  ,{"Y resolution @ TRDin", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y}[#mum]"}
154  ,{"Y pulls @ TRDin", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
155  ,{"Z resolution @ TRDin", "tg(#theta)", "#mu_{z} [#mum]", "#sigma_{z} [#mum]"}
156  ,{"Z pulls @ TRDin", "tg(#theta)", "#mu_{z}", "#sigma_{z}"}
157  ,{"Phi resolution @ TRDin", "tg(#phi)", "#mu_{#phi} [mrad]", "#sigma_{#phi} [mrad]"}
158  ,{"SNP pulls @ TRDin", "tg(#phi)", "#mu_{snp}", "#sigma_{snp}"}
159  ,{"Theta resolution @ TRDin", "tg(#theta)", "#mu_{#theta} [mrad]", "#sigma_{#theta} [mrad]"}
160  ,{"TGL pulls @ TRDin", "tg(#theta)", "#mu_{tgl}", "#sigma_{tgl}"}
161  ,{"P_{t} resolution @ TRDin", "p_{t}^{MC} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "MC: #sigma^{TPC}(#Deltap_{t}/p_{t}^{MC}) [%]"}
162  ,{"1/P_{t} pulls @ TRDin", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC}-1/p_{t}^{MC}", "MC PULL: #sigma_{1/p_{t}}^{TPC}"}
163  ,{"P resolution @ TRDin", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
164  ,{"P pulls @ TRDin", "p^{MC} [GeV/c]", "1/p^{REC}-1/p^{MC}", "MC PULL: #sigma^{TPC}(#Deltap/#sigma_{p})"}
165   // MC track TOF
166  ,{"PosZ", "tg(#theta)", "MC: #mu_{z}^{TOF} [#mum]", "MC: #sigma_{z}^{TOF} [#mum]"}
167  ,{"PullsZ", "tg(#theta)", "MC PULL: #mu_{z}^{TOF}", "MC PULL: #sigma_{z}^{TOF}"}
168   // MC track in TRD
169  ,{"TRD track MC Y resolution", "tg(#phi)", "#mu_{y}^{Trk} [#mum]", "#sigma_{y}^{Trk} [#mum]"}
170  ,{"TRD track MC Y pulls", "tg(#phi)", "#mu_{y}^{Trk}", "#sigma_{y}^{Trk}"}
171  ,{"TRD track MC Z resolution", "tg(#theta)", "#mu_{z}^{Trk} [#mum]", "#sigma_{z}^{Trk} [#mum]"}
172  ,{"TRD track MC Z pulls", "tg(#theta)", "#mu_{z}^{Trk}", "#sigma_{z}^{Trk}"}
173  ,{"TRD track MC Phi resolution", "tg(#phi)", "#mu_{#phi}^{Trk} [mrad]", "#sigma_{#phi}^{Trk} [mrad]"}
174  ,{"TRD track MC SNP pulls", "tg(#phi)", "#mu_{snp}^{Trk}", "#sigma_{snp}^{Trk}"}
175  ,{"TRD track MC Theta resolution", "tg(#theta)", "#mu_{#theta}^{Trk} [mrad]", "#sigma_{#theta}^{Trk} [mrad]"}
176  ,{"TRD track MC TGL pulls", "tg(#theta)", "#mu_{tgl}^{Trk}", "#sigma_{tgl}^{Trk}"}
177  ,{"P_{t} resolution TRD Layer", "p_{t} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "#sigma(#Deltap_{t}/p_{t}^{MC}) [%]"}
178  ,{"1/P_{t} pulls TRD Layer", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC} - 1/p_{t}^{MC}", "#sigma_{1/p_{t}}"}
179  ,{"P resolution TRD Layer", "p [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "#sigma(#Deltap/p^{MC}) [%]"}
180  ,{"[SA] P_{t} resolution TRD Layer", "p_{t}^{MC} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "MC: #sigma^{Trk}(#Deltap_{t}/p_{t}^{MC}) [%]"}
181  ,{"[SA] 1/P_{t} pulls TRD Layer", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC}-1/p_{t}^{MC}", "MC PULL: #sigma_{1/p_{t}}^{Trk}"}
182  ,{"[SA] P resolution TRD Layer", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{Trk}(#Deltap/p^{MC}) [%]"}
183 };
184
185 //________________________________________________________
186 AliTRDresolution::AliTRDresolution()
187   :AliTRDrecoTask()
188   ,fStatus(0)
189   ,fIdxPlot(0)
190   ,fIdxFrame(0)
191   ,fReconstructor(NULL)
192   ,fGeo(NULL)
193   ,fDBPDG(NULL)
194   ,fGraphS(NULL)
195   ,fGraphM(NULL)
196   ,fCl(NULL)
197   ,fTrklt(NULL)
198   ,fMCcl(NULL)
199   ,fMCtrklt(NULL)
200 {
201   //
202   // Default constructor
203   //
204   SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
205 }
206
207 AliTRDresolution::AliTRDresolution(char* name)
208   :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
209   ,fStatus(0)
210   ,fIdxPlot(0)
211   ,fIdxFrame(0)
212   ,fReconstructor(NULL)
213   ,fGeo(NULL)
214   ,fDBPDG(NULL)
215   ,fGraphS(NULL)
216   ,fGraphM(NULL)
217   ,fCl(NULL)
218   ,fTrklt(NULL)
219   ,fMCcl(NULL)
220   ,fMCtrklt(NULL)
221 {
222   //
223   // Default constructor
224   //
225
226   fReconstructor = new AliTRDReconstructor();
227   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
228   fGeo = new AliTRDgeometry();
229
230   InitFunctorList();
231
232   DefineOutput(1, TObjArray::Class()); // cluster2track
233   DefineOutput(2, TObjArray::Class()); // tracklet2track
234   DefineOutput(3, TObjArray::Class()); // cluster2mc
235   DefineOutput(4, TObjArray::Class()); // tracklet2mc
236 }
237
238 //________________________________________________________
239 AliTRDresolution::~AliTRDresolution()
240 {
241   //
242   // Destructor
243   //
244
245   if(fGraphS){fGraphS->Delete(); delete fGraphS;}
246   if(fGraphM){fGraphM->Delete(); delete fGraphM;}
247   delete fGeo;
248   delete fReconstructor;
249   if(gGeoManager) delete gGeoManager;
250   if(fCl){fCl->Delete(); delete fCl;}
251   if(fTrklt){fTrklt->Delete(); delete fTrklt;}
252   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
253   if(fMCtrklt){fMCtrklt->Delete(); delete fMCtrklt;}
254 }
255
256
257 //________________________________________________________
258 void AliTRDresolution::UserCreateOutputObjects()
259 {
260   // spatial resolution
261   OpenFile(1, "RECREATE");
262
263   fContainer = Histos();
264
265   fCl = new TObjArray();
266   fCl->SetOwner(kTRUE);
267   fTrklt = new TObjArray();
268   fTrklt->SetOwner(kTRUE);
269   fMCcl = new TObjArray();
270   fMCcl->SetOwner(kTRUE);
271   fMCtrklt = new TObjArray();
272   fMCtrklt->SetOwner(kTRUE);
273 }
274
275 //________________________________________________________
276 void AliTRDresolution::UserExec(Option_t *opt)
277 {
278   //
279   // Execution part
280   //
281
282   fCl->Delete();
283   fTrklt->Delete();
284   fMCcl->Delete();
285   fMCtrklt->Delete();
286
287   AliTRDrecoTask::UserExec(opt);
288
289   PostData(1, fCl);
290   PostData(2, fTrklt);
291   PostData(3, fMCcl);
292   PostData(4, fMCtrklt);
293 }
294
295 //________________________________________________________
296 TH1* AliTRDresolution::PlotCharge(const AliTRDtrackV1 *track)
297 {
298   //
299   // Plots the charge distribution
300   //
301
302   if(track) fkTrack = track;
303   if(!fkTrack){
304     AliDebug(2, "No Track defined.");
305     return NULL;
306   }
307   TObjArray *arr = NULL;
308   if(!(arr = ((TObjArray*)fContainer->At(kCharge)))){
309     AliWarning("No output container defined.");
310     return NULL;
311   }
312   TH3S* h = NULL;
313
314   AliTRDseedV1 *fTracklet = NULL;  
315   AliTRDcluster *c = NULL;
316   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
317     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
318     if(!fTracklet->IsOK()) continue;
319     Float_t x0 = fTracklet->GetX0();
320     Float_t dq, dl;
321     for(Int_t itb=AliTRDseedV1::kNtb; itb--;){
322       if(!(c = fTracklet->GetClusters(itb))){ 
323         if(!(c = fTracklet->GetClusters(AliTRDseedV1::kNtb+itb))) continue;
324       }
325       dq = fTracklet->GetdQdl(itb, &dl);
326       dl /= 0.15; // dl/dl0, dl0 = 1.5 mm for nominal vd
327       (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dq);
328     }
329
330 //     if(!HasMCdata()) continue;
331 //     UChar_t s;
332 //     Float_t pt0, y0, z0, dydx0, dzdx0;
333 //     if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
334
335   }
336   return h;
337 }
338
339
340 //________________________________________________________
341 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
342 {
343   //
344   // Plot the cluster distributions
345   //
346
347   if(track) fkTrack = track;
348   if(!fkTrack){
349     AliDebug(2, "No Track defined.");
350     return NULL;
351   }
352   TObjArray *arr = NULL;
353   if(!(arr = ((TObjArray*)fContainer->At(kCluster)))){
354     AliWarning("No output container defined.");
355     return NULL;
356   }
357   ULong_t status = fkESD ? fkESD->GetStatus():0;
358
359   Double_t covR[7], cov[3];
360   Float_t x0, y0, z0, dy, dz, dydx, dzdx;
361   AliTRDseedV1 *fTracklet(NULL);  
362   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
363     if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
364     if(!fTracklet->IsOK()) continue;
365     x0 = fTracklet->GetX0();
366
367     // retrive the track angle with the chamber
368     y0   = fTracklet->GetYref(0);
369     z0   = fTracklet->GetZref(0);
370     dydx = fTracklet->GetYref(1);
371     dzdx = fTracklet->GetZref(1);
372     fTracklet->GetCovRef(covR);
373     Float_t tilt(fTracklet->GetTilt());
374     AliTRDcluster *c = NULL;
375     fTracklet->ResetClusterIter(kFALSE);
376     while((c = fTracklet->PrevCluster())){
377       Float_t xc = c->GetX();
378       Float_t yc = c->GetY();
379       Float_t zc = c->GetZ();
380       Float_t dx = x0 - xc; 
381       Float_t yt = y0 - dx*dydx;
382       Float_t zt = z0 - dx*dzdx; 
383       // calculate residuals using tilt correction
384 //       yc -= tilt*(zc-zt); 
385 //       dy = yt - yc;
386
387       // calculate residuals using correct covariance matrix
388       cov[0] = c->GetSigmaY2();
389       cov[1] = c->GetSigmaYZ();
390       cov[2] = c->GetSigmaZ2();
391       // do rotation
392       Double_t sy2(cov[0]), sz2(cov[2]);
393       Double_t t2 = tilt*tilt;
394       Double_t correction = 1./(1. + t2);
395       cov[0] = (sy2+t2*sz2)*correction;
396       cov[1] = tilt*(sz2 - sy2)*correction;
397       cov[2] = (t2*sy2+sz2)*correction;
398       // do inversion
399       Double_t r00=cov[0]+covR[0], r01=cov[1]+covR[1], r11=cov[2]+covR[2];
400       Double_t det=r00*r11 - r01*r01;
401       Double_t tmp=r00; r00=r11/det; r11=tmp/det;
402       dy = (yc - yt)*TMath::Sqrt(r00);
403       dz = (zc - zt)*TMath::Sqrt(r11);
404
405       ((TH2I*)arr->At(0))->Fill(dydx, dy/*, dz*/);
406       ((TH2I*)arr->At(1))->Fill(dydx, dy/TMath::Sqrt(cov[0] /*+ sx2*/ + sy2));
407   
408       if(DebugLevel()>=2){
409         // Get z-position with respect to anode wire
410         Int_t istk = fGeo->GetStack(c->GetDetector());
411         AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
412         Float_t row0 = pp->GetRow0();
413         Float_t d  =  row0 - zt + pp->GetAnodeWireOffset();
414         d -= ((Int_t)(2 * d)) / 2.0;
415         if (d > 0.25) d  = 0.5 - d;
416
417         AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
418         fCl->Add(clInfo);
419         clInfo->SetCluster(c);
420         Float_t covF[] = {cov[0], cov[1], cov[2]};
421         clInfo->SetGlobalPosition(yt, zt, dydx, dzdx, covF);
422         clInfo->SetResolution(dy);
423         clInfo->SetAnisochronity(d);
424         clInfo->SetDriftLength(dx);
425         clInfo->SetTilt(tilt);
426         (*DebugStream()) << "ClusterREC"
427           <<"status="  << status
428           <<"clInfo.=" << clInfo
429           << "\n";
430       }
431     }
432   }
433   return (TH2I*)arr->At(0);
434 }
435
436
437 //________________________________________________________
438 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
439 {
440 // Plot normalized residuals for tracklets to track. 
441 // 
442 // We start from the result that if X=N(|m|, |Cov|)
443 // BEGIN_LATEX
444 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
445 // END_LATEX
446 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial 
447 // reference position. 
448   if(track) fkTrack = track;
449   if(!fkTrack){
450     AliDebug(2, "No Track defined.");
451     return NULL;
452   }
453   TObjArray *arr = NULL;
454   if(!(arr = (TObjArray*)fContainer->At(kTrackTRD ))){
455     AliWarning("No output container defined.");
456     return NULL;
457   }
458
459   Double_t cov[3], covR[7]/*, sqr[3], inv[3]*/;
460   Float_t x, dx, dy, dz;
461   AliTRDseedV1 *fTracklet = NULL;  
462   for(Int_t il=AliTRDgeometry::kNlayer; il--;){
463     if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
464     if(!fTracklet->IsOK()) continue;
465     x    = fTracklet->GetX();
466     dx   = fTracklet->GetX0() - x;
467     // compute dy^2 and dz^2
468     dy   = fTracklet->GetYref(0)-dx*fTracklet->GetYref(1) - fTracklet->GetY();
469     dz   = fTracklet->GetZref(0)-dx*fTracklet->GetZref(1) - fTracklet->GetZ();
470     // compute covariance matrix
471     fTracklet->GetCovAt(x, cov);
472     fTracklet->GetCovRef(covR);
473     cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2]; 
474 /*  // Correct PULL calculation by considering off  
475     // diagonal elements in the covariance matrix
476     // compute square root matrix
477     if(AliTRDseedV1::GetCovInv(cov, inv)==0.) continue;
478     if(AliTRDseedV1::GetCovSqrt(inv, sqr)<0.) continue;
479     Double_t y = sqr[0]*dy+sqr[1]*dz;
480     Double_t z = sqr[1]*dy+sqr[2]*dz;
481     ((TH3*)h)->Fill(y, z, fTracklet->GetYref(1));*/
482
483     ((TH2I*)arr->At(0))->Fill(fTracklet->GetYref(1), dy);
484     ((TH2I*)arr->At(1))->Fill(fTracklet->GetYref(1), dy/TMath::Sqrt(cov[0]));
485     ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), TMath::ATan((fTracklet->GetYref(1)-fTracklet->GetYfit(1))/(1-fTracklet->GetYref(1)*fTracklet->GetYfit(1))));
486     if(!fTracklet->IsRowCross()) continue;
487     ((TH2I*)arr->At(2))->Fill(fTracklet->GetZref(1), dz);
488     ((TH2I*)arr->At(3))->Fill(fTracklet->GetZref(1), dz/TMath::Sqrt(cov[2]));
489   }
490
491
492   return (TH2I*)arr->At(0);
493 }
494
495
496 //________________________________________________________
497 TH1* AliTRDresolution::PlotTrackTPC(const AliTRDtrackV1 *track)
498 {
499 // Store resolution/pulls of Kalman before updating with the TRD information 
500 // at the radial position of the first tracklet. The following points are used 
501 // for comparison  
502 //  - the (y,z,snp) of the first TRD tracklet
503 //  - the (y, z, snp, tgl, pt) of the MC track reference
504 // 
505 // Additionally the momentum resolution/pulls are calculated for usage in the 
506 // PID calculation. 
507
508   if(track) fkTrack = track;
509   if(!fkTrack){
510     AliDebug(2, "No Track defined.");
511     return NULL;
512   }
513   AliExternalTrackParam *tin = NULL;
514   if(!(tin = fkTrack->GetTrackLow())){
515     AliWarning("Track did not entered TRD fiducial volume.");
516     return NULL;
517   }
518   TH1 *h = NULL;
519   
520   Double_t x = tin->GetX();
521   AliTRDseedV1 *tracklet = NULL;  
522   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
523     if(!(tracklet = fkTrack->GetTracklet(ily))) continue;
524     break;
525   }
526   if(!tracklet || TMath::Abs(x-tracklet->GetX())>1.e-3){
527     AliWarning("Tracklet did not match TRD entrance.");
528     return NULL;
529   }
530   const Int_t kNPAR(5);
531   Double_t parR[kNPAR]; memcpy(parR, tin->GetParameter(), kNPAR*sizeof(Double_t));
532   Double_t covR[3*kNPAR]; memcpy(covR, tin->GetCovariance(), 3*kNPAR*sizeof(Double_t));
533   Double_t cov[3]; tracklet->GetCovAt(x, cov);
534
535   // define sum covariances
536   TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
537   Double_t *pc = &covR[0], *pp = &parR[0];
538   for(Int_t ir=0; ir<kNPAR; ir++, pp++){
539     PAR(ir) = (*pp);
540     for(Int_t ic = 0; ic<=ir; ic++,pc++){ 
541       COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
542     }
543   }
544   PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
545   //COV.Print(); PAR.Print();
546
547   //TODO Double_t dydx =  TMath::Sqrt(1.-parR[2]*parR[2])/parR[2]; 
548   Double_t dy = parR[0] - tracklet->GetY(); 
549   TObjArray *arr = (TObjArray*)fContainer->At(kTrackTPC);
550   ((TH2I*)arr->At(0))->Fill(tracklet->GetYref(1), dy);
551   ((TH2I*)arr->At(1))->Fill(tracklet->GetYref(1), dy/TMath::Sqrt(COV(0,0)+cov[0]));
552   if(tracklet->IsRowCross()){
553     Double_t dz = parR[1] - tracklet->GetZ(); 
554     ((TH2I*)arr->At(2))->Fill(tracklet->GetZref(1), dz);
555     ((TH2I*)arr->At(3))->Fill(tracklet->GetZref(1), dz/TMath::Sqrt(COV(1,1)+cov[2]));
556   }
557   Double_t dphi = TMath::ASin(PAR[2])-TMath::ATan(tracklet->GetYfit(1));  ((TH2I*)arr->At(4))->Fill(tracklet->GetYref(1), dphi);
558
559
560   // register reference histo for mini-task
561   h = (TH2I*)arr->At(0);
562
563   if(DebugLevel()>=1){
564     (*DebugStream()) << "trackIn"
565       << "x="       << x
566       << "P="       << &PAR
567       << "C="       << &COV
568       << "\n";
569
570     Double_t y = tracklet->GetY(); 
571     Double_t z = tracklet->GetZ(); 
572     (*DebugStream()) << "trackletIn"
573       << "y="       << y
574       << "z="       << z
575       << "Vy="      << cov[0]
576       << "Cyz="     << cov[1]
577       << "Vz="      << cov[2]
578       << "\n";
579   }
580
581
582   if(!HasMCdata()) return h;
583   UChar_t s;
584   Float_t dx, pt0, x0=tracklet->GetX0(), y0, z0, dydx0, dzdx0;
585   if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
586   Int_t pdg = fkMC->GetPDG(),
587         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
588         sign(0);
589   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
590   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
591   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
592
593   // translate to reference radial position
594   dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
595   Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi))
596   //Fill MC info
597   TVectorD PARMC(kNPAR);
598   PARMC[0]=y0; PARMC[1]=z0;
599   PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
600   PARMC[4]=1./pt0;
601
602 //   TMatrixDSymEigen eigen(COV);
603 //   TVectorD evals = eigen.GetEigenValues();
604 //   TMatrixDSym evalsm(kNPAR);
605 //   for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
606 //   TMatrixD evecs = eigen.GetEigenVectors();
607 //   TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
608   
609   // fill histos
610   arr = (TObjArray*)fContainer->At(kMCtrackTPC);
611   // y resolution/pulls
612   ((TH2I*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0]);
613   ((TH2I*)arr->At(1))->Fill(dydx0, (PARMC[0]-PAR[0])/TMath::Sqrt(COV(0,0)));
614   // z resolution/pulls
615   ((TH2I*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1]);
616   ((TH2I*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)));
617   // phi resolution/snp pulls
618   ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
619   ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
620   // theta resolution/tgl pulls
621   ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
622   ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
623   // pt resolution\\1/pt pulls\\p resolution/pull
624   ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., sign*sIdx);
625   ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), sign*sIdx);
626
627   Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
628   p = TMath::Sqrt(1.+ PAR[3]*PAR[3])/PAR[4];
629   Float_t sp = 
630     p*p*PAR[4]*PAR[4]*COV(4,4)
631    +2.*PAR[3]*COV(3,4)/PAR[4]
632    +PAR[3]*PAR[3]*COV(3,3)/p/p/PAR[4]/PAR[4]/PAR[4]/PAR[4];
633   ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., sign*sIdx);
634   if(sp>0.) ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/TMath::Sqrt(sp), sign*sIdx);
635
636   // fill debug for MC 
637   if(DebugLevel()>=1){
638     (*DebugStream()) << "trackInMC"
639       << "P="   << &PARMC
640       << "\n";
641   }
642   return h;
643 }
644
645 //________________________________________________________
646 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
647 {
648   //
649   // Plot MC distributions
650   //
651
652   if(!HasMCdata()){ 
653     AliWarning("No MC defined. Results will not be available.");
654     return NULL;
655   }
656   if(track) fkTrack = track;
657   if(!fkTrack){
658     AliDebug(2, "No Track defined.");
659     return NULL;
660   }
661   // retriev track characteristics
662   Int_t pdg = fkMC->GetPDG(),
663         sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
664         sign(0),
665         det(-1),
666         label(fkMC->GetLabel());
667   if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
668   TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
669   if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
670   Bool_t kBarrel = Bool_t(fkESD->GetStatus() & AliESDtrack::kTRDin);
671
672   TObjArray *arr(NULL);TH1 *h(NULL);
673   UChar_t s;
674   Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl;
675   Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
676   Double_t covR[7]/*, cov[3]*/;
677
678   if(DebugLevel()>=1){
679     TVectorD dX(12), dY(12), dZ(12), Pt(12), dPt(12), cCOV(12*15);
680     fkMC->PropagateKalman(&dX, &dY, &dZ, &Pt, &dPt, &cCOV);
681     (*DebugStream()) << "MCkalman"
682       << "pdg=" << pdg
683       << "dx="  << &dX
684       << "dy="  << &dY
685       << "dz="  << &dZ
686       << "pt="  << &Pt
687       << "dpt=" << &dPt
688       << "cov=" << &cCOV
689       << "\n";
690   }
691
692   AliTRDReconstructor rec;
693   AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
694   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
695     if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
696        !fTracklet->IsOK())*/ continue;
697
698     det = fTracklet->GetDetector();
699     x0  = fTracklet->GetX0();
700     //radial shift with respect to the MC reference (radial position of the pad plane)
701     x= fTracklet->GetX();
702     if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
703     xAnode  = fTracklet->GetX0();
704
705     // MC track position at reference radial position
706     dx  = x0 - x;
707     if(DebugLevel()>=1){
708       (*DebugStream()) << "MC"
709         << "det="     << det
710         << "pdg="     << pdg
711         << "sgn="     << sign
712         << "barrel="  << kBarrel
713         << "pt="      << pt0
714         << "x="       << x0
715         << "y="       << y0
716         << "z="       << z0
717         << "dydx="    << dydx0
718         << "dzdx="    << dzdx0
719         << "\n";
720     }
721     Float_t ymc = y0 - dx*dydx0;
722     Float_t zmc = z0 - dx*dzdx0;
723     //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
724
725     // Kalman position at reference radial position
726     dx = xAnode - x;
727     dydx = fTracklet->GetYref(1);
728     dzdx = fTracklet->GetZref(1);
729     dzdl = fTracklet->GetTgl();
730     y  = fTracklet->GetYref(0) - dx*dydx;
731     dy = y - ymc;
732     z  = fTracklet->GetZref(0) - dx*dzdx;
733     dz = z - zmc;
734     pt = TMath::Abs(fTracklet->GetPt());
735     fTracklet->GetCovRef(covR);
736
737     arr = (TObjArray*)fContainer->At(kMCtrackTRD);
738     // y resolution/pulls
739     ((TH2I*)arr->At(0))->Fill(dydx0, dy);
740     ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(covR[0]));
741     // z resolution/pulls
742     ((TH2I*)arr->At(2))->Fill(dzdx0, dz);
743     ((TH2I*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]));
744     // phi resolution/ snp pulls
745     Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
746     ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
747     Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0);
748     ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3]));
749     // theta resolution/ tgl pulls
750     Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
751               dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
752     ((TH2I*)arr->At(6))->Fill(dzdl0, 
753     TMath::ATan(dtgl));
754     ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
755     // pt resolution  \\ 1/pt pulls \\ p resolution for PID
756     Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
757              p  = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
758     if(kBarrel){
759       ((TH3S*)((TObjArray*)arr->At(8))->At(ily))->Fill(pt0, pt/pt0-1., sign*sIdx);
760       ((TH3S*)((TObjArray*)arr->At(9))->At(ily))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
761       ((TH3S*)((TObjArray*)arr->At(10))->At(ily))->Fill(p0, p/p0-1., sign*sIdx);
762     } else {
763       ((TH3S*)((TObjArray*)arr->At(11))->At(ily))->Fill(pt0, pt/pt0-1., sign*sIdx);
764       ((TH3S*)((TObjArray*)arr->At(12))->At(ily))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
765       ((TH3S*)((TObjArray*)arr->At(13))->At(ily))->Fill(p0, p/p0-1., sign*sIdx);
766     }
767
768     // Fill Debug stream for Kalman track
769     if(DebugLevel()>=1){
770       (*DebugStream()) << "MCtrack"
771         << "pt="      << pt
772         << "x="       << x
773         << "y="       << y
774         << "z="       << z
775         << "dydx="    << dydx
776         << "dzdx="    << dzdx
777         << "s2y="     << covR[0]
778         << "s2z="     << covR[2]
779         << "\n";
780     }
781
782     // recalculate tracklet based on the MC info
783     AliTRDseedV1 tt(*fTracklet);
784     tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
785     tt.SetZref(1, dzdx0); 
786     tt.SetReconstructor(&rec);
787     tt.Fit(kTRUE, kTRUE);
788     x= tt.GetX();y= tt.GetY();z= tt.GetZ();
789     dydx = tt.GetYfit(1);
790     dx = x0 - x;
791     ymc = y0 - dx*dydx0;
792     zmc = z0 - dx*dzdx0;
793     Bool_t rc = tt.IsRowCross(); 
794     
795     // add tracklet residuals for y and dydx
796     arr = (TObjArray*)fContainer->At(kMCtracklet);
797     if(!rc){
798       dy    = y-ymc;
799
800       Float_t dphi  = (dydx - dydx0);
801       dphi /= (1.- dydx*dydx0);
802
803       ((TH3S*)arr->At(0))->Fill(dydx0, dy, pt0);
804       if(tt.GetS2Y()>0.) ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(tt.GetS2Y()));
805       ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
806     } else {
807       // add tracklet residuals for z
808       dz = z-zmc;
809       ((TH2I*)arr->At(2))->Fill(dzdl0, dz);
810       if(tt.GetS2Z()>0.) ((TH2I*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z()));
811     }
812   
813     // Fill Debug stream for tracklet
814     if(DebugLevel()>=1){
815       Float_t s2y = tt.GetS2Y();
816       Float_t s2z = tt.GetS2Z();
817       (*DebugStream()) << "MCtracklet"
818         << "rc="    << rc
819         << "x="     << x
820         << "y="     << y
821         << "z="     << z
822         << "dydx="  << dydx
823         << "s2y="   << s2y
824         << "s2z="   << s2z
825         << "\n";
826     }
827
828     Int_t istk = AliTRDgeometry::GetStack(det); 
829     AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
830     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
831     Float_t tilt = fTracklet->GetTilt();
832     //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
833
834     arr = (TObjArray*)fContainer->At(kMCcluster);
835     AliTRDcluster *c = NULL;
836     fTracklet->ResetClusterIter(kFALSE);
837     while((c = fTracklet->PrevCluster())){
838       Float_t  q = TMath::Abs(c->GetQ());
839       x = c->GetX(); y = c->GetY();z = c->GetZ();
840       dx = x0 - x; 
841       ymc= y0 - dx*dydx0;
842       zmc= z0 - dx*dzdx0;
843       dy = (y - tilt*(z-zmc)) - ymc;
844       // Fill Histograms
845       if(q>20. && q<250.){ 
846         ((TH3S*)arr->At(0))->Fill(dydx0, dy, pt0);
847         ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(c->GetSigmaY2()));
848       }
849
850       // Fill calibration container
851       Float_t d = zr0 - zmc;
852       d -= ((Int_t)(2 * d)) / 2.0;
853       if (d > 0.25) d  = 0.5 - d;
854       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
855       clInfo->SetCluster(c);
856       clInfo->SetMC(pdg, label);
857       clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
858       clInfo->SetResolution(dy);
859       clInfo->SetAnisochronity(d);
860       clInfo->SetDriftLength(dx-.5*AliTRDgeometry::CamHght());
861       clInfo->SetTilt(tilt);
862       fMCcl->Add(clInfo);
863       if(DebugLevel()>=2){ 
864         if(!clInfoArr) clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
865         clInfoArr->Add(clInfo);
866       }
867     }
868     // Fill Debug Tree
869     if(DebugLevel()>=2 && clInfoArr){
870       (*DebugStream()) << "MCcluster"
871         <<"clInfo.=" << clInfoArr
872         << "\n";
873       clInfoArr->Clear();
874     }
875   }
876   if(clInfoArr) delete clInfoArr;
877   return h;
878 }
879
880
881 //________________________________________________________
882 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
883 {
884   //
885   // Get the reference figures
886   //
887
888   Float_t xy[4] = {0., 0., 0., 0.};
889   if(!gPad){
890     AliWarning("Please provide a canvas to draw results.");
891     return kFALSE;
892   }
893   Int_t first(2),     // first particle species to be drawn
894         nspecies(7); // last particle species to be drawn
895   TList *l = NULL; TVirtualPad *pad=NULL;
896   switch(ifig){
897   case kCharge:
898     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
899     ((TVirtualPad*)l->At(0))->cd();
900     ((TGraphAsymmErrors*)((TObjArray*)fGraphM->At(kCharge))->At(0))->Draw("apl");
901     ((TVirtualPad*)l->At(1))->cd();
902     ((TGraphErrors*)((TObjArray*)fGraphS->At(kCharge))->At(0))->Draw("apl");
903     break;
904   case kCluster:
905     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
906     xy[0] = -.3; xy[1] = -200.; xy[2] = .3; xy[3] = 1000.;
907     ((TVirtualPad*)l->At(0))->cd();
908     if(!GetGraphPlot(&xy[0], kCluster, 0)) break;
909     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
910     ((TVirtualPad*)l->At(1))->cd();
911     if(!GetGraphPlot(&xy[0], kCluster, 1)) break;
912     return kTRUE;
913   case kTrackTRD :
914     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
915     xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
916     ((TVirtualPad*)l->At(0))->cd();
917     if(!GetGraphPlot(&xy[0], kTrackTRD , 0)) break;
918     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
919     ((TVirtualPad*)l->At(1))->cd();
920     if(!GetGraphPlot(&xy[0], kTrackTRD , 1)) break;
921     return kTRUE;
922   case 3: // kTrackTRD  z
923     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
924     xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
925     ((TVirtualPad*)l->At(0))->cd();
926     if(!GetGraphPlot(&xy[0], kTrackTRD , 2)) break;
927     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
928     ((TVirtualPad*)l->At(1))->cd();
929     if(!GetGraphPlot(&xy[0], kTrackTRD , 3)) break;
930     return kTRUE;
931   case 4: // kTrackTRD  phi
932     xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
933     if(GetGraphPlot(&xy[0], kTrackTRD , 4)) return kTRUE;
934     break;
935   case 5: // kTrackTPC y
936     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
937     xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
938     pad = ((TVirtualPad*)l->At(0)); pad->cd();
939     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
940     if(!GetGraphPlot(&xy[0], kTrackTPC, 0)) break;
941     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
942     pad=((TVirtualPad*)l->At(1)); pad->cd();
943     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
944     if(!GetGraphPlot(&xy[0], kTrackTPC, 1)) break;
945     return kTRUE;
946   case 6: // kTrackTPC z
947     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
948     xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
949     pad = ((TVirtualPad*)l->At(0)); pad->cd();
950     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
951     if(!GetGraphPlot(&xy[0], kTrackTPC, 2)) break;
952     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
953     pad = ((TVirtualPad*)l->At(1)); pad->cd();
954     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
955     if(!GetGraphPlot(&xy[0], kTrackTPC, 3)) break;
956     return kTRUE;
957   case 7: // kTrackTPC phi
958     xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
959     if(GetGraphPlot(&xy[0], kTrackTPC, 4)) return kTRUE;
960     break;
961   case 8: // kMCcluster pt <2 GeV/c
962     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
963     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
964     ((TVirtualPad*)l->At(0))->cd();
965     if(!GetGraphPlot(&xy[0], kMCcluster, 0)) break;
966     ((TVirtualPad*)l->At(1))->cd();
967     if(!GetGraphPlot(&xy[0], kMCcluster, 1)) break;
968     return kTRUE;
969   case 9: // kMCcluster pt > 3 GeV/c
970     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
971     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
972     ((TVirtualPad*)l->At(0))->cd();
973     if(!GetGraphPlot(&xy[0], kMCcluster, 2)) break;
974     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
975     ((TVirtualPad*)l->At(1))->cd();
976     if(!GetGraphPlot(&xy[0], kMCcluster, 3)) break;
977     return kTRUE;
978   case 10: //kMCtracklet [y] pt < 3 GeV/c
979     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
980     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =250.;
981     ((TVirtualPad*)l->At(0))->cd();
982     if(!GetGraphPlot(&xy[0], kMCtracklet, 0)) break;
983     ((TVirtualPad*)l->At(1))->cd();
984     if(!GetGraphPlot(&xy[0], kMCtracklet, 1)) break;
985     return kTRUE;
986   case 11: //kMCtracklet [y] pt > 3 GeV/c
987     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
988     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =250.;
989     ((TVirtualPad*)l->At(0))->cd();
990     if(!GetGraphPlot(&xy[0], kMCtracklet, 2)) break;
991     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
992     ((TVirtualPad*)l->At(1))->cd();
993     if(!GetGraphPlot(&xy[0], kMCtracklet, 3)) break;
994     return kTRUE;
995   case 12: //kMCtracklet [z]
996     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
997     xy[0]=-1.; xy[1]=-100.; xy[2]=1.; xy[3] =2500.;
998     ((TVirtualPad*)l->At(0))->cd();
999     if(!GetGraphPlot(&xy[0], kMCtracklet, 4)) break;
1000     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1001     ((TVirtualPad*)l->At(1))->cd();
1002     if(!GetGraphPlot(&xy[0], kMCtracklet, 5)) break;
1003     return kTRUE;
1004   case 13: //kMCtracklet [phi]
1005     xy[0]=-.3; xy[1]=-3.; xy[2]=.3; xy[3] =25.;
1006     if(!GetGraphPlot(&xy[0], kMCtracklet, 6)) break;
1007     return kTRUE;
1008   case 14: //kMCtrackTRD [y]
1009     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1010     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1011     ((TVirtualPad*)l->At(0))->cd();
1012     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 0)) break;
1013     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 3.5;
1014     ((TVirtualPad*)l->At(1))->cd();
1015     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 1)) break;
1016     return kTRUE;
1017   case 15: //kMCtrackTRD [z]
1018     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1019     xy[0]=-1.; xy[1]=-700.; xy[2]=1.; xy[3] =1500.;
1020     ((TVirtualPad*)l->At(0))->cd();
1021     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 2)) break;
1022     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1023     ((TVirtualPad*)l->At(1))->cd();
1024     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 3)) break;
1025     return kTRUE;
1026   case 16: //kMCtrackTRD [phi/snp]
1027     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1028     xy[0]=-.2; xy[1]=-0.2; xy[2]=.2; xy[3] =2.;
1029     ((TVirtualPad*)l->At(0))->cd();
1030     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 4)) break;
1031     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
1032     ((TVirtualPad*)l->At(1))->cd();
1033     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 5)) break;
1034     return kTRUE;
1035   case 17: //kMCtrackTRD [theta/tgl]
1036     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1037     xy[0]=-1.; xy[1]=-0.5; xy[2]=1.; xy[3] =5.;
1038     ((TVirtualPad*)l->At(0))->cd();
1039     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 6)) break;
1040     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
1041     ((TVirtualPad*)l->At(1))->cd();
1042     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 7)) break;
1043     return kTRUE;
1044   case 18: //kMCtrackTRD [pt]
1045     xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 4.;
1046     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1047     pad = (TVirtualPad*)l->At(0); pad->cd();
1048     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1049     if(!GetGraphTrack(&xy[0], 8, first, nspecies)) break;
1050     pad->Modified(); pad->Update(); pad->SetLogx();
1051     pad = (TVirtualPad*)l->At(1); pad->cd();
1052     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1053     if(!GetGraphTrack(&xy[0], 8, 55+first, nspecies, kTRUE)) break;
1054     pad->Modified(); pad->Update(); pad->SetLogx();
1055     return kTRUE;
1056   case 19: //kMCtrackTRD [1/pt] pulls
1057     xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1058     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1059     pad = (TVirtualPad*)l->At(0); pad->cd();
1060     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1061     if(!GetGraphTrack(&xy[0], 9, first, nspecies)) break;
1062     pad = (TVirtualPad*)l->At(1); pad->cd();
1063     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1064     if(!GetGraphTrack(&xy[0], 9, 55+first, nspecies, kTRUE)) break;
1065     return kTRUE;
1066   case 20: //kMCtrackTRD [p]
1067     xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 4.;
1068     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1069     pad = (TVirtualPad*)l->At(0); pad->cd();
1070     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1071     if(!GetGraphTrack(&xy[0], 10, first, nspecies)) break;
1072     pad->Modified(); pad->Update(); pad->SetLogx();
1073     pad = (TVirtualPad*)l->At(1); pad->cd();
1074     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1075     if(!GetGraphTrack(&xy[0], 10, 55+first, nspecies, kTRUE)) break;
1076     pad->Modified(); pad->Update(); pad->SetLogx();
1077     return kTRUE;
1078   case 21: //kMCtrackTRD - SA [pt]
1079     xy[0] = 0.; xy[1] = -5.; xy[2] = 12.; xy[3] = 7.;
1080     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1081     pad = (TVirtualPad*)l->At(0); pad->cd();
1082     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1083     if(!GetGraphTrack(&xy[0], 11, first, nspecies)) break;
1084     pad->SetLogx();
1085     pad = (TVirtualPad*)l->At(1); pad->cd();
1086     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1087     if(!GetGraphTrack(&xy[0], 11, 55+first, nspecies, kTRUE)) break;
1088     pad->SetLogx();
1089     return kTRUE;
1090   case 22: //kMCtrackTRD - SA [1/pt] pulls
1091     xy[0] = 0.; xy[1] = -1.5; xy[2] = 2.; xy[3] = 2.;
1092     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1093     pad = (TVirtualPad*)l->At(0); pad->cd();
1094     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1095     if(!GetGraphTrack(&xy[0], 12, first, nspecies)) break;
1096     pad = (TVirtualPad*)l->At(1); pad->cd();
1097     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1098     if(!GetGraphTrack(&xy[0], 12, 55+first, nspecies, kTRUE)) break;
1099     return kTRUE;
1100   case 23: //kMCtrackTRD - SA [p]
1101     xy[0] = 0.; xy[1] = -7.5; xy[2] = 12.; xy[3] = 10.5;
1102     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1103     pad = (TVirtualPad*)l->At(0); pad->cd();
1104     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1105     if(!GetGraphTrack(&xy[0], 13, first, nspecies)) break;
1106     pad->SetLogx();
1107     pad = (TVirtualPad*)l->At(1); pad->cd();
1108     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1109     if(!GetGraphTrack(&xy[0], 13, 55+first, nspecies, kTRUE)) break;
1110     pad->SetLogx();
1111     return kTRUE;
1112   case 24: // kMCtrackTPC [y]
1113     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1114     xy[0]=-.25; xy[1]=-50.; xy[2]=.25; xy[3] =800.;
1115     ((TVirtualPad*)l->At(0))->cd();
1116     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 0)) break;
1117     xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 2.5;
1118     ((TVirtualPad*)l->At(1))->cd();
1119     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 1)) break;
1120     return kTRUE;
1121   case 25: // kMCtrackTPC [z]
1122     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1123     xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =800.;
1124     ((TVirtualPad*)l->At(0))->cd();
1125     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 2)) break;
1126     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1127     ((TVirtualPad*)l->At(1))->cd();
1128     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 3)) break;
1129     return kTRUE;
1130   case 26: // kMCtrackTPC [phi|snp]
1131     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1132     xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
1133     ((TVirtualPad*)l->At(0))->cd();
1134     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 4)) break;
1135     xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
1136     ((TVirtualPad*)l->At(1))->cd();
1137     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 5)) break;
1138     return kTRUE;
1139   case 27: // kMCtrackTPC [theta|tgl]
1140     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1141     xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.;
1142     ((TVirtualPad*)l->At(0))->cd();
1143     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 6)) break;
1144     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 1.5;
1145     ((TVirtualPad*)l->At(1))->cd();
1146     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 7)) break;
1147     return kTRUE;
1148   case 28: // kMCtrackTPC [pt]
1149     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1150     xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1151     pad=(TVirtualPad*)l->At(0); pad->cd();
1152     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1153     if(!GetGraphTrackTPC(xy, 8, first, nspecies)) break;
1154     pad->SetLogx();
1155     xy[0]=0.; xy[1]=-0.5; xy[2]=2.; xy[3] =3.;
1156     pad = (TVirtualPad*)l->At(1); pad->cd();
1157     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1158     if(!GetGraphTrackTPC(xy, 9, first, nspecies, kTRUE)) break;
1159     return kTRUE;
1160   case 29: // kMCtrackTPC [p]
1161     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1162     xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1163     pad = ((TVirtualPad*)l->At(0));pad->cd();
1164     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1165     if(!GetGraphTrackTPC(xy, 10, first, nspecies)) break;
1166     pad->SetLogx();
1167     xy[0]=0.; xy[1]=-0.5; xy[2]=8.; xy[3] =2.5;
1168     pad = ((TVirtualPad*)l->At(1)); pad->cd();
1169     pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1170     if(!GetGraphTrackTPC(xy, 11, first, nspecies, kTRUE)) break;
1171     return kTRUE;
1172   case 30:  // kMCtrackTOF [z]
1173     return kTRUE;
1174   }
1175   AliWarning(Form("Reference plot [%d] missing result", ifig));
1176   return kFALSE;
1177 }
1178
1179 Char_t const *fgParticle[11]={
1180   " p bar", " K -", " #pi -", " #mu -", " e -",
1181   " No PID",
1182   " e +", " #mu +", " #pi +", " K +", " p",
1183 };
1184 const Color_t fgColorS[11]={
1185 kOrange, kOrange-3, kMagenta+1, kViolet, kRed,
1186 kGray,
1187 kRed, kViolet, kMagenta+1, kOrange-3, kOrange
1188 };
1189 const Color_t fgColorM[11]={
1190 kCyan-5, kAzure-4, kBlue-7, kBlue+2, kViolet+10,
1191 kBlack,
1192 kViolet+10, kBlue+2, kBlue-7, kAzure-4, kCyan-5
1193 };
1194 const Marker_t fgMarker[11]={
1195 30, 30, 26, 25, 24,
1196 28,
1197 20, 21, 22, 29, 29
1198 };
1199 //________________________________________________________
1200 Bool_t AliTRDresolution::PostProcess()
1201 {
1202   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
1203   if (!fContainer) {
1204     AliError("ERROR: list not available");
1205     return kFALSE;
1206   }
1207   Int_t nc(0);
1208   TGraph *gm= NULL, *gs= NULL;
1209   if(!fGraphS && !fGraphM){ 
1210     TObjArray *aM(NULL), *aS(NULL);
1211     Int_t n = fContainer->GetEntriesFast();
1212     fGraphS = new TObjArray(n); fGraphS->SetOwner();
1213     fGraphM = new TObjArray(n); fGraphM->SetOwner();
1214     for(Int_t ig=0; ig<n; ig++){
1215       fGraphM->AddAt(aM = new TObjArray(fgNproj[ig]), ig);
1216       fGraphS->AddAt(aS = new TObjArray(fgNproj[ig]), ig);
1217
1218       for(Int_t ic=0; ic<fgNproj[ig]; ic++){
1219         //printf("g[%d] c[%d] n[%d]\n", ig, ic, fgNcomp[nc]);
1220         if(fgNcomp[nc]>1){
1221           TObjArray *agS(NULL), *agM(NULL);
1222           aS->AddAt(agS = new TObjArray(fgNcomp[nc]), ic); 
1223           aM->AddAt(agM = new TObjArray(fgNcomp[nc]), ic); 
1224           for(Int_t is=fgNcomp[nc]; is--;){
1225             agS->AddAt(gs = new TGraphErrors(), is);
1226             Int_t is0(is%11);
1227             gs->SetMarkerStyle(fgMarker[is0]);
1228             gs->SetMarkerColor(fgColorS[is0]);
1229             gs->SetLineColor(fgColorS[is0]);
1230             gs->SetNameTitle(Form("s_%d%02d%02d", ig, ic, is), fgParticle[is0]);
1231
1232             agM->AddAt(gm = new TGraphErrors(), is);
1233             gm->SetMarkerStyle(fgMarker[is0]);
1234             gm->SetMarkerColor(fgColorM[is0]);
1235             gm->SetLineColor(fgColorM[is0]);
1236             gm->SetNameTitle(Form("m_%d%02d%02d", ig, ic, is), fgParticle[is0]);
1237           }
1238         } else {
1239           aS->AddAt(gs = new TGraphErrors(), ic);
1240           gs->SetMarkerStyle(23);
1241           gs->SetMarkerColor(kRed);
1242           gs->SetLineColor(kRed);
1243           gs->SetNameTitle(Form("s_%d%02d", ig, ic), "");
1244   
1245           aM->AddAt(gm = ig ? (TGraph*)new TGraphErrors() : (TGraph*)new TGraphAsymmErrors(), ic);
1246           gm->SetLineColor(kBlack);
1247           gm->SetMarkerStyle(7);
1248           gm->SetMarkerColor(kBlack);
1249           gm->SetNameTitle(Form("m_%d%02d", ig, ic), "");
1250         }
1251         nc+=fgNcomp[ic];
1252       }
1253     }
1254   }
1255
1256 /*  printf("\n\n\n"); fGraphS->ls();
1257   printf("\n\n\n"); fGraphM->ls();*/
1258   
1259
1260   // DEFINE MODELS
1261   // simple gauss
1262   TF1 fg("fGauss", "gaus", -.5, .5);  
1263   // Landau for charge resolution
1264   TF1 fch("fClCh", "landau", 0., 1000.);  
1265   // Landau for e+- pt resolution
1266   TF1 fpt("fPt", "landau", -0.1, 0.2);  
1267
1268   //PROCESS EXPERIMENTAL DISTRIBUTIONS
1269   // Charge resolution
1270   //Process3DL(kCharge, 0, &fl); 
1271   // Clusters residuals
1272   Process2D(kCluster, 0, &fg, 1.e4); 
1273   Process2D(kCluster, 1, &fg); 
1274   fNRefFigures = 1;
1275   // Tracklet residual/pulls
1276   Process2D(kTrackTRD , 0, &fg, 1.e4); 
1277   Process2D(kTrackTRD , 1, &fg); 
1278   Process2D(kTrackTRD , 2, &fg, 1.e4); 
1279   Process2D(kTrackTRD , 3, &fg); 
1280   Process2D(kTrackTRD , 4, &fg, 1.e3); 
1281   fNRefFigures = 4;
1282   // TPC track residual/pulls
1283   Process2D(kTrackTPC, 0, &fg, 1.e4); 
1284   Process2D(kTrackTPC, 1, &fg); 
1285   Process2D(kTrackTPC, 2, &fg, 1.e4); 
1286   Process2D(kTrackTPC, 3, &fg); 
1287   Process2D(kTrackTPC, 4, &fg, 1.e3); 
1288   fNRefFigures = 7;
1289
1290   if(!HasMCdata()) return kTRUE;
1291
1292
1293   //PROCESS MC RESIDUAL DISTRIBUTIONS
1294
1295   // CLUSTER Y RESOLUTION/PULLS
1296   Process3Drange(kMCcluster, 0, 0, &fg, 1.e4, 1, 1);
1297   Process3Drange(kMCcluster, 0, 1, &fg, 1.e4, 2, 2);
1298   Process3Drange(kMCcluster, 0, 2, &fg, 1.e4, 3, 12);
1299   Process2D(kMCcluster, 1, &fg, 1., 3);
1300   fNRefFigures = 10;
1301
1302   // TRACKLET RESOLUTION/PULLS
1303   Process3Drange(kMCtracklet, 0, 0, &fg, 1.e4, 1, 1); // y
1304   Process3Drange(kMCtracklet, 0, 1, &fg, 1.e4, 2, 2); // y
1305   Process3Drange(kMCtracklet, 0, 2, &fg, 1.e4, 3, 12); // y
1306   Process2D(kMCtracklet, 1, &fg, 1., 3);       // y pulls
1307   Process2D(kMCtracklet, 2, &fg, 1.e4, 4); // z
1308   Process2D(kMCtracklet, 3, &fg, 1., 5);       // z pulls
1309   Process2D(kMCtracklet, 4, &fg, 1.e3, 6); // phi
1310   fNRefFigures = 14;
1311
1312   // TRACK RESOLUTION/PULLS
1313   Process2D(kMCtrackTRD, 0, &fg, 1.e4);   // y
1314   Process2D(kMCtrackTRD, 1, &fg);         // y PULL
1315   Process2D(kMCtrackTRD, 2, &fg, 1.e4);   // z
1316   Process2D(kMCtrackTRD, 3, &fg);         // z PULL
1317   Process2D(kMCtrackTRD, 4, &fg, 1.e3);   // phi
1318   Process2D(kMCtrackTRD, 5, &fg);         // snp PULL
1319   Process2D(kMCtrackTRD, 6, &fg, 1.e3);   // theta
1320   Process2D(kMCtrackTRD, 7, &fg);         // tgl PULL
1321   Process4D(kMCtrackTRD, 8, &fg, 1.e2);   // pt resolution
1322   Process4D(kMCtrackTRD, 8, &fpt, 1.e2, 4);// pt resolution e1- @ L0
1323   Process4D(kMCtrackTRD, 8, &fpt, 1.e2, 6);// pt resolution e1+ @ L0
1324   Process4D(kMCtrackTRD, 8, &fpt, 1.e2, 55+4);// pt resolution e1- @ L5
1325   Process4D(kMCtrackTRD, 8, &fpt, 1.e2, 55+6);// pt resolution e1+ @ L5
1326   Process4D(kMCtrackTRD, 9, &fg);         // 1/pt pulls
1327   Process4D(kMCtrackTRD, 10, &fg, 1.e2);  // p resolution
1328   Process4D(kMCtrackTRD, 10, &fpt, 1.e2, 4);// p resolution e1- @ L0
1329   Process4D(kMCtrackTRD, 10, &fpt, 1.e2, 6);// p resolution e1+ @ L0
1330   Process4D(kMCtrackTRD, 10, &fpt, 1.e2, 55+4);// p resolution e1- @ L5
1331   Process4D(kMCtrackTRD, 10, &fpt, 1.e2, 55+6);// p resolution e1+ @ L5
1332   Process4D(kMCtrackTRD, 11, &fg, 1.e2);   // pt resolution stand alone
1333   Process4D(kMCtrackTRD, 12, &fg);         // 1/pt pulls  stand alone
1334   Process4D(kMCtrackTRD, 13, &fg, 1.e2);  // p resolution  stand alone
1335   fNRefFigures = 24;
1336
1337   // TRACK TPC RESOLUTION/PULLS
1338   Process2D(kMCtrackTPC, 0, &fg, 1.e4);// y resolution
1339   Process2D(kMCtrackTPC, 1, &fg);      // y pulls
1340   Process2D(kMCtrackTPC, 2, &fg, 1.e4);// z resolution
1341   Process2D(kMCtrackTPC, 3, &fg);      // z pulls
1342   Process2D(kMCtrackTPC, 4, &fg, 1.e3);// phi resolution
1343   Process2D(kMCtrackTPC, 5, &fg);      // snp pulls
1344   Process2D(kMCtrackTPC, 6, &fg, 1.e3);// theta resolution
1345   Process2D(kMCtrackTPC, 7, &fg);      // tgl pulls
1346   Process3D(kMCtrackTPC, 8, &fg, 1.e2);// pt resolution
1347   Process3D(kMCtrackTPC, 9, &fg);      // 1/pt pulls
1348   Process3D(kMCtrackTPC, 10, &fg, 1.e2);// p resolution
1349   Process3D(kMCtrackTPC, 11, &fg);      // p pulls
1350   fNRefFigures = 30;
1351
1352   // TRACK HMPID RESOLUTION/PULLS
1353   Process2D(kMCtrackTOF, 0, &fg, 1.e4); // z towards TOF
1354   Process2D(kMCtrackTOF, 1, &fg);       // z towards TOF
1355   fNRefFigures = 31;
1356
1357   return kTRUE;
1358 }
1359
1360
1361 //________________________________________________________
1362 void AliTRDresolution::Terminate(Option_t *opt)
1363 {
1364   AliTRDrecoTask::Terminate(opt);
1365   if(HasPostProcess()) PostProcess();
1366 }
1367
1368 //________________________________________________________
1369 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1370 {
1371 // Helper function to avoid duplication of code
1372 // Make first guesses on the fit parameters
1373
1374   // find the intial parameters of the fit !! (thanks George)
1375   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1376   Double_t sum = 0.;
1377   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1378   f->SetParLimits(0, 0., 3.*sum);
1379   f->SetParameter(0, .9*sum);
1380   Double_t rms = h->GetRMS();
1381   f->SetParLimits(1, -rms, rms);
1382   f->SetParameter(1, h->GetMean());
1383
1384   f->SetParLimits(2, 0., 2.*rms);
1385   f->SetParameter(2, rms);
1386   if(f->GetNpar() <= 4) return;
1387
1388   f->SetParLimits(3, 0., sum);
1389   f->SetParameter(3, .1*sum);
1390
1391   f->SetParLimits(4, -.3, .3);
1392   f->SetParameter(4, 0.);
1393
1394   f->SetParLimits(5, 0., 1.e2);
1395   f->SetParameter(5, 2.e-1);
1396 }
1397
1398 //________________________________________________________
1399 TObjArray* AliTRDresolution::Histos()
1400 {
1401   //
1402   // Define histograms
1403   //
1404
1405   if(fContainer) return fContainer;
1406
1407   fContainer  = new TObjArray(kNviews);
1408   //fContainer->SetOwner(kTRUE);
1409   TH1 *h(NULL);
1410   TObjArray *arr(NULL);
1411
1412   // binnings for plots containing momentum or pt
1413   const Int_t kNpt(14), kNphi(48), kNdy(60);
1414   Float_t Phi=-.48, Dy=-.3, Pt=0.1;
1415   Float_t binsPhi[kNphi+1], binsDy[kNdy+1], binsPt[kNpt+1];
1416   for(Int_t i=0; i<kNphi+1; i++,Phi+=.02) binsPhi[i]=Phi;
1417   for(Int_t i=0; i<kNdy+1; i++,Dy+=.01) binsDy[i]=Dy;
1418   for(Int_t i=0;i<kNpt+1; i++,Pt=TMath::Exp(i*.15)-1.) binsPt[i]=Pt;
1419
1420   // cluster to track residuals/pulls
1421   fContainer->AddAt(arr = new TObjArray(fgNhistos[kCharge]), kCharge);
1422   arr->SetName("Charge");
1423   if(!(h = (TH3S*)gROOT->FindObject("hCharge"))){
1424     h = new TH3S("hCharge", "Charge Resolution", 20, 1., 2., 24, 0., 3.6, 100, 0., 500.);
1425     h->GetXaxis()->SetTitle("dx/dx_{0}");
1426     h->GetYaxis()->SetTitle("x_{d} [cm]");
1427     h->GetZaxis()->SetTitle("dq/dx [ADC/cm]");
1428   } else h->Reset();
1429   arr->AddAt(h, 0);
1430
1431   // cluster to track residuals/pulls
1432   fContainer->AddAt(arr = new TObjArray(fgNhistos[kCluster]), kCluster);
1433   arr->SetName("Cl");
1434   if(!(h = (TH2I*)gROOT->FindObject("hCl"))){
1435     h = new TH2I("hCl", "Cluster Residuals", 21, -.33, .33, 100, -.5, .5);
1436     h->GetXaxis()->SetTitle("tg(#phi)");
1437     h->GetYaxis()->SetTitle("#Delta y [cm]");
1438     h->GetZaxis()->SetTitle("entries");
1439   } else h->Reset();
1440   arr->AddAt(h, 0);
1441   if(!(h = (TH2I*)gROOT->FindObject("hClpull"))){
1442     h = new TH2I("hClpull", "Cluster Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1443     h->GetXaxis()->SetTitle("tg(#phi)");
1444     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1445     h->GetZaxis()->SetTitle("entries");
1446   } else h->Reset();
1447   arr->AddAt(h, 1);
1448
1449   // tracklet to track residuals/pulls in y direction
1450   fContainer->AddAt(arr = new TObjArray(fgNhistos[kTrackTRD ]), kTrackTRD );
1451   arr->SetName("Trklt");
1452   if(!(h = (TH2I*)gROOT->FindObject("hTrkltY"))){
1453     h = new TH2I("hTrkltY", "Tracklet Y Residuals", 21, -.33, .33, 100, -.5, .5);
1454     h->GetXaxis()->SetTitle("#tg(#phi)");
1455     h->GetYaxis()->SetTitle("#Delta y [cm]");
1456     h->GetZaxis()->SetTitle("entries");
1457   } else h->Reset();
1458   arr->AddAt(h, 0);
1459   if(!(h = (TH2I*)gROOT->FindObject("hTrkltYpull"))){
1460     h = new TH2I("hTrkltYpull", "Tracklet Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1461     h->GetXaxis()->SetTitle("#tg(#phi)");
1462     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1463     h->GetZaxis()->SetTitle("entries");
1464   } else h->Reset();
1465   arr->AddAt(h, 1);
1466   // tracklet to track residuals/pulls in z direction
1467   if(!(h = (TH2I*)gROOT->FindObject("hTrkltZ"))){
1468     h = new TH2I("hTrkltZ", "Tracklet Z Residuals", 50, -1., 1., 100, -1.5, 1.5);
1469     h->GetXaxis()->SetTitle("#tg(#theta)");
1470     h->GetYaxis()->SetTitle("#Delta z [cm]");
1471     h->GetZaxis()->SetTitle("entries");
1472   } else h->Reset();
1473   arr->AddAt(h, 2);
1474   if(!(h = (TH2I*)gROOT->FindObject("hTrkltZpull"))){
1475     h = new TH2I("hTrkltZpull", "Tracklet Z Pulls", 50, -1., 1., 100, -5.5, 5.5);
1476     h->GetXaxis()->SetTitle("#tg(#theta)");
1477     h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
1478     h->GetZaxis()->SetTitle("entries");
1479   } else h->Reset();
1480   arr->AddAt(h, 3);
1481   // tracklet to track phi residuals
1482   if(!(h = (TH2I*)gROOT->FindObject("hTrkltPhi"))){
1483     h = new TH2I("hTrkltPhi", "Tracklet #phi Residuals", 21, -.33, .33, 100, -.5, .5);
1484     h->GetXaxis()->SetTitle("tg(#phi)");
1485     h->GetYaxis()->SetTitle("#Delta phi [rad]");
1486     h->GetZaxis()->SetTitle("entries");
1487   } else h->Reset();
1488   arr->AddAt(h, 4);
1489
1490
1491   // tracklet to TPC track residuals/pulls in y direction
1492   fContainer->AddAt(arr = new TObjArray(fgNhistos[kTrackTPC]), kTrackTPC);
1493   arr->SetName("TrkTPC");
1494   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCY"))){
1495     h = new TH2I("hTrkTPCY", "Track[TPC] Y Residuals", 21, -.33, .33, 100, -.5, .5);
1496     h->GetXaxis()->SetTitle("#tg(#phi)");
1497     h->GetYaxis()->SetTitle("#Delta y [cm]");
1498     h->GetZaxis()->SetTitle("entries");
1499   } else h->Reset();
1500   arr->AddAt(h, 0);
1501   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCYpull"))){
1502     h = new TH2I("hTrkTPCYpull", "Track[TPC] Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1503     h->GetXaxis()->SetTitle("#tg(#phi)");
1504     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1505     h->GetZaxis()->SetTitle("entries");
1506   } else h->Reset();
1507   arr->AddAt(h, 1);
1508   // tracklet to TPC track residuals/pulls in z direction
1509   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZ"))){
1510     h = new TH2I("hTrkTPCZ", "Track[TPC] Z Residuals", 50, -1., 1., 100, -1.5, 1.5);
1511     h->GetXaxis()->SetTitle("#tg(#theta)");
1512     h->GetYaxis()->SetTitle("#Delta z [cm]");
1513     h->GetZaxis()->SetTitle("entries");
1514   } else h->Reset();
1515   arr->AddAt(h, 2);
1516   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZpull"))){
1517     h = new TH2I("hTrkTPCZpull", "Track[TPC] Z Pulls", 50, -1., 1., 100, -5.5, 5.5);
1518     h->GetXaxis()->SetTitle("#tg(#theta)");
1519     h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
1520     h->GetZaxis()->SetTitle("entries");
1521   } else h->Reset();
1522   arr->AddAt(h, 3);
1523   // tracklet to TPC track phi residuals
1524   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCPhi"))){
1525     h = new TH2I("hTrkTPCPhi", "Track[TPC] #phi Residuals", 21, -.33, .33, 100, -.5, .5);
1526     h->GetXaxis()->SetTitle("tg(#phi)");
1527     h->GetYaxis()->SetTitle("#Delta phi [rad]");
1528     h->GetZaxis()->SetTitle("entries");
1529   } else h->Reset();
1530   arr->AddAt(h, 4);
1531
1532
1533   // Resolution histos
1534   if(!HasMCdata()) return fContainer;
1535
1536   // cluster y resolution [0]
1537   fContainer->AddAt(arr = new TObjArray(fgNhistos[kMCcluster]), kMCcluster);
1538   arr->SetName("McCl");
1539   if(!(h = (TH3S*)gROOT->FindObject("hMcCl"))){
1540     h = new TH3S("hMcCl", 
1541     "Cluster Resolution;tg(#phi);#Delta y [cm];p_{t} [GeV/c]", 
1542     kNphi, binsPhi, kNdy, binsDy, kNpt, binsPt);
1543   } else h->Reset();
1544   arr->AddAt(h, 0);
1545   if(!(h = (TH2I*)gROOT->FindObject("hMcClPull"))){
1546     h = new TH2I("hMcClPull", "Cluster Pulls", 48, -.48, .48, 100, -4.5, 4.5);
1547     h->GetXaxis()->SetTitle("tg(#phi)");
1548     h->GetYaxis()->SetTitle("#Deltay/#sigma_{y}");
1549     h->GetZaxis()->SetTitle("p_{t} [GeV/c]");
1550   } else h->Reset();
1551   arr->AddAt(h, 1);
1552
1553
1554   // TRACKLET RESOLUTION
1555   fContainer->AddAt(arr = new TObjArray(fgNhistos[kMCtracklet]), kMCtracklet);
1556   arr->SetName("McTrklt");
1557   // tracklet y resolution
1558   if(!(h = (TH3S*)gROOT->FindObject("hMcTrkltY"))){
1559     h = new TH3S("hMcTrkltY", 
1560     "Tracklet Y Resolution;tg(#phi);#Delta y [cm];p_{t} [GeV/c]", 
1561     kNphi, binsPhi, kNdy, binsDy, kNpt, binsPt);
1562   } else h->Reset();
1563   arr->AddAt(h, 0);
1564   // tracklet y pulls
1565   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltYPull"))){
1566     h = new TH2I("hMcTrkltYPull", "Tracklet Pulls (Y)", 48, -.48, .48, 100, -4.5, 4.5);
1567     h->GetXaxis()->SetTitle("tg(#phi)");
1568     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1569     h->GetZaxis()->SetTitle("entries");
1570   } else h->Reset();
1571   arr->AddAt(h, 1);
1572   // tracklet z resolution
1573   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltZ"))){
1574     h = new TH2I("hMcTrkltZ", "Tracklet Resolution (Z)", 100, -1., 1., 100, -1., 1.);
1575     h->GetXaxis()->SetTitle("tg(#theta)");
1576     h->GetYaxis()->SetTitle("#Delta z [cm]");
1577     h->GetZaxis()->SetTitle("entries");
1578   } else h->Reset();
1579   arr->AddAt(h, 2);
1580   // tracklet z pulls
1581   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltZPull"))){
1582     h = new TH2I("hMcTrkltZPull", "Tracklet Pulls (Z)", 100, -1., 1., 100, -3.5, 3.5);
1583     h->GetXaxis()->SetTitle("tg(#theta)");
1584     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1585     h->GetZaxis()->SetTitle("entries");
1586   } else h->Reset();
1587   arr->AddAt(h, 3);
1588   // tracklet phi resolution
1589   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltPhi"))){
1590     h = new TH2I("hMcTrkltPhi", "Tracklet Resolution (#Phi)", 48, -.48, .48, 100, -.15, .15);
1591     h->GetXaxis()->SetTitle("tg(#phi)");
1592     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1593     h->GetZaxis()->SetTitle("entries");
1594   } else h->Reset();
1595   arr->AddAt(h, 4);
1596
1597
1598   // KALMAN TRACK RESOLUTION
1599   fContainer->AddAt(arr = new TObjArray(fgNhistos[kMCtrackTRD]), kMCtrackTRD);
1600   arr->SetName("McTrkTRD");
1601   // Kalman track y resolution
1602   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkY"))){
1603     h = new TH2I("hMcTrkY", "Track Y Resolution", 48, -.48, .48, 100, -.2, .2);
1604     h->GetXaxis()->SetTitle("tg(#phi)");
1605     h->GetYaxis()->SetTitle("#Delta y [cm]");
1606     h->GetZaxis()->SetTitle("entries");
1607   } else h->Reset();
1608   arr->AddAt(h, 0);
1609   // Kalman track y pulls
1610   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkYPull"))){
1611     h = new TH2I("hMcTrkYPull", "Track Y Pulls", 48, -.48, .48, 100, -4., 4.);
1612     h->GetXaxis()->SetTitle("tg(#phi)");
1613     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1614     h->GetZaxis()->SetTitle("entries");
1615   } else h->Reset();
1616   arr->AddAt(h, 1);
1617   // Kalman track Z
1618   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkZ"))){
1619     h = new TH2I("hMcTrkZ", "Track Z Resolution", 100, -1., 1., 100, -1., 1.);
1620     h->GetXaxis()->SetTitle("tg(#theta)");
1621     h->GetYaxis()->SetTitle("#Delta z [cm]");
1622     h->GetZaxis()->SetTitle("entries");
1623   } else h->Reset();
1624   arr->AddAt(h, 2);
1625   // Kalman track Z pulls
1626   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkZPull"))){
1627     h = new TH2I("hMcTrkZPull", "Track Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1628     h->GetXaxis()->SetTitle("tg(#theta)");
1629     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1630     h->GetZaxis()->SetTitle("entries");
1631   } else h->Reset();
1632   arr->AddAt(h, 3);
1633   // Kalman track SNP
1634   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNP"))){
1635     h = new TH2I("hMcTrkSNP", "Track Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3);
1636     h->GetXaxis()->SetTitle("tg(#phi)");
1637     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1638     h->GetZaxis()->SetTitle("entries");
1639   } else h->Reset();
1640   arr->AddAt(h, 4);
1641   // Kalman track SNP pulls
1642   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNPPull"))){
1643     h = new TH2I("hMcTrkSNPPull", "Track SNP Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1644     h->GetXaxis()->SetTitle("tg(#phi)");
1645     h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
1646     h->GetZaxis()->SetTitle("entries");
1647   } else h->Reset();
1648   arr->AddAt(h, 5);
1649   // Kalman track TGL
1650   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGL"))){
1651     h = new TH2I("hMcTrkTGL", "Track Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3);
1652     h->GetXaxis()->SetTitle("tg(#theta)");
1653     h->GetYaxis()->SetTitle("#Delta#theta [rad]");
1654     h->GetZaxis()->SetTitle("entries");
1655   } else h->Reset();
1656   arr->AddAt(h, 6);
1657   // Kalman track TGL pulls
1658   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGLPull"))){
1659     h = new TH2I("hMcTrkTGLPull", "Track TGL  Pulls", 100, -1., 1., 100, -4.5, 4.5);
1660     h->GetXaxis()->SetTitle("tg(#theta)");
1661     h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
1662     h->GetZaxis()->SetTitle("entries");
1663   } else h->Reset();
1664   arr->AddAt(h, 7);
1665
1666   const Int_t kNdpt(150); 
1667   const Int_t kNspc = 2*AliPID::kSPECIES+1;
1668   Float_t DPt=-.1, Spc=-5.5;
1669   Float_t binsSpc[kNspc+1], binsDPt[kNdpt+1];
1670   for(Int_t i=0; i<kNspc+1; i++,Spc+=1.) binsSpc[i]=Spc;
1671   for(Int_t i=0; i<kNdpt+1; i++,DPt+=2.e-3) binsDPt[i]=DPt;
1672   TObjArray *arr2 = NULL; TH3S* h3=NULL;
1673   // Kalman track Pt resolution
1674   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 8);
1675   arr2->SetName("Pt Resolution");
1676   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1677     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkPt%d", il)))){
1678       h3 = new TH3S(Form("hMcTrkPt%d", il), "Track Pt Resolution;p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1679     } else h3->Reset();
1680     arr2->AddAt(h3, il);
1681   }
1682   // Kalman track Pt pulls
1683   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 9);
1684   arr2->SetName("1/Pt Pulls");
1685   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1686     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkPtPulls%d", il)))){
1687       h3 = new TH3S(Form("hMcTrkPtPulls%d", il), 
1688       "Track 1/Pt Pulls;1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", 
1689       kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1690     } else h3->Reset();
1691     arr2->AddAt(h3, il);
1692   }
1693   // Kalman track P resolution
1694   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 10);
1695   arr2->SetName("P Resolution");
1696   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1697     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkP%d", il)))){
1698       h3 = new TH3S(Form("hMcTrkP%d", il), "Track P Resolution;p [GeV/c];#Delta p/p^{MC};SPECIES", kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1699     } else h3->Reset();
1700     arr2->AddAt(h3, il);
1701   }
1702   // TRD stand-alone track Pt resolution
1703   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 11);
1704   arr2->SetName("Pt Resolution [SA]");
1705   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1706     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcSATrkPt%d", il)))){
1707       h3 = new TH3S(Form("hMcSATrkPt%d", il), 
1708       "Track Pt Resolution;p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", 
1709       kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1710     } else h3->Reset();
1711     arr2->AddAt(h3, il);
1712   }
1713   // TRD stand-alone track Pt pulls
1714   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 12);
1715   arr2->SetName("1/Pt Pulls [SA]");
1716   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1717     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcSATrkPtPulls%d", il)))){
1718       h3 = new TH3S(Form("hMcSATrkPtPulls%d", il), 
1719       "Track 1/Pt Pulls;1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", 
1720       kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1721     } else h3->Reset();
1722     arr2->AddAt(h3, il);
1723   }
1724   // TRD stand-alone track P resolution
1725   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 13);
1726   arr2->SetName("P Resolution [SA]");
1727   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1728     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcSATrkP%d", il)))){
1729       h3 = new TH3S(Form("hMcSATrkP%d", il), 
1730       "Track P Resolution;p [GeV/c];#Delta p/p^{MC};SPECIES", 
1731       kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1732     } else h3->Reset();
1733     arr2->AddAt(h3, il);
1734   }
1735
1736   // TPC TRACK RESOLUTION
1737   fContainer->AddAt(arr = new TObjArray(fgNhistos[kMCtrackTPC]), kMCtrackTPC);
1738   arr->SetName("McTrkTPC");
1739   // Kalman track Y
1740   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCY"))){
1741     h = new TH2I("hMcTrkTPCY", "Track[TPC] Y Resolution", 60, -.3, .3, 100, -.5, .5);
1742     h->GetXaxis()->SetTitle("tg(#phi)");
1743     h->GetYaxis()->SetTitle("#Delta y [cm]");
1744     h->GetZaxis()->SetTitle("entries");
1745   } else h->Reset();
1746   arr->AddAt(h, 0);
1747   // Kalman track Y pulls
1748   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCYPull"))){
1749     h = new TH2I("hMcTrkTPCYPull", "Track[TPC] Y Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1750     h->GetXaxis()->SetTitle("tg(#phi)");
1751     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1752     h->GetZaxis()->SetTitle("entries");
1753   } else h->Reset();
1754   arr->AddAt(h, 1);
1755   // Kalman track Z
1756   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCZ"))){
1757     h = new TH2I("hMcTrkTPCZ", "Track[TPC] Z Resolution", 100, -1., 1., 100, -1., 1.);
1758     h->GetXaxis()->SetTitle("tg(#theta)");
1759     h->GetYaxis()->SetTitle("#Delta z [cm]");
1760     h->GetZaxis()->SetTitle("entries");
1761   } else h->Reset();
1762   arr->AddAt(h, 2);
1763   // Kalman track Z pulls
1764   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCZPull"))){
1765     h = new TH2I("hMcTrkTPCZPull", "Track[TPC] Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1766     h->GetXaxis()->SetTitle("tg(#theta)");
1767     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1768     h->GetZaxis()->SetTitle("entries");
1769   } else h->Reset();
1770   arr->AddAt(h, 3);
1771   // Kalman track SNP
1772   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCSNP"))){
1773     h = new TH2I("hMcTrkTPCSNP", "Track[TPC] Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3);
1774     h->GetXaxis()->SetTitle("tg(#phi)");
1775     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1776     h->GetZaxis()->SetTitle("entries");
1777   } else h->Reset();
1778   arr->AddAt(h, 4);
1779   // Kalman track SNP pulls
1780   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCSNPPull"))){
1781     h = new TH2I("hMcTrkTPCSNPPull", "Track[TPC] SNP Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1782     h->GetXaxis()->SetTitle("tg(#phi)");
1783     h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
1784     h->GetZaxis()->SetTitle("entries");
1785   } else h->Reset();
1786   arr->AddAt(h, 5);
1787   // Kalman track TGL
1788   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCTGL"))){
1789     h = new TH2I("hMcTrkTPCTGL", "Track[TPC] Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3);
1790     h->GetXaxis()->SetTitle("tg(#theta)");
1791     h->GetYaxis()->SetTitle("#Delta#theta [rad]");
1792     h->GetZaxis()->SetTitle("entries");
1793   } else h->Reset();
1794   arr->AddAt(h, 6);
1795   // Kalman track TGL pulls
1796   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCTGLPull"))){
1797     h = new TH2I("hMcTrkTPCTGLPull", "Track[TPC] TGL  Pulls", 100, -1., 1., 100, -4.5, 4.5);
1798     h->GetXaxis()->SetTitle("tg(#theta)");
1799     h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
1800     h->GetZaxis()->SetTitle("entries");
1801   } else h->Reset();
1802   arr->AddAt(h, 7);
1803   // Kalman track Pt resolution
1804   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPt"))){
1805     h3 = new TH3S("hMcTrkTPCPt", 
1806     "TRDin Pt Resolution;p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", 
1807     kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1808   } else h3->Reset();
1809   arr->AddAt(h3, 8);
1810   // Kalman track Pt pulls
1811   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPtPulls"))){
1812     h3 = new TH3S("hMcTrkTPCPtPulls", 
1813     "Track[TPC] 1/Pt Pulls;1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", 
1814     kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1815   } else h3->Reset();
1816   arr->AddAt(h3, 9);
1817   // Kalman track P resolution
1818   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCP"))){
1819     h3 = new TH3S("hMcTrkTPCP", 
1820     "TRDin P Resolution;p [GeV/c];#Delta p/p^{MC};SPECIES", 
1821     kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1822   } else h3->Reset();
1823   arr->AddAt(h3, 10);
1824   // Kalman track P pulls
1825   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPPulls"))){
1826     h3 = new TH3S("hMcTrkTPCPPulls", 
1827     "TRDin P Pulls;p^{MC} [GeV/c];#Deltap/#sigma_{p};SPECIES", 
1828     kNpt, 0., 12., 100, -5., 5., kNspc, -5.5, 5.5);
1829   } else h3->Reset();
1830   arr->AddAt(h3, 11);
1831
1832
1833
1834   // Kalman track Z resolution [TOF]
1835   fContainer->AddAt(arr = new TObjArray(fgNhistos[kMCtrackTOF]), kMCtrackTOF);
1836   arr->SetName("McTrkTOF");
1837   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTOFZ"))){
1838     h = new TH2I("hMcTrkTOFZ", "Track[TOF] Z Resolution", 100, -1., 1., 100, -1., 1.);
1839     h->GetXaxis()->SetTitle("tg(#theta)");
1840     h->GetYaxis()->SetTitle("#Delta z [cm]");
1841     h->GetZaxis()->SetTitle("entries");
1842   } else h->Reset();
1843   arr->AddAt(h, 0);
1844   // Kalman track Z pulls
1845   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTOFZPull"))){
1846     h = new TH2I("hMcTrkTOFZPull", "Track[TOF] Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1847     h->GetXaxis()->SetTitle("tg(#theta)");
1848     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1849     h->GetZaxis()->SetTitle("entries");
1850   } else h->Reset();
1851   arr->AddAt(h, 1);
1852
1853   return fContainer;
1854 }
1855
1856 //________________________________________________________
1857 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
1858 {
1859   //
1860   // Do the processing
1861   //
1862
1863   Char_t pn[10]; sprintf(pn, "p%03d", fIdxPlot);
1864   Int_t n = 0;
1865   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1866   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1867
1868   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
1869     Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
1870     TH1D *h = h2->ProjectionY(pn, ibin, ibin);
1871     if(h->GetEntries()<100) continue;
1872     //AdjustF1(h, f);
1873
1874     h->Fit(f, "QN");
1875     
1876     Int_t ip = g[0]->GetN();
1877     g[0]->SetPoint(ip, x, k*f->GetParameter(1));
1878     g[0]->SetPointError(ip, 0., k*f->GetParError(1));
1879     g[1]->SetPoint(ip, x, k*f->GetParameter(2));
1880     g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1881     
1882 /*  
1883     g[0]->SetPoint(ip, x, k*h->GetMean());
1884     g[0]->SetPointError(ip, 0., k*h->GetMeanError());
1885     g[1]->SetPoint(ip, x, k*h->GetRMS());
1886     g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
1887   }
1888   fIdxPlot++;
1889   return kTRUE;
1890 }
1891
1892 //________________________________________________________
1893 Bool_t AliTRDresolution::Process2D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k, Int_t gidx)
1894 {
1895   //
1896   // Do the processing
1897   //
1898
1899   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1900
1901   // retrive containers
1902   TH2I *h2 = idx<0 ? (TH2I*)(fContainer->At(plot)) : (TH2I*)((TObjArray*)(fContainer->At(plot)))->At(idx);
1903   if(!h2) return kFALSE;
1904
1905   TGraphErrors *g[2];
1906   if(gidx<0) gidx=idx;
1907   if(!(g[0] = gidx<0 ? (TGraphErrors*)fGraphM->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(gidx))) return kFALSE;
1908
1909   if(!(g[1] = gidx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(gidx))) return kFALSE;
1910
1911   return Process(h2, f, k, g);
1912 }
1913
1914 //________________________________________________________
1915 Bool_t AliTRDresolution::Process3D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1916 {
1917   //
1918   // Do the processing
1919   //
1920
1921   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1922
1923   // retrive containers
1924   TH3S *h3 = idx<0 ? (TH3S*)(fContainer->At(plot)) : (TH3S*)((TObjArray*)(fContainer->At(plot)))->At(idx);
1925   if(!h3) return kFALSE;
1926
1927   TObjArray *gm, *gs;
1928   if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
1929   if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
1930   TGraphErrors *g[2];
1931
1932   TAxis *az = h3->GetZaxis();
1933   for(Int_t iz=1; iz<=az->GetNbins(); iz++){
1934     if(!(g[0] = (TGraphErrors*)gm->At(iz-1))) return kFALSE;
1935     if(!(g[1] = (TGraphErrors*)gs->At(iz-1))) return kFALSE;
1936     az->SetRange(iz, iz);
1937     if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
1938   }
1939
1940   return kTRUE;
1941 }
1942
1943 //________________________________________________________
1944 Bool_t AliTRDresolution::Process3Drange(ETRDresolutionPlot plot, Int_t hidx, Int_t gidx, TF1 *f, Float_t k, Int_t zbin0, Int_t zbin1)
1945 {
1946   //
1947   // Do the processing
1948   //
1949
1950   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1951
1952   // retrive containers
1953   TH3S *h3 = hidx<0 ? (TH3S*)(fContainer->At(plot)) : (TH3S*)((TObjArray*)(fContainer->At(plot)))->At(hidx);
1954   if(!h3) return kFALSE;
1955
1956   TGraphErrors *g[2];
1957   if(!(g[0] = (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(gidx))) return kFALSE;
1958   if(!(g[1] = (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(gidx))) return kFALSE;
1959
1960   TAxis *az = h3->GetZaxis();
1961   az->SetRange(zbin0, zbin1);
1962   if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
1963   return kTRUE;
1964 }
1965
1966 //________________________________________________________
1967 Bool_t AliTRDresolution::Process3DL(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1968 {
1969   //
1970   // Do the processing
1971   //
1972
1973   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1974
1975   // retrive containers
1976   TH3S *h3 = (TH3S*)((TObjArray*)fContainer->At(plot))->At(idx);
1977   if(!h3) return kFALSE;
1978
1979   TGraphAsymmErrors *gm; 
1980   TGraphErrors *gs;
1981   if(!(gm = (TGraphAsymmErrors*)((TObjArray*)fGraphM->At(plot))->At(0))) return kFALSE;
1982   if(!(gs = (TGraphErrors*)((TObjArray*)fGraphS->At(plot)))) return kFALSE;
1983
1984   Float_t x, r, mpv, xM, xm;
1985   TAxis *ay = h3->GetYaxis();
1986   for(Int_t iy=1; iy<=h3->GetNbinsY(); iy++){
1987     ay->SetRange(iy, iy);
1988     x = ay->GetBinCenter(iy);
1989     TH2F *h2=(TH2F*)h3->Project3D("zx");
1990     TAxis *ax = h2->GetXaxis();
1991     for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
1992       r = ax->GetBinCenter(ix);
1993       TH1D *h1 = h2->ProjectionY("py", ix, ix);
1994       if(h1->Integral()<50) continue;
1995       h1->Fit(f, "QN");
1996
1997       GetLandauMpvFwhm(f, mpv, xm, xM);
1998       Int_t ip = gm->GetN();
1999       gm->SetPoint(ip, x, k*mpv);
2000       gm->SetPointError(ip, 0., 0., k*xm, k*xM);
2001       gs->SetPoint(ip, r, k*(xM-xm)/mpv);
2002       gs->SetPointError(ip, 0., 0.);
2003     }
2004   }
2005
2006   return kTRUE;
2007 }
2008
2009 //________________________________________________________
2010 Bool_t AliTRDresolution::Process4D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k, Int_t n)
2011 {
2012   //
2013   // Do the processing
2014   //
2015
2016   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2017   //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
2018
2019   // retrive containers
2020   TObjArray *arr = (TObjArray*)((TObjArray*)(fContainer->At(plot)))->At(idx);
2021   if(!arr) return kFALSE;
2022
2023   TObjArray *gm, *gs;
2024   if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2025   if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2026
2027   TGraphErrors *g[2];
2028
2029   TH3S *h3(NULL);
2030   for(Int_t ix=0, in=0; ix<arr->GetEntriesFast(); ix++){
2031     if(!(h3 = (TH3S*)arr->At(ix))) return kFALSE;
2032     TAxis *az = h3->GetZaxis();
2033     //printf("  process ix[%d] bins[%d] in[%d]\n", ix, az->GetNbins(), in);
2034     for(Int_t iz=1; iz<=az->GetNbins(); iz++, in++){
2035       if(n>=0 && n!=in) continue;
2036       if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2037       if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2038       //printf("    g0[%s] g1[%s]\n", g[0]->GetName(), g[1]->GetName());
2039       az->SetRange(iz, iz);
2040       if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2041     }
2042   }
2043
2044   return kTRUE;
2045 }
2046
2047 //________________________________________________________
2048 Bool_t AliTRDresolution::GetGraphPlot(Float_t *bb, ETRDresolutionPlot ip, Int_t idx)
2049 {
2050   //
2051   // Get the graphs
2052   //
2053
2054   if(!fGraphS || !fGraphM) return kFALSE;
2055
2056   //printf("plotting task[%d] gidx[%d]\n", ip, idx);
2057   TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphM->At(ip)))->At(idx);
2058   if(!gm) return kFALSE;
2059   TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphS->At(ip)))->At(idx);
2060   if(!gs) return kFALSE;
2061   //printf("gs[%s] gm[%s]\n", gs->GetName(), gm->GetName());
2062   gs->Draw("apl"); gm->Draw("pl");
2063   //return kTRUE;
2064   // titles look up
2065   Int_t nref = 0;
2066   for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fgNproj[jp];
2067   UChar_t jdx = idx<0?0:idx;
2068   for(Int_t jc=0; jc<TMath::Min(jdx,fgNproj[ip]-1); jc++) nref++;
2069   const Char_t **at = fgAxTitle[nref];
2070
2071   Int_t n(0);
2072   if((n=gm->GetN())) {
2073     PutTrendValue(Form("%s_%s", fgPerformanceName[ip], at[0]), gm->GetMean(2));
2074     PutTrendValue(Form("%s_%sRMS", fgPerformanceName[ip], at[0]), gm->GetRMS(2));
2075   }
2076
2077   if((n=gs->GetN())){
2078     gs->Sort(&TGraph::CompareY);
2079     PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[ip], at[0]), gs->GetY()[0]);
2080     PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[ip], at[0]), gs->GetY()[n-1]);
2081     gs->Sort(&TGraph::CompareX); 
2082   }
2083
2084   // axis range
2085   TAxis *ax(NULL); TH1 *hf(NULL);
2086   hf = gs->GetHistogram();
2087   hf->SetTitle(at[0]);
2088   ax = hf->GetXaxis();
2089   ax->SetRangeUser(bb[0], bb[2]);
2090   ax->SetTitle(at[1]);ax->CenterTitle();
2091
2092   ax = hf->GetYaxis();
2093   ax->SetRangeUser(bb[1], bb[3]);
2094   ax->SetTitleOffset(1.1);
2095   ax->SetTitle(at[2]);ax->CenterTitle();
2096
2097   TGaxis *gax = NULL;
2098   gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
2099   gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
2100   //gax->SetVertical();
2101   gax->CenterTitle(); gax->SetTitleOffset(.7);
2102   gax->SetTitle(at[3]); gax->Draw();
2103
2104   // bounding box
2105   TBox *b = new TBox(-.15, bb[1], .15, bb[3]);
2106   b->SetFillStyle(3002);b->SetLineColor(0);
2107   b->SetFillColor(ip<=Int_t(kMCcluster)?kGreen:kBlue);
2108   b->Draw();
2109   return kTRUE;
2110 }
2111
2112 //________________________________________________________
2113 Bool_t AliTRDresolution::GetGraphTrack(Float_t *bb, Int_t idx, Int_t first, Int_t n, Bool_t kLEG)
2114 {
2115   //
2116   // Get the graphs
2117   //
2118
2119   if(!fGraphS || !fGraphM) return kFALSE;
2120
2121   // axis titles look up
2122   Int_t nref(0);
2123   for(Int_t jp=0; jp<Int_t(kMCtrackTRD); jp++) nref+=fgNproj[jp];
2124   nref+=idx;
2125   const Char_t **at = fgAxTitle[nref];
2126   //printf("nref[%d] ax[%s] x[%f %f] y[%f %f]\n", nref, at[0], bb[0], bb[2], bb[1], bb[3]);
2127
2128   TLegend *legM(NULL), *legS(NULL);
2129   if(kLEG){
2130     legM=new TLegend(.35, .6, .65, .9);
2131     legM->SetHeader("Mean");
2132     legM->SetBorderSize(1);
2133     legM->SetFillColor(0);
2134     legS=new TLegend(.65, .6, .95, .9);
2135     legS->SetHeader("Sigma");
2136     legS->SetBorderSize(1);
2137     legS->SetFillColor(0);
2138   }
2139   Int_t layer(first/11);
2140   TH1S *h1(NULL);
2141   h1 = new TH1S(Form("h1TF_%02d", fIdxFrame++), Form("%s %d;%s;%s", at[0], layer, at[1], at[2]), 2, bb[0], bb[2]);
2142   h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2143   h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2); 
2144   // axis range
2145   TAxis *ax = h1->GetXaxis();
2146   ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
2147   ax = h1->GetYaxis();
2148   ax->SetRangeUser(bb[1], bb[3]);
2149   ax->CenterTitle(); ax->SetTitleOffset(1.4);
2150   h1->Draw();
2151 //   TGaxis *gax = NULL;
2152 //   gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
2153 //   gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
2154 //   //gax->SetVertical();
2155 //   gax->CenterTitle(); //gax->SetTitleOffset(.5);gax->SetTitleSize(.06);
2156 //   gax->SetTitle(at[3]); gax->Draw();
2157
2158   TGraphErrors *gm = NULL, *gs = NULL;
2159   TObjArray *a0 = NULL, *a1 = NULL;
2160   a0 = (TObjArray*)((TObjArray*)fGraphM->At(kMCtrackTRD))->At(idx);
2161   a1 = (TObjArray*)((TObjArray*)fGraphS->At(kMCtrackTRD))->At(idx);
2162   Int_t nn(0), m(n/2);
2163   for(Int_t is=first, is0=0; is<first+n; is++, is0++){
2164     if(is0==m) continue; // no PID tracks
2165     if(is0==m-1||is0==m+1) continue; // electron tracks
2166     if(!(gs =  (TGraphErrors*)a1->At(is))) return kFALSE;
2167     if(!(gm =  (TGraphErrors*)a0->At(is))) return kFALSE;
2168
2169     if((nn=gs->GetN())){
2170       gs->Draw("pc"); if(legS) legS->AddEntry(gs, gs->GetTitle(), "pl");
2171       gs->Sort(&TGraph::CompareY);
2172       PutTrendValue(Form("%s_%sSigMin%s", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is0)), gs->GetY()[0]);
2173       PutTrendValue(Form("%s_%sSigMax%s", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is0)), gs->GetY()[nn-1]);
2174       gs->Sort(&TGraph::CompareX); 
2175     }
2176     if(gm->GetN()){
2177       gm->Draw("pc");if(legM) legM->AddEntry(gm, gm->GetTitle(), "pl");
2178       PutTrendValue(Form("%s_%s_%s", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is0)), gm->GetMean(2));
2179       PutTrendValue(Form("%s_%s_%sRMS", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is0)), gm->GetRMS(2));
2180     }
2181   }
2182   if(kLEG){legM->Draw();legS->Draw();}
2183   return kTRUE;
2184 }
2185
2186
2187 //________________________________________________________
2188 Bool_t AliTRDresolution::GetGraphTrackTPC(Float_t *bb, Int_t idx, Int_t ist, Int_t n, Bool_t kLEG)
2189 {
2190   //
2191   // Get the graphs
2192   //
2193
2194   if(!fGraphS || !fGraphM) return kFALSE;
2195
2196   // axis titles look up
2197   Int_t nref = 0;
2198   for(Int_t jp=0; jp<Int_t(kMCtrackTPC); jp++) nref+=fgNproj[jp];
2199   nref+=idx;
2200   const Char_t **at = fgAxTitle[nref];
2201
2202   TLegend *legM(NULL), *legS(NULL);
2203   if(kLEG){
2204     legM=new TLegend(.35, .6, .65, .9);
2205     legM->SetHeader("Mean");
2206     legM->SetBorderSize(1);
2207     legM->SetFillColor(0);
2208     legS=new TLegend(.65, .6, .95, .9);
2209     legS->SetHeader("Sigma");
2210     legS->SetBorderSize(1);
2211     legS->SetFillColor(0);
2212   }
2213   TH1S *h1(NULL);
2214   h1 = new TH1S(Form("h1TF_%02d", fIdxFrame++), at[0], 2, bb[0], bb[2]);
2215   h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2216   h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2); 
2217   // axis range
2218   TAxis *ax = h1->GetXaxis();
2219   ax->SetTitle(at[1]);ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
2220   ax = h1->GetYaxis();
2221   ax->SetRangeUser(bb[1], bb[3]);
2222   ax->SetTitleOffset(1.4);//ax->SetTitleSize(.06);
2223   ax->SetTitle(at[2]);ax->CenterTitle();
2224   h1->Draw();
2225 //   TGaxis *gax = NULL;
2226 //   gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
2227 //   gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
2228 //   //gax->SetVertical();
2229 //   gax->CenterTitle(); //gax->SetTitleOffset(.5);gax->SetTitleSize(.06);
2230 //   gax->SetTitle(at[3]); gax->Draw();
2231
2232   Int_t nn(0), m(n/2);
2233   TGraphErrors *gm = NULL, *gs = NULL;
2234   TObjArray *a0 = NULL, *a1 = NULL;
2235   a0 = (TObjArray*)((TObjArray*)fGraphM->At(kMCtrackTPC))->At(idx);
2236   a1 = (TObjArray*)((TObjArray*)fGraphS->At(kMCtrackTPC))->At(idx);
2237   for(Int_t is=ist, is0=0; is<ist+n; is++, is0++){
2238     if(is0==m) continue; // no PID tracks
2239     //if(is0==m-1||is0==m+1) continue; // electron tracks
2240     if(!(gs =  (TGraphErrors*)a1->At(is))) return kFALSE;
2241     if(!(gm =  (TGraphErrors*)a0->At(is))) return kFALSE;
2242     if((nn=gs->GetN())){
2243       gs->Draw("pl");if(legS) legS->AddEntry(gs, gs->GetTitle(), "pl");
2244       gs->Sort(&TGraph::CompareY);
2245       PutTrendValue(Form("%s_%sSigMin%s", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is0)), gs->GetY()[0]);
2246       PutTrendValue(Form("%s_%sSigMax%s", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is0)), gs->GetY()[nn-1]);
2247       gs->Sort(&TGraph::CompareX); 
2248     }
2249     if(gm->GetN()){
2250       gm->Draw("pl"); if(legM) legM->AddEntry(gm, gm->GetTitle(), "pl");
2251       PutTrendValue(Form("%s_%s_%s", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is0)), gm->GetMean(2));
2252       PutTrendValue(Form("%s_%s_%sRMS", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is0)), gm->GetRMS(2));
2253     }
2254   }
2255   if(kLEG) {legM->Draw(); legS->Draw();}
2256   return kTRUE;
2257 }
2258
2259 //________________________________________________________
2260 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
2261 {
2262   //
2263   // Get the most probable value and the full width half mean 
2264   // of a Landau distribution
2265   //
2266
2267   const Float_t dx = 1.;
2268   mpv = f->GetParameter(1);
2269   Float_t fx, max = f->Eval(mpv);
2270
2271   xm = mpv - dx;
2272   while((fx = f->Eval(xm))>.5*max){
2273     if(fx>max){ 
2274       max = fx;
2275       mpv = xm;
2276     }
2277     xm -= dx;
2278   }
2279
2280   xM += 2*(mpv - xm);
2281   while((fx = f->Eval(xM))>.5*max) xM += dx;
2282 }
2283
2284
2285 //________________________________________________________
2286 void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
2287 {
2288
2289   fReconstructor->SetRecoParam(r);
2290 }