]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/TRD/AliTRDresolution.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
new design for performance plots of the resolution task
[u/mrichter/AliRoot.git] / PWG1 / TRD / AliTRDresolution.cxx
CommitLineData
1ee39b3a 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
92c40e64 51#include <TSystem.h>
52
1ee39b3a 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>
92c40e64 62#include <TLegend.h>
1ee39b3a 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>
92c40e64 70#include <TDatabasePDG.h>
1ee39b3a 71
72#include "AliPID.h"
92d6d80c 73#include "AliLog.h"
92c40e64 74#include "AliESDtrack.h"
1ee39b3a 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"
92c40e64 84#include "AliTRDpidUtil.h"
1ee39b3a 85
86#include "info/AliTRDclusterInfo.h"
87
88ClassImp(AliTRDresolution)
89
92c40e64 90UChar_t const AliTRDresolution::fgNproj[kNviews] = {
1ee39b3a 91 2, 2, 5, 5,
92d6d80c 92 2, 5, 11, 11, 11
1ee39b3a 93};
92c40e64 94Char_t const * AliTRDresolution::fgPerformanceName[kNviews] = {
1ee39b3a 95 "Charge"
96 ,"Cluster2Track"
97 ,"Tracklet2Track"
98 ,"Tracklet2TPC"
99 ,"Cluster2MC"
100 ,"Tracklet2MC"
92d6d80c 101 ,"TRDin2MC"
102 ,"TRDout2MC"
1ee39b3a 103 ,"TRD2MC"
104};
92c40e64 105UChar_t const AliTRDresolution::fgNcomp[kNprojs] = {
92d6d80c 106 1, 1, //2,
107 14, 1, //2,
108 14, 1, 1, 1, 1, //5,
109 14, 1, 1, 1, 1, //5,
92c40e64 110// MC
92d6d80c 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
92c40e64 116};
117Char_t const *AliTRDresolution::fgAxTitle[kNprojs][4] = {
1ee39b3a 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
1295b37f 122 ,{"Cluster2Track residuals", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
123 ,{"Cluster2Track pulls", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
1ee39b3a 124 // TRD tracklet to Kalman fit
1295b37f 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]"}
1ee39b3a 130 // TPC track 2 first TRD tracklet
1295b37f 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]"}
1ee39b3a 136 // MC cluster
92d6d80c 137 ,{"MC Cluster Y resolution", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y} [#mum]"}
1295b37f 138 ,{"MC Cluster Y pulls", "tg(#phi)", "#mu_{y}", "#sigma_{y}"}
1ee39b3a 139 // MC tracklet
92d6d80c 140 ,{"MC Tracklet Y resolution", "tg(#phi)", "#mu_{y} [#mum]", "#sigma_{y}[#mum]"}
1295b37f 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]"}
92d6d80c 145 // MC track TRDin
1295b37f 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}) [%]"}
92d6d80c 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}) [%]"}
1ee39b3a 169 // MC track in TRD
1295b37f 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}) [%]"}
1ee39b3a 181};
182
183//________________________________________________________
184AliTRDresolution::AliTRDresolution()
f2e89a4c 185 :AliTRDrecoTask()
1ee39b3a 186 ,fStatus(0)
187 ,fIdxPlot(0)
92c40e64 188 ,fIdxFrame(0)
b91fdd71 189 ,fReconstructor(NULL)
190 ,fGeo(NULL)
92c40e64 191 ,fDBPDG(NULL)
b91fdd71 192 ,fGraphS(NULL)
193 ,fGraphM(NULL)
194 ,fCl(NULL)
195 ,fTrklt(NULL)
196 ,fMCcl(NULL)
197 ,fMCtrklt(NULL)
f8f46e4d 198{
199 //
200 // Default constructor
201 //
f2e89a4c 202 SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
f8f46e4d 203}
204
705f8b0a 205//________________________________________________________
f8f46e4d 206AliTRDresolution::AliTRDresolution(char* name)
f2e89a4c 207 :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
f8f46e4d 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)
1ee39b3a 220{
221 //
222 // Default constructor
223 //
224
225 fReconstructor = new AliTRDReconstructor();
226 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
227 fGeo = new AliTRDgeometry();
228
229 InitFunctorList();
230
705f8b0a 231 DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
232 DefineOutput(kTrkltToTrk, TObjArray::Class()); // tracklet2track
233 DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
234 DefineOutput(kTrkltToMC, TObjArray::Class()); // tracklet2mc
1ee39b3a 235}
236
237//________________________________________________________
238AliTRDresolution::~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//________________________________________________________
f8f46e4d 257void AliTRDresolution::UserCreateOutputObjects()
1ee39b3a 258{
259 // spatial resolution
f8f46e4d 260 OpenFile(1, "RECREATE");
1ee39b3a 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//________________________________________________________
f8f46e4d 274void AliTRDresolution::UserExec(Option_t *opt)
1ee39b3a 275{
276 //
277 // Execution part
278 //
279
280 fCl->Delete();
281 fTrklt->Delete();
282 fMCcl->Delete();
283 fMCtrklt->Delete();
b4414720 284 AliTRDrecoTask::UserExec(opt);
705f8b0a 285 PostData(kClToTrk, fCl);
286 PostData(kTrkltToTrk, fTrklt);
287 PostData(kClToMC, fMCcl);
288 PostData(kTrkltToMC, fMCtrklt);
1ee39b3a 289}
290
291//________________________________________________________
292TH1* AliTRDresolution::PlotCharge(const AliTRDtrackV1 *track)
293{
294 //
295 // Plots the charge distribution
296 //
297
298 if(track) fkTrack = track;
299 if(!fkTrack){
3d2a3dff 300 AliDebug(4, "No Track defined.");
b91fdd71 301 return NULL;
1ee39b3a 302 }
b91fdd71 303 TObjArray *arr = NULL;
1ee39b3a 304 if(!(arr = ((TObjArray*)fContainer->At(kCharge)))){
305 AliWarning("No output container defined.");
b91fdd71 306 return NULL;
1ee39b3a 307 }
b91fdd71 308 TH3S* h = NULL;
1ee39b3a 309
b91fdd71 310 AliTRDseedV1 *fTracklet = NULL;
311 AliTRDcluster *c = NULL;
1ee39b3a 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//________________________________________________________
337TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
338{
339 //
340 // Plot the cluster distributions
341 //
342
343 if(track) fkTrack = track;
344 if(!fkTrack){
3d2a3dff 345 AliDebug(4, "No Track defined.");
b91fdd71 346 return NULL;
1ee39b3a 347 }
b91fdd71 348 TObjArray *arr = NULL;
1ee39b3a 349 if(!(arr = ((TObjArray*)fContainer->At(kCluster)))){
350 AliWarning("No output container defined.");
b91fdd71 351 return NULL;
1ee39b3a 352 }
05039659 353 ULong_t status = fkESD ? fkESD->GetStatus():0;
1ee39b3a 354
31c8fa8a 355 Double_t covR[7], cov[3];
3d2a3dff 356 Float_t pt, x0, y0, z0, dy, dz, dydx, dzdx;
31c8fa8a 357 AliTRDseedV1 *fTracklet(NULL);
1ee39b3a 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();
3d2a3dff 362 pt = fTracklet->GetPt();
1ee39b3a 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);
31c8fa8a 369 fTracklet->GetCovRef(covR);
370 Float_t tilt(fTracklet->GetTilt());
b91fdd71 371 AliTRDcluster *c = NULL;
1ee39b3a 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;
31c8fa8a 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
3d2a3dff 402 ((TH3S*)arr->At(0))->Fill(dydx, dy/*, dz*/, pt);
1ee39b3a 403 ((TH2I*)arr->At(1))->Fill(dydx, dy/TMath::Sqrt(cov[0] /*+ sx2*/ + sy2));
404
df0514f6 405 if(DebugLevel()>=2){
1ee39b3a 406 // Get z-position with respect to anode wire
1ee39b3a 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);
db99a57a 417 Float_t covF[] = {cov[0], cov[1], cov[2]};
418 clInfo->SetGlobalPosition(yt, zt, dydx, dzdx, covF);
1ee39b3a 419 clInfo->SetResolution(dy);
420 clInfo->SetAnisochronity(d);
421 clInfo->SetDriftLength(dx);
31c8fa8a 422 clInfo->SetTilt(tilt);
df0514f6 423 (*DebugStream()) << "ClusterREC"
05039659 424 <<"status=" << status
1ee39b3a 425 <<"clInfo.=" << clInfo
426 << "\n";
df0514f6 427 }
1ee39b3a 428 }
429 }
430 return (TH2I*)arr->At(0);
431}
432
433
434//________________________________________________________
435TH1* 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){
3d2a3dff 447 AliDebug(4, "No Track defined.");
b91fdd71 448 return NULL;
1ee39b3a 449 }
b91fdd71 450 TObjArray *arr = NULL;
92d6d80c 451 if(!(arr = (TObjArray*)fContainer->At(kTrack ))){
1ee39b3a 452 AliWarning("No output container defined.");
b91fdd71 453 return NULL;
1ee39b3a 454 }
455
456 Double_t cov[3], covR[7]/*, sqr[3], inv[3]*/;
457 Float_t x, dx, dy, dz;
b91fdd71 458 AliTRDseedV1 *fTracklet = NULL;
1ee39b3a 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//________________________________________________________
494TH1* 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){
3d2a3dff 507 AliDebug(4, "No Track defined.");
b91fdd71 508 return NULL;
1ee39b3a 509 }
b91fdd71 510 AliExternalTrackParam *tin = NULL;
1ee39b3a 511 if(!(tin = fkTrack->GetTrackLow())){
512 AliWarning("Track did not entered TRD fiducial volume.");
b91fdd71 513 return NULL;
1ee39b3a 514 }
b91fdd71 515 TH1 *h = NULL;
1ee39b3a 516
517 Double_t x = tin->GetX();
b91fdd71 518 AliTRDseedV1 *tracklet = NULL;
1ee39b3a 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.");
b91fdd71 525 return NULL;
1ee39b3a 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();
92d6d80c 546 TObjArray *arr = (TObjArray*)fContainer->At(kTrackIn);
3d2a3dff 547 ((TH3S*)arr->At(0))->Fill(tracklet->GetYref(1), dy, 1./PAR[4]);
1ee39b3a 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;
92c40e64 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
1ee39b3a 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
92d6d80c 607 arr = (TObjArray*)fContainer->At(kMCtrackIn);
1ee39b3a 608 // y resolution/pulls
3d2a3dff 609 ((TH3S*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0], pt0);
1ee39b3a 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
92c40e64 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;
4226db3e 625 p = TMath::Sqrt(1.+ PAR[3]*PAR[3])/PAR[4];
92c40e64 626 ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., sign*sIdx);
afca20ef 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);
1ee39b3a 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//________________________________________________________
643TH1* 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.");
b91fdd71 651 return NULL;
1ee39b3a 652 }
653 if(track) fkTrack = track;
654 if(!fkTrack){
3d2a3dff 655 AliDebug(4, "No Track defined.");
b91fdd71 656 return NULL;
1ee39b3a 657 }
92c40e64 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;
92c40e64 667
668 TObjArray *arr(NULL);TH1 *h(NULL);
1ee39b3a 669 UChar_t s;
1ee39b3a 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){
5066aa9a 675 TVectorD dX(12), dY(12), dZ(12), Pt(12), dPt(12), cCOV(12*15);
676 fkMC->PropagateKalman(&dX, &dY, &dZ, &Pt, &dPt, &cCOV);
1ee39b3a 677 (*DebugStream()) << "MCkalman"
4226db3e 678 << "pdg=" << pdg
679 << "dx=" << &dX
680 << "dy=" << &dY
681 << "dz=" << &dZ
5066aa9a 682 << "pt=" << &Pt
4226db3e 683 << "dpt=" << &dPt
684 << "cov=" << &cCOV
1ee39b3a 685 << "\n";
686 }
687
688 AliTRDReconstructor rec;
d43e2ad1 689 AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
1ee39b3a 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
92c40e64 707 << "sgn=" << sign
1ee39b3a 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
92d6d80c 732 arr = (TObjArray*)((TObjArray*)fContainer->At(kMCtrack))->At(ily);
1ee39b3a 733 // y resolution/pulls
3d2a3dff 734 ((TH3S*)arr->At(0))->Fill(dydx0, dy, pt0);
1ee39b3a 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
92c40e64 751 Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
752 p = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
afca20ef 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);
1ee39b3a 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
92c40e64 792 ((TH3S*)arr->At(0))->Fill(dydx0, dy, pt0);
1ee39b3a 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);
b91fdd71 824 AliTRDcluster *c = NULL;
1ee39b3a 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;
df0514f6 832 dy = (y - tilt*(z-zmc)) - ymc;
1ee39b3a 833 // Fill Histograms
834 if(q>20. && q<250.){
92c40e64 835 ((TH3S*)arr->At(0))->Fill(dydx0, dy, pt0);
1ee39b3a 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;
1ee39b3a 844 clInfo->SetCluster(c);
845 clInfo->SetMC(pdg, label);
846 clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
847 clInfo->SetResolution(dy);
848 clInfo->SetAnisochronity(d);
d43e2ad1 849 clInfo->SetDriftLength(dx-.5*AliTRDgeometry::CamHght());
1ee39b3a 850 clInfo->SetTilt(tilt);
d43e2ad1 851 fMCcl->Add(clInfo);
852 if(DebugLevel()>=2){
853 if(!clInfoArr) clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
854 clInfoArr->Add(clInfo);
1ee39b3a 855 }
856 }
d43e2ad1 857 // Fill Debug Tree
858 if(DebugLevel()>=2 && clInfoArr){
859 (*DebugStream()) << "MCcluster"
860 <<"clInfo.=" << clInfoArr
861 << "\n";
92c40e64 862 clInfoArr->Clear();
d43e2ad1 863 }
1ee39b3a 864 }
92c40e64 865 if(clInfoArr) delete clInfoArr;
1ee39b3a 866 return h;
867}
868
869
870//________________________________________________________
871Bool_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 }
92d6d80c 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);
1ee39b3a 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();
92d6d80c 891 ga=((TGraphAsymmErrors*)((TObjArray*)fGraphM->At(kCharge))->At(0));
892 if(ga->GetN()) ga->Draw("apl");
1ee39b3a 893 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 894 g = ((TGraphErrors*)((TObjArray*)fGraphS->At(kCharge))->At(0));
895 if(g->GetN()) g->Draw("apl");
1ee39b3a 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();
92d6d80c 901 if(!GetGraphArray(xy, kCluster, 0)) break;
1ee39b3a 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;
92d6d80c 906 case kTrack: // kTrack y
1ee39b3a 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();
92d6d80c 910 if(!GetGraphArray(xy, kTrack , 0)) break;
1ee39b3a 911 xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
912 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 913 if(!GetGraphPlot(&xy[0], kTrack, 1)) break;
1ee39b3a 914 return kTRUE;
92d6d80c 915 case 3: // kTrack z
1ee39b3a 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();
92d6d80c 919 if(!GetGraphPlot(&xy[0], kTrack , 2)) break;
1ee39b3a 920 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
921 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 922 if(!GetGraphPlot(&xy[0], kTrack , 3)) break;
1ee39b3a 923 return kTRUE;
92d6d80c 924 case 4: // kTrack phi
1ee39b3a 925 xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
92d6d80c 926 if(GetGraphPlot(&xy[0], kTrack , 4)) return kTRUE;
1ee39b3a 927 break;
92d6d80c 928 case 5: // kTrackIn y
1ee39b3a 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);
92d6d80c 933 if(!GetGraphArray(xy, kTrackIn, 0)) break;
1ee39b3a 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);
92d6d80c 937 if(!GetGraphPlot(&xy[0], kTrackIn, 1)) break;
1ee39b3a 938 return kTRUE;
92d6d80c 939 case 6: // kTrackIn z
1ee39b3a 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);
92d6d80c 944 if(!GetGraphPlot(&xy[0], kTrackIn, 2)) break;
1ee39b3a 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);
92d6d80c 948 if(!GetGraphPlot(&xy[0], kTrackIn, 3)) break;
1ee39b3a 949 return kTRUE;
92d6d80c 950 case 7: // kTrackIn phi
1ee39b3a 951 xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
92d6d80c 952 if(GetGraphPlot(&xy[0], kTrackIn, 4)) return kTRUE;
1ee39b3a 953 break;
92d6d80c 954 case 8: // kMCcluster
92c40e64 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();
92d6d80c 958 if(!GetGraphArray(xy, kMCcluster, 0)) break;
92c40e64 959 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 960 xy[0]=-.3; xy[1]=-0.5; xy[2]=.3; xy[3]=2.5;
961 if(!GetGraphPlot(xy, kMCcluster, 1)) break;
92c40e64 962 return kTRUE;
92d6d80c 963 case 9: //kMCtracklet [y]
1ee39b3a 964 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
92d6d80c 965 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =500.;
1ee39b3a 966 ((TVirtualPad*)l->At(0))->cd();
92d6d80c 967 if(!GetGraphArray(xy, kMCtracklet, 0)) break;
1ee39b3a 968 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 969 xy[0]=-.3; xy[1]=-0.5; xy[2]=.3; xy[3]=2.5;
970 if(!GetGraphPlot(xy, kMCtracklet, 1)) break;
1ee39b3a 971 return kTRUE;
92d6d80c 972 case 10: //kMCtracklet [z]
92c40e64 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();
92d6d80c 976 if(!GetGraphPlot(&xy[0], kMCtracklet, 2)) break;
92c40e64 977 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
978 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 979 if(!GetGraphPlot(&xy[0], kMCtracklet, 3)) break;
92c40e64 980 return kTRUE;
92d6d80c 981 case 11: //kMCtracklet [phi]
92c40e64 982 xy[0]=-.3; xy[1]=-3.; xy[2]=.3; xy[3] =25.;
92d6d80c 983 if(!GetGraphPlot(&xy[0], kMCtracklet, 4)) break;
1ee39b3a 984 return kTRUE;
92d6d80c 985 case 12: //kMCtrack [y]
1ee39b3a 986 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1295b37f 987 xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1ee39b3a 988 ((TVirtualPad*)l->At(0))->cd();
92d6d80c 989 if(!GetGraphArray(xy, kMCtrack, 0)) break;
1295b37f 990 xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 3.5;
1ee39b3a 991 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 992 if(!GetGraphArray(xy, kMCtrack, 1)) break;
1ee39b3a 993 return kTRUE;
92d6d80c 994 case 13: //kMCtrack [z]
1ee39b3a 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();
92d6d80c 998 if(!GetGraphArray(xy, kMCtrack, 2)) break;
1ee39b3a 999 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1000 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 1001 if(!GetGraphArray(xy, kMCtrack, 3)) break;
1ee39b3a 1002 return kTRUE;
92d6d80c 1003 case 14: //kMCtrack [phi/snp]
1ee39b3a 1004 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
92d6d80c 1005 xy[0]=-.2; xy[1]=-0.5; xy[2]=.2; xy[3] =10.;
1ee39b3a 1006 ((TVirtualPad*)l->At(0))->cd();
92d6d80c 1007 if(!GetGraphArray(xy, kMCtrack, 4)) break;
1ee39b3a 1008 xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
1009 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 1010 if(!GetGraphArray(xy, kMCtrack, 5)) break;
1ee39b3a 1011 return kTRUE;
92d6d80c 1012 case 15: //kMCtrack [theta/tgl]
1ee39b3a 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();
92d6d80c 1016 if(!GetGraphArray(xy, kMCtrack, 6)) break;
1ee39b3a 1017 xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
1018 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 1019 if(!GetGraphArray(xy, kMCtrack, 7)) break;
1ee39b3a 1020 return kTRUE;
92d6d80c 1021 case 16: //kMCtrack [pt]
92c40e64 1022 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1023 pad = (TVirtualPad*)l->At(0); pad->cd();
1295b37f 1024 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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;
31c8fa8a 1034 pad->Modified(); pad->Update(); pad->SetLogx();
92c40e64 1035 pad = (TVirtualPad*)l->At(1); pad->cd();
1295b37f 1036 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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;
31c8fa8a 1044 pad->Modified(); pad->Update(); pad->SetLogx();
92c40e64 1045 return kTRUE;
92d6d80c 1046 case 17: //kMCtrack [pt]
92c40e64 1047 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1048 pad = (TVirtualPad*)l->At(0); pad->cd();
1295b37f 1049 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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();
92c40e64 1060 pad = (TVirtualPad*)l->At(1); pad->cd();
1295b37f 1061 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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();
92c40e64 1072 return kTRUE;
92d6d80c 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
92c40e64 1076 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1077 pad = (TVirtualPad*)l->At(0); pad->cd();
1295b37f 1078 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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;
92c40e64 1086 pad = (TVirtualPad*)l->At(1); pad->cd();
1295b37f 1087 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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;
92c40e64 1095 return kTRUE;
92d6d80c 1096 case 19: //kMCtrack [1/pt] pulls
92c40e64 1097 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1098 pad = (TVirtualPad*)l->At(0); pad->cd();
1295b37f 1099 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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;
92c40e64 1109 pad = (TVirtualPad*)l->At(1); pad->cd();
1295b37f 1110 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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;
1ee39b3a 1119 return kTRUE;
92d6d80c 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.;
92c40e64 1123 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1124 pad = (TVirtualPad*)l->At(0); pad->cd();
1295b37f 1125 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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();
92c40e64 1134 pad = (TVirtualPad*)l->At(1); pad->cd();
1295b37f 1135 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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();
1ee39b3a 1144 return kTRUE;
92d6d80c 1145 case 21: //kMCtrack [p]
92c40e64 1146 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1147 pad = (TVirtualPad*)l->At(0); pad->cd();
1295b37f 1148 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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();
92c40e64 1159 pad = (TVirtualPad*)l->At(1); pad->cd();
1295b37f 1160 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 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();
1ee39b3a 1171 return kTRUE;
92d6d80c 1172 case 22: // kMCtrackIn [y]
1ee39b3a 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();
92d6d80c 1176 if(!GetGraphArray(xy, kMCtrackIn, 0)) break;
1ee39b3a 1177 xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 2.5;
1178 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 1179 if(!GetGraphPlot(&xy[0], kMCtrackIn, 1)) break;
1ee39b3a 1180 return kTRUE;
92d6d80c 1181 case 23: // kMCtrackIn [z]
1ee39b3a 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();
92d6d80c 1185 if(!GetGraphPlot(&xy[0], kMCtrackIn, 2)) break;
1ee39b3a 1186 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1187 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 1188 if(!GetGraphPlot(&xy[0], kMCtrackIn, 3)) break;
1ee39b3a 1189 return kTRUE;
92d6d80c 1190 case 24: // kMCtrackIn [phi|snp]
1ee39b3a 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();
92d6d80c 1194 if(!GetGraphPlot(&xy[0], kMCtrackIn, 4)) break;
1ee39b3a 1195 xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
1196 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 1197 if(!GetGraphPlot(&xy[0], kMCtrackIn, 5)) break;
1ee39b3a 1198 return kTRUE;
92d6d80c 1199 case 25: // kMCtrackIn [theta|tgl]
1ee39b3a 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();
92d6d80c 1203 if(!GetGraphPlot(&xy[0], kMCtrackIn, 6)) break;
1ee39b3a 1204 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 1.5;
1205 ((TVirtualPad*)l->At(1))->cd();
92d6d80c 1206 if(!GetGraphPlot(&xy[0], kMCtrackIn, 7)) break;
1ee39b3a 1207 return kTRUE;
92d6d80c 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]
1ee39b3a 1283 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
db99a57a 1284 xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
92c40e64 1285 pad=(TVirtualPad*)l->At(0); pad->cd();
1295b37f 1286 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 1287 if(!GetGraphArray(xy, kMCtrackOut, 8)) break;
92c40e64 1288 pad->SetLogx();
1295b37f 1289 pad = (TVirtualPad*)l->At(1); pad->cd();
1290 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 1291 if(!GetGraphArray(xy, kMCtrackOut, 8)) break;
1ee39b3a 1292 return kTRUE;
92d6d80c 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]
1ee39b3a 1304 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
db99a57a 1305 xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1ee39b3a 1306 pad = ((TVirtualPad*)l->At(0));pad->cd();
1295b37f 1307 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 1308 if(!GetGraphArray(xy, kMCtrackOut, 10)) break;
92c40e64 1309 pad->SetLogx();
1ee39b3a 1310 pad = ((TVirtualPad*)l->At(1)); pad->cd();
1295b37f 1311 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
92d6d80c 1312 if(!GetGraphArray(xy, kMCtrackOut, 10)) break;
1ee39b3a 1313 return kTRUE;
1314 }
1315 AliWarning(Form("Reference plot [%d] missing result", ifig));
1316 return kFALSE;
1317}
1318
92c40e64 1319Char_t const *fgParticle[11]={
1320 " p bar", " K -", " #pi -", " #mu -", " e -",
1321 " No PID",
1322 " e +", " #mu +", " #pi +", " K +", " p",
1323};
1295b37f 1324const Color_t fgColorS[11]={
1325kOrange, kOrange-3, kMagenta+1, kViolet, kRed,
31c8fa8a 1326kGray,
1295b37f 1327kRed, kViolet, kMagenta+1, kOrange-3, kOrange
1328};
1329const Color_t fgColorM[11]={
1330kCyan-5, kAzure-4, kBlue-7, kBlue+2, kViolet+10,
1331kBlack,
1332kViolet+10, kBlue+2, kBlue-7, kAzure-4, kCyan-5
31c8fa8a 1333};
1334const Marker_t fgMarker[11]={
133530, 30, 26, 25, 24,
133628,
133720, 21, 22, 29, 29
1338};
1ee39b3a 1339//________________________________________________________
1340Bool_t AliTRDresolution::PostProcess()
1341{
1342 //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
1343 if (!fContainer) {
1344 AliError("ERROR: list not available");
1345 return kFALSE;
1346 }
b91fdd71 1347 TGraph *gm= NULL, *gs= NULL;
1ee39b3a 1348 if(!fGraphS && !fGraphM){
31c8fa8a 1349 TObjArray *aM(NULL), *aS(NULL);
1ee39b3a 1350 Int_t n = fContainer->GetEntriesFast();
1351 fGraphS = new TObjArray(n); fGraphS->SetOwner();
1352 fGraphM = new TObjArray(n); fGraphM->SetOwner();
92d6d80c 1353 for(Int_t ig(0), nc(0); ig<n; ig++){
92c40e64 1354 fGraphM->AddAt(aM = new TObjArray(fgNproj[ig]), ig);
1355 fGraphS->AddAt(aS = new TObjArray(fgNproj[ig]), ig);
92d6d80c 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]));
92c40e64 1359 if(fgNcomp[nc]>1){
31c8fa8a 1360 TObjArray *agS(NULL), *agM(NULL);
1361 aS->AddAt(agS = new TObjArray(fgNcomp[nc]), ic);
1362 aM->AddAt(agM = new TObjArray(fgNcomp[nc]), ic);
92c40e64 1363 for(Int_t is=fgNcomp[nc]; is--;){
31c8fa8a 1364 agS->AddAt(gs = new TGraphErrors(), is);
92d6d80c 1365 Int_t is0(is%11), il0(is/11);
31c8fa8a 1366 gs->SetMarkerStyle(fgMarker[is0]);
1295b37f 1367 gs->SetMarkerColor(fgColorS[is0]);
1368 gs->SetLineColor(fgColorS[is0]);
92d6d80c 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));
31c8fa8a 1371
1372 agM->AddAt(gm = new TGraphErrors(), is);
1373 gm->SetMarkerStyle(fgMarker[is0]);
1295b37f 1374 gm->SetMarkerColor(fgColorM[is0]);
1375 gm->SetLineColor(fgColorM[is0]);
92d6d80c 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));
92c40e64 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), "");
1ee39b3a 1391 }
1ee39b3a 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
92c40e64 1404 TF1 fch("fClCh", "landau", 0., 1000.);
1405 // Landau for e+- pt resolution
1406 TF1 fpt("fPt", "landau", -0.1, 0.2);
1ee39b3a 1407
1408 //PROCESS EXPERIMENTAL DISTRIBUTIONS
1409 // Charge resolution
1410 //Process3DL(kCharge, 0, &fl);
1411 // Clusters residuals
92d6d80c 1412 Process3D(kCluster, 0, &fg, 1.e4);
1ee39b3a 1413 Process2D(kCluster, 1, &fg);
92d6d80c 1414 fNRefFigures = 2;
1ee39b3a 1415 // Tracklet residual/pulls
92d6d80c 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;
1ee39b3a 1422 // TPC track residual/pulls
92d6d80c 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;
1ee39b3a 1429
1430 if(!HasMCdata()) return kTRUE;
1431
1432
1433 //PROCESS MC RESIDUAL DISTRIBUTIONS
1434
1435 // CLUSTER Y RESOLUTION/PULLS
92d6d80c 1436 Process3D(kMCcluster, 0, &fg, 1.e4);
1437 Process2D(kMCcluster, 1, &fg, 1.);
1438 fNRefFigures = 9;
1ee39b3a 1439
1440 // TRACKLET RESOLUTION/PULLS
92d6d80c 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;
1ee39b3a 1447
1448 // TRACK RESOLUTION/PULLS
92d6d80c 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;
1ee39b3a 1489
1490 return kTRUE;
1491}
1492
1493
1494//________________________________________________________
1495void AliTRDresolution::Terminate(Option_t *opt)
1496{
1497 AliTRDrecoTask::Terminate(opt);
1498 if(HasPostProcess()) PostProcess();
1499}
1500
1501//________________________________________________________
1502void 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
afca20ef 1531//________________________________________________________
3d2a3dff 1532TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name)
afca20ef 1533{
3d2a3dff 1534// Build performance histograms for AliTRDcluster.vs TRD track or MC
afca20ef 1535// - y reziduals/pulls
3d2a3dff 1536
1537 TObjArray *arr = new TObjArray(2);
afca20ef 1538 arr->SetName(name); arr->SetOwner();
1539 TH1 *h(NULL); char hname[100], htitle[300];
1540
3d2a3dff 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
afca20ef 1548 // tracklet resolution/pull in y direction
3d2a3dff 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);
afca20ef 1554 } else h->Reset();
1555 arr->AddAt(h, 0);
3d2a3dff 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);
afca20ef 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
3d2a3dff 1563 return arr;
1564}
1565
1566//________________________________________________________
1567TObjArray* 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
afca20ef 1577 // tracklet resolution/pull in z direction
3d2a3dff 1578 sprintf(hname, "%s_%s_Z", GetNameId(), name);
1579 sprintf(htitle, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm];entries", GetNameId(), name);
afca20ef 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);
3d2a3dff 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);
afca20ef 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
3d2a3dff 1592 sprintf(hname, "%s_%s_PHI", GetNameId(), name);
1593 sprintf(htitle, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];entries", GetNameId(), name);
afca20ef 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//________________________________________________________
1603TObjArray* 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
afca20ef 1612 TObjArray *arr = BuildMonitorContainerTracklet(name);
1613 arr->Expand(11);
3d2a3dff 1614 TH1 *h(NULL); char hname[100], htitle[300];
1615
afca20ef 1616 // snp pulls
3d2a3dff 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);
afca20ef 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
3d2a3dff 1625 sprintf(hname, "%s_%s_THT", GetNameId(), name);
1626 sprintf(htitle, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
afca20ef 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
3d2a3dff 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);
afca20ef 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
3d2a3dff 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);
afca20ef 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
3d2a3dff 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);
afca20ef 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
3d2a3dff 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);
afca20ef 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
1ee39b3a 1677//________________________________________________________
1678TObjArray* AliTRDresolution::Histos()
1679{
1680 //
1681 // Define histograms
1682 //
1683
1684 if(fContainer) return fContainer;
1685
92c40e64 1686 fContainer = new TObjArray(kNviews);
1ee39b3a 1687 //fContainer->SetOwner(kTRUE);
31c8fa8a 1688 TH1 *h(NULL);
1689 TObjArray *arr(NULL);
1ee39b3a 1690
31c8fa8a 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;
df0514f6 1698
1ee39b3a 1699 // cluster to track residuals/pulls
92d6d80c 1700 fContainer->AddAt(arr = new TObjArray(2), kCharge);
1ee39b3a 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
3d2a3dff 1711 fContainer->AddAt(BuildMonitorContainerCluster("Cl"), kCluster);
afca20ef 1712 // tracklet to TRD track
92d6d80c 1713 fContainer->AddAt(BuildMonitorContainerTracklet("Trk"), kTrack);
afca20ef 1714 // tracklet to TRDin
92d6d80c 1715 fContainer->AddAt(BuildMonitorContainerTracklet("TrkIN"), kTrackIn);
1ee39b3a 1716
1717
1718 // Resolution histos
1719 if(!HasMCdata()) return fContainer;
1720
3d2a3dff 1721 // cluster resolution
1722 fContainer->AddAt(BuildMonitorContainerCluster("MCcl"), kMCcluster);
1ee39b3a 1723
3d2a3dff 1724 // tracklet resolution
1725 fContainer->AddAt(BuildMonitorContainerTracklet("MCtracklet"), kMCtracklet);
1ee39b3a 1726
3d2a3dff 1727 // track resolution
92d6d80c 1728 fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
3d2a3dff 1729 arr->SetName("MCtrk");
1730 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
31c8fa8a 1731
afca20ef 1732 // TRDin TRACK RESOLUTION
92d6d80c 1733 fContainer->AddAt(BuildMonitorContainerTrack("MCtrkIN"), kMCtrackIn);
1ee39b3a 1734
afca20ef 1735 // TRDout TRACK RESOLUTION
92d6d80c 1736 fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
1ee39b3a 1737
1738 return fContainer;
1739}
1740
1741//________________________________________________________
1742Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
1743{
1744 //
1745 // Do the processing
1746 //
1747
92c40e64 1748 Char_t pn[10]; sprintf(pn, "p%03d", fIdxPlot);
1ee39b3a 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);
92d6d80c 1752 AliDebug(4, Form("%s: gm[%s] gs[%s]", pn, g[0]->GetName(), g[1]->GetName()));
1ee39b3a 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);
92d6d80c 1757 if((n=h->GetEntries())<100) continue;
1ee39b3a 1758
1759 h->Fit(f, "QN");
1ee39b3a 1760 Int_t ip = g[0]->GetN();
92d6d80c 1761 AliDebug(4, Form(" x_%d[%f] stat[%d]", ip, x, n));
1ee39b3a 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));
1ee39b3a 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//________________________________________________________
92c40e64 1777Bool_t AliTRDresolution::Process2D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k, Int_t gidx)
1ee39b3a 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;
92d6d80c 1788 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h2->GetName(), h2->GetTitle()));
1789
92c40e64 1790
1ee39b3a 1791 TGraphErrors *g[2];
92c40e64 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;
1ee39b3a 1794
92c40e64 1795 if(!(g[1] = gidx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(gidx))) return kFALSE;
1ee39b3a 1796
1797 return Process(h2, f, k, g);
1798}
1799
1800//________________________________________________________
1801Bool_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;
92d6d80c 1812 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
1ee39b3a 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
92c40e64 1830
1ee39b3a 1831//________________________________________________________
1832Bool_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;
92d6d80c 1843 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
1ee39b3a 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//________________________________________________________
92d6d80c 1876Bool_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//________________________________________________________
1909Bool_t AliTRDresolution::Process3Darray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1ee39b3a 1910{
1911 //
1912 // Do the processing
1913 //
1914
1915 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
92c40e64 1916 //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
1ee39b3a 1917
1918 // retrive containers
92d6d80c 1919 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
1ee39b3a 1920 if(!arr) return kFALSE;
92d6d80c 1921 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
1ee39b3a 1922
92c40e64 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;
1ee39b3a 1926
92d6d80c 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;
1ee39b3a 1931
92d6d80c 1932 if(!(h3 = (TH3S*)a0->At(idx))) return kFALSE;
1933 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h3->GetName(), h3->GetTitle()));
1ee39b3a 1934 TAxis *az = h3->GetZaxis();
92c40e64 1935 for(Int_t iz=1; iz<=az->GetNbins(); iz++, in++){
92c40e64 1936 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
1937 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
1ee39b3a 1938 az->SetRange(iz, iz);
1939 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
1940 }
1941 }
92d6d80c 1942 AliDebug(2, Form("Projections [%d] from [%d]", in, gs->GetEntriesFast()));
1ee39b3a 1943
1944 return kTRUE;
1945}
1946
1947//________________________________________________________
1948Bool_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;
92c40e64 1955
1956 //printf("plotting task[%d] gidx[%d]\n", ip, idx);
1ee39b3a 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;
92c40e64 1961 //printf("gs[%s] gm[%s]\n", gs->GetName(), gm->GetName());
92d6d80c 1962 Int_t n(0);
1963 if(!gs->GetN() || !gm->GetN()) return kFALSE;
1964 gs->Draw("apl");
1965 gm->Draw("pl");
1966
92c40e64 1967 //return kTRUE;
1ee39b3a 1968 // titles look up
1969 Int_t nref = 0;
92c40e64 1970 for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fgNproj[jp];
1ee39b3a 1971 UChar_t jdx = idx<0?0:idx;
92c40e64 1972 for(Int_t jc=0; jc<TMath::Min(jdx,fgNproj[ip]-1); jc++) nref++;
1ee39b3a 1973 const Char_t **at = fgAxTitle[nref];
1974
1ee39b3a 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
1295b37f 1988 TAxis *ax(NULL); TH1 *hf(NULL);
1989 hf = gs->GetHistogram();
1990 hf->SetTitle(at[0]);
1991 ax = hf->GetXaxis();
1ee39b3a 1992 ax->SetRangeUser(bb[0], bb[2]);
1993 ax->SetTitle(at[1]);ax->CenterTitle();
1994
1295b37f 1995 ax = hf->GetYaxis();
1ee39b3a 1996 ax->SetRangeUser(bb[1], bb[3]);
1997 ax->SetTitleOffset(1.1);
1998 ax->SetTitle(at[2]);ax->CenterTitle();
1999
b91fdd71 2000 TGaxis *gax = NULL;
1ee39b3a 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
1ee39b3a 2015//________________________________________________________
92d6d80c 2016Bool_t AliTRDresolution::GetGraphArray(Float_t *bb, ETRDresolutionPlot ip, Int_t idx, Bool_t kLEG, Int_t n, Int_t *sel, const Char_t *explain)
1ee39b3a 2017{
2018 //
2019 // Get the graphs
2020 //
2021
2022 if(!fGraphS || !fGraphM) return kFALSE;
2023
2024 // axis titles look up
92c40e64 2025 Int_t nref(0);
92d6d80c 2026 for(Int_t jp(0); jp<ip; jp++) nref+=fgNproj[jp];
92c40e64 2027 nref+=idx;
1ee39b3a 2028 const Char_t **at = fgAxTitle[nref];
92c40e64 2029
92d6d80c 2030 // build legends if requiered
31c8fa8a 2031 TLegend *legM(NULL), *legS(NULL);
92c40e64 2032 if(kLEG){
1295b37f 2033 legM=new TLegend(.35, .6, .65, .9);
31c8fa8a 2034 legM->SetHeader("Mean");
92d6d80c 2035 legM->SetBorderSize(0);
2036 legM->SetFillStyle(0);
1295b37f 2037 legS=new TLegend(.65, .6, .95, .9);
31c8fa8a 2038 legS->SetHeader("Sigma");
92d6d80c 2039 legS->SetBorderSize(0);
2040 legS->SetFillStyle(0);
1ee39b3a 2041 }
92d6d80c 2042 // build frame
92c40e64 2043 TH1S *h1(NULL);
92d6d80c 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]);
92c40e64 2045 h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2046 h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
1ee39b3a 2047 // axis range
92c40e64 2048 TAxis *ax = h1->GetXaxis();
1295b37f 2049 ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
92c40e64 2050 ax = h1->GetYaxis();
1ee39b3a 2051 ax->SetRangeUser(bb[1], bb[3]);
1295b37f 2052 ax->CenterTitle(); ax->SetTitleOffset(1.4);
92c40e64 2053 h1->Draw();
92d6d80c 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;
1ee39b3a 2066
92c40e64 2067 if((nn=gs->GetN())){
92d6d80c 2068 nPlots++;
31c8fa8a 2069 gs->Draw("pc"); if(legS) legS->AddEntry(gs, gs->GetTitle(), "pl");
92c40e64 2070 gs->Sort(&TGraph::CompareY);
92d6d80c 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]);
92c40e64 2073 gs->Sort(&TGraph::CompareX);
2074 }
2075 if(gm->GetN()){
92d6d80c 2076 nPlots++;
31c8fa8a 2077 gm->Draw("pc");if(legM) legM->AddEntry(gm, gm->GetTitle(), "pl");
92d6d80c 2078 PutTrendValue(Form("%s_%s", fgPerformanceName[kMCtrack], at[0]), gm->GetMean(2));
2079 PutTrendValue(Form("%s_%sRMS", fgPerformanceName[kMCtrack], at[0]), gm->GetRMS(2));
92c40e64 2080 }
2081 }
92d6d80c 2082 if(!nPlots) return kFALSE;
92c40e64 2083 if(kLEG){
92d6d80c 2084 legM->Draw();
2085 legS->Draw();
1ee39b3a 2086 }
1ee39b3a 2087 return kTRUE;
2088}
2089
2090//________________________________________________________
2091void 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//________________________________________________________
2117void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
2118{
2119
2120 fReconstructor->SetRecoParam(r);
2121}