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