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