STEER/AliStrLine.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////
  19 //                                                               //
  20 // A straight line is coded as a point (3 Double_t) and           //
  21 // 3 direction cosines                                           //
  22 //                                                               //
  23 ///////////////////////////////////////////////////////////////////
  24
  25 #include <Riostream.h>
  26 #include <TTree.h>
  27 #include <TMath.h>
  28
  29 #include "AliStrLine.h"
  30
  31 ClassImp(AliStrLine)
  32
  33 //________________________________________________________
  34 AliStrLine::AliStrLine() :
  35   TObject(),
  36   fTpar(0),
  37   fDebug(0)
  38  {
  39   // Default constructor
  40   for(Int_t i=0;i<3;i++) {
  41     fP0[i] = 0.;
  42     fSigma2P0[i] = 0.;
  43     fCd[i] = 0.;
  44   }
  45   for(Int_t i=0;i<9;i++) fWMatrix[i] = 0.;
  46 }
  47
  48 //________________________________________________________
  49 AliStrLine::AliStrLine(Double_t *point, Double_t *cd,Bool_t twopoints) :
  50   TObject(),
  51   fTpar(0),
  52   fDebug(0)
  53 {
  54   // Standard constructor
  55   // if twopoints is true:  point and cd are the 3D coordinates of
  56   //                        two points defininig the straight line
  57   // if twopoint is false: point represents the 3D coordinates of a point
  58   //                       belonging to the straight line and cd is the
  59   //                       direction in space
  60   for(Int_t i=0;i<3;i++){
  61     fSigma2P0[i] = 0.;
  62   }
  63   for(Int_t i=0;i<9;i++) fWMatrix[i] = 0.;
  64   if(twopoints){
  65     InitTwoPoints(point,cd);
  66   }
  67   else {
  68     InitDirection(point,cd);
  69   }
  70 }
  71
  72 //________________________________________________________
  73 AliStrLine::AliStrLine(Float_t *pointf, Float_t *cdf,Bool_t twopoints) :
  74   TObject(),
  75   fTpar(0),
  76   fDebug(0)
  77 {
  78   // Standard constructor - with float arguments
  79   // if twopoints is true:  point and cd are the 3D coordinates of
  80   //                        two points defininig the straight line
  81   // if twopoint is false: point represents the 3D coordinates of a point
  82   //                       belonging to the straight line and cd is the
  83   //                       direction in space
  84   Double_t point[3];
  85   Double_t cd[3];
  86   for(Int_t i=0;i<3;i++){
  87     point[i] = pointf[i];
  88     cd[i] = cdf[i];
  89     fSigma2P0[i] = 0.;
  90   }
  91   for(Int_t i=0;i<9;i++) fWMatrix[i] = 0.;
  92   if(twopoints){
  93     InitTwoPoints(point,cd);
  94   }
  95   else {
  96     InitDirection(point,cd);
  97   }
  98 }
  99
 100 //________________________________________________________
 101 AliStrLine::AliStrLine(Double_t *point, Double_t *sig2point, Double_t *cd,Bool_t twopoints) :
 102   TObject(),
 103   fTpar(0),
 104   fDebug(0)
 105 {
 106   // Standard constructor
 107   // if twopoints is true:  point and cd are the 3D coordinates of
 108   //                        two points defininig the straight line
 109   // if twopoint is false: point represents the 3D coordinates of a point
 110   //                       belonging to the straight line and cd is the
 111   //                       direction in space
 112   for(Int_t i=0;i<3;i++){
 113     fSigma2P0[i] = sig2point[i];
 114   }
 115   for(Int_t i=0;i<9;i++) fWMatrix[i] = 0.;
 116   if(twopoints){
 117     InitTwoPoints(point,cd);
 118   }
 119   else {
 120     InitDirection(point,cd);
 121   }
 122 }
 123
 124 //________________________________________________________
 125 AliStrLine::AliStrLine(Float_t *pointf, Float_t *sig2point, Float_t *cdf,Bool_t twopoints) :
 126   TObject(),
 127   fTpar(0),
 128   fDebug(0)
 129 {
 130   // Standard constructor - with float arguments
 131   // if twopoints is true:  point and cd are the 3D coordinates of
 132   //                        two points defininig the straight line
 133   // if twopoint is false: point represents the 3D coordinates of a point
 134   //                       belonging to the straight line and cd is the
 135   //                       direction in space
 136   Double_t point[3];
 137   Double_t cd[3];
 138   for(Int_t i=0;i<3;i++){
 139     point[i] = pointf[i];
 140     cd[i] = cdf[i];
 141     fSigma2P0[i] = sig2point[i];
 142   }
 143   for(Int_t i=0;i<9;i++) fWMatrix[i] = 0.;
 144   if(twopoints){
 145     InitTwoPoints(point,cd);
 146   }
 147   else {
 148     InitDirection(point,cd);
 149   }
 150 }
 151 //________________________________________________________
 152 AliStrLine::AliStrLine(Double_t *point, Double_t *sig2point, Double_t *wmat, Double_t *cd,Bool_t twopoints) :
 153   TObject(),
 154   fTpar(0),
 155   fDebug(0)
 156 {
 157   // Standard constructor
 158   // if twopoints is true:  point and cd are the 3D coordinates of
 159   //                        two points defininig the straight line
 160   // if twopoint is false: point represents the 3D coordinates of a point
 161   //                       belonging to the straight line and cd is the
 162   //                       direction in space
 163   for(Int_t i=0;i<3;i++){
 164     fSigma2P0[i] = sig2point[i];
 165   }
 166   for(Int_t i=0;i<9;i++) fWMatrix[i] = wmat[i];
 167   if(twopoints){
 168     InitTwoPoints(point,cd);
 169   }
 170   else {
 171     InitDirection(point,cd);
 172   }
 173 }
 174
 175 //________________________________________________________
 176 AliStrLine::AliStrLine(Float_t *pointf, Float_t *sig2point, Float_t *wmat, Float_t *cdf,Bool_t twopoints) :
 177   TObject(),
 178   fTpar(0),
 179   fDebug(0)
 180 {
 181   // Standard constructor - with float arguments
 182   // if twopoints is true:  point and cd are the 3D coordinates of
 183   //                        two points defininig the straight line
 184   // if twopoint is false: point represents the 3D coordinates of a point
 185   //                       belonging to the straight line and cd is the
 186   //                       direction in space
 187   Double_t point[3];
 188   Double_t cd[3];
 189   for(Int_t i=0;i<3;i++){
 190     point[i] = pointf[i];
 191     cd[i] = cdf[i];
 192     fSigma2P0[i] = sig2point[i];
 193   }
 194   for(Int_t i=0;i<9;i++) fWMatrix[i] = wmat[i];
 195   if(twopoints){
 196     InitTwoPoints(point,cd);
 197   }
 198   else {
 199     InitDirection(point,cd);
 200   }
 201 }
 202 //________________________________________________________
 203 void AliStrLine::InitDirection(Double_t *point, Double_t *cd){
 204   // Initialization from a point and a direction
 205   Double_t norm = 0.;
 206   for(Int_t i=0;i<3;i++)norm+=cd[i]*cd[i];
 207   if(norm) {
 208     norm = TMath::Sqrt(norm);
 209     for(Int_t i=0;i<3;i++) cd[i]/=norm;
 210   }
 211   else {
 212     Error("AliStrLine","Null direction cosines!!!");
 213   }
 214   SetP0(point);
 215   SetCd(cd);
 216   fTpar = 0.;
 217   SetDebug();
 218 }
 219
 220 //________________________________________________________
 221 void AliStrLine::InitTwoPoints(Double_t *pA, Double_t *pB){
 222   // Initialization from the coordinates of two
 223   // points in the space
 224   Double_t cd[3];
 225   for(Int_t i=0;i<3;i++)cd[i] = pB[i]-pA[i];
 226   InitDirection(pA,cd);
 227 }
 228
 229 //________________________________________________________
 230 AliStrLine::~AliStrLine() {
 231   // destructor
 232 }
 233
 234 //________________________________________________________
 235 void AliStrLine::PrintStatus() const {
 236   // Print current status
 237   cout <<"=======================================================\n";
 238   cout <<"Direction cosines: ";
 239   for(Int_t i=0;i<3;i++)cout <<fCd[i]<<"; ";
 240   cout <<endl;
 241   cout <<"Known point: ";
 242   for(Int_t i=0;i<3;i++)cout <<fP0[i]<<"; ";
 243   cout <<endl;
 244   cout <<"Error on known point: ";
 245   for(Int_t i=0;i<3;i++)cout <<TMath::Sqrt(fSigma2P0[i])<<"; ";
 246   cout <<endl;
 247   cout <<"Current value for the parameter: "<<fTpar<<endl;
 248   cout <<" Debug flag: "<<fDebug<<endl;
 249 }
 250
 251 //________________________________________________________
 252 Int_t AliStrLine::IsParallelTo(AliStrLine *line) const {
 253   // returns 1 if lines are parallel, 0 if not paralel
 254   Double_t cd2[3];
 255   line->GetCd(cd2);
 256   Double_t vecpx=fCd[1]*cd2[2]-fCd[2]*cd2[1];
 257   if(vecpx!=0) return 0;
 258   Double_t vecpy=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
 259   if(vecpy!=0) return 0;
 260   Double_t vecpz=fCd[0]*cd2[1]-fCd[1]*cd2[0];
 261   if(vecpz!=0) return 0;
 262   return 1;
 263 }
 264 //________________________________________________________
 265 Int_t AliStrLine::Crossrphi(AliStrLine *line){
 266   // Cross 2 lines in the X-Y plane
 267   Double_t p2[3];
 268   Double_t cd2[3];
 269   line->GetP0(p2);
 270   line->GetCd(cd2);
 271   Double_t a=fCd[0];
 272   Double_t b=-cd2[0];
 273   Double_t c=p2[0]-fP0[0];
 274   Double_t d=fCd[1];
 275   Double_t e=-cd2[1];
 276   Double_t f=p2[1]-fP0[1];
 277   Double_t deno = a*e-b*d;
 278   Int_t retcode = 0;
 279   if(deno != 0.) {
 280     fTpar = (c*e-b*f)/deno;
 281   }
 282   else {
 283     retcode = -1;
 284   }
 285   return retcode;
 286 }
 287
 288 //________________________________________________________
 289 Int_t AliStrLine::CrossPoints(AliStrLine *line, Double_t *point1, Double_t *point2){
 290   // Looks for the crossing point estimated starting from the
 291   // DCA segment
 292   Double_t p2[3];
 293   Double_t cd2[3];
 294   line->GetP0(p2);
 295   line->GetCd(cd2);
 296   Int_t i;
 297   Double_t k1 = 0;
 298   for(i=0;i<3;i++)k1+=(fP0[i]-p2[i])*fCd[i];
 299   Double_t k2 = 0;
 300   for(i=0;i<3;i++)k2+=(fP0[i]-p2[i])*cd2[i];
 301   Double_t a11 = 0;
 302   for(i=0;i<3;i++)a11+=fCd[i]*cd2[i];
 303   Double_t a22 = -a11;
 304   Double_t a21 = 0;
 305   for(i=0;i<3;i++)a21+=cd2[i]*cd2[i];
 306   Double_t a12 = 0;
 307   for(i=0;i<3;i++)a12-=fCd[i]*fCd[i];
 308   Double_t deno = a11*a22-a21*a12;
 309   if(deno == 0.) return -1;
 310   fTpar = (a11*k2-a21*k1) / deno;
 311   Double_t par2 = (k1*a22-k2*a12) / deno;
 312   line->SetPar(par2);
 313   GetCurrentPoint(point1);
 314   line->GetCurrentPoint(point2);
 315   return 0;
 316 }
 317 //________________________________________________________________
 318 Int_t AliStrLine::Cross(AliStrLine *line, Double_t *point){
 319
 320   //Finds intersection between lines
 321   Double_t point1[3];
 322   Double_t point2[3];
 323   Int_t retcod=CrossPoints(line,point1,point2);
 324   if(retcod==0){
 325     for(Int_t i=0;i<3;i++)point[i]=(point1[i]+point2[i])/2.;
 326     return 0;
 327   }else{
 328     return -1;
 329   }
 330 }
 331
 332 //___________________________________________________________
 333 Double_t AliStrLine::GetDCA(AliStrLine *line) const{
 334   //Returns the distance of closest approach between two lines
 335   Double_t p2[3];
 336   Double_t cd2[3];
 337   line->GetP0(p2);
 338   line->GetCd(cd2);
 339   Int_t i;
 340   Int_t ispar=IsParallelTo(line);
 341   if(ispar){
 342     Double_t dist1q=0,dist2=0,mod=0;
 343     for(i=0;i<3;i++){
 344       dist1q+=(fP0[i]-p2[i])*(fP0[i]-p2[i]);
 345       dist2+=(fP0[i]-p2[i])*fCd[i];
 346       mod+=fCd[i]*fCd[i];
 347     }
 348     if(mod!=0){
 349       dist2/=mod;
 350       return TMath::Sqrt(dist1q-dist2*dist2);
 351     }else{return -1;}
 352   }else{
 353      Double_t perp[3];
 354      perp[0]=fCd[1]*cd2[2]-fCd[2]*cd2[1];
 355      perp[1]=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
 356      perp[2]=fCd[0]*cd2[1]-fCd[1]*cd2[0];
 357      Double_t mod=0,dist=0;
 358      for(i=0;i<3;i++){
 359        mod+=perp[i]*perp[i];
 360        dist+=(fP0[i]-p2[i])*perp[i];
 361      }
 362      mod=sqrt(mod);
 363      if(mod!=0){
 364        dist/=mod;
 365        return TMath::Abs(dist);
 366      }else{return -1;}
 367   }
 368 }
 369 //________________________________________________________
 370 void AliStrLine::GetCurrentPoint(Double_t *point) const {
 371   // Fills the array point with the current value on the line
 372   for(Int_t i=0;i<3;i++)point[i]=fP0[i]+fCd[i]*fTpar;
 373 }
 374
 375 //________________________________________________________
 376 Double_t AliStrLine::GetDistFromPoint(Double_t *point) const {
 377   // computes distance from point
 378   AliStrLine tmp(point,(Double_t *)fCd,kFALSE);
 379   return this->GetDCA(&tmp);
 380 }