X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliTrackFitterRieman.cxx;h=f4e58101ca2f15b1b9eb2802f4f94400ef07e421;hb=55eb2b56f0af4ddd07785a86d9ac8f106c5fd835;hp=c2aad546c140902e6f372f36f981dacacf38da54;hpb=186e6cedaee51bdc68cc3ed968bea7ba808290be;p=u%2Fmrichter%2FAliRoot.git

diff --git a/STEER/AliTrackFitterRieman.cxx b/STEER/AliTrackFitterRieman.cxx
index c2aad546c14..f4e58101ca2 100644
--- a/STEER/AliTrackFitterRieman.cxx
+++ b/STEER/AliTrackFitterRieman.cxx
@@ -29,58 +29,68 @@
 //
 //////////////////////////////////////////////////////////////////////////////
 
-#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
+  fDebugStream(new TTreeSRedirector("RiemanAlignDebug.root"))
 {
   //
   // default constructor
   //
-  fAlpha = 0.;
-  fNUsed = 0;
-  fConv = kFALSE;
-  fMaxDelta = 3;
-  fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
-  fRieman = new AliRieman(10000);  // allocate rieman
 }
 
 
 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
+  fDebugStream(0)
 {
   //
   // Constructor
   //
-  fAlpha = 0.;
-  fNUsed = 0;
-  fConv = kFALSE;
-  fMaxDelta = 3;
   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))),
+  fDebugStream(0)
 {
   //
   // copy constructor
   //
-  fAlpha = rieman.fAlpha;
-  fNUsed = rieman.fNUsed;
-  fConv = rieman.fConv;
-  fMaxDelta = rieman.fMaxDelta;
-  fRieman = new AliRieman(*(rieman.fRieman));
+  if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
 }
 
 //_____________________________________________________________________________
@@ -92,11 +102,14 @@ AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRiema
   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));
+  fDebugStream = 0;
+  if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
   return *this;
 }
 
@@ -122,8 +135,8 @@ void AliTrackFitterRieman::Reset()
 }
 
 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
@@ -141,10 +154,34 @@ Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
   // 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();
   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);    
@@ -156,7 +193,6 @@ Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
       "fPoints.="<<fPoints<<   // input points
       "\n";
   }
-  if (npoints < 3) return kFALSE;
 
   Bool_t isAlphaCalc = kFALSE;
   AliTrackPoint p,plocal;
@@ -178,9 +214,9 @@ Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
 	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();
@@ -252,30 +288,35 @@ Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
       }
     }  
   
-  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
@@ -284,18 +325,24 @@ Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
     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";
@@ -320,7 +367,7 @@ void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy,
 
 
 
-void AliTrackFitterRieman::Update(){
+Bool_t AliTrackFitterRieman::Update(){
   //
   // 
   //
@@ -335,8 +382,30 @@ void 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
 {
@@ -395,6 +464,12 @@ Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) c
   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+
@@ -411,13 +486,73 @@ Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) c
   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);
+}