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