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