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