[u/mrichter/AliRoot.git] / HLT / src / AliL3ConfMapFit.cxx

// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
//*-- Copyright &copy ASV 

#include <math.h>

#include "AliL3Logging.h"
#include "AliL3ConfMapFit.h"
#include "AliL3Vertex.h"
#include "AliL3ConfMapTrack.h"
#include "AliL3ConfMapPoint.h"

//_____________________________________________________________
// AliL3ConfMapFit
//
// Fit class for conformal mapping tracking

ClassImp(AliL3ConfMapFit)

Double_t AliL3ConfMapFit::pi=3.14159265358979323846;

AliL3ConfMapFit::AliL3ConfMapFit(AliL3ConfMapTrack *track,AliL3Vertex *vertex)
{
  //constructor
  fTrack = track;
  fVertex = vertex;
  BFACT = 0.0029980;
  bField = 0.2;
  
}

Int_t AliL3ConfMapFit::FitHelix()
{
  if(FitCircle())
    {
      LOG(AliL3Log::kError,"AliL3ConfMapFit::FitHelix","TrackFit")<<AliL3Log::kDec<<
	"Problems during circle fit"<<ENDLOG;
      return 1;
    }
  if(FitLine())
    {
      LOG(AliL3Log::kError,"AliL3ConfMapFit::FitHelix","TrackFit")<<AliL3Log::kDec<<
	"Problems during line fit"<<ENDLOG;
      return 1;
    }
  return 0;
}

Int_t AliL3ConfMapFit::FitCircle()
{
  //-----------------------------------------------------------------
  //Fits circle parameters using algorithm
  //described by ChErnov and Oskov in Computer Physics
  //Communications.
  // 
  //Written in FORTRAN by Jawluen Tang, Physics department , UT-Austin 
  //Moved to C by Pablo Yepes
  //Moved to AliROOT by ASV.
  //------------------------------------------------------------------
  
  Double_t wsum  = 0.0 ;
  Double_t xav   = 0.0 ;
  Double_t yav   = 0.0 ;
  
  Int_t num_of_hits = fTrack->GetNumberOfPoints();
  //
  //     Loop over hits calculating average
  Int_t co=0;
  
  for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
    {
      co++;
      AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)fTrack->currentHit;
      cHit->SetXYWeight( 1./ (Double_t)(cHit->GetXerr()*cHit->GetXerr() + cHit->GetYerr()*cHit->GetYerr()) );
      wsum      += cHit->GetXYWeight() ;
      xav       += cHit->GetXYWeight() * cHit->GetX() ;
      yav       += cHit->GetXYWeight() * cHit->GetY() ;
    }
  if(co!=num_of_hits) 
    LOG(AliL3Log::kError,"AliL3ConfMapFit::FitCircle","TrackFit")<<AliL3Log::kDec<<
      "Mismatch of hits. Counter: "<<co<<" nHits: "<<num_of_hits<<ENDLOG;
  if (fTrack->ComesFromMainVertex() == true)
    {    
      wsum += fVertex->GetXYWeight() ;
      xav  += fVertex->GetX() ;
      yav  += fVertex->GetY() ;
    }
  
  xav = xav / wsum ;
  yav = yav / wsum ;
//
//  CALCULATE <X**2>, <XY>, AND <Y**2> WITH <X> = 0, & <Y> = 0
//
  Double_t xxav  = 0.0 ;
  Double_t xyav  = 0.0 ; 
  Double_t yyav  = 0.0 ;
  Double_t xi, yi ;

  for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
    { 
      //AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint *)hits->At(hit_counter);
      AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)fTrack->currentHit;
      xi        = cHit->GetX() - xav ;
      yi        = cHit->GetY() - yav ;
      xxav     += xi * xi * cHit->GetXYWeight() ;
      xyav     += xi * yi * cHit->GetXYWeight() ;
      yyav     += yi * yi * cHit->GetXYWeight() ;
    }
  
  if (fTrack->ComesFromMainVertex() == true)
    {
      xi        = fVertex->GetX() - xav ;
      yi        = fVertex->GetY() - yav ;
      xxav     += xi * xi * fVertex->GetXYWeight() ;
      xyav     += xi * yi * fVertex->GetXYWeight() ;
      yyav     += yi * yi * fVertex->GetXYWeight() ; 
    }
  xxav = xxav / wsum ;
  xyav = xyav / wsum ;
  yyav = yyav / wsum ;
//
//-->  ROTATE COORDINATES SO THAT <XY> = 0
//
//-->  SIGN(C**2 - S**2) = SIGN(XXAV - YYAV) >
//-->  &                                     > ==> NEW : (XXAV-YYAV) > 0
//-->  SIGN(S) = SIGN(XYAV)                  >

  Double_t a = fabs( xxav - yyav ) ;
  Double_t b = 4.0 * xyav * xyav ;

  Double_t asqpb  = a * a + b  ;
  Double_t rasqpb = sqrt ( asqpb) ;

  Double_t splus  = 1.0 + a / rasqpb ;
  Double_t sminus = b / (asqpb * splus) ;

  splus  = sqrt (0.5 * splus ) ;
  sminus = sqrt (0.5 * sminus) ;
//
//->  FIRST REQUIRE : SIGN(C**2 - S**2) = SIGN(XXAV - YYAV)
//
  Double_t sinrot, cosrot ;
  if ( xxav <= yyav ) {
	 cosrot = sminus ;
	 sinrot = splus  ;
  }
  else {
	  cosrot = splus ;
	  sinrot = sminus ;
  }
