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