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