//
//->  REQUIRE : SIGN(S) = SIGN(XYAV) * SIGN(C) (ASSUMING SIGN(C) > 0)
//
  if ( xyav < 0.0 ) sinrot = - sinrot ;
//
//-->  WE NOW HAVE THE SMALLEST ANGLE THAT GUARANTEES <X**2> > <Y**2>
//-->  TO GET THE SIGN OF THE CHARGE RIGHT, THE NEW X-AXIS MUST POINT
//-->  OUTWARD FROM THE ORGIN.  WE ARE FREE TO CHANGE SIGNS OF BOTH
//-->  COSROT AND SINROT SIMULTANEOUSLY TO ACCOMPLISH THIS.
//
//-->  CHOOSE SIGN OF C WISELY TO BE ABLE TO GET THE SIGN OF THE CHARGE
//
  if ( cosrot*xav+sinrot*yav < 0.0 ) {
	  cosrot = -cosrot ;
	  sinrot = -sinrot ;
  }
//
//->  NOW GET <R**2> AND RSCALE= SQRT(<R**2>)
//
  Double_t rrav   = xxav + yyav ;
  Double_t rscale = sqrt(rrav) ;

  xxav   = 0.0 ;
  yyav   = 0.0 ;
  xyav   = 0.0 ;
  Double_t xrrav = 0.0 ;
  Double_t yrrav = 0.0 ;
  Double_t rrrrav  = 0.0 ;

  Double_t xixi, yiyi, riri, wiriri, xold, yold ;
  
  //for (hit_counter=0; hit_counter<num_of_hits; hit_counter++) 
  for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())  
    { 
      //AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)hits->At(hit_counter);  
      AliL3ConfMapPoint* cHit = (AliL3ConfMapPoint*)fTrack->currentHit;

      xold = cHit->GetX() - xav ;
      yold = cHit->GetY() - yav ;
      //
      //-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
      //
      xi = (  cosrot * xold + sinrot * yold ) / rscale ;
      yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
      
      xixi   = xi * xi ;
      yiyi   = yi * yi ;
      riri   = xixi + yiyi ;
      wiriri = cHit->GetXYWeight() * riri ;
      
      xyav   += cHit->GetXYWeight() * xi * yi ;
      xxav   += cHit->GetXYWeight() * xixi ;
      yyav   += cHit->GetXYWeight() * yiyi ;
      
      xrrav  += wiriri * xi ;
      yrrav  += wiriri * yi ;
      rrrrav += wiriri * riri ;
    }
  //
//   Include vertex if required
//
  if (fTrack->ComesFromMainVertex() == true)
    {
	xold = fVertex->GetX() - xav ;
	yold = fVertex->GetY() - yav ;
	//
	//-->  ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
	//
	xi = (  cosrot * xold + sinrot * yold ) / rscale ;
	yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
	
	xixi   = xi * xi ;
	yiyi   = yi * yi ;
	riri   = xixi + yiyi ;
	wiriri = fVertex->GetXYWeight() * riri ;

	xyav   += fVertex->GetXYWeight() * xi * yi ;
	xxav   += fVertex->GetXYWeight() * xixi ;
	yyav   += fVertex->GetXYWeight() * yiyi ;

	xrrav  += wiriri * xi ;
	yrrav  += wiriri * yi ;
	rrrrav += wiriri * riri ;
  }
  //
  //    
  //
  //-->  DIVIDE BY WSUM TO MAKE AVERAGES
  //
  xxav    = xxav   / wsum ;
  yyav    = yyav   / wsum ;
  xrrav   = xrrav  / wsum ;
  yrrav   = yrrav  / wsum ;
  rrrrav  = rrrrav / wsum ;
  xyav    = xyav   / wsum ;

  Int_t const ntry = 5 ;
//
//-->  USE THESE TO GET THE COEFFICIENTS OF THE 4-TH ORDER POLYNIMIAL
//-->  DON'T PANIC - THE THIRD ORDER TERM IS ZERO !
//
  Double_t xrrxrr = xrrav * xrrav ;
  Double_t yrryrr = yrrav * yrrav ;
  Double_t rrrrm1 = rrrrav - 1.0  ;
  Double_t xxyy   = xxav  * yyav  ;        

  Double_t c0  =          rrrrm1*xxyy - xrrxrr*yyav - yrryrr*xxav ;
  Double_t c1  =        - rrrrm1      + xrrxrr      + yrryrr   - 4.0*xxyy ;        
  Double_t c2  =   4.0  + rrrrm1                               - 4.0*xxyy ;           
  Double_t c4  = - 4.0  ;                
//
//-->  COEFFICIENTS OF THE DERIVATIVE - USED IN NEWTON-RAPHSON ITERATIONS
//
  Double_t c2d =   2.0 * c2 ;
  Double_t c4d =   4.0 * c4 ;
//
//-->  0'TH VALUE OF LAMDA - LINEAR INTERPOLATION BETWEEN P(0) & P(YYAV)
//
//   LAMDA = YYAV * C0 / (C0 + YRRSQ * (XXAV-YYAV))
  Double_t lamda  = 0.0 ;
  Double_t dlamda = 0.0 ;
