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