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