[u/mrichter/AliRoot.git] / STEER / AliTrackFitterRieman.cxx

/**************************************************************************
 * 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.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////////////////
//
// Class to the track points on the Riemann sphere. Inputs are
// the set of id's (volids) of the volumes in which residuals are
// calculated to construct a chi2 function to be minimized during 
// the alignment procedures. For the moment the track extrapolation is 
// taken 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.
//
// Internal usage of AliRieman class for minimization
//
//////////////////////////////////////////////////////////////////////////////

#include "TMatrixDSym.h"
#include "TMatrixD.h"
#include "TArrayI.h"
#include "AliTrackFitterRieman.h"
#include "AliLog.h"
#include "TTreeStream.h"
#include "AliRieman.h"
#include "AliLog.h"

ClassImp(AliTrackFitterRieman)


AliTrackFitterRieman::AliTrackFitterRieman():
  AliTrackFitter()
{
  //
  // default constructor
  //
  fAlpha = 0.;
  fNUsed = 0;
  fConv = kFALSE;
  fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
  fRieman = new AliRieman(10000);  // allocate rieman
}


AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
  AliTrackFitter(array,owner)
{
  //
  // Constructor
  //
  fAlpha = 0.;
  fNUsed = 0;
  fConv = kFALSE;
  if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
  fRieman = new AliRieman(10000);  //allocate rieman
}

AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
  AliTrackFitter(rieman)
{
  //
  // copy constructor
  //
  fAlpha = rieman.fAlpha;
  fNUsed = rieman.fNUsed;
  fConv = rieman.fConv;
  fRieman = new AliRieman(*(rieman.fRieman));
}

//_____________________________________________________________________________
AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
{
  //
  // Assignment operator
  //
  if(this==&rieman) return *this;
  ((AliTrackFitter *)this)->operator=(rieman);

  fAlpha  = rieman.fAlpha;
  fNUsed  = rieman.fNUsed;
  fConv   = rieman.fConv;
  fRieman = new AliRieman(*(rieman.fRieman));
  return *this;
}


AliTrackFitterRieman::~AliTrackFitterRieman(){
  //
  //
  //
  delete fRieman;
  delete fDebugStream;
}

void AliTrackFitterRieman::Reset()
{
  // Reset the track parameters and
  // rieman sums
  //
  AliTrackFitter::Reset();
  fRieman->Reset();
  fAlpha = 0.;
  fNUsed = 0;
  fConv =kFALSE;
}

Bool_t AliTrackFitterRieman::Fit(const TArrayI *volIds,const TArrayI *volIdsFit,
				 AliAlignObj::ELayerID layerRangeMin,
			 AliAlignObj::ELayerID layerRangeMax)
{
  // Fit the track points. The method takes as an input
  // the set of id's (volids) of the volumes in which
  // one wants to calculate the residuals.
  // The following parameters are used to define the
  // range of volumes to be used in the fitting
  // As a result two AliTrackPointArray's obects are filled.
  // The first one contains the space points with
  // volume id's from volids list. The second array of points represents
  // the track extrapolations corresponding to the space points
  // in the first array. The two arrays can be used to find
  // the residuals in the volids and consequently construct a
  // chi2 function to be minimized during the alignment
  // procedures. For the moment the track extrapolation is taken
  // 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.

  Reset();
  Int_t npoints = fPoints->GetNPoints();
  if (fPoints && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1){
    Int_t nVol    = volIds->GetSize();
    Int_t nVolFit = volIdsFit->GetSize();
    Int_t volId  = volIds->At(0);    
    (*fDebugStream)<<"PInput"<<
      "NPoints="<<npoints<<    // number of points
      "VolId="<<volId<<        // first vol ID
      "NVol="<<nVol<<          // number of volumes 
      "NvolFit="<<nVolFit<<    // number of volumes to fit
      "fPoints.="<<fPoints<<   // input points
      "\n";
  }
  if (npoints < 3) return kFALSE;

  Bool_t isAlphaCalc = kFALSE;
  AliTrackPoint p,plocal;
//   fPoints->GetPoint(p,0);
//   fAlpha = TMath::ATan2(p.GetY(),p.GetX());

  Int_t npVolId = 0;
  fNUsed = 0;
  Int_t *pindex = new Int_t[npoints];
  for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
    {
      fPoints->GetPoint(p,ipoint);
      UShort_t iVolId = p.GetVolumeID();
      if (FindVolId(volIds,iVolId)) {
	pindex[npVolId] = ipoint;
	npVolId++;
      }
      if (volIdsFit != 0x0) {
	if (!FindVolId(volIdsFit,iVolId)) continue;
      }
      else {
	if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
	    iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,
						AliAlignObj::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){
	printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
	p.Dump();
	plocal.Dump();
	printf("Problematic point\n");	
	printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
      }
      AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
	       TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
      // fNUsed++;  AddPoint should be responsible
    }

  if (npVolId == 0 || fNUsed < fMinNPoints) {
    delete [] pindex;
    return kFALSE;
  }

  Update();
 
  if (!fConv) {
    delete [] pindex;
    return kFALSE;
  }

  if ((fParams[0] == 0) ||
      ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
    delete [] pindex;
    return kFALSE;
  }


  if (fNUsed < fMinNPoints) {
    delete [] pindex;
    return kFALSE;
  }

  fPVolId = new AliTrackPointArray(npVolId);
  fPTrack = new AliTrackPointArray(npVolId);
  AliTrackPoint p2;
  for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
    {
      Int_t index = pindex[ipoint];
      fPoints->GetPoint(p,index);
      if (GetPCA(p,p2)) {
	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]);
	fPVolId->AddPoint(ipoint,&p);
	fPTrack->AddPoint(ipoint,&p2);
      }else{
	// what should be default bahavior -  
	delete [] pindex;
	delete fPVolId;
	delete fPTrack;
	fPVolId =0;
	fPTrack =0;
	return kFALSE;
      }
    }  
  
  if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0){
    //
    // Debug Info
    //
    AliTrackPointArray *lPVolId = new AliTrackPointArray(npVolId);
    AliTrackPointArray *lPTrack = 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);
      }
    }      
    //
    // debug info
    //
    Int_t nVol    = volIds->GetSize();
    Int_t nVolFit = volIdsFit->GetSize();
    Int_t volId   = volIds->At(0);   
    Int_t modId   =0;
    Int_t layer   =  AliAlignObj::VolUIDToLayer(volId,modId);
    
    (*fDebugStream)<<"Fit"<<
      "VolId="<<volId<<        // volume ID
      "Layer="<<layer<<        // layer ID
      "Module="<<modId<<       // module ID
      "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      
      "Rieman.="<<this<<        // original rieman fit
      "Res.="<<residual<<       // residuals of rieman fit  
      "\n";
    delete lPVolId;
    delete lPTrack;
    delete residual;
  }

  delete [] pindex;
  return kTRUE;
}


void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
{
  //
  // add point to rieman fitter
  //
  fRieman->AddPoint(x,y,z,sy,sz);
  fNUsed = fRieman->GetN();
}


void AliTrackFitterRieman::Update(){
  //
  // 
  //
  fRieman->Update();
  fConv = kFALSE;
  if (fRieman->IsValid()){
    for (Int_t ipar=0; ipar<6; ipar++){
      fParams[ipar] = fRieman->GetParam()[ipar];
    }
    fChi2  =  fRieman->GetChi2();
    fNdf   =  fRieman->GetN()- 2;
    fNUsed =  fRieman->GetN();
    fConv = kTRUE;
  }
}

//_____________________________________________________________________________
Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
{
  //
  // Get the closest to a given spacepoint track trajectory point
  // Look for details in the description of the Fit() method
  //
  if (!fConv) return kFALSE;

  // First X and Y coordinates
  Double_t sin = TMath::Sin(fAlpha);
  Double_t cos = TMath::Cos(fAlpha);
  //  fParam[0]  = 1/y0
  //  fParam[1]  = -x0/y0
  //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
  if (fParams[0] == 0) return kFALSE;
  // Track parameters in the global coordinate system
  Double_t x0 = -fParams[1]/fParams[0]*cos -         1./fParams[0]*sin;
  Double_t y0 =          1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
  if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
  Double_t r  = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
                fParams[0];

  // Define space-point refence plane
  Double_t alphap = p.GetAngle();
  Double_t sinp = TMath::Sin(alphap);
  Double_t cosp = TMath::Cos(alphap);
  Double_t x  = p.GetX()*cosp + p.GetY()*sinp;
  Double_t y  = p.GetY()*cosp - p.GetX()*sinp;
  Double_t x0p= x0*cosp + y0*sinp;
  Double_t y0p= y0*cosp - x0*sinp;
  if ((r*r - (x-x0p)*(x-x0p))<0) {
    AliWarning(Form("Track extrapolation failed ! (Track radius = %f, track circle x = %f, space-point x = %f, reference plane angle = %f\n",r,x0p,x,alphap));
    return kFALSE;
  }
  Double_t temp = TMath::Sqrt(r*r - (x-x0p)*(x-x0p));
  Double_t y1 = y0p + temp;
  Double_t y2 = y0p - temp;
  Double_t yprime = y1;
  if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2;

  // Back to the global coordinate system
  Double_t xsecond = x*cosp - yprime*sinp;
  Double_t ysecond = yprime*cosp + x*sinp;

  // Now Z coordinate and track angles
  Double_t x2 = xsecond*cos + ysecond*sin;
  Double_t zsecond = GetZat(x2);
  Double_t dydx = GetDYat(x2);
  Double_t dzdx = GetDZat(x2);

  // Fill the cov matrix of the track extrapolation point
  Double_t cov[6] = {0,0,0,0,0,0};
  Double_t sigmax = 100*100.;
  cov[0] = sigmax;           cov[1] = sigmax*dydx;      cov[2] = sigmax*dzdx;
  cov[3] = sigmax*dydx*dydx; cov[4] = sigmax*dydx*dzdx;
  cov[5] = sigmax*dzdx*dzdx;

  Float_t  newcov[6];
  newcov[0] = cov[0]*cos*cos-
            2*cov[1]*sin*cos+
              cov[3]*sin*sin;
  newcov[1] = cov[1]*(cos*cos-sin*sin)-
             (cov[3]-cov[0])*sin*cos;
  newcov[2] = cov[2]*cos-
              cov[4]*sin;
  newcov[3] = cov[0]*sin*sin+
            2*cov[1]*sin*cos+
              cov[3]*cos*cos;
  newcov[4] = cov[4]*cos+
              cov[2]*sin;
  newcov[5] = cov[5];

  p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
  if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>1){
    AliTrackPoint lp0(p);
    AliTrackPoint lp2(p2);
    if (0)(*fDebugStream)<<"PCA"<<
      "P0.="<<&lp0<<
      "P2.="<<&lp2<<
      "\n";
  }
  return kTRUE;
}
Commit	Line	Data
6b6cba33	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	// Class to the track points on the Riemann sphere. Inputs are
	21	// the set of id's (volids) of the volumes in which residuals are
	22	// calculated to construct a chi2 function to be minimized during
	23	// the alignment procedures. For the moment the track extrapolation is
	24	// taken at the space-point reference plane. The reference plane is
	25	// found using the covariance matrix of the point
	26	// (assuming sigma(x)=0 at the reference coordinate system.
	27	//
8849da04	28	// Internal usage of AliRieman class for minimization
8849da04	29	//
6b6cba33	30	//////////////////////////////////////////////////////////////////////////////
6b6cba33	31
98937d93	32	#include "TMatrixDSym.h"
98937d93	33	#include "TMatrixD.h"
8849da04	34	#include "TArrayI.h"
98937d93	35	#include "AliTrackFitterRieman.h"
24d4520d	36	#include "AliLog.h"
8849da04	37	#include "TTreeStream.h"
	38	#include "AliRieman.h"
	39	#include "AliLog.h"
98937d93	40
	41	ClassImp(AliTrackFitterRieman)
	42
8849da04	43
98937d93	44	AliTrackFitterRieman::AliTrackFitterRieman():
	45	AliTrackFitter()
	46	{
	47	//
	48	// default constructor
	49	//
	50	fAlpha = 0.;
46ae650f	51	fNUsed = 0;
98937d93	52	fConv = kFALSE;
8849da04	53	fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
8849da04	54	fRieman = new AliRieman(10000); // allocate rieman
98937d93	55	}
	56
	57
	58	AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
	59	AliTrackFitter(array,owner)
	60	{
	61	//
	62	// Constructor
	63	//
	64	fAlpha = 0.;
46ae650f	65	fNUsed = 0;
98937d93	66	fConv = kFALSE;
8849da04	67	if (AliLog::GetDebugLevel("","AliTrackFitterRieman")) fDebugStream = new TTreeSRedirector("RiemanAlignDebug.root");
8849da04	68	fRieman = new AliRieman(10000); //allocate rieman
98937d93	69	}
	70
	71	AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
	72	AliTrackFitter(rieman)
	73	{
	74	//
	75	// copy constructor
	76	//
	77	fAlpha = rieman.fAlpha;
46ae650f	78	fNUsed = rieman.fNUsed;
98937d93	79	fConv = rieman.fConv;
8849da04	80	fRieman = new AliRieman(*(rieman.fRieman));
98937d93	81	}
	82
	83	//_____________________________________________________________________________
	84	AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
	85	{
6b6cba33	86	//
6b6cba33	87	// Assignment operator
98937d93	88	//
	89	if(this==&rieman) return *this;
	90	((AliTrackFitter *)this)->operator=(rieman);
	91
8849da04	92	fAlpha = rieman.fAlpha;
	93	fNUsed = rieman.fNUsed;
	94	fConv = rieman.fConv;
	95	fRieman = new AliRieman(*(rieman.fRieman));
98937d93	96	return *this;
	97	}
	98
8849da04	99
	100	AliTrackFitterRieman::~AliTrackFitterRieman(){
	101	//
	102	//
	103	//
	104	delete fRieman;
	105	delete fDebugStream;
	106	}
	107
98937d93	108	void AliTrackFitterRieman::Reset()
	109	{
	110	// Reset the track parameters and
	111	// rieman sums
6b6cba33	112	//
98937d93	113	AliTrackFitter::Reset();
8849da04	114	fRieman->Reset();
98937d93	115	fAlpha = 0.;
46ae650f	116	fNUsed = 0;
98937d93	117	fConv =kFALSE;
	118	}
	119
cc345ce3	120	Bool_t AliTrackFitterRieman::Fit(const TArrayI volIds,const TArrayI volIdsFit,
98937d93	121	AliAlignObj::ELayerID layerRangeMin,
8849da04	122	AliAlignObj::ELayerID layerRangeMax)
98937d93	123	{
98937d93	124	// Fit the track points. The method takes as an input
cc345ce3	125	// the set of id's (volids) of the volumes in which
	126	// one wants to calculate the residuals.
	127	// The following parameters are used to define the
98937d93	128	// range of volumes to be used in the fitting
	129	// As a result two AliTrackPointArray's obects are filled.
	130	// The first one contains the space points with
cc345ce3	131	// volume id's from volids list. The second array of points represents
cc345ce3	132	// the track extrapolations corresponding to the space points
98937d93	133	// in the first array. The two arrays can be used to find
cc345ce3	134	// the residuals in the volids and consequently construct a
98937d93	135	// chi2 function to be minimized during the alignment
98937d93	136	// procedures. For the moment the track extrapolation is taken
cc345ce3	137	// at the space-point reference plane. The reference plane is
	138	// found using the covariance matrix of the point
	139	// (assuming sigma(x)=0 at the reference coordinate system.
98937d93	140
46ae650f	141	Reset();
98937d93	142	Int_t npoints = fPoints->GetNPoints();
8849da04	143	if (fPoints && AliLog::GetDebugLevel("","AliTrackFitterRieman")>1){
	144	Int_t nVol = volIds->GetSize();
	145	Int_t nVolFit = volIdsFit->GetSize();
	146	Int_t volId = volIds->At(0);
	147	(*fDebugStream)<<"PInput"<<
	148	"NPoints="<<npoints<< // number of points
	149	"VolId="<<volId<< // first vol ID
	150	"NVol="<<nVol<< // number of volumes
	151	"NvolFit="<<nVolFit<< // number of volumes to fit
	152	"fPoints.="<<fPoints<< // input points
	153	"\n";
	154	}
98937d93	155	if (npoints < 3) return kFALSE;
98937d93	156
cc345ce3	157	Bool_t isAlphaCalc = kFALSE;
46ae650f	158	AliTrackPoint p,plocal;
cc345ce3	159	// fPoints->GetPoint(p,0);
cc345ce3	160	// fAlpha = TMath::ATan2(p.GetY(),p.GetX());
98937d93	161
98937d93	162	Int_t npVolId = 0;
46ae650f	163	fNUsed = 0;
98937d93	164	Int_t *pindex = new Int_t[npoints];
	165	for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
	166	{
	167	fPoints->GetPoint(p,ipoint);
	168	UShort_t iVolId = p.GetVolumeID();
cc345ce3	169	if (FindVolId(volIds,iVolId)) {
98937d93	170	pindex[npVolId] = ipoint;
	171	npVolId++;
	172	}
cc345ce3	173	if (volIdsFit != 0x0) {
cc345ce3	174	if (!FindVolId(volIdsFit,iVolId)) continue;
46ae650f	175	}
	176	else {
	177	if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) \|\|
	178	iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,
c041444f	179	AliAlignObj::LayerSize(layerRangeMax))) continue;
46ae650f	180	}
cc345ce3	181	if (!isAlphaCalc) {
	182	fAlpha = p.GetAngle();
	183	isAlphaCalc = kTRUE;
	184	}
46ae650f	185	plocal = p.Rotate(fAlpha);
8849da04	186	if (TMath::Abs(plocal.GetX())>500 \|\| TMath::Abs(plocal.GetX())<2){
	187	printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
	188	p.Dump();
	189	plocal.Dump();
	190	printf("Problematic point\n");
	191	printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<</n");
	192	}
46ae650f	193	AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
46ae650f	194	TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
8849da04	195	// fNUsed++; AddPoint should be responsible
98937d93	196	}
98937d93	197
8849da04	198	if (npVolId == 0 \|\| fNUsed < fMinNPoints) {
98937d93	199	delete [] pindex;
	200	return kFALSE;
	201	}
	202
	203	Update();
cc345ce3	204
98937d93	205	if (!fConv) {
	206	delete [] pindex;
	207	return kFALSE;
	208	}
	209
	210	if ((fParams[0] == 0) \|\|
	211	((-fParams[2]fParams[0]+fParams[1]fParams[1]+1) <= 0)) {
	212	delete [] pindex;
	213	return kFALSE;
	214	}
	215
cc345ce3	216
	217	if (fNUsed < fMinNPoints) {
	218	delete [] pindex;
	219	return kFALSE;
	220	}
	221
46ae650f	222	fPVolId = new AliTrackPointArray(npVolId);
46ae650f	223	fPTrack = new AliTrackPointArray(npVolId);
98937d93	224	AliTrackPoint p2;
	225	for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
	226	{
	227	Int_t index = pindex[ipoint];
	228	fPoints->GetPoint(p,index);
	229	if (GetPCA(p,p2)) {
cc345ce3	230	Float_t xyz[3],xyz2[3];
	231	p.GetXYZ(xyz); p2.GetXYZ(xyz2);
	232	// 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]);
46ae650f	233	fPVolId->AddPoint(ipoint,&p);
46ae650f	234	fPTrack->AddPoint(ipoint,&p2);
8849da04	235	}else{
	236	// what should be default bahavior -
	237	delete [] pindex;
	238	delete fPVolId;
	239	delete fPTrack;
	240	fPVolId =0;
	241	fPTrack =0;
	242	return kFALSE;
98937d93	243	}
98937d93	244	}
8849da04	245
	246	if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>0){
	247	//
	248	// Debug Info
	249	//
	250	AliTrackPointArray *lPVolId = new AliTrackPointArray(npVolId);
	251	AliTrackPointArray *lPTrack = new AliTrackPointArray(npVolId);
	252	AliRieman * residual = fRieman->MakeResiduals();
	253	AliTrackPoint p2local;
	254	for (Int_t ipoint = 0; ipoint < npVolId; ipoint++){
	255	Int_t index = pindex[ipoint];
	256	fPoints->GetPoint(p,index);
	257	if (GetPCA(p,p2)) {
	258	plocal = p.Rotate(fAlpha);
	259	// plocal.Rotate(fAlpha);
	260	Float_t xyz[3],xyz2[3];
	261	p2local = plocal;
	262	plocal.GetXYZ(xyz);
	263	xyz2[0] = xyz[0];
	264	xyz2[1] = GetYat(xyz[0]);
	265	xyz2[2] = GetZat(xyz[0]);
	266	p2local.SetXYZ(xyz2);
	267	lPVolId->AddPoint(ipoint,&plocal);
	268	lPTrack->AddPoint(ipoint,&p2local);
	269	}
	270	}
	271	//
	272	// debug info
	273	//
	274	Int_t nVol = volIds->GetSize();
	275	Int_t nVolFit = volIdsFit->GetSize();
	276	Int_t volId = volIds->At(0);
	277	Int_t modId =0;
	278	Int_t layer = AliAlignObj::VolUIDToLayer(volId,modId);
	279
	280	(*fDebugStream)<<"Fit"<<
	281	"VolId="<<volId<< // volume ID
	282	"Layer="<<layer<< // layer ID
	283	"Module="<<modId<< // module ID
	284	"NVol="<<nVol<< // number of volumes
	285	"NvolFit="<<nVolFit<< // number of volumes to fit
	286	"Points0.="<<fPVolId<< // original points
	287	"Points1.="<<fPTrack<< // fitted points
	288	"LPoints0.="<<lPVolId<< // original points - local frame
	289	"LPoints1.="<<lPTrack<< // fitted points - local frame
	290	"Rieman.="<<this<< // original rieman fit
	291	"Res.="<<residual<< // residuals of rieman fit
	292	"\n";
	293	delete lPVolId;
	294	delete lPTrack;
	295	delete residual;
	296	}
98937d93	297
98937d93	298	delete [] pindex;
98937d93	299	return kTRUE;
	300	}
	301
8849da04	302
98937d93	303	void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
	304	{
	305	//
8849da04	306	// add point to rieman fitter
98937d93	307	//
8849da04	308	fRieman->AddPoint(x,y,z,sy,sz);
8849da04	309	fNUsed = fRieman->GetN();
98937d93	310	}
98937d93	311
8849da04	312
8849da04	313
98937d93	314	void AliTrackFitterRieman::Update(){
98937d93	315	//
8849da04	316	//
98937d93	317	//
8849da04	318	fRieman->Update();
	319	fConv = kFALSE;
	320	if (fRieman->IsValid()){
	321	for (Int_t ipar=0; ipar<6; ipar++){
	322	fParams[ipar] = fRieman->GetParam()[ipar];
98937d93	323	}
8849da04	324	fChi2 = fRieman->GetChi2();
	325	fNdf = fRieman->GetN()- 2;
	326	fNUsed = fRieman->GetN();
	327	fConv = kTRUE;
98937d93	328	}
98937d93	329	}
98937d93	330
6b6cba33	331	//_____________________________________________________________________________
98937d93	332	Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
98937d93	333	{
6b6cba33	334	//
98937d93	335	// Get the closest to a given spacepoint track trajectory point
98937d93	336	// Look for details in the description of the Fit() method
6b6cba33	337	//
98937d93	338	if (!fConv) return kFALSE;
	339
	340	// First X and Y coordinates
	341	Double_t sin = TMath::Sin(fAlpha);
	342	Double_t cos = TMath::Cos(fAlpha);
	343	// fParam[0] = 1/y0
	344	// fParam[1] = -x0/y0
	345	// fParam[2] = - (R^2 - x0^2 - y0^2)/y0
	346	if (fParams[0] == 0) return kFALSE;
cc345ce3	347	// Track parameters in the global coordinate system
98937d93	348	Double_t x0 = -fParams[1]/fParams[0]cos - 1./fParams[0]sin;
	349	Double_t y0 = 1./fParams[0]cos - fParams[1]/fParams[0]sin;
	350	if ((-fParams[2]fParams[0]+fParams[1]fParams[1]+1) <= 0) return kFALSE;
6b6cba33	351	Double_t r = TMath::Sqrt(-fParams[2]fParams[0]+fParams[1]fParams[1]+1)/
98937d93	352	fParams[0];
98937d93	353
cc345ce3	354	// Define space-point refence plane
	355	Double_t alphap = p.GetAngle();
	356	Double_t sinp = TMath::Sin(alphap);
	357	Double_t cosp = TMath::Cos(alphap);
	358	Double_t x = p.GetX()cosp + p.GetY()sinp;
	359	Double_t y = p.GetY()cosp - p.GetX()sinp;
	360	Double_t x0p= x0cosp + y0sinp;
	361	Double_t y0p= y0cosp - x0sinp;
6b6cba33	362	if ((rr - (x-x0p)(x-x0p))<0) {
6b6cba33	363	AliWarning(Form("Track extrapolation failed ! (Track radius = %f, track circle x = %f, space-point x = %f, reference plane angle = %f\n",r,x0p,x,alphap));
24d4520d	364	return kFALSE;
24d4520d	365	}
6b6cba33	366	Double_t temp = TMath::Sqrt(rr - (x-x0p)(x-x0p));
cc345ce3	367	Double_t y1 = y0p + temp;
	368	Double_t y2 = y0p - temp;
	369	Double_t yprime = y1;
	370	if(TMath::Abs(y2-y) < TMath::Abs(y1-y)) yprime = y2;
	371
	372	// Back to the global coordinate system
	373	Double_t xsecond = xcosp - yprimesinp;
	374	Double_t ysecond = yprimecosp + xsinp;
	375
	376	// Now Z coordinate and track angles
	377	Double_t x2 = xsecondcos + ysecondsin;
	378	Double_t zsecond = GetZat(x2);
	379	Double_t dydx = GetDYat(x2);
	380	Double_t dzdx = GetDZat(x2);
	381
	382	// Fill the cov matrix of the track extrapolation point
	383	Double_t cov[6] = {0,0,0,0,0,0};
	384	Double_t sigmax = 100*100.;
	385	cov[0] = sigmax; cov[1] = sigmaxdydx; cov[2] = sigmaxdzdx;
	386	cov[3] = sigmaxdydxdydx; cov[4] = sigmaxdydxdzdx;
	387	cov[5] = sigmaxdzdxdzdx;
	388
	389	Float_t newcov[6];
	390	newcov[0] = cov[0]coscos-
	391	2cov[1]sin*cos+
	392	cov[3]sinsin;
	393	newcov[1] = cov[1](coscos-sin*sin)-
	394	(cov[3]-cov[0])sincos;
	395	newcov[2] = cov[2]*cos-
	396	cov[4]*sin;
	397	newcov[3] = cov[0]sinsin+
	398	2cov[1]sin*cos+
	399	cov[3]coscos;
	400	newcov[4] = cov[4]*cos+
	401	cov[2]*sin;
	402	newcov[5] = cov[5];
	403
	404	p2.SetXYZ(xsecond,ysecond,zsecond,newcov);
8849da04	405	if (AliLog::GetDebugLevel("","AliTrackFitterRieman")>1){
	406	AliTrackPoint lp0(p);
	407	AliTrackPoint lp2(p2);
	408	if (0)(*fDebugStream)<<"PCA"<<
	409	"P0.="<<&lp0<<
	410	"P2.="<<&lp2<<
	411	"\n";
	412	}
98937d93	413	return kTRUE;
98937d93	414	}