New class for linear reasiduals minimization based on root TLinearFitter. Possibility...
authorcvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 4 May 2006 12:32:49 +0000 (12:32 +0000)
committercvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 4 May 2006 12:32:49 +0000 (12:32 +0000)
STEER/AliTrackResidualsLinear.cxx [new file with mode: 0644]
STEER/AliTrackResidualsLinear.h [new file with mode: 0644]
STEER/STEERLinkDef.h
STEER/libSTEER.pkg

diff --git a/STEER/AliTrackResidualsLinear.cxx b/STEER/AliTrackResidualsLinear.cxx
new file mode 100644 (file)
index 0000000..ff844ba
--- /dev/null
@@ -0,0 +1,271 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+//   Implementation of the derived class for track residuals
+//   based on linear chi2 minimization 
+//  The minimization relies on the fact that the alignment parameters     
+//   (angles and translations) are small.                                  
+//   TLinearFitter used for minimization
+//   Possibility to fix Paramaters
+//   FixParameter()ReleaseParameter();
+//   Possibility to define fraction of outliers to be skipped
+//
+//   marian.ivanov@cern.ch
+//
+//-----------------------------------------------------------------
+
+#include <TMinuit.h>
+#include <TGeoMatrix.h>
+
+#include "AliLog.h"
+#include "AliAlignObj.h"
+#include "AliTrackPointArray.h"
+#include "AliTrackResidualsLinear.h"
+#include "AliAlignObj.h"
+#include "TLinearFitter.h"
+#include  "TDecompSVD.h"
+
+ClassImp(AliTrackResidualsLinear)
+
+//______________________________________________________________________________
+AliTrackResidualsLinear::AliTrackResidualsLinear():AliTrackResiduals()
+{
+  // Default constructor
+  fFitter = 0;
+  fFraction =-1;
+  fChi2Orig =0;
+  for (Int_t ipar=0; ipar<6; ipar++){
+    fBFixed[ipar] = kFALSE;
+    fFixed[ipar]  = 0;;
+    fParams[ipar]  = 0;
+  }  
+}
+
+//______________________________________________________________________________
+AliTrackResidualsLinear::AliTrackResidualsLinear(Int_t ntracks):
+  AliTrackResiduals(ntracks)
+{
+  // Constructor
+  fFitter = new TLinearFitter(6,"hyp6");
+  fFraction =-1;
+  fChi2Orig =0;
+  for (Int_t ipar=0; ipar<6; ipar++){
+    fBFixed[ipar] = kFALSE;
+    fFixed[ipar]  = 0;
+    fParams[ipar]  = 0;
+  }
+}
+//______________________________________________________________________________
+AliTrackResidualsLinear::AliTrackResidualsLinear(const AliTrackResidualsLinear &res):
+  AliTrackResiduals(res)
+{
+  // Copy constructor
+  //..
+  fFitter = new TLinearFitter(*(res.fFitter));
+  for (Int_t ipar=0; ipar<6; ipar++){
+    fBFixed[ipar]  = res.fBFixed[ipar];
+    fFixed[ipar]   = res.fFixed[ipar];
+    fParams[ipar]  = res.fParams[ipar];
+  }
+  fChi2Orig = res.fChi2Orig;
+}
+
+//______________________________________________________________________________
+AliTrackResidualsLinear &AliTrackResidualsLinear::operator= (const AliTrackResidualsLinear& res)
+{
+  // Assignment operator
+ ((AliTrackResiduals *)this)->operator=(res);
+ return *this;
+}
+//______________________________________________________________________________
+AliTrackResidualsLinear::~AliTrackResidualsLinear()
+{
+  //
+  //
+  //
+  delete fFitter;
+}
+
+
+//______________________________________________________________________________
+Bool_t AliTrackResidualsLinear::Minimize()
+{
+  // Implementation of fast linear Chi2 minimizer
+  // based on TLinear fitter
+  //
+  if (!fFitter) fFitter = new TLinearFitter(6,"hyp6");
+  fFitter->StoreData(kTRUE);
+  fFitter->ClearPoints();
+  fChi2Orig = 0;
+  AliTrackPoint p1,p2;
+  for (Int_t itrack = 0; itrack < fLast; itrack++) {
+    if (!fVolArray[itrack] || !fTrackArray[itrack]) continue;
+    for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) {
+      fVolArray[itrack]->GetPoint(p1,ipoint);
+      fTrackArray[itrack]->GetPoint(p2,ipoint);
+      AddPoints(p1,p2);
+    }
+  }
+  Bool_t isOK = Update();
+  if (!isOK) return isOK;
+  //
+  TGeoHMatrix  matrix;
+  fAlignObj->GetMatrix(matrix);
+  return isOK;
+}
+
+
+
+//______________________________________________________________________________
+void AliTrackResidualsLinear::AddPoints(AliTrackPoint &p, AliTrackPoint &pprime)
+{
+  //
+  // add points to linear fitter - option with correlation betwee measurement in different dimensions
+  // p1      - local point
+  // pprime  - track extrapolation point
+  //
+  Float_t xyz[3],xyzp[3];
+  Float_t cov[6],covp[6];
+  p.GetXYZ(xyz,cov); pprime.GetXYZ(xyzp,covp);
+  //
+  //
+  TMatrixD mcov(3,3);   // local point covariance
+  mcov(0,0) = cov[0]; mcov(0,1) = cov[1]; mcov(0,2) = cov[2];
+  mcov(1,0) = cov[1]; mcov(1,1) = cov[3]; mcov(1,2) = cov[4];
+  mcov(2,0) = cov[2]; mcov(2,1) = cov[4]; mcov(2,2) = cov[5];
+  TMatrixD mcovp(3,3); //  extrapolation point covariance  
+  mcovp(0,0) = covp[0]; mcovp(0,1) = covp[1]; mcovp(0,2) = covp[2];
+  mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4];
+  mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5];
+  mcov+=mcovp;
+  mcov.Invert();
+  if (!mcov.IsValid()) return; 
+  TMatrixD mcovBack = mcov;  // for debug purposes
+  //
+  // decompose matrix
+  //
+  TDecompSVD svd(mcov);              // mcov  = svd.fU * covDiagonal * svd.fV.Invert   
+  if (!svd.Decompose()) return;      // decomposition failed
+  TMatrixD   matrixV = svd.GetV();   // transformation matrix to diagonalize covariance matrix
+  Double_t   covDiagonal[3] = {svd.GetSig()[0],svd.GetSig()[1],svd.GetSig()[2]};    // diagonalized covariance matrix
+  //
+  // residual vector 
+  TMatrixD  deltaR(3,1);
+  deltaR(0,0) = (xyzp[0]-xyz[0]); 
+  deltaR(1,0) = (xyzp[1]-xyz[1]);
+  deltaR(2,0) = (xyzp[2]-xyz[2]);   
+  //
+  // parametrization matrix
+  //
+  TMatrixD        mparam(3,6);
+  mparam(0,0) = 1;      mparam(1,0) = 0;       mparam(2,0) = 0;            // xshift
+  mparam(0,1) = 0;      mparam(1,1) = 1;       mparam(2,1) = 0;            // yshift
+  mparam(0,2) = 0;      mparam(1,2) = 0;       mparam(2,2) = 1;            // zshift
+  mparam(0,3) = 0;      mparam(1,3) =-xyz[2];  mparam(2,3) = xyz[1];       // x rotation
+  mparam(0,4) = xyz[2]; mparam(1,4) = 0;       mparam(2,4) =-xyz[0];       // y rotation
+  mparam(0,5) =-xyz[1]; mparam(1,5) = xyz[0];  mparam(2,5) = 0;            // z rotation
+  //
+  
+  TMatrixD  deltaT(matrixV, TMatrixD::kTransposeMult, deltaR);   // tranformed delta
+  TMatrixD  mparamT(matrixV,TMatrixD::kTransposeMult, mparam);   // tranformed linear transformation
+  if (0){    
+    //
+    // debug part
+    //
+    //   covDiag = U^-1 * mcov * V      -- diagonalization of covariance matrix
+
+    TMatrixD   matrixU = svd.GetU();                      // transformation matrix to diagonalize covariance matrix
+    TMatrixD   matrixUI= svd.GetU(); 
+    matrixUI.Invert();
+    //
+    TMatrixD   test0   = matrixUI*matrixV;                // test matrix - should be unit matrix
+    TMatrixD   test1   = matrixUI*mcovBack*matrixV;       // test matrix - diagonal - should be diagonal with covDiagonal on diag
+    TMatrixD   test2   = matrixU.T()*matrixV;             // test ortogonality - shoul be unit
+    printf("Test matrix 2 - should be unit\n");
+    test2.Print();
+    printf("Test matrix 0 - should be unit\n"); 
+    test0.Print();
+    printf("Test matrix 1 - should be diagonal\n"); 
+    test1.Print();
+    printf("Diagonal matrix\n"); 
+    svd.GetSig().Print();
+    printf("Original param matrix\n"); 
+    mparam.Print();
+    printf("Rotated  param matrix\n"); 
+    mparamT.Print();
+    //
+  }
+  //
+  for (Int_t idim = 0; idim<3; idim++){
+    Double_t yf;     // input points to fit in TLinear fitter
+    Double_t xf[6];  // input points to fit
+    yf = deltaT(idim,0);
+    for (Int_t ipar =0; ipar<6; ipar++) xf[ipar] = mparamT(idim,ipar); 
+    if (covDiagonal[idim]>0.){
+      fFitter->AddPoint(xf,yf, TMath::Sqrt(1/covDiagonal[idim]));
+    }
+    // accumulate chi2
+    fChi2Orig += (yf*yf)*covDiagonal[idim];          
+  }
+  fNdf +=3;
+}
+
+//______________________________________________________________________________
+Bool_t AliTrackResidualsLinear::Update()
+{
+  // Find the alignment parameters
+  // using TLinear fitter  + fill data containers
+  // 
+  //
+  fFitter->Eval();
+  //
+  // TLinear fitter put as first parameter offset - fixing parameter shifted by one
+  //
+  fFitter->FixParameter(0);
+  for (Int_t ipar =0; ipar<6; ipar++){
+    if (fBFixed[ipar])  fFitter->FixParameter(ipar+1,fFixed[ipar]);
+  }
+  if (fFraction>0.5) {
+    fFitter->EvalRobust(fFraction);
+  }else{
+    fFitter->Eval();
+  }
+  //
+  fFitter->ReleaseParameter(0);
+  for (Int_t ipar=0; ipar<7; ipar++) {
+    if (fBFixed[ipar])  fFitter->ReleaseParameter(ipar+1);
+  }
+    
+
+  //
+  fChi2 = fFitter->GetChisquare();
+  fNdf -= 6;
+  TVectorD vector(7);
+  fFitter->GetParameters(vector);
+  fParams[0] = vector[1];
+  fParams[1] = vector[2];
+  fParams[2] = vector[3];  
+  fParams[3] = vector[4];
+  fParams[4] = vector[5];
+  fParams[5] = vector[6];
+  //
+  fAlignObj->SetPars(fParams[0], fParams[1], fParams[2],
+                    TMath::RadToDeg()*fParams[3],
+                    TMath::RadToDeg()*fParams[4],
+                    TMath::RadToDeg()*fParams[5]);
+  return kTRUE;
+}
diff --git a/STEER/AliTrackResidualsLinear.h b/STEER/AliTrackResidualsLinear.h
new file mode 100644 (file)
index 0000000..227c249
--- /dev/null
@@ -0,0 +1,42 @@
+#ifndef ALITRACKRESIDUALSLINEAR_H
+#define ALITRACKRESIDUALSLINEAR_H
+
+//************************************************************************
+// AliTrackResidualsLinear: derived class (from AliTrackResiduals) which   *
+// implements a simple linear minimization of the track residuals chi2.  *
+// The minimization relies on the fact that the alignment parameters     *
+// (angles and translations) are small.                                  *
+//                                                                       *
+//                                                                       *
+//************************************************************************
+
+#include "AliAlignObj.h"
+#include "AliTrackResiduals.h"
+class TLinearFitter;
+
+class AliTrackResidualsLinear : public AliTrackResiduals {
+
+ public:
+  AliTrackResidualsLinear();
+  AliTrackResidualsLinear(Int_t ntracks);
+  AliTrackResidualsLinear(const AliTrackResidualsLinear &res);
+  AliTrackResidualsLinear& operator= (const AliTrackResidualsLinear& res);
+  virtual ~AliTrackResidualsLinear();
+  Bool_t Minimize();
+  void   SetRobust(Float_t fraction){fFraction=fraction;}
+  void   FixParameter(Int_t par, Float_t value) {fBFixed[par]=kTRUE; fFixed[par]= value;}
+  void   ReleaseParameter(Int_t par) {fBFixed[par]=kFALSE;}
+ protected:
+  Bool_t Update();
+  void   AddPoints(AliTrackPoint &p, AliTrackPoint &pprime);
+  TLinearFitter *fFitter;           // !Linear fitter
+  Float_t fFraction;                // fraction of points for robust fitter if less 0 - no robust fitter invoked
+  Float_t fParams[6];               // resulting parameters 
+  Float_t fFixed[6];                // the fixed values of parameters 
+  Float_t fBFixed[6];               // the flag for fixing parameter
+  Double_t fChi2Orig;               // original chi2 
+  ClassDef(AliTrackResidualsLinear,1)
+
+};
+
+#endif
index 1c3a8802ae3e8b17bf1fe993e130edb410120249..fbc4dd5a7db195e003aff283d8a2830ef8b13f4e 100644 (file)
 #pragma link C++ class AliTrackResiduals+;
 #pragma link C++ class AliTrackResidualsChi2+;
 #pragma link C++ class AliTrackResidualsFast+;
 #pragma link C++ class AliTrackResiduals+;
 #pragma link C++ class AliTrackResidualsChi2+;
 #pragma link C++ class AliTrackResidualsFast+;
