]> 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
Fixed size of two arrays (R. Shahoyan)
[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
00a3f7a4 1756//________________________________________________________
1757void AliTRDresolution::MakeSummary()
1758{
1759// Build summary plots
1760
1761 if(!fGraphS || !fGraphM){
1762 AliError("Missing results");
1763 return;
1764 }
1765 Float_t xy[4] = {0., 0., 0., 0.};
1766 Int_t iSumPlot(0);
1767 TCanvas *cOut = new TCanvas(Form("TRDsummary%s_%d", GetName(), iSumPlot), "Cluster & Tracklet Resolution", 1024, 768);
1768 TF1 fg("fg", "gaus", -500., 500.);
1769
1770 if(!HasMCdata()) return;
1771 cOut->Divide(3,2);
1772
1773 TH2 *h2 = new TH2I("h2SF", "", 20, -.2, .2, fgkNresYsegm[fSegmentLevel], -0.5, fgkNresYsegm[fSegmentLevel]-0.5);
1774 h2->GetXaxis()->CenterTitle();
1775 h2->GetYaxis()->CenterTitle();
1776 h2->GetZaxis()->CenterTitle();
1777 TH1 *h1 = new TH1F("h1SF0", "", 120, -200., 1000.);
1778 TGraphErrors *g(NULL);
1779 TAxis *ax(h2->GetXaxis());
1780 Double_t x, y;
1781 TObjArray *a(NULL);
1782
1783 cOut->cd(1);
1784 h2->SetTitle(Form("Cluster R-Phi Resolution;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1785 a=(TObjArray*) ((TObjArray*)fGraphS->At(kMCcluster))->At(0);
1786 for(Int_t iseg(0); iseg<fgkNresYsegm[fSegmentLevel]; iseg++){
1787 g=(TGraphErrors*)a->At(iseg);
1788 for(Int_t in(0); in<g->GetN(); in++){
1789 g->GetPoint(in, x, y);
1790 h2->SetBinContent(ax->FindBin(x), iseg+1, y);
1791 h1->Fill(y);
1792 }
1793 }
1794 Int_t jBinMin(1), jBinMax(h1->GetNbinsX());
1795 for(Int_t ibin(h1->GetMaximumBin()); ibin--;){
1796 if(h1->GetBinContent(ibin)==0){
1797 jBinMin=ibin; break;
1798 }
1799 }
1800 for(Int_t ibin(h1->GetMaximumBin()); ibin++;){
1801 if(h1->GetBinContent(ibin)==0){
1802 jBinMax=ibin; break;
1803 }
1804 }
1805 h2->GetZaxis()->SetRangeUser(h1->GetBinCenter(jBinMin), h1->GetBinCenter(jBinMax));
1806 h2->Draw("col2z");
1807
1808 cOut->cd(2);
1809 h2=(TH2I*)h2->Clone("h2SF"); h2->Reset(); h2->SetTitle(Form("Cluster R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1810 h1=(TH1F*)h1->Clone("h1SF1"); h1->Reset();
1811 //h2->GetZaxis()->SetRangeUser(-100., 100.);
1812 a=(TObjArray*) ((TObjArray*)fGraphM->At(kMCcluster))->At(0);
1813 for(Int_t iseg(0); iseg<fgkNresYsegm[fSegmentLevel]; iseg++){
1814 g=(TGraphErrors*)a->At(iseg);
1815 for(Int_t in(0); in<g->GetN(); in++){
1816 g->GetPoint(in, x, y);
1817 h2->SetBinContent(ax->FindBin(x), iseg+1, y);
1818 h1->Fill(y);
1819 }
1820 }
1821 h1->Fit(&fg, "QN");
1822 Double_t m(fg.GetParameter(1)), s(fg.GetParameter(2));
1823 h2->GetZaxis()->SetRangeUser(m-2.5*s, m+2.5*s);
1824 h2->Draw("col2z");
1825
1826 cOut->cd(3);
1827 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1828 GetGraphArray(xy, kMCcluster, 1, 1);
1829
1830 cOut->cd(4);
1831 h2=(TH2I*)h2->Clone("h2SF"); h2->Reset();
1832 h2->SetTitle(Form("Tracklet R-Phi Resolution;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1833 //h2->GetZaxis()->SetRangeUser(100., 500.);
1834 h1=(TH1F*)h1->Clone("h1SF2"); h1->Reset();
1835 a=(TObjArray*) ((TObjArray*)fGraphS->At(kMCtracklet))->At(0);
1836 for(Int_t iseg(0); iseg<fgkNresYsegm[fSegmentLevel]; iseg++){
1837 g=(TGraphErrors*)a->At(iseg);
1838 for(Int_t in(0); in<g->GetN(); in++){
1839 g->GetPoint(in, x, y);
1840 h2->SetBinContent(ax->FindBin(x), iseg+1, y);
1841 h1->Fill(y);
1842 }
1843 }
1844 jBinMin=1; jBinMax=h1->GetNbinsX();
1845 for(Int_t ibin(h1->GetMaximumBin()); ibin--;){
1846 if(h1->GetBinContent(ibin)==0){
1847 jBinMin=ibin; break;
1848 }
1849 }
1850 for(Int_t ibin(h1->GetMaximumBin()); ibin++;){
1851 if(h1->GetBinContent(ibin)==0){
1852 jBinMax=ibin; break;
1853 }
1854 }
1855 h2->GetZaxis()->SetRangeUser(h1->GetBinCenter(jBinMin), h1->GetBinCenter(jBinMax));
1856 h2->Draw("col2z");
1857
1858 cOut->cd(5);
1859 h2=(TH2I*)h2->Clone("h2SF"); h2->Reset(); h2->SetTitle(Form("Tracklet R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1860 h1=(TH1F*)h1->Clone("h1SF3"); h1->Reset();
1861 a=(TObjArray*) ((TObjArray*)fGraphM->At(kMCtracklet))->At(0);
1862 for(Int_t iseg(0); iseg<fgkNresYsegm[fSegmentLevel]; iseg++){
1863 g=(TGraphErrors*)a->At(iseg);
1864 for(Int_t in(0); in<g->GetN(); in++){
1865 g->GetPoint(in, x, y);
1866 h2->SetBinContent(ax->FindBin(x), iseg+1, y);
1867 h1->Fill(y);
1868 }
1869 }
1870 h1->Fit(&fg, "QN");
1871 m=fg.GetParameter(1); s=fg.GetParameter(2);
1872 h2->GetZaxis()->SetRangeUser(m-2.5*s, m+2.5*s);
1873 h2->Draw("col2z");
1874
1875 cOut->cd(6);
1876 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1877 GetGraphArray(xy, kMCtracklet, 1, 1);
1878
1879 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1880}
1881
1882//________________________________________________________
92c40e64 1883Char_t const *fgParticle[11]={
1884 " p bar", " K -", " #pi -", " #mu -", " e -",
1885 " No PID",
1886 " e +", " #mu +", " #pi +", " K +", " p",
1887};
1295b37f 1888const Color_t fgColorS[11]={
1889kOrange, kOrange-3, kMagenta+1, kViolet, kRed,
31c8fa8a 1890kGray,
1295b37f 1891kRed, kViolet, kMagenta+1, kOrange-3, kOrange
1892};
1893const Color_t fgColorM[11]={
1894kCyan-5, kAzure-4, kBlue-7, kBlue+2, kViolet+10,
1895kBlack,
1896kViolet+10, kBlue+2, kBlue-7, kAzure-4, kCyan-5
31c8fa8a 1897};
1898const Marker_t fgMarker[11]={
189930, 30, 26, 25, 24,
190028,
190120, 21, 22, 29, 29
1902};
1ee39b3a 1903Bool_t AliTRDresolution::PostProcess()
1904{
1905 //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
1906 if (!fContainer) {
1907 AliError("ERROR: list not available");
1908 return kFALSE;
1909 }
b91fdd71 1910 TGraph *gm= NULL, *gs= NULL;
1ee39b3a 1911 if(!fGraphS && !fGraphM){
31c8fa8a 1912 TObjArray *aM(NULL), *aS(NULL);
1ee39b3a 1913 Int_t n = fContainer->GetEntriesFast();
1914 fGraphS = new TObjArray(n); fGraphS->SetOwner();
1915 fGraphM = new TObjArray(n); fGraphM->SetOwner();
92d6d80c 1916 for(Int_t ig(0), nc(0); ig<n; ig++){
92c40e64 1917 fGraphM->AddAt(aM = new TObjArray(fgNproj[ig]), ig);
1918 fGraphS->AddAt(aS = new TObjArray(fgNproj[ig]), ig);
92d6d80c 1919
1920 for(Int_t ic=0; ic<fgNproj[ig]; ic++, nc++){
2589cf64 1921 AliDebug(2, Form("building G[%d] P[%d] N[%d]", ig, ic, fNcomp[nc]));
1922 if(fNcomp[nc]>1){
31c8fa8a 1923 TObjArray *agS(NULL), *agM(NULL);
2589cf64 1924 aS->AddAt(agS = new TObjArray(fNcomp[nc]), ic);
1925 aM->AddAt(agM = new TObjArray(fNcomp[nc]), ic);
1926 for(Int_t is=fNcomp[nc]; is--;){
31c8fa8a 1927 agS->AddAt(gs = new TGraphErrors(), is);
92d6d80c 1928 Int_t is0(is%11), il0(is/11);
31c8fa8a 1929 gs->SetMarkerStyle(fgMarker[is0]);
1295b37f 1930 gs->SetMarkerColor(fgColorS[is0]);
1931 gs->SetLineColor(fgColorS[is0]);
92d6d80c 1932 gs->SetLineStyle(il0);gs->SetLineWidth(2);
2589cf64 1933 gs->SetName(Form("s_%d_%02d_%02d", ig, ic, is));
31c8fa8a 1934
1935 agM->AddAt(gm = new TGraphErrors(), is);
1936 gm->SetMarkerStyle(fgMarker[is0]);
1295b37f 1937 gm->SetMarkerColor(fgColorM[is0]);
1938 gm->SetLineColor(fgColorM[is0]);
92d6d80c 1939 gm->SetLineStyle(il0);gm->SetLineWidth(2);
2589cf64 1940 gm->SetName(Form("m_%d_%02d_%02d", ig, ic, is));
a310e49b 1941 // this is important for labels in the legend
1942 if(ic==0) {
2589cf64 1943 gs->SetTitle(Form("%s %02d", fgkResYsegmName[fSegmentLevel], is%fgkNresYsegm[fSegmentLevel]));
1944 gm->SetTitle(Form("%s %02d", fgkResYsegmName[fSegmentLevel], is%fgkNresYsegm[fSegmentLevel]));
81979445 1945 } else if(ic==1) {
1946 gs->SetTitle(Form("%s Ly[%d]", is%2 ?"z":"y", is/2));
1947 gm->SetTitle(Form("%s Ly[%d]", is%2?"z":"y", is/2));
1948 } else if(ic==2||ic==3) {
1949 gs->SetTitle(Form("%s Ly[%d]", is%2 ?"RC":"no RC", is/2));
1950 gm->SetTitle(Form("%s Ly[%d]", is%2?"RC":"no RC", is/2));
a310e49b 1951 } else if(ic<=7) {
d25116b6 1952 gs->SetTitle(Form("Layer[%d]", is%AliTRDgeometry::kNlayer));
1953 gm->SetTitle(Form("Layer[%d]", is%AliTRDgeometry::kNlayer));
a310e49b 1954 } else {
1955 gs->SetTitle(Form("%s @ ly[%d]", fgParticle[is0], il0));
1956 gm->SetTitle(Form("%s @ ly[%d]", fgParticle[is0], il0));
1957 }
92c40e64 1958 }
1959 } else {
1960 aS->AddAt(gs = new TGraphErrors(), ic);
1961 gs->SetMarkerStyle(23);
1962 gs->SetMarkerColor(kRed);
1963 gs->SetLineColor(kRed);
2589cf64 1964 gs->SetNameTitle(Form("s_%d_%02d", ig, ic), "sigma");
92c40e64 1965
1966 aM->AddAt(gm = ig ? (TGraph*)new TGraphErrors() : (TGraph*)new TGraphAsymmErrors(), ic);
1967 gm->SetLineColor(kBlack);
1968 gm->SetMarkerStyle(7);
1969 gm->SetMarkerColor(kBlack);
2589cf64 1970 gm->SetNameTitle(Form("m_%d_%02d", ig, ic), "mean");
1ee39b3a 1971 }
1ee39b3a 1972 }
1973 }
1974 }
1975
1976/* printf("\n\n\n"); fGraphS->ls();
1977 printf("\n\n\n"); fGraphM->ls();*/
1978
1979
1980 // DEFINE MODELS
1981 // simple gauss
1982 TF1 fg("fGauss", "gaus", -.5, .5);
1983 // Landau for charge resolution
92c40e64 1984 TF1 fch("fClCh", "landau", 0., 1000.);
1985 // Landau for e+- pt resolution
1986 TF1 fpt("fPt", "landau", -0.1, 0.2);
1ee39b3a 1987
1988 //PROCESS EXPERIMENTAL DISTRIBUTIONS
1989 // Charge resolution
1990 //Process3DL(kCharge, 0, &fl);
1991 // Clusters residuals
92d6d80c 1992 Process3D(kCluster, 0, &fg, 1.e4);
81979445 1993 Process3Dlinked(kCluster, 1, &fg);
6558fd69 1994 fNRefFigures = 3;
1ee39b3a 1995 // Tracklet residual/pulls
6558fd69 1996 Process3D(kTrack , 0, &fg, 1.e4);
81979445 1997 Process3Dlinked(kTrack , 1, &fg);
1998 Process3D(kTrack , 2, &fg, 1.e4);
1999 Process3D(kTrack , 3, &fg);
92d6d80c 2000 Process2D(kTrack , 4, &fg, 1.e3);
6558fd69 2001 fNRefFigures = 7;
a310e49b 2002 // TRDin residual/pulls
92d6d80c 2003 Process3D(kTrackIn, 0, &fg, 1.e4);
81979445 2004 Process3Dlinked(kTrackIn, 1, &fg);
2005 Process3D(kTrackIn, 2, &fg, 1.e4);
2006 Process3D(kTrackIn, 3, &fg);
92d6d80c 2007 Process2D(kTrackIn, 4, &fg, 1.e3);
6558fd69 2008 fNRefFigures = 11;
a310e49b 2009 // TRDout residual/pulls
6558fd69 2010 Process3D(kTrackOut, 0, &fg, 1.e3); // scale to fit - see PlotTrackOut
81979445 2011 Process3Dlinked(kTrackOut, 1, &fg);
2012 Process3D(kTrackOut, 2, &fg, 1.e4);
2013 Process3D(kTrackOut, 3, &fg);
a310e49b 2014 Process2D(kTrackOut, 4, &fg, 1.e3);
6558fd69 2015 fNRefFigures = 15;
1ee39b3a 2016
2017 if(!HasMCdata()) return kTRUE;
2018
2019
2020 //PROCESS MC RESIDUAL DISTRIBUTIONS
2021
2022 // CLUSTER Y RESOLUTION/PULLS
92d6d80c 2023 Process3D(kMCcluster, 0, &fg, 1.e4);
81979445 2024 Process3Dlinked(kMCcluster, 1, &fg, 1.);
6558fd69 2025 fNRefFigures = 17;
1ee39b3a 2026
2027 // TRACKLET RESOLUTION/PULLS
92d6d80c 2028 Process3D(kMCtracklet, 0, &fg, 1.e4); // y
81979445 2029 Process3Dlinked(kMCtracklet, 1, &fg, 1.); // y pulls
2030 Process3D(kMCtracklet, 2, &fg, 1.e4); // z
2031 Process3D(kMCtracklet, 3, &fg, 1.); // z pulls
92d6d80c 2032 Process2D(kMCtracklet, 4, &fg, 1.e3); // phi
d25116b6 2033 fNRefFigures = 21;
1ee39b3a 2034
2035 // TRACK RESOLUTION/PULLS
92d6d80c 2036 Process3Darray(kMCtrack, 0, &fg, 1.e4); // y
81979445 2037 Process3DlinkedArray(kMCtrack, 1, &fg); // y PULL
2038 Process3Darray(kMCtrack, 2, &fg, 1.e4); // z
2039 Process3Darray(kMCtrack, 3, &fg); // z PULL
92d6d80c 2040 Process2Darray(kMCtrack, 4, &fg, 1.e3); // phi
2041 Process2Darray(kMCtrack, 5, &fg); // snp PULL
2042 Process2Darray(kMCtrack, 6, &fg, 1.e3); // theta
2043 Process2Darray(kMCtrack, 7, &fg); // tgl PULL
2044 Process3Darray(kMCtrack, 8, &fg, 1.e2); // pt resolution
2045 Process3Darray(kMCtrack, 9, &fg); // 1/pt pulls
2046 Process3Darray(kMCtrack, 10, &fg, 1.e2); // p resolution
d25116b6 2047 fNRefFigures+=16;
92d6d80c 2048
2049 // TRACK TRDin RESOLUTION/PULLS
2050 Process3D(kMCtrackIn, 0, &fg, 1.e4);// y resolution
81979445 2051 Process3Dlinked(kMCtrackIn, 1, &fg); // y pulls
2052 Process3D(kMCtrackIn, 2, &fg, 1.e4);// z resolution
2053 Process3D(kMCtrackIn, 3, &fg); // z pulls
92d6d80c 2054 Process2D(kMCtrackIn, 4, &fg, 1.e3);// phi resolution
2055 Process2D(kMCtrackIn, 5, &fg); // snp pulls
2056 Process2D(kMCtrackIn, 6, &fg, 1.e3);// theta resolution
2057 Process2D(kMCtrackIn, 7, &fg); // tgl pulls
2058 Process3D(kMCtrackIn, 8, &fg, 1.e2);// pt resolution
2059 Process3D(kMCtrackIn, 9, &fg); // 1/pt pulls
2060 Process3D(kMCtrackIn, 10, &fg, 1.e2);// p resolution
d25116b6 2061 fNRefFigures+=8;
92d6d80c 2062
2063 // TRACK TRDout RESOLUTION/PULLS
2064 Process3D(kMCtrackOut, 0, &fg, 1.e4);// y resolution
81979445 2065 Process3Dlinked(kMCtrackOut, 1, &fg); // y pulls
2066 Process3D(kMCtrackOut, 2, &fg, 1.e4);// z resolution
2067 Process3D(kMCtrackOut, 3, &fg); // z pulls
92d6d80c 2068 Process2D(kMCtrackOut, 4, &fg, 1.e3);// phi resolution
2069 Process2D(kMCtrackOut, 5, &fg); // snp pulls
2070 Process2D(kMCtrackOut, 6, &fg, 1.e3);// theta resolution
2071 Process2D(kMCtrackOut, 7, &fg); // tgl pulls
2072 Process3D(kMCtrackOut, 8, &fg, 1.e2);// pt resolution
2073 Process3D(kMCtrackOut, 9, &fg); // 1/pt pulls
2074 Process3D(kMCtrackOut, 10, &fg, 1.e2);// p resolution
d25116b6 2075 fNRefFigures+=8;
1ee39b3a 2076
2077 return kTRUE;
2078}
2079
2080
2081//________________________________________________________
2082void AliTRDresolution::Terminate(Option_t *opt)
2083{
2084 AliTRDrecoTask::Terminate(opt);
2085 if(HasPostProcess()) PostProcess();
2086}
2087
2088//________________________________________________________
2089void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
2090{
2091// Helper function to avoid duplication of code
2092// Make first guesses on the fit parameters
2093
2094 // find the intial parameters of the fit !! (thanks George)
2095 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
2096 Double_t sum = 0.;
2097 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
2098 f->SetParLimits(0, 0., 3.*sum);
2099 f->SetParameter(0, .9*sum);
2100 Double_t rms = h->GetRMS();
2101 f->SetParLimits(1, -rms, rms);
2102 f->SetParameter(1, h->GetMean());
2103
2104 f->SetParLimits(2, 0., 2.*rms);
2105 f->SetParameter(2, rms);
2106 if(f->GetNpar() <= 4) return;
2107
2108 f->SetParLimits(3, 0., sum);
2109 f->SetParameter(3, .1*sum);
2110
2111 f->SetParLimits(4, -.3, .3);
2112 f->SetParameter(4, 0.);
2113
2114 f->SetParLimits(5, 0., 1.e2);
2115 f->SetParameter(5, 2.e-1);
2116}
2117
afca20ef 2118//________________________________________________________
2589cf64 2119TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand)
afca20ef 2120{
3d2a3dff 2121// Build performance histograms for AliTRDcluster.vs TRD track or MC
afca20ef 2122// - y reziduals/pulls
3d2a3dff 2123
2124 TObjArray *arr = new TObjArray(2);
afca20ef 2125 arr->SetName(name); arr->SetOwner();
2126 TH1 *h(NULL); char hname[100], htitle[300];
2127
2128 // tracklet resolution/pull in y direction
3d2a3dff 2129 sprintf(hname, "%s_%s_Y", GetNameId(), name);
2589cf64 2130 sprintf(htitle, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
3d2a3dff 2131 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2589cf64 2132 Int_t nybins=fgkNresYsegm[fSegmentLevel];
2133 if(expand) nybins*=2;
3d2a3dff 2134 h = new TH3S(hname, htitle,
2589cf64 2135 48, -.48, .48, 60, -.15, .15, nybins, -0.5, nybins-0.5);
afca20ef 2136 } else h->Reset();
2137 arr->AddAt(h, 0);
81979445 2138 sprintf(hname, "%s_%s_YZpull", GetNameId(), name);
2589cf64 2139 sprintf(htitle, "YZ pull for \"%s\" @ %s;%s;#Delta y / #sigma_{y};#Delta z / #sigma_{z}", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
81979445 2140 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2589cf64 2141 h = new TH3S(hname, htitle, fgkNresYsegm[fSegmentLevel], -0.5, fgkNresYsegm[fSegmentLevel]-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
afca20ef 2142 } else h->Reset();
2143 arr->AddAt(h, 1);
2144
3d2a3dff 2145 return arr;
2146}
2147
2148//________________________________________________________
2589cf64 2149TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
3d2a3dff 2150{
2151// Build performance histograms for AliExternalTrackParam.vs TRD tracklet
2152// - y reziduals/pulls
2153// - z reziduals/pulls
2154// - phi reziduals
2589cf64 2155 TObjArray *arr = BuildMonitorContainerCluster(name, expand);
3d2a3dff 2156 arr->Expand(5);
2157 TH1 *h(NULL); char hname[100], htitle[300];
2158
afca20ef 2159 // tracklet resolution/pull in z direction
3d2a3dff 2160 sprintf(hname, "%s_%s_Z", GetNameId(), name);
81979445 2161 sprintf(htitle, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm];row cross", GetNameId(), name);
2162 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2163 h = new TH3S(hname, htitle, 50, -1., 1., 100, -1.5, 1.5, 2, -0.5, 1.5);
afca20ef 2164 } else h->Reset();
2165 arr->AddAt(h, 2);
3d2a3dff 2166 sprintf(hname, "%s_%s_Zpull", GetNameId(), name);
81979445 2167 sprintf(htitle, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z / #sigma_{z};row cross", GetNameId(), name);
2168 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2169 h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
dfd7d48b 2170 h->GetZaxis()->SetBinLabel(1, "no RC");
2171 h->GetZaxis()->SetBinLabel(2, "RC");
afca20ef 2172 } else h->Reset();
2173 arr->AddAt(h, 3);
2174
2175 // tracklet to track phi resolution
3d2a3dff 2176 sprintf(hname, "%s_%s_PHI", GetNameId(), name);
2177 sprintf(htitle, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];entries", GetNameId(), name);
afca20ef 2178 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2179 h = new TH2I(hname, htitle, 21, -.33, .33, 100, -.5, .5);
2180 } else h->Reset();
2181 arr->AddAt(h, 4);
2182
2183 return arr;
2184}
2185
2186//________________________________________________________
2187TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
2188{
2189// Build performance histograms for AliExternalTrackParam.vs MC
2190// - y resolution/pulls
2191// - z resolution/pulls
2192// - phi resolution, snp pulls
2193// - theta resolution, tgl pulls
2194// - pt resolution, 1/pt pulls, p resolution
2195
afca20ef 2196 TObjArray *arr = BuildMonitorContainerTracklet(name);
2197 arr->Expand(11);
3d2a3dff 2198 TH1 *h(NULL); char hname[100], htitle[300];
a310e49b 2199 TAxis *ax(NULL);
3d2a3dff 2200
afca20ef 2201 // snp pulls
3d2a3dff 2202 sprintf(hname, "%s_%s_SNPpull", GetNameId(), name);
2203 sprintf(htitle, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp / #sigma_{snp};entries", GetNameId(), name);
afca20ef 2204 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2205 h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
2206 } else h->Reset();
2207 arr->AddAt(h, 5);
2208
2209 // theta resolution
3d2a3dff 2210 sprintf(hname, "%s_%s_THT", GetNameId(), name);
2211 sprintf(htitle, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
afca20ef 2212 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2213 h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
2214 } else h->Reset();
2215 arr->AddAt(h, 6);
2216 // tgl pulls
3d2a3dff 2217 sprintf(hname, "%s_%s_TGLpull", GetNameId(), name);
2218 sprintf(htitle, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl / #sigma_{tgl};entries", GetNameId(), name);
afca20ef 2219 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2220 h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
2221 } else h->Reset();
2222 arr->AddAt(h, 7);
2223
2224 const Int_t kNpt(14);
2225 const Int_t kNdpt(150);
2226 const Int_t kNspc = 2*AliPID::kSPECIES+1;
2227 Float_t Pt=0.1, DPt=-.1, Spc=-5.5;
2228 Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
2229 for(Int_t i=0;i<kNpt+1; i++,Pt=TMath::Exp(i*.15)-1.) binsPt[i]=Pt;
2230 for(Int_t i=0; i<kNspc+1; i++,Spc+=1.) binsSpc[i]=Spc;
2231 for(Int_t i=0; i<kNdpt+1; i++,DPt+=2.e-3) binsDPt[i]=DPt;
2232
2233 // Pt resolution
3d2a3dff 2234 sprintf(hname, "%s_%s_Pt", GetNameId(), name);
2235 sprintf(htitle, "P_{t} res for \"%s\" @ %s;p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
afca20ef 2236 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2237 h = new TH3S(hname, htitle,
2238 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
a310e49b 2239 ax = h->GetZaxis();
2240 for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
afca20ef 2241 } else h->Reset();
2242 arr->AddAt(h, 8);
2243 // 1/Pt pulls
3d2a3dff 2244 sprintf(hname, "%s_%s_1Pt", GetNameId(), name);
2245 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 2246 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2247 h = new TH3S(hname, htitle,
2248 kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
a310e49b 2249 ax = h->GetZaxis();
2250 for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
afca20ef 2251 } else h->Reset();
2252 arr->AddAt(h, 9);
2253 // P resolution
3d2a3dff 2254 sprintf(hname, "%s_%s_P", GetNameId(), name);
2255 sprintf(htitle, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
afca20ef 2256 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2257 h = new TH3S(hname, htitle,
2258 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
a310e49b 2259 ax = h->GetZaxis();
2260 for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
afca20ef 2261 } else h->Reset();
2262 arr->AddAt(h, 10);
2263
2264 return arr;
2265}
2266
2267
1ee39b3a 2268//________________________________________________________
2269TObjArray* AliTRDresolution::Histos()
2270{
2271 //
2272 // Define histograms
2273 //
2274
2275 if(fContainer) return fContainer;
2276
92c40e64 2277 fContainer = new TObjArray(kNviews);
1ee39b3a 2278 //fContainer->SetOwner(kTRUE);
31c8fa8a 2279 TH1 *h(NULL);
2280 TObjArray *arr(NULL);
1ee39b3a 2281
31c8fa8a 2282 // binnings for plots containing momentum or pt
2283 const Int_t kNpt(14), kNphi(48), kNdy(60);
2284 Float_t Phi=-.48, Dy=-.3, Pt=0.1;
2285 Float_t binsPhi[kNphi+1], binsDy[kNdy+1], binsPt[kNpt+1];
2286 for(Int_t i=0; i<kNphi+1; i++,Phi+=.02) binsPhi[i]=Phi;
2287 for(Int_t i=0; i<kNdy+1; i++,Dy+=.01) binsDy[i]=Dy;
2288 for(Int_t i=0;i<kNpt+1; i++,Pt=TMath::Exp(i*.15)-1.) binsPt[i]=Pt;
df0514f6 2289
1ee39b3a 2290 // cluster to track residuals/pulls
92d6d80c 2291 fContainer->AddAt(arr = new TObjArray(2), kCharge);
1ee39b3a 2292 arr->SetName("Charge");
2293 if(!(h = (TH3S*)gROOT->FindObject("hCharge"))){
2294 h = new TH3S("hCharge", "Charge Resolution", 20, 1., 2., 24, 0., 3.6, 100, 0., 500.);
2295 h->GetXaxis()->SetTitle("dx/dx_{0}");
2296 h->GetYaxis()->SetTitle("x_{d} [cm]");
2297 h->GetZaxis()->SetTitle("dq/dx [ADC/cm]");
2298 } else h->Reset();
2299 arr->AddAt(h, 0);
2300
2301 // cluster to track residuals/pulls
3d2a3dff 2302 fContainer->AddAt(BuildMonitorContainerCluster("Cl"), kCluster);
afca20ef 2303 // tracklet to TRD track
2589cf64 2304 fContainer->AddAt(BuildMonitorContainerTracklet("Trk", kTRUE), kTrack);
afca20ef 2305 // tracklet to TRDin
2589cf64 2306 fContainer->AddAt(BuildMonitorContainerTracklet("TrkIN", kTRUE), kTrackIn);
a310e49b 2307 // tracklet to TRDout
2308 fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
1ee39b3a 2309
2310
2311 // Resolution histos
2312 if(!HasMCdata()) return fContainer;
2313
3d2a3dff 2314 // cluster resolution
2315 fContainer->AddAt(BuildMonitorContainerCluster("MCcl"), kMCcluster);
1ee39b3a 2316
3d2a3dff 2317 // tracklet resolution
2318 fContainer->AddAt(BuildMonitorContainerTracklet("MCtracklet"), kMCtracklet);
1ee39b3a 2319
3d2a3dff 2320 // track resolution
92d6d80c 2321 fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
3d2a3dff 2322 arr->SetName("MCtrk");
2323 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
31c8fa8a 2324
afca20ef 2325 // TRDin TRACK RESOLUTION
92d6d80c 2326 fContainer->AddAt(BuildMonitorContainerTrack("MCtrkIN"), kMCtrackIn);
1ee39b3a 2327
afca20ef 2328 // TRDout TRACK RESOLUTION
92d6d80c 2329 fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
1ee39b3a 2330
2331 return fContainer;
2332}
2333
2589cf64 2334//________________________________________________________
fe1d1beb 2335Bool_t AliTRDresolution::Load(const Char_t *file, const Char_t *dir)
2589cf64 2336{
2337// Custom load function. Used to guess the segmentation level of the data.
2338
fe1d1beb 2339 if(!AliTRDrecoTask::Load(file, dir)) return kFALSE;
2589cf64 2340
2341 // look for cluster residual plot - always available
2342 TH3S* h3((TH3S*)((TObjArray*)fContainer->At(kClToTrk))->At(0));
2343 Int_t segmentation(h3->GetNbinsZ()/2);
2344 if(segmentation==fgkNresYsegm[0]){ // default segmentation. Nothing to do
2345 return kTRUE;
2346 } else if(segmentation==fgkNresYsegm[1]){ // stack segmentation.
2347 SetSegmentationLevel(1);
2348 } else if(segmentation==fgkNresYsegm[2]){ // detector segmentation.
2349 SetSegmentationLevel(2);
2350 } else {
2351 AliError(Form("Unknown segmentation [%d].", h3->GetNbinsZ()));
2352 return kFALSE;
2353 }
2354
2355 AliDebug(2, Form("Segmentation set to level \"%s\"", fgkResYsegmName[fSegmentLevel]));
2356 return kTRUE;
2357}
2358
2359
1ee39b3a 2360//________________________________________________________
2361Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
2362{
2363 //
2364 // Do the processing
2365 //
2366
92c40e64 2367 Char_t pn[10]; sprintf(pn, "p%03d", fIdxPlot);
1ee39b3a 2368 Int_t n = 0;
2369 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
2370 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
2589cf64 2371 if(Int_t(h2->GetEntries())){
2372 AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
2373 } else {
2374 AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
2375 fIdxPlot++;
2376 return kTRUE;
2377 }
1ee39b3a 2378
dfd7d48b 2379 const Int_t kINTEGRAL=1;
2380 for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
2381 Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
2382 Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
2383 TH1D *h = h2->ProjectionY(pn, abin, bbin);
2589cf64 2384 if((n=(Int_t)h->GetEntries())<100){
2385 AliDebug(4, Form(" x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
2386 continue;
2387 }
1ee39b3a 2388 h->Fit(f, "QN");
1ee39b3a 2389 Int_t ip = g[0]->GetN();
2589cf64 2390 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 2391 g[0]->SetPoint(ip, x, k*f->GetParameter(1));
2392 g[0]->SetPointError(ip, 0., k*f->GetParError(1));
2393 g[1]->SetPoint(ip, x, k*f->GetParameter(2));
2394 g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1ee39b3a 2395/*
2396 g[0]->SetPoint(ip, x, k*h->GetMean());
2397 g[0]->SetPointError(ip, 0., k*h->GetMeanError());
2398 g[1]->SetPoint(ip, x, k*h->GetRMS());
2399 g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
2400 }
2401 fIdxPlot++;
2402 return kTRUE;
2403}
2404
2405//________________________________________________________
92c40e64 2406Bool_t AliTRDresolution::Process2D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k, Int_t gidx)
1ee39b3a 2407{
2408 //
2409 // Do the processing
2410 //
2411
2412 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2413
2414 // retrive containers
a310e49b 2415 TH2I *h2(NULL);
2416 if(idx<0){
2417 if(!(h2= (TH2I*)(fContainer->At(plot)))) return kFALSE;
2418 } else{
2419 TObjArray *a0(NULL);
2420 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2421 if(!(h2=(TH2I*)a0->At(idx))) return kFALSE;
2422 }
2589cf64 2423 if(Int_t(h2->GetEntries())){
2424 AliDebug(2, Form("p[%d] idx[%d] : h[%s] %s", plot, idx, h2->GetName(), h2->GetTitle()));
2425 } else {
2426 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2427 return kFALSE;
2428 }
92c40e64 2429
1ee39b3a 2430 TGraphErrors *g[2];
92c40e64 2431 if(gidx<0) gidx=idx;
2432 if(!(g[0] = gidx<0 ? (TGraphErrors*)fGraphM->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(gidx))) return kFALSE;
1ee39b3a 2433
92c40e64 2434 if(!(g[1] = gidx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(gidx))) return kFALSE;
1ee39b3a 2435
2436 return Process(h2, f, k, g);
2437}
2438
2439//________________________________________________________
2440Bool_t AliTRDresolution::Process3D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2441{
2442 //
2443 // Do the processing
2444 //
2445
2446 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2447
2448 // retrive containers
a310e49b 2449 TH3S *h3(NULL);
2450 if(idx<0){
2451 if(!(h3= (TH3S*)(fContainer->At(plot)))) return kFALSE;
2452 } else{
2453 TObjArray *a0(NULL);
2454 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2455 if(!(h3=(TH3S*)a0->At(idx))) return kFALSE;
2456 }
2589cf64 2457 if(Int_t(h3->GetEntries())){
2458 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2459 } else {
2460 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2461 return kFALSE;
2462 }
1ee39b3a 2463
2464 TObjArray *gm, *gs;
2465 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2466 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2467 TGraphErrors *g[2];
2468
2469 TAxis *az = h3->GetZaxis();
dfd7d48b 2470 for(Int_t iz(0); iz<gm->GetEntriesFast(); iz++){
2471 if(!(g[0] = (TGraphErrors*)gm->At(iz))) return kFALSE;
2472 if(!(g[1] = (TGraphErrors*)gs->At(iz))) return kFALSE;
2473 az->SetRange(iz+1, iz+1);
1ee39b3a 2474 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2475 }
2476
2477 return kTRUE;
2478}
2479
92c40e64 2480
81979445 2481//________________________________________________________
2482Bool_t AliTRDresolution::Process3Dlinked(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2483{
2484 //
2485 // Do the processing
2486 //
2487
2488 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2489
2490 // retrive containers
2491 TH3S *h3(NULL);
2492 if(idx<0){
2493 if(!(h3= (TH3S*)(fContainer->At(plot)))) return kFALSE;
2494 } else{
2495 TObjArray *a0(NULL);
2496 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2497 if(!(h3=(TH3S*)a0->At(idx))) return kFALSE;
2498 }
2589cf64 2499 if(Int_t(h3->GetEntries())){
2500 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2501 } else {
2502 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2503 return kFALSE;
2504 }
81979445 2505
2506 TObjArray *gm, *gs;
2507 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2508 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2509 TGraphErrors *g[2];
2510
2511 if(!(g[0] = (TGraphErrors*)gm->At(0))) return kFALSE;
2512 if(!(g[1] = (TGraphErrors*)gs->At(0))) return kFALSE;
2513 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2514
2515 if(!(g[0] = (TGraphErrors*)gm->At(1))) return kFALSE;
2516 if(!(g[1] = (TGraphErrors*)gs->At(1))) return kFALSE;
2517 if(!Process((TH2*)h3->Project3D("zx"), f, k, g)) return kFALSE;
2518
2519 return kTRUE;
2520}
2521
2522
1ee39b3a 2523//________________________________________________________
2524Bool_t AliTRDresolution::Process3DL(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2525{
2526 //
2527 // Do the processing
2528 //
2529
2530 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2531
2532 // retrive containers
2533 TH3S *h3 = (TH3S*)((TObjArray*)fContainer->At(plot))->At(idx);
2534 if(!h3) return kFALSE;
92d6d80c 2535 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
1ee39b3a 2536
2537 TGraphAsymmErrors *gm;
2538 TGraphErrors *gs;
2539 if(!(gm = (TGraphAsymmErrors*)((TObjArray*)fGraphM->At(plot))->At(0))) return kFALSE;
2540 if(!(gs = (TGraphErrors*)((TObjArray*)fGraphS->At(plot)))) return kFALSE;
2541
2542 Float_t x, r, mpv, xM, xm;
2543 TAxis *ay = h3->GetYaxis();
2544 for(Int_t iy=1; iy<=h3->GetNbinsY(); iy++){
2545 ay->SetRange(iy, iy);
2546 x = ay->GetBinCenter(iy);
2547 TH2F *h2=(TH2F*)h3->Project3D("zx");
2548 TAxis *ax = h2->GetXaxis();
2549 for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
2550 r = ax->GetBinCenter(ix);
2551 TH1D *h1 = h2->ProjectionY("py", ix, ix);
2552 if(h1->Integral()<50) continue;
2553 h1->Fit(f, "QN");
2554
2555 GetLandauMpvFwhm(f, mpv, xm, xM);
2556 Int_t ip = gm->GetN();
2557 gm->SetPoint(ip, x, k*mpv);
2558 gm->SetPointError(ip, 0., 0., k*xm, k*xM);
2559 gs->SetPoint(ip, r, k*(xM-xm)/mpv);
2560 gs->SetPointError(ip, 0., 0.);
2561 }
2562 }
2563
2564 return kTRUE;
2565}
2566
2567//________________________________________________________
92d6d80c 2568Bool_t AliTRDresolution::Process2Darray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2569{
2570 //
2571 // Do the processing
2572 //
2573
2574 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2575
2576 // retrive containers
2577 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2578 if(!arr) return kFALSE;
2579 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2580
2581 TObjArray *gm, *gs;
2582 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2583 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2584
2585 TGraphErrors *g[2]; TH2I *h2(NULL); TObjArray *a0(NULL);
2586 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2587 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2588
2589 if(!(h2 = (TH2I*)a0->At(idx))) return kFALSE;
2589cf64 2590 if(Int_t(h2->GetEntries())){
2591 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h2->GetName(), h2->GetTitle()));
2592 } else {
2593 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2594 continue;
2595 }
92d6d80c 2596
2597 if(!(g[0] = (TGraphErrors*)gm->At(ia))) return kFALSE;
2598 if(!(g[1] = (TGraphErrors*)gs->At(ia))) return kFALSE;
2599 if(!Process(h2, f, k, g)) return kFALSE;
2600 }
2601
2602 return kTRUE;
2603}
2604
2605//________________________________________________________
2606Bool_t AliTRDresolution::Process3Darray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1ee39b3a 2607{
2608 //
2609 // Do the processing
2610 //
2611
2612 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
92c40e64 2613 //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
1ee39b3a 2614
2615 // retrive containers
92d6d80c 2616 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
1ee39b3a 2617 if(!arr) return kFALSE;
92d6d80c 2618 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
1ee39b3a 2619
92c40e64 2620 TObjArray *gm, *gs;
2621 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2622 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
1ee39b3a 2623
92d6d80c 2624 TGraphErrors *g[2]; TH3S *h3(NULL); TObjArray *a0(NULL);
2625 Int_t in(0);
2626 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2627 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
1ee39b3a 2628
92d6d80c 2629 if(!(h3 = (TH3S*)a0->At(idx))) return kFALSE;
2589cf64 2630 if(Int_t(h3->GetEntries())){
2631 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h3->GetName(), h3->GetTitle()));
2632 } else {
2633 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2634 continue;
2635 }
1ee39b3a 2636 TAxis *az = h3->GetZaxis();
92c40e64 2637 for(Int_t iz=1; iz<=az->GetNbins(); iz++, in++){
dfd7d48b 2638 if(in >= gm->GetEntriesFast()) break;
92c40e64 2639 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2640 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
1ee39b3a 2641 az->SetRange(iz, iz);
2642 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2643 }
2644 }
92d6d80c 2645 AliDebug(2, Form("Projections [%d] from [%d]", in, gs->GetEntriesFast()));
1ee39b3a 2646
2647 return kTRUE;
2648}
2649
81979445 2650//________________________________________________________
2651Bool_t AliTRDresolution::Process3DlinkedArray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2652{
2653 //
2654 // Do the processing
2655 //
2656
2657 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2658 //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
2659
2660 // retrive containers
2661 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2662 if(!arr) return kFALSE;
2663 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2664
2665 TObjArray *gm, *gs;
2666 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2667 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2668
2669 TGraphErrors *g[2]; TH3S *h3(NULL); TObjArray *a0(NULL);
2670 Int_t in(0);
2671 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2672 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2673 if(!(h3 = (TH3S*)a0->At(idx))) return kFALSE;
2589cf64 2674 if(Int_t(h3->GetEntries())){
2675 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h3->GetName(), h3->GetTitle()));
2676 } else {
2677 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2678 continue;
2679 }
81979445 2680 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2681 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2682 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2683 in++;
2684
2685 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2686 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2687 if(!Process((TH2*)h3->Project3D("zx"), f, k, g)) return kFALSE;
2688 in++;
2689 }
2690 AliDebug(2, Form("Projections [%d] from [%d]", in, gs->GetEntriesFast()));
2691
2692 return kTRUE;
2693}
2694
1ee39b3a 2695//________________________________________________________
a310e49b 2696Bool_t AliTRDresolution::GetGraph(Float_t *bb, ETRDresolutionPlot ip, Int_t idx, Bool_t kLEG, const Char_t *explain)
1ee39b3a 2697{
2698 //
2699 // Get the graphs
2700 //
2701
2702 if(!fGraphS || !fGraphM) return kFALSE;
a310e49b 2703 // axis titles look up
1ee39b3a 2704 Int_t nref = 0;
92c40e64 2705 for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fgNproj[jp];
1ee39b3a 2706 UChar_t jdx = idx<0?0:idx;
92c40e64 2707 for(Int_t jc=0; jc<TMath::Min(jdx,fgNproj[ip]-1); jc++) nref++;
2589cf64 2708 Char_t **at = fAxTitle[nref];
1ee39b3a 2709
a310e49b 2710 // build legends if requiered
2711 TLegend *leg(NULL);
2712 if(kLEG){
2713 leg=new TLegend(.65, .6, .95, .9);
2714 leg->SetBorderSize(0);
2715 leg->SetFillStyle(0);
2716 }
2717 // build frame
2718 TH1S *h1(NULL);
2719 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]);
2720 h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2721 h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
2722 // axis range
2723 TAxis *ax = h1->GetXaxis();
2724 ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
2725 ax = h1->GetYaxis();
2726 ax->SetRangeUser(bb[1], bb[3]);
2727 ax->CenterTitle(); ax->SetTitleOffset(1.4);
2728 h1->Draw();
2729 // bounding box
2730 TBox *b = new TBox(-.15, bb[1], .15, bb[3]);
2731 b->SetFillStyle(3002);b->SetLineColor(0);
2732 b->SetFillColor(ip<=Int_t(kMCcluster)?kGreen:kBlue);
2733 b->Draw();
2734
2735 TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphM->At(ip)))->At(idx);
2736 if(!gm) return kFALSE;
2737 TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphS->At(ip)))->At(idx);
2738 if(!gs) return kFALSE;
2739
2740 Int_t n(0), nPlots(0);
1ee39b3a 2741 if((n=gm->GetN())) {
a310e49b 2742 nPlots++;
2743 gm->Draw("pl"); if(leg) leg->AddEntry(gm, gm->GetTitle(), "pl");
1ee39b3a 2744 PutTrendValue(Form("%s_%s", fgPerformanceName[ip], at[0]), gm->GetMean(2));
2745 PutTrendValue(Form("%s_%sRMS", fgPerformanceName[ip], at[0]), gm->GetRMS(2));
2746 }
2747
2748 if((n=gs->GetN())){
a310e49b 2749 nPlots++;
2750 gs->Draw("pl"); if(leg) leg->AddEntry(gs, gs->GetTitle(), "pl");
1ee39b3a 2751 gs->Sort(&TGraph::CompareY);
2752 PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[ip], at[0]), gs->GetY()[0]);
2753 PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[ip], at[0]), gs->GetY()[n-1]);
2754 gs->Sort(&TGraph::CompareX);
2755 }
a310e49b 2756 if(!nPlots) return kFALSE;
2757 if(leg) leg->Draw();
1ee39b3a 2758 return kTRUE;
2759}
2760
1ee39b3a 2761//________________________________________________________
92d6d80c 2762Bool_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 2763{
2764 //
2765 // Get the graphs
2766 //
2767
2768 if(!fGraphS || !fGraphM) return kFALSE;
2769
2770 // axis titles look up
92c40e64 2771 Int_t nref(0);
92d6d80c 2772 for(Int_t jp(0); jp<ip; jp++) nref+=fgNproj[jp];
92c40e64 2773 nref+=idx;
2589cf64 2774 Char_t **at = fAxTitle[nref];
92c40e64 2775
92d6d80c 2776 // build legends if requiered
31c8fa8a 2777 TLegend *legM(NULL), *legS(NULL);
92c40e64 2778 if(kLEG){
1295b37f 2779 legM=new TLegend(.35, .6, .65, .9);
31c8fa8a 2780 legM->SetHeader("Mean");
92d6d80c 2781 legM->SetBorderSize(0);
2782 legM->SetFillStyle(0);
1295b37f 2783 legS=new TLegend(.65, .6, .95, .9);
31c8fa8a 2784 legS->SetHeader("Sigma");
92d6d80c 2785 legS->SetBorderSize(0);
2786 legS->SetFillStyle(0);
1ee39b3a 2787 }
92d6d80c 2788 // build frame
92c40e64 2789 TH1S *h1(NULL);
92d6d80c 2790 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 2791 h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2792 h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
1ee39b3a 2793 // axis range
92c40e64 2794 TAxis *ax = h1->GetXaxis();
1295b37f 2795 ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
92c40e64 2796 ax = h1->GetYaxis();
1ee39b3a 2797 ax->SetRangeUser(bb[1], bb[3]);
1295b37f 2798 ax->CenterTitle(); ax->SetTitleOffset(1.4);
92c40e64 2799 h1->Draw();
92d6d80c 2800
2801 TGraphErrors *gm(NULL), *gs(NULL);
2802 TObjArray *a0(NULL), *a1(NULL);
2803 a0 = (TObjArray*)((TObjArray*)fGraphM->At(ip))->At(idx);
2804 a1 = (TObjArray*)((TObjArray*)fGraphS->At(ip))->At(idx);
2805 if(!n) n=a0->GetEntriesFast();
2806 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'));
2807 Int_t nn(0), nPlots(0);
2808 for(Int_t is(0), is0(0); is<n; is++){
2809 is0 = sel ? sel[is] : is;
2810 if(!(gs = (TGraphErrors*)a1->At(is0))) return kFALSE;
2811 if(!(gm = (TGraphErrors*)a0->At(is0))) return kFALSE;
1ee39b3a 2812
92c40e64 2813 if((nn=gs->GetN())){
92d6d80c 2814 nPlots++;
a310e49b 2815 gs->Draw("pc");
2816 if(legS){
2817 //printf("LegEntry %s [%s]%s\n", at[0], gs->GetName(), gs->GetTitle());
2818 legS->AddEntry(gs, gs->GetTitle(), "pl");
2819 }
92c40e64 2820 gs->Sort(&TGraph::CompareY);
92d6d80c 2821 PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[kMCtrack], at[0]), gs->GetY()[0]);
2822 PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[kMCtrack], at[0]), gs->GetY()[nn-1]);
92c40e64 2823 gs->Sort(&TGraph::CompareX);
2824 }
2825 if(gm->GetN()){
92d6d80c 2826 nPlots++;
a310e49b 2827 gm->Draw("pc");
2828 if(legM){
2829 //printf("LegEntry %s [%s]%s\n", at[0], gm->GetName(), gm->GetTitle());
2830 legM->AddEntry(gm, gm->GetTitle(), "pl");
2831 }
92d6d80c 2832 PutTrendValue(Form("%s_%s", fgPerformanceName[kMCtrack], at[0]), gm->GetMean(2));
2833 PutTrendValue(Form("%s_%sRMS", fgPerformanceName[kMCtrack], at[0]), gm->GetRMS(2));
92c40e64 2834 }
2835 }
92d6d80c 2836 if(!nPlots) return kFALSE;
92c40e64 2837 if(kLEG){
92d6d80c 2838 legM->Draw();
2839 legS->Draw();
1ee39b3a 2840 }
1ee39b3a 2841 return kTRUE;
2842}
2843
2844//________________________________________________________
2845void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
2846{
2847 //
2848 // Get the most probable value and the full width half mean
2849 // of a Landau distribution
2850 //
2851
2852 const Float_t dx = 1.;
2853 mpv = f->GetParameter(1);
2854 Float_t fx, max = f->Eval(mpv);
2855
2856 xm = mpv - dx;
2857 while((fx = f->Eval(xm))>.5*max){
2858 if(fx>max){
2859 max = fx;
2860 mpv = xm;
2861 }
2862 xm -= dx;
2863 }
2864
2865 xM += 2*(mpv - xm);
2866 while((fx = f->Eval(xM))>.5*max) xM += dx;
2867}
2868
2869
2870//________________________________________________________
2871void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
2872{
2873
2874 fReconstructor->SetRecoParam(r);
2875}
2589cf64 2876
2877
2878//________________________________________________________
2879void AliTRDresolution::SetSegmentationLevel(Int_t l)
2880{
2881// Setting the segmentation level to "l"
2882 fSegmentLevel = l;
2883
2884 UShort_t const lNcomp[kNprojs] = {
2885 1, 1, //2,
2886 fgkNresYsegm[fSegmentLevel], 2, //2,
2887 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2888 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2889 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2890 // MC
2891 fgkNresYsegm[fSegmentLevel], 2, //2,
2892 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2893 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, 1, 1, 1, 11, 11, 11, //11
2894 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, 1, 1, 1, 11, 11, 11, //11
2895 6*fgkNresYsegm[fSegmentLevel], 6*2, 6*2, 6*2, 6, 6, 6, 6, 6*11, 6*11, 6*11 //11
2896 };
2897 memcpy(fNcomp, lNcomp, kNprojs*sizeof(UShort_t));
2898
2899 Char_t const *lAxTitle[kNprojs][4] = {
2900 // Charge
2901 {"Impv", "x [cm]", "I_{mpv}", "x/x_{0}"}
2902 ,{"dI/Impv", "x/x_{0}", "#delta I/I_{mpv}", "x[cm]"}
2903 // Clusters to Kalman
2904 ,{"Cluster2Track residuals", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
2905 ,{"Cluster2Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2906 // TRD tracklet to Kalman fit
2907 ,{"Tracklet2Track Y residuals", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
2908 ,{"Tracklet2Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2909 ,{"Tracklet2Track Z residuals", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
2910 ,{"Tracklet2Track Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
2911 ,{"Tracklet2Track Phi residuals", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
2912 // TRDin 2 first TRD tracklet
2913 ,{"Tracklet2Track Y residuals @ TRDin", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
2914 ,{"Tracklet2Track YZ pulls @ TRDin", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2915 ,{"Tracklet2Track Z residuals @ TRDin", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
2916 ,{"Tracklet2Track Z pulls @ TRDin", "tg(#theta)", "z", "#sigma_{z}"}
2917 ,{"Tracklet2Track Phi residuals @ TRDin", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
2918 // TRDout 2 first TRD tracklet
2919 ,{"Tracklet2Track Y residuals @ TRDout", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
2920 ,{"Tracklet2Track YZ pulls @ TRDout", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2921 ,{"Tracklet2Track Z residuals @ TRDout", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
2922 ,{"Tracklet2Track Z pulls @ TRDout", "tg(#theta)", "z", "#sigma_{z}"}
2923 ,{"Tracklet2Track Phi residuals @ TRDout", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
2924 // MC cluster
2925 ,{"MC Cluster Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
2926 ,{"MC Cluster YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2927 // MC tracklet
2928 ,{"MC Tracklet Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
2929 ,{"MC Tracklet YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2930 ,{"MC Tracklet Z resolution", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
2931 ,{"MC Tracklet Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
2932 ,{"MC Tracklet Phi resolution", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
2933 // MC track TRDin
2934 ,{"MC Y resolution @ TRDin", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
2935 ,{"MC YZ pulls @ TRDin", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2936 ,{"MC Z resolution @ TRDin", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
2937 ,{"MC Z pulls @ TRDin", "tg(#theta)", "z", "#sigma_{z}"}
2938 ,{"MC #Phi resolution @ TRDin", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
2939 ,{"MC SNP pulls @ TRDin", "tg(#phi)", "SNP", "#sigma_{snp}"}
2940 ,{"MC #Theta resolution @ TRDin", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
2941 ,{"MC TGL pulls @ TRDin", "tg(#theta)", "TGL", "#sigma_{tgl}"}
2942 ,{"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}) [%]"}
2943 ,{"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}"}
2944 ,{"MC P resolution @ TRDin", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
2945 // MC track TRDout
2946 ,{"MC Y resolution @ TRDout", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
2947 ,{"MC YZ pulls @ TRDout", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2948 ,{"MC Z resolution @ TRDout", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
2949 ,{"MC Z pulls @ TRDout", "tg(#theta)", "z", "#sigma_{z}"}
2950 ,{"MC #Phi resolution @ TRDout", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
2951 ,{"MC SNP pulls @ TRDout", "tg(#phi)", "SNP", "#sigma_{snp}"}
2952 ,{"MC #Theta resolution @ TRDout", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
2953 ,{"MC TGL pulls @ TRDout", "tg(#theta)", "TGL", "#sigma_{tgl}"}
2954 ,{"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}) [%]"}
2955 ,{"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}"}
2956 ,{"MC P resolution @ TRDout", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
2957 // MC track in TRD
2958 ,{"MC Track Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
2959 ,{"MC Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2960 ,{"MC Track Z resolution", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
2961 ,{"MC Track Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
2962 ,{"MC Track #Phi resolution", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
2963 ,{"MC Track SNP pulls", "tg(#phi)", "SNP", "#sigma_{snp}"}
2964 ,{"MC Track #Theta resolution", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
2965 ,{"MC Track TGL pulls", "tg(#theta)", "TGL", "#sigma_{tgl}"}
2966 ,{"MC P_{t} resolution", "p_{t} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "#sigma(#Deltap_{t}/p_{t}^{MC}) [%]"}
2967 ,{"MC 1/P_{t} pulls", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC} - 1/p_{t}^{MC}", "#sigma_{1/p_{t}}"}
2968 ,{"MC P resolution", "p [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "#sigma(#Deltap/p^{MC}) [%]"}
2969 };
2970 memcpy(fAxTitle, lAxTitle, 4*kNprojs*sizeof(Char_t*));
2971}