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 "AliCDBManager.h"
76 #include "AliCDBPath.h"
77 #include "AliCDBEntry.h"
78 #include "AliGeomManager.h"
79 #include "AliMathBase.h"
81 #include "AliTRDresolution.h"
82 #include "AliTRDgeometry.h"
83 #include "AliTRDpadPlane.h"
84 #include "AliTRDcluster.h"
85 #include "AliTRDseedV1.h"
86 #include "AliTRDtrackV1.h"
87 #include "AliTRDReconstructor.h"
88 #include "AliTRDrecoParam.h"
89 #include "AliTRDpidUtil.h"
91 #include "info/AliTRDclusterInfo.h"
93 ClassImp(AliTRDresolution)
95 UChar_t const AliTRDresolution::fgNproj[kNviews] = {
99 Char_t const * AliTRDresolution::fgPerformanceName[kNviews] = {
111 // Configure segmentation for y resolution/residuals
112 Int_t const AliTRDresolution::fgkNresYsegm[3] = {
113 AliTRDgeometry::kNsector
114 ,AliTRDgeometry::kNsector*AliTRDgeometry::kNstack
115 ,AliTRDgeometry::kNdet
117 Char_t const *AliTRDresolution::fgkResYsegmName[3] = {
118 "Sector", "Stack", "Detector"};
121 //________________________________________________________
122 AliTRDresolution::AliTRDresolution()
128 ,fReconstructor(NULL)
139 // Default constructor
141 SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
142 SetSegmentationLevel();
145 //________________________________________________________
146 AliTRDresolution::AliTRDresolution(char* name)
147 :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
152 ,fReconstructor(NULL)
163 // Default constructor
166 fReconstructor = new AliTRDReconstructor();
167 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
168 fGeo = new AliTRDgeometry();
171 SetSegmentationLevel();
173 DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
174 DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
175 /* DefineOutput(kTrkltToTrk, TObjArray::Class()); // tracklet2track
176 DefineOutput(kTrkltToMC, TObjArray::Class()); // tracklet2mc*/
179 //________________________________________________________
180 AliTRDresolution::~AliTRDresolution()
186 if(fGraphS){fGraphS->Delete(); delete fGraphS;}
187 if(fGraphM){fGraphM->Delete(); delete fGraphM;}
189 delete fReconstructor;
190 if(gGeoManager) delete gGeoManager;
191 if(fCl){fCl->Delete(); delete fCl;}
192 if(fMCcl){fMCcl->Delete(); delete fMCcl;}
193 /* if(fTrklt){fTrklt->Delete(); delete fTrklt;}
194 if(fMCtrklt){fMCtrklt->Delete(); delete fMCtrklt;}*/
198 //________________________________________________________
199 void AliTRDresolution::UserCreateOutputObjects()
201 // spatial resolution
204 fReconstructor = new AliTRDReconstructor();
205 fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
208 AliTRDrecoTask::UserCreateOutputObjects();
209 InitExchangeContainers();
212 //________________________________________________________
213 void AliTRDresolution::InitExchangeContainers()
215 fCl = new TObjArray();
216 fCl->SetOwner(kTRUE);
217 fMCcl = new TObjArray();
218 fMCcl->SetOwner(kTRUE);
219 /* fTrklt = new TObjArray();
220 fTrklt->SetOwner(kTRUE);
221 fMCtrklt = new TObjArray();
222 fMCtrklt->SetOwner(kTRUE);*/
223 PostData(kClToTrk, fCl);
224 PostData(kClToMC, fMCcl);
225 /* PostData(kTrkltToTrk, fTrklt);
226 PostData(kTrkltToMC, fMCtrklt);*/
229 //________________________________________________________
230 void AliTRDresolution::UserExec(Option_t *opt)
238 AliCDBManager* ocdb = AliCDBManager::Instance();
239 AliCDBEntry* obj = ocdb->Get(AliCDBPath("GRP", "Geometry", "Data"));
240 AliGeomManager::SetGeometry((TGeoManager*)obj->GetObject());
241 AliGeomManager::GetNalignable("TRD");
242 AliGeomManager::ApplyAlignObjsFromCDB("TRD");
243 fGeo = new AliTRDgeometry();
250 fMCtrklt->Delete();*/
251 AliTRDrecoTask::UserExec(opt);
254 //________________________________________________________
255 Bool_t AliTRDresolution::Pulls(Double_t dyz[2], Double_t cov[3], Double_t tilt)
257 // Helper function to calculate pulls in the yz plane
258 // using proper tilt rotation
259 // Uses functionality defined by AliTRDseedV1.
261 Double_t t2(tilt*tilt);
265 cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
266 cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
267 cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
269 Double_t sqr[3]={0., 0., 0.};
270 if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
271 Double_t invsqr[3]={0., 0., 0.};
272 if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
273 Double_t tmp(dyz[0]);
274 dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
275 dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
279 //________________________________________________________
280 TH1* AliTRDresolution::PlotCharge(const AliTRDtrackV1 *track)
283 // Plots the charge distribution
286 if(track) fkTrack = track;
288 AliDebug(4, "No Track defined.");
291 TObjArray *arr = NULL;
292 if(!fContainer || !(arr = ((TObjArray*)fContainer->At(kCharge)))){
293 AliWarning("No output container defined.");
298 AliTRDseedV1 *fTracklet = NULL;
299 AliTRDcluster *c = NULL;
300 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
301 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
302 if(!fTracklet->IsOK()) continue;
303 Float_t x0 = fTracklet->GetX0();
305 for(Int_t itb=AliTRDseedV1::kNtb; itb--;){
306 if(!(c = fTracklet->GetClusters(itb))){
307 if(!(c = fTracklet->GetClusters(AliTRDseedV1::kNtb+itb))) continue;
309 dq = fTracklet->GetdQdl(itb, &dl);
310 dl /= 0.15; // dl/dl0, dl0 = 1.5 mm for nominal vd
311 (h = (TH3S*)arr->At(0))->Fill(dl, x0-c->GetX(), dq);
314 // if(!HasMCdata()) continue;
316 // Float_t pt0, y0, z0, dydx0, dzdx0;
317 // if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
324 //________________________________________________________
325 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
328 // Plot the cluster distributions
331 if(track) fkTrack = track;
333 AliDebug(4, "No Track defined.");
336 TObjArray *arr = NULL;
337 if(!fContainer || !(arr = ((TObjArray*)fContainer->At(kCluster)))){
338 AliWarning("No output container defined.");
341 ULong_t status = fkESD ? fkESD->GetStatus():0;
344 Double_t covR[7], cov[3], dy[2], dz[2];
345 Float_t pt, x0, y0, z0, dydx, dzdx;
346 AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
347 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
348 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
349 if(!fTracklet->IsOK()) continue;
350 x0 = fTracklet->GetX0();
351 pt = fTracklet->GetPt();
352 sgm[2] = fTracklet->GetDetector();
353 sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
354 sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
356 // retrive the track angle with the chamber
357 y0 = fTracklet->GetYref(0);
358 z0 = fTracklet->GetZref(0);
359 dydx = fTracklet->GetYref(1);
360 dzdx = fTracklet->GetZref(1);
361 fTracklet->GetCovRef(covR);
362 Double_t tilt(fTracklet->GetTilt())
365 ,cost(TMath::Sqrt(corr));
366 AliTRDcluster *c = NULL;
367 fTracklet->ResetClusterIter(kFALSE);
368 while((c = fTracklet->PrevCluster())){
369 Float_t xc = c->GetX();
370 Float_t yc = c->GetY();
371 Float_t zc = c->GetZ();
372 Float_t dx = x0 - xc;
373 Float_t yt = y0 - dx*dydx;
374 Float_t zt = z0 - dx*dzdx;
375 dy[0] = yc-yt; dz[0]= zc-zt;
378 dy[1] = cost*(dy[0] - dz[0]*tilt);
379 dz[1] = cost*(dz[0] + dy[0]*tilt);
380 if(pt>fPtThreshold && c->IsInChamber()) ((TH3S*)arr->At(0))->Fill(dydx, dy[1], sgm[fSegmentLevel]);
382 // tilt rotation of covariance for clusters
383 Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
384 cov[0] = (sy2+t2*sz2)*corr;
385 cov[1] = tilt*(sz2 - sy2)*corr;
386 cov[2] = (t2*sy2+sz2)*corr;
387 // sum with track covariance
388 cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
389 Double_t dyz[2]= {dy[1], dz[1]};
390 Pulls(dyz, cov, tilt);
391 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
393 // Get z-position with respect to anode wire
394 Int_t istk = fGeo->GetStack(c->GetDetector());
395 AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, istk);
396 Float_t row0 = pp->GetRow0();
397 Float_t d = row0 - zt + pp->GetAnodeWireOffset();
398 d -= ((Int_t)(2 * d)) / 2.0;
399 if (d > 0.25) d = 0.5 - d;
401 AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
402 clInfo->SetCluster(c);
403 Float_t covF[] = {cov[0], cov[1], cov[2]};
404 clInfo->SetGlobalPosition(yt, zt, dydx, dzdx, covF);
405 clInfo->SetResolution(dy[1]);
406 clInfo->SetAnisochronity(d);
407 clInfo->SetDriftLength(dx);
408 clInfo->SetTilt(tilt);
409 if(fCl) fCl->Add(clInfo);
410 else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
414 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
415 clInfoArr->SetOwner(kFALSE);
417 clInfoArr->Add(clInfo);
420 if(DebugLevel()>=1 && clInfoArr){
421 (*DebugStream()) << "cluster"
422 <<"status=" << status
423 <<"clInfo.=" << clInfoArr
428 if(clInfoArr) delete clInfoArr;
429 return (TH3S*)arr->At(0);
433 //________________________________________________________
434 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
436 // Plot normalized residuals for tracklets to track.
438 // We start from the result that if X=N(|m|, |Cov|)
440 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
442 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
443 // reference position.
444 if(track) fkTrack = track;
446 AliDebug(4, "No Track defined.");
449 TObjArray *arr = NULL;
450 if(!fContainer || !(arr = (TObjArray*)fContainer->At(kTrack ))){
451 AliWarning("No output container defined.");
456 Double_t cov[3], covR[7]/*, sqr[3], inv[3]*/;
457 Double_t pt, phi, tht, x, dx, dy[2], dz[2];
458 AliTRDseedV1 *fTracklet(NULL);
459 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
460 if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
461 if(!fTracklet->IsOK()) continue;
462 sgm[2] = fTracklet->GetDetector();
463 sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
464 sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
465 x = fTracklet->GetX();
466 dx = fTracklet->GetX0() - x;
467 pt = fTracklet->GetPt();
468 phi = fTracklet->GetYref(1);
469 tht = fTracklet->GetZref(1);
471 dy[0]= fTracklet->GetYref(0)-dx*fTracklet->GetYref(1) - fTracklet->GetY();
472 dz[0]= fTracklet->GetZref(0)-dx*fTracklet->GetZref(1) - fTracklet->GetZ();
473 Double_t tilt(fTracklet->GetTilt())
476 ,cost(TMath::Sqrt(corr));
477 Bool_t rc(fTracklet->IsRowCross());
479 // calculate residuals using tilt rotation
480 dy[1]= cost*(dy[0] - dz[0]*tilt);
481 dz[1]= cost*(dz[0] + dy[0]*tilt);
482 ((TH3S*)arr->At(0))->Fill(phi, dy[1], sgm[fSegmentLevel]+rc*fgkNresYsegm[fSegmentLevel]);
483 ((TH3S*)arr->At(2))->Fill(tht, dz[1], rc);
485 // compute covariance matrix
486 fTracklet->GetCovAt(x, cov);
487 fTracklet->GetCovRef(covR);
488 cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
489 Double_t dyz[2]= {dy[1], dz[1]};
490 Pulls(dyz, cov, tilt);
491 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
492 ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
494 Double_t dphi((phi-fTracklet->GetYfit(1))/(1-phi*fTracklet->GetYfit(1)));
495 Double_t dtht((tht-fTracklet->GetZfit(1))/(1-tht*fTracklet->GetZfit(1)));
496 ((TH2I*)arr->At(4))->Fill(phi, TMath::ATan(dphi));
499 UChar_t err(fTracklet->GetErrorMsg());
500 (*DebugStream()) << "tracklet"
520 return (TH2I*)arr->At(0);
524 //________________________________________________________
525 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
527 // Store resolution/pulls of Kalman before updating with the TRD information
528 // at the radial position of the first tracklet. The following points are used
530 // - the (y,z,snp) of the first TRD tracklet
531 // - the (y, z, snp, tgl, pt) of the MC track reference
533 // Additionally the momentum resolution/pulls are calculated for usage in the
536 if(track) fkTrack = track;
538 AliDebug(4, "No Track defined.");
541 TObjArray *arr = NULL;
542 if(!fContainer || !(arr = (TObjArray*)fContainer->At(kTrackIn))){
543 AliWarning("No output container defined.");
546 AliExternalTrackParam *tin = NULL;
547 if(!(tin = fkTrack->GetTrackIn())){
548 AliWarning("Track did not entered TRD fiducial volume.");
553 Double_t x = tin->GetX();
554 AliTRDseedV1 *fTracklet = NULL;
555 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
556 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
559 if(!fTracklet || TMath::Abs(x-fTracklet->GetX())>1.e-3){
560 AliWarning("Tracklet did not match Track.");
564 sgm[2] = fTracklet->GetDetector();
565 sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
566 sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
567 Double_t tilt(fTracklet->GetTilt())
570 ,cost(TMath::Sqrt(corr));
571 Bool_t rc(fTracklet->IsRowCross());
573 const Int_t kNPAR(5);
574 Double_t parR[kNPAR]; memcpy(parR, tin->GetParameter(), kNPAR*sizeof(Double_t));
575 Double_t covR[3*kNPAR]; memcpy(covR, tin->GetCovariance(), 3*kNPAR*sizeof(Double_t));
576 Double_t cov[3]; fTracklet->GetCovAt(x, cov);
578 // define sum covariances
579 TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
580 Double_t *pc = &covR[0], *pp = &parR[0];
581 for(Int_t ir=0; ir<kNPAR; ir++, pp++){
583 for(Int_t ic = 0; ic<=ir; ic++,pc++){
584 COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
587 PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
588 //COV.Print(); PAR.Print();
590 //TODO Double_t dydx = TMath::Sqrt(1.-parR[2]*parR[2])/parR[2];
591 Double_t dy[2]={parR[0] - fTracklet->GetY(), 0.}
592 ,dz[2]={parR[1] - fTracklet->GetZ(), 0.}
593 ,dphi(TMath::ASin(PAR[2])-TMath::ATan(fTracklet->GetYfit(1)));
594 // calculate residuals using tilt rotation
595 dy[1] = cost*(dy[0] - dz[0]*tilt);
596 dz[1] = cost*(dz[0] + dy[0]*tilt);
598 if(1./PAR[4]>fPtThreshold) ((TH3S*)arr->At(0))->Fill(fTracklet->GetYref(1), dy[1], sgm[fSegmentLevel]+rc*fgkNresYsegm[fSegmentLevel]);
599 ((TH3S*)arr->At(2))->Fill(fTracklet->GetZref(1), dz[1], rc);
600 ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), dphi);
602 Double_t dyz[2] = {dy[1], dz[1]};
603 Double_t cc[3] = {COV(0,0)+cov[0], COV(0,1)+cov[1], COV(1,1)+cov[2]};
604 Pulls(dyz, cc, tilt);
605 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
606 ((TH3S*)arr->At(3))->Fill(fTracklet->GetZref(1), dyz[1], rc);
610 // register reference histo for mini-task
611 h = (TH2I*)arr->At(0);
614 (*DebugStream()) << "trackIn"
620 Double_t y = fTracklet->GetY();
621 Double_t z = fTracklet->GetZ();
622 (*DebugStream()) << "trackletIn"
632 if(!HasMCdata()) return h;
634 Float_t dx, pt0, x0=fTracklet->GetX0(), y0, z0, dydx0, dzdx0;
635 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
636 Int_t pdg = fkMC->GetPDG(),
637 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
639 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
640 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
641 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
643 // translate to reference radial position
644 dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
645 Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi))
647 TVectorD PARMC(kNPAR);
648 PARMC[0]=y0; PARMC[1]=z0;
649 PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
652 // TMatrixDSymEigen eigen(COV);
653 // TVectorD evals = eigen.GetEigenValues();
654 // TMatrixDSym evalsm(kNPAR);
655 // for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
656 // TMatrixD evecs = eigen.GetEigenVectors();
657 // TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
660 if(!(arr = (TObjArray*)fContainer->At(kMCtrackIn))) {
661 AliWarning("No MC container defined.");
665 // y resolution/pulls
666 if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0], sgm[fSegmentLevel]);
667 ((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)));
668 // z resolution/pulls
669 ((TH3S*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1], 0);
670 ((TH3S*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)), 0);
671 // phi resolution/snp pulls
672 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
673 ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
674 // theta resolution/tgl pulls
675 ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
676 ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
677 // pt resolution\\1/pt pulls\\p resolution/pull
678 ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., sign*sIdx);
679 ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), sign*sIdx);
681 Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
682 p = TMath::Sqrt(1.+ PAR[3]*PAR[3])/PAR[4];
683 ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., sign*sIdx);
685 // p*p*PAR[4]*PAR[4]*COV(4,4)
686 // +2.*PAR[3]*COV(3,4)/PAR[4]
687 // +PAR[3]*PAR[3]*COV(3,3)/p/p/PAR[4]/PAR[4]/PAR[4]/PAR[4];
688 // if(sp>0.) ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/TMath::Sqrt(sp), sign*sIdx);
692 (*DebugStream()) << "trackInMC"
699 //________________________________________________________
700 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
702 // Store resolution/pulls of Kalman after last update with the TRD information
703 // at the radial position of the first tracklet. The following points are used
705 // - the (y,z,snp) of the first TRD tracklet
706 // - the (y, z, snp, tgl, pt) of the MC track reference
708 // Additionally the momentum resolution/pulls are calculated for usage in the
711 if(track) fkTrack = track;
713 AliDebug(4, "No Track defined.");
716 TObjArray *arr = NULL;
717 if(!fContainer || !(arr = (TObjArray*)fContainer->At(kTrackOut))){
718 AliWarning("No output container defined.");
721 AliExternalTrackParam *tout = NULL;
722 if(!(tout = fkTrack->GetTrackOut())){
723 AliDebug(2, "Track did not exit TRD.");
728 Double_t x = tout->GetX();
729 AliTRDseedV1 *fTracklet(NULL);
730 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
731 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
734 if(!fTracklet || TMath::Abs(x-fTracklet->GetX())>1.e-3){
735 AliWarning("Tracklet did not match Track position.");
739 sgm[2] = fTracklet->GetDetector();
740 sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
741 sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
742 Double_t tilt(fTracklet->GetTilt())
745 ,cost(TMath::Sqrt(corr));
746 Bool_t rc(fTracklet->IsRowCross());
748 const Int_t kNPAR(5);
749 Double_t parR[kNPAR]; memcpy(parR, tout->GetParameter(), kNPAR*sizeof(Double_t));
750 Double_t covR[3*kNPAR]; memcpy(covR, tout->GetCovariance(), 3*kNPAR*sizeof(Double_t));
751 Double_t cov[3]; fTracklet->GetCovAt(x, cov);
753 // define sum covariances
754 TMatrixDSym COV(kNPAR); TVectorD PAR(kNPAR);
755 Double_t *pc = &covR[0], *pp = &parR[0];
756 for(Int_t ir=0; ir<kNPAR; ir++, pp++){
758 for(Int_t ic = 0; ic<=ir; ic++,pc++){
759 COV(ir,ic) = (*pc); COV(ic,ir) = (*pc);
762 PAR[4] = TMath::Abs(PAR[4]); // remove sign of pt !!
763 //COV.Print(); PAR.Print();
765 //TODO Double_t dydx = TMath::Sqrt(1.-parR[2]*parR[2])/parR[2];
766 Double_t dy[3]={parR[0] - fTracklet->GetY(), 0., 0.}
767 ,dz[3]={parR[1] - fTracklet->GetZ(), 0., 0.}
768 ,dphi(TMath::ASin(PAR[2])-TMath::ATan(fTracklet->GetYfit(1)));
769 // calculate residuals using tilt rotation
770 dy[1] = cost*(dy[0] - dz[0]*tilt);
771 dz[1] = cost*(dz[0] + dy[0]*tilt);
773 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 !!!
774 ((TH3S*)arr->At(2))->Fill(fTracklet->GetZref(1), dz[1], rc);
775 ((TH2I*)arr->At(4))->Fill(fTracklet->GetYref(1), dphi);
777 Double_t dyz[2] = {dy[1], dz[1]};
778 Double_t cc[3] = {COV(0,0)+cov[0], COV(0,1)+cov[1], COV(1,1)+cov[2]};
779 Pulls(dyz, cc, tilt);
780 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
781 ((TH3S*)arr->At(3))->Fill(fTracklet->GetZref(1), dyz[1], rc);
783 // register reference histo for mini-task
784 h = (TH2I*)arr->At(0);
787 (*DebugStream()) << "trackOut"
793 Double_t y = fTracklet->GetY();
794 Double_t z = fTracklet->GetZ();
795 (*DebugStream()) << "trackletOut"
805 if(!HasMCdata()) return h;
807 Float_t dx, pt0, x0=fTracklet->GetX0(), y0, z0, dydx0, dzdx0;
808 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h;
809 Int_t pdg = fkMC->GetPDG(),
810 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
812 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
813 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
814 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
816 // translate to reference radial position
817 dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0;
818 Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi))
820 TVectorD PARMC(kNPAR);
821 PARMC[0]=y0; PARMC[1]=z0;
822 PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm;
825 // TMatrixDSymEigen eigen(COV);
826 // TVectorD evals = eigen.GetEigenValues();
827 // TMatrixDSym evalsm(kNPAR);
828 // for(Int_t ir=0; ir<kNPAR; ir++) for(Int_t ic=0; ic<kNPAR; ic++) evalsm(ir,ic) = (ir==ic ? evals(ir): 0.);
829 // TMatrixD evecs = eigen.GetEigenVectors();
830 // TMatrixD sqrcov(evecs, TMatrixD::kMult, TMatrixD(evalsm, TMatrixD::kMult, evecs.T()));
833 if(!(arr = (TObjArray*)fContainer->At(kMCtrackOut))){
834 AliWarning("No MC container defined.");
837 // y resolution/pulls
838 if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, PARMC[0]-PAR[0], sgm[fSegmentLevel]);
839 ((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)));
840 // z resolution/pulls
841 ((TH3S*)arr->At(2))->Fill(dzdx0, PARMC[1]-PAR[1], 0);
842 ((TH3S*)arr->At(3))->Fill(dzdx0, (PARMC[1]-PAR[1])/TMath::Sqrt(COV(1,1)), 0);
843 // phi resolution/snp pulls
844 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ASin(PARMC[2])-TMath::ASin(PAR[2]));
845 ((TH2I*)arr->At(5))->Fill(dydx0, (PARMC[2]-PAR[2])/TMath::Sqrt(COV(2,2)));
846 // theta resolution/tgl pulls
847 ((TH2I*)arr->At(6))->Fill(dzdx0, TMath::ATan((PARMC[3]-PAR[3])/(1-PARMC[3]*PAR[3])));
848 ((TH2I*)arr->At(7))->Fill(dzdx0, (PARMC[3]-PAR[3])/TMath::Sqrt(COV(3,3)));
849 // pt resolution\\1/pt pulls\\p resolution/pull
850 ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., sign*sIdx);
851 ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), sign*sIdx);
853 Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p;
854 p = TMath::Sqrt(1.+ PAR[3]*PAR[3])/PAR[4];
855 ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., sign*sIdx);
857 // p*p*PAR[4]*PAR[4]*COV(4,4)
858 // +2.*PAR[3]*COV(3,4)/PAR[4]
859 // +PAR[3]*PAR[3]*COV(3,3)/p/p/PAR[4]/PAR[4]/PAR[4]/PAR[4];
860 // if(sp>0.) ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/TMath::Sqrt(sp), sign*sIdx);
864 (*DebugStream()) << "trackOutMC"
871 //________________________________________________________
872 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
875 // Plot MC distributions
879 AliDebug(2, "No MC defined. Results will not be available.");
882 if(track) fkTrack = track;
884 AliDebug(4, "No Track defined.");
888 AliWarning("No output container defined.");
891 // retriev track characteristics
892 Int_t pdg = fkMC->GetPDG(),
893 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
896 label(fkMC->GetLabel());
897 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
898 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
899 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
901 TObjArray *arr(NULL);TH1 *h(NULL);
903 Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl;
904 Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
905 Double_t covR[7]/*, cov[3]*/;
908 TVectorD dX(12), dY(12), dZ(12), Pt(12), dPt(12), cCOV(12*15);
909 fkMC->PropagateKalman(&dX, &dY, &dZ, &Pt, &dPt, &cCOV);
910 (*DebugStream()) << "MCkalman"
921 AliTRDReconstructor rec;
922 AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
923 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
924 if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
925 !fTracklet->IsOK())*/ continue;
927 sgm[2] = fTracklet->GetDetector();
928 sgm[0] = AliTRDgeometry::GetSector(sgm[2]);
929 sgm[1] = sgm[0] * AliTRDgeometry::kNstack + AliTRDgeometry::GetStack(sgm[2]);
930 Double_t tilt(fTracklet->GetTilt())
933 ,cost(TMath::Sqrt(corr));
934 x0 = fTracklet->GetX0();
935 //radial shift with respect to the MC reference (radial position of the pad plane)
936 x= fTracklet->GetX();
937 Bool_t rc(fTracklet->IsRowCross());
938 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
939 xAnode = fTracklet->GetX0();
941 // MC track position at reference radial position
944 (*DebugStream()) << "MC"
956 Float_t ymc = y0 - dx*dydx0;
957 Float_t zmc = z0 - dx*dzdx0;
958 //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p
960 // Kalman position at reference radial position
962 dydx = fTracklet->GetYref(1);
963 dzdx = fTracklet->GetZref(1);
964 dzdl = fTracklet->GetTgl();
965 y = fTracklet->GetYref(0) - dx*dydx;
967 z = fTracklet->GetZref(0) - dx*dzdx;
969 pt = TMath::Abs(fTracklet->GetPt());
970 fTracklet->GetCovRef(covR);
972 arr = (TObjArray*)((TObjArray*)fContainer->At(kMCtrack))->At(ily);
973 // y resolution/pulls
974 if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
975 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(covR[0]), dz/TMath::Sqrt(covR[2]));
976 // z resolution/pulls
977 ((TH3S*)arr->At(2))->Fill(dzdx0, dz, 0);
978 ((TH3S*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(covR[2]), 0);
979 // phi resolution/ snp pulls
980 Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0);
981 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp));
982 Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0);
983 ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3]));
984 // theta resolution/ tgl pulls
985 Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
986 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
987 ((TH2I*)arr->At(6))->Fill(dzdl0,
989 ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
990 // pt resolution \\ 1/pt pulls \\ p resolution for PID
991 Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
992 p = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
993 ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
994 ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
995 ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);
997 // Fill Debug stream for Kalman track
999 (*DebugStream()) << "MCtrack"
1006 << "s2y=" << covR[0]
1007 << "s2z=" << covR[2]
1011 // recalculate tracklet based on the MC info
1012 AliTRDseedV1 tt(*fTracklet);
1013 tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
1014 tt.SetZref(1, dzdx0);
1015 tt.SetReconstructor(&rec);
1017 x= tt.GetX();y= tt.GetY();z= tt.GetZ();
1018 dydx = tt.GetYfit(1);
1020 ymc = y0 - dx*dydx0;
1021 zmc = z0 - dx*dzdx0;
1024 Float_t dphi = (dydx - dydx0);
1025 dphi /= (1.- dydx*dydx0);
1027 // add tracklet residuals for y and dydx
1028 arr = (TObjArray*)fContainer->At(kMCtracklet);
1030 if(pt0>fPtThreshold) ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
1031 if(tt.GetS2Y()>0. && tt.GetS2Z()>0.) ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(tt.GetS2Y()), dz/TMath::Sqrt(tt.GetS2Z()));
1032 ((TH3S*)arr->At(2))->Fill(dzdl0, dz, rc);
1033 if(tt.GetS2Z()>0.) ((TH3S*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z()), rc);
1034 ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi));
1036 // Fill Debug stream for tracklet
1037 if(DebugLevel()>=4){
1038 Float_t s2y = tt.GetS2Y();
1039 Float_t s2z = tt.GetS2Z();
1040 (*DebugStream()) << "MCtracklet"
1051 AliTRDpadPlane *pp = fGeo->GetPadPlane(ily, AliTRDgeometry::GetStack(sgm[2]));
1052 Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
1053 //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
1055 arr = (TObjArray*)fContainer->At(kMCcluster);
1056 AliTRDcluster *c = NULL;
1057 tt.ResetClusterIter(kFALSE);
1058 while((c = tt.PrevCluster())){
1059 Float_t q = TMath::Abs(c->GetQ());
1060 x = c->GetX(); y = c->GetY();z = c->GetZ();
1064 dy = cost*(y - ymc - tilt*(z-zmc));
1065 dz = cost*(z - zmc + tilt*(y-ymc));
1068 if(q>20. && q<250. && pt0>fPtThreshold && c->IsInChamber()){
1069 ((TH3S*)arr->At(0))->Fill(dydx0, dy, sgm[fSegmentLevel]);
1070 ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dy/TMath::Sqrt(c->GetSigmaY2()), dz/TMath::Sqrt(c->GetSigmaZ2()));
1073 // Fill calibration container
1074 Float_t d = zr0 - zmc;
1075 d -= ((Int_t)(2 * d)) / 2.0;
1076 if (d > 0.25) d = 0.5 - d;
1077 AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
1078 clInfo->SetCluster(c);
1079 clInfo->SetMC(pdg, label);
1080 clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
1081 clInfo->SetResolution(dy);
1082 clInfo->SetAnisochronity(d);
1083 clInfo->SetDriftLength(dx);
1084 clInfo->SetTilt(tilt);
1085 if(fMCcl) fMCcl->Add(clInfo);
1086 else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
1087 if(DebugLevel()>=5){
1089 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
1090 clInfoArr->SetOwner(kFALSE);
1092 clInfoArr->Add(clInfo);
1096 if(DebugLevel()>=5 && clInfoArr){
1097 (*DebugStream()) << "MCcluster"
1098 <<"clInfo.=" << clInfoArr
1103 if(clInfoArr) delete clInfoArr;
1108 //________________________________________________________
1109 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
1112 // Get the reference figures
1115 Float_t xy[4] = {0., 0., 0., 0.};
1117 AliWarning("Please provide a canvas to draw results.");
1120 Int_t selection[100], n(0), selStart(0); //
1121 Int_t ly0(0), dly(5);
1122 //Int_t ly0(1), dly(2); // used for SA
1123 TList *l(NULL); TVirtualPad *pad(NULL);
1124 TGraphErrors *g(NULL);TGraphAsymmErrors *ga(NULL);
1126 case 0: // charge resolution
1127 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1128 ((TVirtualPad*)l->At(0))->cd();
1129 ga=((TGraphAsymmErrors*)((TObjArray*)fGraphM->At(kCharge))->At(0));
1130 if(ga->GetN()) ga->Draw("apl");
1131 ((TVirtualPad*)l->At(1))->cd();
1132 g = ((TGraphErrors*)((TObjArray*)fGraphS->At(kCharge))->At(0));
1133 if(g->GetN()) g->Draw("apl");
1135 case 1: // cluster2track residuals
1136 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1137 xy[0] = -.3; xy[1] = -100.; xy[2] = .3; xy[3] = 1000.;
1138 pad = (TVirtualPad*)l->At(0); pad->cd();
1139 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1140 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1141 if(!GetGraphArray(xy, kCluster, 0, 1, n, selection)) break;
1142 pad=(TVirtualPad*)l->At(1); pad->cd();
1143 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1144 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1145 if(!GetGraphArray(xy, kCluster, 0, 1, n, selection)) break;
1147 case 2: // cluster2track residuals
1148 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1149 xy[0] = -.3; xy[1] = -100.; xy[2] = .3; xy[3] = 1000.;
1150 pad = (TVirtualPad*)l->At(0); pad->cd();
1151 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1152 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1153 if(!GetGraphArray(xy, kCluster, 0, 1, n, selection)) break;
1154 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-0.5; xy[3] = 2.5;
1155 pad=(TVirtualPad*)l->At(1); pad->cd();
1156 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1157 if(!GetGraphArray(xy, kCluster, 1, 1)) break;
1160 gPad->Divide(3, 2, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1161 xy[0] = -.3; xy[1] = -20.; xy[2] = .3; xy[3] = 100.;
1162 ((TVirtualPad*)l->At(0))->cd();
1163 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1164 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection)) break;
1166 ((TVirtualPad*)l->At(1))->cd();
1167 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1168 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection)) break;
1170 ((TVirtualPad*)l->At(2))->cd();
1171 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1172 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection)) break;
1174 ((TVirtualPad*)l->At(3))->cd();
1175 selStart=fgkNresYsegm[fSegmentLevel]; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1176 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection, "[RC]")) break;
1178 ((TVirtualPad*)l->At(4))->cd();
1179 selStart=fgkNresYsegm[fSegmentLevel]/3+fgkNresYsegm[fSegmentLevel]; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1180 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection, "[RC]")) break;
1182 ((TVirtualPad*)l->At(5))->cd();
1183 selStart=2*fgkNresYsegm[fSegmentLevel]/3+fgkNresYsegm[fSegmentLevel]; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1184 if(!GetGraphArray(xy, kTrack, 0, 1, n, selection, "[RC]")) break;
1187 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1189 xy[0] = -1.; xy[1] = -150.; xy[2] = 1.; xy[3] = 1000.;
1190 ((TVirtualPad*)l->At(0))->cd();
1192 if(!GetGraphArray(xy, kTrack, 2, 1, 1, selection)) break;
1194 xy[0] = -1.; xy[1] = -1500.; xy[2] = 1.; xy[3] = 10000.;
1195 ((TVirtualPad*)l->At(1))->cd();
1197 if(!GetGraphArray(xy, kTrack, 2, 1, 1, selection)) break;
1200 case 5: // kTrack pulls
1201 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1203 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1204 ((TVirtualPad*)l->At(0))->cd();
1205 if(!GetGraphArray(xy, kTrack, 1, 1)) break;
1207 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1208 ((TVirtualPad*)l->At(1))->cd();
1209 if(!GetGraphArray(xy, kTrack, 3, 1)) break;
1211 case 6: // kTrack phi
1212 xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
1213 if(GetGraph(&xy[0], kTrack , 4)) return kTRUE;
1215 case 7: // kTrackIn y
1216 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1217 xy[0] = -.3; xy[1] = -1500.; xy[2] = .3; xy[3] = 5000.;
1218 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1219 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1220 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1221 if(!GetGraphArray(xy, kTrackIn, 0, 1, n, selection)) break;
1222 pad=((TVirtualPad*)l->At(1)); pad->cd();
1223 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1224 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1225 if(!GetGraphArray(xy, kTrackIn, 0, 1, n, selection)) break;
1227 case 8: // kTrackIn y
1228 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1229 xy[0] = -.3; xy[1] = -1500.; xy[2] = .3; xy[3] = 5000.;
1230 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1231 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1232 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1233 if(!GetGraphArray(xy, kTrackIn, 0, 1, n, selection)) break;
1234 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1235 pad=((TVirtualPad*)l->At(1)); pad->cd();
1236 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1237 if(!GetGraphArray(xy, kTrackIn, 1, 1)) break;
1239 case 9: // kTrackIn z
1240 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1241 xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
1242 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1243 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1245 if(!GetGraphArray(xy, kTrackIn, 2, 1, 1, selection)) break;
1246 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1247 pad = ((TVirtualPad*)l->At(1)); pad->cd();
1248 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1249 if(!GetGraphArray(xy, kTrackIn, 3, 1)) break;
1251 case 10: // kTrackIn phi
1252 xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
1253 if(GetGraph(&xy[0], kTrackIn, 4)) return kTRUE;
1255 case 11: // kTrackOut y
1256 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1257 xy[0] = -.3; xy[1] = -50.; xy[2] = .3; xy[3] = 150.;
1258 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1259 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1260 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1261 if(!GetGraphArray(xy, kTrackOut, 0, 1, n, selection)) break;
1262 pad=((TVirtualPad*)l->At(1)); pad->cd();
1263 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1264 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1265 if(!GetGraphArray(xy, kTrackOut, 0, 1, n, selection)) break;
1267 case 12: // kTrackOut y
1268 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1269 xy[0] = -.3; xy[1] = -50.; xy[2] = .3; xy[3] = 150.;
1270 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1271 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1272 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1273 if(!GetGraphArray(xy, kTrackOut, 0, 1, n, selection)) break;
1274 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1275 pad=((TVirtualPad*)l->At(1)); pad->cd();
1276 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1277 if(!GetGraphArray(xy, kTrackOut, 1, 1)) break;
1279 case 13: // kTrackOut z
1280 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1281 xy[0] = -1.; xy[1] = -1000.; xy[2] = 1.; xy[3] = 4000.;
1282 pad = ((TVirtualPad*)l->At(0)); pad->cd();
1283 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1284 if(!GetGraphArray(xy, kTrackOut, 2, 1)) break;
1285 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1286 pad = ((TVirtualPad*)l->At(1)); pad->cd();
1287 pad->SetMargin(0.1, 0.1, 0.1, 0.01);
1288 if(!GetGraphArray(xy, kTrackOut, 3, 1)) break;
1290 case 14: // kTrackOut phi
1291 xy[0] = -.3; xy[1] = -5.; xy[2] = .3; xy[3] = 50.;
1292 if(GetGraph(&xy[0], kTrackOut, 4)) return kTRUE;
1294 case 15: // kMCcluster
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]=650.;
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, kMCcluster, 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, kMCcluster, 0, 1, n, selection)) break;
1304 case 16: // kMCcluster
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]=650.;
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, kMCcluster, 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, kMCcluster, 1, 1)) break;
1314 case 17: //kMCtracklet [y]
1315 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1316 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =500.;
1317 ((TVirtualPad*)l->At(0))->cd();
1318 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1319 if(!GetGraphArray(xy, kMCtracklet, 0, 1, n, selection)) break;
1320 ((TVirtualPad*)l->At(1))->cd();
1321 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1322 if(!GetGraphArray(xy, kMCtracklet, 0, 1, n, selection)) break;
1324 case 18: //kMCtracklet [y]
1325 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1326 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =500.;
1327 ((TVirtualPad*)l->At(0))->cd();
1328 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1329 if(!GetGraphArray(xy, kMCtracklet, 0, 1, n, selection)) break;
1330 ((TVirtualPad*)l->At(1))->cd();
1331 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1332 if(!GetGraphArray(xy, kMCtracklet, 1, 1)) break;
1334 case 19: //kMCtracklet [z]
1335 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1336 xy[0]=-1.; xy[1]=-100.; xy[2]=1.; xy[3] =2500.;
1337 ((TVirtualPad*)l->At(0))->cd();
1338 if(!GetGraphArray(xy, kMCtracklet, 2)) break;
1339 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1340 ((TVirtualPad*)l->At(1))->cd();
1341 if(!GetGraphArray(xy, kMCtracklet, 3)) break;
1343 case 20: //kMCtracklet [phi]
1344 xy[0]=-.3; xy[1]=-3.; xy[2]=.3; xy[3] =25.;
1345 if(!GetGraph(&xy[0], kMCtracklet, 4)) break;
1347 case 21: //kMCtrack [y] ly [0]
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]*0.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1352 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer1")) break;
1353 ((TVirtualPad*)l->At(1))->cd();
1354 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*0.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1355 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer1")) break;
1357 case 22: //kMCtrack [y] ly [1]
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]*1.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1362 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer2")) break;
1363 ((TVirtualPad*)l->At(1))->cd();
1364 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*1.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1365 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer2")) break;
1367 case 23: //kMCtrack [y] ly [2]
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]*2.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1372 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer3")) break;
1373 ((TVirtualPad*)l->At(1))->cd();
1374 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*2.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1375 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer3")) break;
1377 case 24: //kMCtrack [y] ly [3]
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]*3.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1382 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer4")) break;
1383 ((TVirtualPad*)l->At(1))->cd();
1384 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*3.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1385 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer4")) break;
1387 case 25: //kMCtrack [y] ly [4]
1388 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1389 xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1390 ((TVirtualPad*)l->At(0))->cd();
1391 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*4.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1392 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer5")) break;
1393 ((TVirtualPad*)l->At(1))->cd();
1394 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*4.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1395 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer5")) break;
1397 case 26: //kMCtrack [y] ly [5]
1398 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1399 xy[0]=-.2; xy[1]=-50.; xy[2]=.2; xy[3] =400.;
1400 ((TVirtualPad*)l->At(0))->cd();
1401 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*5.); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1402 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer6")) break;
1403 ((TVirtualPad*)l->At(1))->cd();
1404 selStart=Int_t(fgkNresYsegm[fSegmentLevel]*5.5); for(n=0; n<fgkNresYsegm[fSegmentLevel]/2; n++) selection[n]=selStart+n;
1405 if(!GetGraphArray(xy, kMCtrack, 0, 1, n, selection, "Layer6")) break;
1407 case 27: //kMCtrack [y pulls]
1408 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1409 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 5.5;
1410 ((TVirtualPad*)l->At(0))->cd();
1411 selStart=0; for(n=0; n<6; n++) selection[n]=selStart+n;
1412 if(!GetGraphArray(xy, kMCtrack, 1, 1, n, selection)) break;
1413 ((TVirtualPad*)l->At(1))->cd();
1414 selStart=6; for(n=0; n<6; n++) selection[n]=selStart+n;
1415 if(!GetGraphArray(xy, kMCtrack, 1, 1, n, selection)) break;
1417 case 28: //kMCtrack [z]
1418 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1419 xy[0]=-1.; xy[1]=-1500.; xy[2]=1.; xy[3] =6000.;
1420 ((TVirtualPad*)l->At(0))->cd();
1421 if(!GetGraphArray(xy, kMCtrack, 2)) break;
1422 xy[0] = -1.; xy[1] = -1.5; xy[2] = 1.; xy[3] = 5.;
1423 ((TVirtualPad*)l->At(1))->cd();
1424 if(!GetGraphArray(xy, kMCtrack, 3)) break;
1426 case 29: //kMCtrack [phi/snp]
1427 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1428 xy[0]=-.2; xy[1]=-0.5; xy[2]=.2; xy[3] =10.;
1429 ((TVirtualPad*)l->At(0))->cd();
1430 if(!GetGraphArray(xy, kMCtrack, 4)) break;
1431 xy[0] = -.2; xy[1] = -1.5; xy[2] = .2; xy[3] = 5.;
1432 ((TVirtualPad*)l->At(1))->cd();
1433 if(!GetGraphArray(xy, kMCtrack, 5)) break;
1435 case 30: //kMCtrack [theta/tgl]
1436 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1437 xy[0]=-1.; xy[1]=-0.5; xy[2]=1.; xy[3] =5.;
1438 ((TVirtualPad*)l->At(0))->cd();
1439 if(!GetGraphArray(xy, kMCtrack, 6)) break;
1440 xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5;
1441 ((TVirtualPad*)l->At(1))->cd();
1442 if(!GetGraphArray(xy, kMCtrack, 7)) break;
1444 case 31: //kMCtrack [pt]
1445 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1446 pad = (TVirtualPad*)l->At(0); pad->cd();
1447 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1450 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1451 selection[n++] = il*11 + 2; // pi-
1452 selection[n++] = il*11 + 8; // pi+
1454 xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 4.;
1455 //xy[0] = 0.2; xy[1] = -1.; xy[2] = 7.; xy[3] = 10.; // SA
1456 if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "#pi#pm")) break;
1457 pad->Modified(); pad->Update(); pad->SetLogx();
1458 pad = (TVirtualPad*)l->At(1); pad->cd();
1459 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1462 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1463 selection[n++] = il*11 + 3; // mu-
1464 selection[n++] = il*11 + 7; // mu+
1466 if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "#mu#pm")) break;
1467 pad->Modified(); pad->Update(); pad->SetLogx();
1469 case 32: //kMCtrack [pt]
1470 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1471 pad = (TVirtualPad*)l->At(0); pad->cd();
1472 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1475 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1476 selection[n++] = il*11 + 0; // p bar
1477 selection[n++] = il*11 + 10; // p
1479 xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 8.;
1480 //xy[0] = 0.2; xy[1] = -1.; xy[2] = 7.; xy[3] = 10.; // SA
1481 if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "p&p bar")) break;
1482 pad->Modified(); pad->Update(); pad->SetLogx();
1483 pad = (TVirtualPad*)l->At(1); pad->cd();
1484 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1487 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1488 selection[n++] = il*11 + 4; // e-
1489 selection[n++] = il*11 + 6; // e+
1491 xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 12.;
1492 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 14.; // SA
1493 if(!GetGraphArray(xy, kMCtrack, 8, kTRUE, n, selection, "e#pm")) break;
1494 pad->Modified(); pad->Update(); pad->SetLogx();
1496 case 33: //kMCtrack [1/pt] pulls
1497 xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1498 //xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 4.5; // SA
1499 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1500 pad = (TVirtualPad*)l->At(0); pad->cd();
1501 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1504 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1505 selection[n++] = il*11 + 2; // pi-
1506 selection[n++] = il*11 + 8; // pi+
1508 if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "#pi#pm")) break;
1509 pad = (TVirtualPad*)l->At(1); pad->cd();
1510 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1513 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1514 selection[n++] = il*11 + 3; // mu-
1515 selection[n++] = il*11 + 7; // mu+
1517 if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "#mu#pm")) break;
1519 case 34: //kMCtrack [1/pt] pulls
1520 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1521 pad = (TVirtualPad*)l->At(0); pad->cd();
1522 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1525 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1526 selection[n++] = il*11 + 0; // p bar
1527 selection[n++] = il*11 + 10; // p
1529 xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1530 //xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 6.; // SA
1531 if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "p & p bar")) break;
1532 pad = (TVirtualPad*)l->At(1); pad->cd();
1533 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1536 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1537 selection[n++] = il*11 + 4; // e-
1538 selection[n++] = il*11 + 6; // e+
1540 xy[0] = 0.; xy[1] = -2.; xy[2] = 2.; xy[3] = 4.5;
1541 if(!GetGraphArray(xy, kMCtrack, 9, kTRUE, n, selection, "e#pm")) break;
1543 case 35: //kMCtrack [p]
1544 xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 4.;
1545 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 10.;
1546 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1547 pad = (TVirtualPad*)l->At(0); pad->cd();
1548 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1551 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1552 selection[n++] = il*11 + 2; // pi-
1553 selection[n++] = il*11 + 8; // pi+
1555 if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "#pi#pm")) break;
1556 pad->Modified(); pad->Update(); pad->SetLogx();
1557 pad = (TVirtualPad*)l->At(1); pad->cd();
1558 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1561 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1562 selection[n++] = il*11 + 3; // mu-
1563 selection[n++] = il*11 + 7; // mu+
1565 if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "#mu#pm")) break;
1566 pad->Modified(); pad->Update(); pad->SetLogx();
1568 case 36: //kMCtrack [p]
1569 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1570 pad = (TVirtualPad*)l->At(0); pad->cd();
1571 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1574 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1575 selection[n++] = il*11 + 0; // p bar
1576 selection[n++] = il*11 + 10; // p
1578 xy[0] = 0.2; xy[1] = -.7; xy[2] = 7.; xy[3] = 8.;
1579 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 12.; // SA
1580 if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "p & p bar")) break;
1581 pad->Modified(); pad->Update(); pad->SetLogx();
1582 pad = (TVirtualPad*)l->At(1); pad->cd();
1583 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1586 for(Int_t il(ly0); il<AliTRDgeometry::kNlayer; il+=dly){
1587 selection[n++] = il*11 + 4; // e-
1588 selection[n++] = il*11 + 6; // e+
1590 xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 12.;
1591 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 14.; // SA
1592 if(!GetGraphArray(xy, kMCtrack, 10, kTRUE, n, selection, "e#pm")) break;
1593 pad->Modified(); pad->Update(); pad->SetLogx();
1595 case 37: // kMCtrackIn [y]
1596 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1597 xy[0]=-.25; xy[1]=-1000.; xy[2]=.25; xy[3] =3000.;
1598 ((TVirtualPad*)l->At(0))->cd();
1599 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1600 if(!GetGraphArray(xy, kMCtrackIn, 0, 1, n, selection)) break;
1601 ((TVirtualPad*)l->At(1))->cd();
1602 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1603 if(!GetGraphArray(&xy[0], kMCtrackIn, 0, 1, n, selection)) break;
1605 case 38: // kMCtrackIn [y]
1606 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1607 xy[0]=-.25; xy[1]=-1000.; xy[2]=.25; xy[3] =3000.;
1608 ((TVirtualPad*)l->At(0))->cd();
1609 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1610 if(!GetGraphArray(xy, kMCtrackIn, 0, 1, n, selection)) break;
1611 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1612 ((TVirtualPad*)l->At(1))->cd();
1613 if(!GetGraphArray(xy, kMCtrackIn, 1, 1)) break;
1615 case 39: // kMCtrackIn [z]
1616 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1617 xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =800.;
1618 ((TVirtualPad*)l->At(0))->cd();
1619 if(!GetGraphArray(xy, kMCtrackIn, 2, 1)) break;
1620 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1621 ((TVirtualPad*)l->At(1))->cd();
1622 if(!GetGraphArray(xy, kMCtrackIn, 3, 1)) break;
1624 case 40: // kMCtrackIn [phi|snp]
1625 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1626 xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
1627 ((TVirtualPad*)l->At(0))->cd();
1628 if(!GetGraph(&xy[0], kMCtrackIn, 4)) break;
1629 xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
1630 ((TVirtualPad*)l->At(1))->cd();
1631 if(!GetGraph(&xy[0], kMCtrackIn, 5)) break;
1633 case 41: // kMCtrackIn [theta|tgl]
1634 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1635 xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.;
1636 ((TVirtualPad*)l->At(0))->cd();
1637 if(!GetGraph(&xy[0], kMCtrackIn, 6)) break;
1638 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 1.5;
1639 ((TVirtualPad*)l->At(1))->cd();
1640 if(!GetGraph(&xy[0], kMCtrackIn, 7)) break;
1642 case 42: // kMCtrackIn [pt]
1643 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1644 xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1645 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 10.; // SA
1646 pad=(TVirtualPad*)l->At(0); pad->cd(); pad->SetLogx();
1647 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1648 n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1649 if(!GetGraphArray(xy, kMCtrackIn, 8, 1, n, selection)) break;
1650 pad = (TVirtualPad*)l->At(1); pad->cd(); pad->SetLogx();
1651 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1652 n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1653 if(!GetGraphArray(xy, kMCtrackIn, 8, 1, n, selection)) break;
1655 case 43: //kMCtrackIn [1/pt] pulls
1656 xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1657 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1658 pad = (TVirtualPad*)l->At(0); pad->cd();
1659 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1660 n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1661 if(!GetGraphArray(xy, kMCtrackIn, 9, 1, n, selection)) break;
1662 pad = (TVirtualPad*)l->At(1); pad->cd();
1663 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1664 n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1665 if(!GetGraphArray(xy, kMCtrackIn, 9, 1, n, selection)) break;
1667 case 44: // kMCtrackIn [p]
1668 xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1669 //xy[0] = 0.2; xy[1] = -1.5; xy[2] = 7.; xy[3] = 10.;
1670 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1671 pad = ((TVirtualPad*)l->At(0));pad->cd();pad->SetLogx();
1672 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1673 n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1674 if(!GetGraphArray(xy, kMCtrackIn, 10, 1, n, selection)) break;
1675 pad = ((TVirtualPad*)l->At(1)); pad->cd();pad->SetLogx();
1676 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1677 n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1678 if(!GetGraphArray(xy, kMCtrackIn, 10, 1, n, selection)) break;
1680 case 45: // kMCtrackOut [y]
1681 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1682 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =400.;
1683 ((TVirtualPad*)l->At(0))->cd();
1684 selStart=0; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1685 if(!GetGraphArray(xy, kMCtrackOut, 0, 1, n, selection)) break;
1686 ((TVirtualPad*)l->At(1))->cd();
1687 selStart=fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1688 if(!GetGraphArray(&xy[0], kMCtrackOut, 0, 1, n, selection)) break;
1690 case 46: // kMCtrackOut [y]
1691 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1692 xy[0]=-.3; xy[1]=-50.; xy[2]=.3; xy[3] =400.;
1693 ((TVirtualPad*)l->At(0))->cd();
1694 selStart=2*fgkNresYsegm[fSegmentLevel]/3; for(n=0; n<fgkNresYsegm[fSegmentLevel]/3; n++) selection[n]=selStart+n;
1695 if(!GetGraphArray(xy, kMCtrackOut, 0, 1, n, selection)) break;
1696 xy[0] = -.5; xy[1] = -0.5; xy[2] = fgkNresYsegm[fSegmentLevel]-.5; xy[3] = 2.5;
1697 ((TVirtualPad*)l->At(1))->cd();
1698 if(!GetGraphArray(xy, kMCtrackOut, 1, 1)) break;
1700 case 47: // kMCtrackOut [z]
1701 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1702 xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =1500.;
1703 ((TVirtualPad*)l->At(0))->cd();
1704 if(!GetGraphArray(xy, kMCtrackOut, 2, 1)) break;
1705 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5;
1706 ((TVirtualPad*)l->At(1))->cd();
1707 if(!GetGraphArray(xy, kMCtrackOut, 3, 1)) break;
1709 case 48: // kMCtrackOut [phi|snp]
1710 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1711 xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5;
1712 ((TVirtualPad*)l->At(0))->cd();
1713 if(!GetGraph(&xy[0], kMCtrackOut, 4)) break;
1714 xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5;
1715 ((TVirtualPad*)l->At(1))->cd();
1716 if(!GetGraph(&xy[0], kMCtrackOut, 5)) break;
1718 case 49: // kMCtrackOut [theta|tgl]
1719 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1720 xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.;
1721 ((TVirtualPad*)l->At(0))->cd();
1722 if(!GetGraph(&xy[0], kMCtrackOut, 6)) break;
1723 xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 15.;
1724 ((TVirtualPad*)l->At(1))->cd();
1725 if(!GetGraph(&xy[0], kMCtrackOut, 7)) break;
1727 case 50: // kMCtrackOut [pt]
1728 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1729 xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1730 pad=(TVirtualPad*)l->At(0); pad->cd(); pad->SetLogx();
1731 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1732 n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1733 if(!GetGraphArray(xy, kMCtrackOut, 8, 1, n, selection)) break;
1734 pad = (TVirtualPad*)l->At(1); pad->cd();pad->SetLogx();
1735 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1736 n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1737 if(!GetGraphArray(xy, kMCtrackOut, 8, 1, n, selection)) break;
1739 case 51: //kMCtrackOut [1/pt] pulls
1740 xy[0] = 0.; xy[1] = -1.; xy[2] = 2.; xy[3] = 3.5;
1741 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1742 pad = (TVirtualPad*)l->At(0); pad->cd();
1743 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1744 n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1745 if(!GetGraphArray(xy, kMCtrackOut, 9, 1, n, selection)) break;
1746 pad = (TVirtualPad*)l->At(1); pad->cd();
1747 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1748 n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1749 if(!GetGraphArray(xy, kMCtrackOut, 9, 1, n, selection)) break;
1751 case 52: // kMCtrackOut [p]
1752 gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives();
1753 xy[0] = 0.2; xy[1] = -.8; xy[2] = 7.; xy[3] = 6.;
1754 pad = ((TVirtualPad*)l->At(0));pad->cd();pad->SetLogx();
1755 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1756 n=0; selection[n++]=2; selection[n++]=3; selection[n++]=7; selection[n++]=8;
1757 if(!GetGraphArray(xy, kMCtrackOut, 10, 1, n, selection)) break;
1758 pad = ((TVirtualPad*)l->At(1)); pad->cd();pad->SetLogx();
1759 pad->SetMargin(0.125, 0.015, 0.1, 0.015);
1760 n=0; selection[n++]=0; selection[n++]=4; selection[n++]=6; selection[n++]=10;
1761 if(!GetGraphArray(xy, kMCtrackOut, 10, 1, n, selection)) break;
1764 AliWarning(Form("Reference plot [%d] missing result", ifig));
1768 //________________________________________________________
1769 void AliTRDresolution::MakeSummary()
1771 // Build summary plots
1773 if(!fGraphS || !fGraphM){
1774 AliError("Missing results");
1777 Float_t xy[4] = {0., 0., 0., 0.};
1779 TH2 *h2 = new TH2I("h2SF", "", 20, -.2, .2, fgkNresYsegm[fSegmentLevel], -0.5, fgkNresYsegm[fSegmentLevel]-0.5);
1780 h2->GetXaxis()->CenterTitle();
1781 h2->GetYaxis()->CenterTitle();
1782 h2->GetZaxis()->CenterTitle();h2->GetZaxis()->SetTitleOffset(1.4);
1784 Int_t ih2(0), iSumPlot(0);
1785 TCanvas *cOut = new TCanvas(Form("TRDsummary%s_%d", GetName(), iSumPlot++), "Cluster & Tracklet Resolution", 1024, 768);
1786 cOut->Divide(3,2, 2.e-3, 2.e-3, kYellow-7);
1787 TVirtualPad *p(NULL);
1790 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1791 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1792 h2->SetTitle(Form("Cluster-Track R-Phi Residuals;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1793 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphS->At(kCluster))->At(0), h2);
1794 GetRange(h2, 1, range);
1795 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1800 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1801 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1802 h2->SetTitle(Form("Cluster-Track R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1803 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphM->At(kCluster))->At(0), h2);
1804 GetRange(h2, 0, range);
1805 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1810 p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1811 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1812 GetGraphArray(xy, kCluster, 1, 1);
1815 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1816 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1817 h2->SetTitle(Form("Tracklet-Track R-Phi Residuals;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1818 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphS->At(kTrack))->At(0), h2);
1819 GetRange(h2, 1, range);
1820 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1825 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1826 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1827 h2->SetTitle(Form("Tracklet-Track R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1828 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphM->At(kTrack))->At(0), h2);
1829 GetRange(h2, 0, range);
1830 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1835 p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1836 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1837 GetGraphArray(xy, kTrack, 1, 1);
1839 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1845 cOut->Clear(); cOut->SetName(Form("TRDsummary%s_%d", GetName(), iSumPlot++));
1846 cOut->Divide(3, 2, 2.e-3, 2.e-3, kBlue-10);
1849 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1850 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1851 h2->SetTitle(Form("Cluster-MC R-Phi Resolution;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1852 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphS->At(kMCcluster))->At(0), h2);
1853 GetRange(h2, 1, range);
1854 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1859 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1860 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1862 h2->SetTitle(Form("Cluster-MC R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1863 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphM->At(kMCcluster))->At(0), h2);
1864 GetRange(h2, 0, range);
1865 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1870 p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1871 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1872 GetGraphArray(xy, kMCcluster, 1, 1);
1875 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1876 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1878 h2->SetTitle(Form("Tracklet-MC R-Phi Resolution;tg(#phi);%s;Sigma [#mum]", fgkResYsegmName[fSegmentLevel]));
1879 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphS->At(kMCtracklet))->At(0), h2);
1880 GetRange(h2, 1, range);
1881 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1886 p->SetRightMargin(0.16);p->SetTopMargin(0.06);
1887 h2=(TH2I*)h2->Clone(Form("h2SF_%d", ih2++));
1889 h2->SetTitle(Form("Tracklet-MC R-Phi Systematics;tg(#phi);%s;Mean [#mum]", fgkResYsegmName[fSegmentLevel]));
1890 MakeSummaryPlot((TObjArray*) ((TObjArray*)fGraphM->At(kMCtracklet))->At(0), h2);
1891 GetRange(h2, 0, range);
1892 h2->GetZaxis()->SetRangeUser(range[0], range[1]);
1897 p->SetRightMargin(0.06);p->SetTopMargin(0.06);
1898 xy[0]=-.5; xy[1]=-0.5; xy[2]=fgkNresYsegm[fSegmentLevel]-.5; xy[3]=2.5;
1899 GetGraphArray(xy, kMCtracklet, 1, 1);
1901 cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1905 //________________________________________________________
1906 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1908 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1909 // The "range" vector should be initialized with 2 elements
1910 // Option "mod" can be any of
1911 // - 0 : gaussian like distribution
1912 // - 1 : tailed distribution
1914 Int_t nx(h2->GetNbinsX())
1915 , ny(h2->GetNbinsY())
1917 Double_t *data=new Double_t[n];
1918 for(Int_t ix(1), in(0); ix<=nx; ix++){
1919 for(Int_t iy(1); iy<=ny; iy++)
1920 data[in++] = h2->GetBinContent(ix, iy);
1922 Double_t mean, sigm;
1923 AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1925 range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1926 if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1927 AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1928 TH1S h1("h1SF0", "", 100, range[0], range[1]);
1933 case 0:// gaussian distribution
1935 TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1937 mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1938 range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1939 AliDebug(2, Form(" rangeG[%f %f]", range[0], range[1]));
1942 case 1:// tailed distribution
1944 Int_t bmax(h1.GetMaximumBin());
1945 Int_t jBinMin(1), jBinMax(100);
1946 for(Int_t ibin(bmax); ibin--;){
1947 if(h1.GetBinContent(ibin)==0){
1948 jBinMin=ibin; break;
1951 for(Int_t ibin(bmax); ibin++;){
1952 if(h1.GetBinContent(ibin)==0){
1953 jBinMax=ibin; break;
1956 range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1957 AliDebug(2, Form(" rangeT[%f %f]", range[0], range[1]));
1965 //________________________________________________________
1966 void AliTRDresolution::MakeSummaryPlot(TObjArray *a, TH2 *h2)
1970 TGraphErrors *g(NULL); TAxis *ax(h2->GetXaxis());
1971 for(Int_t iseg(0); iseg<fgkNresYsegm[fSegmentLevel]; iseg++){
1972 g=(TGraphErrors*)a->At(iseg);
1973 for(Int_t in(0); in<g->GetN(); in++){
1974 g->GetPoint(in, x, y);
1975 h2->SetBinContent(ax->FindBin(x), iseg+1, y);
1981 //________________________________________________________
1982 Char_t const *fgParticle[11]={
1983 " p bar", " K -", " #pi -", " #mu -", " e -",
1985 " e +", " #mu +", " #pi +", " K +", " p",
1987 const Color_t fgColorS[11]={
1988 kOrange, kOrange-3, kMagenta+1, kViolet, kRed,
1990 kRed, kViolet, kMagenta+1, kOrange-3, kOrange
1992 const Color_t fgColorM[11]={
1993 kCyan-5, kAzure-4, kBlue-7, kBlue+2, kViolet+10,
1995 kViolet+10, kBlue+2, kBlue-7, kAzure-4, kCyan-5
1997 const Marker_t fgMarker[11]={
2002 Bool_t AliTRDresolution::PostProcess()
2004 //fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
2006 AliError("ERROR: list not available");
2009 TGraph *gm= NULL, *gs= NULL;
2010 if(!fGraphS && !fGraphM){
2011 TObjArray *aM(NULL), *aS(NULL);
2012 Int_t n = fContainer->GetEntriesFast();
2013 fGraphS = new TObjArray(n); fGraphS->SetOwner();
2014 fGraphM = new TObjArray(n); fGraphM->SetOwner();
2015 for(Int_t ig(0), nc(0); ig<n; ig++){
2016 fGraphM->AddAt(aM = new TObjArray(fgNproj[ig]), ig);
2017 fGraphS->AddAt(aS = new TObjArray(fgNproj[ig]), ig);
2019 for(Int_t ic=0; ic<fgNproj[ig]; ic++, nc++){
2020 AliDebug(2, Form("building G[%d] P[%d] N[%d]", ig, ic, fNcomp[nc]));
2022 TObjArray *agS(NULL), *agM(NULL);
2023 aS->AddAt(agS = new TObjArray(fNcomp[nc]), ic);
2024 aM->AddAt(agM = new TObjArray(fNcomp[nc]), ic);
2025 for(Int_t is=fNcomp[nc]; is--;){
2026 agS->AddAt(gs = new TGraphErrors(), is);
2027 Int_t is0(is%11), il0(is/11);
2028 gs->SetMarkerStyle(fgMarker[is0]);
2029 gs->SetMarkerColor(fgColorS[is0]);
2030 gs->SetLineColor(fgColorS[is0]);
2031 gs->SetLineStyle(il0);gs->SetLineWidth(2);
2032 gs->SetName(Form("s_%d_%02d_%02d", ig, ic, is));
2034 agM->AddAt(gm = new TGraphErrors(), is);
2035 gm->SetMarkerStyle(fgMarker[is0]);
2036 gm->SetMarkerColor(fgColorM[is0]);
2037 gm->SetLineColor(fgColorM[is0]);
2038 gm->SetLineStyle(il0);gm->SetLineWidth(2);
2039 gm->SetName(Form("m_%d_%02d_%02d", ig, ic, is));
2040 // this is important for labels in the legend
2042 gs->SetTitle(Form("%s %02d", fgkResYsegmName[fSegmentLevel], is%fgkNresYsegm[fSegmentLevel]));
2043 gm->SetTitle(Form("%s %02d", fgkResYsegmName[fSegmentLevel], is%fgkNresYsegm[fSegmentLevel]));
2045 gs->SetTitle(Form("%s Ly[%d]", is%2 ?"z":"y", is/2));
2046 gm->SetTitle(Form("%s Ly[%d]", is%2?"z":"y", is/2));
2047 } else if(ic==2||ic==3) {
2048 gs->SetTitle(Form("%s Ly[%d]", is%2 ?"RC":"no RC", is/2));
2049 gm->SetTitle(Form("%s Ly[%d]", is%2?"RC":"no RC", is/2));
2051 gs->SetTitle(Form("Layer[%d]", is%AliTRDgeometry::kNlayer));
2052 gm->SetTitle(Form("Layer[%d]", is%AliTRDgeometry::kNlayer));
2054 gs->SetTitle(Form("%s @ ly[%d]", fgParticle[is0], il0));
2055 gm->SetTitle(Form("%s @ ly[%d]", fgParticle[is0], il0));
2059 aS->AddAt(gs = new TGraphErrors(), ic);
2060 gs->SetMarkerStyle(23);
2061 gs->SetMarkerColor(kRed);
2062 gs->SetLineColor(kRed);
2063 gs->SetNameTitle(Form("s_%d_%02d", ig, ic), "sigma");
2065 aM->AddAt(gm = ig ? (TGraph*)new TGraphErrors() : (TGraph*)new TGraphAsymmErrors(), ic);
2066 gm->SetLineColor(kBlack);
2067 gm->SetMarkerStyle(7);
2068 gm->SetMarkerColor(kBlack);
2069 gm->SetNameTitle(Form("m_%d_%02d", ig, ic), "mean");
2075 /* printf("\n\n\n"); fGraphS->ls();
2076 printf("\n\n\n"); fGraphM->ls();*/
2081 TF1 fg("fGauss", "gaus", -.5, .5);
2082 // Landau for charge resolution
2083 TF1 fch("fClCh", "landau", 0., 1000.);
2084 // Landau for e+- pt resolution
2085 TF1 fpt("fPt", "landau", -0.1, 0.2);
2087 //PROCESS EXPERIMENTAL DISTRIBUTIONS
2088 // Charge resolution
2089 //Process3DL(kCharge, 0, &fl);
2090 // Clusters residuals
2091 Process3D(kCluster, 0, &fg, 1.e4);
2092 Process3Dlinked(kCluster, 1, &fg);
2094 // Tracklet residual/pulls
2095 Process3D(kTrack , 0, &fg, 1.e4);
2096 Process3Dlinked(kTrack , 1, &fg);
2097 Process3D(kTrack , 2, &fg, 1.e4);
2098 Process3D(kTrack , 3, &fg);
2099 Process2D(kTrack , 4, &fg, 1.e3);
2101 // TRDin residual/pulls
2102 Process3D(kTrackIn, 0, &fg, 1.e4);
2103 Process3Dlinked(kTrackIn, 1, &fg);
2104 Process3D(kTrackIn, 2, &fg, 1.e4);
2105 Process3D(kTrackIn, 3, &fg);
2106 Process2D(kTrackIn, 4, &fg, 1.e3);
2108 // TRDout residual/pulls
2109 Process3D(kTrackOut, 0, &fg, 1.e3); // scale to fit - see PlotTrackOut
2110 Process3Dlinked(kTrackOut, 1, &fg);
2111 Process3D(kTrackOut, 2, &fg, 1.e4);
2112 Process3D(kTrackOut, 3, &fg);
2113 Process2D(kTrackOut, 4, &fg, 1.e3);
2116 if(!HasMCdata()) return kTRUE;
2119 //PROCESS MC RESIDUAL DISTRIBUTIONS
2121 // CLUSTER Y RESOLUTION/PULLS
2122 Process3D(kMCcluster, 0, &fg, 1.e4);
2123 Process3Dlinked(kMCcluster, 1, &fg, 1.);
2126 // TRACKLET RESOLUTION/PULLS
2127 Process3D(kMCtracklet, 0, &fg, 1.e4); // y
2128 Process3Dlinked(kMCtracklet, 1, &fg, 1.); // y pulls
2129 Process3D(kMCtracklet, 2, &fg, 1.e4); // z
2130 Process3D(kMCtracklet, 3, &fg, 1.); // z pulls
2131 Process2D(kMCtracklet, 4, &fg, 1.e3); // phi
2134 // TRACK RESOLUTION/PULLS
2135 Process3Darray(kMCtrack, 0, &fg, 1.e4); // y
2136 Process3DlinkedArray(kMCtrack, 1, &fg); // y PULL
2137 Process3Darray(kMCtrack, 2, &fg, 1.e4); // z
2138 Process3Darray(kMCtrack, 3, &fg); // z PULL
2139 Process2Darray(kMCtrack, 4, &fg, 1.e3); // phi
2140 Process2Darray(kMCtrack, 5, &fg); // snp PULL
2141 Process2Darray(kMCtrack, 6, &fg, 1.e3); // theta
2142 Process2Darray(kMCtrack, 7, &fg); // tgl PULL
2143 Process3Darray(kMCtrack, 8, &fg, 1.e2); // pt resolution
2144 Process3Darray(kMCtrack, 9, &fg); // 1/pt pulls
2145 Process3Darray(kMCtrack, 10, &fg, 1.e2); // p resolution
2148 // TRACK TRDin RESOLUTION/PULLS
2149 Process3D(kMCtrackIn, 0, &fg, 1.e4);// y resolution
2150 Process3Dlinked(kMCtrackIn, 1, &fg); // y pulls
2151 Process3D(kMCtrackIn, 2, &fg, 1.e4);// z resolution
2152 Process3D(kMCtrackIn, 3, &fg); // z pulls
2153 Process2D(kMCtrackIn, 4, &fg, 1.e3);// phi resolution
2154 Process2D(kMCtrackIn, 5, &fg); // snp pulls
2155 Process2D(kMCtrackIn, 6, &fg, 1.e3);// theta resolution
2156 Process2D(kMCtrackIn, 7, &fg); // tgl pulls
2157 Process3D(kMCtrackIn, 8, &fg, 1.e2);// pt resolution
2158 Process3D(kMCtrackIn, 9, &fg); // 1/pt pulls
2159 Process3D(kMCtrackIn, 10, &fg, 1.e2);// p resolution
2162 // TRACK TRDout RESOLUTION/PULLS
2163 Process3D(kMCtrackOut, 0, &fg, 1.e4);// y resolution
2164 Process3Dlinked(kMCtrackOut, 1, &fg); // y pulls
2165 Process3D(kMCtrackOut, 2, &fg, 1.e4);// z resolution
2166 Process3D(kMCtrackOut, 3, &fg); // z pulls
2167 Process2D(kMCtrackOut, 4, &fg, 1.e3);// phi resolution
2168 Process2D(kMCtrackOut, 5, &fg); // snp pulls
2169 Process2D(kMCtrackOut, 6, &fg, 1.e3);// theta resolution
2170 Process2D(kMCtrackOut, 7, &fg); // tgl pulls
2171 Process3D(kMCtrackOut, 8, &fg, 1.e2);// pt resolution
2172 Process3D(kMCtrackOut, 9, &fg); // 1/pt pulls
2173 Process3D(kMCtrackOut, 10, &fg, 1.e2);// p resolution
2180 //________________________________________________________
2181 void AliTRDresolution::Terminate(Option_t *opt)
2183 AliTRDrecoTask::Terminate(opt);
2184 if(HasPostProcess()) PostProcess();
2187 //________________________________________________________
2188 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
2190 // Helper function to avoid duplication of code
2191 // Make first guesses on the fit parameters
2193 // find the intial parameters of the fit !! (thanks George)
2194 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
2196 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
2197 f->SetParLimits(0, 0., 3.*sum);
2198 f->SetParameter(0, .9*sum);
2199 Double_t rms = h->GetRMS();
2200 f->SetParLimits(1, -rms, rms);
2201 f->SetParameter(1, h->GetMean());
2203 f->SetParLimits(2, 0., 2.*rms);
2204 f->SetParameter(2, rms);
2205 if(f->GetNpar() <= 4) return;
2207 f->SetParLimits(3, 0., sum);
2208 f->SetParameter(3, .1*sum);
2210 f->SetParLimits(4, -.3, .3);
2211 f->SetParameter(4, 0.);
2213 f->SetParLimits(5, 0., 1.e2);
2214 f->SetParameter(5, 2.e-1);
2217 //________________________________________________________
2218 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand)
2220 // Build performance histograms for AliTRDcluster.vs TRD track or MC
2221 // - y reziduals/pulls
2223 TObjArray *arr = new TObjArray(2);
2224 arr->SetName(name); arr->SetOwner();
2225 TH1 *h(NULL); char hname[100], htitle[300];
2227 // tracklet resolution/pull in y direction
2228 sprintf(hname, "%s_%s_Y", GetNameId(), name);
2229 sprintf(htitle, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
2230 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2231 Int_t nybins=fgkNresYsegm[fSegmentLevel];
2232 if(expand) nybins*=2;
2233 h = new TH3S(hname, htitle,
2234 48, -.48, .48, 60, -1.5, 1.5, nybins, -0.5, nybins-0.5);
2237 sprintf(hname, "%s_%s_YZpull", GetNameId(), name);
2238 sprintf(htitle, "YZ pull for \"%s\" @ %s;%s;#Delta y / #sigma_{y};#Delta z / #sigma_{z}", GetNameId(), name, fgkResYsegmName[fSegmentLevel]);
2239 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2240 h = new TH3S(hname, htitle, fgkNresYsegm[fSegmentLevel], -0.5, fgkNresYsegm[fSegmentLevel]-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
2247 //________________________________________________________
2248 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
2250 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
2251 // - y reziduals/pulls
2252 // - z reziduals/pulls
2254 TObjArray *arr = BuildMonitorContainerCluster(name, expand);
2256 TH1 *h(NULL); char hname[100], htitle[300];
2258 // tracklet resolution/pull in z direction
2259 sprintf(hname, "%s_%s_Z", GetNameId(), name);
2260 sprintf(htitle, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm];row cross", GetNameId(), name);
2261 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2262 h = new TH3S(hname, htitle, 50, -1., 1., 100, -1.5, 1.5, 2, -0.5, 1.5);
2265 sprintf(hname, "%s_%s_Zpull", GetNameId(), name);
2266 sprintf(htitle, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z / #sigma_{z};row cross", GetNameId(), name);
2267 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2268 h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
2269 h->GetZaxis()->SetBinLabel(1, "no RC");
2270 h->GetZaxis()->SetBinLabel(2, "RC");
2274 // tracklet to track phi resolution
2275 sprintf(hname, "%s_%s_PHI", GetNameId(), name);
2276 sprintf(htitle, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];entries", GetNameId(), name);
2277 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2278 h = new TH2I(hname, htitle, 21, -.33, .33, 100, -.5, .5);
2285 //________________________________________________________
2286 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
2288 // Build performance histograms for AliExternalTrackParam.vs MC
2289 // - y resolution/pulls
2290 // - z resolution/pulls
2291 // - phi resolution, snp pulls
2292 // - theta resolution, tgl pulls
2293 // - pt resolution, 1/pt pulls, p resolution
2295 TObjArray *arr = BuildMonitorContainerTracklet(name);
2297 TH1 *h(NULL); char hname[100], htitle[300];
2301 sprintf(hname, "%s_%s_SNPpull", GetNameId(), name);
2302 sprintf(htitle, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp / #sigma_{snp};entries", GetNameId(), name);
2303 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2304 h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
2309 sprintf(hname, "%s_%s_THT", GetNameId(), name);
2310 sprintf(htitle, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
2311 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2312 h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
2316 sprintf(hname, "%s_%s_TGLpull", GetNameId(), name);
2317 sprintf(htitle, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl / #sigma_{tgl};entries", GetNameId(), name);
2318 if(!(h = (TH2I*)gROOT->FindObject(hname))){
2319 h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
2323 const Int_t kNpt(14);
2324 const Int_t kNdpt(150);
2325 const Int_t kNspc = 2*AliPID::kSPECIES+1;
2326 Float_t Pt=0.1, DPt=-.1, Spc=-5.5;
2327 Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
2328 for(Int_t i=0;i<kNpt+1; i++,Pt=TMath::Exp(i*.15)-1.) binsPt[i]=Pt;
2329 for(Int_t i=0; i<kNspc+1; i++,Spc+=1.) binsSpc[i]=Spc;
2330 for(Int_t i=0; i<kNdpt+1; i++,DPt+=2.e-3) binsDPt[i]=DPt;
2333 sprintf(hname, "%s_%s_Pt", GetNameId(), name);
2334 sprintf(htitle, "P_{t} res for \"%s\" @ %s;p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
2335 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2336 h = new TH3S(hname, htitle,
2337 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2339 for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2343 sprintf(hname, "%s_%s_1Pt", GetNameId(), name);
2344 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);
2345 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2346 h = new TH3S(hname, htitle,
2347 kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
2349 for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2353 sprintf(hname, "%s_%s_P", GetNameId(), name);
2354 sprintf(htitle, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
2355 if(!(h = (TH3S*)gROOT->FindObject(hname))){
2356 h = new TH3S(hname, htitle,
2357 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
2359 for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
2367 //________________________________________________________
2368 TObjArray* AliTRDresolution::Histos()
2371 // Define histograms
2374 if(fContainer) return fContainer;
2376 fContainer = new TObjArray(kNviews);
2377 //fContainer->SetOwner(kTRUE);
2379 TObjArray *arr(NULL);
2381 // binnings for plots containing momentum or pt
2382 const Int_t kNpt(14), kNphi(48), kNdy(60);
2383 Float_t Phi=-.48, Dy=-.3, Pt=0.1;
2384 Float_t binsPhi[kNphi+1], binsDy[kNdy+1], binsPt[kNpt+1];
2385 for(Int_t i=0; i<kNphi+1; i++,Phi+=.02) binsPhi[i]=Phi;
2386 for(Int_t i=0; i<kNdy+1; i++,Dy+=.01) binsDy[i]=Dy;
2387 for(Int_t i=0;i<kNpt+1; i++,Pt=TMath::Exp(i*.15)-1.) binsPt[i]=Pt;
2389 // cluster to track residuals/pulls
2390 fContainer->AddAt(arr = new TObjArray(2), kCharge);
2391 arr->SetName("Charge");
2392 if(!(h = (TH3S*)gROOT->FindObject("hCharge"))){
2393 h = new TH3S("hCharge", "Charge Resolution", 20, 1., 2., 24, 0., 3.6, 100, 0., 500.);
2394 h->GetXaxis()->SetTitle("dx/dx_{0}");
2395 h->GetYaxis()->SetTitle("x_{d} [cm]");
2396 h->GetZaxis()->SetTitle("dq/dx [ADC/cm]");
2400 // cluster to track residuals/pulls
2401 fContainer->AddAt(BuildMonitorContainerCluster("Cl"), kCluster);
2402 // tracklet to TRD track
2403 fContainer->AddAt(BuildMonitorContainerTracklet("Trk", kTRUE), kTrack);
2404 // tracklet to TRDin
2405 fContainer->AddAt(BuildMonitorContainerTracklet("TrkIN", kTRUE), kTrackIn);
2406 // tracklet to TRDout
2407 fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
2410 // Resolution histos
2411 if(!HasMCdata()) return fContainer;
2413 // cluster resolution
2414 fContainer->AddAt(BuildMonitorContainerCluster("MCcl"), kMCcluster);
2416 // tracklet resolution
2417 fContainer->AddAt(BuildMonitorContainerTracklet("MCtracklet"), kMCtracklet);
2420 fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
2421 arr->SetName("MCtrk");
2422 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
2424 // TRDin TRACK RESOLUTION
2425 fContainer->AddAt(BuildMonitorContainerTrack("MCtrkIN"), kMCtrackIn);
2427 // TRDout TRACK RESOLUTION
2428 fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
2433 //________________________________________________________
2434 Bool_t AliTRDresolution::Load(const Char_t *file, const Char_t *dir)
2436 // Custom load function. Used to guess the segmentation level of the data.
2438 if(!AliTRDrecoTask::Load(file, dir)) return kFALSE;
2440 // look for cluster residual plot - always available
2441 TH3S* h3((TH3S*)((TObjArray*)fContainer->At(kClToTrk))->At(0));
2442 Int_t segmentation(h3->GetNbinsZ()/2);
2443 if(segmentation==fgkNresYsegm[0]){ // default segmentation. Nothing to do
2445 } else if(segmentation==fgkNresYsegm[1]){ // stack segmentation.
2446 SetSegmentationLevel(1);
2447 } else if(segmentation==fgkNresYsegm[2]){ // detector segmentation.
2448 SetSegmentationLevel(2);
2450 AliError(Form("Unknown segmentation [%d].", h3->GetNbinsZ()));
2454 AliDebug(2, Form("Segmentation set to level \"%s\"", fgkResYsegmName[fSegmentLevel]));
2459 //________________________________________________________
2460 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
2463 // Do the processing
2466 Char_t pn[10]; sprintf(pn, "p%03d", fIdxPlot);
2468 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
2469 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
2470 if(Int_t(h2->GetEntries())){
2471 AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
2473 AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
2478 const Int_t kINTEGRAL=1;
2479 for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
2480 Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
2481 Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
2482 TH1D *h = h2->ProjectionY(pn, abin, bbin);
2483 if((n=(Int_t)h->GetEntries())<150){
2484 AliDebug(4, Form(" x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
2488 Int_t ip = g[0]->GetN();
2489 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)));
2490 g[0]->SetPoint(ip, x, k*f->GetParameter(1));
2491 g[0]->SetPointError(ip, 0., k*f->GetParError(1));
2492 g[1]->SetPoint(ip, x, k*f->GetParameter(2));
2493 g[1]->SetPointError(ip, 0., k*f->GetParError(2));
2495 g[0]->SetPoint(ip, x, k*h->GetMean());
2496 g[0]->SetPointError(ip, 0., k*h->GetMeanError());
2497 g[1]->SetPoint(ip, x, k*h->GetRMS());
2498 g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
2504 //________________________________________________________
2505 Bool_t AliTRDresolution::Process2D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k, Int_t gidx)
2508 // Do the processing
2511 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2513 // retrive containers
2516 if(!(h2= (TH2I*)(fContainer->At(plot)))) return kFALSE;
2518 TObjArray *a0(NULL);
2519 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2520 if(!(h2=(TH2I*)a0->At(idx))) return kFALSE;
2522 if(Int_t(h2->GetEntries())){
2523 AliDebug(2, Form("p[%d] idx[%d] : h[%s] %s", plot, idx, h2->GetName(), h2->GetTitle()));
2525 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2530 if(gidx<0) gidx=idx;
2531 if(!(g[0] = gidx<0 ? (TGraphErrors*)fGraphM->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphM->At(plot)))->At(gidx))) return kFALSE;
2533 if(!(g[1] = gidx<0 ? (TGraphErrors*)fGraphS->At(plot) : (TGraphErrors*)((TObjArray*)(fGraphS->At(plot)))->At(gidx))) return kFALSE;
2535 return Process(h2, f, k, g);
2538 //________________________________________________________
2539 Bool_t AliTRDresolution::Process3D(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2542 // Do the processing
2545 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2547 // retrive containers
2550 if(!(h3= (TH3S*)(fContainer->At(plot)))) return kFALSE;
2552 TObjArray *a0(NULL);
2553 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2554 if(!(h3=(TH3S*)a0->At(idx))) return kFALSE;
2556 if(Int_t(h3->GetEntries())){
2557 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2559 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2564 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2565 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2568 TAxis *az = h3->GetZaxis();
2569 for(Int_t iz(0); iz<gm->GetEntriesFast(); iz++){
2570 if(!(g[0] = (TGraphErrors*)gm->At(iz))) return kFALSE;
2571 if(!(g[1] = (TGraphErrors*)gs->At(iz))) return kFALSE;
2572 az->SetRange(iz+1, iz+1);
2573 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2580 //________________________________________________________
2581 Bool_t AliTRDresolution::Process3Dlinked(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2584 // Do the processing
2587 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2589 // retrive containers
2592 if(!(h3= (TH3S*)(fContainer->At(plot)))) return kFALSE;
2594 TObjArray *a0(NULL);
2595 if(!(a0=(TObjArray*)(fContainer->At(plot)))) return kFALSE;
2596 if(!(h3=(TH3S*)a0->At(idx))) return kFALSE;
2598 if(Int_t(h3->GetEntries())){
2599 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2601 AliDebug(2, Form("p[%d] idx[%d] : Missing entries.", plot, idx));
2606 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2607 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2610 if(!(g[0] = (TGraphErrors*)gm->At(0))) return kFALSE;
2611 if(!(g[1] = (TGraphErrors*)gs->At(0))) return kFALSE;
2612 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2614 if(!(g[0] = (TGraphErrors*)gm->At(1))) return kFALSE;
2615 if(!(g[1] = (TGraphErrors*)gs->At(1))) return kFALSE;
2616 if(!Process((TH2*)h3->Project3D("zx"), f, k, g)) return kFALSE;
2622 //________________________________________________________
2623 Bool_t AliTRDresolution::Process3DL(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2626 // Do the processing
2629 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2631 // retrive containers
2632 TH3S *h3 = (TH3S*)((TObjArray*)fContainer->At(plot))->At(idx);
2633 if(!h3) return kFALSE;
2634 AliDebug(2, Form("p[%d] idx[%d] h[%s] %s", plot, idx, h3->GetName(), h3->GetTitle()));
2636 TGraphAsymmErrors *gm;
2638 if(!(gm = (TGraphAsymmErrors*)((TObjArray*)fGraphM->At(plot))->At(0))) return kFALSE;
2639 if(!(gs = (TGraphErrors*)((TObjArray*)fGraphS->At(plot)))) return kFALSE;
2641 Float_t x, r, mpv, xM, xm;
2642 TAxis *ay = h3->GetYaxis();
2643 for(Int_t iy=1; iy<=h3->GetNbinsY(); iy++){
2644 ay->SetRange(iy, iy);
2645 x = ay->GetBinCenter(iy);
2646 TH2F *h2=(TH2F*)h3->Project3D("zx");
2647 TAxis *ax = h2->GetXaxis();
2648 for(Int_t ix=1; ix<=h2->GetNbinsX(); ix++){
2649 r = ax->GetBinCenter(ix);
2650 TH1D *h1 = h2->ProjectionY("py", ix, ix);
2651 if(h1->Integral()<50) continue;
2654 GetLandauMpvFwhm(f, mpv, xm, xM);
2655 Int_t ip = gm->GetN();
2656 gm->SetPoint(ip, x, k*mpv);
2657 gm->SetPointError(ip, 0., 0., k*xm, k*xM);
2658 gs->SetPoint(ip, r, k*(xM-xm)/mpv);
2659 gs->SetPointError(ip, 0., 0.);
2666 //________________________________________________________
2667 Bool_t AliTRDresolution::Process2Darray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2670 // Do the processing
2673 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2675 // retrive containers
2676 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2677 if(!arr) return kFALSE;
2678 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2681 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2682 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2684 TGraphErrors *g[2]; TH2I *h2(NULL); TObjArray *a0(NULL);
2685 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2686 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2688 if(!(h2 = (TH2I*)a0->At(idx))) return kFALSE;
2689 if(Int_t(h2->GetEntries())){
2690 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h2->GetName(), h2->GetTitle()));
2692 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2696 if(!(g[0] = (TGraphErrors*)gm->At(ia))) return kFALSE;
2697 if(!(g[1] = (TGraphErrors*)gs->At(ia))) return kFALSE;
2698 if(!Process(h2, f, k, g)) return kFALSE;
2704 //________________________________________________________
2705 Bool_t AliTRDresolution::Process3Darray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2708 // Do the processing
2711 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2712 //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
2714 // retrive containers
2715 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2716 if(!arr) return kFALSE;
2717 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2720 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2721 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2723 TGraphErrors *g[2]; TH3S *h3(NULL); TObjArray *a0(NULL);
2725 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2726 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2728 if(!(h3 = (TH3S*)a0->At(idx))) return kFALSE;
2729 if(Int_t(h3->GetEntries())){
2730 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h3->GetName(), h3->GetTitle()));
2732 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2735 TAxis *az = h3->GetZaxis();
2736 for(Int_t iz=1; iz<=az->GetNbins(); iz++, in++){
2737 if(in >= gm->GetEntriesFast()) break;
2738 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2739 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2740 az->SetRange(iz, iz);
2741 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2744 AliDebug(2, Form("Projections [%d] from [%d]", in, gs->GetEntriesFast()));
2749 //________________________________________________________
2750 Bool_t AliTRDresolution::Process3DlinkedArray(ETRDresolutionPlot plot, Int_t idx, TF1 *f, Float_t k)
2753 // Do the processing
2756 if(!fContainer || !fGraphS || !fGraphM) return kFALSE;
2757 //printf("Process4D : processing plot[%d] idx[%d]\n", plot, idx);
2759 // retrive containers
2760 TObjArray *arr = (TObjArray*)(fContainer->At(plot));
2761 if(!arr) return kFALSE;
2762 AliDebug(2, Form("p[%d] idx[%d] arr[%s]", plot, idx, arr->GetName()));
2765 if(!(gm = (TObjArray*)((TObjArray*)(fGraphM->At(plot)))->At(idx))) return kFALSE;
2766 if(!(gs = (TObjArray*)((TObjArray*)(fGraphS->At(plot)))->At(idx))) return kFALSE;
2768 TGraphErrors *g[2]; TH3S *h3(NULL); TObjArray *a0(NULL);
2770 for(Int_t ia(0); ia<arr->GetEntriesFast(); ia++){
2771 if(!(a0 = (TObjArray*)arr->At(ia))) continue;
2772 if(!(h3 = (TH3S*)a0->At(idx))) return kFALSE;
2773 if(Int_t(h3->GetEntries())){
2774 AliDebug(4, Form(" idx[%d] h[%s] %s", ia, h3->GetName(), h3->GetTitle()));
2776 AliDebug(2, Form(" idx[%d] : Missing entries.", ia));
2779 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2780 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2781 if(!Process((TH2*)h3->Project3D("yx"), f, k, g)) return kFALSE;
2784 if(!(g[0] = (TGraphErrors*)gm->At(in))) return kFALSE;
2785 if(!(g[1] = (TGraphErrors*)gs->At(in))) return kFALSE;
2786 if(!Process((TH2*)h3->Project3D("zx"), f, k, g)) return kFALSE;
2789 AliDebug(2, Form("Projections [%d] from [%d]", in, gs->GetEntriesFast()));
2794 //________________________________________________________
2795 Bool_t AliTRDresolution::GetGraph(Float_t *bb, ETRDresolutionPlot ip, Int_t idx, Bool_t kLEG, const Char_t *explain)
2801 if(!fGraphS || !fGraphM) return kFALSE;
2802 // axis titles look up
2804 for(Int_t jp=0; jp<(Int_t)ip; jp++) nref+=fgNproj[jp];
2805 UChar_t jdx = idx<0?0:idx;
2806 for(Int_t jc=0; jc<TMath::Min(jdx,fgNproj[ip]-1); jc++) nref++;
2807 Char_t **at = fAxTitle[nref];
2809 // build legends if requiered
2812 leg=new TLegend(.65, .6, .95, .9);
2813 leg->SetBorderSize(0);
2814 leg->SetFillStyle(0);
2818 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]);
2819 h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2820 h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
2822 TAxis *ax = h1->GetXaxis();
2823 ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
2824 ax = h1->GetYaxis();
2825 ax->SetRangeUser(bb[1], bb[3]);
2826 ax->CenterTitle(); ax->SetTitleOffset(1.4);
2829 TBox *b = new TBox(-.15, bb[1], .15, bb[3]);
2830 b->SetFillStyle(3002);b->SetLineColor(0);
2831 b->SetFillColor(ip<=Int_t(kMCcluster)?kGreen:kBlue);
2834 TGraphErrors *gm = idx<0 ? (TGraphErrors*)fGraphM->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphM->At(ip)))->At(idx);
2835 if(!gm) return kFALSE;
2836 TGraphErrors *gs = idx<0 ? (TGraphErrors*)fGraphS->At(ip) : (TGraphErrors*)((TObjArray*)(fGraphS->At(ip)))->At(idx);
2837 if(!gs) return kFALSE;
2839 Int_t n(0), nPlots(0);
2840 if((n=gm->GetN())) {
2842 gm->Draw("pl"); if(leg) leg->AddEntry(gm, gm->GetTitle(), "pl");
2843 PutTrendValue(Form("%s_%s", fgPerformanceName[ip], at[0]), gm->GetMean(2));
2844 PutTrendValue(Form("%s_%sRMS", fgPerformanceName[ip], at[0]), gm->GetRMS(2));
2849 gs->Draw("pl"); if(leg) leg->AddEntry(gs, gs->GetTitle(), "pl");
2850 gs->Sort(&TGraph::CompareY);
2851 PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[ip], at[0]), gs->GetY()[0]);
2852 PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[ip], at[0]), gs->GetY()[n-1]);
2853 gs->Sort(&TGraph::CompareX);
2855 if(!nPlots) return kFALSE;
2856 if(leg) leg->Draw();
2860 //________________________________________________________
2861 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)
2867 if(!fGraphS || !fGraphM) return kFALSE;
2869 // axis titles look up
2871 for(Int_t jp(0); jp<ip; jp++) nref+=fgNproj[jp];
2873 Char_t **at = fAxTitle[nref];
2875 // build legends if requiered
2876 TLegend *legM(NULL), *legS(NULL);
2878 legM=new TLegend(.35, .6, .65, .9);
2879 legM->SetHeader("Mean");
2880 legM->SetBorderSize(0);
2881 legM->SetFillStyle(0);
2882 legS=new TLegend(.65, .6, .95, .9);
2883 legS->SetHeader("Sigma");
2884 legS->SetBorderSize(0);
2885 legS->SetFillStyle(0);
2889 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]);
2890 h1->SetMinimum(bb[1]);h1->SetMaximum(bb[3]);
2891 h1->SetLineColor(kBlack); h1->SetLineWidth(1);h1->SetLineStyle(2);
2893 TAxis *ax = h1->GetXaxis();
2894 ax->CenterTitle();ax->SetMoreLogLabels();ax->SetTitleOffset(1.2);
2895 ax = h1->GetYaxis();
2896 ax->SetRangeUser(bb[1], bb[3]);
2897 ax->CenterTitle(); ax->SetTitleOffset(1.4);
2900 TGraphErrors *gm(NULL), *gs(NULL);
2901 TObjArray *a0(NULL), *a1(NULL);
2902 a0 = (TObjArray*)((TObjArray*)fGraphM->At(ip))->At(idx);
2903 a1 = (TObjArray*)((TObjArray*)fGraphS->At(ip))->At(idx);
2904 if(!n) n=a0->GetEntriesFast();
2905 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'));
2906 Int_t nn(0), nPlots(0);
2907 for(Int_t is(0), is0(0); is<n; is++){
2908 is0 = sel ? sel[is] : is;
2909 if(!(gs = (TGraphErrors*)a1->At(is0))) return kFALSE;
2910 if(!(gm = (TGraphErrors*)a0->At(is0))) return kFALSE;
2912 if((nn=gs->GetN())){
2916 //printf("LegEntry %s [%s]%s\n", at[0], gs->GetName(), gs->GetTitle());
2917 legS->AddEntry(gs, gs->GetTitle(), "pl");
2919 gs->Sort(&TGraph::CompareY);
2920 PutTrendValue(Form("%s_%sSigMin", fgPerformanceName[kMCtrack], at[0]), gs->GetY()[0]);
2921 PutTrendValue(Form("%s_%sSigMax", fgPerformanceName[kMCtrack], at[0]), gs->GetY()[nn-1]);
2922 gs->Sort(&TGraph::CompareX);
2928 //printf("LegEntry %s [%s]%s\n", at[0], gm->GetName(), gm->GetTitle());
2929 legM->AddEntry(gm, gm->GetTitle(), "pl");
2931 PutTrendValue(Form("%s_%s", fgPerformanceName[kMCtrack], at[0]), gm->GetMean(2));
2932 PutTrendValue(Form("%s_%sRMS", fgPerformanceName[kMCtrack], at[0]), gm->GetRMS(2));
2935 if(!nPlots) return kFALSE;
2943 //________________________________________________________
2944 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
2947 // Get the most probable value and the full width half mean
2948 // of a Landau distribution
2951 const Float_t dx = 1.;
2952 mpv = f->GetParameter(1);
2953 Float_t fx, max = f->Eval(mpv);
2956 while((fx = f->Eval(xm))>.5*max){
2965 while((fx = f->Eval(xM))>.5*max) xM += dx;
2969 //________________________________________________________
2970 void AliTRDresolution::SetRecoParam(AliTRDrecoParam *r)
2973 fReconstructor->SetRecoParam(r);
2977 //________________________________________________________
2978 void AliTRDresolution::SetSegmentationLevel(Int_t l)
2980 // Setting the segmentation level to "l"
2983 UShort_t const lNcomp[kNprojs] = {
2985 fgkNresYsegm[fSegmentLevel], 2, //2,
2986 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2987 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2988 2*fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2990 fgkNresYsegm[fSegmentLevel], 2, //2,
2991 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, //5,
2992 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, 1, 1, 1, 11, 11, 11, //11
2993 fgkNresYsegm[fSegmentLevel], 2, 2, 2, 1, 1, 1, 1, 11, 11, 11, //11
2994 6*fgkNresYsegm[fSegmentLevel], 6*2, 6*2, 6*2, 6, 6, 6, 6, 6*11, 6*11, 6*11 //11
2996 memcpy(fNcomp, lNcomp, kNprojs*sizeof(UShort_t));
2998 Char_t const *lAxTitle[kNprojs][4] = {
3000 {"Impv", "x [cm]", "I_{mpv}", "x/x_{0}"}
3001 ,{"dI/Impv", "x/x_{0}", "#delta I/I_{mpv}", "x[cm]"}
3002 // Clusters to Kalman
3003 ,{"Cluster2Track residuals", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3004 ,{"Cluster2Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3005 // TRD tracklet to Kalman fit
3006 ,{"Tracklet2Track Y residuals", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3007 ,{"Tracklet2Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3008 ,{"Tracklet2Track Z residuals", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3009 ,{"Tracklet2Track Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
3010 ,{"Tracklet2Track Phi residuals", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3011 // TRDin 2 first TRD tracklet
3012 ,{"Tracklet2Track Y residuals @ TRDin", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3013 ,{"Tracklet2Track YZ pulls @ TRDin", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3014 ,{"Tracklet2Track Z residuals @ TRDin", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3015 ,{"Tracklet2Track Z pulls @ TRDin", "tg(#theta)", "z", "#sigma_{z}"}
3016 ,{"Tracklet2Track Phi residuals @ TRDin", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3017 // TRDout 2 first TRD tracklet
3018 ,{"Tracklet2Track Y residuals @ TRDout", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3019 ,{"Tracklet2Track YZ pulls @ TRDout", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3020 ,{"Tracklet2Track Z residuals @ TRDout", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3021 ,{"Tracklet2Track Z pulls @ TRDout", "tg(#theta)", "z", "#sigma_{z}"}
3022 ,{"Tracklet2Track Phi residuals @ TRDout", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3024 ,{"MC Cluster Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3025 ,{"MC Cluster YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3027 ,{"MC Tracklet Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
3028 ,{"MC Tracklet YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3029 ,{"MC Tracklet Z resolution", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3030 ,{"MC Tracklet Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
3031 ,{"MC Tracklet Phi resolution", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3033 ,{"MC Y resolution @ TRDin", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
3034 ,{"MC YZ pulls @ TRDin", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3035 ,{"MC Z resolution @ TRDin", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3036 ,{"MC Z pulls @ TRDin", "tg(#theta)", "z", "#sigma_{z}"}
3037 ,{"MC #Phi resolution @ TRDin", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3038 ,{"MC SNP pulls @ TRDin", "tg(#phi)", "SNP", "#sigma_{snp}"}
3039 ,{"MC #Theta resolution @ TRDin", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3040 ,{"MC TGL pulls @ TRDin", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3041 ,{"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}) [%]"}
3042 ,{"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}"}
3043 ,{"MC P resolution @ TRDin", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
3045 ,{"MC Y resolution @ TRDout", "tg(#phi)", "y [#mum]", "#sigma_{y}[#mum]"}
3046 ,{"MC YZ pulls @ TRDout", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3047 ,{"MC Z resolution @ TRDout", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3048 ,{"MC Z pulls @ TRDout", "tg(#theta)", "z", "#sigma_{z}"}
3049 ,{"MC #Phi resolution @ TRDout", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3050 ,{"MC SNP pulls @ TRDout", "tg(#phi)", "SNP", "#sigma_{snp}"}
3051 ,{"MC #Theta resolution @ TRDout", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3052 ,{"MC TGL pulls @ TRDout", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3053 ,{"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}) [%]"}
3054 ,{"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}"}
3055 ,{"MC P resolution @ TRDout", "p^{MC} [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "MC: #sigma^{TPC}(#Deltap/p^{MC}) [%]"}
3057 ,{"MC Track Y resolution", "tg(#phi)", "y [#mum]", "#sigma_{y} [#mum]"}
3058 ,{"MC Track YZ pulls", fgkResYsegmName[fSegmentLevel], "y / z", "#sigma_{y}"}
3059 ,{"MC Track Z resolution", "tg(#theta)", "z [#mum]", "#sigma_{z} [#mum]"}
3060 ,{"MC Track Z pulls", "tg(#theta)", "z", "#sigma_{z}"}
3061 ,{"MC Track #Phi resolution", "tg(#phi)", "#phi [mrad]", "#sigma_{#phi} [mrad]"}
3062 ,{"MC Track SNP pulls", "tg(#phi)", "SNP", "#sigma_{snp}"}
3063 ,{"MC Track #Theta resolution", "tg(#theta)", "#theta [mrad]", "#sigma_{#theta} [mrad]"}
3064 ,{"MC Track TGL pulls", "tg(#theta)", "TGL", "#sigma_{tgl}"}
3065 ,{"MC P_{t} resolution", "p_{t} [GeV/c]", "(p_{t}^{REC}-p_{t}^{MC})/p_{t}^{MC} [%]", "#sigma(#Deltap_{t}/p_{t}^{MC}) [%]"}
3066 ,{"MC 1/P_{t} pulls", "1/p_{t}^{MC} [c/GeV]", "1/p_{t}^{REC} - 1/p_{t}^{MC}", "#sigma_{1/p_{t}}"}
3067 ,{"MC P resolution", "p [GeV/c]", "(p^{REC}-p^{MC})/p^{MC} [%]", "#sigma(#Deltap/p^{MC}) [%]"}
3069 memcpy(fAxTitle, lAxTitle, 4*kNprojs*sizeof(Char_t*));