[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUSeed.cxx

#include <TString.h>
#include <TMath.h>
#include "AliITSUSeed.h"
using namespace TMath;

ClassImp(AliITSUSeed)

//_________________________________________________________________________
AliITSUSeed::AliITSUSeed() 
: fHitsPattern(0)
  ,fClID(0)
  ,fChi2Glo(0)
  ,fChi2Cl(0)
  ,fParent(0)
{
  // def c-tor
  ResetFMatrix();
}

//_________________________________________________________________________
AliITSUSeed::~AliITSUSeed()
{
  // d-rot
}

//_________________________________________________________________________
AliITSUSeed::AliITSUSeed(const AliITSUSeed& src) 
  :AliExternalTrackParam(src)
  ,fHitsPattern(src.fHitsPattern)
  ,fClID(src.fClID)
  ,fChi2Glo(src.fChi2Glo)
  ,fChi2Cl(src.fChi2Cl)
  ,fParent(src.fParent) 
{
  // def c-tor
  for (int i=kNFElem;i--;) fFMatrix[i] = src.fFMatrix[i];
}

//_________________________________________________________________________
AliITSUSeed &AliITSUSeed::operator=(const AliITSUSeed& src) 
{
  // def c-tor
  if (this == &src) return *this;
  fClID        = src.fClID;
  fHitsPattern = src.fHitsPattern;
  fChi2Glo     = src.fChi2Glo;
  fChi2Cl      = src.fChi2Cl;
  fParent      = src.fParent;
  for (int i=kNFElem;i--;) fFMatrix[i] = src.fFMatrix[i];
  AliExternalTrackParam::operator=(src);
  return *this;
}

//_________________________________________________________________________
void AliITSUSeed::Print(Option_t* opt) const
{
  // print seed info
  int lr,cl = GetLrCluster(lr);
  printf("%cLr%d Cl:%4d Chi2Glo:%7.2f(%7.2f) Chi2Cl:",IsKilled() ? '-':' ',
	 lr,cl,GetChi2Glo(),GetChi2GloNrm());
  cl<0 ? printf("   NA  ") : printf("%7.2f",GetChi2Cl());
  printf(" |"); 
  for (int i=0;i<=12;i++) printf("%c",HasClusterOnLayer(i) ? '+':'-'); printf("|\n");
  TString opts = opt; opts.ToLower();
  if (opts.Contains("etp")) AliExternalTrackParam::Print();
  if (opts.Contains("parent") && GetParent()) GetParent()->Print(opt);
}

//______________________________________________________________________________
Float_t AliITSUSeed::GetChi2GloNrm() const
{
  int ndf = 2*GetNLayersHit() - 5;
  return ndf>0 ? fChi2Glo/ndf : fChi2Glo;
}


//______________________________________________________________________________
Int_t AliITSUSeed::Compare(const TObject* obj)  const
{
  // compare clusters accodring to specific mode
  const AliITSUSeed* sd = (const AliITSUSeed*)obj;
  const Float_t kTol = 1e-5;
  if (!IsKilled() && sd->IsKilled()) return -1;
  if ( IsKilled() &&!sd->IsKilled()) return  1;
  //
  if      (GetChi2Glo()+kTol<sd->GetChi2Glo()) return -1;
  else if (GetChi2Glo()-kTol>sd->GetChi2Glo()) return  1;
  return 0;
}

//______________________________________________________________________________
Bool_t AliITSUSeed::IsEqual(const TObject* obj)  const
{
  // compare clusters accodring to specific mode
  const AliITSUSeed* sd = (const AliITSUSeed*)obj;
  const Float_t kTol = 1e-5;
  if (IsKilled() != sd->IsKilled()) return kFALSE;
  return Abs(GetChi2Glo() - sd->GetChi2Glo())<kTol;
}

//______________________________________________________________________________
Bool_t AliITSUSeed::PropagateToX(Double_t xk, Double_t b) 
{
  // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
  Double_t dx=xk-fX;
  if (TMath::Abs(dx)<=kAlmost0)  return kTRUE;

  Double_t crv=GetC(b);
  if (TMath::Abs(b) < kAlmost0Field) crv=0.;
  Double_t x2r = crv*dx;
  Double_t f1=fP[2], f2=f1 + x2r;
  if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
  if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
  if (TMath::Abs(fP[4])< kAlmost0) return kFALSE;

  Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4];
  Double_t 
  &fC00=fC[0],
  &fC10=fC[1],   &fC11=fC[2],  
  &fC20=fC[3],   &fC21=fC[4],   &fC22=fC[5],
  &fC30=fC[6],   &fC31=fC[7],   &fC32=fC[8],   &fC33=fC[9],  
  &fC40=fC[10],  &fC41=fC[11],  &fC42=fC[12],  &fC43=fC[13], &fC44=fC[14];

  Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
  if (TMath::Abs(r1)<kAlmost0)  return kFALSE;
  if (TMath::Abs(r2)<kAlmost0)  return kFALSE;

  fX=xk;
  double dy2dx = (f1+f2)/(r1+r2);
  fP0 += dx*dy2dx;
  if (TMath::Abs(x2r)<0.05) {
    fP1 += dx*(r2 + f2*dy2dx)*fP3;  // Many thanks to P.Hristov !
    fP2 += x2r;
  }
  else { 
    // for small dx/R the linear apporximation of the arc by the segment is OK,
    // but at large dx/R the error is very large and leads to incorrect Z propagation
    // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
    // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
    // Similarly, the rotation angle in linear in dx only for dx<<R
    double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx);   // distance from old position to new one
    double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
    fP1 += rot/crv*fP3;
    fP2  = TMath::Sin(rot + TMath::ASin(fP2));
  }

  //f = F - 1
  double r1i = 1./r1;
  double r2i = 1./r2;
  double tg1 = f1*r1i;
  double tg2 = f2*r2i;
  double v0 = 1. + dy2dx*tg2;
  double v1 = (r1i+r2i)*(dy2dx*(tg1+tg2)+2);
  double v2 = (r1i+r2i)*v0;
  //
  double f24 = dx*crv/fP4;
  double f02 = dx*v1;
  double f04 = dx*v2*f24;
  double f12 = dx*fP3*    (f2*v1+dy2dx-tg2);
  double f13 = dx*r2*v0;
  double f14 = dx*f24*fP3*(f2*v2+dy2dx-tg2);
  //
  //b = C*ft
  Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
  Double_t b02=f24*fC40;
  Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
  Double_t b12=f24*fC41;
  Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
  Double_t b22=f24*fC42;
  Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
  Double_t b42=f24*fC44;
  Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
  Double_t b32=f24*fC43;
  
  //a = f*b = f*C*ft
  Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
  Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
  Double_t a22=f24*b42;

  //F*C*Ft = C + (b + bt + a)
  fC00 += b00 + b00 + a00;
  fC10 += b10 + b01 + a01; 
  fC20 += b20 + b02 + a02;
  fC30 += b30;
  fC40 += b40;
  fC11 += b11 + b11 + a11;
  fC21 += b21 + b12 + a12;
  fC31 += b31; 
  fC41 += b41;
  fC22 += b22 + b22 + a22;
  fC32 += b32;
  fC42 += b42;
  //
  // update stored transformation matrix   F = Fnew*Fold
  fFMatrix[kF04] += f04 + f24*fFMatrix[kF02];
  fFMatrix[kF14] += f14 + f24*fFMatrix[kF12];
  fFMatrix[kF02] += f02;
  fFMatrix[kF12] += f12;
  fFMatrix[kF13] += f13;
  fFMatrix[kF24] += f24;
  //
  CheckCovariance();

  return kTRUE;
}
Commit	Line	Data
f8832015	1	#include <TString.h>
32d38de2	2	#include <TMath.h>
32d38de2	3	#include "AliITSUSeed.h"
32d38de2	4	using namespace TMath;
	5
	6	ClassImp(AliITSUSeed)
	7
	8	//_________________________________________________________________________
	9	AliITSUSeed::AliITSUSeed()
