X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliStrLine.cxx;h=099df02f25aa10050175a7eb37b1be545819e82b;hb=583e63b8cc5bcb315a36593ccc48f486c2907bbf;hp=8eea1aaaaa1c688dda1e041436545050e81f72f6;hpb=2c9641ee47d1932a87a80f71295e0cda4f445c7e;p=u%2Fmrichter%2FAliRoot.git

diff --git a/STEER/AliStrLine.cxx b/STEER/AliStrLine.cxx
index 8eea1aaaaa1..099df02f25a 100644
--- a/STEER/AliStrLine.cxx
+++ b/STEER/AliStrLine.cxx
@@ -12,14 +12,19 @@
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
+
+/* $Id$ */
+
 ///////////////////////////////////////////////////////////////////
 //                                                               //
 // A straight line is coded as a point (3 Double_t) and           //
 // 3 direction cosines                                           //
 //                                                               //
 ///////////////////////////////////////////////////////////////////
+
 #include <Riostream.h>
-#include <TTree.h>
+#include <TMath.h>
+
 #include "AliStrLine.h"
 
 ClassImp(AliStrLine)
@@ -27,21 +32,23 @@ ClassImp(AliStrLine)
 //________________________________________________________
 AliStrLine::AliStrLine() :
   TObject(),
-  fTpar(0),
-  fDebug(0)
+  fWMatrix(0),
+  fTpar(0)
  {
   // Default constructor
   for(Int_t i=0;i<3;i++) {
     fP0[i] = 0.;
+    fSigma2P0[i] = 0.;
     fCd[i] = 0.;
   }
+  SetIdPoints(65535,65535);
 }
 
 //________________________________________________________
-AliStrLine::AliStrLine(Double_t *point, Double_t *cd,Bool_t twopoints) :
+AliStrLine::AliStrLine(const Double_t *const point, const Double_t *const cd, Bool_t twopoints, UShort_t id1, UShort_t id2) :
   TObject(),
-  fTpar(0),
-  fDebug(0)
+  fWMatrix(0),
+  fTpar(0)
 {
   // Standard constructor
   // if twopoints is true:  point and cd are the 3D coordinates of
@@ -49,60 +56,142 @@ AliStrLine::AliStrLine(Double_t *point, Double_t *cd,Bool_t twopoints) :
   // if twopoint is false: point represents the 3D coordinates of a point
   //                       belonging to the straight line and cd is the
   //                       direction in space
-  if(twopoints){
+  for(Int_t i=0;i<3;i++)
+    fSigma2P0[i] = 0.;
+
+  if(twopoints)
     InitTwoPoints(point,cd);
-  }
-  else {
+  else 
     InitDirection(point,cd);
-  }
+
+  SetIdPoints(id1,id2);
 }
 
+
 //________________________________________________________
-AliStrLine::AliStrLine(Float_t *pointf, Float_t *cdf,Bool_t twopoints) :
+AliStrLine::AliStrLine(const Double_t *const point, const Double_t *const sig2point, const Double_t *const cd, Bool_t twopoints, UShort_t id1, UShort_t id2) :
   TObject(),
