From cb29fcbfbe0a46f458d4a397e46c714b52cb204b Mon Sep 17 00:00:00 2001 From: shahoian Date: Wed, 9 Jul 2014 15:05:27 +0200 Subject: [PATCH] Fix from Alex for FPE in TRD: ALIROOT-5503 and tracklet calibration: PWGPP-2 --- PWGPP/TRD/AliTRDcheckTRK.cxx | 8 ++--- TRD/AliTRDrecoParam.cxx | 46 +++++++++++++++++--------- TRD/AliTRDrecoParam.h | 14 ++++---- TRD/AliTRDseedV1.cxx | 63 ++++++++++++++++-------------------- TRD/AliTRDseedV1.h | 15 +++++---- TRD/AliTRDtrackerV1.cxx | 2 +- 6 files changed, 76 insertions(+), 72 deletions(-) diff --git a/PWGPP/TRD/AliTRDcheckTRK.cxx b/PWGPP/TRD/AliTRDcheckTRK.cxx index 507461743f4..11a731dd748 100644 --- a/PWGPP/TRD/AliTRDcheckTRK.cxx +++ b/PWGPP/TRD/AliTRDcheckTRK.cxx @@ -240,11 +240,9 @@ Bool_t AliTRDcheckTRK::PropagateKalman(AliTRDtrackV1 &t, AliExternalTrackParam * } if(HasTrkltRefit()){ // if(!tr->FitRobust(tt.Charge()>0.)) printf("W - AliTRDcheckTRK::PropagateKalman :: FitRobust() failed for Det[%03d]\n", det); - if(!tr->FitRobust(AliTRDgeometry::GetPadPlane(det), prod>0., tt.Charge())) printf("W - AliTRDcheckTRK::PropagateKalman :: FitRobust() failed for Det[%03d]\n", det); - else { - TGeoHMatrix *matrix = AliTRDgeometry::GetClusterMatrix(det); - if (matrix) tr->SetXYZ(matrix); - } + TGeoHMatrix *matrix = AliTRDgeometry::GetClusterMatrix(det); + if(!tr->FitRobust(AliTRDgeometry::GetPadPlane(det), matrix, tt.GetBz(), tt.Charge())) printf("W - AliTRDcheckTRK::PropagateKalman :: FitRobust() failed for Det[%03d]\n", det); + else tr->SetXYZ(matrix); } if(!AliTRDtrackerV1::PropagateToX(tt, tr->GetX0(), fgKalmanStep)) continue; if(!tt.GetTrackIn()) tt.SetTrackIn(); diff --git a/TRD/AliTRDrecoParam.cxx b/TRD/AliTRDrecoParam.cxx index f790fea9c36..294d074d9c4 100644 --- a/TRD/AliTRDrecoParam.cxx +++ b/TRD/AliTRDrecoParam.cxx @@ -156,8 +156,8 @@ AliTRDrecoParam::AliTRDrecoParam(const AliTRDrecoParam &ref) // tracklet params memcpy(fdzdxCorrFactor, ref.fdzdxCorrFactor, 2*sizeof(Double_t)); memcpy(fdzdxCorrRCbias, ref.fdzdxCorrRCbias, 2*sizeof(Double_t)); - memcpy(fYcorrTailCancel, ref.fdzdxCorrRCbias, 6*sizeof(Double_t)); - memcpy(fS2Ycorr, ref.fS2Ycorr, 2*sizeof(Double_t)); + memcpy(fYcorrTailCancel, ref.fdzdxCorrRCbias, 12*sizeof(Double_t)); + memcpy(fS2Ycorr, ref.fS2Ycorr, 4*sizeof(Double_t)); } //______________________________________________________________ @@ -215,8 +215,8 @@ AliTRDrecoParam& AliTRDrecoParam::operator=(const AliTRDrecoParam &ref) // tracklet params memcpy(fdzdxCorrFactor, ref.fdzdxCorrFactor, 2*sizeof(Double_t)); memcpy(fdzdxCorrRCbias, ref.fdzdxCorrRCbias, 2*sizeof(Double_t)); - memcpy(fYcorrTailCancel, ref.fdzdxCorrRCbias, 6*sizeof(Double_t)); - memcpy(fS2Ycorr, ref.fS2Ycorr, 2*sizeof(Double_t)); + memcpy(fYcorrTailCancel, ref.fdzdxCorrRCbias, 12*sizeof(Double_t)); + memcpy(fS2Ycorr, ref.fS2Ycorr, 4*sizeof(Double_t)); return *this; } @@ -351,13 +351,16 @@ void AliTRDrecoParam::SetTrackletParams(Double_t *par) fdzdxCorrRCbias[1] = par[3]; // dz/dx < 0 /// correct x_cross for the bias in dzdx fdzdxXcrossFactor = par[4]; - // y linear q/pt correction due to wrong tail cancellation. - fYcorrTailCancel[0][0] = par[5]; fYcorrTailCancel[0][1] = par[6]; // opposite sign !RC - fYcorrTailCancel[1][0] = par[7]; fYcorrTailCancel[1][1] = par[8]; // same sign !RC - fYcorrTailCancel[2][0] = par[9]; fYcorrTailCancel[2][1] = par[10]; // RC + // y correction due to wrong tail cancellation. + fYcorrTailCancel[0][0] = par[5];fYcorrTailCancel[0][1] = par[6];fYcorrTailCancel[0][2] = par[7]; + fYcorrTailCancel[1][0] = par[8];fYcorrTailCancel[1][1] = par[9];fYcorrTailCancel[1][2] = par[10]; + fYcorrTailCancel[2][0] = par[11];fYcorrTailCancel[2][1] = par[12];fYcorrTailCancel[2][3] = par[13]; + fYcorrTailCancel[3][0] = par[14];fYcorrTailCancel[3][1] = par[15];fYcorrTailCancel[3][3] = par[16]; // inflation factor of error parameterization in r-phi due to wrong estimation of residuals. - fS2Ycorr[0] = par[11]; // opposite sign, - fS2Ycorr[1] = par[12]; // same sign + fS2Ycorr[0] = par[17]; + fS2Ycorr[1] = par[18]; + fS2Ycorr[2] = par[19]; + fS2Ycorr[3] = par[20]; } else { // correct dzdx for the bias in z @@ -368,12 +371,23 @@ void AliTRDrecoParam::SetTrackletParams(Double_t *par) fdzdxCorrRCbias[1] = -0.012; // dz/dx < 0 /// correct x_cross for the bias in dzdx fdzdxXcrossFactor = 0.14; - // y linear q/pt correction due to wrong tail cancellation. - fYcorrTailCancel[0][0] = 0.; fYcorrTailCancel[0][1] = 0.027; // opposite sign !RC - fYcorrTailCancel[1][0] = 0.04; fYcorrTailCancel[1][1] = 0.027; // same sign !RC - fYcorrTailCancel[2][0] = 0.013;fYcorrTailCancel[2][1] = 0.018; // RC + // y correction due to wrong tail cancellation. + // bz<0 && !RC + fYcorrTailCancel[0][0] = 0.04; fYcorrTailCancel[0][1] = 2.151; fYcorrTailCancel[0][2] = 0.013; + // bz>0 && !RC + fYcorrTailCancel[1][0] = 0.034; fYcorrTailCancel[1][1] = 1.817; fYcorrTailCancel[1][2] = -0.01; + // bz<0 && RC + fYcorrTailCancel[2][0] = 0.04; fYcorrTailCancel[2][1] = 2.513; fYcorrTailCancel[2][2] = 0.015; + // bz>0 && RC + fYcorrTailCancel[3][0] = 0.034; fYcorrTailCancel[3][1] = 2.476; fYcorrTailCancel[3][2] = -0.01; // inflation factor of error parameterization in r-phi due to wrong estimation of residuals. - fS2Ycorr[0] = 5.; // opposite sign, - fS2Ycorr[1] = 10; // same sign + // chg<0 && !RC + fS2Ycorr[0] = 5.52; + // chg>0 && !RC + fS2Ycorr[1] = 3.61; + // chg<0 && RC + fS2Ycorr[2] = 4.84; + // chg>0 && RC + fS2Ycorr[3] = 3.24; } } diff --git a/TRD/AliTRDrecoParam.h b/TRD/AliTRDrecoParam.h index d4c9ceefebe..9edca2b9617 100644 --- a/TRD/AliTRDrecoParam.h +++ b/TRD/AliTRDrecoParam.h @@ -94,8 +94,8 @@ public: Double_t GetCorrDZDXbiasRC(Bool_t dzdx) const { return fdzdxCorrRCbias[dzdx];} Double_t GetCorrDZDX(Bool_t rc) const { return fdzdxCorrFactor[rc];} Double_t GetCorrDZDXxcross() const { return fdzdxXcrossFactor;} - void GetYcorrTailCancel(Int_t it, Double_t par[2]) const; - Double_t GetS2Ycorr(Bool_t sgn) const { return fS2Ycorr[sgn];} + void GetYcorrTailCancel(Int_t it, Double_t par[3]) const; + Double_t GetS2Ycorr(Bool_t rc, Bool_t chg) const { return fS2Ycorr[2*rc+chg];} Bool_t IsArgon() const { return TESTBIT(fFlags, kDriftGas); } Bool_t IsCheckTimeConsistency() const { return kCheckTimeConsistency;} @@ -204,8 +204,8 @@ private: Double_t fdzdxCorrFactor[2]; // correction of dzdx estimation due to z bias; [0] for !RC, [1] for RC Double_t fdzdxCorrRCbias[2]; // correction of dzdx estimation bias for RC; [0] for dz/dx>0, [1] for dz/dx<0 Double_t fdzdxXcrossFactor; // bias in dzdx of estimated xcross [RC] - Double_t fYcorrTailCancel[3][2]; // y linear q/pt correction due to wrong tail cancellation. [0] opposite sign !RC, [1] same sign !RC, [2] RC - Double_t fS2Ycorr[2]; // inflation factor of error parameterization in r-phi due to wrong estimation of residuals. [0] opposite sign, [1] same sign + Double_t fYcorrTailCancel[4][3]; // y correction due to wrong tail cancellation. [0] bz<0 && !RC, [1] bz>0 && !RC, [2] bz<0 && RC [3] bz>0 && RC + Double_t fS2Ycorr[4]; // inflation factor of error parameterization in r-phi due to wrong estimation of residuals. // Clusterization parameter Double_t fMinMaxCutSigma; // Threshold sigma noise pad middle @@ -274,10 +274,10 @@ inline void AliTRDrecoParam::SetTCParams(Double_t *par) //___________________________________________________ -inline void AliTRDrecoParam::GetYcorrTailCancel(Int_t it, Double_t par[2]) const +inline void AliTRDrecoParam::GetYcorrTailCancel(Int_t it, Double_t par[3]) const { - if(it<0||it>2) return; - par[0] = fYcorrTailCancel[it][0]; par[1] = fYcorrTailCancel[it][1]; + if(it<0||it>3) return; + memcpy(par, fYcorrTailCancel[it], 3*sizeof(Double_t)); } //___________________________________________________ diff --git a/TRD/AliTRDseedV1.cxx b/TRD/AliTRDseedV1.cxx index b718f5e860a..9b0e50603b4 100644 --- a/TRD/AliTRDseedV1.cxx +++ b/TRD/AliTRDseedV1.cxx @@ -539,7 +539,7 @@ Int_t AliTRDseedV1::GetChargeGaps(Float_t sz[kNtb], Float_t pos[kNtb], Int_t isz //____________________________________________________________________ -Double_t AliTRDseedV1::EstimatedCrossPoint(AliTRDpadPlane *pp) +Double_t AliTRDseedV1::EstimatedCrossPoint(AliTRDpadPlane *pp, Float_t bz) { // Algorithm to estimate cross point in the x-z plane for pad row cross tracklets or the z coordinate of pad row without pad row cross in the local chamber coordinates. // Returns variance of the radial offset from anode wire in case of raw cross or 0 otherwise. @@ -633,7 +633,7 @@ Double_t AliTRDseedV1::EstimatedCrossPoint(AliTRDpadPlane *pp) fZfit[1] = (fZfit[0]+zoff)/dx; // correct dzdx for the bias - UnbiasDZDX(IsRowCross()); + UnbiasDZDX(IsRowCross(), bz); if(IsRowCross()){ // correct x_cross/sigma(x_cross) for the bias in dzdx const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); @@ -652,32 +652,28 @@ Double_t AliTRDseedV1::EstimatedCrossPoint(AliTRDpadPlane *pp) } //____________________________________________________________________ -void AliTRDseedV1::UnbiasDZDX(Bool_t rc) +void AliTRDseedV1::UnbiasDZDX(Bool_t rc, Float_t bz) { // correct dzdx for the bias in z according to MC const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); if(!recoParam) return; - fZfit[1] *= recoParam->GetCorrDZDX(rc); + fZfit[1] *= recoParam->GetCorrDZDX(rc)-(bz>0?0.01:0.); if(rc) fZfit[1] += recoParam->GetCorrDZDXbiasRC(fZfit[1]<0); } //____________________________________________________________________ -Double_t AliTRDseedV1::UnbiasY(Bool_t rc, Bool_t sgn, Int_t chg) +Double_t AliTRDseedV1::UnbiasY(Bool_t rc, Float_t bz) { // correct y coordinate for tail cancellation. This should be fixed by considering TC as a function of q/pt. // rc : TRUE if tracklet crosses rows -// sgn : TRUE if track has same sign with magnetic field -// chg : -1 for negative particles, +1 for the rest +// bz : magnetic field z component const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam(); if(!recoParam) return 0.; - Double_t par[2]={0.}; - if(rc) recoParam->GetYcorrTailCancel(2, par); - else{ - if(sgn && 1./fPt > 1.5) recoParam->GetYcorrTailCancel(1, par); - else if(!sgn) recoParam->GetYcorrTailCancel(0, par); - } - return par[0]+par[1]*chg/fPt; + Double_t par[3]={0.}; + Int_t idx(2*(rc?1:0)+Int_t(bz>0)); + recoParam->GetYcorrTailCancel(idx, par); + return par[0]*TMath::Sin(par[1]*fYref[1])+par[2]; } @@ -2161,7 +2157,7 @@ Bool_t AliTRDseedV1::Fit(UChar_t opt) //____________________________________________________________________ -Bool_t AliTRDseedV1::FitRobust(AliTRDpadPlane *pp, Bool_t sgn, Int_t chg, Int_t opt) +Bool_t AliTRDseedV1::FitRobust(AliTRDpadPlane *pp, TGeoHMatrix *mDet, Float_t bz, Int_t chg, Int_t opt) { // // Linear fit of the clusters attached to the tracklet @@ -2199,7 +2195,7 @@ Bool_t AliTRDseedV1::FitRobust(AliTRDpadPlane *pp, Bool_t sgn, Int_t chg, Int_t Double_t freq(cp?cp->GetSamplingFrequency():10.); // evaluate locally z and dzdx from TRD only information - if(EstimatedCrossPoint(pp)<0.) return kFALSE; + if(EstimatedCrossPoint(pp, bz)<0.) return kFALSE; //printf("D%03d RC[%c] dzdx[%f %f] opt[%d]\n", fDet, IsRowCross()?'y':'n', fZref[1], fZfit[1], opt); Double_t //xchmb = 0.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick(), @@ -2210,13 +2206,11 @@ Bool_t AliTRDseedV1::FitRobust(AliTRDpadPlane *pp, Bool_t sgn, Int_t chg, Int_t Double_t xc[kNclusters], yc[kNclusters], dz(0.), dzdx(0.), s2dz(0.), s2dzdx(0.), sy[kNclusters], s2x((8.33e-2/freq/freq+1.56e-2)*fVD*fVD), // error of 1tb + error of mean time (TRF) - t2(fPad[2]*fPad[2]), - cs(0.); + t2(fPad[2]*fPad[2]), loc[3]={0.}; Int_t n(0), // clusters used in fit - row[]={-1, -1},// pad row spanned by the tracklet - col(-1); // pad column of current cluster - Double_t ycorr(UnbiasY(IsRowCross(), sgn, chg)), - kS2Ycorr(recoParam->GetS2Ycorr(sgn)); + row[]={-1, -1};// pad row spanned by the tracklet + Double_t ycorr(UnbiasY(IsRowCross(), bz)), + kS2Ycorr(recoParam->GetS2Ycorr(IsRowCross(), chg>0)); AliTRDcluster *c(NULL), **jc = &fClusters[0]; for(Int_t ic=0; icGetPadCol()){ - col = c->GetPadCol(); - cs = pp->GetColSize(col); - } - //Use local cluster coordinates - the code should be identical with AliTRDtransform::Transform() !!! - //A.Bercuci 27.11.13 - xc[n] = c->GetXloc(fT0, fVD); // c->GetX(); - yc[n] = c->GetYloc(pp->GetColPos(col) + .5*cs, fS2PRF, cs) - xc[n]*fExB; //c->GetY(); + //Use local cluster coordinates + //A.Bercuci 27.11.13/30.06.14 + Double_t trk[] = {c->GetX(), c->GetY(), c->GetZ()}; + mDet->MasterToLocal(trk, loc); + xc[n] = AliTRDgeometry::AnodePos()-loc[0]; //c->GetXloc(fT0, fVD); // c->GetX(); + yc[n] = loc[1]; //c->GetYloc(pp->GetColPos(col) + .5*cs, fS2PRF, cs) - xc[n]*fExB; //c->GetY(); yc[n]-= fPad[2]*(dz+xc[n]*dzdx); yc[n]-= ycorr; if(IsRowCross()){ // estimate closest distance to anode wire @@ -2293,14 +2285,13 @@ Bool_t AliTRDseedV1::FitRobust(AliTRDpadPlane *pp, Bool_t sgn, Int_t chg, Int_t break; } } - if(col != c->GetPadCol()){ - col = c->GetPadCol(); - cs = pp->GetColSize(col); - } + //Use local cluster coordinates - the code should be identical with AliTRDtransform::Transform() !!! //A.Bercuci 27.11.13 - xc[n] = c->GetXloc(fT0, fVD); // c->GetX(); - yc[n] = c->GetYloc(pp->GetColPos(col) + .5*cs, fS2PRF, cs) - xc[n]*fExB ; + Double_t trk[] = {c->GetX(), c->GetY(), c->GetZ()}; + mDet->MasterToLocal(trk, loc); + xc[n] = AliTRDgeometry::AnodePos()-loc[0]; //c->GetXloc(fT0, fVD); // c->GetX(); + yc[n] = loc[1]; //c->GetYloc(pp->GetColPos(col) + .5*cs, fS2PRF, cs) - xc[n]*fExB; //c->GetY(); yc[n]-= fPad[2]*(dz+xc[n]*dzdx); yc[n]-= ycorr; diff --git a/TRD/AliTRDseedV1.h b/TRD/AliTRDseedV1.h index 80bc7828383..a8dd3b44ab8 100644 --- a/TRD/AliTRDseedV1.h +++ b/TRD/AliTRDseedV1.h @@ -91,8 +91,8 @@ public: void CookLabels(); Bool_t CookPID(); Bool_t Fit(UChar_t opt=0); // OBSOLETE - Bool_t FitRobust(AliTRDpadPlane *pp, Bool_t sgn, Int_t chg, Int_t opt=0); - Double_t EstimatedCrossPoint(AliTRDpadPlane *pp); + Bool_t FitRobust(AliTRDpadPlane *pp, TGeoHMatrix *mdet, Float_t bz, Int_t chg, Int_t opt=0); + Double_t EstimatedCrossPoint(AliTRDpadPlane *pp, Float_t bz); Bool_t Init(const AliTRDtrackV1 *track); void Init(const AliRieman *fit); Bool_t IsEqual(const TObject *inTracklet) const; @@ -216,8 +216,8 @@ public: protected: void Copy(TObject &ref) const; - void UnbiasDZDX(Bool_t rc); - Double_t UnbiasY(Bool_t rc, Bool_t sgn, Int_t chg); + void UnbiasDZDX(Bool_t rc, Float_t bz); + Double_t UnbiasY(Bool_t rc, Float_t bz); private: inline void SetN(Int_t n); @@ -312,11 +312,12 @@ Double_t AliTRDseedV1::GetS2XcrossDZDX(Double_t absdzdx) const //____________________________________________________________ Double_t AliTRDseedV1::GetS2DZDX(Float_t dzdx) const { -// Double_t p0[] = {0.03925, 0.03178}, -// p1[] = {0.06316, 0.06669}; + // Error parametrization for dzdx. + // TODO Should be layer dependent + Double_t p0[] = {0.02835, 0.03925}, p1[] = {0.04746, 0.06316}; - + Double_t s2(p0[IsRowCross()]+p1[IsRowCross()]*dzdx*dzdx); s2*=s2; return s2; diff --git a/TRD/AliTRDtrackerV1.cxx b/TRD/AliTRDtrackerV1.cxx index b4818993733..cca2e272c60 100644 --- a/TRD/AliTRDtrackerV1.cxx +++ b/TRD/AliTRDtrackerV1.cxx @@ -956,7 +956,7 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t) // tilt correction options // 0 : no correction // 2 : pseudo tilt correction - if(!ptrTracklet->FitRobust(fGeom->GetPadPlane(ily, stk), prod>0., t.Charge())){ + if(!ptrTracklet->FitRobust(fGeom->GetPadPlane(ily, stk), matrix, t.GetBz(), t.Charge())){ t.SetErrStat(AliTRDtrackV1::kNoFit, ily); AliDebug(4, "Failed Tracklet Fit"); continue; -- 2.43.0