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