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