From 1fd9389f2c4ab5c90db974b34b282e78e3b1e858 Mon Sep 17 00:00:00 2001 From: abercuci Date: Fri, 5 Jun 2009 14:57:27 +0000 Subject: [PATCH] restructure error calculation : - move all cluster error parameterization to AliTRDcluster - re-implement leading cluster method for radial calculation of pad row cross tracklets - implement pull plots for cluster calibration fix minor errors in performance tasks update documentation --- TRD/AliTRDcluster.cxx | 238 +++++++++++++++++++++-- TRD/AliTRDcluster.h | 5 +- TRD/AliTRDseedV1.cxx | 260 ++++++++++++++------------ TRD/AliTRDseedV1.h | 38 ++-- TRD/AliTRDtracker.cxx | 2 +- TRD/AliTRDtransform.cxx | 16 +- TRD/qaRec/AliTRDcheckDET.cxx | 57 +++--- TRD/qaRec/AliTRDcheckPID.cxx | 33 ++-- TRD/qaRec/AliTRDclusterResolution.cxx | 185 +++++++++++++----- TRD/qaRec/AliTRDclusterResolution.h | 15 +- TRD/qaRec/AliTRDrecoTask.h | 2 +- TRD/qaRec/AliTRDresolution.cxx | 185 +++++++++--------- TRD/qaRec/info/AliTRDclusterInfo.cxx | 7 +- TRD/qaRec/info/AliTRDclusterInfo.h | 2 +- TRD/qaRec/macros/makeResults.C | 3 +- 15 files changed, 684 insertions(+), 364 deletions(-) diff --git a/TRD/AliTRDcluster.cxx b/TRD/AliTRDcluster.cxx index 3a693f6c6b8..18ad6955740 100644 --- a/TRD/AliTRDcluster.cxx +++ b/TRD/AliTRDcluster.cxx @@ -136,7 +136,7 @@ AliTRDcluster::AliTRDcluster(const AliTRDcluster &c) SetY(c.GetY()); SetZ(c.GetZ()); - SetSigmaY2(c.GetSigmaY2()); + AliCluster::SetSigmaY2(c.GetSigmaY2()); SetSigmaZ2(c.GetSigmaZ2()); for (Int_t i = 0; i < 7; i++) { @@ -241,11 +241,11 @@ void AliTRDcluster::Clear(Option_t *) fQ = 0; fCenter = 0; for (Int_t i = 0; i < 3; i++) SetLabel(0,i); - SetX(0); - SetY(0); - SetZ(0); - SetSigmaY2(0); - SetSigmaZ2(0); + SetX(0.); + SetY(0.); + SetZ(0.); + AliCluster::SetSigmaY2(0.); + SetSigmaZ2(0.); SetVolumeId(0); } @@ -268,6 +268,19 @@ Float_t AliTRDcluster::GetSumS() const //___________________________________________________________________________ Double_t AliTRDcluster::GetSX(Int_t tb, Double_t z) { +// Returns the error parameterization in the radial direction for TRD clusters as function of +// the calibrated time bin (tb) and optionally distance to anode wire (z). By default (no z information) +// the largest value over all cluster to wire values is chosen. +// +// The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. +// +//Begin_Html +// +//End_Html +// +// Author +// A.Bercuci + if(tb<1 || tb>=24) return 10.; // return huge [10cm] const Double_t sx[24][10]={ {0.000e+00, 9.352e-01, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 2.309e+00}, @@ -302,10 +315,68 @@ Double_t AliTRDcluster::GetSX(Int_t tb, Double_t z) } //___________________________________________________________________________ -Double_t AliTRDcluster::GetSY(Int_t tb, Double_t z) +Double_t AliTRDcluster::GetSYdrift(Int_t tb, Int_t ly, Double_t/* z*/) { +// Returns the error parameterization for TRD clusters as function of the drift length (here calibrated time bin tb) +// and optionally distance to anode wire (z) for the LUT r-phi cluster shape. By default (no z information) the largest +// value over all cluster to wire values is chosen. +// +// For the LUT method the dependence of s_y with x and d is obtained via a fit to the cluster to MC +// resolution. (see class AliTRDclusterResolution for more details). A normalization to the reference radial position +// x0 = 0.675 (tb=5 for ideal vd) is also applied (see GetSYprf()). The function is *NOT* calibration aware ! +// The result is displayed in the figure below as a 2D plot and also as the projection on the drift axis. A comparison +// with the GAUS parameterization is also given +// +// For the GAUS method the dependence of s_y with x is *analytic* and it is expressed by the relation. +// BEGIN_LATEX +// #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}} +// END_LATEX +// The result is displayed in the figure below as function of the drift time and compared with the LUT parameterization. +//Begin_Html +// +// +//End_Html +// +// Author +// A.Bercuci + + if(tb<1 || tb>=24) return 10.; // return huge [10cm] - const Double_t sy[24][10]={ + const Float_t lSy[6][24] = { + {75.7561, 0.0325, 0.0175, 0.0174, 0.0206, 0.0232, + 0.0253, 0.0262, 0.0265, 0.0264, 0.0266, 0.0257, + 0.0258, 0.0261, 0.0259, 0.0253, 0.0257, 0.0261, + 0.0255, 0.0250, 0.0259, 0.0266, 0.0278, 0.0319 + }, + {49.2252, 0.0371, 0.0204, 0.0189, 0.0230, 0.0261, + 0.0281, 0.0290, 0.0292, 0.0286, 0.0277, 0.0279, + 0.0285, 0.0281, 0.0291, 0.0281, 0.0281, 0.0282, + 0.0272, 0.0282, 0.0282, 0.0284, 0.0310, 0.0334 + }, + {55.1674, 0.0388, 0.0212, 0.0200, 0.0239, 0.0271, + 0.0288, 0.0299, 0.0306, 0.0300, 0.0296, 0.0303, + 0.0293, 0.0290, 0.0291, 0.0294, 0.0295, 0.0290, + 0.0293, 0.0292, 0.0292, 0.0293, 0.0316, 0.0358 + }, + {45.1004, 0.0411, 0.0225, 0.0215, 0.0249, 0.0281, + 0.0301, 0.0315, 0.0320, 0.0308, 0.0318, 0.0321, + 0.0312, 0.0311, 0.0316, 0.0315, 0.0310, 0.0308, + 0.0313, 0.0303, 0.0314, 0.0314, 0.0324, 0.0369 + }, + {43.8614, 0.0420, 0.0239, 0.0224, 0.0268, 0.0296, + 0.0322, 0.0336, 0.0333, 0.0326, 0.0321, 0.0325, + 0.0329, 0.0326, 0.0323, 0.0322, 0.0326, 0.0320, + 0.0329, 0.0319, 0.0314, 0.0329, 0.0341, 0.0373 + }, + {40.5440, 0.0434, 0.0246, 0.0236, 0.0275, 0.0311, + 0.0332, 0.0345, 0.0347, 0.0347, 0.0340, 0.0336, + 0.0339, 0.0344, 0.0339, 0.0341, 0.0341, 0.0342, + 0.0345, 0.0328, 0.0341, 0.0332, 0.0356, 0.0398 + }, + }; + return lSy[ly][tb]; + +/* const Double_t sy[24][10]={ {0.000e+00, 2.610e-01, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 0.000e+00, 4.680e-01}, {3.019e-02, 3.036e-02, 3.131e-02, 3.203e-02, 3.294e-02, 3.407e-02, 3.555e-02, 3.682e-02, 3.766e-02, 3.824e-02}, {1.773e-02, 1.778e-02, 1.772e-02, 1.790e-02, 1.807e-02, 1.833e-02, 1.873e-02, 1.905e-02, 1.958e-02, 2.029e-02}, @@ -331,19 +402,83 @@ Double_t AliTRDcluster::GetSY(Int_t tb, Double_t z) {2.812e-02, 2.786e-02, 2.776e-02, 2.723e-02, 2.695e-02, 2.650e-02, 2.642e-02, 2.617e-02, 2.612e-02, 2.610e-02}, {3.251e-02, 3.267e-02, 3.223e-02, 3.183e-02, 3.125e-02, 3.106e-02, 3.067e-02, 3.010e-02, 2.936e-02, 2.927e-02} }; - if(z>=0. && z<.25) return sy[tb][Int_t(z/.025)]; + if(z>=0. && z<.25) return sy[tb][Int_t(z/.025)] - sy[5][Int_t(z/.025)]; - Double_t m = 1.e-8; for(Int_t id=10; id--;) if(sy[tb][id]>m) m=sy[tb][id]; + Double_t m = -1.e8; for(Int_t id=10; id--;) if((sy[tb][id] - sy[5][id])>m) m=sy[tb][id]-sy[5][id]; - return m; + return m;*/ +} + +//___________________________________________________________________________ +Double_t AliTRDcluster::GetSYcharge(Float_t q) +{ +// Parameterization of the r-phi resolution component due to cluster charge. +// The value is the offset from the nominal cluster resolution defined as the +// cluster resolution at average cluster charge (q0). +// +// BEGIN_LATEX +// #Delta #sigma_{y}(q) = a*(#frac{1}{q} - #frac{1}{q_{0}}) +// q_{0} #approx 50 +// END_LATEX +// The definition is *NOT* robust against gain fluctuations and thus two approaches are possible +// when residual miscalibration are available: +// - determine parameterization with a resolution matching those of the gain +// - define an analytic model which scales with the gain. +// +// For more details please see AliTRDclusterResolution::ProcessCharge() +// +//Begin_Html +// +//End_Html +// +// Author +// A.Bercuci + + const Float_t sq0inv = 0.019962; // [1/q0] + const Float_t sqb = 1.0281564;// [cm] + + return sqb*(1./q - sq0inv); +} + +//___________________________________________________________________________ +Double_t AliTRDcluster::GetSYprf(Int_t ly, Double_t center, Double_t s2) +{ +// Parameterization of the cluster error in the r-phi direction due to charge sharing between +// adiacent pads. Should be identical to what is provided in the OCDB as PRF [TODO] +// +// The parameterization is obtained from fitting cluster resolution at phi=exb and |x-0.675|<0.225. +// For more details see AliTRDclusterResolution::ProcessCenter(). +// +//Begin_Html +// +//End_Html +// +// Author +// A.Bercuci + +/* const Float_t scy[AliTRDgeometry::kNlayer][4] = { + {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02}, + {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02}, + {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02}, + {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02}, + {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02}, + {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02}, + };*/ + const Float_t lPRF[] = {0.438, 0.403, 0.393, 0.382, 0.376, 0.345}; + + return s2*TMath::Gaus(center, 0., lPRF[ly]); } + //___________________________________________________________________________ Double_t AliTRDcluster::GetXcorr(Int_t tb, Double_t z) { - // drift length correction [cm] - // TODO to be parametrized in term of drift velocity - // A.Bercuci (Mar 28 2009) +// Drift length correction [cm]. Due to variation of mean drift velocity along the drift region +// from nominal vd at xd->infinity. For drift velocity determination based on tracking information +// the correction should be negligible. +// +// TODO to be parametrized in term of drift velocity at infinite drift length +// A.Bercuci (Mar 28 2009) if(tb<0 || tb>=24) return 0.; const Int_t nd = 5; @@ -408,8 +543,8 @@ Double_t AliTRDcluster::GetXcorr(Int_t tb, Double_t z) //___________________________________________________________________________ Double_t AliTRDcluster::GetYcorr(Int_t ly, Float_t y) { -// PRF correction TODO to be replaced by the gaussian -// approximation with full error parametrization +// PRF correction for the LUT r-phi cluster shape. + const Float_t cy[AliTRDgeometry::kNlayer][3] = { { 4.014e-04, 8.605e-03, -6.880e+00}, {-3.061e-04, 9.663e-03, -6.789e+00}, @@ -465,7 +600,7 @@ Float_t AliTRDcluster::GetXloc(Double_t t0, Double_t vd, Double_t *const /*q*/, // correction for t0 td -= t0; // time bin corrected for t0 - // BUG in TMath::Nint().root-5.23.02 + // Bug in TMath::Nint().root-5.23.02 // TMath::Nint(3.5) = 4 and TMath::Nint(4.5) = 4 Double_t tmp = td*fFreq; fLocalTimeBin = Char_t(TMath::Floor(tmp)); @@ -483,6 +618,8 @@ Float_t AliTRDcluster::GetXloc(Double_t t0, Double_t vd, Double_t *const /*q*/, return x; /* + // calculate radial posion of clusters in the drift region + // invert drift time function Double_t xM= AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght(), x = vd*td + .5*AliTRDgeometry::CamHght(), @@ -510,7 +647,8 @@ Float_t AliTRDcluster::GetYloc(Double_t y0, Double_t s2, Double_t W, Double_t *c //printf(" s[%3d %3d %3d] w[%f %f] yr[%f %f]\n", fSignals[2], fSignals[3], fSignals[4], w1/(w1+w2), w2/(w1+w2), y1r*W, y2r*W); if(IsRPhiMethod(kCOG)) GetDYcog(); else if(IsRPhiMethod(kLUT)) GetDYlut(); - if(IsRPhiMethod(kGAUS)) GetDYgauss(s2/W/W, y1, y2); + else if(IsRPhiMethod(kGAUS)) GetDYgauss(s2/W/W, y1, y2); + else return 0.; if(y1) (*y1)*=W; if(y2) (*y2)*=W; @@ -518,6 +656,70 @@ Float_t AliTRDcluster::GetYloc(Double_t y0, Double_t s2, Double_t W, Double_t *c return y0+fCenter*W+(IsRPhiMethod(kLUT)?GetYcorr(AliTRDgeometry::GetLayer(fDetector), fCenter):0.); } +//___________________________________________________________________________ +void AliTRDcluster::SetSigmaY2(Float_t s2, Float_t dt, Float_t exb, Float_t x, Float_t z, Float_t tgp) +{ +// Set variance of TRD cluster in the r-phi direction for each method. +// Parameters : +// - s2 - variance due to PRF width for the case of Gauss model. Replaced by parameterization in case of LUT. +// - dt - transversal diffusion coeficient +// - exb - tg of lorentz angle +// - x - drift length - with respect to the anode wire +// - z - offset from the anode wire +// - tgp - local tangent of the track momentum azimuthal angle +// +// The ingredients from which the error is computed are: +// - PRF (charge sharing on adjacent pads) - see AliTRDcluster::GetSYprf() +// - diffusion (dependence with drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSYdrift() +// - charge of the cluster (complex dependence on gain and tail cancellation) - see AliTRDcluster::GetSYcharge() +// - lorentz angle (dependence on the drift length and [2nd order] distance to anode wire) - see AliTRDcluster::GetSX() +// - track angle (superposition of charges on the anode wire) - see AliTRDseedV1::Fit() +// - projection of radial(x) error on r-phi due to fixed value assumed in tracking for x - see AliTRDseedV1::Fit() +// +// The last 2 contributions to cluster error can be estimated only during tracking when the track angle +// is known (tgp). For this reason the errors (and optional position) of TRD clusters are recalculated during +// tracking and thus clusters attached to tracks might differ from bare clusters. +// +// Author: +// A.Bercuci + + Float_t sigmaY2 = 0.; + Int_t ly = AliTRDgeometry::GetLayer(fDetector); + if(IsRPhiMethod(kCOG)) sigmaY2 = 4.e-4; + else if(IsRPhiMethod(kLUT)){ + Float_t sd = GetSYdrift(fLocalTimeBin, ly, z);//printf("drift[%6.2f] ", 1.e4*sd); + sigmaY2 = GetSYprf(ly, fCenter, sd);//printf("PRF[%6.2f] ", 1.e4*sigmaY2); + // add charge contribution TODO scale with respect to s2 + sigmaY2+= GetSYcharge(TMath::Abs(fQ));//printf("Q[%6.2f] ", 1.e4*sigmaY2); + sigmaY2 = TMath::Max(sigmaY2, Float_t(0.)); //!! protection + sigmaY2*= sigmaY2; + } else if(IsRPhiMethod(kGAUS)){ + // PRF contribution + sigmaY2 = s2; + // Diffusion contribution + Double_t sD2 = dt/(1.+exb); sD2 *= sD2; sD2 *= x; + sigmaY2+= sD2; + // add charge contribution TODO scale with respect to s2 + //sigmaY2+= GetSYcharge(TMath::Abs(fQ)); + } + + // store tg^2(phi-a_L) and tg^2(a_L) + Double_t tgg = (tgp-exb)/(1.+tgp*exb); tgg *= tgg; + Double_t exb2= exb*exb; + + // Lorentz angle shift contribution + Float_t sx = GetSX(fLocalTimeBin, z); sx*=sx; + sigmaY2+= exb2*sx;//printf("Al[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2)); + + // Radial contribution due to not measuring x in Kalman model + sigmaY2+= tgg*sx;//printf("x[%6.2f] ", 1.e4*TMath::Sqrt(sigmaY2)); + + // Track angle contribution + sigmaY2+= tgg*x*x*exb2/12.;//printf("angle[%6.2f]\n", 1.e4*TMath::Sqrt(sigmaY2)); + + AliCluster::SetSigmaY2(sigmaY2); +} + //_____________________________________________________________________________ Bool_t AliTRDcluster::IsEqual(const TObject *o) const { diff --git a/TRD/AliTRDcluster.h b/TRD/AliTRDcluster.h index 39c820bc4cb..409cf11e0f8 100644 --- a/TRD/AliTRDcluster.h +++ b/TRD/AliTRDcluster.h @@ -62,7 +62,9 @@ public: Short_t *GetSignals() { return fSignals; } Float_t GetSumS() const; static Double_t GetSX(Int_t tb, Double_t z=-1); - static Double_t GetSY(Int_t tb, Double_t z=-1); + static Double_t GetSYdrift(Int_t tb, Int_t ly=0, Double_t z=-1); + static Double_t GetSYcharge(Float_t q); + static Double_t GetSYprf(Int_t ly, Double_t center, Double_t s2); static Double_t GetXcorr(Int_t tb, Double_t z=-1); static Double_t GetYcorr(Int_t ly, Float_t y); Float_t GetXloc(Double_t t0, Double_t vd, Double_t *const q=0x0, Double_t *const xq=0x0, Double_t z = 0.2); @@ -83,6 +85,7 @@ public: void SetQ(Float_t inQ){ fQ = inQ;} void SetClusterMasking(UChar_t inClusterMasking){ fClusterMasking = inClusterMasking;} void SetShared(Bool_t sh = kTRUE) { SetBit(AliCluster::kShared,sh); } + void SetSigmaY2(Float_t s2, Float_t dt, Float_t exb, Float_t x0, Float_t z=-1., Float_t tgp=0.); void Use(Int_t u = 1) { SetBit(AliCluster::kUsed, u ? kTRUE : kFALSE); } void SetFivePad(Bool_t b = kTRUE) { SetBit(kFivePad,b);} diff --git a/TRD/AliTRDseedV1.cxx b/TRD/AliTRDseedV1.cxx index 8cd43f2b490..be2f38cfe26 100644 --- a/TRD/AliTRDseedV1.cxx +++ b/TRD/AliTRDseedV1.cxx @@ -293,8 +293,9 @@ void AliTRDseedV1::Update(const AliTRDtrackV1 *trk) Double_t fSnp = trk->GetSnp(); Double_t fTgl = trk->GetTgl(); fPt = trk->Pt(); - fYref[1] = fSnp/TMath::Sqrt(1. - fSnp*fSnp); - fZref[1] = fTgl; + Double_t norm =1./TMath::Sqrt(1. - fSnp*fSnp); + fYref[1] = fSnp*norm; + fZref[1] = fTgl*norm; SetCovRef(trk->GetCovariance()); Double_t dx = trk->GetX() - fX0; @@ -474,6 +475,7 @@ Float_t AliTRDseedV1::GetdQdl(Int_t ic, Float_t *dl) const } if(fClusters[ic+kNtb]) dq += TMath::Abs(fClusters[ic+kNtb]->GetQ()); if(dq<1.e-3) return 0.; + Double_t dx = fdX; if(ic-1>=0 && ic+1 + Bool_t kPRINT = kFALSE; if(!fReconstructor->GetRecoParam() ){ AliError("Seed can not be used without a valid RecoParam."); @@ -1007,12 +1016,10 @@ Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt) if(!fClusters[it]) continue; x[irp] = fClusters[it]->GetX(); tb[irp] = fClusters[it]->GetLocalTimeBin(); - printf(" x[%d]=%f t[%d]=%d\n", irp, x[irp], irp, tb[irp]); irp++; } Int_t dtb = tb[1] - tb[0]; fdX = dtb ? (x[0] - x[1]) / dtb : 0.15; - return kTRUE; } @@ -1063,7 +1070,7 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr) // - tilt : switch for tilt pad correction of cluster y position based on // the z, dzdx info from outside [default false]. // - zcorr : switch for using z information to correct for anisochronity -// and a finner error parametrization estimation [default false] +// and a finner error parameterization estimation [default false] // Output : // True if successful // @@ -1071,27 +1078,87 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr) // // Fit in the xy plane // +// The fit is performed to estimate the y position of the tracklet and the track +// angle in the bending plane. The clusters are represented in the chamber coordinate +// system (with respect to the anode wire - see AliTRDtrackerV1::FollowBackProlongation() +// on how this is set). The x and y position of the cluster and also their variances +// are known from clusterizer level (see AliTRDcluster::GetXloc(), AliTRDcluster::GetYloc(), +// AliTRDcluster::GetSX() and AliTRDcluster::GetSY()). +// If gaussian approximation is used to calculate y coordinate of the cluster the position +// is recalculated taking into account the track angle. The general formula to calculate the +// error of cluster position in the gaussian approximation taking into account diffusion and track +// inclination is given for TRD by: +// BEGIN_LATEX +// #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}} + #frac{x^{2}tg^{2}(#phi-#alpha_{L})tg^{2}(#alpha_{L})}{12} +// END_LATEX +// +// Since errors are calculated only in the y directions, radial errors (x direction) are mapped to y +// by projection i.e. +// BEGIN_LATEX +// #sigma_{x|y} = tg(#phi) #sigma_{x} +// END_LATEX +// and also by the lorentz angle correction +// +// Fit in the xz plane +// +// The "fit" is performed to estimate the radial position (x direction) where pad row cross happens. +// If no pad row crossing the z position is taken from geometry and radial position is taken from the xy +// fit (see below). +// +// There are two methods to estimate the radial position of the pad row cross: +// 1. leading cluster radial position : Here the lower part of the tracklet is considered and the last +// cluster registered (at radial x0) on this segment is chosen to mark the pad row crossing. The error +// of the z estimate is given by : +// BEGIN_LATEX +// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12 +// END_LATEX +// The systematic errors for this estimation are generated by the following sources: +// - no charge sharing between pad rows is considered (sharp cross) +// - missing cluster at row cross (noise peak-up, under-threshold signal etc.). +// +// 2. charge fit over the crossing point : Here the full energy deposit along the tracklet is considered +// to estimate the position of the crossing by a fit in the qx plane. The errors in the q directions are +// parameterized as s_q = q^2. The systematic errors for this estimation are generated by the following sources: +// - no general model for the qx dependence +// - physical fluctuations of the charge deposit +// - gain calibration dependence +// +// Estimation of the radial position of the tracklet // +// For pad row cross the radial position is taken from the xz fit (see above). Otherwise it is taken as the +// interpolation point of the tracklet i.e. the point where the error in y of the fit is minimum. The error +// in the y direction of the tracklet is (see AliTRDseedV1::GetCovAt()): +// BEGIN_LATEX +// #sigma_{y} = #sigma^{2}_{y_{0}} + 2xcov(y_{0}, dy/dx) + #sigma^{2}_{dy/dx} +// END_LATEX +// and thus the radial position is: +// BEGIN_LATEX +// x = - cov(y_{0}, dy/dx)/#sigma^{2}_{dy/dx} +// END_LATEX +// +// Estimation of tracklet position error +// +// The error in y direction is the error of the linear fit at the radial position of the tracklet while in the z +// direction is given by the cluster error or pad row cross error. In case of no pad row cross this is given by: +// BEGIN_LATEX +// #sigma_{y} = #sigma^{2}_{y_{0}} - 2cov^{2}(y_{0}, dy/dx)/#sigma^{2}_{dy/dx} + #sigma^{2}_{dy/dx} +// #sigma_{z} = Pad_{length}/12 +// END_LATEX +// For pad row cross the full error is calculated at the radial position of the crossing (see above) and the error +// in z by the width of the crossing region - being a matter of parameterization. +// BEGIN_LATEX +// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12 +// END_LATEX +// In case of no tilt correction (default in the barrel tracking) the tilt is taken into account by the rotation of +// the covariance matrix. See AliTRDseedV1::GetCovAt() for details. +// +// Author +// A.Bercuci if(!IsCalibrated()) Calibrate(); const Int_t kClmin = 8; - - // cluster error parametrization parameters - // 1. sy total charge - const Float_t sq0inv = 0.019962; // [1/q0] - const Float_t sqb = 1.0281564; //[cm] - // 2. sy for the PRF - const Float_t scy[AliTRDgeometry::kNlayer][4] = { - {2.827e-02, 9.600e-04, 4.296e-01, 2.271e-02}, - {2.952e-02,-2.198e-04, 4.146e-01, 2.339e-02}, - {3.090e-02, 1.514e-03, 4.020e-01, 2.402e-02}, - {3.260e-02,-2.037e-03, 3.946e-01, 2.509e-02}, - {3.439e-02,-3.601e-04, 3.883e-01, 2.623e-02}, - {3.510e-02, 2.066e-03, 3.651e-01, 2.588e-02}, - }; - // get track direction Double_t y0 = fYref[0]; Double_t dydx = fYref[1]; @@ -1099,10 +1166,6 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr) Double_t dzdx = fZref[1]; Double_t yt, zt; - // calculation of tg^2(phi - a_L) and tg^2(a_L) - Double_t tgg = (dydx-fExB)/(1.+dydx*fExB); tgg *= tgg; - //Double_t exb2= fExB*fExB; - //AliTRDtrackerV1::AliTRDLeastSquare fitterZ; TLinearFitter fitterY(1, "pol1"); TLinearFitter fitterZ(1, "pol1"); @@ -1110,11 +1173,9 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr) // book cluster information Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters]; - Int_t ily = AliTRDgeometry::GetLayer(fDet); Int_t n = 0; AliTRDcluster *c=0x0, **jc = &fClusters[0]; for (Int_t ic=0; icGetNPads()>4) w = .5; if(c->GetNPads()>5) w = .2; - Int_t tb = c->GetLocalTimeBin(); + // cluster charge qc[n] = TMath::Abs(c->GetQ()); + // pad row of leading + // Radial cluster position //Int_t jc = TMath::Max(fN-3, 0); //xc[fN] = c->GetXloc(fT0, fVD, &qc[jc], &xc[jc]/*, z0 - c->GetX()*dzdx*/); xc[n] = fX0 - c->GetX(); - //Double_t s2 = fS2PRF + fDiffL*fDiffL*xc[n]/(1.+2.*exb2)+tgg*xc[n]*xc[n]*exb2/12.; - //yc[fN] = c->GetYloc(s2, GetPadWidth(), xc[fN], fExB); - yc[n] = c->GetY()-AliTRDcluster::GetYcorr(ily, c->GetCenter()); - zc[n] = c->GetZ(); - - // extrapolated y value for the track + // extrapolated track to cluster position yt = y0 - xc[n]*dydx; - // extrapolated z value for the track zt = z0 - xc[n]*dzdx; - // tilt correction - if(tilt) yc[n] -= GetTilt()*(zc[n] - zt); - - // ELABORATE CLUSTER ERROR - // basic y error (|| to track). - sy[n] = AliTRDcluster::GetSY(tb, zcorr?zt:-1.); - //printf("cluster[%d]\n\tsy[0] = %5.3e [um]\n", fN, sy[fN]*1.e4); - // y error due to total charge - sy[n] += sqb*(1./qc[n] - sq0inv); - //printf("\tsy[1] = %5.3e [um]\n", sy[fN]*1.e4); - // y error due to PRF - sy[n] += scy[ily][0]*TMath::Gaus(c->GetCenter(), scy[ily][1], scy[ily][2]) - scy[ily][3]; - //printf("\tsy[2] = %5.3e [um]\n", sy[fN]*1.e4); - - sy[n] *= sy[n]; - - // ADD ERROR ON x - // error of drift length parallel to the track - Double_t sx = AliTRDcluster::GetSX(tb, zcorr?zt:-1.); // [cm] - //printf("\tsx[0] = %5.3e [um]\n", sx*1.e4); - sx *= sx; // square sx - - // add error from ExB - sy[n] += fExB*fExB*sx; - //printf("\tsy[3] = %5.3e [um^2]\n", sy[fN]*1.e8); - - // global radial error due to misalignment/miscalibration - Double_t sx0 = 0.; sx0 *= sx0; - // add sx contribution to sy due to track angle - sy[n] += tgg*(sx+sx0); - // TODO we should add tilt pad correction here - //printf("\tsy[4] = %5.3e [um^2]\n", sy[fN]*1.e8); - c->SetSigmaY2(sy[n]); - - sy[n] = TMath::Sqrt(sy[n]); - fitterY.AddPoint(&xc[n], yc[n], sy[n]); + + // Recalculate cluster error based on tracking information + c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], zcorr?zt:-1., dydx); + sy[n] = TMath::Sqrt(c->GetSigmaY2()); + + yc[n] = fReconstructor->UseGAUS() ? + c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY(); + zc[n] = c->GetZ(); + //optional tilt correction + if(tilt) yc[n] -= (GetTilt()*(zc[n] - zt)); + + fitterY.AddPoint(&xc[n], yc[n], TMath::Sqrt(sy[n])); fitterZ.AddPoint(&xc[n], qc[n], 1.); n++; } @@ -1199,7 +1232,24 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr) // fit XZ if(IsRowCross()){ + // THE LEADING CLUSTER METHOD + Float_t xMin = fX0; Int_t ic=n=kNclusters-1; jc = &fClusters[ic]; + AliTRDcluster *c0 =0x0, **kc = &fClusters[kNtb-1]; + for(; ic>kNtb; ic--, --jc, --kc){ + if((c0 = (*kc)) && c0->IsInChamber() && (xMin>c0->GetX())) xMin = c0->GetX(); + if(!(c = (*jc))) continue; + if(!c->IsInChamber()) continue; + zc[kNclusters-1] = c->GetZ(); + fX = fX0 - c->GetX(); + } + fZfit[0] = .5*(zc[0]+zc[kNclusters-1]); fZfit[1] = 0.; + // Error parameterization + fS2Z = fdX*fZref[1]; + fS2Z *= fS2Z; fS2Z *= 0.2887; // 1/sqrt(12) +/* + // THE FIT X-Q PLANE METHOD + ic=n=kNclusters-1; jc = &fClusters[ic]; for(; ic>kNtb; ic--, --jc){ if(!(c = (*jc))) continue; if(!c->IsInChamber()) continue; @@ -1224,51 +1274,13 @@ Bool_t AliTRDseedV1::Fit(Bool_t tilt, Bool_t zcorr) fS2Z = 0.05+0.4*TMath::Abs(fZref[1]); fS2Z *= fS2Z; // TODO correct formula //fS2Z = sigma_x*TMath::Abs(fZref[1]); +*/ } else { fZfit[0] = zc[0]; fZfit[1] = 0.; fS2Z = GetPadLength()*GetPadLength()/12.; } fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2]; return kTRUE; -// // determine z offset of the fit -// Float_t zslope = 0.; -// Int_t nchanges = 0, nCross = 0; -// if(nz==2){ // tracklet is crossing pad row -// // Find the break time allowing one chage on pad-rows -// // with maximal number of accepted clusters -// Int_t padRef = zRow[0]; -// for (Int_t ic=1; ic zc[ic] ? 1. : -1.; -// padRef = zRow[ic]; -// nCross = ic; -// nchanges++; -// } -// } -// -// // condition on nCross and reset nchanges TODO -// -// if(nchanges==1){ -// if(dzdx * zslope < 0.){ -// AliInfo("Tracklet-Track mismatch in dzdx. TODO."); -// } -// -// -// //zc[nc] = fitterZ.GetFunctionParameter(0); -// fCross[1] = fYfit[0] - fCross[0] * fYfit[1]; -// fCross[0] = fX0 - fCross[0]; -// } } diff --git a/TRD/AliTRDseedV1.h b/TRD/AliTRDseedV1.h index b7d3d127ef5..6429aa18197 100644 --- a/TRD/AliTRDseedV1.h +++ b/TRD/AliTRDseedV1.h @@ -59,7 +59,7 @@ public: ,kPID = BIT(16) // PID contributor ,kCalib = BIT(17) // calibrated tracklet ,kKink = BIT(18) // kink prolongation tracklet - ,kStandAlone = BIT(19) // stand alone build tracklet + ,kStandAlone = BIT(19) // tracklet build during stand alone track finding }; AliTRDseedV1(Int_t det = -1); @@ -67,11 +67,7 @@ public: AliTRDseedV1(const AliTRDseedV1 &ref); AliTRDseedV1& operator=(const AliTRDseedV1 &ref); -/* Bool_t AttachClustersIter( - AliTRDtrackingChamber *chamber, Float_t quality, - Bool_t kZcorr = kFALSE, AliTRDcluster *c=0x0);*/ - Bool_t AttachClusters( - AliTRDtrackingChamber *chamber, Bool_t tilt = kFALSE); + Bool_t AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt = kFALSE); void Bootstrap(const AliTRDReconstructor *rec); void Calibrate(); void CookdEdx(Int_t nslices); @@ -129,7 +125,7 @@ public: Float_t GetS2Z() const { return fS2Z;} Float_t GetSigmaY() const { return fS2Y > 0. ? TMath::Sqrt(fS2Y) : 0.2;} Float_t GetSnp() const { return fYref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);} - Float_t GetTgl() const { return fZref[1];} + Float_t GetTgl() const { return fZref[1]/TMath::Sqrt(1+fYref[1]*fYref[1]);} Float_t GetTilt() const { return fPad[2];} UInt_t GetTrackletWord() const { return 0;} UShort_t GetVolumeId() const; @@ -137,14 +133,14 @@ public: Float_t GetX() const { return fX0 - fX;} Float_t GetY() const { return fYfit[0] - fYfit[1] * fX;} Double_t GetYat(Double_t x) const { return fYfit[0] - fYfit[1] * (fX0-x);} - Float_t GetYfit(Int_t id) const { return fYfit[id];} - Float_t GetYref(Int_t id) const { return fYref[id];} - Float_t GetZ() const { return fZfit[0] - fZfit[1] * fX;} + Float_t GetYfit(Int_t id) const { return fYfit[id];} + Float_t GetYref(Int_t id) const { return fYref[id];} + Float_t GetZ() const { return fZfit[0] - fZfit[1] * fX;} Double_t GetZat(Double_t x) const { return fZfit[0] - fZfit[1] * (fX0-x);} - Float_t GetZfit(Int_t id) const { return fZfit[id];} - Float_t GetZref(Int_t id) const { return fZref[id];} - Int_t GetYbin() const { return Int_t(GetY()/0.016);} - Int_t GetZbin() const { return Int_t(GetZ()/fPad[0]);} + Float_t GetZfit(Int_t id) const { return fZfit[id];} + Float_t GetZref(Int_t id) const { return fZref[id];} + Int_t GetYbin() const { return Int_t(GetY()/0.016);} + Int_t GetZbin() const { return Int_t(GetZ()/fPad[0]);} inline AliTRDcluster* NextCluster(); inline AliTRDcluster* PrevCluster(); @@ -152,8 +148,8 @@ public: inline void ResetClusterIter(Bool_t forward = kTRUE); void Reset(); - void SetC(Float_t c) { fC = c;} - void SetChi2(Float_t chi2) { fChi2 = chi2;} + void SetC(Float_t c) { fC = c;} + void SetChi2(Float_t chi2) { fChi2 = chi2;} inline void SetCovRef(const Double_t *cov); void SetIndexes(Int_t i, Int_t idx) { fIndexes[i] = idx; } void SetLabels(Int_t *lbls) { memcpy(fLabels, lbls, 3*sizeof(Int_t)); } @@ -199,10 +195,10 @@ private: Short_t fDet; // TRD detector AliTRDcluster *fClusters[kNclusters]; // Clusters Float_t fPad[3]; // local pad definition : length/width/tilt - Float_t fYref[2]; // Reference y - Float_t fZref[2]; // Reference z - Float_t fYfit[2]; // Y fit position +derivation - Float_t fZfit[2]; // Z fit position + Float_t fYref[2]; // Reference y, dydx + Float_t fZref[2]; // Reference z, dz/dx + Float_t fYfit[2]; // Fit y, dy/dx + Float_t fZfit[2]; // Fit z Float_t fPt; // Pt estimate @ tracklet [GeV/c] Float_t fdX; // length of time bin Float_t fX0; // anode wire position @@ -214,7 +210,7 @@ private: Float_t fC; // Curvature Float_t fChi2; // Global chi2 Float_t fdEdx[kNslices]; // dE/dx measurements for tracklet - Float_t fProb[AliPID::kSPECIES]; // PID probabilities + Float_t fProb[AliPID::kSPECIES]; // PID probabilities Int_t fLabels[3]; // most frequent MC labels and total number of different labels Double_t fRefCov[7]; // covariance matrix of the track in the yz plane + the rest of the diagonal elements Double_t fCov[3]; // covariance matrix of the tracklet in the xy plane diff --git a/TRD/AliTRDtracker.cxx b/TRD/AliTRDtracker.cxx index 2ff6df2e365..918f316cfb0 100644 --- a/TRD/AliTRDtracker.cxx +++ b/TRD/AliTRDtracker.cxx @@ -3555,7 +3555,7 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1 } // Set cluster error - cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); + ((AliCluster*)cl[best[bestiter][it]][it])->SetSigmaY2(expectederr); if (!cl[best[bestiter][it]][it]->IsUsed()) { cl[best[bestiter][it]][it]->SetY(cl[best[bestiter][it]][it]->GetY()); } diff --git a/TRD/AliTRDtransform.cxx b/TRD/AliTRDtransform.cxx index ca00a70919d..f5d12ba2573 100644 --- a/TRD/AliTRDtransform.cxx +++ b/TRD/AliTRDtransform.cxx @@ -234,16 +234,16 @@ Bool_t AliTRDtransform::Transform(AliTRDcluster *c) Float_t x = c->GetXloc(t0, vd); - // pad response width with diffusion corrections - Double_t s2 = fCalPRFROC->GetValue(col, row); s2 *= s2; - Float_t dl, dt; - AliTRDCommonParam::Instance()->GetDiffCoeff(dl, dt, vd); - s2 += dl*dl*x/(1.+2.*exb*exb); - s2 -= - 1.5e-1; - // Pad dimensions Double_t rs = fPadPlane->GetRowSize(row); Double_t cs = fPadPlane->GetColSize(col); + + // cluster error with diffusion corrections + Double_t s2 = cs*fCalPRFROC->GetValue(col, row); + s2 *= s2; + Float_t dl, dt; + AliTRDCommonParam::Instance()->GetDiffCoeff(dl, dt, vd); + Double_t y0 = fPadPlane->GetColPos(col) + .5*cs; Double_t loc[] = { kX0shift-x, // Invert the X-position, @@ -257,7 +257,7 @@ Bool_t AliTRDtransform::Transform(AliTRDcluster *c) // store tracking values c->SetX(trk[0]);c->SetY(trk[1]);c->SetZ(trk[2]); - c->SetSigmaY2(s2); + c->SetSigmaY2(s2, dt, exb, x); c->SetSigmaZ2(fPadPlane->GetRowSize(row)*fPadPlane->GetRowSize(row)/12.); return kTRUE; diff --git a/TRD/qaRec/AliTRDcheckDET.cxx b/TRD/qaRec/AliTRDcheckDET.cxx index 8caef1f4a63..3789e55b525 100644 --- a/TRD/qaRec/AliTRDcheckDET.cxx +++ b/TRD/qaRec/AliTRDcheckDET.cxx @@ -183,6 +183,8 @@ Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){ // // Setting Reference Figures // + gPad->SetLogy(0); + gPad->SetLogx(0); switch(ifig){ case kNclustersTrack: ((TH1F*)fContainer->FindObject("hNcls"))->Draw("pl"); @@ -213,6 +215,7 @@ Bool_t AliTRDcheckDET::GetRefFigure(Int_t ifig){ return kTRUE; case kChargeCluster: ((TH1F*)fContainer->FindObject("hQcl"))->Draw("c"); + gPad->SetLogy(1); return kTRUE; case kChargeTracklet: ((TH1F*)fContainer->FindObject("hQtrklt"))->Draw("c"); @@ -904,17 +907,25 @@ void AliTRDcheckDET::MakePlotChi2() TH2S *h2 = (TH2S*)fContainer->At(kChi2); TF1 f("fChi2", "[0]*pow(x, [1]-1)*exp(-0.5*x)", 0., 50.); - + TLegend *leg = new TLegend(.7,.7,.95,.95); + leg->SetBorderSize(1); leg->SetHeader("Tracklets per Track"); Bool_t kFIRST = kTRUE; for(Int_t il=1; il<=h2->GetNbinsX(); il++){ - TH1D *h1 = h2->ProjectionY(Form("py%d", il), il, il); + TH1D *h1 = h2->ProjectionY(Form("pyChi2%d", il), il, il); if(h1->Integral()<50) continue; h1->Scale(1./h1->Integral()); - h1->SetMarkerStyle(7);h1->SetMarkerColor(il);h1->SetLineColor(il); + h1->SetMarkerStyle(7);h1->SetMarkerColor(il); + h1->SetLineColor(il);h1->SetLineStyle(2); f.SetParameter(1, .5*il);f.SetLineColor(il); - h1->Fit(&f, "Q+", kFIRST ? "e1": "e1 same"); + h1->Fit(&f, "QW+", kFIRST ? "pc": "pcsame"); + leg->AddEntry(h1, Form("%d", il), "l"); + if(kFIRST){ + h1->GetXaxis()->SetRangeUser(0., 25.); + } kFIRST = kFALSE; } + leg->Draw(); + gPad->SetLogy(); } @@ -926,41 +937,37 @@ void AliTRDcheckDET::MakePlotNTracklets(){ TH1F *hBAR = (TH1F *)fContainer->FindObject("hNtlsBAR"); TH1F *hSTA = (TH1F *)fContainer->FindObject("hNtlsSTA"); TH1F *hCON = (TH1F *)fContainer->FindObject("hNtls"); + TLegend *leg = new TLegend(0.6, 0.75, 0.89, 0.89); + + Float_t scale = hCON->Integral(); + hCON->Scale(100./scale); + hCON->SetFillColor(kRed);hCON->SetLineColor(kRed); + hCON->SetBarWidth(0.2); + hCON->SetBarOffset(0.6); + hCON->SetStats(kFALSE); + hCON->GetYaxis()->SetRangeUser(0.,40.); + hCON->GetYaxis()->SetTitleOffset(1.2); + hCON->Draw("bar1"); leg->AddEntry(hCON, "Total", "f"); - hBAR->Scale(100./hCON->Integral()); - hBAR->SetFillColor(kRed); + hBAR->Scale(100./scale); + hBAR->SetFillColor(kGreen);hBAR->SetLineColor(kGreen); hBAR->SetBarWidth(0.2); hBAR->SetBarOffset(0.2); hBAR->SetTitle(""); hBAR->SetStats(kFALSE); hBAR->GetYaxis()->SetRangeUser(0.,40.); hBAR->GetYaxis()->SetTitleOffset(1.2); - hBAR->Draw("bar1"); + hBAR->Draw("bar1same"); leg->AddEntry(hBAR, "Barrel", "f"); - hSTA->Scale(100./hCON->Integral()); - hSTA->SetFillColor(kBlue); + hSTA->Scale(100./scale); + hSTA->SetFillColor(kBlue);hSTA->SetLineColor(kBlue); hSTA->SetBarWidth(0.2); hSTA->SetBarOffset(0.4); hSTA->SetTitle(""); hSTA->SetStats(kFALSE); hSTA->GetYaxis()->SetRangeUser(0.,40.); hSTA->GetYaxis()->SetTitleOffset(1.2); - hSTA->Draw("bar1same"); - - hCON->Scale(100./hCON->Integral()); - hCON->SetFillColor(kGreen); - hCON->SetBarWidth(0.2); - hCON->SetBarOffset(0.6); - hCON->SetStats(kFALSE); - hCON->GetYaxis()->SetRangeUser(0.,40.); - hCON->GetYaxis()->SetTitleOffset(1.2); - hCON->Draw("bar1same"); - - TLegend *leg = new TLegend(0.6, 0.75, 0.89, 0.89); - leg->AddEntry(hBAR, "Barrel", "f"); - leg->AddEntry(hSTA, "Stand Alone", "f"); - leg->AddEntry(hCON, "Convoluted", "f"); - + hSTA->Draw("bar1same"); leg->AddEntry(hSTA, "Stand Alone", "f"); leg->Draw(); gPad->Update(); } diff --git a/TRD/qaRec/AliTRDcheckPID.cxx b/TRD/qaRec/AliTRDcheckPID.cxx index 8a2c1f01335..dca4d304d07 100644 --- a/TRD/qaRec/AliTRDcheckPID.cxx +++ b/TRD/qaRec/AliTRDcheckPID.cxx @@ -735,10 +735,10 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig) TList *content = 0x0; switch(ifig){ case kEfficiency:{ - gPad->Divide(2, 1, 0.,0.); + gPad->Divide(2, 1, 1.e-5, 1.e-5); TList *l=gPad->GetListOfPrimitives(); TVirtualPad *pad = ((TVirtualPad*)l->At(0));pad->cd(); - pad->SetMargin(0.07, 0.07, 0.1, 0.001); + pad->SetMargin(0.1, 0.01, 0.1, 0.01); TLegend *legEff = new TLegend(.7, .7, .98, .98); legEff->SetBorderSize(1); @@ -768,7 +768,7 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig) gPad->SetGridx(); pad = ((TVirtualPad*)l->At(1));pad->cd(); - pad->SetMargin(0.07, 0.07, 0.1, 0.001); + pad->SetMargin(0.1, 0.01, 0.1, 0.01); content = (TList *)fGraph->FindObject("Thresholds"); if(!(g = (TGraphErrors*)content->At(AliTRDpidUtil::kLQ))) break; if(!g->GetN()) break; @@ -796,6 +796,7 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig) FIRST = kTRUE; if(!(h2 = (TH2F*)(fContainer->At(kdEdx)))) break; legdEdx->SetHeader("Particle Species"); + gPad->SetMargin(0.1, 0.01, 0.1, 0.01); for(Int_t is = AliPID::kSPECIES-1; is>=0; is--){ Int_t bin = FindBin(is, 2.); h1 = h2->ProjectionY("px", bin, bin); @@ -827,6 +828,7 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig) FIRST = kTRUE; if(!(h2 = (TH2F*)(fContainer->At(kPH)))) break;; legPH->SetHeader("Particle Species"); + gPad->SetMargin(0.1, 0.01, 0.1, 0.01); for(Int_t is=0; isProjectionY("py", bin, bin); @@ -835,8 +837,8 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig) h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is)); h1->SetLineColor(AliTRDCalPID::GetPartColor(is)); if(FIRST){ - h1->GetXaxis()->SetTitle("tb(1/100 ns^{-1})"); - h1->GetYaxis()->SetTitle(" (a.u.)"); + h1->GetXaxis()->SetTitle("t_{drift} [100*ns]"); + h1->GetYaxis()->SetTitle(" [a.u./cm]"); } h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec"); legPH->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "pl"); @@ -859,21 +861,21 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig) legNClus->SetHeader("Particle Species"); for(Int_t is=0; isProjectionY("py", bin, bin); + h1 = h2->ProjectionY("pyNClus", bin, bin); if(!h1->GetEntries()) continue; - h1->Scale(1./h1->Integral()); + h1->Scale(100./h1->Integral()); //h1->SetMarkerStyle(24); //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is)); h1->SetLineColor(AliTRDCalPID::GetPartColor(is)); if(FIRST) h1->GetXaxis()->SetTitle("N^{cl}/tracklet"); - if(FIRST) h1->GetYaxis()->SetTitle(""); + if(FIRST) h1->GetYaxis()->SetTitle("Prob. [%]"); h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec"); legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l"); FIRST = kFALSE; } if(FIRST) break; legNClus->Draw(); - gPad->SetLogy(); + gPad->SetLogy(0); gPad->SetLogx(0); gPad->SetGridy(); gPad->SetGridx(); @@ -890,21 +892,24 @@ Bool_t AliTRDcheckPID::GetRefFigure(Int_t ifig) legNClus->SetHeader("Particle Species"); for(Int_t is=0; isProjectionY("py", bin, bin); + h1 = h2->ProjectionY("pyNTracklets", bin, bin); if(!h1->GetEntries()) continue; - h1->Scale(1./h1->Integral()); + h1->Scale(100./h1->Integral()); //h1->SetMarkerStyle(24); //h1->SetMarkerColor(AliTRDCalPID::GetPartColor(is)); h1->SetLineColor(AliTRDCalPID::GetPartColor(is)); - if(FIRST) h1->GetXaxis()->SetTitle("N^{tl}/track"); - if(FIRST) h1->GetYaxis()->SetTitle(""); + if(FIRST){ + h1->GetXaxis()->SetTitle("N^{trklt}/track"); + h1->GetXaxis()->SetRangeUser(1.,6.); + h1->GetYaxis()->SetTitle("Prob. [%]"); + } h = (TH1F*)h1->DrawClone(FIRST ? "c" : "samec"); legNClus->AddEntry(h, Form("%s", AliTRDCalPID::GetPartName(is)), "l"); FIRST = kFALSE; } if(FIRST) break; legNClus->Draw(); - gPad->SetLogy(); + gPad->SetLogy(0); gPad->SetLogx(0); gPad->SetGridy(); gPad->SetGridx(); diff --git a/TRD/qaRec/AliTRDclusterResolution.cxx b/TRD/qaRec/AliTRDclusterResolution.cxx index 0ff75f02447..ef264300f92 100644 --- a/TRD/qaRec/AliTRDclusterResolution.cxx +++ b/TRD/qaRec/AliTRDclusterResolution.cxx @@ -183,6 +183,7 @@ #include "TObjArray.h" #include "TAxis.h" #include "TF1.h" +#include "TLegend.h" #include "TGraphErrors.h" #include "TLine.h" #include "TH2I.h" @@ -215,6 +216,7 @@ AliTRDclusterResolution::AliTRDclusterResolution(const char *name, const char *t ,fZ(0.) { memset(fR, 0, 4*sizeof(Float_t)); + memset(fP, 0, 4*sizeof(Float_t)); // time drift axis fAt = new TAxis(kNTB, 0., kNTB*fgkTimeBinLength); @@ -254,7 +256,7 @@ void AliTRDclusterResolution::CreateOutputObjects() Bool_t AliTRDclusterResolution::GetRefFigure(Int_t ifig) { if(!fResults) return kFALSE; - + TLegend *leg = 0x0; TList *l = 0x0; TObjArray *arr = 0x0; TH2 *h2 = 0x0;TH1 *h1 = 0x0; @@ -266,20 +268,31 @@ Bool_t AliTRDclusterResolution::GetRefFigure(Int_t ifig) if(!(gs = (TGraphErrors*)arr->At(1))) break; if(!(gp = (TGraphErrors*)arr->At(2))) break; gs->Draw("apl"); - gs->GetHistogram()->GetYaxis()->SetRangeUser(-.01, .6); + gs->GetHistogram()->GetYaxis()->SetRangeUser(-50., 700.); gs->GetHistogram()->SetXTitle("Q [a.u.]"); - gs->GetHistogram()->SetYTitle("#sigma_{y} / #mu_{y} [mm] / freq"); + gs->GetHistogram()->SetYTitle("#sigma_{y} / #mu_{y} [#mum] / freq"); gm->Draw("pl"); gp->Draw("pl"); return kTRUE; case kCenter: if(!(arr = (TObjArray*)fResults->At(kCenter))) break; - gPad->Divide(3, 1); l = gPad->GetListOfPrimitives(); - for(Int_t ipad = 3; ipad--;){ - if(!(h2 = (TH2F*)arr->At(ipad))) return kFALSE; - ((TVirtualPad*)l->At(ipad))->cd(); - h2->Draw("lego2fb"); + gPad->Divide(2, 1); l = gPad->GetListOfPrimitives(); + ((TVirtualPad*)l->At(0))->cd(); + ((TTree*)arr->At(0))->Draw("y:x>>h(23, 0.1, 2.4, 51, -.51, .51)", + "m[0]*(ly==0&&abs(m[0])<1.e-1)", "colz"); + ((TVirtualPad*)l->At(1))->cd(); + leg= new TLegend(.7, .7, .9, .95); + leg->SetBorderSize(0); leg->SetFillColor(0); leg->SetFillStyle(0); + leg->SetHeader("TRD Plane"); + for(Int_t il = 1; il<=AliTRDgeometry::kNlayer; il++){ + if(!(gm = (TGraphErrors*)arr->At(il))) return kFALSE; + gm->Draw(il>1?"pc":"apc"); leg->AddEntry(gm, Form("%d", il-1), "pl"); + if(il>1) continue; + gm->GetHistogram()->SetXTitle("t_{drift} [#mus]"); + gm->GetHistogram()->SetYTitle("#sigma_{y}(x|cen=0) [#mum]"); + gm->GetHistogram()->GetYaxis()->SetRangeUser(150., 500.); } + leg->Draw(); return kTRUE; case kSigm: if(!(arr = (TObjArray*)fResults->At(kSigm))) break; @@ -326,21 +339,16 @@ TObjArray* AliTRDclusterResolution::Histos() fContainer = new TObjArray(kNtasks); //fContainer->SetOwner(kTRUE); - TH2I *h2 = 0x0; TH3S *h3 = 0x0; TObjArray *arr = 0x0; - fContainer->AddAt(h2 = new TH2I("Charge", "dy=f(q)", 50, 2.2, 7.5, 60, -.3, .3), kQRes); - h2->SetXTitle("log(q) [a.u.]"); - h2->SetYTitle("#Delta y[cm]"); - h2->SetZTitle("entries"); - - fContainer->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer), kCenter); + fContainer->AddAt(arr = new TObjArray(2*AliTRDgeometry::kNlayer), kCenter); arr->SetName("Center"); for(Int_t il=0; ilFindObject(Form("hc_l%1d", il)))){ + // add resolution plot for each layer + if(!(h3=(TH3S*)gROOT->FindObject(Form("hCenResLy%d", il)))){ h3 = new TH3S( - Form("hc_l%1d", il), + Form("hCenResLy%d", il), Form(" ly [%d]", il), kNTB, fAt->GetBinLowEdge(1), fAt->GetBinUpEdge(kNTB), // x 51, -.51, .51, // y @@ -350,8 +358,29 @@ TObjArray* AliTRDclusterResolution::Histos() h3->SetZTitle("#Delta y[cm]"); } h3->Reset(); arr->AddAt(h3, il); + // add Pull plot for each layer + if(!(h3=(TH3S*)gROOT->FindObject(Form("hCenPullLy%d", il)))){ + h3 = new TH3S( + Form("hCenPullLy%d", il), + Form(" ly [%d]", il), + kNTB, fAt->GetBinLowEdge(1), fAt->GetBinUpEdge(kNTB), // x + 51, -.51, .51, // y + 60, -4., 4.); // dy + h3->SetXTitle("x [#mus]"); + h3->SetYTitle("y [pw]"); + h3->SetZTitle("#Delta y/#sigma_{y}"); + } h3->Reset(); + arr->AddAt(h3, AliTRDgeometry::kNlayer+il); } + if(!(h3 = (TH3S*)gROOT->FindObject("Charge"))){ + h3 = new TH3S("Charge", "dy=f(q)", 50, 2.2, 7.5, 60, -.3, .3, 60, -4., 4.); + h3->SetXTitle("log(q) [a.u.]"); + h3->SetYTitle("#Delta y[cm]"); + h3->SetZTitle("#Delta y/#sigma_{y}"); + } + fContainer->AddAt(h3, kQRes); + fContainer->AddAt(arr = new TObjArray(kNTB), kSigm); arr->SetName("Resolution"); for(Int_t ix=0; ixAt(kQRes); - h2->Fill(TMath::Log(q), dy); + h3 = (TH3S*)fContainer->At(kQRes); + h3->Fill(TMath::Log(q), dy, dy/TMath::Sqrt(covcl[0])); } // do not use problematic clusters in resolution analysis @@ -429,12 +457,14 @@ void AliTRDclusterResolution::Exec(Option_t *) x = (t+.5)*fgkTimeBinLength; // conservative approach !! // resolution as a function of y displacement from pad center - // only for phi equal exB and clusters close to cathode wire plane - // for ideal simulations time bins 4,5 and 6. - if(TMath::Abs(dydx-fExB) < kDtgPhi && - TMath::Abs(x-0.675)<0.225){ - h3 = (TH3S*)arr0->At(AliTRDgeometry::GetLayer(det)); + // only for phi equal exB + if(TMath::Abs(dydx-fExB) < kDtgPhi/* && + TMath::Abs(x-0.675)<0.225*/){ + Int_t ly(AliTRDgeometry::GetLayer(det)); + h3 = (TH3S*)arr0->At(ly); h3->Fill(x, cli->GetYDisplacement(), dy); + h3 = (TH3S*)arr0->At(AliTRDgeometry::kNlayer+ly); + h3->Fill(x, cli->GetYDisplacement(), dy/TMath::Sqrt(covcl[0])); } Int_t ix = fAt->FindBin(x); @@ -478,12 +508,22 @@ Bool_t AliTRDclusterResolution::PostProcess() g->SetMarkerStyle(7); // pad center dependence - fResults->AddAt(t = new TTree("cent", "dy=f(y,x,ly)"), kCenter); + fResults->AddAt(arr = new TObjArray(AliTRDgeometry::kNlayer+1), kCenter); + arr->SetOwner(); + arr->AddAt( + t = new TTree("cent", "dy=f(y,x,ly)"), 0); t->Branch("ly", &fLy, "ly/B"); t->Branch("x", &fX, "x/F"); t->Branch("y", &fY, "y/F"); t->Branch("m", &fR[0], "m[2]/F"); t->Branch("s", &fR[2], "s[2]/F"); + t->Branch("pm", &fP[0], "pm[2]/F"); + t->Branch("ps", &fP[2], "ps[2]/F"); + for(Int_t il=1; il<=AliTRDgeometry::kNlayer; il++){ + arr->AddAt(g = new TGraphErrors(), il); + g->SetLineColor(il); g->SetLineStyle(il); + g->SetMarkerColor(il);g->SetMarkerStyle(4); + } fResults->AddAt(t = new TTree("sigm", "dy=f(dw,x,dydx)"), kSigm); @@ -567,7 +607,7 @@ void AliTRDclusterResolution::SetVisual() void AliTRDclusterResolution::ProcessCharge() { TH2I *h2 = 0x0; - if((h2 = (TH2I*)fContainer->At(kQRes))) { + if(!(h2 = (TH2I*)fContainer->At(kQRes))) { AliWarning("Missing dy=f(Q) histo"); return; } @@ -593,12 +633,12 @@ void AliTRDclusterResolution::ProcessCharge() // Fill sy^2 = f(q) Int_t ip = gqm->GetN(); - gqm->SetPoint(ip, q, 10.*f.GetParameter(1)); - gqm->SetPointError(ip, 0., 10.*f.GetParError(1)); + gqm->SetPoint(ip, q, 1.e4*f.GetParameter(1)); + gqm->SetPointError(ip, 0., 1.e4*f.GetParError(1)); // correct sigma for ExB effect - gqs->SetPoint(ip, q, 1.e1*f.GetParameter(2)/**f.GetParameter(2)-exb2*sxd2*/); - gqs->SetPointError(ip, 0., 1.e1*f.GetParError(2)/**f.GetParameter(2)*/); + gqs->SetPoint(ip, q, 1.e4*f.GetParameter(2)/**f.GetParameter(2)-exb2*sxd2*/); + gqs->SetPointError(ip, 0., 1.e4*f.GetParError(2)/**f.GetParameter(2)*/); // save probability n += entries; @@ -611,7 +651,7 @@ void AliTRDclusterResolution::ProcessCharge() for(Int_t ip=gqp->GetN(); ip--;){ gqp->GetPoint(ip, q, entries); entries/=n; - gqp->SetPoint(ip, q, entries); + gqp->SetPoint(ip, q, 1.e3*entries); gqs->GetPoint(ip, q, sy); sm += entries*sy; } @@ -619,61 +659,108 @@ void AliTRDclusterResolution::ProcessCharge() // error parametrization s(q) = + b(1/q-1/q0) TF1 fq("fq", "[0] + [1]/x", 20., 250.); gqs->Fit(&fq); - printf("sy(Q) :: sm[%f] b[%f] 1/q0[%f]\n", sm, fq.GetParameter(1), (sm-fq.GetParameter(0))/fq.GetParameter(1)); + //printf("sm=%f [0]=%f [1]=%f\n", 1.e-4*sm, fq.GetParameter(0), fq.GetParameter(1)); + printf(" const Float_t sq0inv = %f; // [1/q0]\n", (sm-fq.GetParameter(0))/fq.GetParameter(1)); + printf(" const Float_t sqb = %f; // [cm]\n", 1.e-4*fq.GetParameter(1)); } //_______________________________________________________ void AliTRDclusterResolution::ProcessCenterPad() { +// Resolution as a function of y displacement from pad center +// only for phi equal exB and clusters close to cathode wire plane +// for ideal simulations time bins 4,5 and 6. + TObjArray *arr = (TObjArray*)fContainer->At(kCenter); if(!arr) { AliWarning("Missing dy=f(y | x, ly) container"); return; } + Float_t s[AliTRDgeometry::kNlayer]; TF1 f("f", "gaus", -.5, .5); - TH1D *h1 = 0x0; TH3S *h3=0x0; - TTree *t = (TTree*)fResults->At(kCenter); + TF1 fp("fp", "gaus", -3.5, 3.5); + + TH1D *h1 = 0x0; TH2F *h2 = 0x0; TH3S *h3r=0x0, *h3p=0x0; + TObjArray *arrRes = (TObjArray*)fResults->At(kCenter); + TTree *t = (TTree*)arrRes->At(0); + TGraphErrors *gs = 0x0; TAxis *ax = 0x0; - for(Int_t il=0; ilGetEntriesFast(); il++){ - if(!(h3 = (TH3S*)arr->At(il))) continue; + printf(" const Float_t lSy[6][24] = {\n {"); + const Int_t nl = AliTRDgeometry::kNlayer; + for(Int_t il=0; ilAt(il))) continue; + if(!(h3p = (TH3S*)arr->At(nl+il))) continue; + gs = (TGraphErrors*)arrRes->At(il+1); fLy = il; - for(Int_t ix=1; ix<=h3->GetXaxis()->GetNbins(); ix++){ - ax = h3->GetXaxis(); - ax->SetRange(ix, ix); + for(Int_t ix=1; ix<=h3r->GetXaxis()->GetNbins(); ix++){ + ax = h3r->GetXaxis(); ax->SetRange(ix, ix); + ax = h3p->GetXaxis(); ax->SetRange(ix, ix); fX = ax->GetBinCenter(ix); //printf(" x[%2d]=%4.2f\n", ix, fX); - for(Int_t iy=1; iy<=h3->GetYaxis()->GetNbins(); iy++){ - ax = h3->GetYaxis(); - ax->SetRange(iy, iy); + for(Int_t iy=1; iy<=h3r->GetYaxis()->GetNbins(); iy++){ + ax = h3r->GetYaxis(); ax->SetRange(iy, iy); + ax = h3p->GetYaxis(); ax->SetRange(iy, iy); fY = ax->GetBinCenter(iy); //printf(" y[%2d]=%5.2f\n", iy, fY); // finish navigation in the HnSparse - h1 = (TH1D*)h3->Project3D("z"); + h1 = (TH1D*)h3r->Project3D("z"); Int_t entries = (Int_t)h1->Integral(); if(entries < 50) continue; //Adjust(&f, h1); h1->Fit(&f, "QN"); - // Fill sy,my=f(y_w,x,ly) fR[0] = f.GetParameter(1); fR[1] = f.GetParError(1); fR[2] = f.GetParameter(2); fR[3] = f.GetParError(2); + h1 = (TH1D*)h3p->Project3D("z"); + h1->Fit(&fp, "QN"); + fP[0] = fp.GetParameter(1); fP[1] = fp.GetParError(1); + fP[2] = fp.GetParameter(2); fP[3] = fp.GetParError(2); + //printf("ly[%d] x[%3.1f] y[%+5.2f] m[%5.3f] s[%5.3f] \n", fLy, fX, fY, fR[0], fR[2]); t->Fill(); + + } } - if(!fCanvas) continue; + t->Draw("y:x>>h(24, 0., 2.4, 51, -.51, .51)", + Form("s[0]*(ly==%d&&abs(m[0])<1.e-1)", fLy), + "goff"); + h2=(TH2F*)gROOT->FindObject("h"); + f.FixParameter(1, 0.); + Int_t n = h2->GetXaxis()->GetNbins(); s[il]=0.; + printf(" {"); + for(Int_t ix=1; ix<=n; ix++){ + ax = h2->GetXaxis(); + fX = ax->GetBinCenter(ix); + h1 = h2->ProjectionY("hCenPy", ix, ix); + //if((Int_t)h1->Integral() < 1.e-10) continue; + h1->Fit(&f, "QN"); + s[il]+=f.GetParameter(2); + printf("%6.4f,%s", f.GetParameter(0), ix%6?" ":"\n "); + Int_t jx = gs->GetN(); + gs->SetPoint(jx, fX, 1.e4*f.GetParameter(0)); + gs->SetPointError(jx, 0., 0./*f.GetParError(0)*/); + } + printf("\b},\n"); + s[il]/=n; + + f.ReleaseParameter(2); - t->Draw("y:x>>h(23, 0.1, 2.4, 51, -.51, .51)", - Form("m[0]*(ly==%d&&abs(m[0])<1.e-1)", fLy), - "lego2fb"); + + if(!fCanvas) continue; + h2->Draw("lego2fb"); fCanvas->Modified(); fCanvas->Update(); if(IsSaveAs()) fCanvas->SaveAs(Form("Figures/ProcessCenter_ly[%d].gif", fLy)); else gSystem->Sleep(100); } + printf(" };\n"); + printf(" const Float_t lPRF[] = {" + "%5.3f, %5.3f, %5.3f, %5.3f, %5.3f, %5.3f};\n", + s[0], s[1], s[2], s[3], s[4], s[5]); } //_______________________________________________________ diff --git a/TRD/qaRec/AliTRDclusterResolution.h b/TRD/qaRec/AliTRDclusterResolution.h index caaa7e7d3ed..2c923f3ade6 100644 --- a/TRD/qaRec/AliTRDclusterResolution.h +++ b/TRD/qaRec/AliTRDclusterResolution.h @@ -18,11 +18,11 @@ public: ,kN = kND*kNTB }; enum EResultContainer { // results container type - kQRes = 0 - ,kCenter = 1 - ,kSigm = 2 - ,kMean = 3 - ,kNtasks = 4 + kCenter = 0 // cluster2center pad calibration + ,kQRes = 1 // resolution on charge dependence + ,kSigm = 2 // sigma cluster as func of x and z + ,kMean = 3 // shift cluster as func of x and z + ,kNtasks = 4 // total number os sub tasks }; enum ECheckBits { // force setting the ExB kSaveAs = BIT(22) @@ -81,9 +81,10 @@ private: Float_t fX; //! local drift length Float_t fY; //! local rphi offset Float_t fZ; //! local anode wire offset - Float_t fR[4]; //! val/error result list + Float_t fR[4]; //! mean/sgm resolution + Float_t fP[4]; //! mean/sgm pulls - ClassDef(AliTRDclusterResolution, 2) // cluster resolution + ClassDef(AliTRDclusterResolution, 3) // cluster resolution }; //___________________________________________________ diff --git a/TRD/qaRec/AliTRDrecoTask.h b/TRD/qaRec/AliTRDrecoTask.h index 656730b85cb..fa3378a5a7f 100644 --- a/TRD/qaRec/AliTRDrecoTask.h +++ b/TRD/qaRec/AliTRDrecoTask.h @@ -45,7 +45,7 @@ public: Bool_t HasPostProcess() const {return TestBit(kPostProcess);}; virtual TObjArray* Histos() {return fContainer;} - virtual Bool_t Load(const Char_t *filename); + virtual Bool_t Load(const Char_t *filename = "TRD.Performance.root"); virtual Bool_t Save(TObjArray *res); virtual Bool_t PostProcess(); virtual Bool_t PutTrendValue(Char_t *name, Double_t val, Double_t err); diff --git a/TRD/qaRec/AliTRDresolution.cxx b/TRD/qaRec/AliTRDresolution.cxx index f9420244ec6..232344f9ff8 100644 --- a/TRD/qaRec/AliTRDresolution.cxx +++ b/TRD/qaRec/AliTRDresolution.cxx @@ -312,10 +312,11 @@ TH1* AliTRDresolution::PlotCluster(const AliTRDtrackV1 *track) yc -= tilt*(zc-zt); // tilt correction dy = yt - yc; - Float_t sx2 = dydx*c->GetSX(c->GetLocalTimeBin()); sx2*=sx2; - Float_t sy2 = c->GetSY(c->GetLocalTimeBin()); sy2*=sy2; + //Float_t sx2 = dydx*c->GetSX(c->GetLocalTimeBin()); sx2*=sx2; + Float_t sy2 = c->GetSigmaY2(); + if(sy2<=0.) continue; ((TH2I*)arr->At(0))->Fill(dydx, dy); - ((TH2I*)arr->At(1))->Fill(dydx, dy/TMath::Sqrt(cov[0] + sx2 + sy2)); + ((TH2I*)arr->At(1))->Fill(dydx, dy/TMath::Sqrt(cov[0] /*+ sx2*/ + sy2)); if(fDebugLevel>=1){ // Get z-position with respect to anode wire @@ -495,10 +496,11 @@ TH1* AliTRDresolution::PlotTrackTPC(const AliTRDtrackV1 *track) if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) return h; // translate to reference radial position dx = x0 - x; y0 -= dx*dydx0; z0 -= dx*dzdx0; + Float_t norm = 1./TMath::Sqrt(1.+dydx0*dydx0); // 1/sqrt(1+tg^2(phi)) //Fill MC info TVectorD PARMC(kNPAR); PARMC[0]=y0; PARMC[1]=z0; - PARMC[2]=dydx0/TMath::Sqrt(1.+dydx0*dydx0); PARMC[3]=dzdx0; + PARMC[2]=dydx0*norm; PARMC[3]=dzdx0*norm; PARMC[4]=1./pt0; // TMatrixDSymEigen eigen(COV); @@ -525,14 +527,13 @@ TH1* AliTRDresolution::PlotTrackTPC(const AliTRDtrackV1 *track) // pt resolution\\1/pt pulls\\p resolution/pull for(Int_t is=AliPID::kSPECIES; is--;){ if(TMath::Abs(fMC->GetPDG())!=AliPID::ParticleCode(is)) continue; - ((TH3S*)arr->At(8))->Fill(pt0, 1.-PARMC[4]/PAR[4], is); + ((TH3S*)arr->At(8))->Fill(pt0, PARMC[4]/PAR[4]-1., is); ((TH3S*)arr->At(9))->Fill(PARMC[4], (PARMC[4]-PAR[4])/TMath::Sqrt(COV(4,4)), is); - Double_t p0 = TMath::Sqrt(1.+ dzdx0*dzdx0)*pt0, p; + Double_t p0 = TMath::Sqrt(1.+ PARMC[3]*PARMC[3])*pt0, p; Float_t sp; p = tracklet->GetMomentum(&sp); - //printf("p[%f] dp[%f] s[%d]\n", p0, 1.-p/p0, is); - ((TH3S*)arr->At(10))->Fill(p0, 1.-p/p0, is); + ((TH3S*)arr->At(10))->Fill(p0, p/p0-1., is); ((TH3S*)arr->At(11))->Fill(p0, (p0-p)/sp, is); break; } @@ -563,7 +564,7 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track) UChar_t s; Int_t pdg = fMC->GetPDG(), det=-1; Int_t label = fMC->GetLabel(); - Double_t xAnode, x, y, z, pt, dydx, dzdx; + Double_t xAnode, x, y, z, pt, dydx, dzdx, dzdl; Float_t pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0; Double_t covR[7]/*, cov[3]*/; @@ -587,6 +588,7 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track) << "\n"; } + AliTRDReconstructor rec; AliTRDseedV1 *fTracklet = 0x0; for(Int_t ily=0; ilyGetTracklet(ily)))/* || @@ -613,18 +615,19 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track) << "dzdx=" << dzdx0 << "\n"; } - Float_t yt = y0 - dx*dydx0; - Float_t zt = z0 - dx*dzdx0; + Float_t ymc = y0 - dx*dydx0; + Float_t zmc = z0 - dx*dzdx0; //p = pt0*TMath::Sqrt(1.+dzdx0*dzdx0); // pt -> p // Kalman position at reference radial position dx = xAnode - x; - y = fTracklet->GetYref(0) - dx*fTracklet->GetYref(1); - dy = yt - y; - z = fTracklet->GetZref(0) - dx*fTracklet->GetZref(1); - dz = zt - z; dydx = fTracklet->GetYref(1); - dzdx = fTracklet->GetTgl(); + dzdx = fTracklet->GetZref(1); + dzdl = fTracklet->GetTgl(); + y = fTracklet->GetYref(0) - dx*dydx; + dy = y - ymc; + z = fTracklet->GetZref(0) - dx*dzdx; + dz = z - zmc; pt = TMath::Abs(fTracklet->GetPt()); fTracklet->GetCovRef(covR); @@ -638,20 +641,22 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track) // phi resolution/ snp pulls Double_t dtgp = (dydx - dydx0)/(1.- dydx*dydx0); ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dtgp)); - //TODO ((TH2I*)arr->At(5))->Fill(dydx0, ); + Double_t dsnp = dydx/TMath::Sqrt(1.+dydx*dydx) - dydx0/TMath::Sqrt(1.+dydx0*dydx0); + ((TH2I*)arr->At(5))->Fill(dydx0, dsnp/TMath::Sqrt(covR[3])); // theta resolution/ tgl pulls - Double_t dtgl = (dzdx - dzdx0)/(1.- dzdx*dzdx0); - ((TH2I*)arr->At(6))->Fill(dzdx0, + Double_t dzdl0 = dzdx0/TMath::Sqrt(1.+dydx0*dydx0), + dtgl = (dzdl - dzdl0)/(1.- dzdl*dzdl0); + ((TH2I*)arr->At(6))->Fill(dzdl0, TMath::ATan(dtgl)); - //TODO ((TH2I*)arr->At(7))->Fill(dydx0, ); + ((TH2I*)arr->At(7))->Fill(dzdl0, (dzdl - dzdl0)/TMath::Sqrt(covR[4])); // pt resolution \\ 1/pt pulls \\ p resolution for PID for(Int_t is=AliPID::kSPECIES; is--;){ if(TMath::Abs(pdg)!=AliPID::ParticleCode(is)) continue; - ((TH3S*)((TObjArray*)arr->At(8))->At(ily))->Fill(pt0, 1.-pt/pt0, is); - //TODO ((TH3S*)((TObjArray*)arr->At(9))->At(ily))->Fill(1./pt0, (pt0/pt-1.)/TMath::Sqrt(covR[4]), is); - Double_t p0 = TMath::Sqrt(1.+ dzdx0*dzdx0)*pt0, - p = TMath::Sqrt(1.+ dzdx*dzdx)*pt; - ((TH3S*)((TObjArray*)arr->At(10))->At(ily))->Fill(p0, 1-p/p0, is); + ((TH3S*)((TObjArray*)arr->At(8))->At(ily))->Fill(pt0, pt/pt0-1., is); + ((TH3S*)((TObjArray*)arr->At(9))->At(ily))->Fill(1./pt0, (1./pt-1./pt0)/TMath::Sqrt(covR[6]), is); + Double_t p0 = TMath::Sqrt(1.+ dzdl0*dzdl0)*pt0, + p = TMath::Sqrt(1.+ dzdl*dzdl)*pt; + ((TH3S*)((TObjArray*)arr->At(10))->At(ily))->Fill(p0, p/p0-1., is); break; } @@ -673,30 +678,31 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track) AliTRDseedV1 tt(*fTracklet); tt.SetZref(0, z0 - (x0-xAnode)*dzdx0); tt.SetZref(1, dzdx0); + tt.SetReconstructor(&rec); tt.Fit(kTRUE, kTRUE); x= tt.GetX();y= tt.GetY();z= tt.GetZ(); dydx = tt.GetYfit(1); dx = x0 - x; - yt = y0 - dx*dydx0; - zt = z0 - dx*dzdx0; + ymc = y0 - dx*dydx0; + zmc = z0 - dx*dzdx0; Bool_t rc = tt.IsRowCross(); // add tracklet residuals for y and dydx arr = (TObjArray*)fContainer->At(kMCtracklet); if(!rc){ - dy = yt-y; + dy = y-ymc; Float_t dphi = (dydx - dydx0); - dphi /= 1.- dydx*dydx0; + dphi /= (1.- dydx*dydx0); ((TH2I*)arr->At(0))->Fill(dydx0, dy); if(tt.GetS2Y()>0.) ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(tt.GetS2Y())); ((TH2I*)arr->At(4))->Fill(dydx0, TMath::ATan(dphi)); } else { // add tracklet residuals for z - dz = zt-z; - ((TH2I*)arr->At(2))->Fill(dzdx0, dz); - if(tt.GetS2Z()>0.) ((TH2I*)arr->At(3))->Fill(dzdx0, dz/TMath::Sqrt(tt.GetS2Z())); + dz = z-zmc; + ((TH2I*)arr->At(2))->Fill(dzdl0, dz); + if(tt.GetS2Z()>0.) ((TH2I*)arr->At(3))->Fill(dzdl0, dz/TMath::Sqrt(tt.GetS2Z())); } // Fill Debug stream for tracklet @@ -725,37 +731,26 @@ TH1* AliTRDresolution::PlotMC(const AliTRDtrackV1 *track) fTracklet->ResetClusterIter(kFALSE); while((c = fTracklet->PrevCluster())){ Float_t q = TMath::Abs(c->GetQ()); - x = c->GetX(); y = c->GetY(); -// Int_t col = c->GetPadCol(); -// Int_t row = c->GetPadRow(); -// Double_t cw = pp->GetColSize(col); -// Double_t y0 = pp->GetColPos(col) + 0.5 * cw; -// Double_t s2 = AliTRDcalibDB::Instance()->GetPRFWidth(det, col, row); s2 *= s2; s2 -= - 1.5e-1; -// y = c->GetYloc(y0, s2, cw); y-=(xAnode-x)*exb; - - z = c->GetZ(); + x = c->GetX(); y = c->GetY();z = c->GetZ(); dx = x0 - x; - yt = y0 - dx*dydx0; - zt = z0 - dx*dzdx0; - dy = yt - (y - tilt*(z-zt)); - Float_t sx2 = dydx0*c->GetSX(c->GetLocalTimeBin()); sx2*=sx2; - Float_t sy2 = c->GetSY(c->GetLocalTimeBin()); sy2*=sy2; - + ymc= y0 - dx*dydx0; + zmc= z0 - dx*dzdx0; + dy = ymc - (y - tilt*(z-zmc)); // Fill Histograms if(q>20. && q<250.){ ((TH2I*)arr->At(0))->Fill(dydx0, dy); - ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(sx2+sy2)); + ((TH2I*)arr->At(1))->Fill(dydx0, dy/TMath::Sqrt(c->GetSigmaY2())); } // Fill calibration container - Float_t d = zr0 - zt; + Float_t d = zr0 - zmc; d -= ((Int_t)(2 * d)) / 2.0; if (d > 0.25) d = 0.5 - d; AliTRDclusterInfo *clInfo = new AliTRDclusterInfo; fMCcl->Add(clInfo); clInfo->SetCluster(c); clInfo->SetMC(pdg, label); - clInfo->SetGlobalPosition(yt, zt, dydx0, dzdx0); + clInfo->SetGlobalPosition(ymc, zmc, dydx0, dzdx0); clInfo->SetResolution(dy); clInfo->SetAnisochronity(d); clInfo->SetDriftLength(((c->GetPadTime()+0.5)*.1)*1.5); @@ -891,7 +886,7 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig) return kTRUE; case 13: //kMCtrackTRD [z] gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); - xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =2000.; + xy[0]=-1.; xy[1]=-700.; xy[2]=1.; xy[3] =1500.; ((TVirtualPad*)l->At(0))->cd(); if(!GetGraphPlot(&xy[0], kMCtrackTRD, 2)) break; xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5; @@ -899,22 +894,22 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig) if(!GetGraphPlot(&xy[0], kMCtrackTRD, 3)) break; return kTRUE; case 14: //kMCtrackTRD [phi/snp] - //gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); - xy[0]=-.2; xy[1]=-2.; xy[2]=.2; xy[3] =5.; - //((TVirtualPad*)l->At(0))->cd(); + gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); + xy[0]=-.2; xy[1]=-0.2; xy[2]=.2; xy[3] =2.; + ((TVirtualPad*)l->At(0))->cd(); if(!GetGraphPlot(&xy[0], kMCtrackTRD, 4)) break; -// xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5; -// ((TVirtualPad*)l->At(1))->cd(); -// if(!GetGraphPlot(&xy[0], kMCtrack, 3)) break; + xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5; + ((TVirtualPad*)l->At(1))->cd(); + if(!GetGraphPlot(&xy[0], kMCtrackTRD, 5)) break; return kTRUE; case 15: //kMCtrackTRD [theta/tgl] - //gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); - xy[0]=-1.; xy[1]=-50.; xy[2]=1.; xy[3] =70.; - //((TVirtualPad*)l->At(0))->cd(); + gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); + xy[0]=-1.; xy[1]=-0.5; xy[2]=1.; xy[3] =5.; + ((TVirtualPad*)l->At(0))->cd(); if(!GetGraphPlot(&xy[0], kMCtrackTRD, 6)) break; -// xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5; -// ((TVirtualPad*)l->At(1))->cd(); -// if(!GetGraphPlot(&xy[0], kMCtrack, 3)) break; + xy[0] = -.2; xy[1] = -0.5; xy[2] = .2; xy[3] = 2.5; + ((TVirtualPad*)l->At(1))->cd(); + if(!GetGraphPlot(&xy[0], kMCtrackTRD, 7)) break; return kTRUE; case 16: //kMCtrackTRD [pt] xy[0] = 0.; xy[1] = -5.; xy[2] = 12.; xy[3] = 7.; @@ -925,7 +920,16 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig) if(!GetGraphTrack(&xy[0], 8, il)) break; } return kTRUE; - case 17: //kMCtrackTRD [p] + case 17: //kMCtrackTRD [1/pt] pulls + xy[0] = 0.; xy[1] = -1.5; xy[2] = 2.; xy[3] = 2.; + gPad->Divide(2, 3, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); + for(Int_t il=0; ilAt(il); pad->cd(); + pad->SetMargin(0.07, 0.07, 0.1, 0.); + if(!GetGraphTrack(&xy[0], 9, il)) break; + } + return kTRUE; + case 18: //kMCtrackTRD [p] xy[0] = 0.; xy[1] = -7.5; xy[2] = 12.; xy[3] = 10.5; gPad->Divide(2, 3, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); for(Int_t il=0; ilDivide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); xy[0]=-.25; xy[1]=-50.; xy[2]=.25; xy[3] =800.; ((TVirtualPad*)l->At(0))->cd(); @@ -943,7 +947,7 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig) ((TVirtualPad*)l->At(1))->cd(); if(!GetGraphPlot(&xy[0], kMCtrackTPC, 1)) break; return kTRUE; - case 19: // kMCtrackTPC [z] + case 20: // kMCtrackTPC [z] gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); xy[0]=-1.; xy[1]=-500.; xy[2]=1.; xy[3] =800.; ((TVirtualPad*)l->At(0))->cd(); @@ -952,25 +956,25 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig) ((TVirtualPad*)l->At(1))->cd(); if(!GetGraphPlot(&xy[0], kMCtrackTPC, 3)) break; return kTRUE; - case 20: // kMCtrackTPC [phi|snp] + case 21: // kMCtrackTPC [phi|snp] gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); xy[0]=-.25; xy[1]=-0.5; xy[2]=.25; xy[3] =2.5; ((TVirtualPad*)l->At(0))->cd(); if(!GetGraphPlot(&xy[0], kMCtrackTPC, 4)) break; - //xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 2.5; + xy[0] = -.25; xy[1] = -0.5; xy[2] = .25; xy[3] = 1.5; ((TVirtualPad*)l->At(1))->cd(); if(!GetGraphPlot(&xy[0], kMCtrackTPC, 5)) break; return kTRUE; - case 21: // kMCtrackTPC [theta|tgl] + case 22: // kMCtrackTPC [theta|tgl] gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); - xy[0]=-1.; xy[1]=-10.5; xy[2]=1.; xy[3] =20.5; + xy[0]=-1.; xy[1]=-1.; xy[2]=1.; xy[3] =4.; ((TVirtualPad*)l->At(0))->cd(); if(!GetGraphPlot(&xy[0], kMCtrackTPC, 6)) break; - xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 2.5; + xy[0] = -1.; xy[1] = -0.5; xy[2] = 1.; xy[3] = 1.5; ((TVirtualPad*)l->At(1))->cd(); if(!GetGraphPlot(&xy[0], kMCtrackTPC, 7)) break; return kTRUE; - case 22: // kMCtrackTPC [pt] + case 23: // kMCtrackTPC [pt] gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); xy[0] = 0.; xy[1] = -.8; xy[2] = 12.; xy[3] = 2.3; ((TVirtualPad*)l->At(0))->cd(); @@ -979,7 +983,7 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig) ((TVirtualPad*)l->At(1))->cd(); if(!GetGraphTrackTPC(xy, 9)) break; return kTRUE; - case 23: // kMCtrackTPC [p] + case 24: // kMCtrackTPC [p] gPad->Divide(2, 1, 1.e-5, 1.e-5); l=gPad->GetListOfPrimitives(); xy[0] = 0.; xy[1] = -.8; xy[2] = 12.; xy[3] = 2.3; pad = ((TVirtualPad*)l->At(0));pad->cd(); @@ -990,10 +994,10 @@ Bool_t AliTRDresolution::GetRefFigure(Int_t ifig) pad->SetMargin(0.12, 0.12, 0.1, 0.04); if(!GetGraphTrackTPC(xy, 11)) break; return kTRUE; - case 24: // kMCtrackTOF [z] + case 25: // kMCtrackTOF [z] return kTRUE; } - AliInfo(Form("Reference plot [%d] missing result", ifig)); + AliWarning(Form("Reference plot [%d] missing result", ifig)); return kFALSE; } @@ -1039,7 +1043,7 @@ Bool_t AliTRDresolution::PostProcess() gm->SetNameTitle(Form("m_%d%02d%d", ig, ic, is), ""); } continue; - } else if(ig==kMCtrackTRD&&(ic==8||ic==10)){ // TRD momentum plot + } else if(ig==kMCtrackTRD&&(ic==8||ic==9||ic==10)){ // TRD momentum plot TObjArray *aaS, *aaM; aS->AddAt(aaS = new TObjArray(AliTRDgeometry::kNlayer), ic); aM->AddAt(aaM = new TObjArray(AliTRDgeometry::kNlayer), ic); @@ -1139,13 +1143,13 @@ Bool_t AliTRDresolution::PostProcess() Process2D(kMCtrackTRD, 2, &fg, 1.e4); // z Process2D(kMCtrackTRD, 3, &fg); // z PULL Process2D(kMCtrackTRD, 4, &fg, 1.e3); // phi - //Process(kMCtrack, 5, &fg); // snp PULL + Process2D(kMCtrackTRD, 5, &fg); // snp PULL Process2D(kMCtrackTRD, 6, &fg, 1.e3); // theta - //Process(kMCtrack, 7, &fg); // tgl PULL + Process2D(kMCtrackTRD, 7, &fg); // tgl PULL Process4D(kMCtrackTRD, 8, &fg, 1.e2); // pt resolution - //Process4D(kMCtrack, 9, &fg); // 1/pt pulls + Process4D(kMCtrackTRD, 9, &fg); // 1/pt pulls Process4D(kMCtrackTRD, 10, &fg, 1.e2); // p resolution - fNRefFigures = 17; + fNRefFigures = 18; // TRACK TPC RESOLUTION/PULLS Process2D(kMCtrackTPC, 0, &fg, 1.e4);// y resolution @@ -1160,12 +1164,12 @@ Bool_t AliTRDresolution::PostProcess() Process3D(kMCtrackTPC, 9, &fg); // 1/pt pulls Process3D(kMCtrackTPC, 10, &fg, 1.e2);// p resolution Process3D(kMCtrackTPC, 11, &fg); // p pulls - fNRefFigures = 23; + fNRefFigures = 24; // TRACK HMPID RESOLUTION/PULLS Process2D(kMCtrackTOF, 0, &fg, 1.e4); // z towards TOF Process2D(kMCtrackTOF, 1, &fg); // z towards TOF - fNRefFigures = 24; + fNRefFigures = 25; return kTRUE; } @@ -1439,7 +1443,7 @@ TObjArray* AliTRDresolution::Histos() arr->AddAt(h, 3); // Kalman track SNP if(!(h = (TH2I*)gROOT->FindObject("hMcTrkSNP"))){ - h = new TH2I("hMcTrkSNP", "Track Phi Resolution", 60, -.3, .3, 100, -.02, .02); + h = new TH2I("hMcTrkSNP", "Track Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3); h->GetXaxis()->SetTitle("tg(#phi)"); h->GetYaxis()->SetTitle("#Delta #phi [rad]"); h->GetZaxis()->SetTitle("entries"); @@ -1455,7 +1459,7 @@ TObjArray* AliTRDresolution::Histos() arr->AddAt(h, 5); // Kalman track TGL if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTGL"))){ - h = new TH2I("hMcTrkTGL", "Track Theta Resolution", 100, -1., 1., 100, -.1, .1); + h = new TH2I("hMcTrkTGL", "Track Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3); h->GetXaxis()->SetTitle("tg(#theta)"); h->GetYaxis()->SetTitle("#Delta#theta [rad]"); h->GetZaxis()->SetTitle("entries"); @@ -1476,7 +1480,7 @@ TObjArray* AliTRDresolution::Histos() arr2->SetName("Track Pt Resolution"); for(Int_t il=0; ilFindObject(Form("hMcTrkPt%d", il)))){ - h3 = new TH3S(Form("hMcTrkPt%d", il), "Track Pt Resolution", 40, 0., 20., 150, -.15, .15, n, -.5, n-.5); + h3 = new TH3S(Form("hMcTrkPt%d", il), "Track Pt Resolution", 40, 0., 20., 150, -.1, .2, n, -.5, n-.5); h3->GetXaxis()->SetTitle("p_{t} [GeV/c]"); h3->GetYaxis()->SetTitle("#Delta p_{t}/p_{t}^{MC}"); h3->GetZaxis()->SetTitle("SPECIES"); @@ -1493,13 +1497,14 @@ TObjArray* AliTRDresolution::Histos() h3->GetYaxis()->SetTitle("#Delta(1/p_{t})/#sigma(1/p_{t}) "); h3->GetZaxis()->SetTitle("SPECIES"); } else h3->Reset(); + arr2->AddAt(h3, il); } // Kalman track P resolution arr->AddAt(arr2 = new TObjArray(AliTRDgeometry::kNlayer), 10); arr2->SetName("Track P Resolution [PID]"); for(Int_t il=0; ilFindObject(Form("hMcTrkP%d", il)))){ - h3 = new TH3S(Form("hMcTrkP%d", il), "Track P Resolution", 40, 0., 20., 150, -.25, .25, n, -.5, n-.5); + h3 = new TH3S(Form("hMcTrkP%d", il), "Track P Resolution", 40, 0., 20., 150, -.15, .35, n, -.5, n-.5); h3->GetXaxis()->SetTitle("p [GeV/c]"); h3->GetYaxis()->SetTitle("#Delta p/p^{MC}"); h3->GetZaxis()->SetTitle("SPECIES"); @@ -1544,7 +1549,7 @@ TObjArray* AliTRDresolution::Histos() arr->AddAt(h, 3); // Kalman track SNP if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCSNP"))){ - h = new TH2I("hMcTrkTPCSNP", "Track[TPC] Phi Resolution", 60, -.3, .3, 100, -.02, .02); + h = new TH2I("hMcTrkTPCSNP", "Track[TPC] Phi Resolution", 60, -.3, .3, 100, -5e-3, 5e-3); h->GetXaxis()->SetTitle("tg(#phi)"); h->GetYaxis()->SetTitle("#Delta #phi [rad]"); h->GetZaxis()->SetTitle("entries"); @@ -1560,7 +1565,7 @@ TObjArray* AliTRDresolution::Histos() arr->AddAt(h, 5); // Kalman track TGL if(!(h = (TH2I*)gROOT->FindObject("hMcTrkTPCTGL"))){ - h = new TH2I("hMcTrkTPCTGL", "Track[TPC] Theta Resolution", 100, -1., 1., 100, -.1, .1); + h = new TH2I("hMcTrkTPCTGL", "Track[TPC] Theta Resolution", 100, -1., 1., 100, -5e-3, 5e-3); h->GetXaxis()->SetTitle("tg(#theta)"); h->GetYaxis()->SetTitle("#Delta#theta [rad]"); h->GetZaxis()->SetTitle("entries"); @@ -1576,7 +1581,7 @@ TObjArray* AliTRDresolution::Histos() arr->AddAt(h, 7); // Kalman track Pt resolution if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCPt"))){ - h3 = new TH3S("hMcTrkTPCPt", "Track[TPC] Pt Resolution", 80, 0., 20., 150, -.15, .15, n, -.5, n-.5); + h3 = new TH3S("hMcTrkTPCPt", "Track[TPC] Pt Resolution", 80, 0., 20., 150, -.1, .2, n, -.5, n-.5); h3->GetXaxis()->SetTitle("p_{t} [GeV/c]"); h3->GetYaxis()->SetTitle("#Delta p_{t}/p_{t}^{MC}"); h3->GetZaxis()->SetTitle("SPECIES"); @@ -1592,7 +1597,7 @@ TObjArray* AliTRDresolution::Histos() arr->AddAt(h3, 9); // Kalman track P resolution if(!(h3 = (TH3S*)gROOT->FindObject("hMcTrkTPCP"))){ - h3 = new TH3S("hMcTrkTPCP", "Track[TPC] P Resolution", 80, 0., 20., 150, -.25, .25, n, -.5, n-.5); + h3 = new TH3S("hMcTrkTPCP", "Track[TPC] P Resolution", 80, 0., 20., 150, -.15, .35, n, -.5, n-.5); h3->GetXaxis()->SetTitle("p [GeV/c]"); h3->GetYaxis()->SetTitle("#Delta p/p^{MC}"); h3->GetZaxis()->SetTitle("SPECIES"); diff --git a/TRD/qaRec/info/AliTRDclusterInfo.cxx b/TRD/qaRec/info/AliTRDclusterInfo.cxx index 93af9963ae5..dbd8d2f2931 100644 --- a/TRD/qaRec/info/AliTRDclusterInfo.cxx +++ b/TRD/qaRec/info/AliTRDclusterInfo.cxx @@ -35,7 +35,7 @@ AliTRDclusterInfo::AliTRDclusterInfo() } //_________________________________________________ -void AliTRDclusterInfo::SetCluster(const AliTRDcluster *c, Float_t *cov) +void AliTRDclusterInfo::SetCluster(const AliTRDcluster *c) { if(!c) return; fDet = c->GetDetector(); @@ -45,8 +45,9 @@ void AliTRDclusterInfo::SetCluster(const AliTRDcluster *c, Float_t *cov) fQ = TMath::Abs(c->GetQ()); fLocalTime = c->GetLocalTimeBin(); fYd = c->GetCenter(); - - if(cov) memcpy(cov, fCovCl, 3*sizeof(Float_t)); + fCovCl[0] = c->GetSigmaY2(); + fCovCl[1] = 0.; + fCovCl[2] = c->GetSigmaZ2(); } //_________________________________________________ diff --git a/TRD/qaRec/info/AliTRDclusterInfo.h b/TRD/qaRec/info/AliTRDclusterInfo.h index 8b1db24ed1f..b946fa38e88 100644 --- a/TRD/qaRec/info/AliTRDclusterInfo.h +++ b/TRD/qaRec/info/AliTRDclusterInfo.h @@ -34,7 +34,7 @@ public: void Print(Option_t *opt="") const; void SetAnisochronity(Float_t d) {fD = d;} - void SetCluster(const AliTRDcluster *c, Float_t *cov=0x0); + void SetCluster(const AliTRDcluster *c); void SetMC(Int_t pdg, Int_t label){ fPdg = pdg; fLbl = label;} diff --git a/TRD/qaRec/macros/makeResults.C b/TRD/qaRec/macros/makeResults.C index 59e357036ea..a75a336d677 100644 --- a/TRD/qaRec/macros/makeResults.C +++ b/TRD/qaRec/macros/makeResults.C @@ -68,7 +68,6 @@ Char_t *libs[] = {"libProofPlayer.so", "libANALYSIS.so", "libTRDqaRec.so"}; // define setup -gStyle->SetOptStat(0); Bool_t mc = kTRUE; Bool_t friends = kTRUE; TCanvas *c = 0x0; @@ -93,6 +92,8 @@ void makeResults(Char_t *opt = "ALL", const Char_t *files=0x0, Bool_t kGRID=kFAL return; } + gStyle->SetOptStat(0); + gStyle->SetOptFit(0); if(files) mergeProd("TRD.Performance.root", files); Int_t fSteerTask = ParseOptions(opt); -- 2.43.0