f8832015	10	: fHitsPattern(0)
	11	,fClID(0)
	12	,fChi2Glo(0)
	13	,fChi2Cl(0)
c61e50c3	14	,fParent(0)
32d38de2	15	{
32d38de2	16	// def c-tor
44785f3e	17	ResetFMatrix();
32d38de2	18	}
	19
	20	//_________________________________________________________________________
	21	AliITSUSeed::~AliITSUSeed()
	22	{
	23	// d-rot
	24	}
	25
	26	//_________________________________________________________________________
	27	AliITSUSeed::AliITSUSeed(const AliITSUSeed& src)
c61e50c3	28	:AliExternalTrackParam(src)
f8832015	29	,fHitsPattern(src.fHitsPattern)
c61e50c3	30	,fClID(src.fClID)
f8832015	31	,fChi2Glo(src.fChi2Glo)
f8832015	32	,fChi2Cl(src.fChi2Cl)
c61e50c3	33	,fParent(src.fParent)
32d38de2	34	{
32d38de2	35	// def c-tor
44785f3e	36	for (int i=kNFElem;i--;) fFMatrix[i] = src.fFMatrix[i];
32d38de2	37	}
	38
	39	//_________________________________________________________________________
	40	AliITSUSeed &AliITSUSeed::operator=(const AliITSUSeed& src)
	41	{
	42	// def c-tor
	43	if (this == &src) return *this;
f8832015	44	fClID = src.fClID;
	45	fHitsPattern = src.fHitsPattern;
	46	fChi2Glo = src.fChi2Glo;
	47	fChi2Cl = src.fChi2Cl;
	48	fParent = src.fParent;
44785f3e	49	for (int i=kNFElem;i--;) fFMatrix[i] = src.fFMatrix[i];
32d38de2	50	AliExternalTrackParam::operator=(src);
	51	return *this;
	52	}