+#pragma link C++ class AliTrackResidualsLinear+;
 #pragma link C++ class AliAlignmentTracks+;
 
 #pragma link C++ class  TTreeDataElement+;
 #pragma link C++ class AliAlignmentTracks+;
 
 #pragma link C++ class  TTreeDataElement+;
index 2fc300577a0c28521502e4d4e6b6464af6afe12a..98068dc9b3eb2d7e0dbe23c81075c491a76c2ed2 100644 (file)
@@ -30,7 +30,7 @@ AliTriggerCondition.cxx \
 AliTriggerDescriptor.cxx \
 AliCentralTrigger.cxx AliRieman.cxx\
 AliTrackFitter.cxx AliTrackFitterRieman.cxx AliTrackFitterStraight.cxx\
 AliTriggerDescriptor.cxx \
 AliCentralTrigger.cxx AliRieman.cxx\
 AliTrackFitter.cxx AliTrackFitterRieman.cxx AliTrackFitterStraight.cxx\
-AliTrackResiduals.cxx AliTrackResidualsChi2.cxx AliTrackResidualsFast.cxx\
+AliTrackResiduals.cxx AliTrackResidualsChi2.cxx AliTrackResidualsFast.cxx AliTrackResidualsLinear.cxx\
 AliAlignmentTracks.cxx \
 AliExpression.cxx \
 AliCTPRawData.cxx AliCTPRawStream.cxx
 AliAlignmentTracks.cxx \
 AliExpression.cxx \
 AliCTPRawData.cxx AliCTPRawStream.cxx