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