f8832015	53
	54	//_________________________________________________________________________
	55	void AliITSUSeed::Print(Option_t* opt) const
	56	{
	57	// print seed info
	58	int lr,cl = GetLrCluster(lr);
3dd9c283	59	printf("%cLr%d Cl:%4d Chi2Glo:%7.2f(%7.2f) Chi2Cl:",IsKilled() ? '-':' ',
	60	lr,cl,GetChi2Glo(),GetChi2GloNrm());
	61	cl<0 ? printf(" NA ") : printf("%7.2f",GetChi2Cl());
f8832015	62	printf(" \|");
0091e9f0	63	for (int i=0;i<=12;i++) printf("%c",HasClusterOnLayer(i) ? '+':'-'); printf("\|\n");
f8832015	64	TString opts = opt; opts.ToLower();
	65	if (opts.Contains("etp")) AliExternalTrackParam::Print();
	66	if (opts.Contains("parent") && GetParent()) GetParent()->Print(opt);
	67	}
3dd9c283	68
	69	//______________________________________________________________________________
	70	Float_t AliITSUSeed::GetChi2GloNrm() const
	71	{
	72	int ndf = 2*GetNLayersHit() - 5;
	73	return ndf>0 ? fChi2Glo/ndf : fChi2Glo;
	74	}
	75
	76
	77	//______________________________________________________________________________
	78	Int_t AliITSUSeed::Compare(const TObject* obj) const
	79	{
	80	// compare clusters accodring to specific mode
f9c7eb32	81	const AliITSUSeed* sd = (const AliITSUSeed*)obj;
3dd9c283	82	const Float_t kTol = 1e-5;
	83	if (!IsKilled() && sd->IsKilled()) return -1;
	84	if ( IsKilled() &&!sd->IsKilled()) return 1;
	85	//
	86	if (GetChi2Glo()+kTol<sd->GetChi2Glo()) return -1;
	87	else if (GetChi2Glo()-kTol>sd->GetChi2Glo()) return 1;
	88	return 0;
	89	}
	90
	91	//______________________________________________________________________________
	92	Bool_t AliITSUSeed::IsEqual(const TObject* obj) const
	93	{
	94	// compare clusters accodring to specific mode
f9c7eb32	95	const AliITSUSeed* sd = (const AliITSUSeed*)obj;
3dd9c283	96	const Float_t kTol = 1e-5;
	97	if (IsKilled() != sd->IsKilled()) return kFALSE;
	98	return Abs(GetChi2Glo() - sd->GetChi2Glo())<kTol;
	99	}