-  fTpar(0),
-  fDebug(0)
+  fWMatrix(0),
+  fTpar(0)
 {
-  // Standard constructor - with float arguments
+  // Standard constructor
   // if twopoints is true:  point and cd are the 3D coordinates of
   //                        two points defininig the straight line
   // if twopoint is false: point represents the 3D coordinates of a point
   //                       belonging to the straight line and cd is the
   //                       direction in space
-  Double_t point[3];
-  Double_t cd[3];
-  for(Int_t i=0;i<3;i++){
-    point[i] = pointf[i];
-    cd[i] = cdf[i];
-  }
-  if(twopoints){
+  for(Int_t i=0;i<3;i++)
+    fSigma2P0[i] = sig2point[i];
+
+  if(twopoints)
     InitTwoPoints(point,cd);
+  else 
+    InitDirection(point,cd);
+
+  SetIdPoints(id1,id2);
+}
+
+//________________________________________________________
+AliStrLine::AliStrLine(const Double_t *const point, const Double_t *const sig2point, const Double_t *const wmat, const Double_t *const cd, Bool_t twopoints, UShort_t id1, UShort_t id2) :
+  TObject(),
+  fWMatrix(0),
+  fTpar(0)
+{
+  // Standard constructor
+  // if twopoints is true:  point and cd are the 3D coordinates of
+  //                        two points defininig the straight line
+  // if twopoint is false: point represents the 3D coordinates of a point
+  //                       belonging to the straight line and cd is the
+  //                       direction in space
+  Int_t k = 0;
+  fWMatrix = new Double_t [6];
+  for(Int_t i=0;i<3;i++){ 
+    fSigma2P0[i] = sig2point[i];
+    for(Int_t j=0;j<3;j++)if(j>=i)fWMatrix[k++]=wmat[3*i+j];
   }
-  else {
+  if(twopoints)
+    InitTwoPoints(point,cd);
+  else 
     InitDirection(point,cd);
+
+  SetIdPoints(id1,id2);
+}
+
+//________________________________________________________
+AliStrLine::AliStrLine(const AliStrLine &source):TObject(source),
+fWMatrix(0),
+fTpar(source.fTpar){
+  // copy constructor
+  for(Int_t i=0;i<3;i++){
+    fP0[i]=source.fP0[i];
+    fSigma2P0[i]=source.fSigma2P0[i];
+    fCd[i]=source.fCd[i];
+  }
+  if(source.fWMatrix){
+    fWMatrix = new Double_t [6];
+    for(Int_t i=0;i<6;i++)fWMatrix[i]=source.fWMatrix[i];
   }
+  for(Int_t i=0;i<2;i++) fIdPoint[i]=source.fIdPoint[i];
 }
 
 //________________________________________________________
-void AliStrLine::InitDirection(Double_t *point, Double_t *cd){
+AliStrLine& AliStrLine::operator=(const AliStrLine& source){
+  // Assignment operator
+  if(this !=&source){
+    this->~AliStrLine();
+    new(this)AliStrLine(source);
+  }
+  return *this;
+}
+
+//________________________________________________________
+void AliStrLine::GetWMatrix(Double_t *wmat)const {
+// Getter for weighting matrix, as a [9] dim. array
+  if(!fWMatrix)return;
+  Int_t k = 0;
+  for(Int_t i=0;i<3;i++){
+    for(Int_t j=0;j<3;j++){
+      if(j>=i){
+	wmat[3*i+j]=fWMatrix[k++];
+      }
+      else{
+	wmat[3*i+j]=wmat[3*j+i];
+      }
+    }
+  }
+} 
+
+//________________________________________________________
+void AliStrLine::SetWMatrix(const Double_t *wmat) {
+// Setter for weighting matrix, strating from a [9] dim. array
+  if(fWMatrix)delete [] fWMatrix;
+  fWMatrix = new Double_t [6];
+  Int_t k = 0;
+  for(Int_t i=0;i<3;i++){
+    for(Int_t j=0;j<3;j++)if(j>=i)fWMatrix[k++]=wmat[3*i+j];
+  }
+}
+
+//________________________________________________________
+void AliStrLine::InitDirection(const Double_t *const point, const Double_t *const cd)
+{
   // Initialization from a point and a direction
-  Double_t norm = 0.;
-  for(Int_t i=0;i<3;i++)norm+=cd[i]*cd[i];
+  Double_t norm = cd[0]*cd[0]+cd[1]*cd[1]+cd[2]*cd[2];
+
   if(norm) {
-    norm = TMath::Sqrt(norm);
-    for(Int_t i=0;i<3;i++) cd[i]/=norm;
-  }
-  else {
-    Error("AliStrLine","Null direction cosines!!!");
+    norm = TMath::Sqrt(1./norm);
+    for(Int_t i=0;i<3;++i) {
+      fP0[i]=point[i];
+      fCd[i]=cd[i]*norm;
+    }
+    fTpar = 0.;
   }
-  SetP0(point);
-  SetCd(cd);
-  fTpar = 0.;
-  SetDebug();
+  else AliFatal("Null direction cosines!!!");
 }
 
 //________________________________________________________
-void AliStrLine::InitTwoPoints(Double_t *pA, Double_t *pB){
+void AliStrLine::InitTwoPoints(const Double_t *const pA, const Double_t *const pB)
+{
   // Initialization from the coordinates of two
   // points in the space
   Double_t cd[3];
@@ -113,6 +202,7 @@ void AliStrLine::InitTwoPoints(Double_t *pA, Double_t *pB){
 //________________________________________________________
 AliStrLine::~AliStrLine() {
   // destructor
+  if(fWMatrix)delete [] fWMatrix;
 }
 
 //________________________________________________________
@@ -125,26 +215,39 @@ void AliStrLine::PrintStatus() const {
   cout <<"Known point: ";
   for(Int_t i=0;i<3;i++)cout <<fP0[i]<<"; ";
   cout <<endl;
+  cout <<"Error on known point: ";
+  for(Int_t i=0;i<3;i++)cout <<TMath::Sqrt(fSigma2P0[i])<<"; ";
+  cout <<endl;
   cout <<"Current value for the parameter: "<<fTpar<<endl;
-  cout <<" Debug flag: "<<fDebug<<endl;
 }
 
 //________________________________________________________
-Int_t AliStrLine::IsParallelTo(AliStrLine *line) const {
+Int_t AliStrLine::IsParallelTo(const AliStrLine *line) const {
   // returns 1 if lines are parallel, 0 if not paralel
+  const Double_t prec=1e-14;
   Double_t cd2[3];
   line->GetCd(cd2);
+
   Double_t vecpx=fCd[1]*cd2[2]-fCd[2]*cd2[1];
-  if(vecpx!=0) return 0;
+  Double_t mod=TMath::Abs(fCd[1]*cd2[2])+TMath::Abs(fCd[2]*cd2[1]);
+  if(TMath::Abs(vecpx) > prec*mod) return 0;
+
   Double_t vecpy=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
-  if(vecpy!=0) return 0;
+  mod=TMath::Abs(fCd[0]*cd2[2])+TMath::Abs(fCd[2]*cd2[0]);
+  if(TMath::Abs(vecpy) > prec*mod) return 0;
+
   Double_t vecpz=fCd[0]*cd2[1]-fCd[1]*cd2[0];
-  if(vecpz!=0) return 0;
+  mod=TMath::Abs(fCd[0]*cd2[1])+TMath::Abs(fCd[1]*cd2[0]);
+  if(TMath::Abs(vecpz) > prec) return 0;
+
   return 1;
 }
 //________________________________________________________
-Int_t AliStrLine::Crossrphi(AliStrLine *line){
+Int_t AliStrLine::Crossrphi(const AliStrLine *line)
+{
   // Cross 2 lines in the X-Y plane
+  const Double_t prec=1e-14;
+  const Double_t big=1e20;
   Double_t p2[3];
   Double_t cd2[3];
   line->GetP0(p2);
@@ -156,11 +259,13 @@ Int_t AliStrLine::Crossrphi(AliStrLine *line){
   Double_t e=-cd2[1];
   Double_t f=p2[1]-fP0[1];
   Double_t deno = a*e-b*d;
+  Double_t mod=TMath::Abs(a*e)+TMath::Abs(b*d);
   Int_t retcode = 0;
-  if(deno != 0.) {
+  if(TMath::Abs(deno) > prec*mod) {
     fTpar = (c*e-b*f)/deno;
   }
   else {
+    fTpar = big;
     retcode = -1;
   }
   return retcode;
@@ -170,24 +275,30 @@ Int_t AliStrLine::Crossrphi(AliStrLine *line){
 Int_t AliStrLine::CrossPoints(AliStrLine *line, Double_t *point1, Double_t *point2){
   // Looks for the crossing point estimated starting from the
   // DCA segment
+  const Double_t prec=1e-14;
   Double_t p2[3];
   Double_t cd2[3];
   line->GetP0(p2);
   line->GetCd(cd2);
   Int_t i;
   Double_t k1 = 0;
-  for(i=0;i<3;i++)k1+=(fP0[i]-p2[i])*fCd[i];
   Double_t k2 = 0;
-  for(i=0;i<3;i++)k2+=(fP0[i]-p2[i])*cd2[i];
   Double_t a11 = 0;
-  for(i=0;i<3;i++)a11+=fCd[i]*cd2[i];
+  for(i=0;i<3;i++){
+    k1+=(fP0[i]-p2[i])*fCd[i];
+    k2+=(fP0[i]-p2[i])*cd2[i];
+    a11+=fCd[i]*cd2[i];
+  }
   Double_t a22 = -a11;
   Double_t a21 = 0;
-  for(i=0;i<3;i++)a21+=cd2[i]*cd2[i];
   Double_t a12 = 0;
-  for(i=0;i<3;i++)a12-=fCd[i]*fCd[i];
+  for(i=0;i<3;i++){
+    a21+=cd2[i]*cd2[i];
+    a12-=fCd[i]*fCd[i];
+  }
   Double_t deno = a11*a22-a21*a12;
-  if(deno == 0.) return -1;
+  Double_t mod = TMath::Abs(a11*a22)+TMath::Abs(a21*a12);
+  if(TMath::Abs(deno) < prec*mod) return -1;
   fTpar = (a11*k2-a21*k1) / deno;
   Double_t par2 = (k1*a22-k2*a12) / deno;
   line->SetPar(par2);
@@ -196,7 +307,8 @@ Int_t AliStrLine::CrossPoints(AliStrLine *line, Double_t *point1, Double_t *poin
   return 0;
 }
 //________________________________________________________________
-Int_t AliStrLine::Cross(AliStrLine *line, Double_t *point){
+Int_t AliStrLine::Cross(AliStrLine *line, Double_t *point)
+{
 
   //Finds intersection between lines
   Double_t point1[3];
@@ -211,8 +323,10 @@ Int_t AliStrLine::Cross(AliStrLine *line, Double_t *point){
 }
 
 //___________________________________________________________
-Double_t AliStrLine::GetDCA(AliStrLine *line) const{
+Double_t AliStrLine::GetDCA(const AliStrLine *line) const
+{
   //Returns the distance of closest approach between two lines
+  const Double_t prec=1e-14;
   Double_t p2[3];
   Double_t cd2[3];
   line->GetP0(p2);
@@ -226,7 +340,7 @@ Double_t AliStrLine::GetDCA(AliStrLine *line) const{
       dist2+=(fP0[i]-p2[i])*fCd[i];
       mod+=fCd[i]*fCd[i];
     }
-    if(mod!=0){
+    if(TMath::Abs(mod) > prec){
       dist2/=mod;
       return TMath::Sqrt(dist1q-dist2*dist2);
     }else{return -1;}
@@ -240,11 +354,9 @@ Double_t AliStrLine::GetDCA(AliStrLine *line) const{
        mod+=perp[i]*perp[i];
        dist+=(fP0[i]-p2[i])*perp[i];
      }
-     mod=sqrt(mod);
-     if(mod!=0){
-       dist/=mod;
-       return TMath::Abs(dist);
-     }else{return -1;}
+     if(TMath::Abs(mod) > prec) {
+       return TMath::Abs(dist/TMath::Sqrt(mod));
+     } else return -1;
   }
 }
 //________________________________________________________
@@ -254,8 +366,9 @@ void AliStrLine::GetCurrentPoint(Double_t *point) const {
 }
 
 //________________________________________________________
-Double_t AliStrLine::GetDistFromPoint(Double_t *point) const {
+Double_t AliStrLine::GetDistFromPoint(const Double_t *point) const 
+{
   // computes distance from point 
-  AliStrLine tmp(point,(Double_t *)fCd,kFALSE);
-  return this->GetDCA(&tmp);
+  AliStrLine tmpline(point, fCd, kFALSE);
+  return GetDCA(&tmpline);
 }