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