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