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