X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=ITS%2FAliITSTPArrayFit.cxx;h=fa132ef406c28a1e331de1b4db601c2a14d27658;hp=e707bdb4b10ebf18c58a8b8e9835fc18058426d6;hb=fcb289102b9b5323c24bf4ec175fb97fcc0260b8;hpb=ef24eb3bd644397bb588db97d293cd8dfa882f27 diff --git a/ITS/AliITSTPArrayFit.cxx b/ITS/AliITSTPArrayFit.cxx index e707bdb4b10..fa132ef406c 100644 --- a/ITS/AliITSTPArrayFit.cxx +++ b/ITS/AliITSTPArrayFit.cxx @@ -124,7 +124,7 @@ AliITSTPArrayFit::AliITSTPArrayFit(const AliITSTPArrayFit &src) : { // copy constructor InitAux(); - memcpy(fCovI,src.fCovI,fNPBooked*6*sizeof(Double_t)); + memcpy(fCovI,src.fCovI,fNPBooked*kNCov*sizeof(Double_t)); for (int i=kMaxParam;i--;) fParams[i] = src.fParams[i]; for (int i=kMaxParamSq;i--;) fParamsCov[i] = src.fParamsCov[i]; memcpy(fCurT,src.fCurT,fNPBooked*sizeof(Double_t)); @@ -148,7 +148,7 @@ AliITSTPArrayFit &AliITSTPArrayFit::operator =(const AliITSTPArrayFit& src) fPntFirst = src.fPntFirst; fPntLast = src.fPntLast; InitAux(); - memcpy(fCovI,src.fCovI,fNPBooked*6*sizeof(Double_t)); + memcpy(fCovI,src.fCovI,fNPBooked*kNCov*sizeof(Double_t)); for (int i=kMaxParam;i--;) fParams[i] = src.fParams[i]; for (int i=kMaxParamSq;i--;) fParamsCov[i] = src.fParamsCov[i]; SetParAxis(src.fParAxis); @@ -211,6 +211,7 @@ void AliITSTPArrayFit::AttachPoints(const AliTrackPointArray* points, Int_t pfir for (int i=kMaxParamSq;i--;) fParamsCov[i] = 0; // InvertPointsCovMat(); + ResetCovIScale(); // } @@ -309,12 +310,13 @@ void AliITSTPArrayFit::InitAux() if (fCovI) delete[] fCovI; if (fCurT) delete[] fCurT; // - fCovI = new Double_t[6*fNPBooked]; + fCovI = new Double_t[kNCov*fNPBooked]; fCurT = new Double_t[fNPBooked+kMaxLrITS]; fElsId = new Int_t[fNPBooked+kMaxLrITS]; fElsDR = new Double_t[fNPBooked+kMaxLrITS]; memset(fElsDR,0,(fNPBooked+kMaxLrITS)*sizeof(Double_t)); - memset(fCovI,0,fNPBooked*6*sizeof(Double_t)); + memset(fCovI,0,fNPBooked*kNCov*sizeof(Double_t)); + ResetCovIScale(); // } @@ -407,10 +409,10 @@ Int_t AliITSTPArrayFit::ChoseParAxis() const } //____________________________________________________ -Double_t AliITSTPArrayFit::GetPosition(Double_t *xyzPCA, const Double_t *xyz, const Double_t *covI) const +Double_t AliITSTPArrayFit::GetPosition(Double_t *xyzPCA, const Double_t *xyz, const Double_t *covI, Double_t sclCovI) const { // calculate the position of the track at PCA to xyz - Double_t t = GetParPCA(xyz,covI); + Double_t t = GetParPCA(xyz,covI,sclCovI); GetPosition(xyzPCA,t); return t; } @@ -444,17 +446,17 @@ void AliITSTPArrayFit::GetResiduals(Double_t *resPCA, const AliTrackPoint *pntCo } //____________________________________________________ -void AliITSTPArrayFit::GetResiduals(Double_t *resPCA, const Double_t *xyz, const Double_t *covI) const +void AliITSTPArrayFit::GetResiduals(Double_t *resPCA, const Double_t *xyz, const Double_t *covI, Double_t sclCovI) const { // calculate the residuals of the track at PCA to xyz - GetPosition(resPCA,xyz,covI); + GetPosition(resPCA,xyz,covI,sclCovI); resPCA[kX] -= xyz[kX]; resPCA[kY] -= xyz[kY]; resPCA[kZ] -= xyz[kZ]; } //____________________________________________________ -Double_t AliITSTPArrayFit::GetParPCALine(const Double_t *xyz, const Double_t *covI) const +Double_t AliITSTPArrayFit::GetParPCALine(const Double_t *xyz, const Double_t *covI/*, Double_t sclCovI*/) const { // get parameter for the point with least weighted distance to the point // @@ -480,13 +482,13 @@ Double_t AliITSTPArrayFit::GetParPCALine(const Double_t *xyz, const Double_t *co } //____________________________________________________ -void AliITSTPArrayFit::GetDResDPosLine(Double_t *dXYZdP, /*const Double_t *xyz,*/ const Double_t *covI) const +void AliITSTPArrayFit::GetDResDPosLine(Double_t *dXYZdP, /*const Double_t *xyz,*/ const Double_t *covI/*,Double_t sclCovI*/) const { // calculate detivative of the PCA residuals vs point position and fill in user provide // array in the format {dXdXp,dY/dXp,dZdXp, ... dXdZp,dYdZp,dZdZp} // Double_t dTdP[3]; - GetDtDPosLine(dTdP, /*xyz,*/ covI); // derivative of the t-param over point position + GetDtDPosLine(dTdP, /*xyz,*/ covI/*,sclCovI*/); // derivative of the t-param over point position // for (int i=3;i--;) { int var = fkAxID[i]; @@ -500,13 +502,13 @@ void AliITSTPArrayFit::GetDResDPosLine(Double_t *dXYZdP, /*const Double_t *xyz,* } //____________________________________________________ -void AliITSTPArrayFit::GetDResDParamsLine(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI) const +void AliITSTPArrayFit::GetDResDParamsLine(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI/*,Double_t sclCovI*/) const { // calculate detivative of the PCA residuals vs line parameters and fill in user provide // array in the format {dXdP0,dYdP0,dZdP0, ... dXdPn,dYdPn,dZdPn} // Double_t dTdP[4]; - Double_t t = GetDtDParamsLine(dTdP, xyz, covI); // derivative of the t-param over line params + Double_t t = GetDtDParamsLine(dTdP, xyz, covI /*,sclCovI*/); // derivative of the t-param over line params // Double_t *curd = dXYZdP + kA0*3; // d/dA0 curd[ fkAxID[kX] ] = fParams[kB0]*dTdP[kA0] + 1.; @@ -642,7 +644,7 @@ void AliITSTPArrayFit::GetDResDPosLine(Double_t *dXYZdP, Int_t ipnt) const AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast)); return; } - GetDResDPosLine(dXYZdP,IsCovIgnored() ? 0 : GetCovI(ipnt)); + GetDResDPosLine(dXYZdP,IsCovIgnored() ? 0 : GetCovI(ipnt)/*,GetCovIScale(ipnt)*/); } //____________________________________________________ @@ -655,7 +657,7 @@ void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, Int_t ipnt) AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast)); return; } - GetDResDParams(dXYZdP, GetPoint(ipnt) , IsCovIgnored() ? 0 : GetCovI(ipnt)); + GetDResDParams(dXYZdP, GetPoint(ipnt) , IsCovIgnored() ? 0 : GetCovI(ipnt),GetCovIScale(ipnt)); } //____________________________________________________ @@ -668,16 +670,16 @@ void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, Int_t ipnt) AliError(Form("Attempt to access the point %d not in the fitted points [%d:%d]",ipnt,fPntFirst,fPntLast)); return; } - GetDResDPos(dXYZdP, GetPoint(ipnt), IsCovIgnored() ? 0 : GetCovI(ipnt)); + GetDResDPos(dXYZdP, GetPoint(ipnt), IsCovIgnored() ? 0 : GetCovI(ipnt), GetCovIScale(ipnt)); } //____________________________________________________ -void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI) +void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI, Double_t sclCovI) { // get residual detivatives over the track parameters for the point with least weighted distance to the point // if (!IsHelix()) { // for the straight line calculate analytically - GetDResDParamsLine(dXYZdP, xyz, covI); + GetDResDParamsLine(dXYZdP, xyz, covI /*,sclCovI*/); return; } // @@ -688,13 +690,13 @@ void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, const Double_t *xyz, con for (int ipar = 5;ipar--;) { double sav = fParams[ipar]; fParams[ipar] -= delta; - GetPosition(xyzVar[0],xyz,covI); + GetPosition(xyzVar[0],xyz,covI,sclCovI); fParams[ipar] += delta/2; - GetPosition(xyzVar[1],xyz,covI); + GetPosition(xyzVar[1],xyz,covI,sclCovI); fParams[ipar] += delta; - GetPosition(xyzVar[2],xyz,covI); + GetPosition(xyzVar[2],xyz,covI,sclCovI); fParams[ipar] += delta/2; - GetPosition(xyzVar[3],xyz,covI); + GetPosition(xyzVar[3],xyz,covI,sclCovI); fParams[ipar] = sav; // restore // double *curd = dXYZdP + 3*ipar; @@ -705,13 +707,13 @@ void AliITSTPArrayFit::GetDResDParams(Double_t *dXYZdP, const Double_t *xyz, con //____________________________________________________ -void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI) +void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, const Double_t *xyz, const Double_t *covI,Double_t sclCovI) const { // get residuals detivative over the point position for the point with least weighted distance to the point // if (!IsHelix()) { // for the straight line calculate analytically - GetDResDPosLine(dXYZdP, /*xyz,*/ covI); + GetDResDPosLine(dXYZdP, /*xyz,*/ covI /*,sclCovI*/); return; } // @@ -723,13 +725,13 @@ void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, const Double_t *xyz, const for (int ipar = 3;ipar--;) { double sav = xyzv[ipar]; xyzv[ipar] -= delta; - GetPosition(xyzVar[0],xyzv,covI); + GetPosition(xyzVar[0],xyzv,covI,sclCovI); xyzv[ipar] += delta/2; - GetPosition(xyzVar[1],xyzv,covI); + GetPosition(xyzVar[1],xyzv,covI,sclCovI); xyzv[ipar] += delta; - GetPosition(xyzVar[2],xyzv,covI); + GetPosition(xyzVar[2],xyzv,covI,sclCovI); xyzv[ipar] += delta/2; - GetPosition(xyzVar[3],xyzv,covI); + GetPosition(xyzVar[3],xyzv,covI,sclCovI); xyzv[ipar] = sav; // restore // double *curd = dXYZdP + 3*ipar; @@ -740,7 +742,7 @@ void AliITSTPArrayFit::GetDResDPos(Double_t *dXYZdP, const Double_t *xyz, const } //________________________________________________________________________________________________________ -Double_t AliITSTPArrayFit::GetParPCAHelix(const Double_t* xyz, const Double_t* covI) const +Double_t AliITSTPArrayFit::GetParPCAHelix(const Double_t* xyz, const Double_t* covI,Double_t sclCovI) const { // find track parameter t (eq.2) corresponding to point of closest approach to xyz // @@ -784,6 +786,7 @@ Double_t AliITSTPArrayFit::GetParPCAHelix(const Double_t* xyz, const Double_t* c dchi2 = dxD*tx + dyD*ty + dzD*tz; ddchi2 = dxDD*tx + dyDD*ty + dxD *ttx + dyD *tty + dzD *ttz; // + if (sclCovI>0) {dchi2 *= sclCovI; ddchi2 *= sclCovI;} } else { // chi2 = dx*dx + dy*dy + dz*dz; @@ -831,9 +834,11 @@ Double_t AliITSTPArrayFit::CalcChi2NDF() const for (int ipnt=fPntFirst;ipnt<=fPntLast;ipnt++) { GetResiduals(dr,ipnt); Double_t* covI = GetCovI(ipnt); - chi2 += dr[kX]*(dr[kX]*covI[ kXX ]+dr[kY]*covI[ kXY ]+dr[kZ]*covI[ kXZ ]) - + dr[kY]*(dr[kX]*covI[ kXY ]+dr[kY]*covI[ kYY ]+dr[kZ]*covI[ kYZ ]) - + dr[kZ]*(dr[kX]*covI[ kXZ ]+dr[kY]*covI[ kYZ ]+dr[kZ]*covI[ kZZ ]); + Double_t chi2p = dr[kX]*(dr[kX]*covI[ kXX ]+dr[kY]*covI[ kXY ]+dr[kZ]*covI[ kXZ ]) + + dr[kY]*(dr[kX]*covI[ kXY ]+dr[kY]*covI[ kYY ]+dr[kZ]*covI[ kYZ ]) + + dr[kZ]*(dr[kX]*covI[ kXZ ]+dr[kY]*covI[ kYZ ]+dr[kZ]*covI[ kZZ ]); + if (covI[kScl]>0) chi2p *= covI[kScl]; // rescaling was requested for this point's errors + chi2 += chi2p; } int ndf = (fPntLast-fPntFirst+1)*3 - GetNParams(); chi2 /= ndf; @@ -962,7 +967,7 @@ Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ) // const float *x=fkPoints->GetX(),*y=fkPoints->GetY(),*z=fkPoints->GetZ(),*cov=fkPoints->GetCov(); // - if (fPntLast>arrU.GetSize()) { + if (fPntLast>=arrU.GetSize()) { arrU.Set(2*fPntLast); arrV.Set(2*fPntLast); arrW.Set(2*fPntLast); @@ -1077,7 +1082,7 @@ Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ) // track direction vector as a - diference between the closest and the next to closest to origin points double v0X = x[minRId1] - x[minRId]; double v0Y = y[minRId1] - y[minRId]; - sqb = (yc*v0X - xc*v0Y)<0 ? -1:1; + sqb = (yc*v0X - xc*v0Y)>0 ? -1:1; } SetCharge( fBz<0 ? -sqb : sqb); } @@ -1194,7 +1199,7 @@ Bool_t AliITSTPArrayFit::FitHelixCrude(Int_t extQ) Double_t normS[3]; // // Positive t-params: along the track direction for negative track, against for positive - Double_t pcur = ptot, ecurr = etot; + // Double_t pcur = ptot, ecurr = etot; for (int ip=fFirstPosT;ipAddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip)); + fParSol->AddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip),GetCovIScale(ip)); } // if (extPT>0) { // add constraint on pt @@ -1682,7 +1687,7 @@ Double_t AliITSTPArrayFit::FitLine() dXYZdLoc[ fkAxID[kY] ] = fParams[kB1]; // dY/dt dXYZdLoc[ fParAxis ] = 1; // - fParSol->AddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip)); + fParSol->AddEquation(dXYZdGlo,dXYZdLoc,xyzRes,GetCovI(ip),GetCovIScale(ip)); } // if (!fParSol->Solve()) { AliError("Failed to fit line"); return -1; } @@ -1857,7 +1862,7 @@ void AliITSTPArrayFit::BuildMaterialLUT(Int_t ntri) for (int i=ntri;i--;) { // if (active) { - int ssID = sID -1 - AliGeomManager::LayerSize(actLrID)*gRandom->Rndm(); + int ssID = sID -1 - (int)(AliGeomManager::LayerSize(actLrID)*gRandom->Rndm()); pg1[0] = pg2[0] = (gRandom->Rndm()-0.5)*tpars[0] + shift; // local X pg2[0] -= 2*shift; pg1[1] = tpars[2]; @@ -1926,10 +1931,13 @@ Bool_t AliITSTPArrayFit::CalcErrorMatrix() GetDResDParams(&dres[0][0],ip); Double_t* covI = GetCovI(ip); for (int i=npar;i--;) - for (int j=i+1;j--;) - cv(i,j) += dres[i][kX]*(dres[j][kX]*covI[ kXX ] + dres[j][kY]*covI[ kXY ] + dres[j][kZ]*covI[ kXZ ]) - + dres[i][kY]*(dres[j][kX]*covI[ kXY ] + dres[j][kY]*covI[ kYY ] + dres[j][kZ]*covI[ kYZ ]) - + dres[i][kZ]*(dres[j][kX]*covI[ kXZ ] + dres[j][kY]*covI[ kYZ ] + dres[j][kZ]*covI[ kZZ ]); + for (int j=i+1;j--;) { + double cvadd = dres[i][kX]*(dres[j][kX]*covI[ kXX ] + dres[j][kY]*covI[ kXY ] + dres[j][kZ]*covI[ kXZ ]) + + dres[i][kY]*(dres[j][kX]*covI[ kXY ] + dres[j][kY]*covI[ kYY ] + dres[j][kZ]*covI[ kYZ ]) + + dres[i][kZ]*(dres[j][kX]*covI[ kXZ ] + dres[j][kY]*covI[ kYZ ] + dres[j][kZ]*covI[ kZZ ]); + if (covI[kScl]>0) cvadd *= covI[kScl]; + cv(i,j) += cvadd; + } } cv.SetSizeUsed(npar); if (cv.InvertChol()) { @@ -1948,8 +1956,8 @@ Double_t AliITSTPArrayFit::CalcParPCA(Int_t ipnt) const // get parameter for the point with least weighted distance to the point const double *xyz = GetPoint(ipnt); const double *covI = GetCovI(ipnt); - if (IsHelix()) return GetParPCAHelix(xyz,covI); - else return GetParPCALine(xyz,covI); + if (IsHelix()) return GetParPCAHelix(xyz,covI,covI[kScl]); + else return GetParPCALine(xyz,covI/*,covI[kScl]*/); } //____________________________________________________