be7e2924f5ded399f54995dacf8f49c3ee89d822
[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDresolution.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercialf purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15
16 /* $Id: AliTRDresolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
17
18 ////////////////////////////////////////////////////////////////////////////
19 //                                                                        //
20 //  TRD tracking resolution                                               //
21 //
22 // The class performs resolution and residual studies 
23 // of the TRD tracks for the following quantities :
24 //   - spatial position (y, [z])
25 //   - angular (phi) tracklet
26 //   - momentum at the track level
27 // 
28 // The class has to be used for regular detector performance checks using the official macros:
29 //   - $ALICE_ROOT/TRD/qaRec/run.C
30 //   - $ALICE_ROOT/TRD/qaRec/makeResults.C
31 // 
32 // For stand alone usage please refer to the following example: 
33 // {  
34 //   gSystem->Load("libANALYSIS.so");
35 //   gSystem->Load("libTRDqaRec.so");
36 //   AliTRDresolution *res = new AliTRDresolution();
37 //   //res->SetMCdata();
38 //   //res->SetVerbose();
39 //   //res->SetVisual();
40 //   res->Load("TRD.TaskResolution.root");
41 //   if(!res->PostProcess()) return;
42 //   res->GetRefFigure(0);
43 // }  
44 //
45 //  Authors:                                                              //
46 //    Alexandru Bercuci <A.Bercuci@gsi.de>                                //
47 //    Markus Fasel <M.Fasel@gsi.de>                                       //
48 //                                                                        //
49 ////////////////////////////////////////////////////////////////////////////
50
51 #include <cstring>
52
53 #include <TROOT.h>
54 #include <TSystem.h>
55 #include <TPDGCode.h>
56 #include <TObjArray.h>
57 #include <TH3.h>
58 #include <TH2.h>
59 #include <TH1.h>
60 #include <TF1.h>
61 #include <TCanvas.h>
62 #include <TGaxis.h>
63 #include <TBox.h>
64 #include <TProfile.h>
65 #include <TGraphErrors.h>
66 #include <TGraphAsymmErrors.h>
67 #include <TMath.h>
68 #include <TMatrixT.h>
69 #include <TVectorT.h>
70 #include "TTreeStream.h"
71 #include "TGeoManager.h"
72
73 #include "AliAnalysisManager.h"
74 #include "AliTrackReference.h"
75 #include "AliTrackPointArray.h"
76 #include "AliCDBManager.h"
77 #include "AliPID.h"
78
79 #include "AliTRDcalibDB.h"
80 #include "AliTRDCommonParam.h"
81 #include "AliTRDSimParam.h"
82 #include "AliTRDgeometry.h"
83 #include "AliTRDpadPlane.h"
84 #include "AliTRDcluster.h"
85 #include "AliTRDseedV1.h"
86 #include "AliTRDtrackV1.h"
87 #include "AliTRDtrackerV1.h"
88 #include "AliTRDReconstructor.h"
89 #include "AliTRDrecoParam.h"
90
91 #include "info/AliTRDclusterInfo.h"
92 #include "info/AliTRDtrackInfo.h"
93 #include "AliTRDresolution.h"
94
95 ClassImp(AliTRDresolution)
96 UChar_t AliTRDresolution::fNElements[kNhistos] = {
97   2, 2, 5, 5,
98   2, 5, 12, 2, 11
99 };
100 Char_t* AliTRDresolution::fPerformanceName[kNhistos] = {
101      "Charge"
102     ,"Cluster2Track"
103     ,"Tracklet2Track"
104     ,"Tracklet2TPC" 
105     ,"Cluster2MC"
106     ,"Tracklet2MC"
107     ,"TPC2MC"
108     ,"TOF/HMPID2MC"
109     ,"TRD2MC"
110 };
111 //Int_t AliTRDresolution::fNPlots = 0; 
112 Char_t *AliTRDresolution::fAxTitle[46][3] = {
113   // Charge
114   {"x [cm]", "I_{mpv}", "x/x_{0}"}
115  ,{"x/x_{0}", "#delta I/I_{mpv}", "x[cm]"}
116   // Clusters to Kalman
117  ,{"tg(#phi)", "#mu_{y}^{cl} [#mum]", "#sigma_{y}^{cl} [#mum]"}
118  ,{"tg(#phi)", "PULL: #mu_{y}^{cl}", "PULL: #sigma_{y}^{cl}"}
119   // TRD tracklet to Kalman fit
120  ,{"tg(#phi)", "#mu_{y}^{trklt} [#mum]", "#sigma_{y}^{trklt} [#mum]"}
121  ,{"tg(#phi)", "PULL: #mu_{y}^{trklt}", "PULL: #sigma_{y}^{trklt}"}
122  ,{"tg(#theta)", "#mu_{z}^{trklt} [#mum]", "#sigma_{z}^{trklt} [#mum]"}
123  ,{"tg(#theta)", "PULL: #mu_{z}^{trklt}", "PULL: #sigma_{z}^{trklt}"}
124  ,{"tg(#phi)", "#mu_{#phi}^{trklt} [mrad]", "#sigma_{#phi}^{trklt} [mrad]"}
125   // TPC track 2 first TRD tracklet
126  ,{"tg(#phi)", "#mu_{y}^{TPC trklt} [#mum]", "#sigma_{y}^{TPC trklt} [#mum]"}
127  ,{"tg(#phi)", "PULL: #mu_{y}^{TPC trklt}", "PULL: #sigma_{y}^{TPC trklt}"}
128  ,{"tg(#theta)", "#mu_{z}^{TPC trklt} [#mum]", "#sigma_{z}^{TPC trklt} [#mum]"}
129  ,{"tg(#theta)", "PULL: #mu_{z}^{TPC trklt}", "PULL: #sigma_{z}^{TPC trklt}"}
130  ,{"tg(#phi)", "#mu_{#phi}^{TPC trklt} [mrad]", "#sigma_{#phi}^{TPC trklt} [mrad]"}
131   // MC cluster
132  ,{"tg(#phi)", "MC: #mu_{y}^{cl} [#mum]", "MC: #sigma_{y}^{cl} [#mum]"}
133  ,{"tg(#phi)", "MC PULL: #mu_{y}^{cl}", "MC PULL: #sigma_{y}^{cl}"}
134   // MC tracklet
135  ,{"tg(#phi)", "MC: #mu_{y}^{trklt} [#mum]", "MC: #sigma_{y}^{trklt} [#mum]"}
136  ,{"tg(#phi)", "MC PULL: #mu_{y}^{trklt}", "MC PULL: #sigma_{y}^{trklt}"}
137  ,{"tg(#theta)", "MC: #mu_{z}^{trklt} [#mum]", "MC: #sigma_{z}^{trklt} [#mum]"}
138  ,{"tg(#theta)", "MC PULL: #mu_{z}^{trklt}", "MC PULL: #sigma_{z}^{trklt}"}
139  ,{"tg(#phi)", "MC: #mu_{#phi}^{trklt} [mrad]", "MC: #sigma_{#phi}^{trklt} [mrad]"}
140   // MC track TPC
141  ,{"tg(#phi)", "MC: #mu_{y}^{TPC} [#mum]", "MC: #sigma_{y}^{TPC} [#mum]"}
142  ,{"tg(#phi)", "MC PULL: #mu_{y}^{TPC}", "MC PULL: #sigma_{y}^{TPC}"}
143  ,{"tg(#theta)", "MC: #mu_{z}^{TPC} [#mum]", "MC: #sigma_{z}^{TPC} [#mum]"}
144  ,{"tg(#theta)", "MC PULL: #mu_{z}^{TPC}", "MC PULL: #sigma_{z}^{TPC}"}
145  ,{"tg(#phi)", "MC: #mu_{#phi}^{TPC} [mrad]", "MC: #sigma_{#phi}^{TPC} [mrad]"}
146  ,{"tg(#phi)", "MC PULL: #mu_{snp}^{TPC}", "MC PULL: #sigma_{snp}^{TPC}"}
147  ,{"tg(#theta)", "MC: #mu_{#theta}^{TPC} [mrad]", "MC: #sigma_{#theta}^{TPC} [mrad]"}
148  ,{"tg(#theta)", "MC PULL: #mu_{tgl}^{TPC}", "MC PULL: #sigma_{tgl}^{TPC}"}
149  ,{"p_{t}^{MC} [GeV/c]", "MC: #mu^{TPC}(#Deltap_{t}/p_{t}^{MC}) [%]", "MC: #sigma^{TPC}(#Deltap_{t}/p_{t}^{MC}) [%]"}
150  ,{"1/p_{t}^{MC} [c/GeV]", "MC PULL: #mu_{1/p_{t}}^{TPC}", "MC PULL: #sigma_{1/p_{t}}^{TPC}"}
151  ,{"p^{MC} [GeV/c]", "MC: #mu^{TPC}(#Deltap/p^{MC}) [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
152  ,{"p^{MC} [GeV/c]", "MC PULL: #mu^{TPC}(#Deltap/#sigma_{p})", "MC PULL: #sigma^{TPC}(#Deltap/#sigma_{p})"}
153   // MC track TOF
154  ,{"tg(#theta)", "MC: #mu_{z}^{TOF} [#mum]", "MC: #sigma_{z}^{TOF} [#mum]"}
155  ,{"tg(#theta)", "MC PULL: #mu_{z}^{TOF}", "MC PULL: #sigma_{z}^{TOF}"}
156   // MC track in TRD
157  ,{"tg(#phi)", "MC: #mu_{y}^{Trk} [#mum]", "MC: #sigma_{y}^{Trk} [#mum]"}
158  ,{"tg(#phi)", "MC PULL: #mu_{y}^{Trk}", "MC PULL: #sigma_{y}^{Trk}"}
159  ,{"tg(#theta)", "MC: #mu_{z}^{Trk} [#mum]", "MC: #sigma_{z}^{Trk} [#mum]"}
160  ,{"tg(#theta)", "MC PULL: #mu_{z}^{Trk}", "MC PULL: #sigma_{z}^{Trk}"}
161  ,{"tg(#phi)", "MC: #mu_{#phi}^{Trk} [mrad]", "MC: #sigma_{#phi}^{Trk} [mrad]"}
162  ,{"tg(#phi)", "MC PULL: #mu_{snp}^{Trk}", "MC PULL: #sigma_{snp}^{Trk}"}
163  ,{"tg(#theta)", "MC: #mu_{#theta}^{Trk} [mrad]", "MC: #sigma_{#theta}^{Trk} [mrad]"}
164  ,{"tg(#theta)", "MC PULL: #mu_{tgl}^{Trk}", "MC PULL: #sigma_{tgl}^{Trk}"}
165  ,{"p_{t}^{MC} [GeV/c]", "MC: #mu^{Trk}(#Deltap_{t}/p_{t}^{MC}) [%]", "MC: #sigma^{Trk}(#Deltap_{t}/p_{t}^{MC}) [%]"}
166  ,{"1/p_{t}^{MC} [c/GeV]", "MC PULL: #mu_{1/p_{t}}^{Trk}", "MC PULL: #sigma_{1/p_{t}}^{Trk}"}
167  ,{"p^{MC} [GeV/c]", "MC: #mu^{Trk}(#Deltap/p^{MC}) [%]", "MC: #sigma^{Trk}(#Deltap/p^{MC}) [%]"}
168 };
169
170 //________________________________________________________
171 AliTRDresolution::AliTRDresolution()
172   :AliTRDrecoTask("resolution", "Spatial and momentum TRD resolution checker")
173   ,fStatus(0)
174   ,fIdxPlot(0)
175   ,fReconstructor(0x0)
176   ,fGeo(0x0)
177   ,fGraphS(0x0)
178   ,fGraphM(0x0)
179   ,fCl(0x0)
180   ,fTrklt(0x0)
181   ,fMCcl(0x0)
182   ,fMCtrklt(0x0)
183 {
184   fReconstructor = new AliTRDReconstructor();
185   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
186   fGeo = new AliTRDgeometry();
187
188   InitFunctorList();
189
190   DefineOutput(1, TObjArray::Class()); // cluster2track
191   DefineOutput(2, TObjArray::Class()); // tracklet2track
192   DefineOutput(3, TObjArray::Class()); // cluster2mc
193   DefineOutput(4, TObjArray::Class()); // tracklet2mc
194 }
195
196 //________________________________________________________
197 AliTRDresolution::~AliTRDresolution()
198 {
199   if(fGraphS){fGraphS->Delete(); delete fGraphS;}
200   if(fGraphM){fGraphM->Delete(); delete fGraphM;}
201   delete fGeo;
202   delete fReconstructor;
203   if(gGeoManager) delete gGeoManager;
204   if(fCl){fCl->Delete(); delete fCl;}
205   if(fTrklt){fTrklt->Delete(); delete fTrklt;}
206   if(fMCcl){fMCcl->Delete(); delete fMCcl;}
207   if(fMCtrklt){fMCtrklt->Delete(); delete fMCtrklt;}
208 }
209
210
211 //________________________________________________________
212 void AliTRDresolution::CreateOutputObjects()
213 {
214   // spatial resolution
215   OpenFile(0, "RECREATE");
216
217   fContainer = Histos();
218
219   fCl = new TObjArray();
220   fCl->SetOwner(kTRUE);
221   fTrklt = new TObjArray();
222   fTrklt->SetOwner(kTRUE);
223   fMCcl = new TObjArray();
224   fMCcl->SetOwner(kTRUE);
225   fMCtrklt = new TObjArray();
226   fMCtrklt->SetOwner(kTRUE);
227 }
228
229 //________________________________________________________
230 void AliTRDresolution::Exec(Option_t *opt)
231 {
232   fCl->Delete();
233   fTrklt->Delete();
234   fMCcl->Delete();
235   fMCtrklt->Delete();
236
237   AliTRDrecoTask::Exec(opt);
238
239   PostData(1, fCl);
240   PostData(2, fTrklt);
241   PostData(3, fMCcl);
242   PostData(4, fMCtrklt);
243 }
244
245 //________________________________________________________
246 TH1* AliTRDresolution::PlotCharge(const AliTRDtrackV1 *track)
247 {
248   if(track) fTrack = track;
249   if(!fTrack){
250     AliWarning("No Track defined.");
251     return 0x0;
252   }
253   TObjArray *arr = 0x0;
254   if(!(arr = ((TObjArray*)fContainer->At(kCharge)))){
255     AliWarning("No output container defined.");
256     return 0x0;
257   }
258   TH3S* h = 0x0;
259
260   AliTRDseedV1 *fTracklet = 0x0;  
261   AliTRDcluster *c = 0x0;
262   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
263     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
264     if(!fTracklet->IsOK()) continue;
265     Float_t x0 = fTracklet->GetX0();
266     Float_t dq, dl;
267     for(Int_t itb=AliTRDseedV1::kNtb; itb--;){
268       if(!(c = fTracklet->GetClusters(itb))){ 
269         if(!(c = fTracklet->GetClusters(AliTRDseedV1::kNtb+itb))) continue;
270       }
271       dq = fTracklet->GetdQdl(itb, &dl);
272       dl /= 0.15; // dl/dl0, dl0 = 1.5 mm for nominal vd
273       (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dq);
274     }
275
276 //     if(!HasMCdata()) continue;
277 //     UChar_t s;
278 //     Float_t pt0, y0, z0, dydx0, dzdx0;
279 //     if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
280
281   }
282   return h;
283 }
284
285
286 //________________________________________________________
287 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
288 {
289   if(track) fTrack = track;
290   if(!fTrack){
291     AliWarning("No Track defined.");
292     return 0x0;
293   }
294   TObjArray *arr = 0x0;
295   if(!(arr = ((TObjArray*)fContainer->At(kCluster)))){
296     AliWarning("No output container defined.");
297     return 0x0;
298   }
299
300   Double_t cov[7];
301   Float_t x0, y0, z0, dy, dydx, dzdx;
302   AliTRDseedV1 *fTracklet = 0x0;  
303   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
304     if(!(fTracklet = fTrack->GetTracklet(ily))) continue;
305     if(!fTracklet->IsOK()) continue;
306     x0 = fTracklet->GetX0();
307
308     // retrive the track angle with the chamber
309     y0   = fTracklet->GetYref(0);
310     z0   = fTracklet->GetZref(0);
311     dydx = fTracklet->GetYref(1);
312     dzdx = fTracklet->GetZref(1);
313     fTracklet->GetCovRef(cov);
314     Float_t tilt = fTracklet->GetTilt();
315     AliTRDcluster *c = 0x0;
316     fTracklet->ResetClusterIter(kFALSE);
317     while((c = fTracklet->PrevCluster())){
318       Float_t xc = c->GetX();
319       Float_t yc = c->GetY();
320       Float_t zc = c->GetZ();
321       Float_t dx = x0 - xc; 
322       Float_t yt = y0 - dx*dydx;
323       Float_t zt = z0 - dx*dzdx; 
324       yc -= tilt*(zc-zt); // tilt correction
325       dy = yt - yc;
326
327       //Float_t sx2 = dydx*c->GetSX(c->GetLocalTimeBin()); sx2*=sx2;
328       Float_t sy2 = c->GetSigmaY2();
329       if(sy2<=0.) continue;
330       ((TH2I*)arr->At(0))->Fill(dydx, dy);
331       ((TH2I*)arr->At(1))->Fill(dydx, dy/TMath::Sqrt(cov[0] /*+ sx2*/ + sy2));
332   
333       if(fDebugLevel>=1){
334         // Get z-position with respect to anode wire
335         //AliTRDSimParam    *simParam    = AliTRDSimParam::Instance();
336         Int_t istk = fGeo->GetStack(c->GetDetector());
337         AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
338         Float_t row0 = pp->GetRow0();
339         Float_t d  =  row0 - zt + pp->GetAnodeWireOffset();
340         d -= ((Int_t)(2 * d)) / 2.0;
341         if (d > 0.25) d  = 0.5 - d;
342
343 /*        AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
344         fCl->Add(clInfo);
345         clInfo->SetCluster(c);
346         clInfo->SetGlobalPosition(yt, zt, dydx, dzdx);
347         clInfo->SetResolution(dy);
348         clInfo->SetAnisochronity(d);
349         clInfo->SetDriftLength(dx);
350         (*fDebugStream) << "ClusterResiduals"
351           <<"clInfo.=" << clInfo
352           << "\n";*/
353       }
354     }
355   }
356   return (TH2I*)arr->At(0);
357 }
358
359
360 //________________________________________________________
361 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
362 {
363 // Plot normalized residuals for tracklets to track. 
364 // 
365 // We start from the result that if X=N(|m|, |Cov|)
366 // BEGIN_LATEX
367 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
368 // END_LATEX
369 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial 
370 // reference position. 
371   if(track) fTrack = track;
372   if(!fTrack){
373     AliWarning("No Track defined.");
374     return 0x0;
375   }
376   TObjArray *arr = 0x0;
377   if(!(arr = (TObjArray*)fContainer->At(kTrackTRD ))){
378     AliWarning("No output container defined.");
379     return 0x0;
380   }
381
382   Double_t cov[3], covR[7]/*, sqr[3], inv[3]*/;
383   Float_t x, dx, dy, dz;
384   AliTRDseedV1 *fTracklet = 0x0;  
385   for(Int_t il=AliTRDgeometry::kNlayer; il--;){
386     if(!(fTracklet = fTrack->GetTracklet(il))) continue;
387     if(!fTracklet->IsOK()) continue;
388     x    = fTracklet->GetX();
389     dx   = fTracklet->GetX0() - x;
390     // compute dy^2 and dz^2
391     dy   = fTracklet->GetYref(0)-dx*fTracklet->GetYref(1) - fTracklet->GetY();
392     dz   = fTracklet->GetZref(0)-dx*fTracklet->GetZref(1) - fTracklet->GetZ();
393     // compute covariance matrix
394     fTracklet->GetCovAt(x, cov);
395     fTracklet->GetCovRef(covR);
396     cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2]; 
397 /*  // Correct PULL calculation by considering off  
398     // diagonal elements in the covariance matrix
399     // compute square root matrix
400     if(AliTRDseedV1::GetCovInv(cov, inv)==0.) continue;
401     if(AliTRDseedV1::GetCovSqrt(inv, sqr)<0.) continue;
402     Double_t y = sqr[0]*dy+sqr[1]*dz;
403     Double_t z = sqr[1]*dy+sqr[2]*dz;
404     ((TH3*)h)->Fill(y, z, fTracklet->GetYref(1));*/
405
406     ((TH2I*)arr->At(0))->Fill(fTracklet->GetYref(1), dy);
407     ((TH2I*)arr->At(1))->Fill(fTracklet->GetYref(1), dy/TMath::Sqrt(cov[0]));
408     ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), TMath::ATan((fTracklet->GetYref(1)-fTracklet->GetYfit(1))/(1-fTracklet->GetYref(1)*fTracklet->GetYfit(1))));
409     if(!fTracklet->IsRowCross()) continue;
410     ((TH2I*)arr->At(2))->Fill(fTracklet->GetZref(1), dz);
411     ((TH2I*)arr->At(3))->Fill(fTracklet->GetZref(1), dz/TMath::Sqrt(cov[2]));
412   }
413
414
415   return (TH2I*)arr->At(0);
416 }
417
418
419 //________________________________________________________
420 TH1* AliTRDresolution::PlotTrackTPC(const AliTRDtrackV1 *track)
421 {
422 // Store resolution/pulls of Kalman before updating with the TRD information 
423 // at the radial position of the first tracklet. The following points are used 
424 // for comparison  
425 //  - the (y,z,snp) of the first TRD tracklet
426 //  - the (y, z, snp, tgl, pt) of the MC track reference
427 // 
428 // Additionally the momentum resolution/pulls are calculated for usage in the 
429 // PID calculation. 
430
431   if(track) fTrack = track;
432   if(!fTrack){
433     AliWarning("No Track defined.");
434     return 0x0;
435   }
436   AliExternalTrackParam *tin = 0x0;
437   if(!(tin = fTrack->GetTrackLow())){
438     AliWarning("Track did not entered TRD fiducial volume.");
439     return 0x0;
440   }
441   TH1 *h = 0x0;
442   
443   Double_t x = tin->GetX();
444   AliTRDseedV1 *tracklet = 0x0;  
445   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
446     if(!(tracklet = fTrack->GetTracklet(ily))) continue;
447     break;
448   }
449   if(!tracklet || TMath::Abs(x-tracklet->GetX())>1.e-3){
450     AliWarning("Tracklet did not match TRD entrance.");
451     return 0x0;
452   }
453   const Int_t kNPAR(5);
454   Double_t parR[kNPAR]; memcpy(parR, tin->GetParameter(), kNPAR*sizeof(Double_t));
455   Double_t covR[3*kNPAR]; memcpy(covR, tin->GetCovariance(), 3*kNPAR*sizeof(Double_t));
456   Double_t cov[3]; tracklet->GetCovAt(x, cov);
457
458   // define sum covariances
459   TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
460   Double_t *pc = &covR[0], *pp = &parR[0];
461   for(Int_t ir=0; ir<kNPAR; ir++, pp++){
462     PAR(ir) = (*pp);
463     for(Int_t ic = 0; ic<=ir; ic++,pc++){ 
464       COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
465     }
466   }
467   PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
468   //COV.Print(); PAR.Print();
469
470   //TODO Double_t dydx =  TMath::Sqrt(1.-parR[2]*parR[2])/parR[2]; 
471   Double_t dy = parR[0] - tracklet->GetY(); 
472   TObjArray *arr = (TObjArray*)fContainer->At(kTrackTPC);
473   ((TH2I*)arr->At(0))->Fill(tracklet->GetYref(1), dy);
474   ((TH2I*)arr->At(1))->Fill(tracklet->GetYref(1), dy/TMath::Sqrt(COV(0,0)+cov[0]));
475   if(tracklet->IsRowCross()){
476     Double_t dz = parR[1] - tracklet->GetZ(); 
477     ((TH2I*)arr->At(2))->Fill(tracklet->GetZref(1), dz);
478     ((TH2I*)arr->At(3))->Fill(tracklet->GetZref(1), dz/TMath::Sqrt(COV(1,1)+cov[2]));
479   }
480   Double_t dphi = TMath::ASin(PAR[2])-TMath::ATan(tracklet->GetYfit(1));  ((TH2I*)arr->At(4))->Fill(tracklet->GetYref(1), dphi);
481
482
483   // register reference histo for mini-task
484   h = (TH2I*)arr->At(0);
485
486   if(fDebugLevel>=1){
487     (*fDebugStream) << "trackIn"
488       << "x="       << x
489       << "P="       << &PAR
490       << "C="       << &COV
491       << "\n";
492
493     Double_t y = tracklet->GetY(); 
494     Double_t z = tracklet->GetZ(); 
495     (*fDebugStream) << "trackletIn"
496       << "y="       << y
497       << "z="       << z
498       << "Vy="      << cov[0]
499       << "Cyz="     << cov[1]
500       << "Vz="      << cov[2]
501       << "\n";
502   }
503
504
505   if(!HasMCdata()) return h;
506   UChar_t s;
507   Float_t dx, pt0, x0=tracklet->GetX0(), y0, z0, dydx0, dzdx0;
508   if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
509   // translate to reference radial position
510   dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
511   Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi))
512   //Fill MC info
513   TVectorD PARMC(kNPAR);
514   PARMC[0]=y0; PARMC[1]=z0;
515   PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
516   PARMC[4]=1./pt0;
517
518 //   TMatrixDSymEigen eigen(COV);
519 //   TVectorD evals = eigen.GetEigenValues();
520 //   TMatrixDSym evalsm(kNPAR);
521 //   for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
522 //   TMatrixD evecs = eigen.GetEigenVectors();
523 //   TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
524   
525   // fill histos
526   arr = (TObjArray*)fContainer->At(kMCtrackTPC);
527   // y resolution/pulls
528   ((TH2I*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0]);
529   ((TH2I*)arr->At(1))->Fill(dydx0, (PARMC[0]-PAR[0])/TMath::Sqrt(COV(0,0)));
530   // z resolution/pulls
531   ((TH2I*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1]);
532   ((TH2I*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)));
533   // phi resolution/snp pulls
534   ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
535   ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
536   // theta resolution/tgl pulls
537   ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
538   ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
539   // pt resolution\\1/pt pulls\\p resolution/pull
540   for(Int_t is=AliPID::kSPECIES; is--;){
541     if(TMath::Abs(fMC->GetPDG())!=AliPID::ParticleCode(is)) continue;
542     ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., is);
543     ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), is);
544
545     Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
546     Float_t sp;
547     p = tracklet->GetMomentum(&sp);
548     ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., is);
549     ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/sp, is);
550     break;
551   }
552
553   // fill debug for MC 
554   if(fDebugLevel>=1){
555     (*fDebugStream) << "trackInMC"
556       << "P="   << &PARMC
557       << "\n";
558   }
559   return h;
560 }
561
562 //________________________________________________________
563 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
564 {
565   if(!HasMCdata()){ 
566     AliWarning("No MC defined. Results will not be available.");
567     return 0x0;
568   }
569   if(track) fTrack = track;
570   if(!fTrack){
571     AliWarning("No Track defined.");
572     return 0x0;
573   }
574   TObjArray *arr = 0x0;
575   TH1 *h = 0x0;
576   UChar_t s;
577   Int_t pdg = fMC->GetPDG(), det=-1;
578   Int_t label = fMC->GetLabel();
579   Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl;
580   Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
581   Double_t covR[7]/*, cov[3]*/;
582
583   if(fDebugLevel>=1){
584     Double_t DX[12], DY[12], DZ[12], DPt[12], COV[12][15];
585     fMC->PropagateKalman(DX, DY, DZ, DPt, COV);
586     (*fDebugStream) << "MCkalman"
587       << "pdg="  << pdg
588       << "dx0="  << DX[0]
589       << "dx1="  << DX[1]
590       << "dx2="  << DX[2]
591       << "dy0="  << DY[0]
592       << "dy1="  << DY[1]
593       << "dy2="  << DY[2]
594       << "dz0="  << DZ[0]
595       << "dz1="  << DZ[1]
596       << "dz2="  << DZ[2]
597       << "dpt0=" << DPt[0]
598       << "dpt1=" << DPt[1]
599       << "dpt2=" << DPt[2]
600       << "\n";
601   }
602
603   AliTRDReconstructor rec;
604   AliTRDseedV1 *fTracklet = 0x0;  
605   for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
606     if(!(fTracklet = fTrack->GetTracklet(ily)))/* ||
607        !fTracklet->IsOK())*/ continue;
608
609     det = fTracklet->GetDetector();
610     x0  = fTracklet->GetX0();
611     //radial shift with respect to the MC reference (radial position of the pad plane)
612     x= fTracklet->GetX();
613     if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
614     xAnode  = fTracklet->GetX0();
615
616     // MC track position at reference radial position
617     dx  = x0 - x;
618     if(fDebugLevel>=1){
619       (*fDebugStream) << "MC"
620         << "det="     << det
621         << "pdg="     << pdg
622         << "pt="      << pt0
623         << "x="       << x0
624         << "y="       << y0
625         << "z="       << z0
626         << "dydx="    << dydx0
627         << "dzdx="    << dzdx0
628         << "\n";
629     }
630     Float_t ymc = y0 - dx*dydx0;
631     Float_t zmc = z0 - dx*dzdx0;
632     //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
633
634     // Kalman position at reference radial position
635     dx = xAnode - x;
636     dydx = fTracklet->GetYref(1);
637     dzdx = fTracklet->GetZref(1);
638     dzdl = fTracklet->GetTgl();
639     y  = fTracklet->GetYref(0) - dx*dydx;
640     dy = y - ymc;
641     z  = fTracklet->GetZref(0) - dx*dzdx;
642     dz = z - zmc;
643     pt = TMath::Abs(fTracklet->GetPt());
644     fTracklet->GetCovRef(covR);
645
646     arr = (TObjArray*)fContainer->At(kMCtrackTRD);
647     // y resolution/pulls
648     ((TH2I*)arr->At(0))->Fill(dydx0, dy);
649     ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(covR[0]));
650     // z resolution/pulls
651     ((TH2I*)arr->At(2))->Fill(dzdx0, dz);
652     ((TH2I*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]));
653     // phi resolution/ snp pulls
654     Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
655     ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
656     Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0);
657     ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3]));
658     // theta resolution/ tgl pulls
659     Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
660               dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
661     ((TH2I*)arr->At(6))->Fill(dzdl0, 
662     TMath::ATan(dtgl));
663     ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
664     // pt resolution  \\ 1/pt pulls \\ p resolution for PID
665     for(Int_t is=AliPID::kSPECIES; is--;){
666       if(TMath::Abs(pdg)!=AliPID::ParticleCode(is)) continue;
667       ((TH3S*)((TObjArray*)arr->At(8))->At(ily))->Fill(pt0, pt/pt0-1., is);
668       ((TH3S*)((TObjArray*)arr->At(9))->At(ily))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), is);
669       Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
670                p  = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
671      ((TH3S*)((TObjArray*)arr->At(10))->At(ily))->Fill(p0, p/p0-1., is);
672       break;
673     }
674
675     // Fill Debug stream for Kalman track
676     if(fDebugLevel>=1){
677       (*fDebugStream) << "MCtrack"
678         << "pt="      << pt
679         << "x="       << x
680         << "y="       << y
681         << "z="       << z
682         << "dydx="    << dydx
683         << "dzdx="    << dzdx
684         << "s2y="     << covR[0]
685         << "s2z="     << covR[2]
686         << "\n";
687     }
688
689     // recalculate tracklet based on the MC info
690     AliTRDseedV1 tt(*fTracklet);
691     tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
692     tt.SetZref(1, dzdx0); 
693     tt.SetReconstructor(&rec);
694     tt.Fit(kTRUE, kTRUE);
695     x= tt.GetX();y= tt.GetY();z= tt.GetZ();
696     dydx = tt.GetYfit(1);
697     dx = x0 - x;
698     ymc = y0 - dx*dydx0;
699     zmc = z0 - dx*dzdx0;
700     Bool_t rc = tt.IsRowCross(); 
701     
702     // add tracklet residuals for y and dydx
703     arr = (TObjArray*)fContainer->At(kMCtracklet);
704     if(!rc){
705       dy    = y-ymc;
706
707       Float_t dphi  = (dydx - dydx0);
708       dphi /= (1.- dydx*dydx0);
709
710       ((TH2I*)arr->At(0))->Fill(dydx0, dy);
711       if(tt.GetS2Y()>0.) ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(tt.GetS2Y()));
712       ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
713     } else {
714       // add tracklet residuals for z
715       dz = z-zmc;
716       ((TH2I*)arr->At(2))->Fill(dzdl0, dz);
717       if(tt.GetS2Z()>0.) ((TH2I*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z()));
718     }
719   
720     // Fill Debug stream for tracklet
721     if(fDebugLevel>=1){
722       Float_t s2y = tt.GetS2Y();
723       Float_t s2z = tt.GetS2Z();
724       (*fDebugStream) << "MCtracklet"
725         << "rc="    << rc
726         << "x="     << x
727         << "y="     << y
728         << "z="     << z
729         << "dydx="  << dydx
730         << "s2y="   << s2y
731         << "s2z="   << s2z
732         << "\n";
733     }
734
735     Int_t istk = AliTRDgeometry::GetStack(det); 
736     AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
737     Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
738     Float_t tilt = fTracklet->GetTilt();
739     //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
740
741     arr = (TObjArray*)fContainer->At(kMCcluster);
742     AliTRDcluster *c = 0x0;
743     fTracklet->ResetClusterIter(kFALSE);
744     while((c = fTracklet->PrevCluster())){
745       Float_t  q = TMath::Abs(c->GetQ());
746       x = c->GetX(); y = c->GetY();z = c->GetZ();
747       dx = x0 - x; 
748       ymc= y0 - dx*dydx0;
749       zmc= z0 - dx*dzdx0;
750       dy = ymc - (y - tilt*(z-zmc));
751       // Fill Histograms
752       if(q>20. && q<250.){ 
753         ((TH2I*)arr->At(0))->Fill(dydx0, dy);
754         ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(c->GetSigmaY2()));
755       }
756
757       // Fill calibration container
758       Float_t d = zr0 - zmc;
759       d -= ((Int_t)(2 * d)) / 2.0;
760       if (d > 0.25) d  = 0.5 - d;
761       AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
762       fMCcl->Add(clInfo);
763       clInfo->SetCluster(c);
764       clInfo->SetMC(pdg, label);
765       clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
766       clInfo->SetResolution(dy);
767       clInfo->SetAnisochronity(d);
768       clInfo->SetDriftLength(((c->GetPadTime()+0.5)*.1)*1.5);
769       //dx-.5*AliTRDgeometry::CamHght());
770       clInfo->SetTilt(tilt);
771
772       // Fill Debug Tree
773       if(fDebugLevel>=2){
774         //clInfo->Print();
775         (*fDebugStream) << "MCcluster"
776           <<"clInfo.=" << clInfo
777           << "\n";
778       }
779     }
780   }
781   return h;
782 }
783
784
785 //________________________________________________________
786 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
787 {
788   Float_t xy[4] = {0., 0., 0., 0.};
789   if(!gPad){
790     AliWarning("Please provide a canvas to draw results.");
791     return kFALSE;
792   }
793   TList *l = 0x0; TVirtualPad *pad=0x0;
794   switch(ifig){
795   case kCharge:
796     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
797     ((TVirtualPad*)l->At(0))->cd();
798     ((TGraphAsymmErrors*)((TObjArray*)fGraphM->At(kCharge))->At(0))->Draw("apl");
799     ((TVirtualPad*)l->At(1))->cd();
800     ((TGraphErrors*)((TObjArray*)fGraphS->At(kCharge))->At(0))->Draw("apl");
801     break;
802   case kCluster:
803     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
804     xy[0] = -.3; xy[1] = -200.; xy[2] = .3; xy[3] = 1000.;
805     ((TVirtualPad*)l->At(0))->cd();
806     if(!GetGraphPlot(&xy[0], kCluster, 0)) break;
807     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
808     ((TVirtualPad*)l->At(1))->cd();
809     if(!GetGraphPlot(&xy[0], kCluster, 1)) break;
810     return kTRUE;
811   case kTrackTRD :
812     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
813     xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
814     ((TVirtualPad*)l->At(0))->cd();
815     if(!GetGraphPlot(&xy[0], kTrackTRD , 0)) break;
816     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
817     ((TVirtualPad*)l->At(1))->cd();
818     if(!GetGraphPlot(&xy[0], kTrackTRD , 1)) break;
819     return kTRUE;
820   case 3: // kTrackTRD  z
821     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
822     xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
823     ((TVirtualPad*)l->At(0))->cd();
824     if(!GetGraphPlot(&xy[0], kTrackTRD , 2)) break;
825     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
826     ((TVirtualPad*)l->At(1))->cd();
827     if(!GetGraphPlot(&xy[0], kTrackTRD , 3)) break;
828     return kTRUE;
829   case 4: // kTrackTRD  phi
830     xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
831     if(GetGraphPlot(&xy[0], kTrackTRD , 4)) return kTRUE;
832     break;
833   case 5: // kTrackTPC y
834     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
835     xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
836     pad = ((TVirtualPad*)l->At(0)); pad->cd();
837     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
838     if(!GetGraphPlot(&xy[0], kTrackTPC, 0)) break;
839     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
840     pad=((TVirtualPad*)l->At(1)); pad->cd();
841     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
842     if(!GetGraphPlot(&xy[0], kTrackTPC, 1)) break;
843     return kTRUE;
844   case 6: // kTrackTPC z
845     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
846     xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
847     pad = ((TVirtualPad*)l->At(0)); pad->cd();
848     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
849     if(!GetGraphPlot(&xy[0], kTrackTPC, 2)) break;
850     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
851     pad = ((TVirtualPad*)l->At(1)); pad->cd();
852     pad->SetMargin(0.1, 0.1, 0.1, 0.01);
853     if(!GetGraphPlot(&xy[0], kTrackTPC, 3)) break;
854     return kTRUE;
855   case 7: // kTrackTPC phi
856     xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
857     if(GetGraphPlot(&xy[0], kTrackTPC, 4)) return kTRUE;
858     break;
859   case 8: // kMCcluster
860     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
861     xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
862     ((TVirtualPad*)l->At(0))->cd();
863     if(!GetGraphPlot(&xy[0], kMCcluster, 0)) break;
864     xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
865     ((TVirtualPad*)l->At(1))->cd();
866     if(!GetGraphPlot(&xy[0], kMCcluster, 1)) break;
867     return kTRUE;
868   case 9: //kMCtracklet [y]
869     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
870     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =250.;
871     ((TVirtualPad*)l->At(0))->cd();
872     if(!GetGraphPlot(&xy[0], kMCtracklet, 0)) break;
873     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
874     ((TVirtualPad*)l->At(1))->cd();
875     if(!GetGraphPlot(&xy[0], kMCtracklet, 1)) break;
876     return kTRUE;
877   case 10: //kMCtracklet [z]
878     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
879     xy[0]=-1.; xy[1]=-100.; xy[2]=1.; xy[3] =2500.;
880     ((TVirtualPad*)l->At(0))->cd();
881     if(!GetGraphPlot(&xy[0], kMCtracklet, 2)) break;
882     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
883     ((TVirtualPad*)l->At(1))->cd();
884     if(!GetGraphPlot(&xy[0], kMCtracklet, 3)) break;
885     return kTRUE;
886   case 11: //kMCtracklet [phi]
887     xy[0]=-.3; xy[1]=-3.; xy[2]=.3; xy[3] =25.;
888     if(!GetGraphPlot(&xy[0], kMCtracklet, 4)) break;
889     return kTRUE;
890   case 12: //kMCtrackTRD [y]
891     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
892     xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =200.;
893     ((TVirtualPad*)l->At(0))->cd();
894     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 0)) break;
895     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
896     ((TVirtualPad*)l->At(1))->cd();
897     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 1)) break;
898     return kTRUE;
899   case 13: //kMCtrackTRD [z]
900     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
901     xy[0]=-1.; xy[1]=-700.; xy[2]=1.; xy[3] =1500.;
902     ((TVirtualPad*)l->At(0))->cd();
903     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 2)) break;
904     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
905     ((TVirtualPad*)l->At(1))->cd();
906     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 3)) break;
907     return kTRUE;
908   case 14: //kMCtrackTRD [phi/snp]
909     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
910     xy[0]=-.2; xy[1]=-0.2; xy[2]=.2; xy[3] =2.;
911     ((TVirtualPad*)l->At(0))->cd();
912     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 4)) break;
913     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
914     ((TVirtualPad*)l->At(1))->cd();
915     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 5)) break;
916     return kTRUE;
917   case 15: //kMCtrackTRD [theta/tgl]
918     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
919     xy[0]=-1.; xy[1]=-0.5; xy[2]=1.; xy[3] =5.;
920     ((TVirtualPad*)l->At(0))->cd();
921     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 6)) break;
922     xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
923     ((TVirtualPad*)l->At(1))->cd();
924     if(!GetGraphPlot(&xy[0], kMCtrackTRD, 7)) break;
925     return kTRUE;
926   case 16: //kMCtrackTRD [pt]
927     xy[0] = 0.; xy[1] = -5.; xy[2] = 12.; xy[3] = 7.;
928     gPad->Divide(2, 3, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
929     for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
930       pad = (TVirtualPad*)l->At(il); pad->cd();
931       pad->SetMargin(0.07, 0.07, 0.1, 0.);
932       if(!GetGraphTrack(&xy[0], 8, il)) break;
933     }
934     return kTRUE;
935   case 17: //kMCtrackTRD [1/pt] pulls
936     xy[0] = 0.; xy[1] = -1.5; xy[2] = 2.; xy[3] = 2.;
937     gPad->Divide(2, 3, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
938     for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
939       pad = (TVirtualPad*)l->At(il); pad->cd();
940       pad->SetMargin(0.07, 0.07, 0.1, 0.);
941       if(!GetGraphTrack(&xy[0], 9, il)) break;
942     }
943     return kTRUE;
944   case 18: //kMCtrackTRD [p]
945     xy[0] = 0.; xy[1] = -7.5; xy[2] = 12.; xy[3] = 10.5;
946     gPad->Divide(2, 3, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
947     for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
948       pad = (TVirtualPad*)l->At(il); pad->cd();
949       pad->SetMargin(0.07, 0.07, 0.1, 0.);
950       if(!GetGraphTrack(&xy[0], 10, il)) break;
951     }
952     return kTRUE;
953   case 19: // kMCtrackTPC [y]
954     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
955     xy[0]=-.25; xy[1]=-50.; xy[2]=.25; xy[3] =800.;
956     ((TVirtualPad*)l->At(0))->cd();
957     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 0)) break;
958     xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 2.5;
959     ((TVirtualPad*)l->At(1))->cd();
960     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 1)) break;
961     return kTRUE;
962   case 20: // kMCtrackTPC [z]
963     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
964     xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =800.;
965     ((TVirtualPad*)l->At(0))->cd();
966     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 2)) break;
967     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
968     ((TVirtualPad*)l->At(1))->cd();
969     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 3)) break;
970     return kTRUE;
971   case 21: // kMCtrackTPC [phi|snp]
972     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
973     xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
974     ((TVirtualPad*)l->At(0))->cd();
975     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 4)) break;
976     xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
977     ((TVirtualPad*)l->At(1))->cd();
978     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 5)) break;
979     return kTRUE;
980   case 22: // kMCtrackTPC [theta|tgl]
981     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
982     xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.;
983     ((TVirtualPad*)l->At(0))->cd();
984     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 6)) break;
985     xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 1.5;
986     ((TVirtualPad*)l->At(1))->cd();
987     if(!GetGraphPlot(&xy[0], kMCtrackTPC, 7)) break;
988     return kTRUE;
989   case 23: // kMCtrackTPC [pt]
990     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
991     xy[0] = 0.; xy[1] = -.8; xy[2] = 12.; xy[3] = 2.3;
992     ((TVirtualPad*)l->At(0))->cd();
993     if(!GetGraphTrackTPC(xy, 8)) break;
994     xy[0]=0.; xy[1]=-0.5; xy[2]=2.; xy[3] =2.5;
995     ((TVirtualPad*)l->At(1))->cd();
996     if(!GetGraphTrackTPC(xy, 9)) break;
997     return kTRUE;
998   case 24: // kMCtrackTPC [p]
999     gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); 
1000     xy[0] = 0.; xy[1] = -.8; xy[2] = 12.; xy[3] = 2.3;
1001     pad = ((TVirtualPad*)l->At(0));pad->cd();
1002     pad->SetMargin(0.12, 0.12, 0.1, 0.04);
1003     if(!GetGraphTrackTPC(xy, 10)) break;
1004     xy[0]=0.; xy[1]=-1.5; xy[2]=12.; xy[3] =2.5;
1005     pad = ((TVirtualPad*)l->At(1)); pad->cd();
1006     pad->SetMargin(0.12, 0.12, 0.1, 0.04);
1007     if(!GetGraphTrackTPC(xy, 11)) break;
1008     return kTRUE;
1009   case 25:  // kMCtrackTOF [z]
1010     return kTRUE;
1011   }
1012   AliWarning(Form("Reference plot [%d] missing result", ifig));
1013   return kFALSE;
1014 }
1015
1016
1017 //________________________________________________________
1018 Bool_t AliTRDresolution::PostProcess()
1019 {
1020   //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
1021   if (!fContainer) {
1022     Printf("ERROR: list not available");
1023     return kFALSE;
1024   }
1025   TGraph *gm= 0x0, *gs= 0x0;
1026   if(!fGraphS && !fGraphM){ 
1027     TObjArray *aM(0x0), *aS(0x0), *a(0x0);
1028     Int_t n = fContainer->GetEntriesFast();
1029     fGraphS = new TObjArray(n); fGraphS->SetOwner();
1030     fGraphM = new TObjArray(n); fGraphM->SetOwner();
1031     for(Int_t ig=0; ig<n; ig++){
1032       fGraphM->AddAt(aM = new TObjArray(fNElements[ig]), ig);
1033       fGraphS->AddAt(aS = new TObjArray(fNElements[ig]), ig);
1034
1035       for(Int_t ic=0; ic<fNElements[ig]; ic++){
1036         if(ig==kMCtrackTPC&&(ic>=8&&ic<=12)){ // TPC momentum plot
1037           aS->AddAt(a = new TObjArray(AliPID::kSPECIES), ic); 
1038           for(Int_t is=AliPID::kSPECIES; is--;){
1039             a->AddAt(gs = new TGraphErrors(), is);
1040             gs->SetMarkerStyle(23);
1041             gs->SetMarkerColor(is ? kRed : kMagenta);
1042             gs->SetLineStyle(is);
1043             gs->SetLineColor(is ? kRed : kMagenta);
1044             gs->SetLineWidth(is ? 1 : 3);
1045             gs->SetNameTitle(Form("s_%d%02d%d", ig, ic, is), "");
1046           }
1047           aM->AddAt(a = new TObjArray(AliPID::kSPECIES), ic); 
1048           for(Int_t is=AliPID::kSPECIES; is--;){
1049             a->AddAt(gm = new TGraphErrors(), is);
1050             gm->SetLineColor(is ? kBlack : kBlue);
1051             gm->SetLineStyle(is);
1052             gm->SetMarkerStyle(7);
1053             gm->SetMarkerColor(is ? kBlack : kBlue);
1054             gm->SetLineWidth(is ? 1 : 3);
1055             gm->SetNameTitle(Form("m_%d%02d%d", ig, ic, is), "");
1056           }
1057           continue;
1058         } else if(ig==kMCtrackTRD&&(ic==8||ic==9||ic==10)){ // TRD momentum plot
1059           TObjArray *aaS, *aaM;
1060           aS->AddAt(aaS = new TObjArray(AliTRDgeometry::kNlayer), ic); 
1061           aM->AddAt(aaM = new TObjArray(AliTRDgeometry::kNlayer), ic);
1062           for(Int_t il=AliTRDgeometry::kNlayer; il--;){
1063             aaS->AddAt(a = new TObjArray(AliPID::kSPECIES), il); 
1064             for(Int_t is=AliPID::kSPECIES; is--;){
1065               a->AddAt(gs = new TGraphErrors(), is);
1066               gs->SetMarkerStyle(23);
1067               gs->SetMarkerColor(is ? kRed : kMagenta);
1068               gs->SetLineStyle(is);
1069               gs->SetLineColor(is ? kRed : kMagenta);
1070               gs->SetLineWidth(is ? 1 : 3);
1071               gs->SetNameTitle(Form("s_%d%02d%d%d", ig, ic, is, il), "");
1072             }
1073             aaM->AddAt(a = new TObjArray(AliPID::kSPECIES), il);
1074             for(Int_t is=AliPID::kSPECIES; is--;){
1075               a->AddAt(gm = new TGraphErrors(), is);
1076               gm->SetMarkerStyle(7);
1077               gm->SetMarkerColor(is ? kBlack : kBlue);
1078               gm->SetLineStyle(is);
1079               gm->SetLineColor(is ? kBlack : kBlue);
1080               gm->SetLineWidth(is ? 1 : 3);
1081               gm->SetNameTitle(Form("m_%d%02d%d%d", ig, ic, is, il), "");
1082             }
1083           } 
1084           continue;
1085         }
1086
1087         aS->AddAt(gs = new TGraphErrors(), ic);
1088         gs->SetMarkerStyle(23);
1089         gs->SetMarkerColor(kRed);
1090         gs->SetLineColor(kRed);
1091         gs->SetNameTitle(Form("s_%d%02d", ig, ic), "");
1092
1093         aM->AddAt(gm = ig ? (TGraph*)new TGraphErrors() : (TGraph*)new TGraphAsymmErrors(), ic);
1094         gm->SetLineColor(kBlack);
1095         gm->SetMarkerStyle(7);
1096         gm->SetMarkerColor(kBlack);
1097         gm->SetNameTitle(Form("m_%d%02d", ig, ic), "");
1098       }
1099     }
1100   }
1101
1102 /*  printf("\n\n\n"); fGraphS->ls();
1103   printf("\n\n\n"); fGraphM->ls();*/
1104   
1105
1106   // DEFINE MODELS
1107   // simple gauss
1108   TF1 fg("fGauss", "gaus", -.5, .5);  
1109   // Landau for charge resolution
1110   TF1 fl("fLandau", "landau", 0., 1000.);  
1111
1112   //PROCESS EXPERIMENTAL DISTRIBUTIONS
1113   // Charge resolution
1114   //Process3DL(kCharge, 0, &fl); 
1115   // Clusters residuals
1116   Process2D(kCluster, 0, &fg, 1.e4); 
1117   Process2D(kCluster, 1, &fg); 
1118   fNRefFigures = 1;
1119   // Tracklet residual/pulls
1120   Process2D(kTrackTRD , 0, &fg, 1.e4); 
1121   Process2D(kTrackTRD , 1, &fg); 
1122   Process2D(kTrackTRD , 2, &fg, 1.e4); 
1123   Process2D(kTrackTRD , 3, &fg); 
1124   Process2D(kTrackTRD , 4, &fg, 1.e3); 
1125   fNRefFigures = 4;
1126   // TPC track residual/pulls
1127   Process2D(kTrackTPC, 0, &fg, 1.e4); 
1128   Process2D(kTrackTPC, 1, &fg); 
1129   Process2D(kTrackTPC, 2, &fg, 1.e4); 
1130   Process2D(kTrackTPC, 3, &fg); 
1131   Process2D(kTrackTPC, 4, &fg, 1.e3); 
1132   fNRefFigures = 7;
1133
1134   if(!HasMCdata()) return kTRUE;
1135
1136
1137   //PROCESS MC RESIDUAL DISTRIBUTIONS
1138
1139   // CLUSTER Y RESOLUTION/PULLS
1140   Process2D(kMCcluster, 0, &fg, 1.e4);
1141   Process2D(kMCcluster, 1, &fg);
1142   fNRefFigures = 8;
1143
1144   // TRACKLET RESOLUTION/PULLS
1145   Process2D(kMCtracklet, 0, &fg, 1.e4); // y
1146   Process2D(kMCtracklet, 1, &fg);       // y pulls
1147   Process2D(kMCtracklet, 2, &fg, 1.e4); // z
1148   Process2D(kMCtracklet, 3, &fg);       // z pulls
1149   Process2D(kMCtracklet, 4, &fg, 1.e3); // phi
1150   fNRefFigures = 11;
1151
1152   // TRACK RESOLUTION/PULLS
1153   Process2D(kMCtrackTRD, 0, &fg, 1.e4);   // y
1154   Process2D(kMCtrackTRD, 1, &fg);         // y PULL
1155   Process2D(kMCtrackTRD, 2, &fg, 1.e4);   // z
1156   Process2D(kMCtrackTRD, 3, &fg);         // z PULL
1157   Process2D(kMCtrackTRD, 4, &fg, 1.e3);   // phi
1158   Process2D(kMCtrackTRD, 5, &fg);         // snp PULL
1159   Process2D(kMCtrackTRD, 6, &fg, 1.e3);   // theta
1160   Process2D(kMCtrackTRD, 7, &fg);         // tgl PULL
1161 //   Process4D(kMCtrackTRD, 8, &fg, 1.e2);   // pt resolution
1162 //   Process4D(kMCtrackTRD, 9, &fg);         // 1/pt pulls
1163 //   Process4D(kMCtrackTRD, 10, &fg, 1.e2);  // p resolution
1164   fNRefFigures = 18;
1165
1166   // TRACK TPC RESOLUTION/PULLS
1167   Process2D(kMCtrackTPC, 0, &fg, 1.e4);// y resolution
1168   Process2D(kMCtrackTPC, 1, &fg);      // y pulls
1169   Process2D(kMCtrackTPC, 2, &fg, 1.e4);// z resolution
1170   Process2D(kMCtrackTPC, 3, &fg);      // z pulls
1171   Process2D(kMCtrackTPC, 4, &fg, 1.e3);// phi resolution
1172   Process2D(kMCtrackTPC, 5, &fg);      // snp pulls
1173   Process2D(kMCtrackTPC, 6, &fg, 1.e3);// theta resolution
1174   Process2D(kMCtrackTPC, 7, &fg);      // tgl pulls
1175 //   Process3D(kMCtrackTPC, 8, &fg, 1.e2);// pt resolution
1176 //   Process3D(kMCtrackTPC, 9, &fg);      // 1/pt pulls
1177 //   Process3D(kMCtrackTPC, 10, &fg, 1.e2);// p resolution
1178 //   Process3D(kMCtrackTPC, 11, &fg);      // p pulls
1179   fNRefFigures = 24;
1180
1181   // TRACK HMPID RESOLUTION/PULLS
1182   Process2D(kMCtrackTOF, 0, &fg, 1.e4); // z towards TOF
1183   Process2D(kMCtrackTOF, 1, &fg);       // z towards TOF
1184   fNRefFigures = 25;
1185
1186   return kTRUE;
1187 }
1188
1189
1190 //________________________________________________________
1191 void AliTRDresolution::Terminate(Option_t *)
1192 {
1193   if(fDebugStream){ 
1194     delete fDebugStream;
1195     fDebugStream = 0x0;
1196     fDebugLevel = 0;
1197   }
1198   if(HasPostProcess()) PostProcess();
1199 }
1200
1201 //________________________________________________________
1202 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1203 {
1204 // Helper function to avoid duplication of code
1205 // Make first guesses on the fit parameters
1206
1207   // find the intial parameters of the fit !! (thanks George)
1208   Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1209   Double_t sum = 0.;
1210   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1211   f->SetParLimits(0, 0., 3.*sum);
1212   f->SetParameter(0, .9*sum);
1213   Double_t rms = h->GetRMS();
1214   f->SetParLimits(1, -rms, rms);
1215   f->SetParameter(1, h->GetMean());
1216
1217   f->SetParLimits(2, 0., 2.*rms);
1218   f->SetParameter(2, rms);
1219   if(f->GetNpar() <= 4) return;
1220
1221   f->SetParLimits(3, 0., sum);
1222   f->SetParameter(3, .1*sum);
1223
1224   f->SetParLimits(4, -.3, .3);
1225   f->SetParameter(4, 0.);
1226
1227   f->SetParLimits(5, 0., 1.e2);
1228   f->SetParameter(5, 2.e-1);
1229 }
1230
1231 //________________________________________________________
1232 TObjArray* AliTRDresolution::Histos()
1233 {
1234   if(fContainer) return fContainer;
1235
1236   fContainer  = new TObjArray(kNhistos);
1237   //fContainer->SetOwner(kTRUE);
1238   TH1 *h = 0x0;
1239   TObjArray *arr = 0x0;
1240
1241   // cluster to track residuals/pulls
1242   fContainer->AddAt(arr = new TObjArray(fNElements[kCharge]), kCharge);
1243   arr->SetName("Charge");
1244   if(!(h = (TH3S*)gROOT->FindObject("hCharge"))){
1245     h = new TH3S("hCharge", "Charge Resolution", 20, 1., 2., 24, 0., 3.6, 100, 0., 500.);
1246     h->GetXaxis()->SetTitle("dx/dx_{0}");
1247     h->GetYaxis()->SetTitle("x_{d} [cm]");
1248     h->GetZaxis()->SetTitle("dq/dx [ADC/cm]");
1249   } else h->Reset();
1250   arr->AddAt(h, 0);
1251
1252   // cluster to track residuals/pulls
1253   fContainer->AddAt(arr = new TObjArray(fNElements[kCluster]), kCluster);
1254   arr->SetName("Cl");
1255   if(!(h = (TH2I*)gROOT->FindObject("hCl"))){
1256     h = new TH2I("hCl", "Cluster Residuals", 21, -.33, .33, 100, -.5, .5);
1257     h->GetXaxis()->SetTitle("tg(#phi)");
1258     h->GetYaxis()->SetTitle("#Delta y [cm]");
1259     h->GetZaxis()->SetTitle("entries");
1260   } else h->Reset();
1261   arr->AddAt(h, 0);
1262   if(!(h = (TH2I*)gROOT->FindObject("hClpull"))){
1263     h = new TH2I("hClpull", "Cluster Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1264     h->GetXaxis()->SetTitle("tg(#phi)");
1265     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1266     h->GetZaxis()->SetTitle("entries");
1267   } else h->Reset();
1268   arr->AddAt(h, 1);
1269
1270   // tracklet to track residuals/pulls in y direction
1271   fContainer->AddAt(arr = new TObjArray(fNElements[kTrackTRD ]), kTrackTRD );
1272   arr->SetName("Trklt");
1273   if(!(h = (TH2I*)gROOT->FindObject("hTrkltY"))){
1274     h = new TH2I("hTrkltY", "Tracklet Y Residuals", 21, -.33, .33, 100, -.5, .5);
1275     h->GetXaxis()->SetTitle("#tg(#phi)");
1276     h->GetYaxis()->SetTitle("#Delta y [cm]");
1277     h->GetZaxis()->SetTitle("entries");
1278   } else h->Reset();
1279   arr->AddAt(h, 0);
1280   if(!(h = (TH2I*)gROOT->FindObject("hTrkltYpull"))){
1281     h = new TH2I("hTrkltYpull", "Tracklet Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1282     h->GetXaxis()->SetTitle("#tg(#phi)");
1283     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1284     h->GetZaxis()->SetTitle("entries");
1285   } else h->Reset();
1286   arr->AddAt(h, 1);
1287   // tracklet to track residuals/pulls in z direction
1288   if(!(h = (TH2I*)gROOT->FindObject("hTrkltZ"))){
1289     h = new TH2I("hTrkltZ", "Tracklet Z Residuals", 50, -1., 1., 100, -1.5, 1.5);
1290     h->GetXaxis()->SetTitle("#tg(#theta)");
1291     h->GetYaxis()->SetTitle("#Delta z [cm]");
1292     h->GetZaxis()->SetTitle("entries");
1293   } else h->Reset();
1294   arr->AddAt(h, 2);
1295   if(!(h = (TH2I*)gROOT->FindObject("hTrkltZpull"))){
1296     h = new TH2I("hTrkltZpull", "Tracklet Z Pulls", 50, -1., 1., 100, -5.5, 5.5);
1297     h->GetXaxis()->SetTitle("#tg(#theta)");
1298     h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
1299     h->GetZaxis()->SetTitle("entries");
1300   } else h->Reset();
1301   arr->AddAt(h, 3);
1302   // tracklet to track phi residuals
1303   if(!(h = (TH2I*)gROOT->FindObject("hTrkltPhi"))){
1304     h = new TH2I("hTrkltPhi", "Tracklet #phi Residuals", 21, -.33, .33, 100, -.5, .5);
1305     h->GetXaxis()->SetTitle("tg(#phi)");
1306     h->GetYaxis()->SetTitle("#Delta phi [rad]");
1307     h->GetZaxis()->SetTitle("entries");
1308   } else h->Reset();
1309   arr->AddAt(h, 4);
1310
1311
1312   // tracklet to TPC track residuals/pulls in y direction
1313   fContainer->AddAt(arr = new TObjArray(fNElements[kTrackTPC]), kTrackTPC);
1314   arr->SetName("TrkTPC");
1315   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCY"))){
1316     h = new TH2I("hTrkTPCY", "Track[TPC] Y Residuals", 21, -.33, .33, 100, -.5, .5);
1317     h->GetXaxis()->SetTitle("#tg(#phi)");
1318     h->GetYaxis()->SetTitle("#Delta y [cm]");
1319     h->GetZaxis()->SetTitle("entries");
1320   } else h->Reset();
1321   arr->AddAt(h, 0);
1322   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCYpull"))){
1323     h = new TH2I("hTrkTPCYpull", "Track[TPC] Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1324     h->GetXaxis()->SetTitle("#tg(#phi)");
1325     h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1326     h->GetZaxis()->SetTitle("entries");
1327   } else h->Reset();
1328   arr->AddAt(h, 1);
1329   // tracklet to TPC track residuals/pulls in z direction
1330   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZ"))){
1331     h = new TH2I("hTrkTPCZ", "Track[TPC] Z Residuals", 50, -1., 1., 100, -1.5, 1.5);
1332     h->GetXaxis()->SetTitle("#tg(#theta)");
1333     h->GetYaxis()->SetTitle("#Delta z [cm]");
1334     h->GetZaxis()->SetTitle("entries");
1335   } else h->Reset();
1336   arr->AddAt(h, 2);
1337   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZpull"))){
1338     h = new TH2I("hTrkTPCZpull", "Track[TPC] Z Pulls", 50, -1., 1., 100, -5.5, 5.5);
1339     h->GetXaxis()->SetTitle("#tg(#theta)");
1340     h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
1341     h->GetZaxis()->SetTitle("entries");
1342   } else h->Reset();
1343   arr->AddAt(h, 3);
1344   // tracklet to TPC track phi residuals
1345   if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCPhi"))){
1346     h = new TH2I("hTrkTPCPhi", "Track[TPC] #phi Residuals", 21, -.33, .33, 100, -.5, .5);
1347     h->GetXaxis()->SetTitle("tg(#phi)");
1348     h->GetYaxis()->SetTitle("#Delta phi [rad]");
1349     h->GetZaxis()->SetTitle("entries");
1350   } else h->Reset();
1351   arr->AddAt(h, 4);
1352
1353
1354   // Resolution histos
1355   if(!HasMCdata()) return fContainer;
1356
1357   // cluster y resolution [0]
1358   fContainer->AddAt(arr = new TObjArray(fNElements[kMCcluster]), kMCcluster);
1359   arr->SetName("McCl");
1360   if(!(h = (TH2I*)gROOT->FindObject("hMcCl"))){
1361     h = new TH2I("hMcCl", "Cluster Resolution", 48, -.48, .48, 100, -.3, .3);
1362     h->GetXaxis()->SetTitle("tg(#phi)");
1363     h->GetYaxis()->SetTitle("#Delta y [cm]");
1364     h->GetZaxis()->SetTitle("entries");
1365   } else h->Reset();
1366   arr->AddAt(h, 0);
1367   if(!(h = (TH2I*)gROOT->FindObject("hMcClPull"))){
1368     h = new TH2I("hMcClPull", "Cluster Pulls", 48, -.48, .48, 100, -4.5, 4.5);
1369     h->GetXaxis()->SetTitle("tg(#phi)");
1370     h->GetYaxis()->SetTitle("#Deltay/#sigma_{y}");
1371     h->GetZaxis()->SetTitle("entries");
1372   } else h->Reset();
1373   arr->AddAt(h, 1);
1374
1375
1376   // TRACKLET RESOLUTION
1377   fContainer->AddAt(arr = new TObjArray(fNElements[kMCtracklet]), kMCtracklet);
1378   arr->SetName("McTrklt");
1379   // tracklet y resolution
1380   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltY"))){
1381     h = new TH2I("hMcTrkltY", "Tracklet Resolution (Y)", 48, -.48, .48, 100, -.2, .2);
1382     h->GetXaxis()->SetTitle("tg(#phi)");
1383     h->GetYaxis()->SetTitle("#Delta y [cm]");
1384     h->GetZaxis()->SetTitle("entries");
1385   } else h->Reset();
1386   arr->AddAt(h, 0);
1387   // tracklet y pulls
1388   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltYPull"))){
1389     h = new TH2I("hMcTrkltYPull", "Tracklet Pulls (Y)", 48, -.48, .48, 100, -4.5, 4.5);
1390     h->GetXaxis()->SetTitle("tg(#phi)");
1391     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1392     h->GetZaxis()->SetTitle("entries");
1393   } else h->Reset();
1394   arr->AddAt(h, 1);
1395   // tracklet z resolution
1396   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltZ"))){
1397     h = new TH2I("hMcTrkltZ", "Tracklet Resolution (Z)", 100, -1., 1., 100, -1., 1.);
1398     h->GetXaxis()->SetTitle("tg(#theta)");
1399     h->GetYaxis()->SetTitle("#Delta z [cm]");
1400     h->GetZaxis()->SetTitle("entries");
1401   } else h->Reset();
1402   arr->AddAt(h, 2);
1403   // tracklet z pulls
1404   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltZPull"))){
1405     h = new TH2I("hMcTrkltZPull", "Tracklet Pulls (Z)", 100, -1., 1., 100, -3.5, 3.5);
1406     h->GetXaxis()->SetTitle("tg(#theta)");
1407     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1408     h->GetZaxis()->SetTitle("entries");
1409   } else h->Reset();
1410   arr->AddAt(h, 3);
1411   // tracklet phi resolution
1412   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltPhi"))){
1413     h = new TH2I("hMcTrkltPhi", "Tracklet Resolution (#Phi)", 48, -.48, .48, 100, -.15, .15);
1414     h->GetXaxis()->SetTitle("tg(#phi)");
1415     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1416     h->GetZaxis()->SetTitle("entries");
1417   } else h->Reset();
1418   arr->AddAt(h, 4);
1419
1420
1421   // KALMAN TRACK RESOLUTION
1422   fContainer->AddAt(arr = new TObjArray(fNElements[kMCtrackTRD]), kMCtrackTRD);
1423   arr->SetName("McTrkTRD");
1424   // Kalman track y resolution
1425   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkY"))){
1426     h = new TH2I("hMcTrkY", "Track Y Resolution", 48, -.48, .48, 100, -.2, .2);
1427     h->GetXaxis()->SetTitle("tg(#phi)");
1428     h->GetYaxis()->SetTitle("#Delta y [cm]");
1429     h->GetZaxis()->SetTitle("entries");
1430   } else h->Reset();
1431   arr->AddAt(h, 0);
1432   // Kalman track y pulls
1433   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkYPull"))){
1434     h = new TH2I("hMcTrkYPull", "Track Y Pulls", 48, -.48, .48, 100, -4., 4.);
1435     h->GetXaxis()->SetTitle("tg(#phi)");
1436     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1437     h->GetZaxis()->SetTitle("entries");
1438   } else h->Reset();
1439   arr->AddAt(h, 1);
1440   // Kalman track Z
1441   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkZ"))){
1442     h = new TH2I("hMcTrkZ", "Track Z Resolution", 100, -1., 1., 100, -1., 1.);
1443     h->GetXaxis()->SetTitle("tg(#theta)");
1444     h->GetYaxis()->SetTitle("#Delta z [cm]");
1445     h->GetZaxis()->SetTitle("entries");
1446   } else h->Reset();
1447   arr->AddAt(h, 2);
1448   // Kalman track Z pulls
1449   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkZPull"))){
1450     h = new TH2I("hMcTrkZPull", "Track Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1451     h->GetXaxis()->SetTitle("tg(#theta)");
1452     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1453     h->GetZaxis()->SetTitle("entries");
1454   } else h->Reset();
1455   arr->AddAt(h, 3);
1456   // Kalman track SNP
1457   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNP"))){
1458     h = new TH2I("hMcTrkSNP", "Track Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3);
1459     h->GetXaxis()->SetTitle("tg(#phi)");
1460     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1461     h->GetZaxis()->SetTitle("entries");
1462   } else h->Reset();
1463   arr->AddAt(h, 4);
1464   // Kalman track SNP pulls
1465   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNPPull"))){
1466     h = new TH2I("hMcTrkSNPPull", "Track SNP Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1467     h->GetXaxis()->SetTitle("tg(#phi)");
1468     h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
1469     h->GetZaxis()->SetTitle("entries");
1470   } else h->Reset();
1471   arr->AddAt(h, 5);
1472   // Kalman track TGL
1473   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGL"))){
1474     h = new TH2I("hMcTrkTGL", "Track Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3);
1475     h->GetXaxis()->SetTitle("tg(#theta)");
1476     h->GetYaxis()->SetTitle("#Delta#theta [rad]");
1477     h->GetZaxis()->SetTitle("entries");
1478   } else h->Reset();
1479   arr->AddAt(h, 6);
1480   // Kalman track TGL pulls
1481   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGLPull"))){
1482     h = new TH2I("hMcTrkTGLPull", "Track TGL  Pulls", 100, -1., 1., 100, -4.5, 4.5);
1483     h->GetXaxis()->SetTitle("tg(#theta)");
1484     h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
1485     h->GetZaxis()->SetTitle("entries");
1486   } else h->Reset();
1487   arr->AddAt(h, 7);
1488   // Kalman track Pt resolution
1489   const Int_t n = AliPID::kSPECIES;
1490   TObjArray *arr2 = 0x0; TH3S* h3=0x0;
1491   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 8);
1492   arr2->SetName("Track Pt Resolution");
1493   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1494     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkPt%d", il)))){
1495       h3 = new TH3S(Form("hMcTrkPt%d", il), "Track Pt Resolution", 40, 0., 20., 150, -.1, .2, n, -.5, n-.5);
1496       h3->GetXaxis()->SetTitle("p_{t} [GeV/c]");
1497       h3->GetYaxis()->SetTitle("#Delta p_{t}/p_{t}^{MC}");
1498       h3->GetZaxis()->SetTitle("SPECIES");
1499     } else h3->Reset();
1500     arr2->AddAt(h3, il);
1501   }
1502   // Kalman track Pt pulls
1503   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 9);
1504   arr2->SetName("Track 1/Pt Pulls");
1505   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1506     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkPtPulls%d", il)))){
1507       h3 = new TH3S(Form("hMcTrkPtPulls%d", il), "Track 1/Pt Pulls", 40, 0., 2., 100, -4., 4., n, -.5, n-.5);
1508       h3->GetXaxis()->SetTitle("1/p_{t}^{MC} [c/GeV]");
1509       h3->GetYaxis()->SetTitle("#Delta(1/p_{t})/#sigma(1/p_{t}) ");
1510       h3->GetZaxis()->SetTitle("SPECIES");
1511     } else h3->Reset();
1512     arr2->AddAt(h3, il);
1513   }
1514   // Kalman track P resolution
1515   arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 10);
1516   arr2->SetName("Track P Resolution [PID]");
1517   for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1518     if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkP%d", il)))){
1519       h3 = new TH3S(Form("hMcTrkP%d", il), "Track P Resolution", 40, 0., 20., 150, -.15, .35, n, -.5, n-.5);
1520       h3->GetXaxis()->SetTitle("p [GeV/c]");
1521       h3->GetYaxis()->SetTitle("#Delta p/p^{MC}");
1522       h3->GetZaxis()->SetTitle("SPECIES");
1523     } else h3->Reset();
1524     arr2->AddAt(h3, il);
1525   }
1526
1527   // TPC TRACK RESOLUTION
1528   fContainer->AddAt(arr = new TObjArray(fNElements[kMCtrackTPC]), kMCtrackTPC);
1529   arr->SetName("McTrkTPC");
1530   // Kalman track Y
1531   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCY"))){
1532     h = new TH2I("hMcTrkTPCY", "Track[TPC] Y Resolution", 60, -.3, .3, 100, -.5, .5);
1533     h->GetXaxis()->SetTitle("tg(#phi)");
1534     h->GetYaxis()->SetTitle("#Delta y [cm]");
1535     h->GetZaxis()->SetTitle("entries");
1536   } else h->Reset();
1537   arr->AddAt(h, 0);
1538   // Kalman track Y pulls
1539   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCYPull"))){
1540     h = new TH2I("hMcTrkTPCYPull", "Track[TPC] Y Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1541     h->GetXaxis()->SetTitle("tg(#phi)");
1542     h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1543     h->GetZaxis()->SetTitle("entries");
1544   } else h->Reset();
1545   arr->AddAt(h, 1);
1546   // Kalman track Z
1547   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCZ"))){
1548     h = new TH2I("hMcTrkTPCZ", "Track[TPC] Z Resolution", 100, -1., 1., 100, -1., 1.);
1549     h->GetXaxis()->SetTitle("tg(#theta)");
1550     h->GetYaxis()->SetTitle("#Delta z [cm]");
1551     h->GetZaxis()->SetTitle("entries");
1552   } else h->Reset();
1553   arr->AddAt(h, 2);
1554   // Kalman track Z pulls
1555   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCZPull"))){
1556     h = new TH2I("hMcTrkTPCZPull", "Track[TPC] Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1557     h->GetXaxis()->SetTitle("tg(#theta)");
1558     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1559     h->GetZaxis()->SetTitle("entries");
1560   } else h->Reset();
1561   arr->AddAt(h, 3);
1562   // Kalman track SNP
1563   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCSNP"))){
1564     h = new TH2I("hMcTrkTPCSNP", "Track[TPC] Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3);
1565     h->GetXaxis()->SetTitle("tg(#phi)");
1566     h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1567     h->GetZaxis()->SetTitle("entries");
1568   } else h->Reset();
1569   arr->AddAt(h, 4);
1570   // Kalman track SNP pulls
1571   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCSNPPull"))){
1572     h = new TH2I("hMcTrkTPCSNPPull", "Track[TPC] SNP Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1573     h->GetXaxis()->SetTitle("tg(#phi)");
1574     h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
1575     h->GetZaxis()->SetTitle("entries");
1576   } else h->Reset();
1577   arr->AddAt(h, 5);
1578   // Kalman track TGL
1579   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCTGL"))){
1580     h = new TH2I("hMcTrkTPCTGL", "Track[TPC] Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3);
1581     h->GetXaxis()->SetTitle("tg(#theta)");
1582     h->GetYaxis()->SetTitle("#Delta#theta [rad]");
1583     h->GetZaxis()->SetTitle("entries");
1584   } else h->Reset();
1585   arr->AddAt(h, 6);
1586   // Kalman track TGL pulls
1587   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCTGLPull"))){
1588     h = new TH2I("hMcTrkTPCTGLPull", "Track[TPC] TGL  Pulls", 100, -1., 1., 100, -4.5, 4.5);
1589     h->GetXaxis()->SetTitle("tg(#theta)");
1590     h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
1591     h->GetZaxis()->SetTitle("entries");
1592   } else h->Reset();
1593   arr->AddAt(h, 7);
1594   // Kalman track Pt resolution
1595   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPt"))){
1596     h3 = new TH3S("hMcTrkTPCPt", "Track[TPC] Pt Resolution", 80, 0., 20., 150, -.1, .2, n, -.5, n-.5);
1597     h3->GetXaxis()->SetTitle("p_{t} [GeV/c]");
1598     h3->GetYaxis()->SetTitle("#Delta p_{t}/p_{t}^{MC}");
1599     h3->GetZaxis()->SetTitle("SPECIES");
1600   } else h3->Reset();
1601   arr->AddAt(h3, 8);
1602   // Kalman track Pt pulls
1603   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPtPulls"))){
1604     h3 = new TH3S("hMcTrkTPCPtPulls", "Track[TPC] 1/Pt Pulls", 80, 0., 2., 100, -4., 4., n, -.5, n-.5);
1605     h3->GetXaxis()->SetTitle("1/p_{t}^{MC} [c/GeV]");
1606     h3->GetYaxis()->SetTitle("#Delta(1/p_{t})/#sigma(1/p_{t}) ");
1607     h3->GetZaxis()->SetTitle("SPECIES");
1608   } else h3->Reset();
1609   arr->AddAt(h3, 9);
1610   // Kalman track P resolution
1611   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCP"))){
1612     h3 = new TH3S("hMcTrkTPCP", "Track[TPC] P Resolution", 80, 0., 20., 150, -.15, .35, n, -.5, n-.5);
1613     h3->GetXaxis()->SetTitle("p [GeV/c]");
1614     h3->GetYaxis()->SetTitle("#Delta p/p^{MC}");
1615     h3->GetZaxis()->SetTitle("SPECIES");
1616   } else h3->Reset();
1617   arr->AddAt(h3, 10);
1618   // Kalman track Pt pulls
1619   if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPPulls"))){
1620     h3 = new TH3S("hMcTrkTPCPPulls", "Track[TPC] P Pulls", 80, 0., 20., 100, -5., 5., n, -.5, n-.5);
1621     h3->GetXaxis()->SetTitle("p^{MC} [GeV/c]");
1622     h3->GetYaxis()->SetTitle("#Deltap/#sigma_{p}");
1623     h3->GetZaxis()->SetTitle("SPECIES");
1624   } else h3->Reset();
1625   arr->AddAt(h3, 11);
1626
1627
1628
1629   // Kalman track Z resolution [TOF]
1630   fContainer->AddAt(arr = new TObjArray(fNElements[kMCtrackTOF]), kMCtrackTOF);
1631   arr->SetName("McTrkTOF");
1632   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTOFZ"))){
1633     h = new TH2I("hMcTrkTOFZ", "Track[TOF] Z Resolution", 100, -1., 1., 100, -1., 1.);
1634     h->GetXaxis()->SetTitle("tg(#theta)");
1635     h->GetYaxis()->SetTitle("#Delta z [cm]");
1636     h->GetZaxis()->SetTitle("entries");
1637   } else h->Reset();
1638   arr->AddAt(h, 0);
1639   // Kalman track Z pulls
1640   if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTOFZPull"))){
1641     h = new TH2I("hMcTrkTOFZPull", "Track[TOF] Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1642     h->GetXaxis()->SetTitle("tg(#theta)");
1643     h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1644     h->GetZaxis()->SetTitle("entries");
1645   } else h->Reset();
1646   arr->AddAt(h, 1);
1647
1648   return fContainer;
1649 }
1650
1651 //________________________________________________________
1652 Bool_t AliTRDresolution::Process(TH2* h2, TF1 *f, Float_t k, TGraphErrors **g)
1653 {
1654   Char_t pn[10]; sprintf(pn, "p%02d", fIdxPlot);
1655   Int_t n = 0;
1656   if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1657   if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1658
1659   for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
1660     Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
1661     TH1D *h = h2->ProjectionY(pn, ibin, ibin);
1662     if(h->GetEntries()<100) continue;
1663     //AdjustF1(h, f);
1664
1665     h->Fit(f, "QN");
1666     
1667     Int_t ip = g[0]->GetN();
1668     g[0]->SetPoint(ip, x, k*f->GetParameter(1));
1669     g[0]->SetPointError(ip, 0., k*f->GetParError(1));
1670     g[1]->SetPoint(ip, x, k*f->GetParameter(2));
1671     g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1672     
1673 /*  
1674     g[0]->SetPoint(ip, x, k*h->GetMean());
1675     g[0]->SetPointError(ip, 0., k*h->GetMeanError());
1676     g[1]->SetPoint(ip, x, k*h->GetRMS());
1677     g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
1678   }
1679   fIdxPlot++;
1680   return kTRUE;
1681 }
1682
1683 //________________________________________________________
1684 Bool_t AliTRDresolution::Process2D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1685 {
1686   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1687
1688   // retrive containers
1689   TH2I *h2 = idx<0 ? (TH2I*)(fContainer->At(plot)) : (TH2I*)((TObjArray*)(fContainer->At(plot)))->At(idx);
1690   if(!h2) return kFALSE;
1691   TGraphErrors *g[2];
1692   if(!(g[0] = idx<0 ? (TGraphErrors*)fGraphM->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
1693
1694   if(!(g[1] = idx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
1695
1696   return Process(h2, f, k, g);
1697 }
1698
1699 //________________________________________________________
1700 Bool_t AliTRDresolution::Process3D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1701 {
1702   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1703
1704   // retrive containers
1705   TH3S *h3 = idx<0 ? (TH3S*)(fContainer->At(plot)) : (TH3S*)((TObjArray*)(fContainer->At(plot)))->At(idx);
1706   if(!h3) return kFALSE;
1707
1708   TObjArray *gm, *gs;
1709   if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
1710   if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
1711   TGraphErrors *g[2];
1712
1713   TAxis *az = h3->GetZaxis();
1714   for(Int_t iz=1; iz<=az->GetNbins(); iz++){
1715     if(!(g[0] = (TGraphErrors*)gm->At(iz-1))) return kFALSE;
1716     if(!(g[1] = (TGraphErrors*)gs->At(iz-1))) return kFALSE;
1717     az->SetRange(iz, iz);
1718     if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
1719   }
1720
1721   return kTRUE;
1722 }
1723
1724 //________________________________________________________
1725 Bool_t AliTRDresolution::Process3DL(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1726 {
1727   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1728
1729   // retrive containers
1730   TH3S *h3 = (TH3S*)((TObjArray*)fContainer->At(plot))->At(idx);
1731   if(!h3) return kFALSE;
1732
1733   TGraphAsymmErrors *gm; 
1734   TGraphErrors *gs;
1735   if(!(gm = (TGraphAsymmErrors*)((TObjArray*)fGraphM->At(plot))->At(0))) return kFALSE;
1736   if(!(gs = (TGraphErrors*)((TObjArray*)fGraphS->At(plot)))) return kFALSE;
1737
1738   Float_t x, r, mpv, xM, xm;
1739   TAxis *ay = h3->GetYaxis();
1740   for(Int_t iy=1; iy<=h3->GetNbinsY(); iy++){
1741     ay->SetRange(iy, iy);
1742     x = ay->GetBinCenter(iy);
1743     TH2F *h2=(TH2F*)h3->Project3D("zx");
1744     TAxis *ax = h2->GetXaxis();
1745     for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
1746       r = ax->GetBinCenter(ix);
1747       TH1D *h1 = h2->ProjectionY("py", ix, ix);
1748       if(h1->Integral()<50) continue;
1749       h1->Fit(f, "QN");
1750
1751       GetLandauMpvFwhm(f, mpv, xm, xM);
1752       Int_t ip = gm->GetN();
1753       gm->SetPoint(ip, x, k*mpv);
1754       gm->SetPointError(ip, 0., 0., k*xm, k*xM);
1755       gs->SetPoint(ip, r, k*(xM-xm)/mpv);
1756       gs->SetPointError(ip, 0., 0.);
1757     }
1758   }
1759
1760   return kTRUE;
1761 }
1762
1763 //________________________________________________________
1764 Bool_t AliTRDresolution::Process4D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1765 {
1766   if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1767
1768   // retrive containers
1769   TObjArray *arr = (TObjArray*)((TObjArray*)(fContainer->At(plot)))->At(idx);
1770   if(!arr) return kFALSE;
1771
1772
1773   TObjArray *gm[2], *gs[2];
1774   if(!(gm[0] = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
1775   if(!(gs[0] = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
1776
1777   TGraphErrors *g[2];
1778
1779   TH3S *h3 = 0x0;
1780   for(Int_t ix=0; ix<arr->GetEntriesFast(); ix++){
1781     if(!(h3 = (TH3S*)arr->At(ix))) return kFALSE;
1782     if(!(gm[1] = (TObjArray*)gm[0]->At(ix))) return kFALSE;
1783     if(!(gs[1] = (TObjArray*)gs[0]->At(ix))) return kFALSE;
1784     TAxis *az = h3->GetZaxis();
1785     for(Int_t iz=1; iz<=az->GetNbins(); iz++){
1786       if(!(g[0] = (TGraphErrors*)gm[1]->At(iz-1))) return kFALSE;
1787       if(!(g[1] = (TGraphErrors*)gs[1]->At(iz-1))) return kFALSE;
1788       az->SetRange(iz, iz);
1789       if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
1790     }
1791   }
1792
1793   return kTRUE;
1794 }
1795
1796 //________________________________________________________
1797 Bool_t AliTRDresolution::GetGraphPlot(Float_t *bb, ETRDresolutionPlot ip, Int_t idx)
1798 {
1799   if(!fGraphS || !fGraphM) return kFALSE;
1800   TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphM->At(ip)))->At(idx);
1801   if(!gm) return kFALSE;
1802   TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphS->At(ip)))->At(idx);
1803   if(!gs) return kFALSE;
1804   gs->Draw("apl"); gm->Draw("pl");
1805
1806   PutTrendValue(Form("%s_[%d]M", fPerformanceName[ip], idx), gm->GetMean(2), gm->GetRMS(2));
1807   gs->Sort(&TGraph::CompareY); Int_t n = gs->GetN();
1808   PutTrendValue(Form("%s_[%d]Sm", fPerformanceName[ip], idx), gs->GetY()[0], 0.);
1809   PutTrendValue(Form("%s_[%d]SM", fPerformanceName[ip], idx), gs->GetY()[n-1], 0.);
1810
1811   // axis titles look up
1812   Int_t nref = 0;
1813   for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fNElements[jp];
1814   UChar_t jdx = idx<0?0:idx;
1815   for(Int_t jc=0; jc<TMath::Min(jdx,fNElements[ip]-1); jc++) nref++;
1816   Char_t **at = fAxTitle[nref];
1817
1818   // axis range
1819   TAxis *ax = 0x0;
1820   ax = gs->GetHistogram()->GetXaxis();
1821   ax->SetRangeUser(bb[0], bb[2]);
1822   ax->SetTitle(*at);ax->CenterTitle();
1823
1824   ax = gs->GetHistogram()->GetYaxis();
1825   ax->SetRangeUser(bb[1], bb[3]);
1826   ax->SetTitleOffset(1.1);
1827   ax->SetTitle(*(++at));ax->CenterTitle();
1828
1829   TGaxis *gax = 0x0;
1830   gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
1831   gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
1832   //gax->SetVertical();
1833   gax->CenterTitle(); gax->SetTitleOffset(.7);
1834   gax->SetTitle(*(++at)); gax->Draw();
1835
1836   // bounding box
1837   TBox *b = new TBox(-.15, bb[1], .15, bb[3]);
1838   b->SetFillStyle(3002);b->SetLineColor(0);
1839   b->SetFillColor(ip<=Int_t(kMCcluster)?kGreen:kBlue);
1840   b->Draw();
1841   return kTRUE;
1842 }
1843
1844
1845 //________________________________________________________
1846 Bool_t AliTRDresolution::GetGraphTrack(Float_t *bb, Int_t idx, Int_t il)
1847 {
1848   if(!fGraphS || !fGraphM) return kFALSE;
1849
1850   TGraphErrors *gm = 0x0, *gs = 0x0;
1851   TObjArray *a0 = fGraphS, *a1 = 0x0;
1852   a1 = (TObjArray*)a0->At(kMCtrackTRD); a0 = a1;
1853   a1 = (TObjArray*)a0->At(idx); a0 = a1;
1854   a1 = (TObjArray*)a0->At(il); a0 = a1;
1855   for(Int_t is=0; is<AliPID::kSPECIES; is++){
1856     if(!(gs =  (TGraphErrors*)a0->At(is))) return kFALSE;
1857     if(!gs->GetN()) continue;
1858     gs->Draw(is ? "pl" : "apl");
1859   }
1860   gs =  (TGraphErrors*)a0->At(0);
1861   // axis titles look up
1862   Int_t nref = 0;
1863   for(Int_t jp=0; jp<Int_t(kMCtrackTRD); jp++) nref+=fNElements[jp];
1864   for(Int_t jc=0; jc<idx; jc++) nref++;
1865   Char_t **at = fAxTitle[nref];
1866   // axis range
1867   TAxis *ax = gs->GetHistogram()->GetXaxis();
1868   ax->SetRangeUser(bb[0], bb[2]);
1869   ax->SetTitle(*at);ax->CenterTitle();
1870
1871   ax = gs->GetHistogram()->GetYaxis();
1872   ax->SetRangeUser(bb[1], bb[3]);
1873   ax->SetTitleOffset(.5);ax->SetTitleSize(.06);
1874   ax->SetTitle(*(++at));ax->CenterTitle();
1875
1876   TGaxis *gax = 0x0;
1877   gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
1878   gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
1879   //gax->SetVertical();
1880   gax->CenterTitle(); gax->SetTitleOffset(.5);gax->SetTitleSize(.06);
1881   gax->SetTitle(*(++at)); gax->Draw();
1882
1883
1884   a0 = fGraphM;
1885   a1 = (TObjArray*)a0->At(kMCtrackTRD); a0 = a1;
1886   a1 = (TObjArray*)a0->At(idx); a0 = a1;
1887   a1 = (TObjArray*)a0->At(il); a0 = a1;
1888   for(Int_t is=0; is<AliPID::kSPECIES; is++){
1889     if(!(gm =  (TGraphErrors*)a0->At(is))) return kFALSE;
1890     if(!gm->GetN()) continue;
1891     gm->Draw("pl");
1892   }
1893
1894   return kTRUE;
1895 }
1896
1897
1898 //________________________________________________________
1899 Bool_t AliTRDresolution::GetGraphTrackTPC(Float_t *bb, Int_t sel)
1900 {
1901   if(!fGraphS || !fGraphM) return kFALSE;
1902
1903   TGraphErrors *gm = 0x0, *gs = 0x0;
1904   TObjArray *a0 = fGraphS, *a1 = 0x0;
1905   a1 = (TObjArray*)a0->At(kMCtrackTPC); a0 = a1;
1906   a1 = (TObjArray*)a0->At(sel); a0 = a1;
1907   for(Int_t is=0; is<AliPID::kSPECIES; is++){
1908     if(!(gs =  (TGraphErrors*)a0->At(is))) return kFALSE;
1909     if(!gs->GetN()) continue;
1910     gs->Draw(is ? "pl" : "apl");
1911   }
1912   gs =  (TGraphErrors*)a0->At(0);
1913   // axis titles look up
1914   Int_t nref = 0;
1915   for(Int_t jp=0; jp<Int_t(kMCtrackTPC); jp++) nref+=fNElements[jp];
1916   for(Int_t jc=0; jc<sel; jc++) nref++;
1917   Char_t **at = fAxTitle[nref];
1918   // axis range
1919   TAxis *ax = gs->GetHistogram()->GetXaxis();
1920   ax->SetRangeUser(bb[0], bb[2]);
1921   ax->SetTitle(*at);ax->CenterTitle();
1922
1923   ax = gs->GetHistogram()->GetYaxis();
1924   ax->SetRangeUser(bb[1], bb[3]);
1925   ax->SetTitleOffset(1.);ax->SetTitleSize(0.05);
1926   ax->SetTitle(*(++at));ax->CenterTitle();
1927
1928   TGaxis *gax = 0x0;
1929   gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
1930   gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
1931   //gax->SetVertical();
1932   gax->CenterTitle(); gax->SetTitleOffset(.7);gax->SetTitleSize(0.05);
1933   gax->SetTitle(*(++at)); gax->Draw();
1934
1935
1936   a0 = fGraphM;
1937   a1 = (TObjArray*)a0->At(kMCtrackTPC); a0 = a1;
1938   a1 = (TObjArray*)a0->At(sel); a0 = a1;
1939   for(Int_t is=0; is<AliPID::kSPECIES; is++){
1940     if(!(gm =  (TGraphErrors*)a0->At(is))) return kFALSE;
1941     if(!gm->GetN()) continue;
1942     gm->Draw("pl");
1943   }
1944
1945   return kTRUE;
1946 }
1947
1948 //________________________________________________________
1949 void AliTRDresolution::GetLandauMpvFwhm(TF1 *f, Float_t &mpv, Float_t &xm, Float_t &xM)
1950 {
1951   const Float_t dx = 1.;
1952   mpv = f->GetParameter(1);
1953   Float_t fx, max = f->Eval(mpv);
1954
1955   xm = mpv - dx;
1956   while((fx = f->Eval(xm))>.5*max){
1957     if(fx>max){ 
1958       max = fx;
1959       mpv = xm;
1960     }
1961     xm -= dx;
1962   }
1963
1964   xM += 2*(mpv - xm);
1965   while((fx = f->Eval(xM))>.5*max) xM += dx;
1966 }
1967
1968
1969 //________________________________________________________
1970 void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
1971 {
1972
1973   fReconstructor->SetRecoParam(r);
1974 }