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