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