//
//////////////////////////////////////////////////////////////////////////////
-#include "TMatrixDSym.h"
-#include "TMatrixD.h"
-#include "TArrayI.h"
-#include "AliTrackFitterRieman.h"
+#include <TArrayI.h>
+#include <TLinearFitter.h>
+#include <TMath.h>
+#include <TMatrixD.h>
+#include <TMatrixDSym.h>
+#include <TRandom.h>
+#include <TTreeStream.h>
+
#include "AliLog.h"
-#include "TTreeStream.h"
-#include "AliRieman.h"
#include "AliLog.h"
+#include "AliRieman.h"
+#include "AliTrackFitterRieman.h"
ClassImp(AliTrackFitterRieman)
AliTrackFitterRieman::AliTrackFitterRieman():
- AliTrackFitter()
+ AliTrackFitter(),
+ fBCorrection(kFALSE),
+ fAlpha(0.),
+ fNUsed(0),
+ fConv(kFALSE),
+ fMaxDelta(3),
+ fRieman(new AliRieman(10000)), // allocate rieman
+ fMinPointRadius(2.),
+ fMaxPointRadius(500.),
+ fDebugStream(new TTreeSRedirector("RiemanAlignDebug.root"))
{
//
// default constructor
//
- fAlpha = 0.;
- fNUsed = 0;
- fConv = kFALSE;
- fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
- fRieman = new AliRieman(10000); // allocate rieman
+ fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
+ fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
}
AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
- AliTrackFitter(array,owner)
+ AliTrackFitter(array,owner),
+ fBCorrection(kFALSE),
+ fAlpha(0.),
+ fNUsed(0),
+ fConv(kFALSE),
+ fMaxDelta(3),
+ fRieman(new AliRieman(10000)), //allocate rieman
+ fMinPointRadius(2.),
+ fMaxPointRadius(500.),
+ fDebugStream(0)
{
//
// Constructor
//
- fAlpha = 0.;
- fNUsed = 0;
- fConv = kFALSE;
+ fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
+ fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
- fRieman = new AliRieman(10000); //allocate rieman
}
AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
- AliTrackFitter(rieman)
+ AliTrackFitter(rieman),
+ fBCorrection(rieman.fBCorrection),
+ fAlpha(rieman.fAlpha),
+ fNUsed(rieman.fNUsed),
+ fConv(rieman.fConv),
+ fMaxDelta(rieman.fMaxDelta),
+ fRieman(new AliRieman(*(rieman.fRieman))),
+ fMinPointRadius(rieman.fMinPointRadius),
+ fMaxPointRadius(rieman.fMaxPointRadius),
+ fDebugStream(0)
{
//
// copy constructor
//
- fAlpha = rieman.fAlpha;
- fNUsed = rieman.fNUsed;
- fConv = rieman.fConv;
- fRieman = new AliRieman(*(rieman.fRieman));
+ fCorrY[0] = fCorrY[1] = fCorrY[2] = fCorrY[3] = 0;
+ fCorrZ[0] = fCorrZ[1] = fCorrZ[2] = fCorrZ[3] = 0;
+ if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
}
//_____________________________________________________________________________
if(this==&rieman) return *this;
((AliTrackFitter *)this)->operator=(rieman);
+ fBCorrection = rieman.fBCorrection;
fAlpha = rieman.fAlpha;
fNUsed = rieman.fNUsed;
fConv = rieman.fConv;
+ fMaxDelta = rieman.fMaxDelta;
fRieman = new AliRieman(*(rieman.fRieman));
+ fMinPointRadius = rieman.fMinPointRadius;
+ fMaxPointRadius = rieman.fMaxPointRadius;
+ fDebugStream = 0;
+ if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
return *this;
}
}
Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
- AliAlignObj::ELayerID layerRangeMin,
- AliAlignObj::ELayerID layerRangeMax)
+ AliGeomManager::ELayerID layerRangeMin,
+ AliGeomManager::ELayerID layerRangeMax)
{
// Fit the track points. The method takes as an input
// the set of id's (volids) of the volumes in which
// at the space-point reference plane. The reference plane is
// found using the covariance matrix of the point
// (assuming sigma(x)=0 at the reference coordinate system.
+ Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
+
+ // Float_t debugRatio = 1./(1.+debugLevel);
+ Float_t debugRatio = debugLevel? 1.0/debugLevel : 1.0;
- Reset();
+ if (!fPoints) {
+ AliError("Track points array not available! Exiting...");
+ return kFALSE;
+ }
Int_t npoints = fPoints->GetNPoints();
- if (fPoints && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1){
+ if ( npoints<fMinNPoints) return kFALSE;
+ //
+ // fast count points
+ if (volIdsFit != 0x0) {
+ Int_t countFit = 0;
+ Int_t countPoint = 0;
+ for (Int_t ipoint = 0; ipoint < npoints; ipoint++) {
+ if (FindVolId(volIds,fPoints->GetVolumeID()[ipoint]))
+ countPoint++;
+ if (volIdsFit != 0x0) {
+ if (FindVolId(volIdsFit,fPoints->GetVolumeID()[ipoint]))
+ countFit++;
+ }
+ }
+ if (countPoint==0) return kFALSE;
+ if ((countFit<fMinNPoints) && (volIdsFit != 0x0)) return kFALSE;
+ }
+ //
+
+ Reset();
+
+ if (fPoints && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1&& gRandom->Rndm()<debugRatio){
Int_t nVol = volIds->GetSize();
Int_t nVolFit = volIdsFit->GetSize();
Int_t volId = volIds->At(0);
"fPoints.="<<fPoints<< // input points
"\n";
}
- if (npoints < 3) return kFALSE;
Bool_t isAlphaCalc = kFALSE;
AliTrackPoint p,plocal;
if (!FindVolId(volIdsFit,iVolId)) continue;
}
else {
- if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
- iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,
- AliAlignObj::LayerSize(layerRangeMax))) continue;
+ if (iVolId < AliGeomManager::LayerToVolUID(layerRangeMin,0) ||
+ iVolId > AliGeomManager::LayerToVolUID(layerRangeMax,
+ AliGeomManager::LayerSize(layerRangeMax))) continue;
}
if (!isAlphaCalc) {
fAlpha = p.GetAngle();
isAlphaCalc = kTRUE;
}
plocal = p.Rotate(fAlpha);
- if (TMath::Abs(plocal.GetX())>500 || TMath::Abs(plocal.GetX())<2){
+ if (TMath::Abs(plocal.GetX())>fMaxPointRadius || TMath::Abs(plocal.GetX())<fMinPointRadius || plocal.GetCov()[3]<=0 ||plocal.GetCov()[5]<=0 ){
printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
p.Dump();
plocal.Dump();
printf("Problematic point\n");
printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
+ }else{
+ AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
+ TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
}
- AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
- TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
// fNUsed++; AddPoint should be responsible
}
{
Int_t index = pindex[ipoint];
fPoints->GetPoint(p,index);
- if (GetPCA(p,p2)) {
+ if (GetPCA(p,p2) && (
+ TMath::Abs(p.GetX()-p2.GetX())<fMaxDelta &&
+ TMath::Abs(p.GetY()-p2.GetY())<fMaxDelta &&
+ TMath::Abs(p.GetZ()-p2.GetZ())<fMaxDelta
+ )) {
Float_t xyz[3],xyz2[3];
p.GetXYZ(xyz); p2.GetXYZ(xyz2);
- // printf("residuals %f %d %d %f %f %f %f %f %f\n",fChi2,fNUsed,fConv,xyz[0],xyz[1],xyz[2],xyz2[0]-xyz[0],xyz2[1]-xyz[1],xyz2[2]-xyz[2]);
+ // printf("residuals %f %d %d %f %f %f %f %f %f\n",fChi2,fNUsed,fConv,xyz[0],xyz[1],xyz[2],xyz2[0]-xyz[0],xyz2[1]-xyz[1],xyz2[2]-xyz[2]);
fPVolId->AddPoint(ipoint,&p);
fPTrack->AddPoint(ipoint,&p2);
}else{
}
}
- if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0){
+ if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
//
// Debug Info
//
AliTrackPointArray *lPVolId = new AliTrackPointArray(npVolId);
AliTrackPointArray *lPTrack = new AliTrackPointArray(npVolId);
+ AliTrackPointArray *lPTrackE = new AliTrackPointArray(npVolId);
AliRieman * residual = fRieman->MakeResiduals();
- AliTrackPoint p2local;
for (Int_t ipoint = 0; ipoint < npVolId; ipoint++){
- Int_t index = pindex[ipoint];
- fPoints->GetPoint(p,index);
- if (GetPCA(p,p2)) {
- plocal = p.Rotate(fAlpha);
- // plocal.Rotate(fAlpha);
- Float_t xyz[3],xyz2[3];
- p2local = plocal;
- plocal.GetXYZ(xyz);
- xyz2[0] = xyz[0];
- xyz2[1] = GetYat(xyz[0]);
- xyz2[2] = GetZat(xyz[0]);
- p2local.SetXYZ(xyz2);
- lPVolId->AddPoint(ipoint,&plocal);
- lPTrack->AddPoint(ipoint,&p2local);
- }
+ AliTrackPoint p0, p0local;
+ AliTrackPoint pFit, pFitlocal, pFitLocalE;
+ fPVolId->GetPoint(p0,ipoint);
+ Float_t lAngle = p0.GetAngle();
+ p0local= p0.MasterToLocal();
+ fPTrack->GetPoint(pFit,ipoint);
+ pFitlocal= pFit.Rotate(lAngle);
+ //
+ Float_t xyz[3], cov[6];
+ xyz[0] = pFitlocal.GetX();
+ xyz[1] = pFitlocal.GetY();
+ xyz[2] = pFitlocal.GetZ();
+ for (Int_t icov=0; icov<6; icov++) cov[icov]=0;
+ cov[3] = GetErrY2at(xyz[0]);
+ cov[5] = GetErrZ2at(xyz[0]);
+ pFitLocalE.SetXYZ(xyz,cov);
+ //
+ lPVolId->AddPoint(ipoint,&p0local);
+ lPTrack->AddPoint(ipoint,&pFitlocal);
+ lPTrackE->AddPoint(ipoint,&pFitLocalE);
}
//
// debug info
Int_t nVolFit = volIdsFit->GetSize();
Int_t volId = volIds->At(0);
Int_t modId =0;
- Int_t layer = AliAlignObj::VolUIDToLayer(volId,modId);
+ Int_t layer = AliGeomManager::VolUIDToLayer(volId,modId);
+ Int_t volIdFit = volIdsFit->At(0);
+ Int_t modIdFit =0;
+ Int_t layerFit = AliGeomManager::VolUIDToLayer(volIdFit,modIdFit);
(*fDebugStream)<<"Fit"<<
"VolId="<<volId<< // volume ID
"Layer="<<layer<< // layer ID
"Module="<<modId<< // module ID
+ "LayerFit="<<layerFit<< // layer ID fit
+ "ModuleFit="<<modIdFit<< // module ID fit
"NVol="<<nVol<< // number of volumes
"NvolFit="<<nVolFit<< // number of volumes to fit
"Points0.="<<fPVolId<< // original points
"Points1.="<<fPTrack<< // fitted points
"LPoints0.="<<lPVolId<< // original points - local frame
"LPoints1.="<<lPTrack<< // fitted points - local frame
+ "LPointsE.="<<lPTrackE<< // fitted points with ext error - local frame
"Rieman.="<<this<< // original rieman fit
"Res.="<<residual<< // residuals of rieman fit
"\n";
-void AliTrackFitterRieman::Update(){
+Bool_t AliTrackFitterRieman::Update(){
//
//
//
fNUsed = fRieman->GetN();
fConv = kTRUE;
}
+ //
+ //
+ //
+ TLinearFitter fitY(3,"pol2");
+ TLinearFitter fitZ(3,"pol2");
+ for (Int_t ip=0; ip<fRieman->GetN();ip++){
+ Double_t x = fRieman->GetX()[ip];
+ fitY.AddPoint(&x,fRieman->GetY()[ip]-fRieman->GetYat(x),1);
+ fitZ.AddPoint(&x,fRieman->GetZ()[ip]-fRieman->GetZat(x),1);
+ }
+ fitY.Eval();
+ fitZ.Eval();
+ for (Int_t iparam=0; iparam<3; iparam++){
+ fCorrY[iparam]=fitY.GetParameter(iparam);
+ fCorrZ[iparam]=fitZ.GetParameter(iparam);
+ }
+ fCorrY[3]=fitY.GetChisquare()/Float_t(fRieman->GetN()-3);
+ fCorrZ[3]=fitZ.GetChisquare()/Float_t(fRieman->GetN()-3);
+
+ return kTRUE;
}
+
+
//_____________________________________________________________________________
Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
{
cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx;
cov[5] = sigmax*dzdx*dzdx;
+ Double_t sigmay2 = GetErrY2at(x2);
+ Double_t sigmaz2 = GetErrZ2at(x2);
+ cov[3] += sigmay2;
+ cov[5] += sigmaz2;
+
+
Float_t newcov[6];
newcov[0] = cov[0]*cos*cos-
2*cov[1]*sin*cos+
newcov[5] = cov[5];
p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
- if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>1){
+ Int_t debugLevel = AliLog::GetDebugLevel("","AliTrackFitterRieman");
+ Float_t debugRatio = 1./(1.+debugLevel);
+ if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0 && gRandom->Rndm()<debugRatio){
AliTrackPoint lp0(p);
AliTrackPoint lp2(p2);
- if (0)(*fDebugStream)<<"PCA"<<
- "P0.="<<&lp0<<
+ AliTrackPoint localp0(p);
+ AliTrackPoint localp2(p2);
+ Float_t lAngle = lp0.GetAngle();
+ localp0 = localp0.Rotate(lAngle);
+ localp2 = localp2.Rotate(lAngle);
+
+ (*fDebugStream)<<"PCA"<<
+ "P0.="<<&lp0<< //global position
"P2.="<<&lp2<<
+ "LP0.="<<&localp0<< //local position
+ "LP2.="<<&localp2<<
"\n";
}
return kTRUE;
}
+
+Double_t AliTrackFitterRieman::GetYat(Double_t x) const {
+ //
+ // get y position at given point
+ //
+ Double_t correction=0;
+ if (fBCorrection){ // systematic effect correction
+ correction = fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
+ }
+ return fRieman->GetYat(x)+correction;
+}
+
+Double_t AliTrackFitterRieman::GetZat(Double_t x) const {
+ //
+ // get z position at given point
+ //
+ Double_t correction=0;
+ if (fBCorrection){ // systematic effect correction
+ correction = fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
+ }
+ return fRieman->GetZat(x)+correction;
+}
+
+Double_t AliTrackFitterRieman::GetErrY2at(Double_t x) const {
+ //
+ // get estimate of extrapolation error
+ //
+ Double_t error = fRieman->GetErrY(x);
+ Double_t correction=0;
+ if (fBCorrection){ // everestimate error due systematic effect
+ error *=fCorrY[3];
+ correction = fCorrY[0]+fCorrY[1]*x +fCorrY[2]*x*x;
+ correction *=correction;
+ }
+ return TMath::Sqrt(error+correction);
+}
+
+Double_t AliTrackFitterRieman::GetErrZ2at(Double_t x) const {
+ //
+ // get estimate of extrapolation error
+ //
+ Double_t error = fRieman->GetErrZ(x)*fCorrZ[3];
+ Double_t correction=0;
+ if (fBCorrection){
+ error *= fCorrZ[3];
+ correction = fCorrZ[0]+fCorrZ[1]*x +fCorrZ[2]*x*x;
+ correction*= correction;
+ }
+ return TMath::Sqrt(error+correction);
+}
+
+void AliTrackFitterRieman::SetParam(Int_t i, Double_t par) {
+ if (i<0 || i>5) return;
+ fParams[i]=par;
+ fRieman->GetParam()[i]=par;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff