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