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