//
  Double_t chiscl = wsum * rscale * rscale ;
  Double_t dlamax = 0.001 / chiscl ;   
   
  Double_t p, pd ;
  for ( int itry = 1 ; itry <= ntry ; itry++ ) {
     p      = c0 + lamda * (c1 + lamda * (c2 + lamda * lamda * c4 )) ;
     pd     = (c1 + lamda * (c2d + lamda * lamda * c4d)) ;
     dlamda = -p / pd ;
     lamda  = lamda + dlamda ;
     if (fabs(dlamda)<   dlamax) break ;
  }

  Double_t chi2 = (Double_t)(chiscl * lamda) ;
 
  fTrack->SetChiSq1(chi2);
  // Double_t dchisq = chiscl * dlamda ;	     
//
//-->  NOW CALCULATE THE MATRIX ELEMENTS FOR ALPHA, BETA & KAPPA
//
  Double_t h11   = xxav  -     lamda ;
  Double_t h14   = xrrav ;
  Double_t h22   = yyav  -     lamda ; 
  Double_t h24   = yrrav ;
  Double_t h34   = 1.0   + 2.0*lamda ;
  if ( h11 == 0.0 || h22 == 0.0 ){
    LOG(AliL3Log::kError,"AliL3ConfMapFit::FitCircle","TrackFit")<<AliL3Log::kDec<<
      "Problems fitting circle"<<ENDLOG;
    return 1 ;
  }
  Double_t rootsq = (h14*h14)/(h11*h11) + 4.0*h34 ;

  Double_t ratio, kappa, beta ;
  if ( fabs(h22) > fabs(h24) ) {
     ratio  = h24 / h22 ;
     rootsq = ratio * ratio + rootsq ;
     kappa = 1.0 / sqrt(rootsq) ;
     beta  = - ratio * kappa ;
  }
  else {
     ratio  = h22 / h24 ;
     rootsq = 1.0 + ratio * ratio * rootsq ;
     beta  = 1.0 / sqrt(rootsq) ;
     if ( h24 > 0 ) beta = - beta ;
     kappa = -ratio * beta ;
  }            
  Double_t alpha = - (h14/h11) * kappa ;
//
//-->  transform these into the lab coordinate system
//-->  first get kappa and back to real dimensions
//
  Double_t kappa1 = kappa / rscale ;
  Double_t dbro   = 0.5   / kappa1 ;
//
//-->  next rotate alpha and beta and scale
//
  Double_t alphar = (cosrot * alpha - sinrot * beta)* dbro ;
  Double_t betar  = (sinrot * alpha + cosrot * beta)* dbro ;
//
//-->  then translate by (xav,yav)
//
  Double_t acent  = (double)(xav - alphar) ;
  Double_t bcent  = (double)(yav - betar ) ;
  Double_t radius = (double)dbro ;
//
//   Get charge
//
  Int_t q = ( ( yrrav < 0 ) ? 1 : -1 ) ;

  fTrack->SetCharge(q);
  
//
//    Get other track parameters
//
  Double_t x0, y0,phi0,r0,psi,pt ;
  if ( fTrack->ComesFromMainVertex() == true ) 
    {
      //flag = 1 ; // primary track flag
      x0   = fVertex->GetX() ;
      y0   = fVertex->GetY() ;
      phi0 = fVertex->GetPhi() ;
      r0   = fVertex->GetR() ;
      fTrack->SetPhi0(phi0);
      fTrack->SetR0(r0);
    } 
  else 
    {
      //AliL3ConfMapPoint *lHit = (AliL3ConfMapPoint*)hits->Last();
      AliL3ConfMapPoint *lHit = (AliL3ConfMapPoint*)fTrack->lastHit;
      //flag =  0 ; // primary track flag
      x0   =  lHit->GetX()  ;
      y0   =  lHit->GetY()  ;
      phi0 =  atan2(lHit->GetY(),lHit->GetX());
      if ( phi0 < 0 ) phi0 += 2*pi;
      r0   =  sqrt ( lHit->GetX() * lHit->GetX() + lHit->GetY() * lHit->GetY() )  ;
      fTrack->SetPhi0(phi0);
      fTrack->SetR0(r0);
    }
  //
  psi  = (Double_t)atan2(bcent-y0,acent-x0) ;
  psi  = psi + q * 0.5F * pi ;
  if ( psi < 0 ) psi = psi + 2*pi;
  
  pt   = (Double_t)(BFACT * bField * radius ) ;
  fTrack->SetPsi(psi);
  fTrack->SetPt(pt);

  //
//    Get errors from fast fit
//
  //if ( getPara()->getErrors ) getErrorsCircleFit ( acent, bcent, radius ) ;
