1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliTRDresolution.cxx 27496 2008-07-22 08:35:45Z cblume $ */
18 ////////////////////////////////////////////////////////////////////////////
20 // TRD tracking resolution //
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
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
32 // For stand alone usage please refer to the following example:
34 // gSystem->Load("libANALYSIS.so");
35 // gSystem->Load("libTRDqaRec.so");
36 // AliTRDresolution *res = new AliTRDresolution();
37 // //res->SetMCdata();
38 // //res->SetVerbose();
39 // //res->SetVisual();
41 // if(!res->PostProcess()) return;
42 // res->GetRefFigure(0);
46 // Alexandru Bercuci <A.Bercuci@gsi.de> //
47 // Markus Fasel <M.Fasel@gsi.de> //
49 ////////////////////////////////////////////////////////////////////////////
52 #include <TObjArray.h>
60 #include <TGraphErrors.h>
61 #include <TGraphAsymmErrors.h>
65 #include <TTreeStream.h>
66 #include <TGeoManager.h>
70 #include "AliTRDresolution.h"
71 #include "AliTRDgeometry.h"
72 #include "AliTRDpadPlane.h"
73 #include "AliTRDcluster.h"
74 #include "AliTRDseedV1.h"
75 #include "AliTRDtrackV1.h"
76 #include "AliTRDReconstructor.h"
77 #include "AliTRDrecoParam.h"
79 #include "info/AliTRDclusterInfo.h"
81 ClassImp(AliTRDresolution)
83 UChar_t const AliTRDresolution::fgNElements[kNhistos] = {
87 Char_t const * AliTRDresolution::fgPerformanceName[kNhistos] = {
98 Char_t const *AliTRDresolution::fgAxTitle[46][4] = {
100 {"Impv", "x [cm]", "I_{mpv}", "x/x_{0}"}
101 ,{"dI/Impv", "x/x_{0}", "#delta I/I_{mpv}", "x[cm]"}
102 // Clusters to Kalman
103 ,{"Pos", "tg(#phi)", "#mu_{y}^{cl} [#mum]", "#sigma_{y}^{cl} [#mum]"}
104 ,{"Pulls", "tg(#phi)", "PULL: #mu_{y}^{cl}", "PULL: #sigma_{y}^{cl}"}
105 // TRD tracklet to Kalman fit
106 ,{"PosY", "tg(#phi)", "#mu_{y}^{trklt} [#mum]", "#sigma_{y}^{trklt} [#mum]"}
107 ,{"PullsY", "tg(#phi)", "PULL: #mu_{y}^{trklt}", "PULL: #sigma_{y}^{trklt}"}
108 ,{"PosZ", "tg(#theta)", "#mu_{z}^{trklt} [#mum]", "#sigma_{z}^{trklt} [#mum]"}
109 ,{"PullsZ", "tg(#theta)", "PULL: #mu_{z}^{trklt}", "PULL: #sigma_{z}^{trklt}"}
110 ,{"Phi", "tg(#phi)", "#mu_{#phi}^{trklt} [mrad]", "#sigma_{#phi}^{trklt} [mrad]"}
111 // TPC track 2 first TRD tracklet
112 ,{"PosY", "tg(#phi)", "#mu_{y}^{TPC trklt} [#mum]", "#sigma_{y}^{TPC trklt} [#mum]"}
113 ,{"PullsY", "tg(#phi)", "PULL: #mu_{y}^{TPC trklt}", "PULL: #sigma_{y}^{TPC trklt}"}
114 ,{"PosZ", "tg(#theta)", "#mu_{z}^{TPC trklt} [#mum]", "#sigma_{z}^{TPC trklt} [#mum]"}
115 ,{"PullsZ", "tg(#theta)", "PULL: #mu_{z}^{TPC trklt}", "PULL: #sigma_{z}^{TPC trklt}"}
116 ,{"Phi", "tg(#phi)", "#mu_{#phi}^{TPC trklt} [mrad]", "#sigma_{#phi}^{TPC trklt} [mrad]"}
118 ,{"Pos", "tg(#phi)", "MC: #mu_{y}^{cl} [#mum]", "MC: #sigma_{y}^{cl} [#mum]"}
119 ,{"Pulls", "tg(#phi)", "MC PULL: #mu_{y}^{cl}", "MC PULL: #sigma_{y}^{cl}"}
121 ,{"PosY", "tg(#phi)", "MC: #mu_{y}^{trklt} [#mum]", "MC: #sigma_{y}^{trklt} [#mum]"}
122 ,{"PullsY", "tg(#phi)", "MC PULL: #mu_{y}^{trklt}", "MC PULL: #sigma_{y}^{trklt}"}
123 ,{"PosZ", "tg(#theta)", "MC: #mu_{z}^{trklt} [#mum]", "MC: #sigma_{z}^{trklt} [#mum]"}
124 ,{"PullsZ", "tg(#theta)", "MC PULL: #mu_{z}^{trklt}", "MC PULL: #sigma_{z}^{trklt}"}
125 ,{"Phi", "tg(#phi)", "MC: #mu_{#phi}^{trklt} [mrad]", "MC: #sigma_{#phi}^{trklt} [mrad]"}
127 ,{"PosY", "tg(#phi)", "MC: #mu_{y}^{TPC} [#mum]", "MC: #sigma_{y}^{TPC} [#mum]"}
128 ,{"PullsY", "tg(#phi)", "MC PULL: #mu_{y}^{TPC}", "MC PULL: #sigma_{y}^{TPC}"}
129 ,{"PosZ", "tg(#theta)", "MC: #mu_{z}^{TPC} [#mum]", "MC: #sigma_{z}^{TPC} [#mum]"}
130 ,{"PullsZ", "tg(#theta)", "MC PULL: #mu_{z}^{TPC}", "MC PULL: #sigma_{z}^{TPC}"}
131 ,{"Phi", "tg(#phi)", "MC: #mu_{#phi}^{TPC} [mrad]", "MC: #sigma_{#phi}^{TPC} [mrad]"}
132 ,{"PullsSNP", "tg(#phi)", "MC PULL: #mu_{snp}^{TPC}", "MC PULL: #sigma_{snp}^{TPC}"}
133 ,{"Theta", "tg(#theta)", "MC: #mu_{#theta}^{TPC} [mrad]", "MC: #sigma_{#theta}^{TPC} [mrad]"}
134 ,{"PullsTGL", "tg(#theta)", "MC PULL: #mu_{tgl}^{TPC}", "MC PULL: #sigma_{tgl}^{TPC}"}
135 ,{"Pt", "p_{t}^{MC} [GeV/c]", "MC: #mu^{TPC}(#Deltap_{t}/p_{t}^{MC}) [%]", "MC: #sigma^{TPC}(#Deltap_{t}/p_{t}^{MC}) [%]"}
136 ,{"Pulls1Pt", "1/p_{t}^{MC} [c/GeV]", "MC PULL: #mu_{1/p_{t}}^{TPC}", "MC PULL: #sigma_{1/p_{t}}^{TPC}"}
137 ,{"P", "p^{MC} [GeV/c]", "MC: #mu^{TPC}(#Deltap/p^{MC}) [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
138 ,{"PullsP", "p^{MC} [GeV/c]", "MC PULL: #mu^{TPC}(#Deltap/#sigma_{p})", "MC PULL: #sigma^{TPC}(#Deltap/#sigma_{p})"}
140 ,{"PosZ", "tg(#theta)", "MC: #mu_{z}^{TOF} [#mum]", "MC: #sigma_{z}^{TOF} [#mum]"}
141 ,{"PullsZ", "tg(#theta)", "MC PULL: #mu_{z}^{TOF}", "MC PULL: #sigma_{z}^{TOF}"}
143 ,{"PosY", "tg(#phi)", "MC: #mu_{y}^{Trk} [#mum]", "MC: #sigma_{y}^{Trk} [#mum]"}
144 ,{"PullsY", "tg(#phi)", "MC PULL: #mu_{y}^{Trk}", "MC PULL: #sigma_{y}^{Trk}"}
145 ,{"PosZ", "tg(#theta)", "MC: #mu_{z}^{Trk} [#mum]", "MC: #sigma_{z}^{Trk} [#mum]"}
146 ,{"PullsZ", "tg(#theta)", "MC PULL: #mu_{z}^{Trk}", "MC PULL: #sigma_{z}^{Trk}"}
147 ,{"Phi", "tg(#phi)", "MC: #mu_{#phi}^{Trk} [mrad]", "MC: #sigma_{#phi}^{Trk} [mrad]"}
148 ,{"PullsSNP", "tg(#phi)", "MC PULL: #mu_{snp}^{Trk}", "MC PULL: #sigma_{snp}^{Trk}"}
149 ,{"Theta", "tg(#theta)", "MC: #mu_{#theta}^{Trk} [mrad]", "MC: #sigma_{#theta}^{Trk} [mrad]"}
150 ,{"PullsTGL", "tg(#theta)", "MC PULL: #mu_{tgl}^{Trk}", "MC PULL: #sigma_{tgl}^{Trk}"}
151 ,{"Pt", "p_{t}^{MC} [GeV/c]", "MC: #mu^{Trk}(#Deltap_{t}/p_{t}^{MC}) [%]", "MC: #sigma^{Trk}(#Deltap_{t}/p_{t}^{MC}) [%]"}
152 ,{"Pulls1Pt", "1/p_{t}^{MC} [c/GeV]", "MC PULL: #mu_{1/p_{t}}^{Trk}", "MC PULL: #sigma_{1/p_{t}}^{Trk}"}
153 ,{"P", "p^{MC} [GeV/c]", "MC: #mu^{Trk}(#Deltap/p^{MC}) [%]", "MC: #sigma^{Trk}(#Deltap/p^{MC}) [%]"}
156 //________________________________________________________
157 AliTRDresolution::AliTRDresolution()
158 :AliTRDrecoTask("resolution", "Spatial and momentum TRD resolution checker")
161 ,fReconstructor(NULL)
171 // Default constructor
174 fReconstructor = new AliTRDReconstructor();
175 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
176 fGeo = new AliTRDgeometry();
180 DefineOutput(1, TObjArray::Class()); // cluster2track
181 DefineOutput(2, TObjArray::Class()); // tracklet2track
182 DefineOutput(3, TObjArray::Class()); // cluster2mc
183 DefineOutput(4, TObjArray::Class()); // tracklet2mc
186 //________________________________________________________
187 AliTRDresolution::~AliTRDresolution()
193 if(fGraphS){fGraphS->Delete(); delete fGraphS;}
194 if(fGraphM){fGraphM->Delete(); delete fGraphM;}
196 delete fReconstructor;
197 if(gGeoManager) delete gGeoManager;
198 if(fCl){fCl->Delete(); delete fCl;}
199 if(fTrklt){fTrklt->Delete(); delete fTrklt;}
200 if(fMCcl){fMCcl->Delete(); delete fMCcl;}
201 if(fMCtrklt){fMCtrklt->Delete(); delete fMCtrklt;}
205 //________________________________________________________
206 void AliTRDresolution::CreateOutputObjects()
208 // spatial resolution
209 OpenFile(0, "RECREATE");
211 fContainer = Histos();
213 fCl = new TObjArray();
214 fCl->SetOwner(kTRUE);
215 fTrklt = new TObjArray();
216 fTrklt->SetOwner(kTRUE);
217 fMCcl = new TObjArray();
218 fMCcl->SetOwner(kTRUE);
219 fMCtrklt = new TObjArray();
220 fMCtrklt->SetOwner(kTRUE);
223 //________________________________________________________
224 void AliTRDresolution::Exec(Option_t *opt)
235 AliTRDrecoTask::Exec(opt);
240 PostData(4, fMCtrklt);
243 //________________________________________________________
244 TH1* AliTRDresolution::PlotCharge(const AliTRDtrackV1 *track)
247 // Plots the charge distribution
250 if(track) fkTrack = track;
252 AliWarning("No Track defined.");
255 TObjArray *arr = NULL;
256 if(!(arr = ((TObjArray*)fContainer->At(kCharge)))){
257 AliWarning("No output container defined.");
262 AliTRDseedV1 *fTracklet = NULL;
263 AliTRDcluster *c = NULL;
264 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
265 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
266 if(!fTracklet->IsOK()) continue;
267 Float_t x0 = fTracklet->GetX0();
269 for(Int_t itb=AliTRDseedV1::kNtb; itb--;){
270 if(!(c = fTracklet->GetClusters(itb))){
271 if(!(c = fTracklet->GetClusters(AliTRDseedV1::kNtb+itb))) continue;
273 dq = fTracklet->GetdQdl(itb, &dl);
274 dl /= 0.15; // dl/dl0, dl0 = 1.5 mm for nominal vd
275 (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dq);
278 // if(!HasMCdata()) continue;
280 // Float_t pt0, y0, z0, dydx0, dzdx0;
281 // if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
288 //________________________________________________________
289 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
292 // Plot the cluster distributions
295 if(track) fkTrack = track;
297 AliWarning("No Track defined.");
300 TObjArray *arr = NULL;
301 if(!(arr = ((TObjArray*)fContainer->At(kCluster)))){
302 AliWarning("No output container defined.");
305 ULong_t status = fkESD ? fkESD->GetStatus():0;
308 Float_t x0, y0, z0, dy, dydx, dzdx;
309 AliTRDseedV1 *fTracklet = NULL;
310 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
311 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
312 if(!fTracklet->IsOK()) continue;
313 x0 = fTracklet->GetX0();
315 // retrive the track angle with the chamber
316 y0 = fTracklet->GetYref(0);
317 z0 = fTracklet->GetZref(0);
318 dydx = fTracklet->GetYref(1);
319 dzdx = fTracklet->GetZref(1);
320 fTracklet->GetCovRef(cov);
321 Float_t tilt = fTracklet->GetTilt();
322 AliTRDcluster *c = NULL;
323 fTracklet->ResetClusterIter(kFALSE);
324 while((c = fTracklet->PrevCluster())){
325 Float_t xc = c->GetX();
326 Float_t yc = c->GetY();
327 Float_t zc = c->GetZ();
328 Float_t dx = x0 - xc;
329 Float_t yt = y0 - dx*dydx;
330 Float_t zt = z0 - dx*dzdx;
331 yc -= tilt*(zc-zt); // tilt correction
334 //Float_t sx2 = dydx*c->GetSX(c->GetLocalTimeBin()); sx2*=sx2;
335 Float_t sy2 = c->GetSigmaY2();
336 if(sy2<=0.) continue;
337 ((TH2I*)arr->At(0))->Fill(dydx, dy);
338 ((TH2I*)arr->At(1))->Fill(dydx, dy/TMath::Sqrt(cov[0] /*+ sx2*/ + sy2));
341 // Get z-position with respect to anode wire
342 Int_t istk = fGeo->GetStack(c->GetDetector());
343 AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
344 Float_t row0 = pp->GetRow0();
345 Float_t d = row0 - zt + pp->GetAnodeWireOffset();
346 d -= ((Int_t)(2 * d)) / 2.0;
347 if (d > 0.25) d = 0.5 - d;
349 AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
351 clInfo->SetCluster(c);
352 Float_t covR[] = {cov[0], cov[1], cov[2]};
353 clInfo->SetGlobalPosition(yt, zt, dydx, dzdx, covR);
354 clInfo->SetResolution(dy);
355 clInfo->SetAnisochronity(d);
356 clInfo->SetDriftLength(dx);
357 (*DebugStream()) << "ClusterREC"
358 <<"status=" << status
359 <<"clInfo.=" << clInfo
364 return (TH2I*)arr->At(0);
368 //________________________________________________________
369 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
371 // Plot normalized residuals for tracklets to track.
373 // We start from the result that if X=N(|m|, |Cov|)
375 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
377 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
378 // reference position.
379 if(track) fkTrack = track;
381 AliWarning("No Track defined.");
384 TObjArray *arr = NULL;
385 if(!(arr = (TObjArray*)fContainer->At(kTrackTRD ))){
386 AliWarning("No output container defined.");
390 Double_t cov[3], covR[7]/*, sqr[3], inv[3]*/;
391 Float_t x, dx, dy, dz;
392 AliTRDseedV1 *fTracklet = NULL;
393 for(Int_t il=AliTRDgeometry::kNlayer; il--;){
394 if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
395 if(!fTracklet->IsOK()) continue;
396 x = fTracklet->GetX();
397 dx = fTracklet->GetX0() - x;
398 // compute dy^2 and dz^2
399 dy = fTracklet->GetYref(0)-dx*fTracklet->GetYref(1) - fTracklet->GetY();
400 dz = fTracklet->GetZref(0)-dx*fTracklet->GetZref(1) - fTracklet->GetZ();
401 // compute covariance matrix
402 fTracklet->GetCovAt(x, cov);
403 fTracklet->GetCovRef(covR);
404 cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
405 /* // Correct PULL calculation by considering off
406 // diagonal elements in the covariance matrix
407 // compute square root matrix
408 if(AliTRDseedV1::GetCovInv(cov, inv)==0.) continue;
409 if(AliTRDseedV1::GetCovSqrt(inv, sqr)<0.) continue;
410 Double_t y = sqr[0]*dy+sqr[1]*dz;
411 Double_t z = sqr[1]*dy+sqr[2]*dz;
412 ((TH3*)h)->Fill(y, z, fTracklet->GetYref(1));*/
414 ((TH2I*)arr->At(0))->Fill(fTracklet->GetYref(1), dy);
415 ((TH2I*)arr->At(1))->Fill(fTracklet->GetYref(1), dy/TMath::Sqrt(cov[0]));
416 ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), TMath::ATan((fTracklet->GetYref(1)-fTracklet->GetYfit(1))/(1-fTracklet->GetYref(1)*fTracklet->GetYfit(1))));
417 if(!fTracklet->IsRowCross()) continue;
418 ((TH2I*)arr->At(2))->Fill(fTracklet->GetZref(1), dz);
419 ((TH2I*)arr->At(3))->Fill(fTracklet->GetZref(1), dz/TMath::Sqrt(cov[2]));
423 return (TH2I*)arr->At(0);
427 //________________________________________________________
428 TH1* AliTRDresolution::PlotTrackTPC(const AliTRDtrackV1 *track)
430 // Store resolution/pulls of Kalman before updating with the TRD information
431 // at the radial position of the first tracklet. The following points are used
433 // - the (y,z,snp) of the first TRD tracklet
434 // - the (y, z, snp, tgl, pt) of the MC track reference
436 // Additionally the momentum resolution/pulls are calculated for usage in the
439 if(track) fkTrack = track;
441 AliWarning("No Track defined.");
444 AliExternalTrackParam *tin = NULL;
445 if(!(tin = fkTrack->GetTrackLow())){
446 AliWarning("Track did not entered TRD fiducial volume.");
451 Double_t x = tin->GetX();
452 AliTRDseedV1 *tracklet = NULL;
453 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
454 if(!(tracklet = fkTrack->GetTracklet(ily))) continue;
457 if(!tracklet || TMath::Abs(x-tracklet->GetX())>1.e-3){
458 AliWarning("Tracklet did not match TRD entrance.");
461 const Int_t kNPAR(5);
462 Double_t parR[kNPAR]; memcpy(parR, tin->GetParameter(), kNPAR*sizeof(Double_t));
463 Double_t covR[3*kNPAR]; memcpy(covR, tin->GetCovariance(), 3*kNPAR*sizeof(Double_t));
464 Double_t cov[3]; tracklet->GetCovAt(x, cov);
466 // define sum covariances
467 TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
468 Double_t *pc = &covR[0], *pp = &parR[0];
469 for(Int_t ir=0; ir<kNPAR; ir++, pp++){
471 for(Int_t ic = 0; ic<=ir; ic++,pc++){
472 COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
475 PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
476 //COV.Print(); PAR.Print();
478 //TODO Double_t dydx = TMath::Sqrt(1.-parR[2]*parR[2])/parR[2];
479 Double_t dy = parR[0] - tracklet->GetY();
480 TObjArray *arr = (TObjArray*)fContainer->At(kTrackTPC);
481 ((TH2I*)arr->At(0))->Fill(tracklet->GetYref(1), dy);
482 ((TH2I*)arr->At(1))->Fill(tracklet->GetYref(1), dy/TMath::Sqrt(COV(0,0)+cov[0]));
483 if(tracklet->IsRowCross()){
484 Double_t dz = parR[1] - tracklet->GetZ();
485 ((TH2I*)arr->At(2))->Fill(tracklet->GetZref(1), dz);
486 ((TH2I*)arr->At(3))->Fill(tracklet->GetZref(1), dz/TMath::Sqrt(COV(1,1)+cov[2]));
488 Double_t dphi = TMath::ASin(PAR[2])-TMath::ATan(tracklet->GetYfit(1)); ((TH2I*)arr->At(4))->Fill(tracklet->GetYref(1), dphi);
491 // register reference histo for mini-task
492 h = (TH2I*)arr->At(0);
495 (*DebugStream()) << "trackIn"
501 Double_t y = tracklet->GetY();
502 Double_t z = tracklet->GetZ();
503 (*DebugStream()) << "trackletIn"
513 if(!HasMCdata()) return h;
515 Float_t dx, pt0, x0=tracklet->GetX0(), y0, z0, dydx0, dzdx0;
516 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
517 // translate to reference radial position
518 dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
519 Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi))
521 TVectorD PARMC(kNPAR);
522 PARMC[0]=y0; PARMC[1]=z0;
523 PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
526 // TMatrixDSymEigen eigen(COV);
527 // TVectorD evals = eigen.GetEigenValues();
528 // TMatrixDSym evalsm(kNPAR);
529 // for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
530 // TMatrixD evecs = eigen.GetEigenVectors();
531 // TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
534 arr = (TObjArray*)fContainer->At(kMCtrackTPC);
535 // y resolution/pulls
536 ((TH2I*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0]);
537 ((TH2I*)arr->At(1))->Fill(dydx0, (PARMC[0]-PAR[0])/TMath::Sqrt(COV(0,0)));
538 // z resolution/pulls
539 ((TH2I*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1]);
540 ((TH2I*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)));
541 // phi resolution/snp pulls
542 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
543 ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
544 // theta resolution/tgl pulls
545 ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
546 ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
547 // pt resolution\\1/pt pulls\\p resolution/pull
548 for(Int_t is=AliPID::kSPECIES; is--;){
549 if(TMath::Abs(fkMC->GetPDG())!=AliPID::ParticleCode(is)) continue;
550 ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., is);
551 ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), is);
553 Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
555 p = tracklet->GetMomentum(&sp);
556 ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., is);
557 ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/sp, is);
563 (*DebugStream()) << "trackInMC"
570 //________________________________________________________
571 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
574 // Plot MC distributions
578 AliWarning("No MC defined. Results will not be available.");
581 if(track) fkTrack = track;
583 AliWarning("No Track defined.");
586 TObjArray *arr = NULL;
589 Int_t pdg = fkMC->GetPDG(), det=-1;
590 Int_t label = fkMC->GetLabel();
591 Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl;
592 Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
593 Double_t covR[7]/*, cov[3]*/;
596 Double_t dX[12], dY[12], dZ[12], dPt[12], cCOV[12][15];
597 fkMC->PropagateKalman(dX, dY, dZ, dPt, cCOV);
598 (*DebugStream()) << "MCkalman"
615 AliTRDReconstructor rec;
616 AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
617 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
618 if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
619 !fTracklet->IsOK())*/ continue;
621 det = fTracklet->GetDetector();
622 x0 = fTracklet->GetX0();
623 //radial shift with respect to the MC reference (radial position of the pad plane)
624 x= fTracklet->GetX();
625 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
626 xAnode = fTracklet->GetX0();
628 // MC track position at reference radial position
631 (*DebugStream()) << "MC"
642 Float_t ymc = y0 - dx*dydx0;
643 Float_t zmc = z0 - dx*dzdx0;
644 //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
646 // Kalman position at reference radial position
648 dydx = fTracklet->GetYref(1);
649 dzdx = fTracklet->GetZref(1);
650 dzdl = fTracklet->GetTgl();
651 y = fTracklet->GetYref(0) - dx*dydx;
653 z = fTracklet->GetZref(0) - dx*dzdx;
655 pt = TMath::Abs(fTracklet->GetPt());
656 fTracklet->GetCovRef(covR);
658 arr = (TObjArray*)fContainer->At(kMCtrackTRD);
659 // y resolution/pulls
660 ((TH2I*)arr->At(0))->Fill(dydx0, dy);
661 ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(covR[0]));
662 // z resolution/pulls
663 ((TH2I*)arr->At(2))->Fill(dzdx0, dz);
664 ((TH2I*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]));
665 // phi resolution/ snp pulls
666 Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
667 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
668 Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0);
669 ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3]));
670 // theta resolution/ tgl pulls
671 Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
672 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
673 ((TH2I*)arr->At(6))->Fill(dzdl0,
675 ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
676 // pt resolution \\ 1/pt pulls \\ p resolution for PID
677 for(Int_t is=AliPID::kSPECIES; is--;){
678 if(TMath::Abs(pdg)!=AliPID::ParticleCode(is)) continue;
679 ((TH3S*)((TObjArray*)arr->At(8))->At(ily))->Fill(pt0, pt/pt0-1., is);
680 ((TH3S*)((TObjArray*)arr->At(9))->At(ily))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), is);
681 Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
682 p = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
683 ((TH3S*)((TObjArray*)arr->At(10))->At(ily))->Fill(p0, p/p0-1., is);
687 // Fill Debug stream for Kalman track
689 (*DebugStream()) << "MCtrack"
701 // recalculate tracklet based on the MC info
702 AliTRDseedV1 tt(*fTracklet);
703 tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
704 tt.SetZref(1, dzdx0);
705 tt.SetReconstructor(&rec);
706 tt.Fit(kTRUE, kTRUE);
707 x= tt.GetX();y= tt.GetY();z= tt.GetZ();
708 dydx = tt.GetYfit(1);
712 Bool_t rc = tt.IsRowCross();
714 // add tracklet residuals for y and dydx
715 arr = (TObjArray*)fContainer->At(kMCtracklet);
719 Float_t dphi = (dydx - dydx0);
720 dphi /= (1.- dydx*dydx0);
722 ((TH2I*)arr->At(0))->Fill(dydx0, dy);
723 if(tt.GetS2Y()>0.) ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(tt.GetS2Y()));
724 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
726 // add tracklet residuals for z
728 ((TH2I*)arr->At(2))->Fill(dzdl0, dz);
729 if(tt.GetS2Z()>0.) ((TH2I*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z()));
732 // Fill Debug stream for tracklet
734 Float_t s2y = tt.GetS2Y();
735 Float_t s2z = tt.GetS2Z();
736 (*DebugStream()) << "MCtracklet"
747 Int_t istk = AliTRDgeometry::GetStack(det);
748 AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
749 Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
750 Float_t tilt = fTracklet->GetTilt();
751 //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
753 arr = (TObjArray*)fContainer->At(kMCcluster);
754 AliTRDcluster *c = NULL;
755 fTracklet->ResetClusterIter(kFALSE);
756 while((c = fTracklet->PrevCluster())){
757 Float_t q = TMath::Abs(c->GetQ());
758 x = c->GetX(); y = c->GetY();z = c->GetZ();
762 dy = (y - tilt*(z-zmc)) - ymc;
765 ((TH2I*)arr->At(0))->Fill(dydx0, dy);
766 ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(c->GetSigmaY2()));
769 // Fill calibration container
770 Float_t d = zr0 - zmc;
771 d -= ((Int_t)(2 * d)) / 2.0;
772 if (d > 0.25) d = 0.5 - d;
773 AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
774 clInfo->SetCluster(c);
775 clInfo->SetMC(pdg, label);
776 clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
777 clInfo->SetResolution(dy);
778 clInfo->SetAnisochronity(d);
779 clInfo->SetDriftLength(dx-.5*AliTRDgeometry::CamHght());
780 clInfo->SetTilt(tilt);
783 if(!clInfoArr) clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
784 clInfoArr->Add(clInfo);
788 if(DebugLevel()>=2 && clInfoArr){
789 (*DebugStream()) << "MCcluster"
790 <<"clInfo.=" << clInfoArr
792 delete clInfoArr; clInfoArr=NULL;
799 //________________________________________________________
800 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
803 // Get the reference figures
806 Float_t xy[4] = {0., 0., 0., 0.};
808 AliWarning("Please provide a canvas to draw results.");
811 TList *l = NULL; TVirtualPad *pad=NULL;
814 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
815 ((TVirtualPad*)l->At(0))->cd();
816 ((TGraphAsymmErrors*)((TObjArray*)fGraphM->At(kCharge))->At(0))->Draw("apl");
817 ((TVirtualPad*)l->At(1))->cd();
818 ((TGraphErrors*)((TObjArray*)fGraphS->At(kCharge))->At(0))->Draw("apl");
821 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
822 xy[0] = -.3; xy[1] = -200.; xy[2] = .3; xy[3] = 1000.;
823 ((TVirtualPad*)l->At(0))->cd();
824 if(!GetGraphPlot(&xy[0], kCluster, 0)) break;
825 xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
826 ((TVirtualPad*)l->At(1))->cd();
827 if(!GetGraphPlot(&xy[0], kCluster, 1)) break;
830 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
831 xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
832 ((TVirtualPad*)l->At(0))->cd();
833 if(!GetGraphPlot(&xy[0], kTrackTRD , 0)) break;
834 xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
835 ((TVirtualPad*)l->At(1))->cd();
836 if(!GetGraphPlot(&xy[0], kTrackTRD , 1)) break;
838 case 3: // kTrackTRD z
839 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
840 xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
841 ((TVirtualPad*)l->At(0))->cd();
842 if(!GetGraphPlot(&xy[0], kTrackTRD , 2)) break;
843 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
844 ((TVirtualPad*)l->At(1))->cd();
845 if(!GetGraphPlot(&xy[0], kTrackTRD , 3)) break;
847 case 4: // kTrackTRD phi
848 xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
849 if(GetGraphPlot(&xy[0], kTrackTRD , 4)) return kTRUE;
851 case 5: // kTrackTPC y
852 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
853 xy[0] = -.3; xy[1] = -500.; xy[2] = .3; xy[3] = 1500.;
854 pad = ((TVirtualPad*)l->At(0)); pad->cd();
855 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
856 if(!GetGraphPlot(&xy[0], kTrackTPC, 0)) break;
857 xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
858 pad=((TVirtualPad*)l->At(1)); pad->cd();
859 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
860 if(!GetGraphPlot(&xy[0], kTrackTPC, 1)) break;
862 case 6: // kTrackTPC z
863 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
864 xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
865 pad = ((TVirtualPad*)l->At(0)); pad->cd();
866 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
867 if(!GetGraphPlot(&xy[0], kTrackTPC, 2)) break;
868 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
869 pad = ((TVirtualPad*)l->At(1)); pad->cd();
870 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
871 if(!GetGraphPlot(&xy[0], kTrackTPC, 3)) break;
873 case 7: // kTrackTPC phi
874 xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
875 if(GetGraphPlot(&xy[0], kTrackTPC, 4)) return kTRUE;
877 case 8: // kMCcluster
878 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
879 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
880 ((TVirtualPad*)l->At(0))->cd();
881 if(!GetGraphPlot(&xy[0], kMCcluster, 0)) break;
882 xy[0] = -.3; xy[1] = -0.5; xy[2] = .3; xy[3] = 2.5;
883 ((TVirtualPad*)l->At(1))->cd();
884 if(!GetGraphPlot(&xy[0], kMCcluster, 1)) break;
886 case 9: //kMCtracklet [y]
887 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
888 xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =250.;
889 ((TVirtualPad*)l->At(0))->cd();
890 if(!GetGraphPlot(&xy[0], kMCtracklet, 0)) break;
891 xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
892 ((TVirtualPad*)l->At(1))->cd();
893 if(!GetGraphPlot(&xy[0], kMCtracklet, 1)) break;
895 case 10: //kMCtracklet [z]
896 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
897 xy[0]=-1.; xy[1]=-100.; xy[2]=1.; xy[3] =2500.;
898 ((TVirtualPad*)l->At(0))->cd();
899 if(!GetGraphPlot(&xy[0], kMCtracklet, 2)) break;
900 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
901 ((TVirtualPad*)l->At(1))->cd();
902 if(!GetGraphPlot(&xy[0], kMCtracklet, 3)) break;
904 case 11: //kMCtracklet [phi]
905 xy[0]=-.3; xy[1]=-3.; xy[2]=.3; xy[3] =25.;
906 if(!GetGraphPlot(&xy[0], kMCtracklet, 4)) break;
908 case 12: //kMCtrackTRD [y]
909 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
910 xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =200.;
911 ((TVirtualPad*)l->At(0))->cd();
912 if(!GetGraphPlot(&xy[0], kMCtrackTRD, 0)) break;
913 xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
914 ((TVirtualPad*)l->At(1))->cd();
915 if(!GetGraphPlot(&xy[0], kMCtrackTRD, 1)) break;
917 case 13: //kMCtrackTRD [z]
918 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
919 xy[0]=-1.; xy[1]=-700.; xy[2]=1.; xy[3] =1500.;
920 ((TVirtualPad*)l->At(0))->cd();
921 if(!GetGraphPlot(&xy[0], kMCtrackTRD, 2)) break;
922 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
923 ((TVirtualPad*)l->At(1))->cd();
924 if(!GetGraphPlot(&xy[0], kMCtrackTRD, 3)) break;
926 case 14: //kMCtrackTRD [phi/snp]
927 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
928 xy[0]=-.2; xy[1]=-0.2; xy[2]=.2; xy[3] =2.;
929 ((TVirtualPad*)l->At(0))->cd();
930 if(!GetGraphPlot(&xy[0], kMCtrackTRD, 4)) break;
931 xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
932 ((TVirtualPad*)l->At(1))->cd();
933 if(!GetGraphPlot(&xy[0], kMCtrackTRD, 5)) break;
935 case 15: //kMCtrackTRD [theta/tgl]
936 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
937 xy[0]=-1.; xy[1]=-0.5; xy[2]=1.; xy[3] =5.;
938 ((TVirtualPad*)l->At(0))->cd();
939 if(!GetGraphPlot(&xy[0], kMCtrackTRD, 6)) break;
940 xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
941 ((TVirtualPad*)l->At(1))->cd();
942 if(!GetGraphPlot(&xy[0], kMCtrackTRD, 7)) break;
944 case 16: //kMCtrackTRD [pt]
945 xy[0] = 0.; xy[1] = -5.; xy[2] = 12.; xy[3] = 7.;
946 gPad->Divide(2, 3, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
947 for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
948 pad = (TVirtualPad*)l->At(il); pad->cd();
949 pad->SetMargin(0.07, 0.07, 0.1, 0.);
950 if(!GetGraphTrack(&xy[0], 8, il)) break;
953 case 17: //kMCtrackTRD [1/pt] pulls
954 xy[0] = 0.; xy[1] = -1.5; xy[2] = 2.; xy[3] = 2.;
955 gPad->Divide(2, 3, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
956 for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
957 pad = (TVirtualPad*)l->At(il); pad->cd();
958 pad->SetMargin(0.07, 0.07, 0.1, 0.);
959 if(!GetGraphTrack(&xy[0], 9, il)) break;
962 case 18: //kMCtrackTRD [p]
963 xy[0] = 0.; xy[1] = -7.5; xy[2] = 12.; xy[3] = 10.5;
964 gPad->Divide(2, 3, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
965 for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
966 pad = (TVirtualPad*)l->At(il); pad->cd();
967 pad->SetMargin(0.07, 0.07, 0.1, 0.);
968 if(!GetGraphTrack(&xy[0], 10, il)) break;
971 case 19: // kMCtrackTPC [y]
972 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
973 xy[0]=-.25; xy[1]=-50.; xy[2]=.25; xy[3] =800.;
974 ((TVirtualPad*)l->At(0))->cd();
975 if(!GetGraphPlot(&xy[0], kMCtrackTPC, 0)) break;
976 xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 2.5;
977 ((TVirtualPad*)l->At(1))->cd();
978 if(!GetGraphPlot(&xy[0], kMCtrackTPC, 1)) break;
980 case 20: // kMCtrackTPC [z]
981 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
982 xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =800.;
983 ((TVirtualPad*)l->At(0))->cd();
984 if(!GetGraphPlot(&xy[0], kMCtrackTPC, 2)) break;
985 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
986 ((TVirtualPad*)l->At(1))->cd();
987 if(!GetGraphPlot(&xy[0], kMCtrackTPC, 3)) break;
989 case 21: // kMCtrackTPC [phi|snp]
990 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
991 xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
992 ((TVirtualPad*)l->At(0))->cd();
993 if(!GetGraphPlot(&xy[0], kMCtrackTPC, 4)) break;
994 xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
995 ((TVirtualPad*)l->At(1))->cd();
996 if(!GetGraphPlot(&xy[0], kMCtrackTPC, 5)) break;
998 case 22: // kMCtrackTPC [theta|tgl]
999 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1000 xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.;
1001 ((TVirtualPad*)l->At(0))->cd();
1002 if(!GetGraphPlot(&xy[0], kMCtrackTPC, 6)) break;
1003 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 1.5;
1004 ((TVirtualPad*)l->At(1))->cd();
1005 if(!GetGraphPlot(&xy[0], kMCtrackTPC, 7)) break;
1007 case 23: // kMCtrackTPC [pt]
1008 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1009 xy[0] = 0.; xy[1] = -.8; xy[2] = 12.; xy[3] = 2.3;
1010 ((TVirtualPad*)l->At(0))->cd();
1011 if(!GetGraphTrackTPC(xy, 8)) break;
1012 xy[0]=0.; xy[1]=-0.5; xy[2]=2.; xy[3] =2.5;
1013 ((TVirtualPad*)l->At(1))->cd();
1014 if(!GetGraphTrackTPC(xy, 9)) break;
1016 case 24: // kMCtrackTPC [p]
1017 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1018 xy[0] = 0.; xy[1] = -.8; xy[2] = 12.; xy[3] = 2.3;
1019 pad = ((TVirtualPad*)l->At(0));pad->cd();
1020 pad->SetMargin(0.12, 0.12, 0.1, 0.04);
1021 if(!GetGraphTrackTPC(xy, 10)) break;
1022 xy[0]=0.; xy[1]=-1.5; xy[2]=12.; xy[3] =2.5;
1023 pad = ((TVirtualPad*)l->At(1)); pad->cd();
1024 pad->SetMargin(0.12, 0.12, 0.1, 0.04);
1025 if(!GetGraphTrackTPC(xy, 11)) break;
1027 case 25: // kMCtrackTOF [z]
1030 AliWarning(Form("Reference plot [%d] missing result", ifig));
1035 //________________________________________________________
1036 Bool_t AliTRDresolution::PostProcess()
1038 //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
1040 AliError("ERROR: list not available");
1043 TGraph *gm= NULL, *gs= NULL;
1044 if(!fGraphS && !fGraphM){
1045 TObjArray *aM(NULL), *aS(NULL), *a(NULL);
1046 Int_t n = fContainer->GetEntriesFast();
1047 fGraphS = new TObjArray(n); fGraphS->SetOwner();
1048 fGraphM = new TObjArray(n); fGraphM->SetOwner();
1049 for(Int_t ig=0; ig<n; ig++){
1050 fGraphM->AddAt(aM = new TObjArray(fgNElements[ig]), ig);
1051 fGraphS->AddAt(aS = new TObjArray(fgNElements[ig]), ig);
1053 for(Int_t ic=0; ic<fgNElements[ig]; ic++){
1054 if(ig==kMCtrackTPC&&(ic>=8&&ic<=12)){ // TPC momentum plot
1055 aS->AddAt(a = new TObjArray(AliPID::kSPECIES), ic);
1056 for(Int_t is=AliPID::kSPECIES; is--;){
1057 a->AddAt(gs = new TGraphErrors(), is);
1058 gs->SetMarkerStyle(23);
1059 gs->SetMarkerColor(is ? kRed : kMagenta);
1060 gs->SetLineStyle(is);
1061 gs->SetLineColor(is ? kRed : kMagenta);
1062 gs->SetLineWidth(is ? 1 : 3);
1063 gs->SetNameTitle(Form("s_%d%02d%d", ig, ic, is), "");
1065 aM->AddAt(a = new TObjArray(AliPID::kSPECIES), ic);
1066 for(Int_t is=AliPID::kSPECIES; is--;){
1067 a->AddAt(gm = new TGraphErrors(), is);
1068 gm->SetLineColor(is ? kBlack : kBlue);
1069 gm->SetLineStyle(is);
1070 gm->SetMarkerStyle(7);
1071 gm->SetMarkerColor(is ? kBlack : kBlue);
1072 gm->SetLineWidth(is ? 1 : 3);
1073 gm->SetNameTitle(Form("m_%d%02d%d", ig, ic, is), "");
1076 } else if(ig==kMCtrackTRD&&(ic==8||ic==9||ic==10)){ // TRD momentum plot
1077 TObjArray *aaS, *aaM;
1078 aS->AddAt(aaS = new TObjArray(AliTRDgeometry::kNlayer), ic);
1079 aM->AddAt(aaM = new TObjArray(AliTRDgeometry::kNlayer), ic);
1080 for(Int_t il=AliTRDgeometry::kNlayer; il--;){
1081 aaS->AddAt(a = new TObjArray(AliPID::kSPECIES), il);
1082 for(Int_t is=AliPID::kSPECIES; is--;){
1083 a->AddAt(gs = new TGraphErrors(), is);
1084 gs->SetMarkerStyle(23);
1085 gs->SetMarkerColor(is ? kRed : kMagenta);
1086 gs->SetLineStyle(is);
1087 gs->SetLineColor(is ? kRed : kMagenta);
1088 gs->SetLineWidth(is ? 1 : 3);
1089 gs->SetNameTitle(Form("s_%d%02d%d%d", ig, ic, is, il), "");
1091 aaM->AddAt(a = new TObjArray(AliPID::kSPECIES), il);
1092 for(Int_t is=AliPID::kSPECIES; is--;){
1093 a->AddAt(gm = new TGraphErrors(), is);
1094 gm->SetMarkerStyle(7);
1095 gm->SetMarkerColor(is ? kBlack : kBlue);
1096 gm->SetLineStyle(is);
1097 gm->SetLineColor(is ? kBlack : kBlue);
1098 gm->SetLineWidth(is ? 1 : 3);
1099 gm->SetNameTitle(Form("m_%d%02d%d%d", ig, ic, is, il), "");
1105 aS->AddAt(gs = new TGraphErrors(), ic);
1106 gs->SetMarkerStyle(23);
1107 gs->SetMarkerColor(kRed);
1108 gs->SetLineColor(kRed);
1109 gs->SetNameTitle(Form("s_%d%02d", ig, ic), "");
1111 aM->AddAt(gm = ig ? (TGraph*)new TGraphErrors() : (TGraph*)new TGraphAsymmErrors(), ic);
1112 gm->SetLineColor(kBlack);
1113 gm->SetMarkerStyle(7);
1114 gm->SetMarkerColor(kBlack);
1115 gm->SetNameTitle(Form("m_%d%02d", ig, ic), "");
1120 /* printf("\n\n\n"); fGraphS->ls();
1121 printf("\n\n\n"); fGraphM->ls();*/
1126 TF1 fg("fGauss", "gaus", -.5, .5);
1127 // Landau for charge resolution
1128 TF1 fl("fLandau", "landau", 0., 1000.);
1130 //PROCESS EXPERIMENTAL DISTRIBUTIONS
1131 // Charge resolution
1132 //Process3DL(kCharge, 0, &fl);
1133 // Clusters residuals
1134 Process2D(kCluster, 0, &fg, 1.e4);
1135 Process2D(kCluster, 1, &fg);
1137 // Tracklet residual/pulls
1138 Process2D(kTrackTRD , 0, &fg, 1.e4);
1139 Process2D(kTrackTRD , 1, &fg);
1140 Process2D(kTrackTRD , 2, &fg, 1.e4);
1141 Process2D(kTrackTRD , 3, &fg);
1142 Process2D(kTrackTRD , 4, &fg, 1.e3);
1144 // TPC track residual/pulls
1145 Process2D(kTrackTPC, 0, &fg, 1.e4);
1146 Process2D(kTrackTPC, 1, &fg);
1147 Process2D(kTrackTPC, 2, &fg, 1.e4);
1148 Process2D(kTrackTPC, 3, &fg);
1149 Process2D(kTrackTPC, 4, &fg, 1.e3);
1152 if(!HasMCdata()) return kTRUE;
1155 //PROCESS MC RESIDUAL DISTRIBUTIONS
1157 // CLUSTER Y RESOLUTION/PULLS
1158 Process2D(kMCcluster, 0, &fg, 1.e4);
1159 Process2D(kMCcluster, 1, &fg);
1162 // TRACKLET RESOLUTION/PULLS
1163 Process2D(kMCtracklet, 0, &fg, 1.e4); // y
1164 Process2D(kMCtracklet, 1, &fg); // y pulls
1165 Process2D(kMCtracklet, 2, &fg, 1.e4); // z
1166 Process2D(kMCtracklet, 3, &fg); // z pulls
1167 Process2D(kMCtracklet, 4, &fg, 1.e3); // phi
1170 // TRACK RESOLUTION/PULLS
1171 Process2D(kMCtrackTRD, 0, &fg, 1.e4); // y
1172 Process2D(kMCtrackTRD, 1, &fg); // y PULL
1173 Process2D(kMCtrackTRD, 2, &fg, 1.e4); // z
1174 Process2D(kMCtrackTRD, 3, &fg); // z PULL
1175 Process2D(kMCtrackTRD, 4, &fg, 1.e3); // phi
1176 Process2D(kMCtrackTRD, 5, &fg); // snp PULL
1177 Process2D(kMCtrackTRD, 6, &fg, 1.e3); // theta
1178 Process2D(kMCtrackTRD, 7, &fg); // tgl PULL
1179 Process4D(kMCtrackTRD, 8, &fg, 1.e2); // pt resolution
1180 Process4D(kMCtrackTRD, 9, &fg); // 1/pt pulls
1181 Process4D(kMCtrackTRD, 10, &fg, 1.e2); // p resolution
1184 // TRACK TPC RESOLUTION/PULLS
1185 Process2D(kMCtrackTPC, 0, &fg, 1.e4);// y resolution
1186 Process2D(kMCtrackTPC, 1, &fg); // y pulls
1187 Process2D(kMCtrackTPC, 2, &fg, 1.e4);// z resolution
1188 Process2D(kMCtrackTPC, 3, &fg); // z pulls
1189 Process2D(kMCtrackTPC, 4, &fg, 1.e3);// phi resolution
1190 Process2D(kMCtrackTPC, 5, &fg); // snp pulls
1191 Process2D(kMCtrackTPC, 6, &fg, 1.e3);// theta resolution
1192 Process2D(kMCtrackTPC, 7, &fg); // tgl pulls
1193 Process3D(kMCtrackTPC, 8, &fg, 1.e2);// pt resolution
1194 Process3D(kMCtrackTPC, 9, &fg); // 1/pt pulls
1195 Process3D(kMCtrackTPC, 10, &fg, 1.e2);// p resolution
1196 Process3D(kMCtrackTPC, 11, &fg); // p pulls
1199 // TRACK HMPID RESOLUTION/PULLS
1200 Process2D(kMCtrackTOF, 0, &fg, 1.e4); // z towards TOF
1201 Process2D(kMCtrackTOF, 1, &fg); // z towards TOF
1208 //________________________________________________________
1209 void AliTRDresolution::Terminate(Option_t *opt)
1211 AliTRDrecoTask::Terminate(opt);
1212 if(HasPostProcess()) PostProcess();
1215 //________________________________________________________
1216 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1218 // Helper function to avoid duplication of code
1219 // Make first guesses on the fit parameters
1221 // find the intial parameters of the fit !! (thanks George)
1222 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1224 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1225 f->SetParLimits(0, 0., 3.*sum);
1226 f->SetParameter(0, .9*sum);
1227 Double_t rms = h->GetRMS();
1228 f->SetParLimits(1, -rms, rms);
1229 f->SetParameter(1, h->GetMean());
1231 f->SetParLimits(2, 0., 2.*rms);
1232 f->SetParameter(2, rms);
1233 if(f->GetNpar() <= 4) return;
1235 f->SetParLimits(3, 0., sum);
1236 f->SetParameter(3, .1*sum);
1238 f->SetParLimits(4, -.3, .3);
1239 f->SetParameter(4, 0.);
1241 f->SetParLimits(5, 0., 1.e2);
1242 f->SetParameter(5, 2.e-1);
1245 //________________________________________________________
1246 TObjArray* AliTRDresolution::Histos()
1249 // Define histograms
1252 if(fContainer) return fContainer;
1254 fContainer = new TObjArray(kNhistos);
1255 //fContainer->SetOwner(kTRUE);
1257 TObjArray *arr = NULL;
1259 const Int_t kPbins(12); // binning in momentum range should depend on the statistics analyzed
1261 // cluster to track residuals/pulls
1262 fContainer->AddAt(arr = new TObjArray(fgNElements[kCharge]), kCharge);
1263 arr->SetName("Charge");
1264 if(!(h = (TH3S*)gROOT->FindObject("hCharge"))){
1265 h = new TH3S("hCharge", "Charge Resolution", 20, 1., 2., 24, 0., 3.6, 100, 0., 500.);
1266 h->GetXaxis()->SetTitle("dx/dx_{0}");
1267 h->GetYaxis()->SetTitle("x_{d} [cm]");
1268 h->GetZaxis()->SetTitle("dq/dx [ADC/cm]");
1272 // cluster to track residuals/pulls
1273 fContainer->AddAt(arr = new TObjArray(fgNElements[kCluster]), kCluster);
1275 if(!(h = (TH2I*)gROOT->FindObject("hCl"))){
1276 h = new TH2I("hCl", "Cluster Residuals", 21, -.33, .33, 100, -.5, .5);
1277 h->GetXaxis()->SetTitle("tg(#phi)");
1278 h->GetYaxis()->SetTitle("#Delta y [cm]");
1279 h->GetZaxis()->SetTitle("entries");
1282 if(!(h = (TH2I*)gROOT->FindObject("hClpull"))){
1283 h = new TH2I("hClpull", "Cluster Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1284 h->GetXaxis()->SetTitle("tg(#phi)");
1285 h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1286 h->GetZaxis()->SetTitle("entries");
1290 // tracklet to track residuals/pulls in y direction
1291 fContainer->AddAt(arr = new TObjArray(fgNElements[kTrackTRD ]), kTrackTRD );
1292 arr->SetName("Trklt");
1293 if(!(h = (TH2I*)gROOT->FindObject("hTrkltY"))){
1294 h = new TH2I("hTrkltY", "Tracklet Y Residuals", 21, -.33, .33, 100, -.5, .5);
1295 h->GetXaxis()->SetTitle("#tg(#phi)");
1296 h->GetYaxis()->SetTitle("#Delta y [cm]");
1297 h->GetZaxis()->SetTitle("entries");
1300 if(!(h = (TH2I*)gROOT->FindObject("hTrkltYpull"))){
1301 h = new TH2I("hTrkltYpull", "Tracklet Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1302 h->GetXaxis()->SetTitle("#tg(#phi)");
1303 h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1304 h->GetZaxis()->SetTitle("entries");
1307 // tracklet to track residuals/pulls in z direction
1308 if(!(h = (TH2I*)gROOT->FindObject("hTrkltZ"))){
1309 h = new TH2I("hTrkltZ", "Tracklet Z Residuals", 50, -1., 1., 100, -1.5, 1.5);
1310 h->GetXaxis()->SetTitle("#tg(#theta)");
1311 h->GetYaxis()->SetTitle("#Delta z [cm]");
1312 h->GetZaxis()->SetTitle("entries");
1315 if(!(h = (TH2I*)gROOT->FindObject("hTrkltZpull"))){
1316 h = new TH2I("hTrkltZpull", "Tracklet Z Pulls", 50, -1., 1., 100, -5.5, 5.5);
1317 h->GetXaxis()->SetTitle("#tg(#theta)");
1318 h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
1319 h->GetZaxis()->SetTitle("entries");
1322 // tracklet to track phi residuals
1323 if(!(h = (TH2I*)gROOT->FindObject("hTrkltPhi"))){
1324 h = new TH2I("hTrkltPhi", "Tracklet #phi Residuals", 21, -.33, .33, 100, -.5, .5);
1325 h->GetXaxis()->SetTitle("tg(#phi)");
1326 h->GetYaxis()->SetTitle("#Delta phi [rad]");
1327 h->GetZaxis()->SetTitle("entries");
1332 // tracklet to TPC track residuals/pulls in y direction
1333 fContainer->AddAt(arr = new TObjArray(fgNElements[kTrackTPC]), kTrackTPC);
1334 arr->SetName("TrkTPC");
1335 if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCY"))){
1336 h = new TH2I("hTrkTPCY", "Track[TPC] Y Residuals", 21, -.33, .33, 100, -.5, .5);
1337 h->GetXaxis()->SetTitle("#tg(#phi)");
1338 h->GetYaxis()->SetTitle("#Delta y [cm]");
1339 h->GetZaxis()->SetTitle("entries");
1342 if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCYpull"))){
1343 h = new TH2I("hTrkTPCYpull", "Track[TPC] Y Pulls", 21, -.33, .33, 100, -4.5, 4.5);
1344 h->GetXaxis()->SetTitle("#tg(#phi)");
1345 h->GetYaxis()->SetTitle("#Delta y/#sigma_{y}");
1346 h->GetZaxis()->SetTitle("entries");
1349 // tracklet to TPC track residuals/pulls in z direction
1350 if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZ"))){
1351 h = new TH2I("hTrkTPCZ", "Track[TPC] Z Residuals", 50, -1., 1., 100, -1.5, 1.5);
1352 h->GetXaxis()->SetTitle("#tg(#theta)");
1353 h->GetYaxis()->SetTitle("#Delta z [cm]");
1354 h->GetZaxis()->SetTitle("entries");
1357 if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCZpull"))){
1358 h = new TH2I("hTrkTPCZpull", "Track[TPC] Z Pulls", 50, -1., 1., 100, -5.5, 5.5);
1359 h->GetXaxis()->SetTitle("#tg(#theta)");
1360 h->GetYaxis()->SetTitle("#Delta z/#sigma_{z}");
1361 h->GetZaxis()->SetTitle("entries");
1364 // tracklet to TPC track phi residuals
1365 if(!(h = (TH2I*)gROOT->FindObject("hTrkTPCPhi"))){
1366 h = new TH2I("hTrkTPCPhi", "Track[TPC] #phi Residuals", 21, -.33, .33, 100, -.5, .5);
1367 h->GetXaxis()->SetTitle("tg(#phi)");
1368 h->GetYaxis()->SetTitle("#Delta phi [rad]");
1369 h->GetZaxis()->SetTitle("entries");
1374 // Resolution histos
1375 if(!HasMCdata()) return fContainer;
1377 // cluster y resolution [0]
1378 fContainer->AddAt(arr = new TObjArray(fgNElements[kMCcluster]), kMCcluster);
1379 arr->SetName("McCl");
1380 if(!(h = (TH2I*)gROOT->FindObject("hMcCl"))){
1381 h = new TH2I("hMcCl", "Cluster Resolution", 48, -.48, .48, 100, -.3, .3);
1382 h->GetXaxis()->SetTitle("tg(#phi)");
1383 h->GetYaxis()->SetTitle("#Delta y [cm]");
1384 h->GetZaxis()->SetTitle("entries");
1387 if(!(h = (TH2I*)gROOT->FindObject("hMcClPull"))){
1388 h = new TH2I("hMcClPull", "Cluster Pulls", 48, -.48, .48, 100, -4.5, 4.5);
1389 h->GetXaxis()->SetTitle("tg(#phi)");
1390 h->GetYaxis()->SetTitle("#Deltay/#sigma_{y}");
1391 h->GetZaxis()->SetTitle("entries");
1396 // TRACKLET RESOLUTION
1397 fContainer->AddAt(arr = new TObjArray(fgNElements[kMCtracklet]), kMCtracklet);
1398 arr->SetName("McTrklt");
1399 // tracklet y resolution
1400 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltY"))){
1401 h = new TH2I("hMcTrkltY", "Tracklet Resolution (Y)", 48, -.48, .48, 100, -.2, .2);
1402 h->GetXaxis()->SetTitle("tg(#phi)");
1403 h->GetYaxis()->SetTitle("#Delta y [cm]");
1404 h->GetZaxis()->SetTitle("entries");
1408 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltYPull"))){
1409 h = new TH2I("hMcTrkltYPull", "Tracklet Pulls (Y)", 48, -.48, .48, 100, -4.5, 4.5);
1410 h->GetXaxis()->SetTitle("tg(#phi)");
1411 h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1412 h->GetZaxis()->SetTitle("entries");
1415 // tracklet z resolution
1416 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltZ"))){
1417 h = new TH2I("hMcTrkltZ", "Tracklet Resolution (Z)", 100, -1., 1., 100, -1., 1.);
1418 h->GetXaxis()->SetTitle("tg(#theta)");
1419 h->GetYaxis()->SetTitle("#Delta z [cm]");
1420 h->GetZaxis()->SetTitle("entries");
1424 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltZPull"))){
1425 h = new TH2I("hMcTrkltZPull", "Tracklet Pulls (Z)", 100, -1., 1., 100, -3.5, 3.5);
1426 h->GetXaxis()->SetTitle("tg(#theta)");
1427 h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1428 h->GetZaxis()->SetTitle("entries");
1431 // tracklet phi resolution
1432 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkltPhi"))){
1433 h = new TH2I("hMcTrkltPhi", "Tracklet Resolution (#Phi)", 48, -.48, .48, 100, -.15, .15);
1434 h->GetXaxis()->SetTitle("tg(#phi)");
1435 h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1436 h->GetZaxis()->SetTitle("entries");
1441 // KALMAN TRACK RESOLUTION
1442 fContainer->AddAt(arr = new TObjArray(fgNElements[kMCtrackTRD]), kMCtrackTRD);
1443 arr->SetName("McTrkTRD");
1444 // Kalman track y resolution
1445 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkY"))){
1446 h = new TH2I("hMcTrkY", "Track Y Resolution", 48, -.48, .48, 100, -.2, .2);
1447 h->GetXaxis()->SetTitle("tg(#phi)");
1448 h->GetYaxis()->SetTitle("#Delta y [cm]");
1449 h->GetZaxis()->SetTitle("entries");
1452 // Kalman track y pulls
1453 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkYPull"))){
1454 h = new TH2I("hMcTrkYPull", "Track Y Pulls", 48, -.48, .48, 100, -4., 4.);
1455 h->GetXaxis()->SetTitle("tg(#phi)");
1456 h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1457 h->GetZaxis()->SetTitle("entries");
1461 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkZ"))){
1462 h = new TH2I("hMcTrkZ", "Track Z Resolution", 100, -1., 1., 100, -1., 1.);
1463 h->GetXaxis()->SetTitle("tg(#theta)");
1464 h->GetYaxis()->SetTitle("#Delta z [cm]");
1465 h->GetZaxis()->SetTitle("entries");
1468 // Kalman track Z pulls
1469 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkZPull"))){
1470 h = new TH2I("hMcTrkZPull", "Track Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1471 h->GetXaxis()->SetTitle("tg(#theta)");
1472 h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1473 h->GetZaxis()->SetTitle("entries");
1477 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNP"))){
1478 h = new TH2I("hMcTrkSNP", "Track Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3);
1479 h->GetXaxis()->SetTitle("tg(#phi)");
1480 h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1481 h->GetZaxis()->SetTitle("entries");
1484 // Kalman track SNP pulls
1485 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNPPull"))){
1486 h = new TH2I("hMcTrkSNPPull", "Track SNP Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1487 h->GetXaxis()->SetTitle("tg(#phi)");
1488 h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
1489 h->GetZaxis()->SetTitle("entries");
1493 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGL"))){
1494 h = new TH2I("hMcTrkTGL", "Track Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3);
1495 h->GetXaxis()->SetTitle("tg(#theta)");
1496 h->GetYaxis()->SetTitle("#Delta#theta [rad]");
1497 h->GetZaxis()->SetTitle("entries");
1500 // Kalman track TGL pulls
1501 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGLPull"))){
1502 h = new TH2I("hMcTrkTGLPull", "Track TGL Pulls", 100, -1., 1., 100, -4.5, 4.5);
1503 h->GetXaxis()->SetTitle("tg(#theta)");
1504 h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
1505 h->GetZaxis()->SetTitle("entries");
1508 // Kalman track Pt resolution
1509 const Int_t n = AliPID::kSPECIES;
1510 TObjArray *arr2 = NULL; TH3S* h3=NULL;
1511 arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 8);
1512 arr2->SetName("Track Pt Resolution");
1513 for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1514 if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkPt%d", il)))){
1515 h3 = new TH3S(Form("hMcTrkPt%d", il), "Track Pt Resolution", kPbins, 0., 12., 150, -.1, .2, n, -.5, n-.5);
1516 h3->GetXaxis()->SetTitle("p_{t} [GeV/c]");
1517 h3->GetYaxis()->SetTitle("#Delta p_{t}/p_{t}^{MC}");
1518 h3->GetZaxis()->SetTitle("SPECIES");
1520 arr2->AddAt(h3, il);
1522 // Kalman track Pt pulls
1523 arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 9);
1524 arr2->SetName("Track 1/Pt Pulls");
1525 for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1526 if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkPtPulls%d", il)))){
1527 h3 = new TH3S(Form("hMcTrkPtPulls%d", il), "Track 1/Pt Pulls", kPbins, 0., 2., 100, -4., 4., n, -.5, n-.5);
1528 h3->GetXaxis()->SetTitle("1/p_{t}^{MC} [c/GeV]");
1529 h3->GetYaxis()->SetTitle("#Delta(1/p_{t})/#sigma(1/p_{t}) ");
1530 h3->GetZaxis()->SetTitle("SPECIES");
1532 arr2->AddAt(h3, il);
1534 // Kalman track P resolution
1535 arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 10);
1536 arr2->SetName("Track P Resolution [PID]");
1537 for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++){
1538 if(!(h3 = (TH3S*)gROOT->FindObject(Form("hMcTrkP%d", il)))){
1539 h3 = new TH3S(Form("hMcTrkP%d", il), "Track P Resolution", kPbins, 0., 12., 150, -.15, .35, n, -.5, n-.5);
1540 h3->GetXaxis()->SetTitle("p [GeV/c]");
1541 h3->GetYaxis()->SetTitle("#Delta p/p^{MC}");
1542 h3->GetZaxis()->SetTitle("SPECIES");
1544 arr2->AddAt(h3, il);
1547 // TPC TRACK RESOLUTION
1548 fContainer->AddAt(arr = new TObjArray(fgNElements[kMCtrackTPC]), kMCtrackTPC);
1549 arr->SetName("McTrkTPC");
1551 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCY"))){
1552 h = new TH2I("hMcTrkTPCY", "Track[TPC] Y Resolution", 60, -.3, .3, 100, -.5, .5);
1553 h->GetXaxis()->SetTitle("tg(#phi)");
1554 h->GetYaxis()->SetTitle("#Delta y [cm]");
1555 h->GetZaxis()->SetTitle("entries");
1558 // Kalman track Y pulls
1559 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCYPull"))){
1560 h = new TH2I("hMcTrkTPCYPull", "Track[TPC] Y Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1561 h->GetXaxis()->SetTitle("tg(#phi)");
1562 h->GetYaxis()->SetTitle("#Delta y / #sigma_{y}");
1563 h->GetZaxis()->SetTitle("entries");
1567 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCZ"))){
1568 h = new TH2I("hMcTrkTPCZ", "Track[TPC] Z Resolution", 100, -1., 1., 100, -1., 1.);
1569 h->GetXaxis()->SetTitle("tg(#theta)");
1570 h->GetYaxis()->SetTitle("#Delta z [cm]");
1571 h->GetZaxis()->SetTitle("entries");
1574 // Kalman track Z pulls
1575 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCZPull"))){
1576 h = new TH2I("hMcTrkTPCZPull", "Track[TPC] Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1577 h->GetXaxis()->SetTitle("tg(#theta)");
1578 h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1579 h->GetZaxis()->SetTitle("entries");
1583 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCSNP"))){
1584 h = new TH2I("hMcTrkTPCSNP", "Track[TPC] Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3);
1585 h->GetXaxis()->SetTitle("tg(#phi)");
1586 h->GetYaxis()->SetTitle("#Delta #phi [rad]");
1587 h->GetZaxis()->SetTitle("entries");
1590 // Kalman track SNP pulls
1591 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCSNPPull"))){
1592 h = new TH2I("hMcTrkTPCSNPPull", "Track[TPC] SNP Pulls", 60, -.3, .3, 100, -4.5, 4.5);
1593 h->GetXaxis()->SetTitle("tg(#phi)");
1594 h->GetYaxis()->SetTitle("#Delta(sin(#phi)) / #sigma_{sin(#phi)}");
1595 h->GetZaxis()->SetTitle("entries");
1599 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCTGL"))){
1600 h = new TH2I("hMcTrkTPCTGL", "Track[TPC] Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3);
1601 h->GetXaxis()->SetTitle("tg(#theta)");
1602 h->GetYaxis()->SetTitle("#Delta#theta [rad]");
1603 h->GetZaxis()->SetTitle("entries");
1606 // Kalman track TGL pulls
1607 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCTGLPull"))){
1608 h = new TH2I("hMcTrkTPCTGLPull", "Track[TPC] TGL Pulls", 100, -1., 1., 100, -4.5, 4.5);
1609 h->GetXaxis()->SetTitle("tg(#theta)");
1610 h->GetYaxis()->SetTitle("#Delta(tg(#theta)) / #sigma_{tg(#theta)}");
1611 h->GetZaxis()->SetTitle("entries");
1614 // Kalman track Pt resolution
1615 if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPt"))){
1616 h3 = new TH3S("hMcTrkTPCPt", "Track[TPC] Pt Resolution", kPbins, 0., 12., 100, -.1, .2, n, -.5, n-.5);
1617 h3->GetXaxis()->SetTitle("p_{t} [GeV/c]");
1618 h3->GetYaxis()->SetTitle("#Delta p_{t}/p_{t}^{MC}");
1619 h3->GetZaxis()->SetTitle("SPECIES");
1622 // Kalman track Pt pulls
1623 if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPtPulls"))){
1624 h3 = new TH3S("hMcTrkTPCPtPulls", "Track[TPC] 1/Pt Pulls", kPbins, 0., 2., 100, -4., 4., n, -.5, n-.5);
1625 h3->GetXaxis()->SetTitle("1/p_{t}^{MC} [c/GeV]");
1626 h3->GetYaxis()->SetTitle("#Delta(1/p_{t})/#sigma(1/p_{t}) ");
1627 h3->GetZaxis()->SetTitle("SPECIES");
1630 // Kalman track P resolution
1631 if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCP"))){
1632 h3 = new TH3S("hMcTrkTPCP", "Track[TPC] P Resolution", kPbins, 0., 12., 100, -.15, .35, n, -.5, n-.5);
1633 h3->GetXaxis()->SetTitle("p [GeV/c]");
1634 h3->GetYaxis()->SetTitle("#Delta p/p^{MC}");
1635 h3->GetZaxis()->SetTitle("SPECIES");
1638 // Kalman track Pt pulls
1639 if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPPulls"))){
1640 h3 = new TH3S("hMcTrkTPCPPulls", "Track[TPC] P Pulls", kPbins, 0., 12., 100, -5., 5., n, -.5, n-.5);
1641 h3->GetXaxis()->SetTitle("p^{MC} [GeV/c]");
1642 h3->GetYaxis()->SetTitle("#Deltap/#sigma_{p}");
1643 h3->GetZaxis()->SetTitle("SPECIES");
1649 // Kalman track Z resolution [TOF]
1650 fContainer->AddAt(arr = new TObjArray(fgNElements[kMCtrackTOF]), kMCtrackTOF);
1651 arr->SetName("McTrkTOF");
1652 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTOFZ"))){
1653 h = new TH2I("hMcTrkTOFZ", "Track[TOF] Z Resolution", 100, -1., 1., 100, -1., 1.);
1654 h->GetXaxis()->SetTitle("tg(#theta)");
1655 h->GetYaxis()->SetTitle("#Delta z [cm]");
1656 h->GetZaxis()->SetTitle("entries");
1659 // Kalman track Z pulls
1660 if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTOFZPull"))){
1661 h = new TH2I("hMcTrkTOFZPull", "Track[TOF] Z Pulls", 100, -1., 1., 100, -4.5, 4.5);
1662 h->GetXaxis()->SetTitle("tg(#theta)");
1663 h->GetYaxis()->SetTitle("#Delta z / #sigma_{z}");
1664 h->GetZaxis()->SetTitle("entries");
1671 //________________________________________________________
1672 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
1675 // Do the processing
1678 Char_t pn[10]; sprintf(pn, "p%02d", fIdxPlot);
1680 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
1681 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
1683 for(Int_t ibin = 1; ibin <= h2->GetNbinsX(); ibin++){
1684 Double_t x = h2->GetXaxis()->GetBinCenter(ibin);
1685 TH1D *h = h2->ProjectionY(pn, ibin, ibin);
1686 if(h->GetEntries()<100) continue;
1691 Int_t ip = g[0]->GetN();
1692 g[0]->SetPoint(ip, x, k*f->GetParameter(1));
1693 g[0]->SetPointError(ip, 0., k*f->GetParError(1));
1694 g[1]->SetPoint(ip, x, k*f->GetParameter(2));
1695 g[1]->SetPointError(ip, 0., k*f->GetParError(2));
1698 g[0]->SetPoint(ip, x, k*h->GetMean());
1699 g[0]->SetPointError(ip, 0., k*h->GetMeanError());
1700 g[1]->SetPoint(ip, x, k*h->GetRMS());
1701 g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
1707 //________________________________________________________
1708 Bool_t AliTRDresolution::Process2D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1711 // Do the processing
1714 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1716 // retrive containers
1717 TH2I *h2 = idx<0 ? (TH2I*)(fContainer->At(plot)) : (TH2I*)((TObjArray*)(fContainer->At(plot)))->At(idx);
1718 if(!h2) return kFALSE;
1720 if(!(g[0] = idx<0 ? (TGraphErrors*)fGraphM->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
1722 if(!(g[1] = idx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
1724 return Process(h2, f, k, g);
1727 //________________________________________________________
1728 Bool_t AliTRDresolution::Process3D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1731 // Do the processing
1734 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1736 // retrive containers
1737 TH3S *h3 = idx<0 ? (TH3S*)(fContainer->At(plot)) : (TH3S*)((TObjArray*)(fContainer->At(plot)))->At(idx);
1738 if(!h3) return kFALSE;
1741 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
1742 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
1745 TAxis *az = h3->GetZaxis();
1746 for(Int_t iz=1; iz<=az->GetNbins(); iz++){
1747 if(!(g[0] = (TGraphErrors*)gm->At(iz-1))) return kFALSE;
1748 if(!(g[1] = (TGraphErrors*)gs->At(iz-1))) return kFALSE;
1749 az->SetRange(iz, iz);
1750 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
1756 //________________________________________________________
1757 Bool_t AliTRDresolution::Process3DL(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1760 // Do the processing
1763 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1765 // retrive containers
1766 TH3S *h3 = (TH3S*)((TObjArray*)fContainer->At(plot))->At(idx);
1767 if(!h3) return kFALSE;
1769 TGraphAsymmErrors *gm;
1771 if(!(gm = (TGraphAsymmErrors*)((TObjArray*)fGraphM->At(plot))->At(0))) return kFALSE;
1772 if(!(gs = (TGraphErrors*)((TObjArray*)fGraphS->At(plot)))) return kFALSE;
1774 Float_t x, r, mpv, xM, xm;
1775 TAxis *ay = h3->GetYaxis();
1776 for(Int_t iy=1; iy<=h3->GetNbinsY(); iy++){
1777 ay->SetRange(iy, iy);
1778 x = ay->GetBinCenter(iy);
1779 TH2F *h2=(TH2F*)h3->Project3D("zx");
1780 TAxis *ax = h2->GetXaxis();
1781 for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
1782 r = ax->GetBinCenter(ix);
1783 TH1D *h1 = h2->ProjectionY("py", ix, ix);
1784 if(h1->Integral()<50) continue;
1787 GetLandauMpvFwhm(f, mpv, xm, xM);
1788 Int_t ip = gm->GetN();
1789 gm->SetPoint(ip, x, k*mpv);
1790 gm->SetPointError(ip, 0., 0., k*xm, k*xM);
1791 gs->SetPoint(ip, r, k*(xM-xm)/mpv);
1792 gs->SetPointError(ip, 0., 0.);
1799 //________________________________________________________
1800 Bool_t AliTRDresolution::Process4D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
1803 // Do the processing
1806 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
1808 // retrive containers
1809 TObjArray *arr = (TObjArray*)((TObjArray*)(fContainer->At(plot)))->At(idx);
1810 if(!arr) return kFALSE;
1813 TObjArray *gm[2], *gs[2];
1814 if(!(gm[0] = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
1815 if(!(gs[0] = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
1820 for(Int_t ix=0; ix<arr->GetEntriesFast(); ix++){
1821 if(!(h3 = (TH3S*)arr->At(ix))) return kFALSE;
1822 if(!(gm[1] = (TObjArray*)gm[0]->At(ix))) return kFALSE;
1823 if(!(gs[1] = (TObjArray*)gs[0]->At(ix))) return kFALSE;
1824 TAxis *az = h3->GetZaxis();
1825 for(Int_t iz=1; iz<=az->GetNbins(); iz++){
1826 if(!(g[0] = (TGraphErrors*)gm[1]->At(iz-1))) return kFALSE;
1827 if(!(g[1] = (TGraphErrors*)gs[1]->At(iz-1))) return kFALSE;
1828 az->SetRange(iz, iz);
1829 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
1836 //________________________________________________________
1837 Bool_t AliTRDresolution::GetGraphPlot(Float_t *bb, ETRDresolutionPlot ip, Int_t idx)
1843 if(!fGraphS || !fGraphM) return kFALSE;
1844 TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphM->At(ip)))->At(idx);
1845 if(!gm) return kFALSE;
1846 TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphS->At(ip)))->At(idx);
1847 if(!gs) return kFALSE;
1848 gs->Draw("apl"); gm->Draw("pl");
1852 for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fgNElements[jp];
1853 UChar_t jdx = idx<0?0:idx;
1854 for(Int_t jc=0; jc<TMath::Min(jdx,fgNElements[ip]-1); jc++) nref++;
1855 const Char_t **at = fgAxTitle[nref];
1858 if((n=gm->GetN())) {
1859 PutTrendValue(Form("%s_%s", fgPerformanceName[ip], at[0]), gm->GetMean(2));
1860 PutTrendValue(Form("%s_%sRMS", fgPerformanceName[ip], at[0]), gm->GetRMS(2));
1864 gs->Sort(&TGraph::CompareY);
1865 PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[ip], at[0]), gs->GetY()[0]);
1866 PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[ip], at[0]), gs->GetY()[n-1]);
1867 gs->Sort(&TGraph::CompareX);
1872 ax = gs->GetHistogram()->GetXaxis();
1873 ax->SetRangeUser(bb[0], bb[2]);
1874 ax->SetTitle(at[1]);ax->CenterTitle();
1876 ax = gs->GetHistogram()->GetYaxis();
1877 ax->SetRangeUser(bb[1], bb[3]);
1878 ax->SetTitleOffset(1.1);
1879 ax->SetTitle(at[2]);ax->CenterTitle();
1882 gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
1883 gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
1884 //gax->SetVertical();
1885 gax->CenterTitle(); gax->SetTitleOffset(.7);
1886 gax->SetTitle(at[3]); gax->Draw();
1889 TBox *b = new TBox(-.15, bb[1], .15, bb[3]);
1890 b->SetFillStyle(3002);b->SetLineColor(0);
1891 b->SetFillColor(ip<=Int_t(kMCcluster)?kGreen:kBlue);
1897 //________________________________________________________
1898 Bool_t AliTRDresolution::GetGraphTrack(Float_t *bb, Int_t idx, Int_t il)
1904 if(!fGraphS || !fGraphM) return kFALSE;
1906 // axis titles look up
1908 for(Int_t jp=0; jp<Int_t(kMCtrackTRD); jp++) nref+=fgNElements[jp];
1909 for(Int_t jc=0; jc<idx; jc++) nref++;
1910 const Char_t **at = fgAxTitle[nref];
1912 TGraphErrors *gm = NULL, *gs = NULL;
1913 TObjArray *a0 = fGraphS, *a1 = NULL;
1914 a1 = (TObjArray*)a0->At(kMCtrackTRD); a0 = a1;
1915 a1 = (TObjArray*)a0->At(idx); a0 = a1;
1916 a1 = (TObjArray*)a0->At(il); a0 = a1;
1917 for(Int_t is=0; is<AliPID::kSPECIES; is++){
1918 if(!(gs = (TGraphErrors*)a0->At(is))) return kFALSE;
1919 if(!gs->GetN()) continue;
1920 gs->Draw(is ? "pl" : "apl");
1921 gs->Sort(&TGraph::CompareY); Int_t n = gs->GetN();
1922 PutTrendValue(Form("%s_%sSigMin%s", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is)), gs->GetY()[0]);
1923 PutTrendValue(Form("%s_%sSigMax%s", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is)), gs->GetY()[n-1]);
1924 gs->Sort(&TGraph::CompareX);
1926 gs = (TGraphErrors*)a0->At(0);
1929 TAxis *ax = gs->GetHistogram()->GetXaxis();
1930 ax->SetRangeUser(bb[0], bb[2]);
1931 ax->SetTitle(at[1]);ax->CenterTitle();
1933 ax = gs->GetHistogram()->GetYaxis();
1934 ax->SetRangeUser(bb[1], bb[3]);
1935 ax->SetTitleOffset(.5);ax->SetTitleSize(.06);
1936 ax->SetTitle(at[2]);ax->CenterTitle();
1939 gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
1940 gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
1941 //gax->SetVertical();
1942 gax->CenterTitle(); gax->SetTitleOffset(.5);gax->SetTitleSize(.06);
1943 gax->SetTitle(at[3]); gax->Draw();
1947 a1 = (TObjArray*)a0->At(kMCtrackTRD); a0 = a1;
1948 a1 = (TObjArray*)a0->At(idx); a0 = a1;
1949 a1 = (TObjArray*)a0->At(il); a0 = a1;
1950 for(Int_t is=0; is<AliPID::kSPECIES; is++){
1951 if(!(gm = (TGraphErrors*)a0->At(is))) return kFALSE;
1952 if(!gm->GetN()) continue;
1954 PutTrendValue(Form("%s_%s_%s", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is)), gm->GetMean(2));
1955 PutTrendValue(Form("%s_%s_%sRMS", fgPerformanceName[kMCtrackTRD], at[0], AliPID::ParticleShortName(is)), gm->GetRMS(2));
1962 //________________________________________________________
1963 Bool_t AliTRDresolution::GetGraphTrackTPC(Float_t *bb, Int_t sel)
1969 if(!fGraphS || !fGraphM) return kFALSE;
1971 // axis titles look up
1973 for(Int_t jp=0; jp<Int_t(kMCtrackTPC); jp++) nref+=fgNElements[jp];
1974 for(Int_t jc=0; jc<sel; jc++) nref++;
1975 const Char_t **at = fgAxTitle[nref];
1977 TGraphErrors *gm = NULL, *gs = NULL;
1978 TObjArray *a0 = fGraphS, *a1 = NULL;
1979 a1 = (TObjArray*)a0->At(kMCtrackTPC); a0 = a1;
1980 a1 = (TObjArray*)a0->At(sel); a0 = a1;
1981 for(Int_t is=AliPID::kSPECIES; is--;){
1982 if(!(gs = (TGraphErrors*)a0->At(is))) return kFALSE;
1983 if(!gs->GetN()) continue;
1984 gs->Draw(is ? "pl" : "apl");
1985 gs->Sort(&TGraph::CompareY); Int_t n = gs->GetN();
1986 PutTrendValue(Form("%s_%sSigMin%s", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is)), gs->GetY()[0]);
1987 PutTrendValue(Form("%s_%sSigMax%s", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is)), gs->GetY()[n-1]);
1988 gs->Sort(&TGraph::CompareX);
1990 gs = (TGraphErrors*)a0->At(0);
1992 TAxis *ax = gs->GetHistogram()->GetXaxis();
1993 ax->SetRangeUser(bb[0], bb[2]);
1994 ax->SetTitle(at[1]);ax->CenterTitle();
1996 ax = gs->GetHistogram()->GetYaxis();
1997 ax->SetRangeUser(bb[1], bb[3]);
1998 ax->SetTitleOffset(1.);ax->SetTitleSize(0.05);
1999 ax->SetTitle(at[2]);ax->CenterTitle();
2002 gax = new TGaxis(bb[2], bb[1], bb[2], bb[3], bb[1], bb[3], 510, "+U");
2003 gax->SetLineColor(kRed);gax->SetLineWidth(2);gax->SetTextColor(kRed);
2004 //gax->SetVertical();
2005 gax->CenterTitle(); gax->SetTitleOffset(.7);gax->SetTitleSize(0.05);
2006 gax->SetTitle(at[3]); gax->Draw();
2010 a1 = (TObjArray*)a0->At(kMCtrackTPC); a0 = a1;
2011 a1 = (TObjArray*)a0->At(sel); a0 = a1;
2012 for(Int_t is=AliPID::kSPECIES; is--;){
2013 if(!(gm = (TGraphErrors*)a0->At(is))) return kFALSE;
2014 if(!gm->GetN()) continue;
2016 PutTrendValue(Form("%s_%s_%s", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is)), gm->GetMean(2));
2017 PutTrendValue(Form("%s_%s_%sRMS", fgPerformanceName[kMCtrackTPC], at[0], AliPID::ParticleShortName(is)), gm->GetRMS(2));
2023 //________________________________________________________
2024 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
2027 // Get the most probable value and the full width half mean
2028 // of a Landau distribution
2031 const Float_t dx = 1.;
2032 mpv = f->GetParameter(1);
2033 Float_t fx, max = f->Eval(mpv);
2036 while((fx = f->Eval(xm))>.5*max){
2045 while((fx = f->Eval(xM))>.5*max) xM += dx;
2049 //________________________________________________________
2050 void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
2053 fReconstructor->SetRecoParam(r);