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>
58 #include <THnSparse.h>
64 #include <TGraphErrors.h>
65 #include <TGraphAsymmErrors.h>
66 #include <TLinearFitter.h>
70 #include <TTreeStream.h>
71 #include <TGeoManager.h>
72 #include <TDatabasePDG.h>
76 #include "AliESDtrack.h"
77 #include "AliMathBase.h"
78 #include "AliTrackPointArray.h"
80 #include "AliTRDresolution.h"
81 #include "AliTRDgeometry.h"
82 #include "AliTRDtransform.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"
90 #include "AliTRDinfoGen.h"
92 #include "info/AliTRDclusterInfo.h"
94 ClassImp(AliTRDresolution)
95 //ClassImp(AliTRDresolution::AliTRDresolutionProjection)
97 Int_t const AliTRDresolution::fgkNbins[kNdim] = {
98 Int_t(kNbunchCross)/*bc*/,
104 Int_t(kNcharge)*AliPID::kSPECIES+1/*chg*species*/,
106 }; //! no of bins/projection
107 Double_t const AliTRDresolution::fgkMin[kNdim] = {
114 -AliPID::kSPECIES-0.5,
116 }; //! low limits for projections
117 Double_t const AliTRDresolution::fgkMax[kNdim] = {
118 Int_t(kNbunchCross)-0.5,
124 AliPID::kSPECIES+0.5,
126 }; //! high limits for projections
127 Char_t const *AliTRDresolution::fgkTitle[kNdim] = {
136 }; //! title of projection
138 Char_t const * AliTRDresolution::fgPerformanceName[kNclasses] = {
146 // ,"Tracklet2TRDout"
149 Float_t AliTRDresolution::fgPtBin[kNpt+1];
151 //________________________________________________________
152 AliTRDresolution::AliTRDresolution()
164 // Default constructor
166 SetNameTitle("TRDresolution", "TRD spatial and momentum resolution");
167 MakePtSegmentation();
170 //________________________________________________________
171 AliTRDresolution::AliTRDresolution(char* name, Bool_t xchange)
172 :AliTRDrecoTask(name, "TRD spatial and momentum resolution")
183 // Default constructor
187 MakePtSegmentation();
189 SetUseExchangeContainers();
190 DefineOutput(kClToTrk, TObjArray::Class()); // cluster2track
191 DefineOutput(kClToMC, TObjArray::Class()); // cluster2mc
195 //________________________________________________________
196 AliTRDresolution::~AliTRDresolution()
202 if(fProj){fProj->Delete(); delete fProj;}
203 if(fCl){fCl->Delete(); delete fCl;}
204 if(fMCcl){fMCcl->Delete(); delete fMCcl;}
208 //________________________________________________________
209 void AliTRDresolution::UserCreateOutputObjects()
211 // spatial resolution
213 AliTRDrecoTask::UserCreateOutputObjects();
214 if(UseExchangeContainers()) InitExchangeContainers();
217 //________________________________________________________
218 void AliTRDresolution::InitExchangeContainers()
220 // Init containers for subsequent tasks (AliTRDclusterResolution)
222 fCl = new TObjArray(200); fCl->SetOwner(kTRUE);
223 fMCcl = new TObjArray(); fMCcl->SetOwner(kTRUE);
224 PostData(kClToTrk, fCl);
225 PostData(kClToMC, fMCcl);
228 //________________________________________________________
229 void AliTRDresolution::UserExec(Option_t *opt)
235 if(fCl) fCl->Delete();
236 if(fMCcl) fMCcl->Delete();
237 AliTRDrecoTask::UserExec(opt);
240 //________________________________________________________
241 Bool_t AliTRDresolution::Pulls(Double_t* /*dyz[2]*/, Double_t* /*cov[3]*/, Double_t /*tilt*/) const
243 // Helper function to calculate pulls in the yz plane
244 // using proper tilt rotation
245 // Uses functionality defined by AliTRDseedV1.
249 Double_t t2(tilt*tilt);
250 // exit door until a bug fix is found for AliTRDseedV1::GetCovSqrt
254 cc[0] = cov[0] - 2.*tilt*cov[1] + t2*cov[2];
255 cc[1] = cov[1]*(1.-t2) + tilt*(cov[0] - cov[2]);
256 cc[2] = t2*cov[0] + 2.*tilt*cov[1] + cov[2];
258 Double_t sqr[3]={0., 0., 0.};
259 if(AliTRDseedV1::GetCovSqrt(cc, sqr)) return kFALSE;
260 Double_t invsqr[3]={0., 0., 0.};
261 if(AliTRDseedV1::GetCovInv(sqr, invsqr)<1.e-40) return kFALSE;
262 Double_t tmp(dyz[0]);
263 dyz[0] = invsqr[0]*tmp + invsqr[1]*dyz[1];
264 dyz[1] = invsqr[1]*tmp + invsqr[2]*dyz[1];
269 //________________________________________________________
270 TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track)
273 // Plot the cluster distributions
276 if(track) fkTrack = track;
278 AliDebug(4, "No Track defined.");
281 if(TMath::Abs(fkESD->GetTOFbc()) > 1){
282 AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
285 if(fPt<fPtThreshold){
286 AliDebug(4, Form("Track with pt[%6.4f] under threshold.", fPt));
290 if(!fContainer || !(H = ((THnSparse*)fContainer->At(kCluster)))){
291 AliWarning("No output container defined.");
295 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
296 Double_t val[kNdim]; //Float_t exb, vd, t0, s2, dl, dt;
297 TObjArray *clInfoArr(NULL);
298 AliTRDseedV1 *fTracklet(NULL);
299 AliTRDcluster *c(NULL), *cc(NULL);
300 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
301 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
302 if(!fTracklet->IsOK()) continue;
303 //fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
307 val[kPt] = TMath::ATan(fTracklet->GetYref(1))*TMath::RadToDeg();
308 Float_t corr = 1./TMath::Sqrt(1.+fTracklet->GetYref(1)*fTracklet->GetYref(1)+fTracklet->GetZref(1)*fTracklet->GetZref(1));
310 Float_t padCorr(fTracklet->GetTilt()*fTracklet->GetPadLength());
311 fTracklet->ResetClusterIter(kTRUE);
312 while((c = fTracklet->NextCluster())){
313 Float_t xc(c->GetX()),
314 q(TMath::Abs(c->GetQ()));
315 if(row0<0) row0 = c->GetPadRow();
317 val[kYrez] = c->GetY() + padCorr*(c->GetPadRow() - row0) -fTracklet->GetYat(xc);
318 val[kPrez] = fTracklet->GetX0()-xc;
319 val[kZrez] = 0.; Int_t ic(0), tb(c->GetLocalTimeBin());;
320 if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
321 if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
322 if(ic) val[kZrez] /= (ic*q);
323 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0.:(TMath::Max(q*corr, Float_t(3.)));
325 /* // tilt rotation of covariance for clusters
326 Double_t sy2(c->GetSigmaY2()), sz2(c->GetSigmaZ2());
327 cov[0] = (sy2+t2*sz2)*corr;
328 cov[1] = tilt*(sz2 - sy2)*corr;
329 cov[2] = (t2*sy2+sz2)*corr;
330 // sum with track covariance
331 cov[0]+=covR[0]; cov[1]+=covR[1]; cov[2]+=covR[2];
332 Double_t dyz[2]= {dy[1], dz[1]};
333 Pulls(dyz, cov, tilt);*/
335 // Get z-position with respect to anode wire
336 Float_t yt(fTracklet->GetYref(0)-val[kZrez]*fTracklet->GetYref(1)),
337 zt(fTracklet->GetZref(0)-val[kZrez]*fTracklet->GetZref(1));
338 Int_t istk = geo->GetStack(c->GetDetector());
339 AliTRDpadPlane *pp = geo->GetPadPlane(ily, istk);
340 Float_t rowZ = pp->GetRow0();
341 Float_t d = rowZ - zt + pp->GetAnodeWireOffset();
342 d -= ((Int_t)(2 * d)) / 2.0;
343 if (d > 0.25) d = 0.5 - d;
345 AliTRDclusterInfo *clInfo(NULL);
346 clInfo = new AliTRDclusterInfo;
347 clInfo->SetCluster(c);
348 //Float_t covF[] = {cov[0], cov[1], cov[2]};
349 clInfo->SetGlobalPosition(yt, zt, fTracklet->GetYref(1), fTracklet->GetZref(1)/*, covF*/);
350 clInfo->SetResolution(val[kYrez]);
351 clInfo->SetAnisochronity(d);
352 clInfo->SetDriftLength(val[kZrez]);
353 clInfo->SetTilt(fTracklet->GetTilt());
354 if(fCl) fCl->Add(clInfo);
355 //else AliDebug(1, "Cl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
359 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
360 clInfoArr->SetOwner(kFALSE);
362 clInfoArr->Add(clInfo);
365 if(DebugLevel()>=2 && clInfoArr){
366 ULong_t status = fkESD->GetStatus();
367 (*DebugStream()) << "cluster"
368 <<"status=" << status
369 <<"clInfo.=" << clInfoArr
374 if(clInfoArr) delete clInfoArr;
376 return NULL;//H->Projection(kEta, kPhi);
380 //________________________________________________________
381 TH1* AliTRDresolution::PlotTracklet(const AliTRDtrackV1 *track)
383 // Plot normalized residuals for tracklets to track.
385 // We start from the result that if X=N(|m|, |Cov|)
387 // (Cov^{-1})^{1/2}X = N((Cov^{-1})^{1/2}*|m|, |1|)
389 // in our case X=(y_trklt - y_trk z_trklt - z_trk) and |Cov| = |Cov_trklt| + |Cov_trk| at the radial
390 // reference position.
391 if(track) fkTrack = track;
393 AliDebug(4, "No Track defined.");
396 if(TMath::Abs(fkESD->GetTOFbc())>1){
397 AliDebug(4, Form("Track with BC_index[%d] not used.", fkESD->GetTOFbc()));
401 if(!fContainer || !(H = (THnSparse*)fContainer->At(kTracklet))){
402 AliWarning("No output container defined.");
406 Double_t val[kNdim+1];
407 AliTRDseedV1 *fTracklet(NULL);
408 for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++){
409 if(!(fTracklet = fkTrack->GetTracklet(il))) continue;
410 if(!fTracklet->IsOK()) continue;
414 val[kSpeciesChgRC]= fSpecies;
415 val[kPt] = GetPtBin(fTracklet->GetMomentum());
416 Double_t dyt(fTracklet->GetYfit(0) - fTracklet->GetYref(0)),
417 dzt(fTracklet->GetZfit(0) - fTracklet->GetZref(0)),
418 dydx(fTracklet->GetYfit(1)),
419 tilt(fTracklet->GetTilt());
420 // correct for tilt rotation
421 val[kYrez] = dyt - dzt*tilt;
422 val[kZrez] = dzt + dyt*tilt;
423 dydx+= tilt*fTracklet->GetZref(1);
424 val[kPrez] = TMath::ATan((dydx - fTracklet->GetYref(1))/(1.+ fTracklet->GetYref(1)*dydx)) * TMath::RadToDeg();
425 if(fTracklet->IsRowCross()){
426 val[kSpeciesChgRC]= 0.;
427 // val[kPrez] = fkTrack->Charge(); // may be better defined
429 Float_t exb, vd, t0, s2, dl, dt;
430 fTracklet->GetCalibParam(exb, vd, t0, s2, dl, dt);
431 val[kZrez] = TMath::ATan((fTracklet->GetYref(1) - exb)/(1+fTracklet->GetYref(1)*exb));
433 val[kNdim] = fTracklet->GetdQdl();
434 if(DebugLevel()>=1) H->Fill(val);
436 // // compute covariance matrix
437 // fTracklet->GetCovAt(x, cov);
438 // fTracklet->GetCovRef(covR);
439 // cov[0] += covR[0]; cov[1] += covR[1]; cov[2] += covR[2];
440 // Double_t dyz[2]= {dy[1], dz[1]};
441 // Pulls(dyz, cov, tilt);
442 // ((TH3S*)arr->At(1))->Fill(sgm[fSegmentLevel], dyz[0], dyz[1]);
443 // ((TH3S*)arr->At(3))->Fill(tht, dyz[1], rc);
446 Bool_t rc(fTracklet->IsRowCross());
447 UChar_t err(fTracklet->GetErrorMsg());
448 Double_t x(fTracklet->GetX()),
449 pt(fTracklet->GetPt()),
450 yt(fTracklet->GetYref(0)),
451 zt(fTracklet->GetZref(0)),
452 phi(fTracklet->GetYref(1)),
453 tht(fTracklet->GetZref(1));
454 Int_t ncl(fTracklet->GetN()),
455 det(fTracklet->GetDetector());
456 (*DebugStream()) << "tracklet"
467 <<"dy=" << val[kYrez]
468 <<"dz=" << val[kZrez]
469 <<"dphi="<< val[kPrez]
470 <<"dQ ="<< val[kNdim]
476 return NULL;//H->Projection(kEta, kPhi);
480 //________________________________________________________
481 TH1* AliTRDresolution::PlotTrackIn(const AliTRDtrackV1 *track)
483 // Store resolution/pulls of Kalman before updating with the TRD information
484 // at the radial position of the first tracklet. The following points are used
486 // - the (y,z,snp) of the first TRD tracklet
487 // - the (y, z, snp, tgl, pt) of the MC track reference
489 // Additionally the momentum resolution/pulls are calculated for usage in the
491 //printf("AliTRDresolution::PlotTrackIn() :: track[%p]\n", (void*)track);
493 if(track) fkTrack = track;
495 AliDebug(4, "No Track defined.");
500 THnSparseI *H=(THnSparseI*)fContainer->At(kTrackIn);
502 AliError(Form("Missing container @ %d", Int_t(kTrackIn)));
505 // check input track status
506 AliExternalTrackParam *tin(NULL);
507 if(!(tin = fkTrack->GetTrackIn())){
508 AliError("Track did not entered TRD fiducial volume.");
511 // check first tracklet
512 AliTRDseedV1 *fTracklet(fkTrack->GetTracklet(0));
514 AliDebug(3, "No Tracklet in ly[0]. Skip track.");
517 // check radial position
518 Double_t x = tin->GetX();
519 if(TMath::Abs(x-fTracklet->GetX())>1.e-3){
520 AliDebug(1, Form("Tracklet did not match Track. dx[cm]=%+4.1f", x-fTracklet->GetX()));
523 //printf("USE y[%+f] dydx[%+f]\n", fTracklet->GetYfit(0), fTracklet->GetYfit(1));
525 Int_t bc(TMath::Abs(fkESD->GetTOFbc())%2);
526 const Double_t *parR(tin->GetParameter());
527 Double_t dyt(fTracklet->GetYfit(0)-parR[0]), dzt(fTracklet->GetZfit(0)-parR[1]),
528 phit(fTracklet->GetYfit(1)),
529 tilt(fTracklet->GetTilt()),
530 norm(1./TMath::Sqrt((1.-parR[2])*(1.+parR[2])));
532 // correct for tilt rotation
533 Double_t dy = dyt - dzt*tilt,
535 dx = dy/(parR[2]*norm-parR[3]*norm*tilt);
536 phit += tilt*parR[3];
537 Double_t dphi = TMath::ATan(phit) - TMath::ASin(parR[2]);
539 Double_t val[kNdim+2];
543 val[kSpeciesChgRC]= fTracklet->IsRowCross()?0:fkTrack->Charge();
544 val[kPt] = fPt<0.8?0:(fPt<1.5?1:2);//GetPtBin(fPt);
547 val[kPrez] = dphi*TMath::RadToDeg();
548 val[kNdim] = fTracklet->GetDetector();
552 (*DebugStream()) << "trackIn"
553 <<"tracklet.=" << fTracklet
558 return NULL; // H->Projection(kEta, kPhi);
562 //________________________________________________________
563 TH1* AliTRDresolution::PlotTrackOut(const AliTRDtrackV1 *track)
565 // Store resolution/pulls of Kalman after last update with the TRD information
566 // at the radial position of the first tracklet. The following points are used
568 // - the (y,z,snp) of the first TRD tracklet
569 // - the (y, z, snp, tgl, pt) of the MC track reference
571 // Additionally the momentum resolution/pulls are calculated for usage in the
574 if(track) fkTrack = track;
578 //________________________________________________________
579 TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track)
582 // Plot MC distributions
586 AliDebug(2, "No MC defined. Results will not be available.");
589 if(track) fkTrack = track;
591 AliDebug(4, "No Track defined.");
594 Int_t bc(TMath::Abs(fkESD->GetTOFbc()));
598 AliWarning("No output container defined.");
601 // retriev track characteristics
602 Int_t pdg = fkMC->GetPDG(),
603 sIdx(AliTRDpidUtil::Pdg2Pid(TMath::Abs(pdg))+1), // species index
606 label(fkMC->GetLabel());
608 if(!fDBPDG) fDBPDG=TDatabasePDG::Instance();
609 TParticlePDG *ppdg(fDBPDG->GetParticle(pdg));
610 if(ppdg) sign = ppdg->Charge() > 0. ? 1 : -1;
613 AliTRDgeometry *geo(AliTRDinfoGen::Geometry());
614 AliTRDseedV1 *fTracklet(NULL); TObjArray *clInfoArr(NULL);
616 Double_t x, y, z, pt, dydx, dzdx, dzdl;
617 Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
618 Double_t covR[7]/*, cov[3]*/;
620 /* if(DebugLevel()>=3){
621 // get first detector
623 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
624 if(!(fTracklet = fkTrack->GetTracklet(ily))) continue;
625 det = fTracklet->GetDetector();
629 TVectorD X(12), Y(12), Z(12), dX(12), dY(12), dZ(12), vPt(12), dPt(12), budget(12), cCOV(12*15);
631 m = fkTrack->GetMass();
632 if(fkMC->PropagateKalman(&X, &Y, &Z, &dX, &dY, &dZ, &vPt, &dPt, &budget, &cCOV, m)){
633 (*DebugStream()) << "MCkalman"
650 AliTRDcluster *c(NULL);
651 Double_t val[kNdim+1];
652 for(Int_t ily=0; ily<AliTRDgeometry::kNlayer; ily++){
653 if(!(fTracklet = fkTrack->GetTracklet(ily)))/* ||
654 !fTracklet->IsOK())*/ continue;
656 x= x0 = fTracklet->GetX();
657 Bool_t rc(fTracklet->IsRowCross()); Float_t eta, phi;
658 if(!fkMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, eta, phi, s)) continue;
660 // MC track position at reference radial position
662 Float_t ymc = y0 - dx*dydx0;
663 Float_t zmc = z0 - dx*dzdx0;
669 val[kSpeciesChgRC]= rc?0.:sign*sIdx;
670 val[kPt] = pt0<0.8?0:1;//GetPtBin(pt0);
671 Double_t tilt(fTracklet->GetTilt());
673 // ,corr(1./(1. + t2))
674 // ,cost(TMath::Sqrt(corr));
676 AliExternalTrackParam *tin(fkTrack->GetTrackIn());
677 if(ily==0 && tin){ // trackIn residuals
678 // check radial position
679 if(TMath::Abs(tin->GetX()-x)>1.e-3) AliDebug(1, Form("TrackIn radial mismatch. dx[cm]=%+4.1f", tin->GetX()-x));
681 // Float_t phi = TMath::ATan2(y0, x0);
682 val[kBC] = (bc>=kNbunchCross)?(kNbunchCross-1):bc;
683 val[kYrez] = tin->GetY()-ymc;
684 val[kZrez] = tin->GetZ()-zmc;
685 val[kPrez] = (TMath::ASin(tin->GetSnp())-TMath::ATan(dydx0))*TMath::RadToDeg();
686 if((H = (THnSparseI*)fContainer->At(kMCtrackIn))) H->Fill(val);
689 if(bc>1) break; // do nothing for the rest of TRD objects if satellite bunch
692 dydx = fTracklet->GetYref(1);
693 dzdx = fTracklet->GetZref(1);
694 dzdl = fTracklet->GetTgl();
695 y = fTracklet->GetYref(0);
697 z = fTracklet->GetZref(0);
699 pt = TMath::Abs(fTracklet->GetPt());
700 fTracklet->GetCovRef(covR);
703 val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
705 val[kNdim] = 1.e2*(pt/pt0-1.);
706 if((H = (THnSparse*)fContainer->At(kMCtrack))) H->Fill(val);
707 /* // theta resolution/ tgl pulls
708 Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0),
709 dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0);
710 ((TH2I*)arr->At(6))->Fill(dzdl0, TMath::ATan(dtgl));
711 ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4]));
712 // pt resolution \\ 1/pt pulls \\ p resolution for PID
713 Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0,
714 p = TMath::Sqrt(1.+ dzdl*dzdl)*pt;
715 ((TH3S*)((TObjArray*)arr->At(8)))->Fill(pt0, pt/pt0-1., sign*sIdx);
716 ((TH3S*)((TObjArray*)arr->At(9)))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), sign*sIdx);
717 ((TH3S*)((TObjArray*)arr->At(10)))->Fill(p0, p/p0-1., sign*sIdx);*/
719 // Fill Debug stream for MC track
721 Int_t det(fTracklet->GetDetector());
722 (*DebugStream()) << "MC"
734 // Fill Debug stream for Kalman track
735 (*DebugStream()) << "MCtrack"
747 // tracklet residuals
748 dydx = fTracklet->GetYfit(1) + tilt*dzdx0;
749 dzdx = fTracklet->GetZfit(1);
750 y = fTracklet->GetYfit(0);
752 z = fTracklet->GetZfit(0);
754 val[kYrez] = dy - dz*tilt;
755 val[kPrez] = TMath::ATan((dydx - dydx0)/(1.+ dydx*dydx0))*TMath::RadToDeg();
756 val[kZrez] = dz + dy*tilt;
757 // val[kNdim] = pt/pt0-1.;
758 if((H = (THnSparse*)fContainer->At(kMCtracklet))) H->Fill(val);
761 // Fill Debug stream for tracklet
763 Float_t s2y = fTracklet->GetS2Y();
764 Float_t s2z = fTracklet->GetS2Z();
765 (*DebugStream()) << "MCtracklet"
776 AliTRDpadPlane *pp = geo->GetPadPlane(ily, AliTRDgeometry::GetStack(fTracklet->GetDetector()));
777 Float_t zr0 = pp->GetRow0() + pp->GetAnodeWireOffset();
778 //Double_t exb = AliTRDCommonParam::Instance()->GetOmegaTau(1.5);
780 H = (THnSparse*)fContainer->At(kMCcluster);
781 val[kPt] = TMath::ATan(dydx0)*TMath::RadToDeg();
782 //Float_t corr = 1./TMath::Sqrt(1.+dydx0*dydx0+dzdx0*dzdx0);
784 Float_t padCorr(tilt*fTracklet->GetPadLength());
785 fTracklet->ResetClusterIter(kTRUE);
786 while((c = fTracklet->NextCluster())){
787 if(row0<0) row0 = c->GetPadRow();
788 x = c->GetX();//+fXcorr[c->GetDetector()][c->GetLocalTimeBin()];
789 y = c->GetY() + padCorr*(c->GetPadRow() - row0);
796 val[kYrez] = dy - dz*tilt;
798 val[kZrez] = 0.; AliTRDcluster *cc(NULL); Int_t ic(0), tb(c->GetLocalTimeBin()); Float_t q(TMath::Abs(c->GetQ()));
799 if((cc = fTracklet->GetClusters(tb-1))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
800 if((cc = fTracklet->GetClusters(tb-2))) {val[kZrez] += TMath::Abs(cc->GetQ()); ic++;}
801 if(ic) val[kZrez] /= (ic*q);
805 // Fill calibration container
806 Float_t d = zr0 - zmc;
807 d -= ((Int_t)(2 * d)) / 2.0;
808 if (d > 0.25) d = 0.5 - d;
809 AliTRDclusterInfo *clInfo = new AliTRDclusterInfo;
810 clInfo->SetCluster(c);
811 clInfo->SetMC(pdg, label);
812 clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0);
813 clInfo->SetResolution(dy);
814 clInfo->SetAnisochronity(d);
815 clInfo->SetDriftLength(dx);
816 clInfo->SetTilt(tilt);
817 if(fMCcl) fMCcl->Add(clInfo);
818 else AliDebug(1, "MCcl exchange container missing. Activate by calling \"InitExchangeContainers()\"");
821 clInfoArr=new TObjArray(AliTRDseedV1::kNclusters);
822 clInfoArr->SetOwner(kFALSE);
824 clInfoArr->Add(clInfo);
828 if(DebugLevel()>=5 && clInfoArr){
829 (*DebugStream()) << "MCcluster"
830 <<"clInfo.=" << clInfoArr
835 if(clInfoArr) delete clInfoArr;
840 //__________________________________________________________________________
841 Int_t AliTRDresolution::GetPtBin(Float_t pt)
843 // Find pt bin according to local pt segmentation
845 while(ipt<AliTRDresolution::kNpt){
846 if(pt<fgPtBin[ipt+1]) break;
852 //________________________________________________________
853 Float_t AliTRDresolution::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt)
855 // return mean number of entries/bin of histogram "h"
856 // if option "opt" is given the following values are accepted:
857 // "<" : consider only entries less than "cut"
858 // ">" : consider only entries greater than "cut"
860 //Int_t dim(h->GetDimension());
861 Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ());
862 Double_t sum(0.); Int_t n(0);
863 for(Int_t ix(1); ix<=nbx; ix++)
864 for(Int_t iy(1); iy<=nby; iy++)
865 for(Int_t iz(1); iz<=nbz; iz++){
866 if(strcmp(opt, "")==0){sum += h->GetBinContent(ix, iy, iz); n++;}
868 if(strcmp(opt, "<")==0) {
869 if(h->GetBinContent(ix, iy, iz)<cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
870 } else if(strcmp(opt, ">")==0){
871 if(h->GetBinContent(ix, iy, iz)>cut) {sum += h->GetBinContent(ix, iy, iz); n++;}
872 } else {sum += h->GetBinContent(ix, iy, iz); n++;}
878 //________________________________________________________
879 Bool_t AliTRDresolution::GetRefFigure(Int_t ifig)
882 // Get the reference figures
886 AliWarning("Please provide a canvas to draw results.");
889 /* Int_t selection[100], n(0), selStart(0); //
890 Int_t ly0(0), dly(5);
891 TList *l(NULL); TVirtualPad *pad(NULL); */
896 AliWarning(Form("Reference plot [%d] missing result", ifig));
901 //________________________________________________________
902 void AliTRDresolution::MakePtSegmentation(Float_t pt0, Float_t dpt)
905 for(Int_t j(0); j<=kNpt; j++){
906 pt0+=(TMath::Exp(j*j*dpt)-1.);
911 //________________________________________________________
912 void AliTRDresolution::MakeSummary()
914 // Build summary plots
917 AliError("Missing results");
920 TVirtualPad *p(NULL); TCanvas *cOut(NULL);
921 TObjArray *arr(NULL); TH2 *h2(NULL);
923 // cluster resolution
925 gStyle->SetPalette(1);
926 const Int_t nClViews(11);
927 const Char_t *vClName[nClViews] = {"ClY", "ClYn", "ClYp", "ClQn", "ClQp", "ClYXTCp", "ClYXTCn", "ClYXPh", "ClYXPh", "ClY", "ClYn"};
928 const UChar_t vClOpt[nClViews] = {1, 1, 1, 0, 0, 0, 0, 0, 1, 0, 0};
929 const Int_t nTrkltViews(10);
930 const Char_t *vTrkltName[nTrkltViews] = {"TrkltY", "TrkltYn", "TrkltYp", "TrkltPhn", "TrkltPhp", "TrkltZ", "TrkltQn", "TrkltQp", "TrkltPn", "TrkltPp"};
931 const Int_t nTrkInViews(4);
932 const Char_t *vTrkInName[nTrkInViews][6] = {
933 {"TrkInY", "TrkInYn", "TrkInYp", "TrkInZ", "TrkInPhn", "TrkInPhp"},
934 {"TrkInYnl", "TrkInYni", "TrkInYnh", "TrkInYpl", "TrkInYpi", "TrkInYph"},
935 {"TrkInXnl", "TrkInZn", "TrkInXl", "TrkInXpl", "TrkInZp", "TrkInYh"},
936 {"TrkInPhnl", "TrkInPhni", "TrkInPhnh", "TrkInPhpl", "TrkInPhpi", "TrkInPhph"}};
937 const Int_t nTrkViews(10);
938 const Char_t *vTrkName[nTrkViews] = {"TrkY", "TrkYn", "TrkYp", "TrkPhn", "TrkPhp", "TrkZ", "TrkQn", "TrkQp", "TrkPn", "TrkPp"};
939 const Char_t *typName[] = {"", "MC"};
941 for(Int_t ityp(0); ityp<(HasMCdata()?2:1); ityp++){
942 if((arr = (TObjArray*)fProj->At(ityp?kMCcluster:kCluster))){
943 for(Int_t iview(0); iview<nClViews; iview++){
944 cOut = new TCanvas(Form("TRDsummary%s_%sCl%02d", GetName(), typName[ityp], iview), "Cluster Resolution", 1024, 768);
945 cOut->Divide(3,2, 1.e-5, 1.e-5);
947 for(Int_t iplot(0); iplot<6; iplot++){
948 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
949 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vClName[iview], iplot)))) continue;
951 if(vClOpt[iview]==0) h2->Draw("colz");
952 else if(vClOpt[iview]==1) DrawSigma(h2, 1.e4, 2.e2, 6.5e2, "#sigma(#Deltay) [#mum]");
954 if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
958 // tracklet systematic
959 if((arr = (TObjArray*)fProj->At(ityp?kMCtracklet:kTracklet))){
960 for(Int_t iview(0); iview<nTrkltViews; iview++){
961 cOut = new TCanvas(Form("TRDsummary%s_%sTrklt%02d", GetName(), typName[ityp], iview), "Tracklet Resolution", 1024, 768);
962 cOut->Divide(3,2, 1.e-5, 1.e-5);
964 for(Int_t iplot(0); iplot<6; iplot++){
965 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
966 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s%d_2D", typName[ityp], vTrkltName[iview], iplot)))) continue;
967 h2->Draw("colz"); nplot++;
969 if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
973 // trackIn systematic
974 if((arr = (TObjArray*)fProj->At(ityp?kMCtrackIn:kTrackIn))){
975 for(Int_t iview(0); iview<nTrkInViews; iview++){
976 cOut = new TCanvas(Form("TRDsummary%s_%sTrkIn%02d", GetName(), typName[ityp], iview), "Track IN Resolution", 1024, 768);
977 cOut->Divide(3,2, 1.e-5, 1.e-5);
979 for(Int_t iplot(0); iplot<6; iplot++){
980 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581);p->SetTopMargin(0.08262712);
981 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%s_2D", typName[ityp], vTrkInName[iview][iplot])))) continue;
982 h2->Draw("colz"); nplot++;
984 if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
989 // track MC systematic
990 if((arr = (TObjArray*)fProj->At(kMCtrack))) {
991 for(Int_t iview(0); iview<nTrkViews; iview++){
992 cOut = new TCanvas(Form("TRDsummary%s_MCTrk%02d", GetName(), iview), "Track Resolution", 1024, 768);
993 cOut->Divide(3,2, 1.e-5, 1.e-5);
995 for(Int_t iplot(0); iplot<6; iplot++){
996 p=cOut->cd(iplot+1); p->SetRightMargin(0.1572581); p->SetTopMargin(0.08262712);
997 if(!(h2 = (TH2*)arr->FindObject(Form("H%s%d_2D", vTrkName[iview], iplot)))) continue;
998 h2->Draw("colz"); nplot++;
1000 if(nplot) cOut->SaveAs(Form("%s.gif", cOut->GetName()));
1006 gStyle->SetPalette(1);
1009 //________________________________________________________
1010 void AliTRDresolution::DrawSigma(TH2 *h2, Float_t scale, Float_t m, Float_t M, const Char_t *title)
1012 // Draw error bars scaled with "scale" instead of content values
1013 //use range [m,M] if limits are specified
1016 TH2 *h2e = (TH2F*)h2->Clone(Form("%s_E", h2->GetName()));
1017 h2e->SetContour(10);
1018 if(M>m) h2e->GetZaxis()->SetRangeUser(m, M);
1019 if(title) h2e->GetZaxis()->SetTitle(title);
1021 for(Int_t ix(1); ix<=h2->GetNbinsX(); ix++){
1022 for(Int_t iy(1); iy<=h2->GetNbinsY(); iy++){
1023 if(h2->GetBinContent(ix, iy)<-100.) continue;
1024 Float_t v(scale*h2->GetBinError(ix, iy));
1025 if(M>m && v<m) v=m+TMath::Abs((M-m)*1.e-3);
1026 h2e->SetBinContent(ix, iy, v);
1032 //________________________________________________________
1033 void AliTRDresolution::GetRange(TH2 *h2, Char_t mod, Float_t *range)
1035 // Returns the range of the bulk of data in histogram h2. Removes outliers.
1036 // The "range" vector should be initialized with 2 elements
1037 // Option "mod" can be any of
1038 // - 0 : gaussian like distribution
1039 // - 1 : tailed distribution
1041 Int_t nx(h2->GetNbinsX())
1042 , ny(h2->GetNbinsY())
1044 Double_t *data=new Double_t[n];
1045 for(Int_t ix(1), in(0); ix<=nx; ix++){
1046 for(Int_t iy(1); iy<=ny; iy++)
1047 data[in++] = h2->GetBinContent(ix, iy);
1049 Double_t mean, sigm;
1050 AliMathBase::EvaluateUni(n, data, mean, sigm, Int_t(n*.8));
1052 range[0]=mean-3.*sigm; range[1]=mean+3.*sigm;
1053 if(mod==1) range[0]=TMath::Max(Float_t(1.e-3), range[0]);
1054 AliDebug(2, Form("h[%s] range0[%f %f]", h2->GetName(), range[0], range[1]));
1055 TH1S h1("h1SF0", "", 100, range[0], range[1]);
1060 case 0:// gaussian distribution
1062 TF1 fg("fg", "gaus", mean-3.*sigm, mean+3.*sigm);
1064 mean=fg.GetParameter(1); sigm=fg.GetParameter(2);
1065 range[0] = mean-2.5*sigm;range[1] = mean+2.5*sigm;
1066 AliDebug(2, Form(" rangeG[%f %f]", range[0], range[1]));
1069 case 1:// tailed distribution
1071 Int_t bmax(h1.GetMaximumBin());
1072 Int_t jBinMin(1), jBinMax(100);
1073 for(Int_t ibin(bmax); ibin--;){
1074 if(h1.GetBinContent(ibin)<1.){
1075 jBinMin=ibin; break;
1078 for(Int_t ibin(bmax); ibin++;){
1079 if(h1.GetBinContent(ibin)<1.){
1080 jBinMax=ibin; break;
1083 range[0]=h1.GetBinCenter(jBinMin); range[1]=h1.GetBinCenter(jBinMax);
1084 AliDebug(2, Form(" rangeT[%f %f]", range[0], range[1]));
1093 //________________________________________________________
1094 Bool_t AliTRDresolution::MakeProjectionCluster(Bool_t mc)
1097 const Int_t kNcontours(9);
1098 const Int_t kNstat(300);
1099 Int_t cidx=mc?kMCcluster:kCluster;
1100 if(fProj && fProj->At(cidx)) return kTRUE;
1102 AliError("Missing data container.");
1106 if(!(H = (THnSparse*)fContainer->At(cidx))){
1107 AliError(Form("Missing/Wrong data @ %d.", cidx));
1110 Int_t ndim(H->GetNdimensions());
1111 Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1112 TAxis *aa[kNdim], *as(NULL), *apt(NULL); memset(aa, 0, sizeof(TAxis*) * kNdim);
1113 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1114 if(ndim > kPt) apt = H->GetAxis(kPt);
1115 if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1116 // build list of projections
1117 const Int_t nsel(12), npsel(5);
1118 // define rebinning strategy
1119 const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1120 AliTRDresolutionProjection hp[kClNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1121 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1122 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1123 isel++; // new selection
1124 hp[ih].Build(Form("H%sClY%d", mc?"MC":"", ily), Form("Clusters :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1125 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1126 php[isel][np[isel]++] = &hp[ih++];
1127 hp[ih].Build(Form("H%sClYn%d", mc?"MC":"", ily), Form("Clusters[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1128 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1129 php[isel][np[isel]++] = &hp[ih++];
1130 hp[ih].Build(Form("H%sClQn%d", mc?"MC":"", ily), Form("Clusters[-]:: Charge distribution ly%d", ily), kEta, kPhi, kSpeciesChgRC, aa);
1131 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1132 hp[ih].SetShowRange(20., 40.);
1133 php[isel][np[isel]++] = &hp[ih++];
1134 hp[ih].Build(Form("H%sClYXTCn%d", mc?"MC":"", ily), Form("Clusters[-]:: r-#phi(x,TC) residuals ly%d", ily), kPrez, kZrez, kYrez, aa);
1135 // hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1136 php[isel][np[isel]++] = &hp[ih++];
1137 hp[ih].Build(Form("H%sClYXPh%d", mc?"MC":"", ily), Form("Clusters :: r-#phi(x,#Phi) residuals ly%d", ily), kPrez, kPt, kYrez, aa);
1138 // hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1139 php[isel][np[isel]++] = &hp[ih++];
1140 isel++; // new selection
1141 php[isel][np[isel]++] = &hp[ih-5]; // relink HClY
1142 hp[ih].Build(Form("H%sClYp%d", mc?"MC":"", ily), Form("Clusters[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1143 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1144 php[isel][np[isel]++] = &hp[ih++];
1145 hp[ih].Build(Form("H%sClQp%d", mc?"MC":"", ily), Form("Clusters[+]:: Charge distribution ly%d", ily), kEta, kPhi, kSpeciesChgRC, aa);
1146 hp[ih].SetShowRange(20., 40.);
1147 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1148 php[isel][np[isel]++] = &hp[ih++];
1149 hp[ih].Build(Form("H%sClYXTCp%d", mc?"MC":"", ily), Form("Clusters[+]:: r-#phi(x,TC) residuals ly%d", ily), kPrez, kZrez, kYrez, aa);
1150 // hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1151 php[isel][np[isel]++] = &hp[ih++];
1152 php[isel][np[isel]++] = &hp[ih-4]; // relink HClYXPh
1155 Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1), chBin(apt?apt->FindBin(0.):-1);
1156 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1157 v = H->GetBinContent(ib, coord); if(v<1.) continue;
1160 if(rcBin>0 && coord[kSpeciesChgRC] == rcBin) continue;
1163 ch = 0; // [-] track
1164 if(chBin>0 && coord[kPt] > chBin) ch = 1; // [+] track
1167 for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v);
1169 TObjArray *arr(NULL);
1170 fProj->AddAt(arr = new TObjArray(kClNproj), cidx);
1174 if(!hp[ih].fH) continue;
1175 Int_t mid(1), nstat(kNstat);
1176 if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=300;}
1177 if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1184 //________________________________________________________
1185 Bool_t AliTRDresolution::MakeProjectionTracklet(Bool_t mc)
1188 const Int_t kNcontours(9);
1189 const Int_t kNstat(100);
1190 Int_t cidx=mc?kMCtracklet:kTracklet;
1191 if(fProj && fProj->At(cidx)) return kTRUE;
1193 AliError("Missing data container.");
1197 if(!(H = (THnSparse*)fContainer->At(cidx))){
1198 AliError(Form("Missing/Wrong data @ %d.", cidx));
1201 Int_t ndim(H->GetNdimensions());
1202 Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
1203 TAxis *aa[kNdim+1], *as(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
1204 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1205 if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1206 // build list of projections
1207 const Int_t nsel(18), npsel(6);
1208 // define rebinning strategy
1209 const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1210 AliTRDresolutionProjection hp[kTrkltNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1211 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1212 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1213 isel++; // new selection
1214 hp[ih].Build(Form("H%sTrkltY%d", mc?"MC":"", ily), Form("Tracklets :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1215 hp[ih].SetShowRange(-0.03,0.03);
1216 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1217 php[isel][np[isel]++] = &hp[ih++];
1218 hp[ih].Build(Form("H%sTrkltYn%d", mc?"MC":"", ily), Form("Tracklets[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1219 hp[ih].SetShowRange(-0.03,0.03);
1220 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1221 php[isel][np[isel]++] = &hp[ih++];
1222 hp[ih].Build(Form("H%sTrkltPhn%d", mc?"MC":"", ily), Form("Tracklets[-]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1223 hp[ih].SetShowRange(-0.5,0.5);
1224 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1225 php[isel][np[isel]++] = &hp[ih++];
1226 hp[ih].Build(Form("H%sTrkltPn%d", mc?"MC":"", ily), Form("Tracklets[-]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1227 hp[ih].SetShowRange(6.,12.);
1228 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1229 php[isel][np[isel]++] = &hp[ih++];
1230 hp[ih].Build(Form("H%sTrkltYPn%d", mc?"MC":"", ily), Form("Tracklets[-]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1231 php[isel][np[isel]++] = &hp[ih++];
1232 hp[ih].Build(Form("H%sTrkltQn%d", mc?"MC":"", ily), Form("Tracklets[-]:: dQdl ly%d", ily), kEta, kPhi, kNdim, aa);
1233 hp[ih].SetShowRange(700.,1100.);
1234 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1235 php[isel][np[isel]++] = &hp[ih++];
1236 isel++; // new selection
1237 php[isel][np[isel]++] = &hp[ih-6]; // relink first histo
1238 hp[ih].Build(Form("H%sTrkltYp%d", mc?"MC":"", ily), Form("Tracklets[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1239 hp[ih].SetShowRange(-0.03,0.03);
1240 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1241 php[isel][np[isel]++] = &hp[ih++];
1242 hp[ih].Build(Form("H%sTrkltPhp%d", mc?"MC":"", ily), Form("Tracklets[+]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1243 hp[ih].SetShowRange(-0.5,0.5);
1244 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1245 php[isel][np[isel]++] = &hp[ih++];
1246 hp[ih].Build(Form("H%sTrkltPp%d", mc?"MC":"", ily), Form("Tracklets[+]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1247 hp[ih].SetShowRange(6.,12.);
1248 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1249 php[isel][np[isel]++] = &hp[ih++];
1250 hp[ih].Build(Form("H%sTrkltYPp%d", mc?"MC":"", ily), Form("Tracklets[+]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1251 php[isel][np[isel]++] = &hp[ih++];
1252 hp[ih].Build(Form("H%sTrkltQp%d", mc?"MC":"", ily), Form("Tracklets[+]:: dQdl ly%d", ily), kEta, kPhi, kNdim, aa);
1253 hp[ih].SetShowRange(700.,1100.);
1254 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1255 php[isel][np[isel]++] = &hp[ih++];
1256 isel++; // new selection
1257 hp[ih].Build(Form("H%sTrkltZ%d", mc?"MC":"", ily), Form("Tracklets[RC]:: z residuals ly%d", ily), kEta, kPhi, kZrez, aa);
1258 hp[ih].SetShowRange(-0.1,0.1);
1259 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1260 php[isel][np[isel]++] = &hp[ih++];
1263 Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1);
1264 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1265 v = H->GetBinContent(ib, coord);
1267 ly = coord[kBC]-1; // layer selection
1269 ch = 0; // [-] track
1270 if(rcBin>0){ // debug mode in which species are also saved
1271 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
1272 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
1275 for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v);
1277 TObjArray *arr(NULL);
1278 fProj->AddAt(arr = new TObjArray(kTrkltNproj), cidx);
1282 if(!hp[ih].fH) continue;
1283 Int_t mid(0), nstat(kNstat);
1284 if(strchr(hp[ih].fH->GetName(), 'Q')){ mid=2; nstat=200;}
1285 if(!(h2 = hp[ih].Projection2D(nstat, kNcontours, mid))) continue;
1291 //________________________________________________________
1292 Bool_t AliTRDresolution::MakeProjectionTrackIn(Bool_t mc)
1296 const Int_t kNcontours(9);
1297 const Int_t kNstat(30);
1298 Int_t cidx=mc?kMCtrackIn:kTrackIn;
1299 if(fProj && fProj->At(cidx)) return kTRUE;
1301 AliError("Missing data container.");
1305 if(!(H = (THnSparse*)fContainer->At(cidx))){
1306 AliError(Form("Missing/Wrong data @ %d.", Int_t(cidx)));
1310 Int_t coord[kNdim]; memset(coord, 0, sizeof(Int_t) * kNdim); Double_t v = 0.;
1311 Int_t ndim(H->GetNdimensions());
1312 TAxis *aa[kNdim+1], *as(NULL), *ap(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
1313 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1314 if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1315 if(ndim > kPt) ap = H->GetAxis(kPt);
1316 // build list of projections
1317 const Int_t nsel(15), npsel(3);
1318 // define rebinning strategy
1319 const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1320 AliTRDresolutionProjection hp[kMCTrkInNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1321 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1322 // define list of projections
1323 isel++; // negative low pt tracks
1324 hp[ih].Build(Form("H%sTrkInYnl", mc?"MC":""), "TrackIn[-]:: #Deltay{p_{t}[GeV/c]<0.8}", kEta, kPhi, kYrez, aa);
1325 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1326 php[isel][np[isel]++] = &hp[ih++];
1327 hp[ih].Build(Form("H%sTrkInPhnl", mc?"MC":""), "TrackIn[-]:: #Delta#phi{p_{t}[GeV/c]<0.8}", kEta, kPhi, kPrez, aa);
1328 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1329 php[isel][np[isel]++] = &hp[ih++];
1330 hp[ih].Build(Form("H%sTrkInXnl", mc?"MC":""), "TrackIn[-]:: #Deltax{p_{t}[GeV/c]<0.8}", kEta, kPhi, kNdim+1, aa);
1331 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1332 php[isel][np[isel]++] = &hp[ih++];
1333 isel++; // negative intermediate pt tracks
1334 hp[ih].Build(Form("H%sTrkInYni", mc?"MC":""), "TrackIn[-]:: #Deltay{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kYrez, aa);
1335 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1336 php[isel][np[isel]++] = &hp[ih++];
1337 hp[ih].Build(Form("H%sTrkInPhni", mc?"MC":""), "TrackIn[-]:: #Delta#phi{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kPrez, aa);
1338 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1339 php[isel][np[isel]++] = &hp[ih++];
1340 hp[ih].Build(Form("H%sTrkInXni", mc?"MC":""), "TrackIn[-]:: #Deltax{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kNdim+1, aa);
1341 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1342 php[isel][np[isel]++] = &hp[ih++];
1343 isel++; // negative high pt tracks
1344 hp[ih].Build(Form("H%sTrkInYnh", mc?"MC":""), "TrackIn[-]:: #Deltay{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kYrez, aa);
1345 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1346 php[isel][np[isel]++] = &hp[ih++];
1347 hp[ih].Build(Form("H%sTrkInPhnh", mc?"MC":""), "TrackIn[-]:: #Delta#phi{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kPrez, aa);
1348 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1349 php[isel][np[isel]++] = &hp[ih++];
1350 hp[ih].Build(Form("H%sTrkInXnh", mc?"MC":""), "TrackIn[-]:: #Deltax{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kNdim+1, aa);
1351 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1352 php[isel][np[isel]++] = &hp[ih++];
1353 isel++; // positive low pt tracks
1354 hp[ih].Build(Form("H%sTrkInYpl", mc?"MC":""), "TrackIn[+]:: #Deltay{p_{t}[GeV/c]<0.8}", kEta, kPhi, kYrez, aa);
1355 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1356 php[isel][np[isel]++] = &hp[ih++];
1357 hp[ih].Build(Form("H%sTrkInPhpl", mc?"MC":""), "TrackIn[+]:: #Delta#phi{p_{t}[GeV/c]<0.8}", kEta, kPhi, kPrez, aa);
1358 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1359 php[isel][np[isel]++] = &hp[ih++];
1360 hp[ih].Build(Form("H%sTrkInXpl", mc?"MC":""), "TrackIn[+]:: #Deltax{p_{t}[GeV/c]<0.8}", kEta, kPhi, kNdim+1, aa);
1361 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1362 php[isel][np[isel]++] = &hp[ih++];
1363 isel++; // positive intermediate pt tracks
1364 hp[ih].Build(Form("H%sTrkInYpi", mc?"MC":""), "TrackIn[+]:: #Deltay{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kYrez, aa);
1365 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1366 php[isel][np[isel]++] = &hp[ih++];
1367 hp[ih].Build(Form("H%sTrkInPhpi", mc?"MC":""), "TrackIn[+]:: #Delta#phi{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kPrez, aa);
1368 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1369 php[isel][np[isel]++] = &hp[ih++];
1370 hp[ih].Build(Form("H%sTrkInXpi", mc?"MC":""), "TrackIn[+]:: #Deltax{0.8<=p_{t}[GeV/c]<1.5}", kEta, kPhi, kNdim+1, aa);
1371 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1372 php[isel][np[isel]++] = &hp[ih++];
1373 isel++; // positive high pt tracks
1374 hp[ih].Build(Form("H%sTrkInYph", mc?"MC":""), "TrackIn[+]:: #Deltay{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kYrez, aa);
1375 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1376 php[isel][np[isel]++] = &hp[ih++];
1377 hp[ih].Build(Form("H%sTrkInPhph", mc?"MC":""), "TrackIn[+]:: #Delta#phi{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kPrez, aa);
1378 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1379 php[isel][np[isel]++] = &hp[ih++];
1380 hp[ih].Build(Form("H%sTrkInXph", mc?"MC":""), "TrackIn[+]:: #Deltax{p_{t}[GeV/c]>=1.5}", kEta, kPhi, kNdim+1, aa);
1381 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1382 php[isel][np[isel]++] = &hp[ih++];
1383 isel++; // negative RC tracks
1384 hp[ih].Build(Form("H%sTrkInZn", mc?"MC":""), "TrackIn[RC-]:: #Deltaz", kEta, kPhi, kZrez, aa);
1385 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1386 php[isel][np[isel]++] = &hp[ih++];
1387 isel++; // positive RC tracks
1388 hp[ih].Build(Form("H%sTrkInZp", mc?"MC":""), "TrackIn[RC+]:: #Deltaz", kEta, kPhi, kZrez, aa);
1389 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1390 php[isel][np[isel]++] = &hp[ih++];
1392 for(Int_t is(0); is<AliPID::kSPECIES; is++){
1393 isel++; // negative MC tracks
1394 hp[ih].Build(Form("HMCTrkInYn%s", AliPID::ParticleShortName(is)), Form("TrackIn[%s-]:: #Deltay", AliPID::ParticleShortName(is)), kEta, kPhi, kYrez, aa);
1395 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1396 php[isel][np[isel]++] = &hp[ih++];
1397 hp[ih].Build(Form("HMCTrkInPhn%s", AliPID::ParticleShortName(is)), Form("TrackIn[%s-]:: #Delta#phi", AliPID::ParticleShortName(is)), kEta, kPhi, kPrez, aa);
1398 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1399 php[isel][np[isel]++] = &hp[ih++];
1401 for(Int_t is(0); is<AliPID::kSPECIES; is++){
1402 isel++; // positive MC tracks
1403 hp[ih].Build(Form("HMCTrkInYp%s", AliPID::ParticleShortName(is)), Form("TrackIn[%s+]:: #Deltay", AliPID::ParticleShortName(is)), kEta, kPhi, kYrez, aa);
1404 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1405 php[isel][np[isel]++] = &hp[ih++];
1406 hp[ih].Build(Form("HMCTrkInPhp%s", AliPID::ParticleShortName(is)), Form("TrackIn[%s+]:: #Delta#phi", AliPID::ParticleShortName(is)), kEta, kPhi, kPrez, aa);
1407 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1408 php[isel][np[isel]++] = &hp[ih++];
1411 // for(Int_t jsel(0); jsel<=isel; jsel++){
1412 // printf("Selection [%2d]\n", jsel);
1413 // for(Int_t n(0); n<np[jsel]; n++){
1414 // printf(" %d %s\n", n, php[jsel][n]->fH->GetName());
1419 Int_t ch(0), pt(0), rcBin(as?as->FindBin(0.):-1);
1420 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1421 v = H->GetBinContent(ib, coord);
1423 if(coord[kBC]>1) continue; // bunch cross cut
1425 ch = 0; // [-] track
1426 if(rcBin>0){ // debug mode in which species are also saved
1427 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
1428 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
1432 if(ap) pt = coord[kPt]-1;
1434 Int_t selection = ch*3+pt;
1435 for(Int_t jh(0); jh<np[selection]; jh++) php[selection][jh]->Increment(coord, v);
1436 if(!mc || rcBin<0 || ch==2) continue;
1437 Int_t spec = Int_t(TMath::Abs(as->GetBinCenter(coord[kSpeciesChgRC])))-1;
1438 selection+=(ch*AliPID::kSPECIES+spec);
1439 for(Int_t jh(0); jh<np[selection]; jh++) php[selection][jh]->Increment(coord, v);
1441 TObjArray *arr(NULL);
1442 fProj->AddAt(arr = new TObjArray(mc?kMCTrkInNproj:kTrkInNproj), cidx);
1446 if(!hp[ih].fH) continue;
1447 if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
1450 // build combined performance plots
1451 /*!dy negative tracks all momenta*/
1453 hp[iproj]+=hp[npsel+iproj]; hp[iproj]+=hp[npsel*2+iproj]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInYn", mc?"MC":""), "TrackIn[-]:: #Deltay");
1454 if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, arr->GetEntries());
1455 /*!dy positive tracks all momenta*/
1457 hp[iproj]+=hp[npsel+iproj]; hp[iproj]+=hp[npsel*2+iproj]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInYp", mc?"MC":""), "TrackIn[+]:: #Deltay");
1458 if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, arr->GetEntries());
1459 /*!dy all tracks all momenta*/
1460 hp[0]+=hp[9];hp[0].fH->SetNameTitle(Form("H%sTrkInY", mc?"MC":""), "TrackIn :: #Deltay");
1461 if((h2 = hp[0].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, arr->GetEntries());
1462 /*!dy all tracks high momenta*/
1464 hp[iproj]+=hp[iproj+9];hp[iproj].fH->SetNameTitle(Form("H%sTrkInYh", mc?"MC":""), "TrackIn :: #Deltay{p_{t}[GeV/c]>=1.5}");
1465 if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, arr->GetEntries());
1466 /*!dx all tracks low momenta*/
1469 hp[iproj]+=hp[iproj+9];hp[iproj].fH->SetNameTitle(Form("H%sTrkInXl", mc?"MC":""), "TrackIn :: #Deltax{p_{t}[GeV/c]<0.8}");
1470 if((h2 = hp[iproj].Projection2D(kNstat, kNcontours, 1))) arr->AddAt(h2, arr->GetEntries());
1472 /*!dphi negative tracks all momenta*/
1474 hp[iproj]+=hp[npsel+iproj]; hp[iproj]+=hp[npsel*2+iproj]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInPhn", mc?"MC":""), "TrackIn[-]:: #Delta#phi");
1475 if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, arr->GetEntries());
1476 /*!dphi positive tracks all momenta*/
1478 hp[iproj]+=hp[npsel+iproj]; hp[iproj]+=hp[npsel*2+iproj]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInPhp", mc?"MC":""), "TrackIn[+]:: #Delta#phi");
1479 if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, arr->GetEntries());
1480 /*!dy[RC] tracks all charges*/
1482 hp[iproj]+=hp[iproj+1]; hp[iproj].fH->SetNameTitle(Form("H%sTrkInZ", mc?"MC":""), "TrackIn[RC]:: #Deltaz");
1483 if((h2 = hp[iproj].Projection2D(kNstat, kNcontours))) arr->AddAt(h2, arr->GetEntries());
1489 //________________________________________________________
1490 Bool_t AliTRDresolution::MakeProjectionTrack()
1493 const Int_t kNcontours(9);
1494 const Int_t kNstat(100);
1495 Int_t cidx(kMCtrack);
1496 if(fProj && fProj->At(cidx)) return kTRUE;
1498 AliError("Missing data container.");
1502 if(!(H = (THnSparse*)fContainer->At(cidx))){
1503 AliError(Form("Missing/Wrong data @ %d.", cidx));
1506 Int_t ndim(H->GetNdimensions());
1507 Int_t coord[kNdim+1]; memset(coord, 0, sizeof(Int_t) * (kNdim+1)); Double_t v = 0.;
1508 TAxis *aa[kNdim+1], *as(NULL); memset(aa, 0, sizeof(TAxis*) * (kNdim+1));
1509 for(Int_t id(0); id<ndim; id++) aa[id] = H->GetAxis(id);
1510 if(ndim > kSpeciesChgRC) as = H->GetAxis(kSpeciesChgRC);
1511 // build list of projections
1512 const Int_t nsel(18), npsel(6);
1513 // define rebinning strategy
1514 const Int_t nEtaPhi(4); Int_t rebinEtaPhiX[nEtaPhi] = {1, 2, 5, 1}, rebinEtaPhiY[nEtaPhi] = {2, 1, 1, 5};
1515 AliTRDresolutionProjection hp[kTrkNproj], *php[nsel][npsel]; memset(php, 0, nsel*npsel*sizeof(AliTRDresolutionProjection*));
1516 Int_t ih(0), isel(-1), np[nsel]; memset(np, 0, nsel*sizeof(Int_t));
1517 for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
1518 isel++; // new selection
1519 hp[ih].Build(Form("HTrkY%d", ily), Form("Tracks :: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1520 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1521 php[isel][np[isel]++] = &hp[ih++];
1522 hp[ih].Build(Form("HTrkYn%d", ily), Form("Tracks[-]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1523 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1524 php[isel][np[isel]++] = &hp[ih++];
1525 hp[ih].Build(Form("HTrkPhn%d", ily), Form("Tracks[-]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1526 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1527 php[isel][np[isel]++] = &hp[ih++];
1528 hp[ih].Build(Form("HTrkPn%d", ily), Form("Tracks[-]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1529 hp[ih].SetShowRange(6.,12.);
1530 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1531 php[isel][np[isel]++] = &hp[ih++];
1532 hp[ih].Build(Form("HTrkYPn%d", ily), Form("Tracks[-]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1533 php[isel][np[isel]++] = &hp[ih++];
1534 hp[ih].Build(Form("HTrkQn%d", ily), Form("Tracks[-]:: #Deltap_{t}/p_{t} ly%d", ily), kEta, kPhi, kNdim, aa);
1535 //hp[ih].SetShowRange(700.,1100.);
1536 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1537 php[isel][np[isel]++] = &hp[ih++];
1538 isel++; // new selection
1539 php[isel][np[isel]++] = &hp[ih-6]; // relink first histo
1540 hp[ih].Build(Form("HTrkYp%d", ily), Form("Tracks[+]:: r-#phi residuals ly%d", ily), kEta, kPhi, kYrez, aa);
1541 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1542 php[isel][np[isel]++] = &hp[ih++];
1543 hp[ih].Build(Form("HTrkPhp%d", ily), Form("Tracks[+]:: #Delta#phi residuals ly%d", ily), kEta, kPhi, kPrez, aa);
1544 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1545 php[isel][np[isel]++] = &hp[ih++];
1546 hp[ih].Build(Form("HTrkPp%d", ily), Form("Tracks[+]:: Momentum distribution ly%d", ily), kEta, kPhi, kPt, aa);
1547 hp[ih].SetShowRange(6.,12.);
1548 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1549 php[isel][np[isel]++] = &hp[ih++];
1550 hp[ih].Build(Form("HTrkYPp%d", ily), Form("Tracks[+]:: r-#phi/p_{t} residuals ly%d", ily), kPt, kPhi, kYrez, aa);
1551 php[isel][np[isel]++] = &hp[ih++];
1552 hp[ih].Build(Form("HTrkQp%d", ily), Form("Tracks[+]:: #Deltap_{t}/p_{t} ly%d", ily), kEta, kPhi, kNdim, aa);
1553 //hp[ih].SetShowRange(700.,1100.);
1554 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1555 php[isel][np[isel]++] = &hp[ih++];
1556 isel++; // new selection
1557 hp[ih].Build(Form("HTrkZ%d", ily), Form("Tracks[RC]:: z residuals ly%d", ily), kEta, kPhi, kZrez, aa);
1558 hp[ih].SetRebinStrategy(nEtaPhi, rebinEtaPhiX, rebinEtaPhiY);
1559 php[isel][np[isel]++] = &hp[ih++];
1562 Int_t ly(0), ch(0), rcBin(as?as->FindBin(0.):-1);
1563 for (Long64_t ib(0); ib < H->GetNbins(); ib++) {
1564 v = H->GetBinContent(ib, coord);
1566 ly = coord[kBC]-1; // layer selection
1568 ch = 0; // [-] track
1569 if(rcBin>0){ // debug mode in which species are also saved
1570 if(coord[kSpeciesChgRC] > rcBin) ch = 1; // [+] track
1571 else if(coord[kSpeciesChgRC] == rcBin) ch = 2; // [RC] track
1574 for(Int_t jh(0); jh<np[isel]; jh++) php[isel][jh]->Increment(coord, v);
1576 TObjArray *arr(NULL);
1577 fProj->AddAt(arr = new TObjArray(kTrkNproj), cidx);
1581 if(!hp[ih].fH) continue;
1582 if(!(h2 = hp[ih].Projection2D(kNstat, kNcontours))) continue;
1588 //________________________________________________________
1589 Bool_t AliTRDresolution::PostProcess()
1591 // Fit, Project, Combine, Extract values from the containers filled during execution
1594 AliError("ERROR: list not available");
1598 AliInfo("Building array of projections ...");
1599 fProj = new TObjArray(kNclasses); fProj->SetOwner(kTRUE);
1602 //PROCESS EXPERIMENTAL DISTRIBUTIONS
1603 // Clusters residuals
1604 if(!MakeProjectionCluster()) return kFALSE;
1606 // Tracklet residual/pulls
1607 if(!MakeProjectionTracklet()) return kFALSE;
1609 // TRDin residual/pulls
1610 if(!MakeProjectionTrackIn()) return kFALSE;
1613 if(!HasMCdata()) return kTRUE;
1614 //PROCESS MC RESIDUAL DISTRIBUTIONS
1616 // CLUSTER Y RESOLUTION/PULLS
1617 if(!MakeProjectionCluster(kTRUE)) return kFALSE;
1620 // TRACKLET RESOLUTION/PULLS
1621 if(!MakeProjectionTracklet(kTRUE)) return kFALSE;
1624 // TRACK RESOLUTION/PULLS
1625 if(!MakeProjectionTrack()) return kFALSE;
1628 // TRACK TRDin RESOLUTION/PULLS
1629 if(!MakeProjectionTrackIn(kTRUE)) return kFALSE;
1636 //________________________________________________________
1637 void AliTRDresolution::Terminate(Option_t *opt)
1639 AliTRDrecoTask::Terminate(opt);
1640 if(HasPostProcess()) PostProcess();
1643 //________________________________________________________
1644 void AliTRDresolution::AdjustF1(TH1 *h, TF1 *f)
1646 // Helper function to avoid duplication of code
1647 // Make first guesses on the fit parameters
1649 // find the intial parameters of the fit !! (thanks George)
1650 Int_t nbinsy = Int_t(.5*h->GetNbinsX());
1652 for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;
1653 f->SetParLimits(0, 0., 3.*sum);
1654 f->SetParameter(0, .9*sum);
1655 Double_t rms = h->GetRMS();
1656 f->SetParLimits(1, -rms, rms);
1657 f->SetParameter(1, h->GetMean());
1659 f->SetParLimits(2, 0., 2.*rms);
1660 f->SetParameter(2, rms);
1661 if(f->GetNpar() <= 4) return;
1663 f->SetParLimits(3, 0., sum);
1664 f->SetParameter(3, .1*sum);
1666 f->SetParLimits(4, -.3, .3);
1667 f->SetParameter(4, 0.);
1669 f->SetParLimits(5, 0., 1.e2);
1670 f->SetParameter(5, 2.e-1);
1673 //________________________________________________________
1674 TObjArray* AliTRDresolution::BuildMonitorContainerCluster(const char* name, Bool_t expand, Float_t range)
1676 // Build performance histograms for AliTRDcluster.vs TRD track or MC
1677 // - y reziduals/pulls
1679 TObjArray *arr = new TObjArray(2);
1680 arr->SetName(name); arr->SetOwner();
1681 TH1 *h(NULL); char hname[100], htitle[300];
1683 // tracklet resolution/pull in y direction
1684 snprintf(hname, 100, "%s_%s_Y", GetNameId(), name);
1685 snprintf(htitle, 300, "Y res for \"%s\" @ %s;tg(#phi);#Delta y [cm];%s", GetNameId(), name, "Detector");
1686 Float_t rr = range<0.?fDyRange:range;
1687 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1689 if(expand) nybins*=2;
1690 h = new TH3S(hname, htitle,
1691 48, -.48, .48, // phi
1693 nybins, -0.5, nybins-0.5);// segment
1696 snprintf(hname, 100, "%s_%s_YZpull", GetNameId(), name);
1697 snprintf(htitle, 300, "YZ pull for \"%s\" @ %s;%s;#Delta y / #sigma_{y};#Delta z / #sigma_{z}", GetNameId(), name, "Detector");
1698 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1699 h = new TH3S(hname, htitle, 540, -0.5, 540-0.5, 100, -4.5, 4.5, 100, -4.5, 4.5);
1706 //________________________________________________________
1707 TObjArray* AliTRDresolution::BuildMonitorContainerTracklet(const char* name, Bool_t expand)
1709 // Build performance histograms for AliExternalTrackParam.vs TRD tracklet
1710 // - y reziduals/pulls
1711 // - z reziduals/pulls
1713 TObjArray *arr = BuildMonitorContainerCluster(name, expand, 0.05);
1715 TH1 *h(NULL); char hname[100], htitle[300];
1717 // tracklet resolution/pull in z direction
1718 snprintf(hname, 100, "%s_%s_Z", GetNameId(), name);
1719 snprintf(htitle, 300, "Z res for \"%s\" @ %s;tg(#theta);#Delta z [cm]", GetNameId(), name);
1720 if(!(h = (TH2S*)gROOT->FindObject(hname))){
1721 h = new TH2S(hname, htitle, 50, -1., 1., 100, -.05, .05);
1724 snprintf(hname, 100, "%s_%s_Zpull", GetNameId(), name);
1725 snprintf(htitle, 300, "Z pull for \"%s\" @ %s;tg(#theta);#Delta z / #sigma_{z};row cross", GetNameId(), name);
1726 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1727 h = new TH3S(hname, htitle, 50, -1., 1., 100, -5.5, 5.5, 2, -0.5, 1.5);
1728 h->GetZaxis()->SetBinLabel(1, "no RC");
1729 h->GetZaxis()->SetBinLabel(2, "RC");
1733 // tracklet to track phi resolution
1734 snprintf(hname, 100, "%s_%s_PHI", GetNameId(), name);
1735 snprintf(htitle, 300, "#Phi res for \"%s\" @ %s;tg(#phi);#Delta #phi [rad];%s", GetNameId(), name, "Detector");
1737 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1738 h = new TH3S(hname, htitle, 48, -.48, .48, 100, -.5, .5, nsgms, -0.5, nsgms-0.5);
1745 //________________________________________________________
1746 TObjArray* AliTRDresolution::BuildMonitorContainerTrack(const char* name)
1748 // Build performance histograms for AliExternalTrackParam.vs MC
1749 // - y resolution/pulls
1750 // - z resolution/pulls
1751 // - phi resolution, snp pulls
1752 // - theta resolution, tgl pulls
1753 // - pt resolution, 1/pt pulls, p resolution
1755 TObjArray *arr = BuildMonitorContainerTracklet(name);
1757 TH1 *h(NULL); char hname[100], htitle[300];
1761 snprintf(hname, 100, "%s_%s_SNPpull", GetNameId(), name);
1762 snprintf(htitle, 300, "SNP pull for \"%s\" @ %s;tg(#phi);#Delta snp / #sigma_{snp};entries", GetNameId(), name);
1763 if(!(h = (TH2I*)gROOT->FindObject(hname))){
1764 h = new TH2I(hname, htitle, 60, -.3, .3, 100, -4.5, 4.5);
1769 snprintf(hname, 100, "%s_%s_THT", GetNameId(), name);
1770 snprintf(htitle, 300, "#Theta res for \"%s\" @ %s;tg(#theta);#Delta #theta [rad];entries", GetNameId(), name);
1771 if(!(h = (TH2I*)gROOT->FindObject(hname))){
1772 h = new TH2I(hname, htitle, 100, -1., 1., 100, -5e-3, 5e-3);
1776 snprintf(hname, 100, "%s_%s_TGLpull", GetNameId(), name);
1777 snprintf(htitle, 300, "TGL pull for \"%s\" @ %s;tg(#theta);#Delta tgl / #sigma_{tgl};entries", GetNameId(), name);
1778 if(!(h = (TH2I*)gROOT->FindObject(hname))){
1779 h = new TH2I(hname, htitle, 100, -1., 1., 100, -4.5, 4.5);
1783 const Int_t kNdpt(150);
1784 const Int_t kNspc = 2*AliPID::kSPECIES+1;
1785 Float_t lPt=0.1, lDPt=-.1, lSpc=-5.5;
1786 Float_t binsPt[kNpt+1], binsSpc[kNspc+1], binsDPt[kNdpt+1];
1787 for(Int_t i=0;i<kNpt+1; i++,lPt=TMath::Exp(i*.15)-1.) binsPt[i]=lPt;
1788 for(Int_t i=0; i<kNspc+1; i++,lSpc+=1.) binsSpc[i]=lSpc;
1789 for(Int_t i=0; i<kNdpt+1; i++,lDPt+=2.e-3) binsDPt[i]=lDPt;
1792 snprintf(hname, 100, "%s_%s_Pt", GetNameId(), name);
1793 snprintf(htitle, 300, "#splitline{P_{t} res for}{\"%s\" @ %s};p_{t} [GeV/c];#Delta p_{t}/p_{t}^{MC};SPECIES", GetNameId(), name);
1794 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1795 h = new TH3S(hname, htitle,
1796 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1797 //ax = h->GetZaxis();
1798 //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1802 snprintf(hname, 100, "%s_%s_1Pt", GetNameId(), name);
1803 snprintf(htitle, 300, "#splitline{1/P_{t} pull for}{\"%s\" @ %s};1/p_{t}^{MC} [c/GeV];#Delta(1/p_{t})/#sigma(1/p_{t});SPECIES", GetNameId(), name);
1804 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1805 h = new TH3S(hname, htitle,
1806 kNpt, 0., 2., 100, -4., 4., kNspc, -5.5, 5.5);
1807 //ax = h->GetZaxis();
1808 //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1812 snprintf(hname, 100, "%s_%s_P", GetNameId(), name);
1813 snprintf(htitle, 300, "P res for \"%s\" @ %s;p [GeV/c];#Delta p/p^{MC};SPECIES", GetNameId(), name);
1814 if(!(h = (TH3S*)gROOT->FindObject(hname))){
1815 h = new TH3S(hname, htitle,
1816 kNpt, binsPt, kNdpt, binsDPt, kNspc, binsSpc);
1817 //ax = h->GetZaxis();
1818 //for(Int_t ib(1); ib<=ax->GetNbins(); ib++) ax->SetBinLabel(ib, fgParticle[ib-1]);
1826 //________________________________________________________
1827 TObjArray* AliTRDresolution::Histos()
1830 // Define histograms
1833 if(fContainer) return fContainer;
1835 fContainer = new TObjArray(kNclasses); fContainer->SetOwner(kTRUE);
1837 const Int_t nhn(100); Char_t hn[nhn]; TString st;
1839 //++++++++++++++++++++++
1840 // cluster to tracklet residuals/pulls
1841 snprintf(hn, nhn, "h%s", fgPerformanceName[kCluster]);
1842 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1843 const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", "Q/angle", "#Phi [deg]"};
1844 const Int_t clNbins[kNdim] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 45, 10, 30, 15};
1845 const Double_t clMin[kNdim] = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., -.5, 0.1, -2., -45},
1846 clMax[kNdim] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, 118., 45};
1847 st = "cluster spatial&charge resolution;";
1848 // define minimum info to be saved in non debug mode
1849 Int_t ndim=DebugLevel()>=1?kNdim:4;
1850 for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
1851 H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
1853 fContainer->AddAt(H, kCluster);
1854 //++++++++++++++++++++++
1855 // tracklet to TRD track
1856 snprintf(hn, nhn, "h%s", fgPerformanceName[kTracklet]);
1857 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1858 Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
1859 Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
1860 Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1861 Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1862 // set specific fields
1863 trMin[kYrez] = -0.45; trMax[kYrez] = 0.45;
1864 trMin[kPrez] = -4.5; trMax[kPrez] = 4.5;
1865 trMin[kZrez] = -1.5; trMax[kZrez] = 1.5;
1866 // trNbins[kSpeciesChgRC] = 3;trMin[kSpeciesChgRC] = -1.5; trMax[kSpeciesChgRC] = 1.5;
1867 trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
1868 trTitle[kNdim]=StrDup("dq/dl [a.u.]"); trNbins[kNdim] = 30; trMin[kNdim] = 100.; trMax[kNdim] = 3100;
1870 st = "tracklet spatial&charge resolution;";
1871 // define minimum info to be saved in non debug mode
1872 Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
1873 for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
1874 H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
1876 fContainer->AddAt(H, kTracklet);
1877 //++++++++++++++++++++++
1878 // tracklet to TRDin
1879 snprintf(hn, nhn, "h%s", fgPerformanceName[kTrackIn]);
1880 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1881 // set specific fields
1882 const Int_t mdim(kNdim+2);
1883 Char_t *trinTitle[mdim]; memcpy(trinTitle, fgkTitle, kNdim*sizeof(Char_t*));
1884 Int_t trinNbins[mdim]; memcpy(trinNbins, fgkNbins, kNdim*sizeof(Int_t));
1885 Double_t trinMin[mdim]; memcpy(trinMin, fgkMin, kNdim*sizeof(Double_t));
1886 Double_t trinMax[mdim]; memcpy(trinMax, fgkMax, kNdim*sizeof(Double_t));
1887 trinNbins[kSpeciesChgRC] = 3;trinMin[kSpeciesChgRC] = -1.5; trinMax[kSpeciesChgRC] = 1.5;
1888 trinTitle[kNdim]=StrDup("detector"); trinNbins[kNdim] = 540; trinMin[kNdim] = -0.5; trinMax[kNdim] = 539.5;
1889 trinTitle[kNdim+1]=StrDup("dx [cm]"); trinNbins[kNdim+1]=48; trinMin[kNdim+1]=-2.4; trinMax[kNdim+1]=2.4;
1890 st = "r-#phi/z/angular residuals @ TRD entry;";
1891 // define minimum info to be saved in non debug mode
1892 Int_t ndim=DebugLevel()>=1?mdim:7;
1893 for(Int_t idim(0); idim<ndim; idim++){st+=trinTitle[idim]; st+=";";}
1894 H = new THnSparseI(hn, st.Data(), ndim, trinNbins, trinMin, trinMax);
1896 fContainer->AddAt(H, kTrackIn);
1897 // tracklet to TRDout
1898 // fContainer->AddAt(BuildMonitorContainerTracklet("TrkOUT"), kTrackOut);
1901 // Resolution histos
1902 if(!HasMCdata()) return fContainer;
1904 //++++++++++++++++++++++
1905 // cluster to TrackRef residuals/pulls
1906 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCcluster]);
1907 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1908 const Char_t *clTitle[kNdim] = {"layer", fgkTitle[kPhi], fgkTitle[kEta], fgkTitle[kYrez], "#Deltax [cm]", "Q</Q", fgkTitle[kSpeciesChgRC], "#Phi [deg]"};
1909 const Int_t clNbins[kNdim] = {AliTRDgeometry::kNlayer, fgkNbins[kPhi], fgkNbins[kEta], fgkNbins[kYrez], 20, 10, fgkNbins[kSpeciesChgRC], 15};
1910 const Double_t clMin[kNdim] = {-0.5, fgkMin[kPhi], fgkMin[kEta], fgkMin[kYrez]/10., 0., 0.1, fgkMin[kSpeciesChgRC], -45},
1911 clMax[kNdim] = {AliTRDgeometry::kNlayer-0.5, fgkMax[kPhi], fgkMax[kEta], fgkMax[kYrez]/10., 4., 2.1, fgkMax[kSpeciesChgRC], 45};
1912 st = "MC cluster spatial resolution;";
1913 // define minimum info to be saved in non debug mode
1914 Int_t ndim=DebugLevel()>=1?kNdim:4;
1915 for(Int_t idim(0); idim<ndim; idim++){ st += clTitle[idim]; st+=";";}
1916 H = new THnSparseI(hn, st.Data(), ndim, clNbins, clMin, clMax);
1918 fContainer->AddAt(H, kMCcluster);
1919 //++++++++++++++++++++++
1920 // tracklet to TrackRef
1921 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtracklet]);
1922 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1923 Char_t *trTitle[kNdim]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
1924 Int_t trNbins[kNdim]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
1925 Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1926 Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1927 // set specific fields
1928 trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
1929 trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
1930 trMin[kPrez] = -4.5; trMax[kPrez] = -trMin[kPrez];
1931 trMin[kZrez] = -1.5; trMax[kZrez] = -trMin[kZrez];
1933 st = "MC tracklet spatial resolution;";
1934 // define minimum info to be saved in non debug mode
1935 Int_t ndim=DebugLevel()>=1?kNdim:4;
1936 for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
1937 H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
1939 fContainer->AddAt(H, kMCtracklet);
1940 //++++++++++++++++++++++
1941 // TRDin to TrackRef
1942 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrackIn]);
1943 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1944 st = "MC r-#phi/z/angular residuals @ TRD entry;";
1945 // set specific fields
1946 Double_t trMin[kNdim]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1947 Double_t trMax[kNdim]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1948 trMin[kYrez] = -0.54; trMax[kYrez] = -trMin[kYrez];
1949 trMin[kPrez] = -2.4; trMax[kPrez] = -trMin[kPrez];
1950 trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
1951 // define minimum info to be saved in non debug mode
1952 Int_t ndim=DebugLevel()>=1?kNdim:7;
1953 for(Int_t idim(0); idim<ndim; idim++){ st += fgkTitle[idim]; st+=";";}
1954 H = new THnSparseI(hn, st.Data(), ndim, fgkNbins, trMin, trMax);
1956 fContainer->AddAt(H, kMCtrackIn);
1957 //++++++++++++++++++++++
1958 // track to TrackRef
1959 snprintf(hn, nhn, "h%s", fgPerformanceName[kMCtrack]);
1960 if(!(H = (THnSparseI*)gROOT->FindObject(hn))){
1961 Char_t *trTitle[kNdim+1]; memcpy(trTitle, fgkTitle, kNdim*sizeof(Char_t*));
1962 Int_t trNbins[kNdim+1]; memcpy(trNbins, fgkNbins, kNdim*sizeof(Int_t));
1963 Double_t trMin[kNdim+1]; memcpy(trMin, fgkMin, kNdim*sizeof(Double_t));
1964 Double_t trMax[kNdim+1]; memcpy(trMax, fgkMax, kNdim*sizeof(Double_t));
1965 // set specific fields
1966 trTitle[kBC]=StrDup("layer"); trNbins[kBC] = AliTRDgeometry::kNlayer; trMin[kBC] = -0.5; trMax[kBC] = AliTRDgeometry::kNlayer-0.5;
1967 trTitle[kNdim]=StrDup("#Deltap_{t}/p_{t} [%]"); trNbins[kNdim] = 25; trMin[kNdim] = -4.5; trMax[kNdim] = 20.5;
1968 trMin[kYrez] = -0.9; trMax[kYrez] = -trMin[kYrez];
1969 trMin[kPrez] = -1.5; trMax[kPrez] = -trMin[kPrez];
1970 trMin[kZrez] = -0.9; trMax[kZrez] = -trMin[kZrez];
1972 st = "MC track spatial&p_{t} resolution;";
1973 // define minimum info to be saved in non debug mode
1974 Int_t ndim=DebugLevel()>=1?(kNdim+1):4;
1975 for(Int_t idim(0); idim<ndim; idim++){ st += trTitle[idim]; st+=";";}
1976 H = new THnSparseI(hn, st.Data(), ndim, trNbins, trMin, trMax);
1978 fContainer->AddAt(H, kMCtrack);
1980 // // cluster resolution
1981 // fContainer->AddAt(BuildMonitorContainerCluster("MCcl"), kMCcluster);
1982 // // track resolution
1983 // TObjArray *arr(NULL);
1984 // fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kMCtrack);
1985 // arr->SetName("MCtrk");
1986 // for(Int_t il(0); il<AliTRDgeometry::kNlayer; il++) arr->AddAt(BuildMonitorContainerTrack(Form("MCtrk_Ly%d", il)), il);
1987 // // TRDin TRACK RESOLUTION
1988 // fContainer->AddAt(H, kMCtrackIn);
1989 // // TRDout TRACK RESOLUTION
1990 // fContainer->AddAt(BuildMonitorContainerTrack("MCtrkOUT"), kMCtrackOut);
1995 //________________________________________________________
1996 Bool_t AliTRDresolution::Process(TH2* const h2, TGraphErrors **g, Int_t stat)
1998 // Robust function to process sigma/mean for 2D plot dy(x)
1999 // For each x bin a gauss fit is performed on the y projection and the range
2000 // with the minimum chi2/ndf is choosen
2003 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer input container.\n");
2006 if(!Int_t(h2->GetEntries())){
2007 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : Empty h[%s - %s].\n", h2->GetName(), h2->GetTitle());
2010 if(!g || !g[0]|| !g[1]) {
2011 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>0) printf("D-AliTRDresolution::Process() : NULL pointer output container.\n");
2015 TAxis *ax(h2->GetXaxis()), *ay(h2->GetYaxis());
2016 Float_t ymin(ay->GetXmin()), ymax(ay->GetXmax()), dy(ay->GetBinWidth(1)), y0(0.), y1(0.);
2017 TF1 f("f", "gaus", ymin, ymax);
2019 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
2020 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
2022 if((h=(TH1D*)gROOT->FindObject("py"))) delete h;
2023 Double_t x(0.), y(0.), ex(0.), ey(0.), sy(0.), esy(0.);
2027 Float_t chi2OverNdf(0.);
2028 for(Int_t ix = 1, np=0; ix<=ax->GetNbins(); ix++){
2029 x = ax->GetBinCenter(ix); ex= ax->GetBinWidth(ix)*0.288; // w/sqrt(12)
2030 ymin = ay->GetXmin(); ymax = ay->GetXmax();
2032 h = h2->ProjectionY("py", ix, ix);
2033 if((n=(Int_t)h->GetEntries())<stat){
2034 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("I-AliTRDresolution::Process() : Low statistics @ x[%f] stat[%d]=%d [%d].\n", x, ix, n, stat);
2037 // looking for a first order mean value
2038 f.SetParameter(1, 0.); f.SetParameter(2, 3.e-2);
2040 chi2OverNdf = f.GetChisquare()/f.GetNDF();
2041 printf("x[%f] range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", x, ymin, ymax, f.GetChisquare(),f.GetNDF(),chi2OverNdf);
2042 y = f.GetParameter(1); y0 = y-4*dy; y1 = y+4*dy;
2043 ey = f.GetParError(1);
2044 sy = f.GetParameter(2); esy = f.GetParError(2);
2045 // // looking for the best chi2/ndf value
2046 // while(ymin<y0 && ymax>y1){
2047 // f.SetParameter(1, y);
2048 // f.SetParameter(2, sy);
2049 // h->Fit(&f, "QNW", "", y0, y1);
2050 // printf(" range[%f %f] chi2[%f] ndf[%d] chi2/ndf[%f]\n", y0, y1, f.GetChisquare(),f.GetNDF(),f.GetChisquare()/f.GetNDF());
2051 // if(f.GetChisquare()/f.GetNDF() < Chi2OverNdf){
2052 // chi2OverNdf = f.GetChisquare()/f.GetNDF();
2053 // y = f.GetParameter(1); ey = f.GetParError(1);
2054 // sy = f.GetParameter(2); esy = f.GetParError(2);
2055 // printf(" set y[%f] sy[%f] chi2/ndf[%f]\n", y, sy, chi2OverNdf);
2059 g[0]->SetPoint(np, x, y);
2060 g[0]->SetPointError(np, ex, ey);
2061 g[1]->SetPoint(np, x, sy);
2062 g[1]->SetPointError(np, ex, esy);
2069 //________________________________________________________
2070 Bool_t AliTRDresolution::Process(TH2 * const h2, TF1 *f, Float_t k, TGraphErrors **g)
2073 // Do the processing
2076 Char_t pn[10]; snprintf(pn, 10, "p%03d", fIdxPlot);
2078 if((n=g[0]->GetN())) for(;n--;) g[0]->RemovePoint(n);
2079 if((n=g[1]->GetN())) for(;n--;) g[1]->RemovePoint(n);
2080 if(Int_t(h2->GetEntries())){
2081 AliDebug(4, Form("%s: g[%s %s]", pn, g[0]->GetName(), g[0]->GetTitle()));
2083 AliDebug(2, Form("%s: g[%s %s]: Missing entries.", pn, g[0]->GetName(), g[0]->GetTitle()));
2088 const Int_t kINTEGRAL=1;
2089 for(Int_t ibin = 0; ibin < Int_t(h2->GetNbinsX()/kINTEGRAL); ibin++){
2090 Int_t abin(ibin*kINTEGRAL+1),bbin(abin+kINTEGRAL-1),mbin(abin+Int_t(kINTEGRAL/2));
2091 Double_t x = h2->GetXaxis()->GetBinCenter(mbin);
2092 TH1D *h = h2->ProjectionY(pn, abin, bbin);
2093 if((n=(Int_t)h->GetEntries())<150){
2094 AliDebug(4, Form(" x[%f] range[%d %d] stat[%d] low statistics !", x, abin, bbin, n));
2098 Int_t ip = g[0]->GetN();
2099 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)));
2100 g[0]->SetPoint(ip, x, k*f->GetParameter(1));
2101 g[0]->SetPointError(ip, 0., k*f->GetParError(1));
2102 g[1]->SetPoint(ip, x, k*f->GetParameter(2));
2103 g[1]->SetPointError(ip, 0., k*f->GetParError(2));
2105 g[0]->SetPoint(ip, x, k*h->GetMean());
2106 g[0]->SetPointError(ip, 0., k*h->GetMeanError());
2107 g[1]->SetPoint(ip, x, k*h->GetRMS());
2108 g[1]->SetPointError(ip, 0., k*h->GetRMSError());*/
2115 //____________________________________________________________________
2116 Bool_t AliTRDresolution::FitTrack(const Int_t np, AliTrackPoint *points, Float_t param[10])
2119 // Fit track with a staight line using the "np" clusters stored in the array "points".
2120 // The following particularities are stored in the clusters from points:
2121 // 1. pad tilt as cluster charge
2122 // 2. pad row cross or vertex constrain as fake cluster with cluster type 1
2123 // The parameters of the straight line fit are stored in the array "param" in the following order :
2124 // param[0] - x0 reference radial position
2125 // param[1] - y0 reference r-phi position @ x0
2126 // param[2] - z0 reference z position @ x0
2127 // param[3] - slope dy/dx
2128 // param[4] - slope dz/dx
2131 // Function should be used to refit tracks for B=0T
2135 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTrack: Not enough clusters to fit a track [%d].\n", np);
2138 TLinearFitter yfitter(2, "pol1"), zfitter(2, "pol1");
2141 for(Int_t ip(0); ip<np; ip++) x0+=points[ip].GetX();
2144 Double_t x, y, z, dx, tilt(0.);
2145 for(Int_t ip(0); ip<np; ip++){
2146 x = points[ip].GetX(); z = points[ip].GetZ();
2148 zfitter.AddPoint(&dx, z, points[ip].GetClusterType()?1.e-3:1.);
2150 if(zfitter.Eval() != 0) return kFALSE;
2152 Double_t z0 = zfitter.GetParameter(0);
2153 Double_t dzdx = zfitter.GetParameter(1);
2154 for(Int_t ip(0); ip<np; ip++){
2155 if(points[ip].GetClusterType()) continue;
2156 x = points[ip].GetX();
2158 y = points[ip].GetY();
2159 z = points[ip].GetZ();
2160 tilt = points[ip].GetCharge();
2161 y -= tilt*(-dzdx*dx + z - z0);
2162 Float_t xyz[3] = {x, y, z}; points[ip].SetXYZ(xyz);
2163 yfitter.AddPoint(&dx, y, 1.);
2165 if(yfitter.Eval() != 0) return kFALSE;
2166 Double_t y0 = yfitter.GetParameter(0);
2167 Double_t dydx = yfitter.GetParameter(1);
2169 param[0] = x0; param[1] = y0; param[2] = z0; param[3] = dydx; param[4] = dzdx;
2170 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>3) printf("D-AliTRDresolution::FitTrack: x0[%f] y0[%f] z0[%f] dydx[%f] dzdx[%f].\n", x0, y0, z0, dydx, dzdx);
2174 //____________________________________________________________________
2175 Bool_t AliTRDresolution::FitTracklet(const Int_t ly, const Int_t np, const AliTrackPoint *points, const Float_t param[10], Float_t par[3])
2178 // Fit tracklet with a staight line using the coresponding subset of clusters out of the total "np" clusters stored in the array "points".
2179 // See function FitTrack for the data stored in the "clusters" array
2181 // The parameters of the straight line fit are stored in the array "param" in the following order :
2182 // par[0] - x0 reference radial position
2183 // par[1] - y0 reference r-phi position @ x0
2184 // par[2] - slope dy/dx
2187 // Function should be used to refit tracks for B=0T
2190 TLinearFitter yfitter(2, "pol1");
2192 // grep data for tracklet
2193 Double_t x0(0.), x[60], y[60], dy[60];
2195 for(Int_t ip(0); ip<np; ip++){
2196 if(points[ip].GetClusterType()) continue;
2197 if(points[ip].GetVolumeID() != ly) continue;
2198 Float_t xt(points[ip].GetX())
2199 ,yt(param[1] + param[3] * (xt - param[0]));
2201 y[nly] = points[ip].GetY();
2207 if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>1) printf("D-AliTRDresolution::FitTracklet: Not enough clusters to fit a tracklet [%d].\n", nly);
2210 // set radial reference for fit
2213 // find tracklet core
2214 Double_t mean(0.), sig(1.e3);
2215 AliMathBase::EvaluateUni(nly, dy, mean, sig, 0);
2217 // simple cluster error parameterization
2218 Float_t kSigCut = TMath::Sqrt(5.e-4 + param[3]*param[3]*0.018);
2220 // fit tracklet core
2221 for(Int_t jly(0); jly<nly; jly++){
2222 if(TMath::Abs(dy[jly]-mean)>kSigCut) continue;
2223 Double_t dx(x[jly]-x0);
2224 yfitter.AddPoint(&dx, y[jly], 1.);
2226 if(yfitter.Eval() != 0) return kFALSE;
2228 par[1] = yfitter.GetParameter(0);
2229 par[2] = yfitter.GetParameter(1);
2233 //____________________________________________________________________
2234 Bool_t AliTRDresolution::UseTrack(const Int_t np, const AliTrackPoint *points, Float_t param[10])
2237 // Global selection mechanism of tracksbased on cluster to fit residuals
2238 // The parameters are the same as used ni function FitTrack().
2240 const Float_t kS(0.6), kM(0.2);
2241 TH1S h("h1", "", 100, -5.*kS, 5.*kS);
2242 Float_t dy, dz, s, m;
2243 for(Int_t ip(0); ip<np; ip++){
2244 if(points[ip].GetClusterType()) continue;
2245 Float_t x0(points[ip].GetX())
2246 ,y0(param[1] + param[3] * (x0 - param[0]))
2247 ,z0(param[2] + param[4] * (x0 - param[0]));
2248 dy=points[ip].GetY() - y0; h.Fill(dy);
2249 dz=points[ip].GetZ() - z0;
2251 TF1 fg("fg", "gaus", -5.*kS, 5.*kS);
2252 fg.SetParameter(1, 0.);
2253 fg.SetParameter(2, 2.e-2);
2255 m=fg.GetParameter(1); s=fg.GetParameter(2);
2256 if(s>kS || TMath::Abs(m)>kM) return kFALSE;
2260 //________________________________________________________
2261 void AliTRDresolution::GetLandauMpvFwhm(TF1 * const f, Float_t &mpv, Float_t &xm, Float_t &xM)
2264 // Get the most probable value and the full width half mean
2265 // of a Landau distribution
2268 const Float_t dx = 1.;
2269 mpv = f->GetParameter(1);
2270 Float_t fx, max = f->Eval(mpv);
2273 while((fx = f->Eval(xm))>.5*max){
2282 while((fx = f->Eval(xM))>.5*max) xM += dx;
2286 // #include "TFile.h"
2287 // //________________________________________________________
2288 // Bool_t AliTRDresolution::LoadCorrection(const Char_t *file)
2291 // AliWarning("Use cluster position as in reconstruction.");
2292 // SetLoadCorrection();
2295 // TDirectory *cwd(gDirectory);
2296 // TString fileList;
2297 // FILE *filePtr = fopen(file, "rt");
2299 // AliWarning(Form("Couldn't open correction list \"%s\". Use cluster position as in reconstruction.", file));
2300 // SetLoadCorrection();
2303 // TH2 *h2 = new TH2F("h2", ";time [time bins];detector;dx [#mum]", 30, -0.5, 29.5, AliTRDgeometry::kNdet, -0.5, AliTRDgeometry::kNdet-0.5);
2304 // while(fileList.Gets(filePtr)){
2305 // if(!TFile::Open(fileList.Data())) {
2306 // AliWarning(Form("Couldn't open \"%s\"", fileList.Data()));
2308 // } else AliInfo(Form("\"%s\"", fileList.Data()));
2310 // TTree *tSys = (TTree*)gFile->Get("tSys");
2311 // h2->SetDirectory(gDirectory); h2->Reset("ICE");
2312 // tSys->Draw("det:t>>h2", "dx", "goff");
2313 // for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
2314 // for(Int_t it(0); it<30; it++) fXcorr[idet][it]+=(1.e-4*h2->GetBinContent(it+1, idet+1));
2316 // h2->SetDirectory(cwd);
2321 // if(AliLog::GetDebugLevel("PWG1", "AliTRDresolution")>=2){
2322 // for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
2323 // printf("DET|");for(Int_t it(0); it<30; it++) printf(" tb%02d|", it); printf("\n");
2324 // for(Int_t il(0); il<184; il++) printf("-"); printf("\n");
2325 // FILE *fout = fopen("TRD.ClusterCorrection.txt", "at");
2326 // fprintf(fout, " static const Double_t dx[AliTRDgeometry::kNdet][30]={\n");
2327 // for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){
2328 // printf("%03d|", idet);
2329 // fprintf(fout, " {");
2330 // for(Int_t it(0); it<30; it++){
2331 // printf("%+5.0f|", 1.e4*fXcorr[idet][it]);
2332 // fprintf(fout, "%+6.4f%c", fXcorr[idet][it], it==29?' ':',');
2335 // fprintf(fout, "}%c\n", idet==AliTRDgeometry::kNdet-1?' ':',');
2337 // fprintf(fout, " };\n");
2339 // SetLoadCorrection();
2343 //________________________________________________________
2344 AliTRDresolution::AliTRDresolutionProjection::AliTRDresolutionProjection()
2351 memset(fAx, 0, 3*sizeof(Int_t));
2352 memset(fRange, 0, 4*sizeof(Float_t));
2355 //________________________________________________________
2356 AliTRDresolution::AliTRDresolutionProjection::~AliTRDresolutionProjection()
2362 //________________________________________________________
2363 void AliTRDresolution::AliTRDresolutionProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[])
2365 // check and build (if neccessary) projection determined by axis "ix", "iy" and "iz"
2366 if(!aa[ix] || !aa[iy] || !aa[iz]) return;
2367 TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]);
2368 fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()),
2369 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),
2370 ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),
2371 az->GetNbins(), az->GetXmin(), az->GetXmax());
2372 fAx[0] = ix; fAx[1] = iy; fAx[2] = iz;
2373 fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.;
2376 //________________________________________________________
2377 AliTRDresolution::AliTRDresolutionProjection& AliTRDresolution::AliTRDresolutionProjection::operator+=(const AliTRDresolutionProjection& other)
2379 // increment projections
2380 if(!fH || !other.fH) return *this;
2385 //________________________________________________________
2386 void AliTRDresolution::AliTRDresolutionProjection::Increment(Int_t bin[], Double_t v)
2388 // increment bin with value "v" pointed by general coord in "bin"
2390 fH->AddBinContent(fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), Int_t(v));
2393 //________________________________________________________
2394 TH2* AliTRDresolution::AliTRDresolutionProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid)
2396 // build the 2D projection and adjust binning
2398 const Char_t *title[] = {"Mean", "#mu", "MPV"};
2399 if(!fH) return NULL;
2400 TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());
2402 if(!(h2s = (TH2*)fH->Project3D("yx"))) return NULL;
2403 Int_t irebin(0), dxBin(1), dyBin(1);
2404 while(irebin<fNrebin && (AliTRDresolution::GetMeanStat(h2s, .5, ">")<nstat)){
2405 h2s->Rebin2D(fRebinX[irebin], fRebinY[irebin]);
2406 dxBin*=fRebinX[irebin];dyBin*=fRebinY[irebin];
2409 Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());
2414 Float_t dz=(fRange[1]-fRange[1])/ncol;
2415 TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" "));
2416 TH2 *h2 = new TH2F(Form("%s_2D", fH->GetName()),
2417 Form("%s;%s;%s;%s(%s) %s", fH->GetTitle(), ax->GetTitle(), ay->GetTitle(), title[mid], (*tokenTitle)[0]->GetName(), tokenTitle->GetEntriesFast()>1?(*tokenTitle)[1]->GetName():""),
2418 nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax());
2419 h2->SetContour(ncol);
2420 h2->GetZaxis()->CenterTitle();
2421 h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]);
2422 //printf("%s[%s] nx[%d] ny[%d]\n", h2->GetName(), h2->GetTitle(), nx, ny);
2423 for(Int_t iy(0); iy<ny; iy++){
2424 for(Int_t ix(0); ix<nx; ix++){
2425 h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin+1, iy*dyBin+1, (iy+1)*dyBin+1);
2426 Int_t ne((Int_t)h->Integral());
2428 h2->SetBinContent(ix+1, iy+1, -999);
2429 h2->SetBinError(ix+1, iy+1, 1.);
2431 Float_t v(h->GetMean()), ve(h->GetRMS());
2433 TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());
2434 fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve);
2436 v = fg.GetParameter(1); ve = fg.GetParameter(2);
2437 } else if (mid==2) {
2438 TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());
2439 fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, v); fl.SetParameter(2, ve);
2441 v = fl.GetMaximumX(); ve = fl.GetParameter(2);
2442 /* TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax());
2443 fgle.SetParameter(0, fl.GetParameter(0));
2444 fgle.SetParameter(1, fl.GetParameter(1));
2445 fgle.SetParameter(2, fl.GetParameter(2));
2446 fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.);
2447 h->Fit(&fgle, "WQ");
2448 v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/
2450 if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz);
2451 else h2->SetBinContent(ix+1, iy+1, v);
2452 h2->SetBinError(ix+1, iy+1, ve);
2460 void AliTRDresolution::AliTRDresolutionProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[])
2462 // define rebinning strategy for this projection
2464 fRebinX = new Int_t[n]; memcpy(fRebinX, rebx, n*sizeof(Int_t));
2465 fRebinY = new Int_t[n]; memcpy(fRebinY, reby, n*sizeof(Int_t));