//
  return 0 ;
  
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//    Fit Line in s-z plane
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Int_t AliL3ConfMapFit::FitLine ( )
{
  //
  //Initialization 
  //
  Double_t sum = 0.F ;
  Double_t ss  = 0.F ;
  Double_t sz  = 0.F ;
  Double_t sss = 0.F ;
  Double_t ssz = 0.F ;
  //
  //find sum , sums ,sumz, sumss 
  // 
  Double_t dx, dy ;
  Double_t radius = (Double_t)(fTrack->GetPt() / ( BFACT * bField ) ) ;

  //TObjArray *hits = fTrack->GetHits();
  //Int_t num_of_hits = fTrack->GetNumberOfPoints();

  if ( fTrack->ComesFromMainVertex() == true ) 
    {
      dx = ((AliL3ConfMapPoint*)fTrack->firstHit)->GetX() - fVertex->GetX();
      dy = ((AliL3ConfMapPoint*)fTrack->firstHit)->GetY() - fVertex->GetY() ;
    }
  else 
    {
      dx = ((AliL3ConfMapPoint *)fTrack->firstHit)->GetX() - ((AliL3ConfMapPoint *)fTrack->lastHit)->GetX() ;
      dy = ((AliL3ConfMapPoint *)fTrack->firstHit)->GetY() - ((AliL3ConfMapPoint *)fTrack->lastHit)->GetY() ;
      //dx = ((AliL3ConfMapPoint *)hits->First())->GetX() - ((AliL3ConfMapPoint *)hits->Last())->GetX() ;
      //dy = ((AliL3ConfMapPoint *)hits->First())->GetY() - ((AliL3ConfMapPoint *)hits->Last())->GetY() ;
    }
  
  Double_t localPsi = 0.5F * sqrt ( dx*dx + dy*dy ) / radius ;
  Double_t total_s ;
  
  if ( fabs(localPsi) < 1. ) 
    {
      total_s = 2.0 * radius * asin ( localPsi ) ;
    } 
  else 
    { 
      total_s = 2.0 * radius * pi ;
    } 
  
  AliL3ConfMapPoint *previousHit = NULL;

  // FtfBaseHit *previousHit = 0  ;
  
  //for ( startLoop() ; done() ; nextHit() ) {
  Double_t dpsi,s;

  //  for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
  for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())  
    {
      // AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)hits->At(hit_counter);
      AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)fTrack->currentHit;
      // if ( currentHit != firstHit ) 
      if(cHit != fTrack->firstHit)//  hits->First())
	{
	  dx   = cHit->GetX() - previousHit->GetX() ;
	  dy   = cHit->GetY() - previousHit->GetY() ;
	  dpsi = 0.5 * (Double_t)sqrt ( dx*dx + dy*dy ) / radius ;
	  fTrack->SetPsierr(dpsi);
	  s = previousHit->GetS() - 2.0 * radius * (Double_t)asin ( dpsi ) ;
	  cHit->SetS(s);
	}
      else
	cHit->SetS(total_s);
      //	cHit->s = total_s ;
    
      sum += cHit->GetZWeight() ;
      ss  += cHit->GetZWeight() * cHit->GetS() ;
      sz  += cHit->GetZWeight() * cHit->GetZ() ;
      sss += cHit->GetZWeight() * cHit->GetS() * cHit->GetS() ;
      ssz += cHit->GetZWeight() * cHit->GetS() * cHit->GetZ() ;
      previousHit = cHit ;
    }
  
  Double_t chi2,det = sum * sss - ss * ss;
  if ( fabs(det) < 1e-20)
    { 
      chi2 = 99999.F ;
      fTrack->SetChiSq2(chi2);
      return 0 ;
    }
  
  //Compute the best fitted parameters A,B
  Double_t tanl,z0,dtanl,dz0;

  tanl = (Double_t)((sum * ssz - ss * sz ) / det );
  z0   = (Double_t)((sz * sss - ssz * ss ) / det );

  fTrack->SetTgl(tanl);
  fTrack->SetZ0(z0);
  
  //     calculate chi-square 
  
  chi2 = 0.;
  Double_t r1 ;
  
  //for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
  for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())  
    {
      //AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)hits->At(hit_counter);
      AliL3ConfMapPoint *cHit = (AliL3ConfMapPoint*)fTrack->currentHit;
      r1   = cHit->GetZ() - tanl * cHit->GetS() - z0 ;
      chi2 += (Double_t) ( (Double_t)cHit->GetZWeight() * (r1 * r1) );
    }
  fTrack->SetChiSq2(chi2);
  //
  //     calculate estimated variance
  //      varsq=chi/(double(n)-2.) 
  //     calculate covariance matrix 
  //      siga=sqrt(varsq*sxx/det) 
  //      sigb=sqrt(varsq*sum/det) 
  //
  dtanl = (Double_t) ( sum / det );
  dz0   = (Double_t) ( sss / det );
  
  fTrack->SetTglerr(dtanl);
  fTrack->SetZ0err(dz0);
  
  return 0 ;
}
Commit	Line	Data
b661165c	1	// Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
	2	//*-- Copyright &copy ASV
	3
108615fc	4	#include <math.h>
	5
	6	#include "AliL3Logging.h"
	7	#include "AliL3ConfMapFit.h"
	8	#include "AliL3Vertex.h"
	9	#include "AliL3ConfMapTrack.h"
	10	#include "AliL3ConfMapPoint.h"
	11
