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 ////////////////////////////////////////////////////////////////////////////
54 #include <TObjArray.h>
63 #include <TGraphErrors.h>
64 #include <TGraphAsymmErrors.h>
68 #include <TTreeStream.h>
69 #include <TGeoManager.h>
70 #include <TDatabasePDG.h>
74 #include "AliESDtrack.h"
75 #include "AliMathBase.h"
77 #include "AliTRDresolution.h"
78 #include "AliTRDgeometry.h"
79 #include "AliTRDpadPlane.h"
80 #include "AliTRDcluster.h"
81 #include "AliTRDseedV1.h"
82 #include "AliTRDtrackV1.h"
83 #include "AliTRDReconstructor.h"
84 #include "AliTRDrecoParam.h"
85 #include "AliTRDpidUtil.h"
86 #include "AliTRDinfoGen.h"
88 #include "info/AliTRDclusterInfo.h"
90 ClassImp(AliTRDresolution)
92 UChar_t const AliTRDresolution::fgNproj[kNviews] = {
96 Char_t const * AliTRDresolution::fgPerformanceName[kNviews] = {
108 Char_t const * AliTRDresolution::fgParticle[11]={
109 " p bar", " K -", " #pi -", " #mu -", " e -",
111 " e +", " #mu +", " #pi +", " K +", " p",
114 // Configure segmentation for y resolution/residuals
115 Int_t const AliTRDresolution::fgkNresYsegm[3] = {
116 AliTRDgeometry::kNsector
117 ,AliTRDgeometry::kNsector*AliTRDgeometry::kNstack
118 ,AliTRDgeometry::kNdet
120 Char_t const *AliTRDresolution::fgkResYsegmName[3] = {
121 "Sector", "Stack", "Detector"};
124 //________________________________________________________
125 AliTRDresolution::AliTRDresolution()
140 // Default constructor
142 SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
143 SetSegmentationLevel();
146 //________________________________________________________
147 AliTRDresolution::AliTRDresolution(char* name)
148 :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
162 // Default constructor
166 SetSegmentationLevel();
168 DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
169 DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
170 /* DefineOutput(kTrkltToTrk, TObjArray::Class()); // tracklet2track
171 DefineOutput(kTrkltToMC, TObjArray::Class()); // tracklet2mc*/
174 //________________________________________________________
175 AliTRDresolution::~AliTRDresolution()
181 if(fGraphS){fGraphS->Delete(); delete fGraphS;}
182 if(fGraphM){fGraphM->Delete(); delete fGraphM;}
183 if(fCl){fCl->Delete(); delete fCl;}
184 if(fMCcl){fMCcl->Delete(); delete fMCcl;}
185 /* if(fTrklt){fTrklt->Delete(); delete fTrklt;}
186 if(fMCtrklt){fMCtrklt->Delete(); delete fMCtrklt;}*/
190 //________________________________________________________
191 void AliTRDresolution::UserCreateOutputObjects()
193 // spatial resolution
195 AliTRDrecoTask::UserCreateOutputObjects();
196 InitExchangeContainers();
199 //________________________________________________________
200 void AliTRDresolution::InitExchangeContainers()
202 // Init containers for subsequent tasks (AliTRDclusterResolution)
204 fCl = new TObjArray(200);
205 fCl->SetOwner(kTRUE);
206 fMCcl = new TObjArray();
207 fMCcl->SetOwner(kTRUE);
208 /* fTrklt = new TObjArray();
209 fTrklt->SetOwner(kTRUE);
210 fMCtrklt = new TObjArray();
211 fMCtrklt->SetOwner(kTRUE);*/
212 PostData(kClToTrk, fCl);
213 PostData(kClToMC, fMCcl);
214 /* PostData(kTrkltToTrk, fTrklt);
215 PostData(kTrkltToMC, fMCtrklt);*/
218 //________________________________________________________
219 void AliTRDresolution::UserExec(Option_t *opt)
228 fMCtrklt->Delete();*/
229 AliTRDrecoTask::UserExec(opt);
232 //________________________________________________________
233 Bool_t AliTRDresolution::Pulls(Double_t dyz[2], Double_t cov[3], Double_t tilt) const
235 // Helper function to calculate pulls in the yz plane
236 // using proper tilt rotation
237 // Uses functionality defined by AliTRDseedV1.
239 Double_t t2(tilt*tilt);
243 cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
244 cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
245 cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
247 Double_t sqr[3]={0., 0., 0.};
248 if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
249 Double_t invsqr[3]={0., 0., 0.};
250 if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
251 Double_t tmp(dyz[0]);
252 dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
253 dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
257 //________________________________________________________
258 TH1* AliTRDresolution::PlotCharge(const AliTRDtrackV1 *track)
261 // Plots the charge distribution
264 if(track) fkTrack = track;
266 AliDebug(4, "No Track defined.");
269 TObjArray *arr = NULL;
270 if(!fContainer || !(arr = ((TObjArray*)fContainer->At(kCharge)))){
271 AliWarning("No output container defined.");
276 AliTRDseedV1 *fTracklet = NULL;
277 AliTRDcluster *c = NULL;
278 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
279 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
280 if(!fTracklet->IsOK()) continue;
281 Float_t x0 = fTracklet->GetX0();
283 for(Int_t itb=AliTRDseedV1::kNtb; itb--;){
284 if(!(c = fTracklet->GetClusters(itb))){
285 if(!(c = fTracklet->GetClusters(AliTRDseedV1::kNtb+itb))) continue;
287 dqdl = fTracklet->GetdQdl(itb, &dl);
288 if(dqdl<1.e-5) continue;
289 dl /= 0.15; // dl/dl0, dl0 = 1.5 mm for nominal vd
290 (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dqdl);
293 // if(!HasMCdata()) continue;
295 // Float_t pt0, y0, z0, dydx0, dzdx0;
296 // if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
303 //________________________________________________________
304 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
307 // Plot the cluster distributions
310 if(track) fkTrack = track;
312 AliDebug(4, "No Track defined.");
315 TObjArray *arr = NULL;
316 if(!fContainer || !(arr = ((TObjArray*)fContainer->At(kCluster)))){
317 AliWarning("No output container defined.");
320 ULong_t status = fkESD ? fkESD->GetStatus():0;
323 Double_t covR[7], cov[3], dy[2], dz[2];
324 Float_t pt, x0, y0, z0, dydx, dzdx;
325 const AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
326 AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
327 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
328 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
329 if(!fTracklet->IsOK()) continue;
330 x0 = fTracklet->GetX0();
331 pt = fTracklet->GetPt();
332 sgm[2] = fTracklet->GetDetector();
333 sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
334 sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
336 // retrive the track angle with the chamber
337 y0 = fTracklet->GetYref(0);
338 z0 = fTracklet->GetZref(0);
339 dydx = fTracklet->GetYref(1);
340 dzdx = fTracklet->GetZref(1);
341 fTracklet->GetCovRef(covR);
342 Double_t tilt(fTracklet->GetTilt())
345 ,cost(TMath::Sqrt(corr));
346 AliTRDcluster *c = NULL;
347 fTracklet->ResetClusterIter(kFALSE);
348 while((c = fTracklet->PrevCluster())){
349 Float_t xc = c->GetX();
350 Float_t yc = c->GetY();
351 Float_t zc = c->GetZ();
352 Float_t dx = x0 - xc;
353 Float_t yt = y0 - dx*dydx;
354 Float_t zt = z0 - dx*dzdx;
355 dy[0] = yc-yt; dz[0]= zc-zt;
358 dy[1] = cost*(dy[0] - dz[0]*tilt);
359 dz[1] = cost*(dz[0] + dy[0]*tilt);
360 if(pt>fPtThreshold && c->IsInChamber()) ((TH3S*)arr->At(0))->Fill(dydx, dy[1], sgm[fSegmentLevel]);
362 // tilt rotation of covariance for clusters
363 Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
364 cov[0] = (sy2+t2*sz2)*corr;
365 cov[1] = tilt*(sz2 - sy2)*corr;
366 cov[2] = (t2*sy2+sz2)*corr;
367 // sum with track covariance
368 cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
369 Double_t dyz[2]= {dy[1], dz[1]};
370 Pulls(dyz, cov, tilt);
371 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
373 // Get z-position with respect to anode wire
374 Int_t istk = geo->GetStack(c->GetDetector());
375 AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
376 Float_t row0 = pp->GetRow0();
377 Float_t d = row0 - zt + pp->GetAnodeWireOffset();
378 d -= ((Int_t)(2 * d)) / 2.0;
379 if (d > 0.25) d = 0.5 - d;
381 AliTRDclusterInfo *clInfo(NULL);
382 clInfo = new AliTRDclusterInfo;
383 clInfo->SetCluster(c);
384 Float_t covF[] = {cov[0], cov[1], cov[2]};
385 clInfo->SetGlobalPosition(yt, zt, dydx, dzdx, covF);
386 clInfo->SetResolution(dy[1]);
387 clInfo->SetAnisochronity(d);
388 clInfo->SetDriftLength(dx);
389 clInfo->SetTilt(tilt);
390 if(fCl) fCl->Add(clInfo);
391 else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
395 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
396 clInfoArr->SetOwner(kFALSE);
398 clInfoArr->Add(clInfo);
401 if(DebugLevel()>=1 && clInfoArr){
402 (*DebugStream()) << "cluster"
403 <<"status=" << status
404 <<"clInfo.=" << clInfoArr
409 if(clInfoArr) delete clInfoArr;
410 return (TH3S*)arr->At(0);
414 //________________________________________________________
415 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
417 // Plot normalized residuals for tracklets to track.
419 // We start from the result that if X=N(|m|, |Cov|)
421 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
423 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
424 // reference position.
425 if(track) fkTrack = track;
427 AliDebug(4, "No Track defined.");
430 TObjArray *arr = NULL;
431 if(!fContainer || !(arr = (TObjArray*)fContainer->At(kTrack ))){
432 AliWarning("No output container defined.");
437 Double_t cov[3], covR[7]/*, sqr[3], inv[3]*/;
438 Double_t pt, phi, tht, x, dx, dy[2], dz[2];
439 AliTRDseedV1 *fTracklet(NULL);
440 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
441 if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
442 if(!fTracklet->IsOK()) continue;
443 sgm[2] = fTracklet->GetDetector();
444 sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
445 sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
446 x = fTracklet->GetX();
447 dx = fTracklet->GetX0() - x;
448 pt = fTracklet->GetPt();
449 phi = fTracklet->GetYref(1);
450 tht = fTracklet->GetZref(1);
452 dy[0]= fTracklet->GetYref(0)-dx*fTracklet->GetYref(1) - fTracklet->GetY();
453 dz[0]= fTracklet->GetZref(0)-dx*fTracklet->GetZref(1) - fTracklet->GetZ();
454 Double_t tilt(fTracklet->GetTilt())
457 ,cost(TMath::Sqrt(corr));
458 Bool_t rc(fTracklet->IsRowCross());
460 // calculate residuals using tilt rotation
461 dy[1]= cost*(dy[0] - dz[0]*tilt);
462 dz[1]= cost*(dz[0] + dy[0]*tilt);
463 ((TH3S*)arr->At(0))->Fill(phi, dy[1], sgm[fSegmentLevel]+rc*fgkNresYsegm[fSegmentLevel]);
464 ((TH3S*)arr->At(2))->Fill(tht, dz[1], rc);
466 // compute covariance matrix
467 fTracklet->GetCovAt(x, cov);
468 fTracklet->GetCovRef(covR);
469 cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
470 Double_t dyz[2]= {dy[1], dz[1]};
471 Pulls(dyz, cov, tilt);
472 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
473 ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
475 Double_t dphi((phi-fTracklet->GetYfit(1))/(1-phi*fTracklet->GetYfit(1)));
476 Double_t dtht((tht-fTracklet->GetZfit(1))/(1-tht*fTracklet->GetZfit(1)));
477 ((TH2I*)arr->At(4))->Fill(phi, TMath::ATan(dphi));
480 UChar_t err(fTracklet->GetErrorMsg());
481 (*DebugStream()) << "tracklet"
501 return (TH2I*)arr->At(0);
505 //________________________________________________________
506 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
508 // Store resolution/pulls of Kalman before updating with the TRD information
509 // at the radial position of the first tracklet. The following points are used
511 // - the (y,z,snp) of the first TRD tracklet
512 // - the (y, z, snp, tgl, pt) of the MC track reference
514 // Additionally the momentum resolution/pulls are calculated for usage in the
517 if(track) fkTrack = track;
519 AliDebug(4, "No Track defined.");
522 TObjArray *arr = NULL;
523 if(!fContainer || !(arr = (TObjArray*)fContainer->At(kTrackIn))){
524 AliWarning("No output container defined.");
527 AliExternalTrackParam *tin = NULL;
528 if(!(tin = fkTrack->GetTrackIn())){
529 AliWarning("Track did not entered TRD fiducial volume.");
534 Double_t x = tin->GetX();
535 AliTRDseedV1 *fTracklet = NULL;
536 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
537 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
540 if(!fTracklet || TMath::Abs(x-fTracklet->GetX())>1.e-3){
541 AliWarning("Tracklet did not match Track.");
545 sgm[2] = fTracklet->GetDetector();
546 sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
547 sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
548 Double_t tilt(fTracklet->GetTilt())
551 ,cost(TMath::Sqrt(corr));
552 Bool_t rc(fTracklet->IsRowCross());
554 const Int_t kNPAR(5);
555 Double_t parR[kNPAR]; memcpy(parR, tin->GetParameter(), kNPAR*sizeof(Double_t));
556 Double_t covR[3*kNPAR]; memcpy(covR, tin->GetCovariance(), 3*kNPAR*sizeof(Double_t));
557 Double_t cov[3]; fTracklet->GetCovAt(x, cov);
559 // define sum covariances
560 TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
561 Double_t *pc = &covR[0], *pp = &parR[0];
562 for(Int_t ir=0; ir<kNPAR; ir++, pp++){
564 for(Int_t ic = 0; ic<=ir; ic++,pc++){
565 COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
568 PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
569 //COV.Print(); PAR.Print();
571 //TODO Double_t dydx = TMath::Sqrt(1.-parR[2]*parR[2])/parR[2];
572 Double_t dy[2]={parR[0] - fTracklet->GetY(), 0.}
573 ,dz[2]={parR[1] - fTracklet->GetZ(), 0.}
574 ,dphi(TMath::ASin(PAR[2])-TMath::ATan(fTracklet->GetYfit(1)));
575 // calculate residuals using tilt rotation
576 dy[1] = cost*(dy[0] - dz[0]*tilt);
577 dz[1] = cost*(dz[0] + dy[0]*tilt);
579 if(1./PAR[4]>fPtThreshold) ((TH3S*)arr->At(0))->Fill(fTracklet->GetYref(1), dy[1], sgm[fSegmentLevel]+rc*fgkNresYsegm[fSegmentLevel]);
580 ((TH3S*)arr->At(2))->Fill(fTracklet->GetZref(1), dz[1], rc);
581 ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), dphi);
583 Double_t dyz[2] = {dy[1], dz[1]};
584 Double_t cc[3] = {COV(0,0)+cov[0], COV(0,1)+cov[1], COV(1,1)+cov[2]};
585 Pulls(dyz, cc, tilt);
586 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
587 ((TH3S*)arr->At(3))->Fill(fTracklet->GetZref(1), dyz[1], rc);
591 // register reference histo for mini-task
592 h = (TH2I*)arr->At(0);
595 (*DebugStream()) << "trackIn"
601 Double_t y = fTracklet->GetY();
602 Double_t z = fTracklet->GetZ();
603 (*DebugStream()) << "trackletIn"
613 if(!HasMCdata()) return h;
615 Float_t dx, pt0, x0=fTracklet->GetX0(), y0, z0, dydx0, dzdx0;
616 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
617 Int_t pdg = fkMC->GetPDG(),
618 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
620 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
621 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
622 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
624 // translate to reference radial position
625 dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
626 Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi))
628 TVectorD PARMC(kNPAR);
629 PARMC[0]=y0; PARMC[1]=z0;
630 PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
633 // TMatrixDSymEigen eigen(COV);
634 // TVectorD evals = eigen.GetEigenValues();
635 // TMatrixDSym evalsm(kNPAR);
636 // for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
637 // TMatrixD evecs = eigen.GetEigenVectors();
638 // TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
641 if(!(arr = (TObjArray*)fContainer->At(kMCtrackIn))) {
642 AliWarning("No MC container defined.");
646 // y resolution/pulls
647 if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0], sgm[fSegmentLevel]);
648 ((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)));
649 // z resolution/pulls
650 ((TH3S*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1], 0);
651 ((TH3S*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)), 0);
652 // phi resolution/snp pulls
653 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
654 ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
655 // theta resolution/tgl pulls
656 ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
657 ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
658 // pt resolution\\1/pt pulls\\p resolution/pull
659 ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., sign*sIdx);
660 ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), sign*sIdx);
662 Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
663 p = TMath::Sqrt(1.+ PAR[3]*PAR[3])/PAR[4];
664 ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., sign*sIdx);
666 // p*p*PAR[4]*PAR[4]*COV(4,4)
667 // +2.*PAR[3]*COV(3,4)/PAR[4]
668 // +PAR[3]*PAR[3]*COV(3,3)/p/p/PAR[4]/PAR[4]/PAR[4]/PAR[4];
669 // if(sp>0.) ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/TMath::Sqrt(sp), sign*sIdx);
673 (*DebugStream()) << "trackInMC"
680 //________________________________________________________
681 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
683 // Store resolution/pulls of Kalman after last update with the TRD information
684 // at the radial position of the first tracklet. The following points are used
686 // - the (y,z,snp) of the first TRD tracklet
687 // - the (y, z, snp, tgl, pt) of the MC track reference
689 // Additionally the momentum resolution/pulls are calculated for usage in the
692 if(track) fkTrack = track;
694 AliDebug(4, "No Track defined.");
697 TObjArray *arr = NULL;
698 if(!fContainer || !(arr = (TObjArray*)fContainer->At(kTrackOut))){
699 AliWarning("No output container defined.");
702 AliExternalTrackParam *tout = NULL;
703 if(!(tout = fkTrack->GetTrackOut())){
704 AliDebug(2, "Track did not exit TRD.");
709 Double_t x = tout->GetX();
710 AliTRDseedV1 *fTracklet(NULL);
711 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
712 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
715 if(!fTracklet || TMath::Abs(x-fTracklet->GetX())>1.e-3){
716 AliWarning("Tracklet did not match Track position.");
720 sgm[2] = fTracklet->GetDetector();
721 sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
722 sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
723 Double_t tilt(fTracklet->GetTilt())
726 ,cost(TMath::Sqrt(corr));
727 Bool_t rc(fTracklet->IsRowCross());
729 const Int_t kNPAR(5);
730 Double_t parR[kNPAR]; memcpy(parR, tout->GetParameter(), kNPAR*sizeof(Double_t));
731 Double_t covR[3*kNPAR]; memcpy(covR, tout->GetCovariance(), 3*kNPAR*sizeof(Double_t));
732 Double_t cov[3]; fTracklet->GetCovAt(x, cov);
734 // define sum covariances
735 TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
736 Double_t *pc = &covR[0], *pp = &parR[0];
737 for(Int_t ir=0; ir<kNPAR; ir++, pp++){
739 for(Int_t ic = 0; ic<=ir; ic++,pc++){
740 COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
743 PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
744 //COV.Print(); PAR.Print();
746 //TODO Double_t dydx = TMath::Sqrt(1.-parR[2]*parR[2])/parR[2];
747 Double_t dy[3]={parR[0] - fTracklet->GetY(), 0., 0.}
748 ,dz[3]={parR[1] - fTracklet->GetZ(), 0., 0.}
749 ,dphi(TMath::ASin(PAR[2])-TMath::ATan(fTracklet->GetYfit(1)));
750 // calculate residuals using tilt rotation
751 dy[1] = cost*(dy[0] - dz[0]*tilt);
752 dz[1] = cost*(dz[0] + dy[0]*tilt);
754 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 !!!
755 ((TH3S*)arr->At(2))->Fill(fTracklet->GetZref(1), dz[1], rc);
756 ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), dphi);
758 Double_t dyz[2] = {dy[1], dz[1]};
759 Double_t cc[3] = {COV(0,0)+cov[0], COV(0,1)+cov[1], COV(1,1)+cov[2]};
760 Pulls(dyz, cc, tilt);
761 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
762 ((TH3S*)arr->At(3))->Fill(fTracklet->GetZref(1), dyz[1], rc);
764 // register reference histo for mini-task
765 h = (TH2I*)arr->At(0);
768 (*DebugStream()) << "trackOut"
774 Double_t y = fTracklet->GetY();
775 Double_t z = fTracklet->GetZ();
776 (*DebugStream()) << "trackletOut"
786 if(!HasMCdata()) return h;
788 Float_t dx, pt0, x0=fTracklet->GetX0(), y0, z0, dydx0, dzdx0;
789 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
790 Int_t pdg = fkMC->GetPDG(),
791 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
793 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
794 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
795 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
797 // translate to reference radial position
798 dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
799 Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi))
801 TVectorD PARMC(kNPAR);
802 PARMC[0]=y0; PARMC[1]=z0;
803 PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
806 // TMatrixDSymEigen eigen(COV);
807 // TVectorD evals = eigen.GetEigenValues();
808 // TMatrixDSym evalsm(kNPAR);
809 // for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
810 // TMatrixD evecs = eigen.GetEigenVectors();
811 // TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
814 if(!(arr = (TObjArray*)fContainer->At(kMCtrackOut))){
815 AliWarning("No MC container defined.");
818 // y resolution/pulls
819 if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0], sgm[fSegmentLevel]);
820 ((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)));
821 // z resolution/pulls
822 ((TH3S*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1], 0);
823 ((TH3S*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)), 0);
824 // phi resolution/snp pulls
825 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
826 ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
827 // theta resolution/tgl pulls
828 ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
829 ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
830 // pt resolution\\1/pt pulls\\p resolution/pull
831 ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., sign*sIdx);
832 ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), sign*sIdx);
834 Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
835 p = TMath::Sqrt(1.+ PAR[3]*PAR[3])/PAR[4];
836 ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., sign*sIdx);
838 // p*p*PAR[4]*PAR[4]*COV(4,4)
839 // +2.*PAR[3]*COV(3,4)/PAR[4]
840 // +PAR[3]*PAR[3]*COV(3,3)/p/p/PAR[4]/PAR[4]/PAR[4]/PAR[4];
841 // if(sp>0.) ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/TMath::Sqrt(sp), sign*sIdx);
845 (*DebugStream()) << "trackOutMC"
852 //________________________________________________________
853 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
856 // Plot MC distributions
860 AliDebug(2, "No MC defined. Results will not be available.");
863 if(track) fkTrack = track;
865 AliDebug(4, "No Track defined.");
869 AliWarning("No output container defined.");
872 // retriev track characteristics
873 Int_t pdg = fkMC->GetPDG(),
874 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
877 label(fkMC->GetLabel());
878 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
879 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
880 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
882 TObjArray *arr(NULL);TH1 *h(NULL);
884 Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl;
885 Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
886 Double_t covR[7]/*, cov[3]*/;
889 TVectorD dX(12), dY(12), dZ(12), vPt(12), dPt(12), cCOV(12*15);
890 fkMC->PropagateKalman(&dX, &dY, &dZ, &vPt, &dPt, &cCOV);
891 (*DebugStream()) << "MCkalman"
901 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
902 AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
903 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
904 if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
905 !fTracklet->IsOK())*/ continue;
907 sgm[2] = fTracklet->GetDetector();
908 sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
909 sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
910 Double_t tilt(fTracklet->GetTilt())
913 ,cost(TMath::Sqrt(corr));
914 x0 = fTracklet->GetX0();
915 //radial shift with respect to the MC reference (radial position of the pad plane)
916 x= fTracklet->GetX();
917 Bool_t rc(fTracklet->IsRowCross());
918 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
919 xAnode = fTracklet->GetX0();
921 // MC track position at reference radial position
924 (*DebugStream()) << "MC"
936 Float_t ymc = y0 - dx*dydx0;
937 Float_t zmc = z0 - dx*dzdx0;
938 //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
940 // Kalman position at reference radial position
942 dydx = fTracklet->GetYref(1);
943 dzdx = fTracklet->GetZref(1);
944 dzdl = fTracklet->GetTgl();
945 y = fTracklet->GetYref(0) - dx*dydx;
947 z = fTracklet->GetZref(0) - dx*dzdx;
949 pt = TMath::Abs(fTracklet->GetPt());
950 fTracklet->GetCovRef(covR);
952 arr = (TObjArray*)((TObjArray*)fContainer->At(kMCtrack))->At(ily);
953 // y resolution/pulls
954 if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
955 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(covR[0]), dz/TMath::Sqrt(covR[2]));
956 // z resolution/pulls
957 ((TH3S*)arr->At(2))->Fill(dzdx0, dz, 0);
958 ((TH3S*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]), 0);
959 // phi resolution/ snp pulls
960 Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
961 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
962 Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0);
963 ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3]));
964 // theta resolution/ tgl pulls
965 Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
966 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
967 ((TH2I*)arr->At(6))->Fill(dzdl0,
969 ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
970 // pt resolution \\ 1/pt pulls \\ p resolution for PID
971 Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
972 p = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
973 ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
974 ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
975 ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);
977 // Fill Debug stream for Kalman track
979 (*DebugStream()) << "MCtrack"
991 // recalculate tracklet based on the MC info
992 AliTRDseedV1 tt(*fTracklet);
993 tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
994 tt.SetZref(1, dzdx0);
995 tt.SetReconstructor(AliTRDinfoGen::Reconstructor());
997 x= tt.GetX();y= tt.GetY();z= tt.GetZ();
998 dydx = tt.GetYfit(1);
1000 ymc = y0 - dx*dydx0;
1001 zmc = z0 - dx*dzdx0;
1004 Float_t dphi = (dydx - dydx0);
1005 dphi /= (1.- dydx*dydx0);
1007 // add tracklet residuals for y and dydx
1008 arr = (TObjArray*)fContainer->At(kMCtracklet);
1010 if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
1011 if(tt.GetS2Y()>0. && tt.GetS2Z()>0.) ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(tt.GetS2Y()), dz/TMath::Sqrt(tt.GetS2Z()));
1012 ((TH3S*)arr->At(2))->Fill(dzdl0, dz, rc);
1013 if(tt.GetS2Z()>0.) ((TH3S*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z()), rc);
1014 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
1016 // Fill Debug stream for tracklet
1017 if(DebugLevel()>=4){
1018 Float_t s2y = tt.GetS2Y();
1019 Float_t s2z = tt.GetS2Z();
1020 (*DebugStream()) << "MCtracklet"
1031 AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(sgm[2]));
1032 Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
1033 //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
1035 arr = (TObjArray*)fContainer->At(kMCcluster);
1036 AliTRDcluster *c = NULL;
1037 tt.ResetClusterIter(kFALSE);
1038 while((c = tt.PrevCluster())){
1039 Float_t q = TMath::Abs(c->GetQ());
1040 x = c->GetX(); y = c->GetY();z = c->GetZ();
1044 dy = cost*(y - ymc - tilt*(z-zmc));
1045 dz = cost*(z - zmc + tilt*(y-ymc));
1048 if(q>20. && q<250. && pt0>fPtThreshold && c->IsInChamber()){
1049 ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
1050 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(c->GetSigmaY2()), dz/TMath::Sqrt(c->GetSigmaZ2()));
1053 // Fill calibration container
1054 Float_t d = zr0 - zmc;
1055 d -= ((Int_t)(2 * d)) / 2.0;
1056 if (d > 0.25) d = 0.5 - d;
1057 AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
1058 clInfo->SetCluster(c);
1059 clInfo->SetMC(pdg, label);
1060 clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
1061 clInfo->SetResolution(dy);
1062 clInfo->SetAnisochronity(d);
1063 clInfo->SetDriftLength(dx);
1064 clInfo->SetTilt(tilt);
1065 if(fMCcl) fMCcl->Add(clInfo);
1066 else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
1067 if(DebugLevel()>=5){
1069 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
1070 clInfoArr->SetOwner(kFALSE);
1072 clInfoArr->Add(clInfo);
1076 if(DebugLevel()>=5 && clInfoArr){
1077 (*DebugStream()) << "MCcluster"
1078 <<"clInfo.=" << clInfoArr
1083 if(clInfoArr) delete clInfoArr;
1088 //________________________________________________________
1089 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
1092 // Get the reference figures
1095 Float_t xy[4] = {0., 0., 0., 0.};
1097 AliWarning("Please provide a canvas to draw results.");
1100 Int_t selection[100], n(0), selStart(0); //
1101 Int_t ly0(0), dly(5);
1102 //Int_t ly0(1), dly(2); // used for SA
1103 TList *l(NULL); TVirtualPad *pad(NULL);
1104 TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);
1106 case 0: // charge resolution
1107 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1108 ((TVirtualPad*)l->At(0))->cd();
1109 ga=((TGraphAsymmErrors*)((TObjArray*)fGraphM->At(kCharge))->At(0));
1110 if(ga->GetN()) ga->Draw("apl");
1111 ((TVirtualPad*)l->At(1))->cd();
1112 g = ((TGraphErrors*)((TObjArray*)fGraphS->At(kCharge))->At(0));
1113 if(g->GetN()) g->Draw("apl");
1115 case 1: // cluster2track residuals
1116 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1117 xy[0] = -.3; xy[1] = -100.; xy[2] = .3; xy[3] = 1000.;
1118 pad = (TVirtualPad*)l->At(0); pad->cd();
1119 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1120 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1121 if(!GetGraphArray(xy, kCluster, 0, 1, n, selection)) break;
1122 pad=(TVirtualPad*)l->At(1); pad->cd();
1123 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1124 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1125 if(!GetGraphArray(xy, kCluster, 0, 1, n, selection)) break;
1127 case 2: // cluster2track residuals
1128 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1129 xy[0] = -.3; xy[1] = -100.; xy[2] = .3; xy[3] = 1000.;
1130 pad = (TVirtualPad*)l->At(0); pad->cd();
1131 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1132 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1133 if(!GetGraphArray(xy, kCluster, 0, 1, n, selection)) break;
1134 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-0.5; xy[3] = 2.5;
1135 pad=(TVirtualPad*)l->At(1); pad->cd();
1136 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1137 if(!GetGraphArray(xy, kCluster, 1, 1)) break;
1140 gPad->Divide(3, 2, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1141 xy[0] = -.3; xy[1] = -20.; xy[2] = .3; xy[3] = 100.;
1142 ((TVirtualPad*)l->At(0))->cd();
1143 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1144 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection)) break;
1146 ((TVirtualPad*)l->At(1))->cd();
1147 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1148 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection)) break;
1150 ((TVirtualPad*)l->At(2))->cd();
1151 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1152 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection)) break;
1154 ((TVirtualPad*)l->At(3))->cd();
1155 selStart=fgkNresYsegm[fSegmentLevel]; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1156 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection, "[RC]")) break;
1158 ((TVirtualPad*)l->At(4))->cd();
1159 selStart=fgkNresYsegm[fSegmentLevel]/3+fgkNresYsegm[fSegmentLevel]; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1160 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection, "[RC]")) break;
1162 ((TVirtualPad*)l->At(5))->cd();
1163 selStart=2*fgkNresYsegm[fSegmentLevel]/3+fgkNresYsegm[fSegmentLevel]; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1164 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection, "[RC]")) break;
1167 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1169 xy[0] = -1.; xy[1] = -150.; xy[2] = 1.; xy[3] = 1000.;
1170 ((TVirtualPad*)l->At(0))->cd();
1172 if(!GetGraphArray(xy, kTrack, 2, 1, 1, selection)) break;
1174 xy[0] = -1.; xy[1] = -1500.; xy[2] = 1.; xy[3] = 10000.;
1175 ((TVirtualPad*)l->At(1))->cd();
1177 if(!GetGraphArray(xy, kTrack, 2, 1, 1, selection)) break;
1180 case 5: // kTrack pulls
1181 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1183 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1184 ((TVirtualPad*)l->At(0))->cd();
1185 if(!GetGraphArray(xy, kTrack, 1, 1)) break;
1187 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1188 ((TVirtualPad*)l->At(1))->cd();
1189 if(!GetGraphArray(xy, kTrack, 3, 1)) break;
1191 case 6: // kTrack phi
1192 xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
1193 if(GetGraph(&xy[0], kTrack , 4)) return kTRUE;
1195 case 7: // kTrackIn y
1196 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1197 xy[0] = -.3; xy[1] = -1500.; xy[2] = .3; xy[3] = 5000.;
1198 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1199 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1200 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1201 if(!GetGraphArray(xy, kTrackIn, 0, 1, n, selection)) break;
1202 pad=((TVirtualPad*)l->At(1)); pad->cd();
1203 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1204 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1205 if(!GetGraphArray(xy, kTrackIn, 0, 1, n, selection)) break;
1207 case 8: // kTrackIn y
1208 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1209 xy[0] = -.3; xy[1] = -1500.; xy[2] = .3; xy[3] = 5000.;
1210 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1211 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1212 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1213 if(!GetGraphArray(xy, kTrackIn, 0, 1, n, selection)) break;
1214 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1215 pad=((TVirtualPad*)l->At(1)); pad->cd();
1216 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1217 if(!GetGraphArray(xy, kTrackIn, 1, 1)) break;
1219 case 9: // kTrackIn z
1220 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1221 xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
1222 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1223 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1225 if(!GetGraphArray(xy, kTrackIn, 2, 1, 1, selection)) break;
1226 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1227 pad = ((TVirtualPad*)l->At(1)); pad->cd();
1228 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1229 if(!GetGraphArray(xy, kTrackIn, 3, 1)) break;
1231 case 10: // kTrackIn phi
1232 xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
1233 if(GetGraph(&xy[0], kTrackIn, 4)) return kTRUE;
1235 case 11: // kTrackOut y
1236 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1237 xy[0] = -.3; xy[1] = -50.; xy[2] = .3; xy[3] = 150.;
1238 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1239 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1240 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1241 if(!GetGraphArray(xy, kTrackOut, 0, 1, n, selection)) break;
1242 pad=((TVirtualPad*)l->At(1)); pad->cd();
1243 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1244 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1245 if(!GetGraphArray(xy, kTrackOut, 0, 1, n, selection)) break;
1247 case 12: // kTrackOut y
1248 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1249 xy[0] = -.3; xy[1] = -50.; xy[2] = .3; xy[3] = 150.;
1250 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1251 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1252 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1253 if(!GetGraphArray(xy, kTrackOut, 0, 1, n, selection)) break;
1254 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1255 pad=((TVirtualPad*)l->At(1)); pad->cd();
1256 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1257 if(!GetGraphArray(xy, kTrackOut, 1, 1)) break;
1259 case 13: // kTrackOut z
1260 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1261 xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
1262 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1263 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1264 if(!GetGraphArray(xy, kTrackOut, 2, 1)) break;
1265 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1266 pad = ((TVirtualPad*)l->At(1)); pad->cd();
1267 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1268 if(!GetGraphArray(xy, kTrackOut, 3, 1)) break;
1270 case 14: // kTrackOut phi
1271 xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
1272 if(GetGraph(&xy[0], kTrackOut, 4)) return kTRUE;
1274 case 15: // kMCcluster
1275 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1276 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
1277 ((TVirtualPad*)l->At(0))->cd();
1278 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1279 if(!GetGraphArray(xy, kMCcluster, 0, 1, n, selection)) break;
1280 ((TVirtualPad*)l->At(1))->cd();
1281 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1282 if(!GetGraphArray(xy, kMCcluster, 0, 1, n, selection)) break;
1284 case 16: // kMCcluster
1285 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1286 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3]=650.;
1287 ((TVirtualPad*)l->At(0))->cd();
1288 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1289 if(!GetGraphArray(xy, kMCcluster, 0, 1, n, selection)) break;
1290 ((TVirtualPad*)l->At(1))->cd();
1291 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1292 if(!GetGraphArray(xy, kMCcluster, 1, 1)) break;
1294 case 17: //kMCtracklet [y]
1295 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1296 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =500.;
1297 ((TVirtualPad*)l->At(0))->cd();
1298 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1299 if(!GetGraphArray(xy, kMCtracklet, 0, 1, n, selection)) break;
1300 ((TVirtualPad*)l->At(1))->cd();
1301 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1302 if(!GetGraphArray(xy, kMCtracklet, 0, 1, n, selection)) break;
1304 case 18: //kMCtracklet [y]
1305 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1306 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =500.;
1307 ((TVirtualPad*)l->At(0))->cd();
1308 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1309 if(!GetGraphArray(xy, kMCtracklet, 0, 1, n, selection)) break;
1310 ((TVirtualPad*)l->At(1))->cd();
1311 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1312 if(!GetGraphArray(xy, kMCtracklet, 1, 1)) break;
1314 case 19: //kMCtracklet [z]
1315 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1316 xy[0]=-1.; xy[1]=-100.; xy[2]=1.; xy[3] =2500.;
1317 ((TVirtualPad*)l->At(0))->cd();
1318 if(!GetGraphArray(xy, kMCtracklet, 2)) break;
1319 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1320 ((TVirtualPad*)l->At(1))->cd();
1321 if(!GetGraphArray(xy, kMCtracklet, 3)) break;
1323 case 20: //kMCtracklet [phi]
1324 xy[0]=-.3; xy[1]=-3.; xy[2]=.3; xy[3] =25.;
1325 if(!GetGraph(&xy[0], kMCtracklet, 4)) break;
1327 case 21: //kMCtrack [y] ly [0]
1328 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1329 xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1330 ((TVirtualPad*)l->At(0))->cd();
1331 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*0.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1332 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer1")) break;
1333 ((TVirtualPad*)l->At(1))->cd();
1334 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*0.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1335 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer1")) break;
1337 case 22: //kMCtrack [y] ly [1]
1338 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1339 xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1340 ((TVirtualPad*)l->At(0))->cd();
1341 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*1.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1342 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer2")) break;
1343 ((TVirtualPad*)l->At(1))->cd();
1344 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*1.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1345 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer2")) break;
1347 case 23: //kMCtrack [y] ly [2]
1348 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1349 xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1350 ((TVirtualPad*)l->At(0))->cd();
1351 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*2.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1352 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer3")) break;
1353 ((TVirtualPad*)l->At(1))->cd();
1354 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*2.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1355 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer3")) break;
1357 case 24: //kMCtrack [y] ly [3]
1358 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1359 xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1360 ((TVirtualPad*)l->At(0))->cd();
1361 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*3.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1362 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer4")) break;
1363 ((TVirtualPad*)l->At(1))->cd();
1364 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*3.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1365 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer4")) break;
1367 case 25: //kMCtrack [y] ly [4]
1368 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1369 xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1370 ((TVirtualPad*)l->At(0))->cd();
1371 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*4.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1372 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer5")) break;
1373 ((TVirtualPad*)l->At(1))->cd();
1374 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*4.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1375 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer5")) break;
1377 case 26: //kMCtrack [y] ly [5]
1378 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1379 xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1380 ((TVirtualPad*)l->At(0))->cd();
1381 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*5.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1382 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer6")) break;
1383 ((TVirtualPad*)l->At(1))->cd();
1384 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*5.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1385 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer6")) break;
1387 case 27: //kMCtrack [y pulls]
1388 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1389 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 5.5;
1390 ((TVirtualPad*)l->At(0))->cd();
1391 selStart=0; for(n=0; n<6; n++) selection[n]=selStart+n;
1392 if(!GetGraphArray(xy, kMCtrack, 1, 1, n, selection)) break;
1393 ((TVirtualPad*)l->At(1))->cd();
1394 selStart=6; for(n=0; n<6; n++) selection[n]=selStart+n;
1395 if(!GetGraphArray(xy, kMCtrack, 1, 1, n, selection)) break;
1397 case 28: //kMCtrack [z]
1398 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1399 xy[0]=-1.; xy[1]=-1500.; xy[2]=1.; xy[3] =6000.;
1400 ((TVirtualPad*)l->At(0))->cd();
1401 if(!GetGraphArray(xy, kMCtrack, 2)) break;
1402 xy[0] = -1.; xy[1] = -1.5; xy[2] = 1.; xy[3] = 5.;
1403 ((TVirtualPad*)l->At(1))->cd();
1404 if(!GetGraphArray(xy, kMCtrack, 3)) break;
1406 case 29: //kMCtrack [phi/snp]
1407 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1408 xy[0]=-.2; xy[1]=-0.5; xy[2]=.2; xy[3] =10.;
1409 ((TVirtualPad*)l->At(0))->cd();
1410 if(!GetGraphArray(xy, kMCtrack, 4)) break;
1411 xy[0] = -.2; xy[1] = -1.5; xy[2] = .2; xy[3] = 5.;
1412 ((TVirtualPad*)l->At(1))->cd();
1413 if(!GetGraphArray(xy, kMCtrack, 5)) break;
1415 case 30: //kMCtrack [theta/tgl]
1416 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1417 xy[0]=-1.; xy[1]=-0.5; xy[2]=1.; xy[3] =5.;
1418 ((TVirtualPad*)l->At(0))->cd();
1419 if(!GetGraphArray(xy, kMCtrack, 6)) break;
1420 xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
1421 ((TVirtualPad*)l->At(1))->cd();
1422 if(!GetGraphArray(xy, kMCtrack, 7)) break;
1424 case 31: //kMCtrack [pt]
1425 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1426 pad = (TVirtualPad*)l->At(0); pad->cd();
1427 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1430 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1431 selection[n++] = il*11 + 2; // pi-
1432 selection[n++] = il*11 + 8; // pi+
1434 xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 4.;
1435 //xy[0] = 0.2; xy[1] = -1.; xy[2] = 7.; xy[3] = 10.; // SA
1436 if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "#pi#pm")) break;
1437 pad->Modified(); pad->Update(); pad->SetLogx();
1438 pad = (TVirtualPad*)l->At(1); pad->cd();
1439 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1442 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1443 selection[n++] = il*11 + 3; // mu-
1444 selection[n++] = il*11 + 7; // mu+
1446 if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "#mu#pm")) break;
1447 pad->Modified(); pad->Update(); pad->SetLogx();
1449 case 32: //kMCtrack [pt]
1450 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1451 pad = (TVirtualPad*)l->At(0); pad->cd();
1452 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1455 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1456 selection[n++] = il*11 + 0; // p bar
1457 selection[n++] = il*11 + 10; // p
1459 xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 8.;
1460 //xy[0] = 0.2; xy[1] = -1.; xy[2] = 7.; xy[3] = 10.; // SA
1461 if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "p&p bar")) break;
1462 pad->Modified(); pad->Update(); pad->SetLogx();
1463 pad = (TVirtualPad*)l->At(1); pad->cd();
1464 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1467 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1468 selection[n++] = il*11 + 4; // e-
1469 selection[n++] = il*11 + 6; // e+
1471 xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 12.;
1472 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 14.; // SA
1473 if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "e#pm")) break;
1474 pad->Modified(); pad->Update(); pad->SetLogx();
1476 case 33: //kMCtrack [1/pt] pulls
1477 xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1478 //xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 4.5; // SA
1479 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1480 pad = (TVirtualPad*)l->At(0); pad->cd();
1481 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1484 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1485 selection[n++] = il*11 + 2; // pi-
1486 selection[n++] = il*11 + 8; // pi+
1488 if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "#pi#pm")) break;
1489 pad = (TVirtualPad*)l->At(1); pad->cd();
1490 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1493 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1494 selection[n++] = il*11 + 3; // mu-
1495 selection[n++] = il*11 + 7; // mu+
1497 if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "#mu#pm")) break;
1499 case 34: //kMCtrack [1/pt] pulls
1500 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1501 pad = (TVirtualPad*)l->At(0); pad->cd();
1502 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1505 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1506 selection[n++] = il*11 + 0; // p bar
1507 selection[n++] = il*11 + 10; // p
1509 xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1510 //xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 6.; // SA
1511 if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "p & p bar")) break;
1512 pad = (TVirtualPad*)l->At(1); pad->cd();
1513 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1516 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1517 selection[n++] = il*11 + 4; // e-
1518 selection[n++] = il*11 + 6; // e+
1520 xy[0] = 0.; xy[1] = -2.; xy[2] = 2.; xy[3] = 4.5;
1521 if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "e#pm")) break;
1523 case 35: //kMCtrack [p]
1524 xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 4.;
1525 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 10.;
1526 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1527 pad = (TVirtualPad*)l->At(0); pad->cd();
1528 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1531 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1532 selection[n++] = il*11 + 2; // pi-
1533 selection[n++] = il*11 + 8; // pi+
1535 if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "#pi#pm")) break;
1536 pad->Modified(); pad->Update(); pad->SetLogx();
1537 pad = (TVirtualPad*)l->At(1); pad->cd();
1538 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1541 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1542 selection[n++] = il*11 + 3; // mu-
1543 selection[n++] = il*11 + 7; // mu+
1545 if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "#mu#pm")) break;
1546 pad->Modified(); pad->Update(); pad->SetLogx();
1548 case 36: //kMCtrack [p]
1549 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1550 pad = (TVirtualPad*)l->At(0); pad->cd();
1551 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1554 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1555 selection[n++] = il*11 + 0; // p bar
1556 selection[n++] = il*11 + 10; // p
1558 xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 8.;
1559 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 12.; // SA
1560 if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "p & p bar")) break;
1561 pad->Modified(); pad->Update(); pad->SetLogx();
1562 pad = (TVirtualPad*)l->At(1); pad->cd();
1563 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1566 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1567 selection[n++] = il*11 + 4; // e-
1568 selection[n++] = il*11 + 6; // e+
1570 xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 12.;
1571 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 14.; // SA
1572 if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "e#pm")) break;
1573 pad->Modified(); pad->Update(); pad->SetLogx();
1575 case 37: // kMCtrackIn [y]
1576 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1577 xy[0]=-.25; xy[1]=-1000.; xy[2]=.25; xy[3] =3000.;
1578 ((TVirtualPad*)l->At(0))->cd();
1579 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1580 if(!GetGraphArray(xy, kMCtrackIn, 0, 1, n, selection)) break;
1581 ((TVirtualPad*)l->At(1))->cd();
1582 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1583 if(!GetGraphArray(&xy[0], kMCtrackIn, 0, 1, n, selection)) break;
1585 case 38: // kMCtrackIn [y]
1586 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1587 xy[0]=-.25; xy[1]=-1000.; xy[2]=.25; xy[3] =3000.;
1588 ((TVirtualPad*)l->At(0))->cd();
1589 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1590 if(!GetGraphArray(xy, kMCtrackIn, 0, 1, n, selection)) break;
1591 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1592 ((TVirtualPad*)l->At(1))->cd();
1593 if(!GetGraphArray(xy, kMCtrackIn, 1, 1)) break;
1595 case 39: // kMCtrackIn [z]
1596 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1597 xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =800.;
1598 ((TVirtualPad*)l->At(0))->cd();
1599 if(!GetGraphArray(xy, kMCtrackIn, 2, 1)) break;
1600 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1601 ((TVirtualPad*)l->At(1))->cd();
1602 if(!GetGraphArray(xy, kMCtrackIn, 3, 1)) break;
1604 case 40: // kMCtrackIn [phi|snp]
1605 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1606 xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
1607 ((TVirtualPad*)l->At(0))->cd();
1608 if(!GetGraph(&xy[0], kMCtrackIn, 4)) break;
1609 xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
1610 ((TVirtualPad*)l->At(1))->cd();
1611 if(!GetGraph(&xy[0], kMCtrackIn, 5)) break;
1613 case 41: // kMCtrackIn [theta|tgl]
1614 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1615 xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.;
1616 ((TVirtualPad*)l->At(0))->cd();
1617 if(!GetGraph(&xy[0], kMCtrackIn, 6)) break;
1618 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 1.5;
1619 ((TVirtualPad*)l->At(1))->cd();
1620 if(!GetGraph(&xy[0], kMCtrackIn, 7)) break;
1622 case 42: // kMCtrackIn [pt]
1623 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1624 xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1625 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 10.; // SA
1626 pad=(TVirtualPad*)l->At(0); pad->cd(); pad->SetLogx();
1627 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1628 n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1629 if(!GetGraphArray(xy, kMCtrackIn, 8, 1, n, selection)) break;
1630 pad = (TVirtualPad*)l->At(1); pad->cd(); pad->SetLogx();
1631 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1632 n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1633 if(!GetGraphArray(xy, kMCtrackIn, 8, 1, n, selection)) break;
1635 case 43: //kMCtrackIn [1/pt] pulls
1636 xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1637 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1638 pad = (TVirtualPad*)l->At(0); pad->cd();
1639 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1640 n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1641 if(!GetGraphArray(xy, kMCtrackIn, 9, 1, n, selection)) break;
1642 pad = (TVirtualPad*)l->At(1); pad->cd();
1643 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1644 n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1645 if(!GetGraphArray(xy, kMCtrackIn, 9, 1, n, selection)) break;
1647 case 44: // kMCtrackIn [p]
1648 xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1649 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 10.;
1650 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1651 pad = ((TVirtualPad*)l->At(0));pad->cd();pad->SetLogx();
1652 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1653 n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1654 if(!GetGraphArray(xy, kMCtrackIn, 10, 1, n, selection)) break;
1655 pad = ((TVirtualPad*)l->At(1)); pad->cd();pad->SetLogx();
1656 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1657 n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1658 if(!GetGraphArray(xy, kMCtrackIn, 10, 1, n, selection)) break;
1660 case 45: // kMCtrackOut [y]
1661 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1662 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =400.;
1663 ((TVirtualPad*)l->At(0))->cd();
1664 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1665 if(!GetGraphArray(xy, kMCtrackOut, 0, 1, n, selection)) break;
1666 ((TVirtualPad*)l->At(1))->cd();
1667 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1668 if(!GetGraphArray(&xy[0], kMCtrackOut, 0, 1, n, selection)) break;
1670 case 46: // kMCtrackOut [y]
1671 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1672 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =400.;
1673 ((TVirtualPad*)l->At(0))->cd();
1674 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1675 if(!GetGraphArray(xy, kMCtrackOut, 0, 1, n, selection)) break;
1676 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1677 ((TVirtualPad*)l->At(1))->cd();
1678 if(!GetGraphArray(xy, kMCtrackOut, 1, 1)) break;
1680 case 47: // kMCtrackOut [z]
1681 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1682 xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =1500.;
1683 ((TVirtualPad*)l->At(0))->cd();
1684 if(!GetGraphArray(xy, kMCtrackOut, 2, 1)) break;
1685 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1686 ((TVirtualPad*)l->At(1))->cd();
1687 if(!GetGraphArray(xy, kMCtrackOut, 3, 1)) break;
1689 case 48: // kMCtrackOut [phi|snp]
1690 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1691 xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
1692 ((TVirtualPad*)l->At(0))->cd();
1693 if(!GetGraph(&xy[0], kMCtrackOut, 4)) break;
1694 xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
1695 ((TVirtualPad*)l->At(1))->cd();
1696 if(!GetGraph(&xy[0], kMCtrackOut, 5)) break;
1698 case 49: // kMCtrackOut [theta|tgl]
1699 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1700 xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.;
1701 ((TVirtualPad*)l->At(0))->cd();
1702 if(!GetGraph(&xy[0], kMCtrackOut, 6)) break;
1703 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 15.;
1704 ((TVirtualPad*)l->At(1))->cd();
1705 if(!GetGraph(&xy[0], kMCtrackOut, 7)) break;
1707 case 50: // kMCtrackOut [pt]
1708 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1709 xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1710 pad=(TVirtualPad*)l->At(0); pad->cd(); pad->SetLogx();
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, kMCtrackOut, 8, 1, n, selection)) break;
1714 pad = (TVirtualPad*)l->At(1); pad->cd();pad->SetLogx();
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, kMCtrackOut, 8, 1, n, selection)) break;
1719 case 51: //kMCtrackOut [1/pt] pulls
1720 xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1721 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1722 pad = (TVirtualPad*)l->At(0); pad->cd();
1723 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1724 n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1725 if(!GetGraphArray(xy, kMCtrackOut, 9, 1, n, selection)) break;
1726 pad = (TVirtualPad*)l->At(1); pad->cd();
1727 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1728 n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1729 if(!GetGraphArray(xy, kMCtrackOut, 9, 1, n, selection)) break;
1731 case 52: // kMCtrackOut [p]
1732 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1733 xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1734 pad = ((TVirtualPad*)l->At(0));pad->cd();pad->SetLogx();
1735 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1736 n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1737 if(!GetGraphArray(xy, kMCtrackOut, 10, 1, n, selection)) break;
1738 pad = ((TVirtualPad*)l->At(1)); pad->cd();pad->SetLogx();
1739 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1740 n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1741 if(!GetGraphArray(xy, kMCtrackOut, 10, 1, n, selection)) break;
1744 AliWarning(Form("Reference plot [%d] missing result", ifig));
1748 //________________________________________________________
1749 void AliTRDresolution::MakeSummary()
1751 // Build summary plots
1753 if(!fGraphS || !fGraphM){
1754 AliError("Missing results");
1757 Float_t xy[4] = {0., 0., 0., 0.};
1759 TH2 *h2 = new TH2I("h2SF", "", 20, -.2, .2, fgkNresYsegm[fSegmentLevel], -0.5, fgkNresYsegm[fSegmentLevel]-0.5);
1760 h2->GetXaxis()->CenterTitle();
1761 h2->GetYaxis()->CenterTitle();
1762 h2->GetZaxis()->CenterTitle();h2->GetZaxis()->SetTitleOffset(1.4);
1764 Int_t ih2(0), iSumPlot(0);
1765 TCanvas *cOut = new TCanvas(Form("TRDsummary%s_%d", GetName(), iSumPlot++), "Cluster & Tracklet Resolution", 1024, 768);
1766 cOut->Divide(3,2, 2.e-3, 2.e-3, kYellow-7);
1767 TVirtualPad *p(NULL);
1770 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1771 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1772 h2->SetTitle(Form("Cluster-Track R-Phi Residuals;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1773 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphS->At(kCluster))->At(0), h2);
1774 GetRange(h2, 1, range);
1775 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1780 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1781 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1782 h2->SetTitle(Form("Cluster-Track R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1783 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphM->At(kCluster))->At(0), h2);
1784 GetRange(h2, 0, range);
1785 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1790 p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1791 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1792 GetGraphArray(xy, kCluster, 1, 1);
1795 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1796 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1797 h2->SetTitle(Form("Tracklet-Track R-Phi Residuals;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1798 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphS->At(kTrack))->At(0), h2);
1799 GetRange(h2, 1, range);
1800 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1805 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1806 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1807 h2->SetTitle(Form("Tracklet-Track R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1808 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphM->At(kTrack))->At(0), h2);
1809 GetRange(h2, 0, range);
1810 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1815 p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1816 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1817 GetGraphArray(xy, kTrack, 1, 1);
1819 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1825 cOut->Clear(); cOut->SetName(Form("TRDsummary%s_%d", GetName(), iSumPlot++));
1826 cOut->Divide(3, 2, 2.e-3, 2.e-3, kBlue-10);
1829 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1830 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1831 h2->SetTitle(Form("Cluster-MC R-Phi Resolution;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1832 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphS->At(kMCcluster))->At(0), h2);
1833 GetRange(h2, 1, range);
1834 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1839 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1840 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1842 h2->SetTitle(Form("Cluster-MC R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1843 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphM->At(kMCcluster))->At(0), h2);
1844 GetRange(h2, 0, range);
1845 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1850 p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1851 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1852 GetGraphArray(xy, kMCcluster, 1, 1);
1855 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1856 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1858 h2->SetTitle(Form("Tracklet-MC R-Phi Resolution;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1859 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphS->At(kMCtracklet))->At(0), h2);
1860 GetRange(h2, 1, range);
1861 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1866 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1867 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1869 h2->SetTitle(Form("Tracklet-MC R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1870 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphM->At(kMCtracklet))->At(0), h2);
1871 GetRange(h2, 0, range);
1872 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1877 p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1878 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1879 GetGraphArray(xy, kMCtracklet, 1, 1);
1881 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1885 //________________________________________________________
1886 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1888 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1889 // The "range" vector should be initialized with 2 elements
1890 // Option "mod" can be any of
1891 // - 0 : gaussian like distribution
1892 // - 1 : tailed distribution
1894 Int_t nx(h2->GetNbinsX())
1895 , ny(h2->GetNbinsY())
1897 Double_t *data=new Double_t[n];
1898 for(Int_t ix(1), in(0); ix<=nx; ix++){
1899 for(Int_t iy(1); iy<=ny; iy++)
1900 data[in++] = h2->GetBinContent(ix, iy);
1902 Double_t mean, sigm;
1903 AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1905 range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1906 if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1907 AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1908 TH1S h1("h1SF0", "", 100, range[0], range[1]);
1913 case 0:// gaussian distribution
1915 TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1917 mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1918 range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1919 AliDebug(2, Form(" rangeG[%f %f]", range[0], range[1]));
1922 case 1:// tailed distribution
1924 Int_t bmax(h1.GetMaximumBin());
1925 Int_t jBinMin(1), jBinMax(100);
1926 for(Int_t ibin(bmax); ibin--;){
1927 if(h1.GetBinContent(ibin)<1.){
1928 jBinMin=ibin; break;
1931 for(Int_t ibin(bmax); ibin++;){
1932 if(h1.GetBinContent(ibin)<1.){
1933 jBinMax=ibin; break;
1936 range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1937 AliDebug(2, Form(" rangeT[%f %f]", range[0], range[1]));
1945 //________________________________________________________
1946 void AliTRDresolution::MakeSummaryPlot(TObjArray *a, TH2 *h2)
1948 // Core functionality for MakeSummary function.
1952 TGraphErrors *g(NULL); TAxis *ax(h2->GetXaxis());
1953 for(Int_t iseg(0); iseg<fgkNresYsegm[fSegmentLevel]; iseg++){
1954 g=(TGraphErrors*)a->At(iseg);
1955 for(Int_t in(0); in<g->GetN(); in++){
1956 g->GetPoint(in, x, y);
1957 h2->SetBinContent(ax->FindBin(x), iseg+1, y);
1963 //________________________________________________________
1964 Bool_t AliTRDresolution::PostProcess()
1966 // Fit, Project, Combine, Extract values from the containers filled during execution
1968 /*fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));*/
1970 AliError("ERROR: list not available");
1974 // define general behavior parameters
1975 const Color_t fgColorS[11]={
1976 kOrange, kOrange-3, kMagenta+1, kViolet, kRed,
1978 kRed, kViolet, kMagenta+1, kOrange-3, kOrange
1980 const Color_t fgColorM[11]={
1981 kCyan-5, kAzure-4, kBlue-7, kBlue+2, kViolet+10,
1983 kViolet+10, kBlue+2, kBlue-7, kAzure-4, kCyan-5
1985 const Marker_t fgMarker[11]={
1991 TGraph *gm= NULL, *gs= NULL;
1992 if(!fGraphS && !fGraphM){
1993 TObjArray *aM(NULL), *aS(NULL);
1994 Int_t n = fContainer->GetEntriesFast();
1995 fGraphS = new TObjArray(n); fGraphS->SetOwner();
1996 fGraphM = new TObjArray(n); fGraphM->SetOwner();
1997 for(Int_t ig(0), nc(0); ig<n; ig++){
1998 fGraphM->AddAt(aM = new TObjArray(fgNproj[ig]), ig);
1999 fGraphS->AddAt(aS = new TObjArray(fgNproj[ig]), ig);
2001 for(Int_t ic=0; ic<fgNproj[ig]; ic++, nc++){
2002 AliDebug(2, Form("building G[%d] P[%d] N[%d]", ig, ic, fNcomp[nc]));
2004 TObjArray *agS(NULL), *agM(NULL);
2005 aS->AddAt(agS = new TObjArray(fNcomp[nc]), ic);
2006 aM->AddAt(agM = new TObjArray(fNcomp[nc]), ic);
2007 for(Int_t is=fNcomp[nc]; is--;){
2008 agS->AddAt(gs = new TGraphErrors(), is);
2009 Int_t is0(is%11), il0(is/11);
2010 gs->SetMarkerStyle(fgMarker[is0]);
2011 gs->SetMarkerColor(fgColorS[is0]);
2012 gs->SetLineColor(fgColorS[is0]);
2013 gs->SetLineStyle(il0);gs->SetLineWidth(2);
2014 gs->SetName(Form("s_%d_%02d_%02d", ig, ic, is));
2016 agM->AddAt(gm = new TGraphErrors(), is);
2017 gm->SetMarkerStyle(fgMarker[is0]);
2018 gm->SetMarkerColor(fgColorM[is0]);
2019 gm->SetLineColor(fgColorM[is0]);
2020 gm->SetLineStyle(il0);gm->SetLineWidth(2);
2021 gm->SetName(Form("m_%d_%02d_%02d", ig, ic, is));
2022 // this is important for labels in the legend
2024 gs->SetTitle(Form("%s %02d", fgkResYsegmName[fSegmentLevel], is%fgkNresYsegm[fSegmentLevel]));
2025 gm->SetTitle(Form("%s %02d", fgkResYsegmName[fSegmentLevel], is%fgkNresYsegm[fSegmentLevel]));
2027 gs->SetTitle(Form("%s Ly[%d]", is%2 ?"z":"y", is/2));
2028 gm->SetTitle(Form("%s Ly[%d]", is%2?"z":"y", is/2));
2029 } else if(ic==2||ic==3) {
2030 gs->SetTitle(Form("%s Ly[%d]", is%2 ?"RC":"no RC", is/2));
2031 gm->SetTitle(Form("%s Ly[%d]", is%2?"RC":"no RC", is/2));
2033 gs->SetTitle(Form("Layer[%d]", is%AliTRDgeometry::kNlayer));
2034 gm->SetTitle(Form("Layer[%d]", is%AliTRDgeometry::kNlayer));
2036 gs->SetTitle(Form("%s @ ly[%d]", fgParticle[is0], il0));
2037 gm->SetTitle(Form("%s @ ly[%d]", fgParticle[is0], il0));
2041 aS->AddAt(gs = new TGraphErrors(), ic);
2042 gs->SetMarkerStyle(23);
2043 gs->SetMarkerColor(kRed);
2044 gs->SetLineColor(kRed);
2045 gs->SetNameTitle(Form("s_%d_%02d", ig, ic), "sigma");
2047 aM->AddAt(gm = ig ? (TGraph*)new TGraphErrors() : (TGraph*)new TGraphAsymmErrors(), ic);
2048 gm->SetLineColor(kBlack);
2049 gm->SetMarkerStyle(7);
2050 gm->SetMarkerColor(kBlack);
2051 gm->SetNameTitle(Form("m_%d_%02d", ig, ic), "mean");
2057 /* printf("\n\n\n"); fGraphS->ls();
2058 printf("\n\n\n"); fGraphM->ls();*/
2063 TF1 fg("fGauss", "gaus", -.5, .5);
2064 // Landau for charge resolution
2065 TF1 fch("fClCh", "landau", 0., 1000.);
2066 // Landau for e+- pt resolution
2067 TF1 fpt("fPt", "landau", -0.1, 0.2);
2069 //PROCESS EXPERIMENTAL DISTRIBUTIONS
2070 // Charge resolution
2071 //Process3DL(kCharge, 0, &fl);
2072 // Clusters residuals
2073 Process3D(kCluster, 0, &fg, 1.e4);
2074 Process3Dlinked(kCluster, 1, &fg);
2076 // Tracklet residual/pulls
2077 Process3D(kTrack , 0, &fg, 1.e4);
2078 Process3Dlinked(kTrack , 1, &fg);
2079 Process3D(kTrack , 2, &fg, 1.e4);
2080 Process3D(kTrack , 3, &fg);
2081 Process2D(kTrack , 4, &fg, 1.e3);
2083 // TRDin residual/pulls
2084 Process3D(kTrackIn, 0, &fg, 1.e4);
2085 Process3Dlinked(kTrackIn, 1, &fg);
2086 Process3D(kTrackIn, 2, &fg, 1.e4);
2087 Process3D(kTrackIn, 3, &fg);
2088 Process2D(kTrackIn, 4, &fg, 1.e3);
2090 // TRDout residual/pulls
2091 Process3D(kTrackOut, 0, &fg, 1.e3); // scale to fit - see PlotTrackOut
2092 Process3Dlinked(kTrackOut, 1, &fg);
2093 Process3D(kTrackOut, 2, &fg, 1.e4);
2094 Process3D(kTrackOut, 3, &fg);
2095 Process2D(kTrackOut, 4, &fg, 1.e3);
2098 if(!HasMCdata()) return kTRUE;
2101 //PROCESS MC RESIDUAL DISTRIBUTIONS
2103 // CLUSTER Y RESOLUTION/PULLS
2104 Process3D(kMCcluster, 0, &fg, 1.e4);
2105 Process3Dlinked(kMCcluster, 1, &fg, 1.);
2108 // TRACKLET RESOLUTION/PULLS
2109 Process3D(kMCtracklet, 0, &fg, 1.e4); // y
2110 Process3Dlinked(kMCtracklet, 1, &fg, 1.); // y pulls
2111 Process3D(kMCtracklet, 2, &fg, 1.e4); // z
2112 Process3D(kMCtracklet, 3, &fg, 1.); // z pulls
2113 Process2D(kMCtracklet, 4, &fg, 1.e3); // phi
2116 // TRACK RESOLUTION/PULLS
2117 Process3Darray(kMCtrack, 0, &fg, 1.e4); // y
2118 Process3DlinkedArray(kMCtrack, 1, &fg); // y PULL
2119 Process3Darray(kMCtrack, 2, &fg, 1.e4); // z
2120 Process3Darray(kMCtrack, 3, &fg); // z PULL
2121 Process2Darray(kMCtrack, 4, &fg, 1.e3); // phi
2122 Process2Darray(kMCtrack, 5, &fg); // snp PULL
2123 Process2Darray(kMCtrack, 6, &fg, 1.e3); // theta
2124 Process2Darray(kMCtrack, 7, &fg); // tgl PULL
2125 Process3Darray(kMCtrack, 8, &fg, 1.e2); // pt resolution
2126 Process3Darray(kMCtrack, 9, &fg); // 1/pt pulls
2127 Process3Darray(kMCtrack, 10, &fg, 1.e2); // p resolution
2130 // TRACK TRDin RESOLUTION/PULLS
2131 Process3D(kMCtrackIn, 0, &fg, 1.e4);// y resolution
2132 Process3Dlinked(kMCtrackIn, 1, &fg); // y pulls
2133 Process3D(kMCtrackIn, 2, &fg, 1.e4);// z resolution
2134 Process3D(kMCtrackIn, 3, &fg); // z pulls
2135 Process2D(kMCtrackIn, 4, &fg, 1.e3);// phi resolution
2136 Process2D(kMCtrackIn, 5, &fg); // snp pulls
2137 Process2D(kMCtrackIn, 6, &fg, 1.e3);// theta resolution
2138 Process2D(kMCtrackIn, 7, &fg); // tgl pulls
2139 Process3D(kMCtrackIn, 8, &fg, 1.e2);// pt resolution
2140 Process3D(kMCtrackIn, 9, &fg); // 1/pt pulls
2141 Process3D(kMCtrackIn, 10, &fg, 1.e2);// p resolution
2144 // TRACK TRDout RESOLUTION/PULLS
2145 Process3D(kMCtrackOut, 0, &fg, 1.e4);// y resolution
2146 Process3Dlinked(kMCtrackOut, 1, &fg); // y pulls
2147 Process3D(kMCtrackOut, 2, &fg, 1.e4);// z resolution
2148 Process3D(kMCtrackOut, 3, &fg); // z pulls
2149 Process2D(kMCtrackOut, 4, &fg, 1.e3);// phi resolution
2150 Process2D(kMCtrackOut, 5, &fg); // snp pulls
2151 Process2D(kMCtrackOut, 6, &fg, 1.e3);// theta resolution
2152 Process2D(kMCtrackOut, 7, &fg); // tgl pulls
2153 Process3D(kMCtrackOut, 8, &fg, 1.e2);// pt resolution
2154 Process3D(kMCtrackOut, 9, &fg); // 1/pt pulls
2155 Process3D(kMCtrackOut, 10, &fg, 1.e2);// p resolution
2162 //________________________________________________________
2163 void AliTRDresolution::Terminate(Option_t *opt)
2165 AliTRDrecoTask::Terminate(opt);
2166 if(HasPostProcess()) PostProcess();
2169 //________________________________________________________
2170 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
2172 // Helper function to avoid duplication of code
2173 // Make first guesses on the fit parameters
2175 // find the intial parameters of the fit !! (thanks George)
2176 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
2178 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
2179 f->SetParLimits(0, 0., 3.*sum);
2180 f->SetParameter(0, .9*sum);
2181 Double_t rms = h->GetRMS();
2182 f->SetParLimits(1, -rms, rms);
2183 f->SetParameter(1, h->GetMean());
2185 f->SetParLimits(2, 0., 2.*rms);
2186 f->SetParameter(2, rms);
2187 if(f->GetNpar() <= 4) return;
2189 f->SetParLimits(3, 0., sum);
2190 f->SetParameter(3, .1*sum);
2192 f->SetParLimits(4, -.3, .3);
2193 f->SetParameter(4, 0.);
2195 f->SetParLimits(5, 0., 1.e2);
2196 f->SetParameter(5, 2.e-1);
2199 //________________________________________________________
2200 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand)
2202 // Build performance histograms for AliTRDcluster.vs TRD track or MC
2203 // - y reziduals/pulls
2205 TObjArray *arr = new TObjArray(2);
2206 arr->SetName(name); arr->SetOwner();
2207 TH1 *h(NULL); char hname[100], htitle[300];
2209 // tracklet resolution/pull in y direction
2210 sprintf(hname, "%s_%s_Y", GetNameId(), name);
2211 sprintf(htitle, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
2212 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2213 Int_t nybins=fgkNresYsegm[fSegmentLevel];
2214 if(expand) nybins*=2;
2215 h = new TH3S(hname, htitle,
2216 48, -.48, .48, 60, -1.5, 1.5, nybins, -0.5, nybins-0.5);
2219 sprintf(hname, "%s_%s_YZpull", GetNameId(), name);
2220 sprintf(htitle, "YZ pull for \"%s\" @ %s;%s;#Delta y / #sigma_{y};#Delta z / #sigma_{z}", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
2221 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2222 h = new TH3S(hname, htitle, fgkNresYsegm[fSegmentLevel], -0.5, fgkNresYsegm[fSegmentLevel]-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
2229 //________________________________________________________
2230 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
2232 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
2233 // - y reziduals/pulls
2234 // - z reziduals/pulls
2236 TObjArray *arr = BuildMonitorContainerCluster(name, expand);
2238 TH1 *h(NULL); char hname[100], htitle[300];
2240 // tracklet resolution/pull in z direction
2241 sprintf(hname, "%s_%s_Z", GetNameId(), name);
2242 sprintf(htitle, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm];row cross", GetNameId(), name);
2243 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2244 h = new TH3S(hname, htitle, 50, -1., 1., 100, -1.5, 1.5, 2, -0.5, 1.5);
2247 sprintf(hname, "%s_%s_Zpull", GetNameId(), name);
2248 sprintf(htitle, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z / #sigma_{z};row cross", GetNameId(), name);
2249 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2250 h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
2251 h->GetZaxis()->SetBinLabel(1, "no RC");
2252 h->GetZaxis()->SetBinLabel(2, "RC");
2256 // tracklet to track phi resolution
2257 sprintf(hname, "%s_%s_PHI", GetNameId(), name);
2258 sprintf(htitle, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];entries", GetNameId(), name);
2259 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2260 h = new TH2I(hname, htitle, 21, -.33, .33, 100, -.5, .5);
2267 //________________________________________________________
2268 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
2270 // Build performance histograms for AliExternalTrackParam.vs MC
2271 // - y resolution/pulls
2272 // - z resolution/pulls
2273 // - phi resolution, snp pulls
2274 // - theta resolution, tgl pulls
2275 // - pt resolution, 1/pt pulls, p resolution
2277 TObjArray *arr = BuildMonitorContainerTracklet(name);
2279 TH1 *h(NULL); char hname[100], htitle[300];
2283 sprintf(hname, "%s_%s_SNPpull", GetNameId(), name);
2284 sprintf(htitle, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp / #sigma_{snp};entries", GetNameId(), name);
2285 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2286 h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
2291 sprintf(hname, "%s_%s_THT", GetNameId(), name);
2292 sprintf(htitle, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
2293 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2294 h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
2298 sprintf(hname, "%s_%s_TGLpull", GetNameId(), name);
2299 sprintf(htitle, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl / #sigma_{tgl};entries", GetNameId(), name);
2300 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2301 h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
2305 const Int_t kNpt(14);
2306 const Int_t kNdpt(150);
2307 const Int_t kNspc = 2*AliPID::kSPECIES+1;
2308 Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
2309 Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
2310 for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
2311 for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
2312 for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
2315 sprintf(hname, "%s_%s_Pt", GetNameId(), name);
2316 sprintf(htitle, "P_{t} res for \"%s\" @ %s;p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
2317 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2318 h = new TH3S(hname, htitle,
2319 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2321 for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2325 sprintf(hname, "%s_%s_1Pt", GetNameId(), name);
2326 sprintf(htitle, "1/P_{t} pull for \"%s\" @ %s;1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", GetNameId(), name);
2327 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2328 h = new TH3S(hname, htitle,
2329 kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
2331 for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2335 sprintf(hname, "%s_%s_P", GetNameId(), name);
2336 sprintf(htitle, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
2337 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2338 h = new TH3S(hname, htitle,
2339 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2341 for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2349 //________________________________________________________
2350 TObjArray* AliTRDresolution::Histos()
2353 // Define histograms
2356 if(fContainer) return fContainer;
2358 fContainer = new TObjArray(kNviews);
2359 //fContainer->SetOwner(kTRUE);
2361 TObjArray *arr(NULL);
2363 // binnings for plots containing momentum or pt
2364 const Int_t kNpt(14), kNphi(48), kNdy(60);
2365 Float_t lPhi=-.48, lDy=-.3, lPt=0.1;
2366 Float_t binsPhi[kNphi+1], binsDy[kNdy+1], binsPt[kNpt+1];
2367 for(Int_t i=0; i<kNphi+1; i++,lPhi+=.02) binsPhi[i]=lPhi;
2368 for(Int_t i=0; i<kNdy+1; i++,lDy+=.01) binsDy[i]=lDy;
2369 for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
2371 // cluster to track residuals/pulls
2372 fContainer->AddAt(arr = new TObjArray(2), kCharge);
2373 arr->SetName("Charge");
2374 if(!(h = (TH3S*)gROOT->FindObject("hCharge"))){
2375 h = new TH3S("hCharge", "Charge Resolution", 20, 1., 2., 24, 0., 3.6, 100, 0., 500.);
2376 h->GetXaxis()->SetTitle("dx/dx_{0}");
2377 h->GetYaxis()->SetTitle("x_{d} [cm]");
2378 h->GetZaxis()->SetTitle("dq/dx [ADC/cm]");
2382 // cluster to track residuals/pulls
2383 fContainer->AddAt(BuildMonitorContainerCluster("Cl"), kCluster);
2384 // tracklet to TRD track
2385 fContainer->AddAt(BuildMonitorContainerTracklet("Trk", kTRUE), kTrack);
2386 // tracklet to TRDin
2387 fContainer->AddAt(BuildMonitorContainerTracklet("TrkIN", kTRUE), kTrackIn);
2388 // tracklet to TRDout
2389 fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
2392 // Resolution histos
2393 if(!HasMCdata()) return fContainer;
2395 // cluster resolution
2396 fContainer->AddAt(BuildMonitorContainerCluster("MCcl"), kMCcluster);
2398 // tracklet resolution
2399 fContainer->AddAt(BuildMonitorContainerTracklet("MCtracklet"), kMCtracklet);
2402 fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
2403 arr->SetName("MCtrk");
2404 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
2406 // TRDin TRACK RESOLUTION
2407 fContainer->AddAt(BuildMonitorContainerTrack("MCtrkIN"), kMCtrackIn);
2409 // TRDout TRACK RESOLUTION
2410 fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
2415 //________________________________________________________
2416 Bool_t AliTRDresolution::Load(const Char_t *file, const Char_t *dir)
2418 // Custom load function. Used to guess the segmentation level of the data.
2420 if(!AliTRDrecoTask::Load(file, dir)) return kFALSE;
2422 // look for cluster residual plot - always available
2423 TH3S* h3((TH3S*)((TObjArray*)fContainer->At(kClToTrk))->At(0));
2424 Int_t segmentation(h3->GetNbinsZ()/2);
2425 if(segmentation==fgkNresYsegm[0]){ // default segmentation. Nothing to do
2427 } else if(segmentation==fgkNresYsegm[1]){ // stack segmentation.
2428 SetSegmentationLevel(1);
2429 } else if(segmentation==fgkNresYsegm[2]){ // detector segmentation.
2430 SetSegmentationLevel(2);
2432 AliError(Form("Unknown segmentation [%d].", h3->GetNbinsZ()));
2436 AliDebug(2, Form("Segmentation set to level \"%s\"", fgkResYsegmName[fSegmentLevel]));
2441 //________________________________________________________
2442 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
2445 // Do the processing
2448 Char_t pn[10]; sprintf(pn, "p%03d", fIdxPlot);
2450 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
2451 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
2452 if(Int_t(h2->GetEntries())){
2453 AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
2455 AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
2460 const Int_t kINTEGRAL=1;
2461 for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
2462 Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
2463 Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
2464 TH1D *h = h2->ProjectionY(pn, abin, bbin);
2465 if((n=(Int_t)h->GetEntries())<150){
2466 AliDebug(4, Form(" x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
2470 Int_t ip = g[0]->GetN();
2471 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)));
2472 g[0]->SetPoint(ip, x, k*f->GetParameter(1));
2473 g[0]->SetPointError(ip, 0., k*f->GetParError(1));
2474 g[1]->SetPoint(ip, x, k*f->GetParameter(2));
2475 g[1]->SetPointError(ip, 0., k*f->GetParError(2));
2477 g[0]->SetPoint(ip, x, k*h->GetMean());
2478 g[0]->SetPointError(ip, 0., k*h->GetMeanError());
2479 g[1]->SetPoint(ip, x, k*h->GetRMS());
2480 g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
2486 //________________________________________________________
2487 Bool_t AliTRDresolution::Process2D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k, Int_t gidx)
2490 // Do the processing
2493 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2495 // retrive containers
2498 if(!(h2= (TH2I*)(fContainer->At(plot)))) return kFALSE;
2500 TObjArray *a0(NULL);
2501 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2502 if(!(h2=(TH2I*)a0->At(idx))) return kFALSE;
2504 if(Int_t(h2->GetEntries())){
2505 AliDebug(2, Form("p[%d] idx[%d] : h[%s] %s", plot, idx, h2->GetName(), h2->GetTitle()));
2507 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2512 if(gidx<0) gidx=idx;
2513 if(!(g[0] = gidx<0 ? (TGraphErrors*)fGraphM->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(gidx))) return kFALSE;
2515 if(!(g[1] = gidx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(gidx))) return kFALSE;
2517 return Process(h2, f, k, g);
2520 //________________________________________________________
2521 Bool_t AliTRDresolution::Process3D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2524 // Do the processing
2527 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2529 // retrive containers
2532 if(!(h3= (TH3S*)(fContainer->At(plot)))) return kFALSE;
2534 TObjArray *a0(NULL);
2535 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2536 if(!(h3=(TH3S*)a0->At(idx))) return kFALSE;
2538 if(Int_t(h3->GetEntries())){
2539 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2541 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2546 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2547 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2550 TAxis *az = h3->GetZaxis();
2551 for(Int_t iz(0); iz<gm->GetEntriesFast(); iz++){
2552 if(!(g[0] = (TGraphErrors*)gm->At(iz))) return kFALSE;
2553 if(!(g[1] = (TGraphErrors*)gs->At(iz))) return kFALSE;
2554 az->SetRange(iz+1, iz+1);
2555 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2562 //________________________________________________________
2563 Bool_t AliTRDresolution::Process3Dlinked(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2566 // Do the processing
2569 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2571 // retrive containers
2574 if(!(h3= (TH3S*)(fContainer->At(plot)))) return kFALSE;
2576 TObjArray *a0(NULL);
2577 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2578 if(!(h3=(TH3S*)a0->At(idx))) return kFALSE;
2580 if(Int_t(h3->GetEntries())){
2581 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2583 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2588 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2589 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2592 if(!(g[0] = (TGraphErrors*)gm->At(0))) return kFALSE;
2593 if(!(g[1] = (TGraphErrors*)gs->At(0))) return kFALSE;
2594 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2596 if(!(g[0] = (TGraphErrors*)gm->At(1))) return kFALSE;
2597 if(!(g[1] = (TGraphErrors*)gs->At(1))) return kFALSE;
2598 if(!Process((TH2*)h3->Project3D("zx"), f, k, g)) return kFALSE;
2604 //________________________________________________________
2605 Bool_t AliTRDresolution::Process3DL(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2608 // Do the processing
2611 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2613 // retrive containers
2614 TH3S *h3 = (TH3S*)((TObjArray*)fContainer->At(plot))->At(idx);
2615 if(!h3) return kFALSE;
2616 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2618 TGraphAsymmErrors *gm;
2620 if(!(gm = (TGraphAsymmErrors*)((TObjArray*)fGraphM->At(plot))->At(0))) return kFALSE;
2621 if(!(gs = (TGraphErrors*)((TObjArray*)fGraphS->At(plot)))) return kFALSE;
2623 Float_t x, r, mpv, xM, xm;
2624 TAxis *ay = h3->GetYaxis();
2625 for(Int_t iy=1; iy<=h3->GetNbinsY(); iy++){
2626 ay->SetRange(iy, iy);
2627 x = ay->GetBinCenter(iy);
2628 TH2F *h2=(TH2F*)h3->Project3D("zx");
2629 TAxis *ax = h2->GetXaxis();
2630 for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
2631 r = ax->GetBinCenter(ix);
2632 TH1D *h1 = h2->ProjectionY("py", ix, ix);
2633 if(h1->Integral()<50) continue;
2636 GetLandauMpvFwhm(f, mpv, xm, xM);
2637 Int_t ip = gm->GetN();
2638 gm->SetPoint(ip, x, k*mpv);
2639 gm->SetPointError(ip, 0., 0., k*xm, k*xM);
2640 gs->SetPoint(ip, r, k*(xM-xm)/mpv);
2641 gs->SetPointError(ip, 0., 0.);
2648 //________________________________________________________
2649 Bool_t AliTRDresolution::Process2Darray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2652 // Do the processing
2655 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2657 // retrive containers
2658 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2659 if(!arr) return kFALSE;
2660 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2663 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2664 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2666 TGraphErrors *g[2]; TH2I *h2(NULL); TObjArray *a0(NULL);
2667 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2668 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2670 if(!(h2 = (TH2I*)a0->At(idx))) return kFALSE;
2671 if(Int_t(h2->GetEntries())){
2672 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h2->GetName(), h2->GetTitle()));
2674 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2678 if(!(g[0] = (TGraphErrors*)gm->At(ia))) return kFALSE;
2679 if(!(g[1] = (TGraphErrors*)gs->At(ia))) return kFALSE;
2680 if(!Process(h2, f, k, g)) return kFALSE;
2686 //________________________________________________________
2687 Bool_t AliTRDresolution::Process3Darray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2690 // Do the processing
2693 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2694 //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
2696 // retrive containers
2697 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2698 if(!arr) return kFALSE;
2699 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2702 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2703 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2705 TGraphErrors *g[2]; TH3S *h3(NULL); TObjArray *a0(NULL);
2707 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2708 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2710 if(!(h3 = (TH3S*)a0->At(idx))) return kFALSE;
2711 if(Int_t(h3->GetEntries())){
2712 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h3->GetName(), h3->GetTitle()));
2714 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2717 TAxis *az = h3->GetZaxis();
2718 for(Int_t iz=1; iz<=az->GetNbins(); iz++, in++){
2719 if(in >= gm->GetEntriesFast()) break;
2720 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2721 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2722 az->SetRange(iz, iz);
2723 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2726 AliDebug(2, Form("Projections [%d] from [%d]", in, gs->GetEntriesFast()));
2731 //________________________________________________________
2732 Bool_t AliTRDresolution::Process3DlinkedArray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2735 // Do the processing
2738 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2739 //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
2741 // retrive containers
2742 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2743 if(!arr) return kFALSE;
2744 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2747 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2748 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2750 TGraphErrors *g[2]; TH3S *h3(NULL); TObjArray *a0(NULL);
2752 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2753 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2754 if(!(h3 = (TH3S*)a0->At(idx))) return kFALSE;
2755 if(Int_t(h3->GetEntries())){
2756 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h3->GetName(), h3->GetTitle()));
2758 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2761 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2762 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2763 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2766 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2767 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2768 if(!Process((TH2*)h3->Project3D("zx"), f, k, g)) return kFALSE;
2771 AliDebug(2, Form("Projections [%d] from [%d]", in, gs->GetEntriesFast()));
2776 //________________________________________________________
2777 Bool_t AliTRDresolution::GetGraph(Float_t *bb, ETRDresolutionPlot ip, Int_t idx, Bool_t kLEG, const Char_t *explain)
2783 if(!fGraphS || !fGraphM) return kFALSE;
2784 // axis titles look up
2786 for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fgNproj[jp];
2787 UChar_t jdx = idx<0?0:idx;
2788 for(Int_t jc=0; jc<TMath::Min(jdx,fgNproj[ip]-1); jc++) nref++;
2789 Char_t **at = fAxTitle[nref];
2791 // build legends if requiered
2794 leg=new TLegend(.65, .6, .95, .9);
2795 leg->SetBorderSize(0);
2796 leg->SetFillStyle(0);
2800 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]);
2801 h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2802 h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
2804 TAxis *ax = h1->GetXaxis();
2805 ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
2806 ax = h1->GetYaxis();
2807 ax->SetRangeUser(bb[1], bb[3]);
2808 ax->CenterTitle(); ax->SetTitleOffset(1.4);
2811 TBox *b = new TBox(-.15, bb[1], .15, bb[3]);
2812 b->SetFillStyle(3002);b->SetLineColor(0);
2813 b->SetFillColor(ip<=Int_t(kMCcluster)?kGreen:kBlue);
2816 TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphM->At(ip)))->At(idx);
2817 if(!gm) return kFALSE;
2818 TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphS->At(ip)))->At(idx);
2819 if(!gs) return kFALSE;
2821 Int_t n(0), nPlots(0);
2822 if((n=gm->GetN())) {
2824 gm->Draw("pl"); if(leg) leg->AddEntry(gm, gm->GetTitle(), "pl");
2825 PutTrendValue(Form("%s_%s", fgPerformanceName[ip], at[0]), gm->GetMean(2));
2826 PutTrendValue(Form("%s_%sRMS", fgPerformanceName[ip], at[0]), gm->GetRMS(2));
2831 gs->Draw("pl"); if(leg) leg->AddEntry(gs, gs->GetTitle(), "pl");
2832 gs->Sort(&TGraph::CompareY);
2833 PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[ip], at[0]), gs->GetY()[0]);
2834 PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[ip], at[0]), gs->GetY()[n-1]);
2835 gs->Sort(&TGraph::CompareX);
2837 if(!nPlots) return kFALSE;
2838 if(leg) leg->Draw();
2842 //________________________________________________________
2843 Bool_t AliTRDresolution::GetGraphArray(Float_t *bb, ETRDresolutionPlot ip, Int_t idx, Bool_t kLEG, Int_t n, Int_t *sel, const Char_t *explain)
2849 if(!fGraphS || !fGraphM) return kFALSE;
2851 // axis titles look up
2853 for(Int_t jp(0); jp<ip; jp++) nref+=fgNproj[jp];
2855 Char_t **at = fAxTitle[nref];
2857 // build legends if requiered
2858 TLegend *legM(NULL), *legS(NULL);
2860 legM=new TLegend(.35, .6, .65, .9);
2861 legM->SetHeader("Mean");
2862 legM->SetBorderSize(0);
2863 legM->SetFillStyle(0);
2864 legS=new TLegend(.65, .6, .95, .9);
2865 legS->SetHeader("Sigma");
2866 legS->SetBorderSize(0);
2867 legS->SetFillStyle(0);
2871 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]);
2872 h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2873 h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
2875 TAxis *ax = h1->GetXaxis();
2876 ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
2877 ax = h1->GetYaxis();
2878 ax->SetRangeUser(bb[1], bb[3]);
2879 ax->CenterTitle(); ax->SetTitleOffset(1.4);
2882 TGraphErrors *gm(NULL), *gs(NULL);
2883 TObjArray *a0(NULL), *a1(NULL);
2884 a0 = (TObjArray*)((TObjArray*)fGraphM->At(ip))->At(idx);
2885 a1 = (TObjArray*)((TObjArray*)fGraphS->At(ip))->At(idx);
2886 if(!n) n=a0->GetEntriesFast();
2887 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'));
2888 Int_t nn(0), nPlots(0);
2889 for(Int_t is(0), is0(0); is<n; is++){
2890 is0 = sel ? sel[is] : is;
2891 if(!(gs = (TGraphErrors*)a1->At(is0))) return kFALSE;
2892 if(!(gm = (TGraphErrors*)a0->At(is0))) return kFALSE;
2894 if((nn=gs->GetN())){
2898 //printf("LegEntry %s [%s]%s\n", at[0], gs->GetName(), gs->GetTitle());
2899 legS->AddEntry(gs, gs->GetTitle(), "pl");
2901 gs->Sort(&TGraph::CompareY);
2902 PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[kMCtrack], at[0]), gs->GetY()[0]);
2903 PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[kMCtrack], at[0]), gs->GetY()[nn-1]);
2904 gs->Sort(&TGraph::CompareX);
2910 //printf("LegEntry %s [%s]%s\n", at[0], gm->GetName(), gm->GetTitle());
2911 legM->AddEntry(gm, gm->GetTitle(), "pl");
2913 PutTrendValue(Form("%s_%s", fgPerformanceName[kMCtrack], at[0]), gm->GetMean(2));
2914 PutTrendValue(Form("%s_%sRMS", fgPerformanceName[kMCtrack], at[0]), gm->GetRMS(2));
2917 if(!nPlots) return kFALSE;
2925 //________________________________________________________
2926 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
2929 // Get the most probable value and the full width half mean
2930 // of a Landau distribution
2933 const Float_t dx = 1.;
2934 mpv = f->GetParameter(1);
2935 Float_t fx, max = f->Eval(mpv);
2938 while((fx = f->Eval(xm))>.5*max){
2947 while((fx = f->Eval(xM))>.5*max) xM += dx;
2951 //________________________________________________________
2952 void AliTRDresolution::SetSegmentationLevel(Int_t l)
2954 // Setting the segmentation level to "l"
2957 UShort_t const lNcomp[kNprojs] = {
2959 fgkNresYsegm[fSegmentLevel], 2, //2,
2960 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2961 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2962 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2964 fgkNresYsegm[fSegmentLevel], 2, //2,
2965 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2966 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, 1, 1, 1, 11, 11, 11, //11
2967 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, 1, 1, 1, 11, 11, 11, //11
2968 6*fgkNresYsegm[fSegmentLevel], 6*2, 6*2, 6*2, 6, 6, 6, 6, 6*11, 6*11, 6*11 //11
2970 memcpy(fNcomp, lNcomp, kNprojs*sizeof(UShort_t));
2972 Char_t const *lAxTitle[kNprojs][4] = {
2974 {"Impv", "x [cm]", "I_{mpv}", "x/x_{0}"}
2975 ,{"dI/Impv", "x/x_{0}", "#delta I/I_{mpv}", "x[cm]"}
2976 // Clusters to Kalman
2977 ,{"Cluster2Track residuals", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
2978 ,{"Cluster2Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2979 // TRD tracklet to Kalman fit
2980 ,{"Tracklet2Track Y residuals", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
2981 ,{"Tracklet2Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2982 ,{"Tracklet2Track Z residuals", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
2983 ,{"Tracklet2Track Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
2984 ,{"Tracklet2Track Phi residuals", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
2985 // TRDin 2 first TRD tracklet
2986 ,{"Tracklet2Track Y residuals @ TRDin", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
2987 ,{"Tracklet2Track YZ pulls @ TRDin", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2988 ,{"Tracklet2Track Z residuals @ TRDin", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
2989 ,{"Tracklet2Track Z pulls @ TRDin", "tg(#theta)", "z", "#sigma_{z}"}
2990 ,{"Tracklet2Track Phi residuals @ TRDin", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
2991 // TRDout 2 first TRD tracklet
2992 ,{"Tracklet2Track Y residuals @ TRDout", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
2993 ,{"Tracklet2Track YZ pulls @ TRDout", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
2994 ,{"Tracklet2Track Z residuals @ TRDout", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
2995 ,{"Tracklet2Track Z pulls @ TRDout", "tg(#theta)", "z", "#sigma_{z}"}
2996 ,{"Tracklet2Track Phi residuals @ TRDout", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
2998 ,{"MC Cluster Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
2999 ,{"MC Cluster YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3001 ,{"MC Tracklet Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
3002 ,{"MC Tracklet YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3003 ,{"MC Tracklet Z resolution", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3004 ,{"MC Tracklet Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
3005 ,{"MC Tracklet Phi resolution", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3007 ,{"MC Y resolution @ TRDin", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
3008 ,{"MC YZ pulls @ TRDin", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3009 ,{"MC Z resolution @ TRDin", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3010 ,{"MC Z pulls @ TRDin", "tg(#theta)", "z", "#sigma_{z}"}
3011 ,{"MC #Phi resolution @ TRDin", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3012 ,{"MC SNP pulls @ TRDin", "tg(#phi)", "SNP", "#sigma_{snp}"}
3013 ,{"MC #Theta resolution @ TRDin", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3014 ,{"MC TGL pulls @ TRDin", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3015 ,{"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}) [%]"}
3016 ,{"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}"}
3017 ,{"MC P resolution @ TRDin", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
3019 ,{"MC Y resolution @ TRDout", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
3020 ,{"MC YZ pulls @ TRDout", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3021 ,{"MC Z resolution @ TRDout", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3022 ,{"MC Z pulls @ TRDout", "tg(#theta)", "z", "#sigma_{z}"}
3023 ,{"MC #Phi resolution @ TRDout", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3024 ,{"MC SNP pulls @ TRDout", "tg(#phi)", "SNP", "#sigma_{snp}"}
3025 ,{"MC #Theta resolution @ TRDout", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3026 ,{"MC TGL pulls @ TRDout", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3027 ,{"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}) [%]"}
3028 ,{"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}"}
3029 ,{"MC P resolution @ TRDout", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
3031 ,{"MC Track Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3032 ,{"MC Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3033 ,{"MC Track Z resolution", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3034 ,{"MC Track Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
3035 ,{"MC Track #Phi resolution", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3036 ,{"MC Track SNP pulls", "tg(#phi)", "SNP", "#sigma_{snp}"}
3037 ,{"MC Track #Theta resolution", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3038 ,{"MC Track TGL pulls", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3039 ,{"MC P_{t} resolution", "p_{t} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "#sigma(#Deltap_{t}/p_{t}^{MC}) [%]"}
3040 ,{"MC 1/P_{t} pulls", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC} - 1/p_{t}^{MC}", "#sigma_{1/p_{t}}"}
3041 ,{"MC P resolution", "p [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "#sigma(#Deltap/p^{MC}) [%]"}
3043 memcpy(fAxTitle, lAxTitle, 4*kNprojs*sizeof(Char_t*));