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