44785f3e	100
	101	//______________________________________________________________________________
	102	Bool_t AliITSUSeed::PropagateToX(Double_t xk, Double_t b)
	103	{
	104	// Propagate this track to the plane X=xk (cm) in the field "b" (kG)
	105	Double_t dx=xk-fX;
	106	if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
	107
	108	Double_t crv=GetC(b);
	109	if (TMath::Abs(b) < kAlmost0Field) crv=0.;
	110	Double_t x2r = crv*dx;
	111	Double_t f1=fP[2], f2=f1 + x2r;
	112	if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
	113	if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
	114	if (TMath::Abs(fP[4])< kAlmost0) return kFALSE;
	115
	116	Double_t &fP0=fP[0], &fP1=fP[1], &fP2=fP[2], &fP3=fP[3], &fP4=fP[4];
	117	Double_t
	118	&fC00=fC[0],
	119	&fC10=fC[1], &fC11=fC[2],
	120	&fC20=fC[3], &fC21=fC[4], &fC22=fC[5],
	121	&fC30=fC[6], &fC31=fC[7], &fC32=fC[8], &fC33=fC[9],
	122	&fC40=fC[10], &fC41=fC[11], &fC42=fC[12], &fC43=fC[13], &fC44=fC[14];
	123
	124	Double_t r1=TMath::Sqrt((1.-f1)(1.+f1)), r2=TMath::Sqrt((1.-f2)(1.+f2));
	125	if (TMath::Abs(r1)<kAlmost0) return kFALSE;
	126	if (TMath::Abs(r2)<kAlmost0) return kFALSE;
	127
	128	fX=xk;
	129	double dy2dx = (f1+f2)/(r1+r2);
	130	fP0 += dx*dy2dx;
	131	if (TMath::Abs(x2r)<0.05) {
	132	fP1 += dx(r2 + f2dy2dx)*fP3; // Many thanks to P.Hristov !
	133	fP2 += x2r;
	134	}
	135	else {
	136	// for small dx/R the linear apporximation of the arc by the segment is OK,
	137	// but at large dx/R the error is very large and leads to incorrect Z propagation
	138	// angle traversed delta = 2asin(dist_start_end / R / 2), hence the arc is: RdeltaPhi
	139	// The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)sqrt(2+f1f2+r1*r2)
	140	// Similarly, the rotation angle in linear in dx only for dx<<R
	141	double chord = dxTMath::Sqrt(1+dy2dxdy2dx); // distance from old position to new one
	142	double rot = 2TMath::ASin(0.5chord*crv); // angular difference seen from the circle center
	143	fP1 += rot/crv*fP3;
	144	fP2 = TMath::Sin(rot + TMath::ASin(fP2));
	145	}
	146
	147	//f = F - 1
	148	double r1i = 1./r1;
	149	double r2i = 1./r2;
	150	double tg1 = f1*r1i;
	151	double tg2 = f2*r2i;
	152	double v0 = 1. + dy2dx*tg2;
	153	double v1 = (r1i+r2i)(dy2dx(tg1+tg2)+2);
	154	double v2 = (r1i+r2i)*v0;
	155	//
	156	double f24 = dx*crv/fP4;
	157	double f02 = dx*v1;
	158	double f04 = dxv2f24;
	159	double f12 = dxfP3 (f2*v1+dy2dx-tg2);
	160	double f13 = dxr2v0;
	161	double f14 = dxf24fP3(f2v2+dy2dx-tg2);
	162	//
	163	//b = C*ft
164	Double_t b00=f02fC20 + f04fC40, b01=f12fC20 + f14fC40 + f13*fC30;
165	Double_t b02=f24*fC40;
166	Double_t b10=f02fC21 + f04fC41, b11=f12fC21 + f14fC41 + f13*fC31;
167	Double_t b12=f24*fC41;
168	Double_t b20=f02fC22 + f04fC42, b21=f12fC22 + f14fC42 + f13*fC32;
169	Double_t b22=f24*fC42;
170	Double_t b40=f02fC42 + f04fC44, b41=f12fC42 + f14fC44 + f13*fC43;
171	Double_t b42=f24*fC44;
172	Double_t b30=f02fC32 + f04fC43, b31=f12fC32 + f14fC43 + f13*fC33;
173	Double_t b32=f24*fC43;
174
175	//a = fb = fC*ft
176	Double_t a00=f02b20+f04b40,a01=f02b21+f04b41,a02=f02b22+f04b42;
177	Double_t a11=f12b21+f14b41+f13b31,a12=f12b22+f14b42+f13b32;
178	Double_t a22=f24*b42;
179
180	//FCFt = C + (b + bt + a)
181	fC00 += b00 + b00 + a00;
182	fC10 += b10 + b01 + a01;
183	fC20 += b20 + b02 + a02;
184	fC30 += b30;
185	fC40 += b40;
186	fC11 += b11 + b11 + a11;
187	fC21 += b21 + b12 + a12;
188	fC31 += b31;
189	fC41 += b41;
190	fC22 += b22 + b22 + a22;
191	fC32 += b32;
192	fC42 += b42;
193	//
194	// update stored transformation matrix F = Fnew*Fold
195	fFMatrix[kF04] += f04 + f24*fFMatrix[kF02];
196	fFMatrix[kF14] += f14 + f24*fFMatrix[kF12];
197	fFMatrix[kF02] += f02;
198	fFMatrix[kF12] += f12;
199	fFMatrix[kF13] += f13;
200	fFMatrix[kF24] += f24;
201	//
202	CheckCovariance();
203
204	return kTRUE;
205	}