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