]> 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
remove left over print statement
[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
2513
1ee39b3a 2514//________________________________________________________
2515Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
2516{
2517 //
2518 // Do the processing
2519 //
2520
92c40e64 2521 Char_t pn[10]; sprintf(pn, "p%03d", fIdxPlot);
1ee39b3a 2522 Int_t n = 0;
2523 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
2524 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
2589cf64 2525 if(Int_t(h2->GetEntries())){
2526 AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
2527 } else {
2528 AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
2529 fIdxPlot++;
2530 return kTRUE;
2531 }
1ee39b3a 2532
dfd7d48b 2533 const Int_t kINTEGRAL=1;
2534 for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
2535 Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
2536 Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
2537 TH1D *h = h2->ProjectionY(pn, abin, bbin);
068e2c00 2538 if((n=(Int_t)h->GetEntries())<150){
2589cf64 2539 AliDebug(4, Form(" x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
2540 continue;
2541 }
1ee39b3a 2542 h->Fit(f, "QN");
1ee39b3a 2543 Int_t ip = g[0]->GetN();
2589cf64 2544 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 2545 g[0]->SetPoint(ip, x, k*f->GetParameter(1));
2546 g[0]->SetPointError(ip, 0., k*f->GetParError(1));
2547 g[1]->SetPoint(ip, x, k*f->GetParameter(2));
2548 g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1ee39b3a 2549/*
2550 g[0]->SetPoint(ip, x, k*h->GetMean());
2551 g[0]->SetPointError(ip, 0., k*h->GetMeanError());
2552 g[1]->SetPoint(ip, x, k*h->GetRMS());
2553 g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
2554 }
2555 fIdxPlot++;
2556 return kTRUE;
2557}
2558
2559//________________________________________________________
92c40e64 2560Bool_t AliTRDresolution::Process2D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k, Int_t gidx)
1ee39b3a 2561{
2562 //
2563 // Do the processing
2564 //
2565
2566 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2567
2568 // retrive containers
a310e49b 2569 TH2I *h2(NULL);
2570 if(idx<0){
2571 if(!(h2= (TH2I*)(fContainer->At(plot)))) return kFALSE;
2572 } else{
2573 TObjArray *a0(NULL);
2574 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2575 if(!(h2=(TH2I*)a0->At(idx))) return kFALSE;
2576 }
2589cf64 2577 if(Int_t(h2->GetEntries())){
2578 AliDebug(2, Form("p[%d] idx[%d] : h[%s] %s", plot, idx, h2->GetName(), h2->GetTitle()));
2579 } else {
2580 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2581 return kFALSE;
2582 }
92c40e64 2583
1ee39b3a 2584 TGraphErrors *g[2];
92c40e64 2585 if(gidx<0) gidx=idx;
2586 if(!(g[0] = gidx<0 ? (TGraphErrors*)fGraphM->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(gidx))) return kFALSE;
1ee39b3a 2587
92c40e64 2588 if(!(g[1] = gidx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(gidx))) return kFALSE;
1ee39b3a 2589
2590 return Process(h2, f, k, g);
2591}
2592
2593//________________________________________________________
2594Bool_t AliTRDresolution::Process3D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2595{
2596 //
2597 // Do the processing
2598 //
2599
2600 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2601
2602 // retrive containers
a310e49b 2603 TH3S *h3(NULL);
2604 if(idx<0){
2605 if(!(h3= (TH3S*)(fContainer->At(plot)))) return kFALSE;
2606 } else{
2607 TObjArray *a0(NULL);
2608 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2609 if(!(h3=(TH3S*)a0->At(idx))) return kFALSE;
2610 }
2589cf64 2611 if(Int_t(h3->GetEntries())){
2612 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2613 } else {
2614 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2615 return kFALSE;
2616 }
1ee39b3a 2617
2618 TObjArray *gm, *gs;
2619 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2620 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2621 TGraphErrors *g[2];
2622
2623 TAxis *az = h3->GetZaxis();
dfd7d48b 2624 for(Int_t iz(0); iz<gm->GetEntriesFast(); iz++){
2625 if(!(g[0] = (TGraphErrors*)gm->At(iz))) return kFALSE;
2626 if(!(g[1] = (TGraphErrors*)gs->At(iz))) return kFALSE;
2627 az->SetRange(iz+1, iz+1);
1ee39b3a 2628 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2629 }
2630
2631 return kTRUE;
2632}
2633
92c40e64 2634
81979445 2635//________________________________________________________
2636Bool_t AliTRDresolution::Process3Dlinked(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2637{
2638 //
2639 // Do the processing
2640 //
2641
2642 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2643
2644 // retrive containers
2645 TH3S *h3(NULL);
2646 if(idx<0){
2647 if(!(h3= (TH3S*)(fContainer->At(plot)))) return kFALSE;
2648 } else{
2649 TObjArray *a0(NULL);
2650 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2651 if(!(h3=(TH3S*)a0->At(idx))) return kFALSE;
2652 }
2589cf64 2653 if(Int_t(h3->GetEntries())){
2654 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2655 } else {
2656 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2657 return kFALSE;
2658 }
81979445 2659
2660 TObjArray *gm, *gs;
2661 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2662 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2663 TGraphErrors *g[2];
2664
2665 if(!(g[0] = (TGraphErrors*)gm->At(0))) return kFALSE;
2666 if(!(g[1] = (TGraphErrors*)gs->At(0))) return kFALSE;
2667 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2668
2669 if(!(g[0] = (TGraphErrors*)gm->At(1))) return kFALSE;
2670 if(!(g[1] = (TGraphErrors*)gs->At(1))) return kFALSE;
2671 if(!Process((TH2*)h3->Project3D("zx"), f, k, g)) return kFALSE;
2672
2673 return kTRUE;
2674}
2675
2676
1ee39b3a 2677//________________________________________________________
2678Bool_t AliTRDresolution::Process3DL(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2679{
2680 //
2681 // Do the processing
2682 //
2683
2684 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2685
2686 // retrive containers
2687 TH3S *h3 = (TH3S*)((TObjArray*)fContainer->At(plot))->At(idx);
2688 if(!h3) return kFALSE;
92d6d80c 2689 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
1ee39b3a 2690
2691 TGraphAsymmErrors *gm;
2692 TGraphErrors *gs;
2693 if(!(gm = (TGraphAsymmErrors*)((TObjArray*)fGraphM->At(plot))->At(0))) return kFALSE;
2694 if(!(gs = (TGraphErrors*)((TObjArray*)fGraphS->At(plot)))) return kFALSE;
2695
2696 Float_t x, r, mpv, xM, xm;
2697 TAxis *ay = h3->GetYaxis();
2698 for(Int_t iy=1; iy<=h3->GetNbinsY(); iy++){
2699 ay->SetRange(iy, iy);
2700 x = ay->GetBinCenter(iy);
2701 TH2F *h2=(TH2F*)h3->Project3D("zx");
2702 TAxis *ax = h2->GetXaxis();
2703 for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
2704 r = ax->GetBinCenter(ix);
2705 TH1D *h1 = h2->ProjectionY("py", ix, ix);
2706 if(h1->Integral()<50) continue;
2707 h1->Fit(f, "QN");
2708
2709 GetLandauMpvFwhm(f, mpv, xm, xM);
2710 Int_t ip = gm->GetN();
2711 gm->SetPoint(ip, x, k*mpv);
2712 gm->SetPointError(ip, 0., 0., k*xm, k*xM);
2713 gs->SetPoint(ip, r, k*(xM-xm)/mpv);
2714 gs->SetPointError(ip, 0., 0.);
2715 }
2716 }
2717
2718 return kTRUE;
2719}
2720
2721//________________________________________________________
92d6d80c 2722Bool_t AliTRDresolution::Process2Darray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2723{
2724 //
2725 // Do the processing
2726 //
2727
2728 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2729
2730 // retrive containers
2731 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2732 if(!arr) return kFALSE;
2733 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2734
2735 TObjArray *gm, *gs;
2736 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2737 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2738
2739 TGraphErrors *g[2]; TH2I *h2(NULL); TObjArray *a0(NULL);
2740 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2741 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2742
2743 if(!(h2 = (TH2I*)a0->At(idx))) return kFALSE;
2589cf64 2744 if(Int_t(h2->GetEntries())){
2745 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h2->GetName(), h2->GetTitle()));
2746 } else {
2747 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2748 continue;
2749 }
92d6d80c 2750
2751 if(!(g[0] = (TGraphErrors*)gm->At(ia))) return kFALSE;
2752 if(!(g[1] = (TGraphErrors*)gs->At(ia))) return kFALSE;
2753 if(!Process(h2, f, k, g)) return kFALSE;
2754 }
2755
2756 return kTRUE;
2757}
2758
2759//________________________________________________________
2760Bool_t AliTRDresolution::Process3Darray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1ee39b3a 2761{
2762 //
2763 // Do the processing
2764 //
2765
2766 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
92c40e64 2767 //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
1ee39b3a 2768
2769 // retrive containers
92d6d80c 2770 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
1ee39b3a 2771 if(!arr) return kFALSE;
92d6d80c 2772 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
1ee39b3a 2773
92c40e64 2774 TObjArray *gm, *gs;
2775 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2776 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
1ee39b3a 2777
92d6d80c 2778 TGraphErrors *g[2]; TH3S *h3(NULL); TObjArray *a0(NULL);
2779 Int_t in(0);
2780 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2781 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
1ee39b3a 2782
92d6d80c 2783 if(!(h3 = (TH3S*)a0->At(idx))) return kFALSE;
2589cf64 2784 if(Int_t(h3->GetEntries())){
2785 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h3->GetName(), h3->GetTitle()));
2786 } else {
2787 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2788 continue;
2789 }
1ee39b3a 2790 TAxis *az = h3->GetZaxis();
92c40e64 2791 for(Int_t iz=1; iz<=az->GetNbins(); iz++, in++){
dfd7d48b 2792 if(in >= gm->GetEntriesFast()) break;
92c40e64 2793 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2794 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
1ee39b3a 2795 az->SetRange(iz, iz);
2796 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2797 }
2798 }
92d6d80c 2799 AliDebug(2, Form("Projections [%d] from [%d]", in, gs->GetEntriesFast()));
1ee39b3a 2800
2801 return kTRUE;
2802}
2803
81979445 2804//________________________________________________________
2805Bool_t AliTRDresolution::Process3DlinkedArray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2806{
2807 //
2808 // Do the processing
2809 //
2810
2811 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2812 //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
2813
2814 // retrive containers
2815 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2816 if(!arr) return kFALSE;
2817 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2818
2819 TObjArray *gm, *gs;
2820 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2821 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2822
2823 TGraphErrors *g[2]; TH3S *h3(NULL); TObjArray *a0(NULL);
2824 Int_t in(0);
2825 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2826 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2827 if(!(h3 = (TH3S*)a0->At(idx))) return kFALSE;
2589cf64 2828 if(Int_t(h3->GetEntries())){
2829 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h3->GetName(), h3->GetTitle()));
2830 } else {
2831 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2832 continue;
2833 }
81979445 2834 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2835 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2836 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2837 in++;
2838
2839 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2840 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2841 if(!Process((TH2*)h3->Project3D("zx"), f, k, g)) return kFALSE;
2842 in++;
2843 }
2844 AliDebug(2, Form("Projections [%d] from [%d]", in, gs->GetEntriesFast()));
2845
2846 return kTRUE;
2847}
2848
1ee39b3a 2849//________________________________________________________
a310e49b 2850Bool_t AliTRDresolution::GetGraph(Float_t *bb, ETRDresolutionPlot ip, Int_t idx, Bool_t kLEG, const Char_t *explain)
1ee39b3a 2851{
2852 //
2853 // Get the graphs
2854 //
2855
2856 if(!fGraphS || !fGraphM) return kFALSE;
a310e49b 2857 // axis titles look up
1ee39b3a 2858 Int_t nref = 0;
92c40e64 2859 for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fgNproj[jp];
1ee39b3a 2860 UChar_t jdx = idx<0?0:idx;
92c40e64 2861 for(Int_t jc=0; jc<TMath::Min(jdx,fgNproj[ip]-1); jc++) nref++;
2589cf64 2862 Char_t **at = fAxTitle[nref];
1ee39b3a 2863
a310e49b 2864 // build legends if requiered
2865 TLegend *leg(NULL);
2866 if(kLEG){
2867 leg=new TLegend(.65, .6, .95, .9);
2868 leg->SetBorderSize(0);
2869 leg->SetFillStyle(0);
2870 }
2871 // build frame
2872 TH1S *h1(NULL);
2873 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]);
2874 h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2875 h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
2876 // axis range
2877 TAxis *ax = h1->GetXaxis();
2878 ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
2879 ax = h1->GetYaxis();
2880 ax->SetRangeUser(bb[1], bb[3]);
2881 ax->CenterTitle(); ax->SetTitleOffset(1.4);
2882 h1->Draw();
2883 // bounding box
2884 TBox *b = new TBox(-.15, bb[1], .15, bb[3]);
2885 b->SetFillStyle(3002);b->SetLineColor(0);
2886 b->SetFillColor(ip<=Int_t(kMCcluster)?kGreen:kBlue);
2887 b->Draw();
2888
2889 TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphM->At(ip)))->At(idx);
2890 if(!gm) return kFALSE;
2891 TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphS->At(ip)))->At(idx);
2892 if(!gs) return kFALSE;
2893
2894 Int_t n(0), nPlots(0);
1ee39b3a 2895 if((n=gm->GetN())) {
a310e49b 2896 nPlots++;
2897 gm->Draw("pl"); if(leg) leg->AddEntry(gm, gm->GetTitle(), "pl");
1ee39b3a 2898 PutTrendValue(Form("%s_%s", fgPerformanceName[ip], at[0]), gm->GetMean(2));
2899 PutTrendValue(Form("%s_%sRMS", fgPerformanceName[ip], at[0]), gm->GetRMS(2));
2900 }
2901
2902 if((n=gs->GetN())){
a310e49b 2903 nPlots++;
2904 gs->Draw("pl"); if(leg) leg->AddEntry(gs, gs->GetTitle(), "pl");
1ee39b3a 2905 gs->Sort(&TGraph::CompareY);
2906 PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[ip], at[0]), gs->GetY()[0]);
2907 PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[ip], at[0]), gs->GetY()[n-1]);
2908 gs->Sort(&TGraph::CompareX);
2909 }
a310e49b 2910 if(!nPlots) return kFALSE;
2911 if(leg) leg->Draw();
1ee39b3a 2912 return kTRUE;
2913}
2914
1ee39b3a 2915//________________________________________________________
92d6d80c 2916Bool_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 2917{
2918 //
2919 // Get the graphs
2920 //
2921
2922 if(!fGraphS || !fGraphM) return kFALSE;
2923
2924 // axis titles look up
92c40e64 2925 Int_t nref(0);
92d6d80c 2926 for(Int_t jp(0); jp<ip; jp++) nref+=fgNproj[jp];
92c40e64 2927 nref+=idx;
2589cf64 2928 Char_t **at = fAxTitle[nref];
92c40e64 2929
92d6d80c 2930 // build legends if requiered
31c8fa8a 2931 TLegend *legM(NULL), *legS(NULL);
92c40e64 2932 if(kLEG){
1295b37f 2933 legM=new TLegend(.35, .6, .65, .9);
31c8fa8a 2934 legM->SetHeader("Mean");
92d6d80c 2935 legM->SetBorderSize(0);
2936 legM->SetFillStyle(0);
1295b37f 2937 legS=new TLegend(.65, .6, .95, .9);
31c8fa8a 2938 legS->SetHeader("Sigma");
92d6d80c 2939 legS->SetBorderSize(0);
2940 legS->SetFillStyle(0);
1ee39b3a 2941 }
92d6d80c 2942 // build frame
92c40e64 2943 TH1S *h1(NULL);
92d6d80c 2944 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 2945 h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2946 h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
1ee39b3a 2947 // axis range
92c40e64 2948 TAxis *ax = h1->GetXaxis();
1295b37f 2949 ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
92c40e64 2950 ax = h1->GetYaxis();
1ee39b3a 2951 ax->SetRangeUser(bb[1], bb[3]);
1295b37f 2952 ax->CenterTitle(); ax->SetTitleOffset(1.4);
92c40e64 2953 h1->Draw();
92d6d80c 2954
2955 TGraphErrors *gm(NULL), *gs(NULL);
2956 TObjArray *a0(NULL), *a1(NULL);
2957 a0 = (TObjArray*)((TObjArray*)fGraphM->At(ip))->At(idx);
2958 a1 = (TObjArray*)((TObjArray*)fGraphS->At(ip))->At(idx);
2959 if(!n) n=a0->GetEntriesFast();
2960 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'));
2961 Int_t nn(0), nPlots(0);
2962 for(Int_t is(0), is0(0); is<n; is++){
2963 is0 = sel ? sel[is] : is;
2964 if(!(gs = (TGraphErrors*)a1->At(is0))) return kFALSE;
2965 if(!(gm = (TGraphErrors*)a0->At(is0))) return kFALSE;
1ee39b3a 2966
92c40e64 2967 if((nn=gs->GetN())){
92d6d80c 2968 nPlots++;
a310e49b 2969 gs->Draw("pc");
2970 if(legS){
2971 //printf("LegEntry %s [%s]%s\n", at[0], gs->GetName(), gs->GetTitle());
2972 legS->AddEntry(gs, gs->GetTitle(), "pl");
2973 }
92c40e64 2974 gs->Sort(&TGraph::CompareY);
92d6d80c 2975 PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[kMCtrack], at[0]), gs->GetY()[0]);
2976 PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[kMCtrack], at[0]), gs->GetY()[nn-1]);
92c40e64 2977 gs->Sort(&TGraph::CompareX);
2978 }
2979 if(gm->GetN()){
92d6d80c 2980 nPlots++;
a310e49b 2981 gm->Draw("pc");
2982 if(legM){
2983 //printf("LegEntry %s [%s]%s\n", at[0], gm->GetName(), gm->GetTitle());
2984 legM->AddEntry(gm, gm->GetTitle(), "pl");
2985 }
92d6d80c 2986 PutTrendValue(Form("%s_%s", fgPerformanceName[kMCtrack], at[0]), gm->GetMean(2));
2987 PutTrendValue(Form("%s_%sRMS", fgPerformanceName[kMCtrack], at[0]), gm->GetRMS(2));
92c40e64 2988 }
2989 }
92d6d80c 2990 if(!nPlots) return kFALSE;
92c40e64 2991 if(kLEG){
92d6d80c 2992 legM->Draw();
2993 legS->Draw();
1ee39b3a 2994 }
1ee39b3a 2995 return kTRUE;
2996}
2997
08bdd9d1 2998//____________________________________________________________________
2999Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
3000{
3001//
3002// Fit track with a staight line using the "np" clusters stored in the array "points".
3003// The following particularities are stored in the clusters from points:
3004// 1. pad tilt as cluster charge
3005// 2. pad row cross or vertex constrain as fake cluster with cluster type 1
3006// The parameters of the straight line fit are stored in the array "param" in the following order :
3007// param[0] - x0 reference radial position
3008// param[1] - y0 reference r-phi position @ x0
3009// param[2] - z0 reference z position @ x0
3010// param[3] - slope dy/dx
3011// param[4] - slope dz/dx
3012//
3013// Attention :
3014// Function should be used to refit tracks for B=0T
3015//
3016
3017 if(np<40){
3018 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].", np);
3019 return kFALSE;
3020 }
3021 TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
3022
3023 Double_t x0(0.);
3024 for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
3025 x0/=Float_t(np);
3026
00a56108 3027 Double_t x, y, z, dx, tilt(0.);
08bdd9d1 3028 for(Int_t ip(0); ip<np; ip++){
3029 x = points[ip].GetX(); z = points[ip].GetZ();
3030 dx = x - x0;
3031 zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
3032 }
3033 if(zfitter.Eval() != 0) return kFALSE;
3034
3035 Double_t z0 = zfitter.GetParameter(0);
3036 Double_t dzdx = zfitter.GetParameter(1);
3037 for(Int_t ip(0); ip<np; ip++){
3038 if(points[ip].GetClusterType()) continue;
3039 x = points[ip].GetX();
3040 dx = x - x0;
3041 y = points[ip].GetY();
3042 z = points[ip].GetZ();
3043 tilt = points[ip].GetCharge();
3044 y -= tilt*(-dzdx*dx + z - z0);
3045 Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
3046 yfitter.AddPoint(&dx, y, 1.);
3047 }
3048 if(yfitter.Eval() != 0) return kFALSE;
3049 Double_t y0 = yfitter.GetParameter(0);
3050 Double_t dydx = yfitter.GetParameter(1);
3051
3052 param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
3053 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);
3054 return kTRUE;
3055}
3056
5f7f1c07 3057//____________________________________________________________________
3058Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, AliTrackPoint *points, const Float_t param[10], Float_t par[3])
3059{
3060//
3061// Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
3062// See function FitTrack for the data stored in the "clusters" array
3063
3064// The parameters of the straight line fit are stored in the array "param" in the following order :
3065// par[0] - x0 reference radial position
3066// par[1] - y0 reference r-phi position @ x0
3067// par[2] - slope dy/dx
3068//
3069// Attention :
3070// Function should be used to refit tracks for B=0T
3071//
3072
3073 TLinearFitter yfitter(2, "pol1");
3074
3075 // grep data for tracklet
3076 Double_t x0(0.), x[60], y[60], dy[60];
3077 Int_t nly(0);
3078 for(Int_t ip(0); ip<np; ip++){
3079 if(points[ip].GetClusterType()) continue;
3080 if(points[ip].GetVolumeID() != ly) continue;
3081 Float_t xt(points[ip].GetX())
3082 ,yt(param[1] + param[3] * (xt - param[0]));
3083 x[nly] = xt;
3084 y[nly] = points[ip].GetY();
3085 dy[nly]= y[nly]-yt;
3086 x0 += xt;
3087 nly++;
3088 }
3089 if(nly<10){
3090 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].", nly);
3091 return kFALSE;
3092 }
3093 // set radial reference for fit
3094 x0 /= Float_t(nly);
3095
3096 // find tracklet core
3097 Double_t mean(0.), sig(1.e3);
3098 AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
3099
3100 // simple cluster error parameterization
3101 Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
3102
3103 // fit tracklet core
3104 for(Int_t jly(0); jly<nly; jly++){
3105 if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
3106 Double_t dx(x[jly]-x0);
3107 yfitter.AddPoint(&dx, y[jly], 1.);
3108 }
3109 if(yfitter.Eval() != 0) return kFALSE;
3110 par[0] = x0;
3111 par[1] = yfitter.GetParameter(0);
3112 par[2] = yfitter.GetParameter(1);
3113 return kTRUE;
3114}
3115
00a56108 3116//____________________________________________________________________
3117Bool_t AliTRDresolution::UseTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
3118{
3119//
3120// Global selection mechanism of tracksbased on cluster to fit residuals
3121// The parameters are the same as used ni function FitTrack().
3122
3123 const Float_t kS(0.6), kM(0.2);
3124 TH1S h("h1", "", 100, -5.*kS, 5.*kS);
3125 Float_t dy, dz, s, m;
3126 for(Int_t ip(0); ip<np; ip++){
3127 if(points[ip].GetClusterType()) continue;
3128 Float_t x0(points[ip].GetX())
3129 ,y0(param[1] + param[3] * (x0 - param[0]))
3130 ,z0(param[2] + param[4] * (x0 - param[0]));
3131 dy=points[ip].GetY() - y0; h.Fill(dy);
3132 dz=points[ip].GetZ() - z0;
3133 }
3134 TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
3135 fg.SetParameter(1, 0.);
3136 fg.SetParameter(2, 2.e-2);
3137 h.Fit(&fg, "QN");
3138 m=fg.GetParameter(1); s=fg.GetParameter(2);
3139 if(s>kS || TMath::Abs(m)>kM) return kFALSE;
3140 return kTRUE;
3141}
3142
1ee39b3a 3143//________________________________________________________
3144void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
3145{
3146 //
3147 // Get the most probable value and the full width half mean
3148 // of a Landau distribution
3149 //
3150
3151 const Float_t dx = 1.;
3152 mpv = f->GetParameter(1);
3153 Float_t fx, max = f->Eval(mpv);
3154
3155 xm = mpv - dx;
3156 while((fx = f->Eval(xm))>.5*max){
3157 if(fx>max){
3158 max = fx;
3159 mpv = xm;
3160 }
3161 xm -= dx;
3162 }
3163
3164 xM += 2*(mpv - xm);
3165 while((fx = f->Eval(xM))>.5*max) xM += dx;
3166}
3167
3168
2589cf64 3169//________________________________________________________
3170void AliTRDresolution::SetSegmentationLevel(Int_t l)
3171{
3172// Setting the segmentation level to "l"
3173 fSegmentLevel = l;
3174
3175 UShort_t const lNcomp[kNprojs] = {
3176 1, 1, //2,
3177 fgkNresYsegm[fSegmentLevel], 2, //2,
3178 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
3179 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
3180 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
3181 // MC
3182 fgkNresYsegm[fSegmentLevel], 2, //2,
3183 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
3184 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, 1, 1, 1, 11, 11, 11, //11
3185 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, 1, 1, 1, 11, 11, 11, //11
3186 6*fgkNresYsegm[fSegmentLevel], 6*2, 6*2, 6*2, 6, 6, 6, 6, 6*11, 6*11, 6*11 //11
3187 };
3188 memcpy(fNcomp, lNcomp, kNprojs*sizeof(UShort_t));
3189
3190 Char_t const *lAxTitle[kNprojs][4] = {
3191 // Charge
3192 {"Impv", "x [cm]", "I_{mpv}", "x/x_{0}"}
3193 ,{"dI/Impv", "x/x_{0}", "#delta I/I_{mpv}", "x[cm]"}
3194 // Clusters to Kalman
3195 ,{"Cluster2Track residuals", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3196 ,{"Cluster2Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3197 // TRD tracklet to Kalman fit
3198 ,{"Tracklet2Track Y residuals", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3199 ,{"Tracklet2Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3200 ,{"Tracklet2Track Z residuals", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3201 ,{"Tracklet2Track Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
3202 ,{"Tracklet2Track Phi residuals", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3203 // TRDin 2 first TRD tracklet
3204 ,{"Tracklet2Track Y residuals @ TRDin", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3205 ,{"Tracklet2Track YZ pulls @ TRDin", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3206 ,{"Tracklet2Track Z residuals @ TRDin", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3207 ,{"Tracklet2Track Z pulls @ TRDin", "tg(#theta)", "z", "#sigma_{z}"}
3208 ,{"Tracklet2Track Phi residuals @ TRDin", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3209 // TRDout 2 first TRD tracklet
3210 ,{"Tracklet2Track Y residuals @ TRDout", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3211 ,{"Tracklet2Track YZ pulls @ TRDout", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3212 ,{"Tracklet2Track Z residuals @ TRDout", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3213 ,{"Tracklet2Track Z pulls @ TRDout", "tg(#theta)", "z", "#sigma_{z}"}
3214 ,{"Tracklet2Track Phi residuals @ TRDout", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3215 // MC cluster
3216 ,{"MC Cluster Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3217 ,{"MC Cluster YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3218 // MC tracklet
3219 ,{"MC Tracklet Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
3220 ,{"MC Tracklet YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3221 ,{"MC Tracklet Z resolution", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3222 ,{"MC Tracklet Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
3223 ,{"MC Tracklet Phi resolution", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3224 // MC track TRDin
3225 ,{"MC Y resolution @ TRDin", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
3226 ,{"MC YZ pulls @ TRDin", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3227 ,{"MC Z resolution @ TRDin", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3228 ,{"MC Z pulls @ TRDin", "tg(#theta)", "z", "#sigma_{z}"}
3229 ,{"MC #Phi resolution @ TRDin", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3230 ,{"MC SNP pulls @ TRDin", "tg(#phi)", "SNP", "#sigma_{snp}"}
3231 ,{"MC #Theta resolution @ TRDin", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3232 ,{"MC TGL pulls @ TRDin", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3233 ,{"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}) [%]"}
3234 ,{"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}"}
3235 ,{"MC P resolution @ TRDin", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
3236 // MC track TRDout
3237 ,{"MC Y resolution @ TRDout", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
3238 ,{"MC YZ pulls @ TRDout", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3239 ,{"MC Z resolution @ TRDout", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3240 ,{"MC Z pulls @ TRDout", "tg(#theta)", "z", "#sigma_{z}"}
3241 ,{"MC #Phi resolution @ TRDout", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3242 ,{"MC SNP pulls @ TRDout", "tg(#phi)", "SNP", "#sigma_{snp}"}
3243 ,{"MC #Theta resolution @ TRDout", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3244 ,{"MC TGL pulls @ TRDout", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3245 ,{"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}) [%]"}
3246 ,{"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}"}
3247 ,{"MC P resolution @ TRDout", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
3248 // MC track in TRD
3249 ,{"MC Track Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3250 ,{"MC Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3251 ,{"MC Track Z resolution", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3252 ,{"MC Track Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
3253 ,{"MC Track #Phi resolution", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3254 ,{"MC Track SNP pulls", "tg(#phi)", "SNP", "#sigma_{snp}"}
3255 ,{"MC Track #Theta resolution", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3256 ,{"MC Track TGL pulls", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3257 ,{"MC P_{t} resolution", "p_{t} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "#sigma(#Deltap_{t}/p_{t}^{MC}) [%]"}
3258 ,{"MC 1/P_{t} pulls", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC} - 1/p_{t}^{MC}", "#sigma_{1/p_{t}}"}
3259 ,{"MC P resolution", "p [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "#sigma(#Deltap/p^{MC}) [%]"}
3260 };
3261 memcpy(fAxTitle, lAxTitle, 4*kNprojs*sizeof(Char_t*));
3262}