b661165c	12	//_____________________________________________________________
108615fc	13	// AliL3ConfMapFit
	14	//
	15	// Fit class for conformal mapping tracking
	16
	17	ClassImp(AliL3ConfMapFit)
	18
	19	Double_t AliL3ConfMapFit::pi=3.14159265358979323846;
	20
	21	AliL3ConfMapFit::AliL3ConfMapFit(AliL3ConfMapTrack track,AliL3Vertex vertex)
	22	{
	23	//constructor
	24	fTrack = track;
	25	fVertex = vertex;
	26	BFACT = 0.0029980;
	27	bField = 0.2;
	28
	29	}
	30
	31	Int_t AliL3ConfMapFit::FitHelix()
	32	{
	33	if(FitCircle())
	34	{
	35	LOG(AliL3Log::kError,"AliL3ConfMapFit::FitHelix","TrackFit")<<AliL3Log::kDec<<
	36	"Problems during circle fit"<<ENDLOG;
	37	return 1;
	38	}
	39	if(FitLine())
	40	{
	41	LOG(AliL3Log::kError,"AliL3ConfMapFit::FitHelix","TrackFit")<<AliL3Log::kDec<<
	42	"Problems during line fit"<<ENDLOG;
	43	return 1;
	44	}
	45	return 0;
	46	}
	47
	48	Int_t AliL3ConfMapFit::FitCircle()
	49	{
	50	//-----------------------------------------------------------------
	51	//Fits circle parameters using algorithm
	52	//described by ChErnov and Oskov in Computer Physics
	53	//Communications.
	54	//
	55	//Written in FORTRAN by Jawluen Tang, Physics department , UT-Austin
	56	//Moved to C by Pablo Yepes
	57	//Moved to AliROOT by ASV.
	58	//------------------------------------------------------------------
	59
	60	Double_t wsum = 0.0 ;
	61	Double_t xav = 0.0 ;
	62	Double_t yav = 0.0 ;
	63
	64	Int_t num_of_hits = fTrack->GetNumberOfPoints();
	65	//
	66	// Loop over hits calculating average
	67	Int_t co=0;
	68
	69	for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
	70	{
	71	co++;
	72	AliL3ConfMapPoint cHit = (AliL3ConfMapPoint)fTrack->currentHit;
	73	cHit->SetXYWeight( 1./ (Double_t)(cHit->GetXerr()cHit->GetXerr() + cHit->GetYerr()cHit->GetYerr()) );
	74	wsum += cHit->GetXYWeight() ;
	75	xav += cHit->GetXYWeight() * cHit->GetX() ;
	76	yav += cHit->GetXYWeight() * cHit->GetY() ;
77	}
78	if(co!=num_of_hits)
79	LOG(AliL3Log::kError,"AliL3ConfMapFit::FitCircle","TrackFit")<<AliL3Log::kDec<<
80	"Mismatch of hits. Counter: "<<co<<" nHits: "<<num_of_hits<<ENDLOG;
81	if (fTrack->ComesFromMainVertex() == true)
82	{
83	wsum += fVertex->GetXYWeight() ;
84	xav += fVertex->GetX() ;
85	yav += fVertex->GetY() ;
86	}
87
88	xav = xav / wsum ;
89	yav = yav / wsum ;
90	//
91	// CALCULATE <X2>, <XY>, AND <Y2> WITH <X> = 0, & <Y> = 0
92	//
93	Double_t xxav = 0.0 ;
94	Double_t xyav = 0.0 ;
95	Double_t yyav = 0.0 ;
96	Double_t xi, yi ;
97
98	for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
99	{
100	//AliL3ConfMapPoint cHit = (AliL3ConfMapPoint )hits->At(hit_counter);
101	AliL3ConfMapPoint cHit = (AliL3ConfMapPoint)fTrack->currentHit;
102	xi = cHit->GetX() - xav ;
103	yi = cHit->GetY() - yav ;
104	xxav += xi * xi * cHit->GetXYWeight() ;
105	xyav += xi * yi * cHit->GetXYWeight() ;
106	yyav += yi * yi * cHit->GetXYWeight() ;
107	}
108
109	if (fTrack->ComesFromMainVertex() == true)
110	{
111	xi = fVertex->GetX() - xav ;
112	yi = fVertex->GetY() - yav ;
113	xxav += xi * xi * fVertex->GetXYWeight() ;
114	xyav += xi * yi * fVertex->GetXYWeight() ;
115	yyav += yi * yi * fVertex->GetXYWeight() ;
116	}
117	xxav = xxav / wsum ;
118	xyav = xyav / wsum ;
119	yyav = yyav / wsum ;
120	//
121	//--> ROTATE COORDINATES SO THAT <XY> = 0
122	//
123	//--> SIGN(C2 - S2) = SIGN(XXAV - YYAV) >
124	//--> & > ==> NEW : (XXAV-YYAV) > 0
125	//--> SIGN(S) = SIGN(XYAV) >
126
127	Double_t a = fabs( xxav - yyav ) ;
128	Double_t b = 4.0 * xyav * xyav ;
129
130	Double_t asqpb = a * a + b ;
131	Double_t rasqpb = sqrt ( asqpb) ;
132
133	Double_t splus = 1.0 + a / rasqpb ;
134	Double_t sminus = b / (asqpb * splus) ;
135
136	splus = sqrt (0.5 * splus ) ;
137	sminus = sqrt (0.5 * sminus) ;
138	//
139	//-> FIRST REQUIRE : SIGN(C2 - S2) = SIGN(XXAV - YYAV)
140	//
141	Double_t sinrot, cosrot ;
142	if ( xxav <= yyav ) {
143	cosrot = sminus ;
144	sinrot = splus ;
145	}
146	else {
147	cosrot = splus ;
148	sinrot = sminus ;
149	}
150	//
151	//-> REQUIRE : SIGN(S) = SIGN(XYAV) * SIGN(C) (ASSUMING SIGN(C) > 0)
152	//
153	if ( xyav < 0.0 ) sinrot = - sinrot ;
154	//
155	//--> WE NOW HAVE THE SMALLEST ANGLE THAT GUARANTEES <X2> > <Y2>
156	//--> TO GET THE SIGN OF THE CHARGE RIGHT, THE NEW X-AXIS MUST POINT
157	//--> OUTWARD FROM THE ORGIN. WE ARE FREE TO CHANGE SIGNS OF BOTH
158	//--> COSROT AND SINROT SIMULTANEOUSLY TO ACCOMPLISH THIS.
159	//
160	//--> CHOOSE SIGN OF C WISELY TO BE ABLE TO GET THE SIGN OF THE CHARGE
161	//
162	if ( cosrotxav+sinrotyav < 0.0 ) {
163	cosrot = -cosrot ;
164	sinrot = -sinrot ;
165	}
166	//
167	//-> NOW GET <R2> AND RSCALE= SQRT(<R2>)
168	//
169	Double_t rrav = xxav + yyav ;
170	Double_t rscale = sqrt(rrav) ;
171
172	xxav = 0.0 ;
173	yyav = 0.0 ;
174	xyav = 0.0 ;
175	Double_t xrrav = 0.0 ;
176	Double_t yrrav = 0.0 ;
177	Double_t rrrrav = 0.0 ;
178
179	Double_t xixi, yiyi, riri, wiriri, xold, yold ;
180
181	//for (hit_counter=0; hit_counter<num_of_hits; hit_counter++)
182	for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
183	{
184	//AliL3ConfMapPoint cHit = (AliL3ConfMapPoint)hits->At(hit_counter);
185	AliL3ConfMapPoint* cHit = (AliL3ConfMapPoint*)fTrack->currentHit;
186
187	xold = cHit->GetX() - xav ;
188	yold = cHit->GetY() - yav ;
189	//
190	//--> ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
191	//
192	xi = ( cosrot * xold + sinrot * yold ) / rscale ;
193	yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
194
195	xixi = xi * xi ;
196	yiyi = yi * yi ;
197	riri = xixi + yiyi ;
198	wiriri = cHit->GetXYWeight() * riri ;
199
200	xyav += cHit->GetXYWeight() * xi * yi ;
201	xxav += cHit->GetXYWeight() * xixi ;
202	yyav += cHit->GetXYWeight() * yiyi ;
203
204	xrrav += wiriri * xi ;
205	yrrav += wiriri * yi ;
206	rrrrav += wiriri * riri ;
207	}
208	//
209	// Include vertex if required
210	//
211	if (fTrack->ComesFromMainVertex() == true)
212	{
213	xold = fVertex->GetX() - xav ;
214	yold = fVertex->GetY() - yav ;
215	//
216	//--> ROTATE SO THAT <XY> = 0 & DIVIDE BY RSCALE SO THAT <R**2> = 1
217	//
218	xi = ( cosrot * xold + sinrot * yold ) / rscale ;
219	yi = ( -sinrot * xold + cosrot * yold ) / rscale ;
220
221	xixi = xi * xi ;
222	yiyi = yi * yi ;
223	riri = xixi + yiyi ;
224	wiriri = fVertex->GetXYWeight() * riri ;
225
226	xyav += fVertex->GetXYWeight() * xi * yi ;
227	xxav += fVertex->GetXYWeight() * xixi ;
228	yyav += fVertex->GetXYWeight() * yiyi ;
229
230	xrrav += wiriri * xi ;
231	yrrav += wiriri * yi ;
232	rrrrav += wiriri * riri ;
233	}
234	//
235	//
236	//
237	//--> DIVIDE BY WSUM TO MAKE AVERAGES
238	//
239	xxav = xxav / wsum ;
240	yyav = yyav / wsum ;
241	xrrav = xrrav / wsum ;
242	yrrav = yrrav / wsum ;
243	rrrrav = rrrrav / wsum ;
244	xyav = xyav / wsum ;
245
246	Int_t const ntry = 5 ;
247	//
248	//--> USE THESE TO GET THE COEFFICIENTS OF THE 4-TH ORDER POLYNIMIAL
249	//--> DON'T PANIC - THE THIRD ORDER TERM IS ZERO !
250	//
251	Double_t xrrxrr = xrrav * xrrav ;
252	Double_t yrryrr = yrrav * yrrav ;
253	Double_t rrrrm1 = rrrrav - 1.0 ;
254	Double_t xxyy = xxav * yyav ;
255
256	Double_t c0 = rrrrm1xxyy - xrrxrryyav - yrryrr*xxav ;
257	Double_t c1 = - rrrrm1 + xrrxrr + yrryrr - 4.0*xxyy ;
258	Double_t c2 = 4.0 + rrrrm1 - 4.0*xxyy ;
259	Double_t c4 = - 4.0 ;
260	//
261	//--> COEFFICIENTS OF THE DERIVATIVE - USED IN NEWTON-RAPHSON ITERATIONS
262	//
263	Double_t c2d = 2.0 * c2 ;
264	Double_t c4d = 4.0 * c4 ;
265	//
266	//--> 0'TH VALUE OF LAMDA - LINEAR INTERPOLATION BETWEEN P(0) & P(YYAV)
267	//
268	// LAMDA = YYAV * C0 / (C0 + YRRSQ * (XXAV-YYAV))
269	Double_t lamda = 0.0 ;
270	Double_t dlamda = 0.0 ;
271	//
272	Double_t chiscl = wsum * rscale * rscale ;
273	Double_t dlamax = 0.001 / chiscl ;
274
275	Double_t p, pd ;
276	for ( int itry = 1 ; itry <= ntry ; itry++ ) {
277	p = c0 + lamda * (c1 + lamda * (c2 + lamda * lamda * c4 )) ;
278	pd = (c1 + lamda * (c2d + lamda * lamda * c4d)) ;
279	dlamda = -p / pd ;
280	lamda = lamda + dlamda ;
281	if (fabs(dlamda)< dlamax) break ;
282	}
283
284	Double_t chi2 = (Double_t)(chiscl * lamda) ;
285
286	fTrack->SetChiSq1(chi2);
287	// Double_t dchisq = chiscl * dlamda ;
288	//
289	//--> NOW CALCULATE THE MATRIX ELEMENTS FOR ALPHA, BETA & KAPPA
290	//
291	Double_t h11 = xxav - lamda ;
292	Double_t h14 = xrrav ;
293	Double_t h22 = yyav - lamda ;
294	Double_t h24 = yrrav ;
295	Double_t h34 = 1.0 + 2.0*lamda ;
296	if ( h11 == 0.0 \|\| h22 == 0.0 ){
297	LOG(AliL3Log::kError,"AliL3ConfMapFit::FitCircle","TrackFit")<<AliL3Log::kDec<<
298	"Problems fitting circle"<<ENDLOG;
299	return 1 ;
300	}
301	Double_t rootsq = (h14h14)/(h11h11) + 4.0*h34 ;
302
303	Double_t ratio, kappa, beta ;
304	if ( fabs(h22) > fabs(h24) ) {
305	ratio = h24 / h22 ;
306	rootsq = ratio * ratio + rootsq ;
307	kappa = 1.0 / sqrt(rootsq) ;
308	beta = - ratio * kappa ;
309	}
310	else {
311	ratio = h22 / h24 ;
312	rootsq = 1.0 + ratio * ratio * rootsq ;
313	beta = 1.0 / sqrt(rootsq) ;
314	if ( h24 > 0 ) beta = - beta ;
315	kappa = -ratio * beta ;
316	}
317	Double_t alpha = - (h14/h11) * kappa ;
318	//
319	//--> transform these into the lab coordinate system
320	//--> first get kappa and back to real dimensions
321	//
322	Double_t kappa1 = kappa / rscale ;
323	Double_t dbro = 0.5 / kappa1 ;
324	//
325	//--> next rotate alpha and beta and scale
326	//
327	Double_t alphar = (cosrot * alpha - sinrot * beta)* dbro ;
328	Double_t betar = (sinrot * alpha + cosrot * beta)* dbro ;
329	//
330	//--> then translate by (xav,yav)
331	//
332	Double_t acent = (double)(xav - alphar) ;
333	Double_t bcent = (double)(yav - betar ) ;
334	Double_t radius = (double)dbro ;
335	//
336	// Get charge
337	//
5dd30052	338	Int_t q = ( ( yrrav < 0 ) ? 1 : -1 ) ;
108615fc	339
	340	fTrack->SetCharge(q);
	341
	342	//
	343	// Get other track parameters
	344	//
	345	Double_t x0, y0,phi0,r0,psi,pt ;
	346	if ( fTrack->ComesFromMainVertex() == true )
	347	{
	348	//flag = 1 ; // primary track flag
	349	x0 = fVertex->GetX() ;
	350	y0 = fVertex->GetY() ;
	351	phi0 = fVertex->GetPhi() ;
	352	r0 = fVertex->GetR() ;
	353	fTrack->SetPhi0(phi0);
	354	fTrack->SetR0(r0);
	355	}
	356	else
	357	{
	358	//AliL3ConfMapPoint lHit = (AliL3ConfMapPoint)hits->Last();
	359	AliL3ConfMapPoint lHit = (AliL3ConfMapPoint)fTrack->lastHit;
	360	//flag = 0 ; // primary track flag
	361	x0 = lHit->GetX() ;
	362	y0 = lHit->GetY() ;
	363	phi0 = atan2(lHit->GetY(),lHit->GetX());
	364	if ( phi0 < 0 ) phi0 += 2*pi;
	365	r0 = sqrt ( lHit->GetX() * lHit->GetX() + lHit->GetY() * lHit->GetY() ) ;
	366	fTrack->SetPhi0(phi0);
	367	fTrack->SetR0(r0);
	368	}
	369	//
	370	psi = (Double_t)atan2(bcent-y0,acent-x0) ;
	371	psi = psi + q * 0.5F * pi ;
	372	if ( psi < 0 ) psi = psi + 2*pi;
	373
	374	pt = (Double_t)(BFACT * bField * radius ) ;
	375	fTrack->SetPsi(psi);
	376	fTrack->SetPt(pt);
	377
	378	//
	379	// Get errors from fast fit
	380	//
	381	//if ( getPara()->getErrors ) getErrorsCircleFit ( acent, bcent, radius ) ;
	382	//
	383	return 0 ;
	384
	385	}
	386
	387	//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	388	// Fit Line in s-z plane
	389	//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	390	Int_t AliL3ConfMapFit::FitLine ( )
	391	{
	392	//
	393	//Initialization
	394	//
	395	Double_t sum = 0.F ;
	396	Double_t ss = 0.F ;
	397	Double_t sz = 0.F ;
	398	Double_t sss = 0.F ;
	399	Double_t ssz = 0.F ;
	400	//
	401	//find sum , sums ,sumz, sumss
	402	//
403	Double_t dx, dy ;
404	Double_t radius = (Double_t)(fTrack->GetPt() / ( BFACT * bField ) ) ;
405
406	//TObjArray *hits = fTrack->GetHits();
407	//Int_t num_of_hits = fTrack->GetNumberOfPoints();
408
409	if ( fTrack->ComesFromMainVertex() == true )
410	{
411	dx = ((AliL3ConfMapPoint*)fTrack->firstHit)->GetX() - fVertex->GetX();
412	dy = ((AliL3ConfMapPoint*)fTrack->firstHit)->GetY() - fVertex->GetY() ;
413	}
414	else
415	{
416	dx = ((AliL3ConfMapPoint )fTrack->firstHit)->GetX() - ((AliL3ConfMapPoint )fTrack->lastHit)->GetX() ;
417	dy = ((AliL3ConfMapPoint )fTrack->firstHit)->GetY() - ((AliL3ConfMapPoint )fTrack->lastHit)->GetY() ;
418	//dx = ((AliL3ConfMapPoint )hits->First())->GetX() - ((AliL3ConfMapPoint )hits->Last())->GetX() ;
419	//dy = ((AliL3ConfMapPoint )hits->First())->GetY() - ((AliL3ConfMapPoint )hits->Last())->GetY() ;
420	}
421
422	Double_t localPsi = 0.5F * sqrt ( dxdx + dydy ) / radius ;
423	Double_t total_s ;
424
425	if ( fabs(localPsi) < 1. )
426	{
427	total_s = 2.0 * radius * asin ( localPsi ) ;
428	}
429	else
430	{
431	total_s = 2.0 * radius * pi ;
432	}
433
434	AliL3ConfMapPoint *previousHit = NULL;
435
436	// FtfBaseHit *previousHit = 0 ;
437
438	//for ( startLoop() ; done() ; nextHit() ) {
439	Double_t dpsi,s;
440
441	// for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
442	for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
443	{
444	// AliL3ConfMapPoint cHit = (AliL3ConfMapPoint)hits->At(hit_counter);
445	AliL3ConfMapPoint cHit = (AliL3ConfMapPoint)fTrack->currentHit;
446	// if ( currentHit != firstHit )
447	if(cHit != fTrack->firstHit)// hits->First())
448	{
449	dx = cHit->GetX() - previousHit->GetX() ;
450	dy = cHit->GetY() - previousHit->GetY() ;
451	dpsi = 0.5 * (Double_t)sqrt ( dxdx + dydy ) / radius ;
452	fTrack->SetPsierr(dpsi);
453	s = previousHit->GetS() - 2.0 * radius * (Double_t)asin ( dpsi ) ;
454	cHit->SetS(s);
455	}
456	else
457	cHit->SetS(total_s);
458	// cHit->s = total_s ;
459
460	sum += cHit->GetZWeight() ;
461	ss += cHit->GetZWeight() * cHit->GetS() ;
462	sz += cHit->GetZWeight() * cHit->GetZ() ;
463	sss += cHit->GetZWeight() * cHit->GetS() * cHit->GetS() ;
464	ssz += cHit->GetZWeight() * cHit->GetS() * cHit->GetZ() ;
465	previousHit = cHit ;
466	}
467
468	Double_t chi2,det = sum * sss - ss * ss;
469	if ( fabs(det) < 1e-20)
470	{
471	chi2 = 99999.F ;
472	fTrack->SetChiSq2(chi2);
473	return 0 ;
474	}
475
476	//Compute the best fitted parameters A,B
477	Double_t tanl,z0,dtanl,dz0;
478
479	tanl = (Double_t)((sum * ssz - ss * sz ) / det );
480	z0 = (Double_t)((sz * sss - ssz * ss ) / det );
481
482	fTrack->SetTgl(tanl);
483	fTrack->SetZ0(z0);
484
485	// calculate chi-square
486
487	chi2 = 0.;
488	Double_t r1 ;
489
490	//for(hit_counter=0; hit_counter<num_of_hits; hit_counter++)
491	for(fTrack->StartLoop(); fTrack->LoopDone(); fTrack->GetNextHit())
492	{
493	//AliL3ConfMapPoint cHit = (AliL3ConfMapPoint)hits->At(hit_counter);
494	AliL3ConfMapPoint cHit = (AliL3ConfMapPoint)fTrack->currentHit;
495	r1 = cHit->GetZ() - tanl * cHit->GetS() - z0 ;
496	chi2 += (Double_t) ( (Double_t)cHit->GetZWeight() * (r1 * r1) );
497	}
498	fTrack->SetChiSq2(chi2);
499	//
500	// calculate estimated variance
501	// varsq=chi/(double(n)-2.)
502	// calculate covariance matrix
503	// siga=sqrt(varsq*sxx/det)
504	// sigb=sqrt(varsq*sum/det)
505	//
506	dtanl = (Double_t) ( sum / det );
507	dz0 = (Double_t) ( sss / det );
508
509	fTrack->SetTglerr(dtanl);
510	fTrack->SetZ0err(dz0);
511
512	return 0 ;
513	}