CheckCovariance();
}
+//_____________________________________________________________________________
+void AliExternalTrackParam::CopyFromVTrack(const AliVTrack *vTrack)
+{
+ //
+ // Recreate TrackParams from VTrack
+ // This is not a copy contructor !
+ //
+ if (!vTrack) {
+ AliError("Source VTrack is NULL");
+ return;
+ }
+ if (this==vTrack) {
+ AliError("Copy of itself is requested");
+ return;
+ }
+ //
+ if (vTrack->InheritsFrom(AliExternalTrackParam::Class())) {
+ AliDebug(1,"Source itself is AliExternalTrackParam, using assignment operator");
+ *this = *(AliExternalTrackParam*)vTrack;
+ return;
+ }
+ //
+ AliVTrack::operator=( *vTrack );
+ //
+ Double_t xyz[3],pxpypz[3],cv[21];
+ vTrack->GetXYZ(xyz);
+ pxpypz[0]=vTrack->Px();
+ pxpypz[1]=vTrack->Py();
+ pxpypz[2]=vTrack->Pz();
+ vTrack->GetCovarianceXYZPxPyPz(cv);
+ Short_t sign = (Short_t)vTrack->Charge();
+ Set(xyz,pxpypz,cv,sign);
+}
+
//_____________________________________________________________________________
AliExternalTrackParam::AliExternalTrackParam(const AliVTrack *vTrack) :
AliVTrack(),
Set(xyz,pxpypz,cv,sign);
}
+/*
//_____________________________________________________________________________
void AliExternalTrackParam::Set(Double_t xyz[3],Double_t pxpypz[3],
Double_t cv[21],Short_t sign)
//
Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
// protection: avoid alpha being too close to 0 or +-pi/2
- if (TMath::Abs(sn)<kSafe) {
- fAlpha = kSafe;
+ if (TMath::Abs(sn)<2*kSafe) {
+ if (fAlpha>0) fAlpha += fAlpha< TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
+ else fAlpha += fAlpha>-TMath::Pi()/2. ? -2*kSafe : 2*kSafe;
cs=TMath::Cos(fAlpha);
sn=TMath::Sin(fAlpha);
}
- else if (cs<kSafe) {
- fAlpha -= TMath::Sign(kSafe, fAlpha);
+ else if (TMath::Abs(cs)<2*kSafe) {
+ if (fAlpha>0) fAlpha += fAlpha> TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
+ else fAlpha += fAlpha>-TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
cs=TMath::Cos(fAlpha);
- sn=TMath::Sin(fAlpha);
+ sn=TMath::Sin(fAlpha);
}
// Get the vertex of origin and the momentum
TVector3 ver(xyz[0],xyz[1],xyz[2]);
ver.RotateZ(-fAlpha);
mom.RotateZ(-fAlpha);
+ //
// x of the reference plane
fX = ver.X();
fP[2] = TMath::Sin(mom.Phi());
fP[3] = mom.Pz()/mom.Pt();
fP[4] = TMath::Sign(1/mom.Pt(),charge);
-
+ //
+ if (TMath::Abs( 1-fP[2]) < 3*kSafe) fP[2] = 1.- 3*kSafe; //Protection
+ else if (TMath::Abs(-1-fP[2]) < 3*kSafe) fP[2] =-1.+ 3*kSafe; //Protection
+ //
// Covariance matrix (formulas to be simplified)
-
- if (TMath::Abs( 1-fP[2]) < kSafe) fP[2] = 1.- kSafe; //Protection
- else if (TMath::Abs(-1-fP[2]) < kSafe) fP[2] =-1.+ kSafe; //Protection
-
Double_t pt=1./TMath::Abs(fP[4]);
- Double_t r=TMath::Sqrt((1.-fP[2])*(1.+fP[2]));
-
+ // avoid alpha+phi being to close to +-pi/2 in the cov.matrix evaluation
+ double fp2 = fP[2];
+ Double_t r=TMath::Sqrt((1.-fp2)*(1.+fp2));
+ //
Double_t m00=-sn;// m10=cs;
- Double_t m23=-pt*(sn + fP[2]*cs/r), m43=-pt*pt*(r*cs - fP[2]*sn);
- Double_t m24= pt*(cs - fP[2]*sn/r), m44=-pt*pt*(r*sn + fP[2]*cs);
+ Double_t m23=-pt*(sn + fp2*cs/r), m43=-pt*pt*(r*cs - fp2*sn);
+ Double_t m24= pt*(cs - fp2*sn/r), m44=-pt*pt*(r*sn + fp2*cs);
Double_t m35=pt, m45=-pt*pt*fP[3];
m43*=GetSign();
Double_t a1=cv[13]-cv[9]*(m23*m44+m43*m24)/m23/m43;
Double_t a2=m23*m24-m23*(m23*m44+m43*m24)/m43;
Double_t a3=m43*m44-m43*(m23*m44+m43*m24)/m23;
- Double_t a4=cv[14]-2.*cv[9]*m24*m44/m23/m43;
+ Double_t a4=cv[14]+2.*cv[9]; //cv[14]-2.*cv[9]*m24*m44/m23/m43;
Double_t a5=m24*m24-2.*m24*m44*m23/m43;
Double_t a6=m44*m44-2.*m24*m44*m43/m23;
fC[0 ] = cv[0 ]+cv[2 ];
fC[1 ] = TMath::Sign(cv34,cv[3 ]/m00);
fC[2 ] = cv[5 ];
- fC[3 ] = (cv[10]/m44-cv[6]/m43)/(m24/m44-m23/m43)/m00;
+ fC[3 ] = (cv[10]*m43-cv[6]*m44)/(m24*m43-m23*m44)/m00;
fC[10] = (cv[6]/m00-fC[3 ]*m23)/m43;
fC[6 ] = (cv[15]/m00-fC[10]*m45)/m35;
- fC[4 ] = (cv[12]-cv[8]*m44/m43)/(m24-m23*m44/m43);
+ fC[4 ] = (cv[12]*m43-cv[8]*m44)/(m24*m43-m23*m44);
fC[11] = (cv[8]-fC[4]*m23)/m43;
fC[7 ] = cv[17]/m35-fC[11]*m45/m35;
- fC[5 ] = TMath::Abs((a4-a6*a1/a3)/(a5-a6*a2/a3));
- fC[14] = TMath::Abs(a1/a3-a2*fC[5]/a3);
+ fC[5 ] = TMath::Abs((a4*a3-a6*a1)/(a5*a3-a6*a2));
+ fC[14] = TMath::Abs((a1-a2*fC[5])/a3);
fC[12] = (cv[9]-fC[5]*m23*m23-fC[14]*m43*m43)/m23/m43;
Double_t b1=cv[18]-fC[12]*m23*m45-fC[14]*m43*m45;
Double_t b2=m23*m35;
return;
}
+*/
+
+//_____________________________________________________________________________
+void AliExternalTrackParam::Set(Double_t xyz[3],Double_t pxpypz[3],
+ Double_t cv[21],Short_t sign)
+{
+ //
+ // create external track parameters from the global parameters
+ // x,y,z,px,py,pz and their 6x6 covariance matrix
+ // A.Dainese 10.10.08
+
+ // Calculate alpha: the rotation angle of the corresponding local system.
+ //
+ // For global radial position inside the beam pipe, alpha is the
+ // azimuthal angle of the momentum projected on (x,y).
+ //
+ // For global radial position outside the ITS, alpha is the
+ // azimuthal angle of the centre of the TPC sector in which the point
+ // xyz lies
+ //
+ const double kSafe = 1e-5;
+ Double_t radPos2 = xyz[0]*xyz[0]+xyz[1]*xyz[1];
+ Double_t radMax = 45.; // approximately ITS outer radius
+ if (radPos2 < radMax*radMax) { // inside the ITS
+ fAlpha = TMath::ATan2(pxpypz[1],pxpypz[0]);
+ } else { // outside the ITS
+ Float_t phiPos = TMath::Pi()+TMath::ATan2(-xyz[1], -xyz[0]);
+ fAlpha =
+ TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10);
+ }
+ //
+ Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
+ // protection: avoid alpha being too close to 0 or +-pi/2
+ if (TMath::Abs(sn)<2*kSafe) {
+ if (fAlpha>0) fAlpha += fAlpha< TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
+ else fAlpha += fAlpha>-TMath::Pi()/2. ? -2*kSafe : 2*kSafe;
+ cs=TMath::Cos(fAlpha);
+ sn=TMath::Sin(fAlpha);
+ }
+ else if (TMath::Abs(cs)<2*kSafe) {
+ if (fAlpha>0) fAlpha += fAlpha> TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
+ else fAlpha += fAlpha>-TMath::Pi()/2. ? 2*kSafe : -2*kSafe;
+ cs=TMath::Cos(fAlpha);
+ sn=TMath::Sin(fAlpha);
+ }
+ // Get the vertex of origin and the momentum
+ TVector3 ver(xyz[0],xyz[1],xyz[2]);
+ TVector3 mom(pxpypz[0],pxpypz[1],pxpypz[2]);
+ //
+ // Rotate to the local coordinate system
+ ver.RotateZ(-fAlpha);
+ mom.RotateZ(-fAlpha);
+
+ //
+ // x of the reference plane
+ fX = ver.X();
+
+ Double_t charge = (Double_t)sign;
+
+ fP[0] = ver.Y();
+ fP[1] = ver.Z();
+ fP[2] = TMath::Sin(mom.Phi());
+ fP[3] = mom.Pz()/mom.Pt();
+ fP[4] = TMath::Sign(1/mom.Pt(),charge);
+ //
+ if (TMath::Abs( 1-fP[2]) < kSafe) fP[2] = 1.- kSafe; //Protection
+ else if (TMath::Abs(-1-fP[2]) < kSafe) fP[2] =-1.+ kSafe; //Protection
+ //
+ // Covariance matrix (formulas to be simplified)
+ Double_t pt=1./TMath::Abs(fP[4]);
+ Double_t r=TMath::Sqrt((1.-fP[2])*(1.+fP[2]));
+ //
+ Double_t cv34 = TMath::Sqrt(cv[3 ]*cv[3 ]+cv[4 ]*cv[4 ]);
+ //
+ Int_t special = 0;
+ double sgcheck = r*sn + fP[2]*cs;
+ if (TMath::Abs(sgcheck)>=1-kSafe) { // special case: lab phi is +-pi/2
+ special = 1;
+ sgcheck = TMath::Sign(1.0,sgcheck);
+ }
+ else if (TMath::Abs(sgcheck)<kSafe) {
+ sgcheck = TMath::Sign(1.0,cs);
+ special = 2; // special case: lab phi is 0
+ }
+ //
+ fC[0 ] = cv[0 ]+cv[2 ];
+ fC[1 ] = TMath::Sign(cv34,-cv[3 ]*sn);
+ fC[2 ] = cv[5 ];
+ //
+ if (special==1) {
+ double pti = 1./pt;
+ double pti2 = pti*pti;
+ int q = GetSign();
+ fC[3 ] = cv[6]*pti;
+ fC[4 ] = -sgcheck*cv[8]*r*pti;
+ fC[5 ] = TMath::Abs(cv[9]*r*r*pti2);
+ fC[6 ] = (cv[10]*fP[3]-sgcheck*cv[15])*pti/r;
+ fC[7 ] = (cv[17]-sgcheck*cv[12]*fP[3])*pti;
+ fC[8 ] = (-sgcheck*cv[18]+cv[13]*fP[3])*r*pti2;
+ fC[9 ] = TMath::Abs( cv[20]-2*sgcheck*cv[19]*fP[3]+cv[14]*fP[3]*fP[3])*pti2;
+ fC[10] = cv[10]*pti2/r*q;
+ fC[11] = -sgcheck*cv[12]*pti2*q;
+ fC[12] = cv[13]*r*pti*pti2*q;
+ fC[13] = (-sgcheck*cv[19]+cv[14]*fP[3])*r*pti2*pti;
+ fC[14] = TMath::Abs(cv[14]*pti2*pti2);
+ } else if (special==2) {
+ double pti = 1./pt;
+ double pti2 = pti*pti;
+ int q = GetSign();
+ fC[3 ] = -cv[10]*pti*cs/sn;
+ fC[4 ] = cv[12]*cs*pti;
+ fC[5 ] = TMath::Abs(cv[14]*cs*cs*pti2);
+ fC[6 ] = (sgcheck*cv[6]*fP[3]-cv[15])*pti/sn;
+ fC[7 ] = (cv[17]-sgcheck*cv[8]*fP[3])*pti;
+ fC[8 ] = (cv[19]-sgcheck*cv[13]*fP[3])*cs*pti2;
+ fC[9 ] = TMath::Abs( cv[20]-2*sgcheck*cv[18]*fP[3]+cv[9]*fP[3]*fP[3])*pti2;
+ fC[10] = sgcheck*cv[6]*pti2/sn*q;
+ fC[11] = -sgcheck*cv[8]*pti2*q;
+ fC[12] = -sgcheck*cv[13]*cs*pti*pti2*q;
+ fC[13] = (-sgcheck*cv[18]+cv[9]*fP[3])*pti2*pti*q;
+ fC[14] = TMath::Abs(cv[9]*pti2*pti2);
+ }
+ else {
+ Double_t m00=-sn;// m10=cs;
+ Double_t m23=-pt*(sn + fP[2]*cs/r), m43=-pt*pt*(r*cs - fP[2]*sn);
+ Double_t m24= pt*(cs - fP[2]*sn/r), m44=-pt*pt*(r*sn + fP[2]*cs);
+ Double_t m35=pt, m45=-pt*pt*fP[3];
+ //
+ m43*=GetSign();
+ m44*=GetSign();
+ m45*=GetSign();
+ //
+ Double_t a1=cv[13]-cv[9]*(m23*m44+m43*m24)/m23/m43;
+ Double_t a2=m23*m24-m23*(m23*m44+m43*m24)/m43;
+ Double_t a3=m43*m44-m43*(m23*m44+m43*m24)/m23;
+ Double_t a4=cv[14]+2.*cv[9]; //cv[14]-2.*cv[9]*m24*m44/m23/m43;
+ Double_t a5=m24*m24-2.*m24*m44*m23/m43;
+ Double_t a6=m44*m44-2.*m24*m44*m43/m23;
+ //
+ fC[3 ] = (cv[10]*m43-cv[6]*m44)/(m24*m43-m23*m44)/m00;
+ fC[10] = (cv[6]/m00-fC[3 ]*m23)/m43;
+ fC[6 ] = (cv[15]/m00-fC[10]*m45)/m35;
+ fC[4 ] = (cv[12]*m43-cv[8]*m44)/(m24*m43-m23*m44);
+ fC[11] = (cv[8]-fC[4]*m23)/m43;
+ fC[7 ] = cv[17]/m35-fC[11]*m45/m35;
+ fC[5 ] = TMath::Abs((a4*a3-a6*a1)/(a5*a3-a6*a2));
+ fC[14] = TMath::Abs((a1-a2*fC[5])/a3);
+ fC[12] = (cv[9]-fC[5]*m23*m23-fC[14]*m43*m43)/m23/m43;
+ Double_t b1=cv[18]-fC[12]*m23*m45-fC[14]*m43*m45;
+ Double_t b2=m23*m35;
+ Double_t b3=m43*m35;
+ Double_t b4=cv[19]-fC[12]*m24*m45-fC[14]*m44*m45;
+ Double_t b5=m24*m35;
+ Double_t b6=m44*m35;
+ fC[8 ] = (b4-b6*b1/b3)/(b5-b6*b2/b3);
+ fC[13] = b1/b3-b2*fC[8]/b3;
+ fC[9 ] = TMath::Abs((cv[20]-fC[14]*(m45*m45)-fC[13]*2.*m35*m45)/(m35*m35));
+ }
+ CheckCovariance();
+
+ return;
+}
//_____________________________________________________________________________
void AliExternalTrackParam::Reset() {
// "xTimesRho" - is the product length*density (g/cm^2).
// It should be passed as negative when propagating tracks
// from the intreaction point to the outside of the central barrel.
- // "mass" - the mass of this particle (GeV/c^2).
+ // "mass" - the mass of this particle (GeV/c^2). Negative mass means charge=2 particle
// "dEdx" - mean enery loss (GeV/(g/cm^2)
// "anglecorr" - switch for the angular correction
//------------------------------------------------------------------
}
Double_t p=GetP();
+ if (mass<0) p += p; // q=2 particle
Double_t p2=p*p;
Double_t beta2=p2/(p2 + mass*mass);
double lt = 1+0.038*TMath::Log(TMath::Abs(xOverX0));
if (lt>0) theta2 *= lt*lt;
}
+ if (mass<0) theta2 *= 4; // q=2 particle
if(theta2>TMath::Pi()*TMath::Pi()) return kFALSE;
cC22 = theta2*((1.-fP2)*(1.+fP2))*(1. + fP3*fP3);
cC33 = theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
Double_t dE=dEdx*xTimesRho;
Double_t e=TMath::Sqrt(p2 + mass*mass);
if ( TMath::Abs(dE) > 0.3*e ) return kFALSE; //30% energy loss is too much!
- //cP4 = (1.- e/p2*dE);
if ( (1.+ dE/p2*(dE + 2*e)) < 0. ) return kFALSE;
cP4 = 1./TMath::Sqrt(1.+ dE/p2*(dE + 2*e)); //A precise formula by Ruben !
if (TMath::Abs(fP4*cP4)>100.) return kFALSE; //Do not track below 10 MeV/c
// "xTimesRho" - is the product length*density (g/cm^2).
// It should be passed as negative when propagating tracks
// from the intreaction point to the outside of the central barrel.
- // "mass" - the mass of this particle (GeV/c^2).
+ // "mass" - the mass of this particle (GeV/c^2). mass<0 means charge=2
// "anglecorr" - switch for the angular correction
// "Bethe" - function calculating the energy loss (GeV/(g/cm^2))
//------------------------------------------------------------------
-
+
Double_t bg=GetP()/mass;
+ if (mass<0) {
+ if (mass<-990) {
+ AliDebug(2,Form("Mass %f corresponds to unknown PID particle",mass));
+ return kFALSE;
+ }
+ bg = -2*bg;
+ }
Double_t dEdx=Bethe(bg);
+ if (mass<0) dEdx *= 4;
return CorrectForMeanMaterialdEdx(xOverX0,xTimesRho,mass,dEdx,anglecorr);
}
// "xTimesRho" - is the product length*density (g/cm^2).
// It should be passed as negative when propagating tracks
// from the intreaction point to the outside of the central barrel.
- // "mass" - the mass of this particle (GeV/c^2).
+ // "mass" - the mass of this particle (GeV/c^2). mass<0 means charge=2 particle
// "density" - mean density (g/cm^3)
// "zOverA" - mean Z/A
// "exEnergy" - mean exitation energy (GeV)
//------------------------------------------------------------------
Double_t bg=GetP()/mass;
+ if (mass<0) {
+ if (mass<-990) {
+ AliDebug(2,Form("Mass %f corresponds to unknown PID particle",mass));
+ return kFALSE;
+ }
+ bg = -2*bg;
+ }
Double_t dEdx=BetheBlochGeant(bg,density,jp1,jp2,exEnergy,zOverA);
+ if (mass<0) dEdx *= 4;
return CorrectForMeanMaterialdEdx(xOverX0,xTimesRho,mass,dEdx,anglecorr);
}
Double_t x=fX;
Double_t ca=TMath::Cos(alpha-fAlpha), sa=TMath::Sin(alpha-fAlpha);
Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2)); // Improve precision
-
// RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
// direction in local frame is along the X axis
if ((cf*ca+sf*sa)<0) {
- AliWarning(Form("Rotation failed: local cos(phi) would become %.2f",cf*ca+sf*sa));
+ AliDebug(1,Form("Rotation failed: local cos(phi) would become %.2f",cf*ca+sf*sa));
return kFALSE;
}
//
Double_t tmp=sf*ca - cf*sa;
+
if (TMath::Abs(tmp) >= kAlmost1) {
if (TMath::Abs(tmp) > 1.+ Double_t(FLT_EPSILON))
AliWarning(Form("Rotation failed ! %.10e",tmp));
return kFALSE;
}
-
fAlpha = alpha;
fX = x*ca + fP0*sa;
fP0= -x*sa + fP0*ca;
return kTRUE;
}
+//______________________________________________________
+Bool_t AliExternalTrackParam::RotateParamOnly(Double_t alpha)
+{
+ // rotate to new frame, ignore covariance
+ if (TMath::Abs(fP[2]) >= kAlmost1) {
+ AliError(Form("Precondition is not satisfied: |sin(phi)|>1 ! %f",fP[2]));
+ return kFALSE;
+ }
+ //
+ if (alpha < -TMath::Pi()) alpha += 2*TMath::Pi();
+ else if (alpha >= TMath::Pi()) alpha -= 2*TMath::Pi();
+ //
+ Double_t &fP0=fP[0];
+ Double_t &fP2=fP[2];
+ //
+ Double_t x=fX;
+ Double_t ca=TMath::Cos(alpha-fAlpha), sa=TMath::Sin(alpha-fAlpha);
+ Double_t sf=fP2, cf=TMath::Sqrt((1.- fP2)*(1.+fP2)); // Improve precision
+ // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
+ // direction in local frame is along the X axis
+ if ((cf*ca+sf*sa)<0) {
+ AliDebug(1,Form("Rotation failed: local cos(phi) would become %.2f",cf*ca+sf*sa));
+ return kFALSE;
+ }
+ //
+ Double_t tmp=sf*ca - cf*sa;
+
+ if (TMath::Abs(tmp) >= kAlmost1) {
+ if (TMath::Abs(tmp) > 1.+ Double_t(FLT_EPSILON))
+ AliWarning(Form("Rotation failed ! %.10e",tmp));
+ return kFALSE;
+ }
+ fAlpha = alpha;
+ fX = x*ca + fP0*sa;
+ fP0= -x*sa + fP0*ca;
+ fP2= tmp;
+ return kTRUE;
+}
+
+Bool_t AliExternalTrackParam::Invert() {
+ //------------------------------------------------------------------
+ // Transform this track to the local coord. system rotated by 180 deg.
+ //------------------------------------------------------------------
+ fX = -fX;
+ fAlpha += TMath::Pi();
+ while (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
+ while (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
+ //
+ fP[0] = -fP[0];
+ //fP[2] = -fP[2];
+ fP[3] = -fP[3];
+ fP[4] = -fP[4];
+ //
+ fC[1] = -fC[1]; // since the fP1 and fP2 are not inverted, their covariances with others change sign
+ fC[3] = -fC[3];
+ fC[7] = -fC[7];
+ fC[8] = -fC[8];
+ fC[11] = -fC[11];
+ fC[12] = -fC[12];
+ //
+ return kTRUE;
+}
+
Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t b) {
//----------------------------------------------------------------
// Propagate this track to the plane X=xk (cm) in the field "b" (kG)
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;
- }
+ fP2 += x2r;
+ if (TMath::Abs(x2r)<0.05) fP1 += dx*(r2 + f2*dy2dx)*fP3; // Many thanks to P.Hristov !
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));
+ // 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;
+ //
+ double rot = TMath::ASin(r1*f2 - r2*f1); // more economic version from Yura.
+ if (f1*f1+f2*f2>1 && f1*f2<0) { // special cases of large rotations or large abs angles
+ if (f2>0) rot = TMath::Pi() - rot; //
+ else rot = -TMath::Pi() - rot;
+ }
+ fP1 += fP3/crv*rot;
}
//f = F - 1
-
+ /*
Double_t f02= dx/(r1*r1*r1); Double_t cc=crv/fP4;
Double_t f04=0.5*dx*dx/(r1*r1*r1); f04*=cc;
Double_t f12= dx*fP3*f1/(r1*r1*r1);
Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1); f14*=cc;
Double_t f13= dx/r1;
Double_t f24= dx; f24*=cc;
-
+ */
+ Double_t rinv = 1./r1;
+ Double_t r3inv = rinv*rinv*rinv;
+ Double_t f24= x2r/fP4;
+ Double_t f02= dx*r3inv;
+ Double_t f04=0.5*f24*f02;
+ Double_t f12= f02*fP3*f1;
+ Double_t f14=0.5*f24*f02*fP3*f1;
+ Double_t f13= dx*rinv;
+
//b = C*ft
Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
Double_t b02=f24*fC40;
return kTRUE;
}
+Bool_t AliExternalTrackParam::PropagateParamOnlyTo(Double_t xk, Double_t b) {
+ //----------------------------------------------------------------
+ // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
+ // Only parameters are propagated, not the matrix. To be used for small
+ // distances only (<mm, i.e. misalignment)
+ //----------------------------------------------------------------
+ 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 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);
+ fP[0] += dx*dy2dx;
+ fP[2] += x2r;
+ if (TMath::Abs(x2r)<0.05) fP[1] += dx*(r2 + f2*dy2dx)*fP[3]; // Many thanks to P.Hristov !
+ 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)
+ // 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;
+ //
+ double rot = TMath::ASin(r1*f2 - r2*f1); // more economic version from Yura.
+ if (f1*f1+f2*f2>1 && f1*f2<0) { // special cases of large rotations or large abs angles
+ if (f2>0) rot = TMath::Pi() - rot; //
+ else rot = -TMath::Pi() - rot;
+ }
+ fP[1] += fP[3]/crv*rot;
+ }
+ return kTRUE;
+}
+
Bool_t
AliExternalTrackParam::Propagate(Double_t alpha, Double_t x, Double_t b) {
//------------------------------------------------------------------
}
Double_t
-AliExternalTrackParam::GetPredictedChi2(Double_t p[2],Double_t cov[3]) const {
+AliExternalTrackParam::GetPredictedChi2(const Double_t p[2],const Double_t cov[3]) const {
//----------------------------------------------------------------
// Estimate the chi2 of the space point "p" with the cov. matrix "cov"
//----------------------------------------------------------------
}
Double_t AliExternalTrackParam::
-GetPredictedChi2(Double_t p[3],Double_t covyz[3],Double_t covxyz[3]) const {
+GetPredictedChi2(const Double_t p[3],const Double_t covyz[3],const Double_t covxyz[3]) const {
//----------------------------------------------------------------
// Estimate the chi2 of the 3D space point "p" and
// the full covariance matrix "covyz" and "covxyz"
// and estimated at the same reference plane.
//----------------------------------------------------------------
- if (TMath::Abs(1. - t->GetAlpha()/GetAlpha()) > FLT_EPSILON) {
+ if (TMath::Abs(t->GetAlpha()-GetAlpha()) > FLT_EPSILON) {
AliError("The reference systems of the tracks differ !");
return kVeryBig;
}
- if (TMath::Abs(1. - t->GetX()/GetX()) > FLT_EPSILON) {
+ if (TMath::Abs(t->GetX()-GetX()) > FLT_EPSILON) {
AliError("The reference of the tracks planes differ !");
return kVeryBig;
}
return res;
}
-Bool_t AliExternalTrackParam::Update(Double_t p[2], Double_t cov[3]) {
+Bool_t AliExternalTrackParam::Update(const Double_t p[2], const Double_t cov[3]) {
//------------------------------------------------------------------
// Update the track parameters with the space point "p" having
// the covariance matrix "cov"
fC33-=k30*c03+k31*c13;
fC43-=k30*c04+k31*c14;
-
+
fC44-=k40*c04+k41*c14;
CheckCovariance();
return phi;
}
+Double_t AliExternalTrackParam::PhiPos() const {
+ //---------------------------------------------------------------------
+ // Returns the azimuthal angle of position
+ // 0 <= phi < 2*pi
+ //---------------------------------------------------------------------
+ Double_t r[3]={0.,0.,0.};
+ GetXYZ(r);
+ Double_t phi=TMath::ATan2(r[1],r[0]);
+ if (phi<0.) phi+=2.*TMath::Pi();
+
+ return phi;
+}
+
Double_t AliExternalTrackParam::M() const {
// return particle mass
Double_t dx=x-fX;
if(TMath::Abs(dx)<=kAlmost0) {z=fP[1]; return kTRUE;}
- Double_t f1=fP[2], f2=f1 + dx*GetC(b);
+ Double_t crv=GetC(b);
+ 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;
Double_t r1=sqrt((1.-f1)*(1.+f1)), r2=sqrt((1.-f2)*(1.+f2));
- z = fP[1] + dx*(r2 + f2*(f1+f2)/(r1+r2))*fP[3]; // Many thanks to P.Hristov !
+ double dy2dx = (f1+f2)/(r1+r2);
+ if (TMath::Abs(x2r)<0.05) {
+ z = fP[1] + dx*(r2 + f2*dy2dx)*fP[3]; // Many thanks to P.Hristov !
+ }
+ 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
+ z = fP[1] + rot/crv*fP[3];
+ }
return kTRUE;
}
Double_t dx=x-fX;
if(TMath::Abs(dx)<=kAlmost0) return GetXYZ(r);
- Double_t f1=fP[2], f2=f1 + dx*GetC(b);
+ Double_t crv=GetC(b);
+ Double_t x2r = crv*dx;
+ Double_t f1=fP[2], f2=f1 + dx*crv;
if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
+ double dy2dx = (f1+f2)/(r1+r2);
r[0] = x;
- r[1] = fP[0] + dx*(f1+f2)/(r1+r2);
- r[2] = fP[1] + dx*(r2 + f2*(f1+f2)/(r1+r2))*fP[3];//Thanks to Andrea & Peter
+ r[1] = fP[0] + dx*dy2dx;
+ if (TMath::Abs(x2r)<0.05) {
+ r[2] = fP[1] + dx*(r2 + f2*dy2dx)*fP[3];//Thanks to Andrea & Peter
+ }
+ 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
+ r[2] = fP[1] + rot/crv*fP[3];
+ }
return Local2GlobalPosition(r,fAlpha);
}
Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
//f = F - 1
+ /*
Double_t f02= dx/(r1*r1*r1); Double_t cc=crv/fP4;
Double_t f04=0.5*dx*dx/(r1*r1*r1); f04*=cc;
Double_t f12= dx*fP3*f1/(r1*r1*r1);
Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1); f14*=cc;
Double_t f13= dx/r1;
Double_t f24= dx; f24*=cc;
-
+ */
+ Double_t rinv = 1./r1;
+ Double_t r3inv = rinv*rinv*rinv;
+ Double_t f24= x2r/fP4;
+ Double_t f02= dx*r3inv;
+ Double_t f04=0.5*f24*f02;
+ Double_t f12= f02*fP3*f1;
+ Double_t f14=0.5*f24*f02*fP3*f1;
+ Double_t f13= dx*rinv;
+
//b = C*ft
Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
Double_t b02=f24*fC40;
return kTRUE;
}
+Bool_t AliExternalTrackParam::PropagateParamOnlyBxByBzTo(Double_t xk, const Double_t b[3]) {
+ //----------------------------------------------------------------
+ // Extrapolate this track params (w/o cov matrix) to the plane X=xk in the field b[].
+ //
+ // X [cm] is in the "tracking coordinate system" of this track.
+ // b[]={Bx,By,Bz} [kG] is in the Global coordidate system.
+ //----------------------------------------------------------------
+
+ Double_t dx=xk-fX;
+ if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
+ if (TMath::Abs(fP[4])<=kAlmost0) return kFALSE;
+ // Do not propagate tracks outside the ALICE detector
+ if (TMath::Abs(dx)>1e5 ||
+ TMath::Abs(GetY())>1e5 ||
+ TMath::Abs(GetZ())>1e5) {
+ AliWarning(Form("Anomalous track, target X:%f",xk));
+ Print();
+ return kFALSE;
+ }
+
+ Double_t crv=GetC(b[2]);
+ if (TMath::Abs(b[2]) < 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;
+ //
+ Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
+ //
+ // Appoximate step length
+ double dy2dx = (f1+f2)/(r1+r2);
+ Double_t step = (TMath::Abs(x2r)<0.05) ? dx*TMath::Abs(r2 + f2*dy2dx) // chord
+ : 2.*TMath::ASin(0.5*dx*TMath::Sqrt(1.+dy2dx*dy2dx)*crv)/crv; // arc
+ step *= TMath::Sqrt(1.+ GetTgl()*GetTgl());
+
+ // Get the track's (x,y,z) and (px,py,pz) in the Global System
+ Double_t r[3]; GetXYZ(r);
+ Double_t p[3]; GetPxPyPz(p);
+ Double_t pp=GetP();
+ p[0] /= pp;
+ p[1] /= pp;
+ p[2] /= pp;
+
+ // Rotate to the system where Bx=By=0.
+ Double_t bt=TMath::Sqrt(b[0]*b[0] + b[1]*b[1]);
+ Double_t cosphi=1., sinphi=0.;
+ if (bt > kAlmost0) {cosphi=b[0]/bt; sinphi=b[1]/bt;}
+ Double_t bb=TMath::Sqrt(b[0]*b[0] + b[1]*b[1] + b[2]*b[2]);
+ Double_t costet=1., sintet=0.;
+ if (bb > kAlmost0) {costet=b[2]/bb; sintet=bt/bb;}
+ Double_t vect[7];
+
+ vect[0] = costet*cosphi*r[0] + costet*sinphi*r[1] - sintet*r[2];
+ vect[1] = -sinphi*r[0] + cosphi*r[1];
+ vect[2] = sintet*cosphi*r[0] + sintet*sinphi*r[1] + costet*r[2];
+
+ vect[3] = costet*cosphi*p[0] + costet*sinphi*p[1] - sintet*p[2];
+ vect[4] = -sinphi*p[0] + cosphi*p[1];
+ vect[5] = sintet*cosphi*p[0] + sintet*sinphi*p[1] + costet*p[2];
+
+ vect[6] = pp;
+
+ // Do the helix step
+ g3helx3(GetSign()*bb,step,vect);
+
+ // Rotate back to the Global System
+ r[0] = cosphi*costet*vect[0] - sinphi*vect[1] + cosphi*sintet*vect[2];
+ r[1] = sinphi*costet*vect[0] + cosphi*vect[1] + sinphi*sintet*vect[2];
+ r[2] = -sintet*vect[0] + costet*vect[2];
+
+ p[0] = cosphi*costet*vect[3] - sinphi*vect[4] + cosphi*sintet*vect[5];
+ p[1] = sinphi*costet*vect[3] + cosphi*vect[4] + sinphi*sintet*vect[5];
+ p[2] = -sintet*vect[3] + costet*vect[5];
+
+ // Rotate back to the Tracking System
+ Double_t cosalp = TMath::Cos(fAlpha);
+ Double_t sinalp =-TMath::Sin(fAlpha);
+
+ Double_t
+ t = cosalp*r[0] - sinalp*r[1];
+ r[1] = sinalp*r[0] + cosalp*r[1];
+ r[0] = t;
+
+ t = cosalp*p[0] - sinalp*p[1];
+ p[1] = sinalp*p[0] + cosalp*p[1];
+ p[0] = t;
+
+ // Do the final correcting step to the target plane (linear approximation)
+ Double_t x=r[0], y=r[1], z=r[2];
+ if (TMath::Abs(dx) > kAlmost0) {
+ if (TMath::Abs(p[0]) < kAlmost0) return kFALSE;
+ dx = xk - r[0];
+ x += dx;
+ y += p[1]/p[0]*dx;
+ z += p[2]/p[0]*dx;
+ }
+
+
+ // Calculate the track parameters
+ t=TMath::Sqrt(p[0]*p[0] + p[1]*p[1]);
+ fX = x;
+ fP[0] = y;
+ fP[1] = z;
+ fP[2] = p[1]/t;
+ fP[3] = p[2]/t;
+ fP[4] = GetSign()/(t*pp);
+
+ return kTRUE;
+}
+
+
Bool_t AliExternalTrackParam::Translate(Double_t *vTrasl,Double_t *covV){
//
//Translation: in the event mixing, the tracks can be shifted
// In case the diagonal element is bigger than the maximal allowed value, it is set to
// the limit and the off-diagonal elements that correspond to it are set to zero.
- fC[0] = TMath::Abs(fC[0]);
- if (fC[0]>kC0max) {
- fC[0] = kC0max;
- fC[1] = 0;
- fC[3] = 0;
- fC[6] = 0;
- fC[10] = 0;
- }
- fC[2] = TMath::Abs(fC[2]);
- if (fC[2]>kC2max) {
- fC[2] = kC2max;
- fC[1] = 0;
- fC[4] = 0;
- fC[7] = 0;
- fC[11] = 0;
- }
- fC[5] = TMath::Abs(fC[5]);
- if (fC[5]>kC5max) {
- fC[5] = kC5max;
- fC[3] = 0;
- fC[4] = 0;
- fC[8] = 0;
- fC[12] = 0;
- }
- fC[9] = TMath::Abs(fC[9]);
- if (fC[9]>kC9max) {
- fC[9] = kC9max;
- fC[6] = 0;
- fC[7] = 0;
- fC[8] = 0;
- fC[13] = 0;
- }
- fC[14] = TMath::Abs(fC[14]);
- if (fC[14]>kC14max) {
- fC[14] = kC14max;
- fC[10] = 0;
- fC[11] = 0;
- fC[12] = 0;
- fC[13] = 0;
- }
-
+ fC[0] = TMath::Abs(fC[0]);
+ if (fC[0]>kC0max) {
+ double scl = TMath::Sqrt(kC0max/fC[0]);
+ fC[0] = kC0max;
+ fC[1] *= scl;
+ fC[3] *= scl;
+ fC[6] *= scl;
+ fC[10] *= scl;
+ }
+ fC[2] = TMath::Abs(fC[2]);
+ if (fC[2]>kC2max) {
+ double scl = TMath::Sqrt(kC2max/fC[2]);
+ fC[2] = kC2max;
+ fC[1] *= scl;
+ fC[4] *= scl;
+ fC[7] *= scl;
+ fC[11] *= scl;
+ }
+ fC[5] = TMath::Abs(fC[5]);
+ if (fC[5]>kC5max) {
+ double scl = TMath::Sqrt(kC5max/fC[5]);
+ fC[5] = kC5max;
+ fC[3] *= scl;
+ fC[4] *= scl;
+ fC[8] *= scl;
+ fC[12] *= scl;
+ }
+ fC[9] = TMath::Abs(fC[9]);
+ if (fC[9]>kC9max) {
+ double scl = TMath::Sqrt(kC9max/fC[9]);
+ fC[9] = kC9max;
+ fC[6] *= scl;
+ fC[7] *= scl;
+ fC[8] *= scl;
+ fC[13] *= scl;
+ }
+ fC[14] = TMath::Abs(fC[14]);
+ if (fC[14]>kC14max) {
+ double scl = TMath::Sqrt(kC14max/fC[14]);
+ fC[14] = kC14max;
+ fC[10] *= scl;
+ fC[11] *= scl;
+ fC[12] *= scl;
+ fC[13] *= scl;
+ }
+
// The part below is used for tests and normally is commented out
// TMatrixDSym m(5);
// TVectorD eig(5);
// eig.Print();
// }
}
+
+Bool_t AliExternalTrackParam::ConstrainToVertex(const AliVVertex* vtx, Double_t b[3])
+{
+ // Constrain TPC inner params constrained
+ //
+ if (!vtx)
+ return kFALSE;
+
+ Double_t dz[2], cov[3];
+ if (!PropagateToDCABxByBz(vtx, b, 3, dz, cov))
+ return kFALSE;
+
+ Double_t covar[6];
+ vtx->GetCovarianceMatrix(covar);
+
+ Double_t p[2]= { fP[0] - dz[0], fP[1] - dz[1] };
+ Double_t c[3]= { covar[2], 0., covar[5] };
+
+ Double_t chi2C = GetPredictedChi2(p,c);
+ if (chi2C>kVeryBig)
+ return kFALSE;
+
+ if (!Update(p,c))
+ return kFALSE;
+
+ return kTRUE;
+}
+
+//___________________________________________________________________________________________
+Bool_t AliExternalTrackParam::GetXatLabR(Double_t r,Double_t &x, Double_t bz, Int_t dir) const
+{
+ // Get local X of the track position estimated at the radius lab radius r.
+ // The track curvature is accounted exactly
+ //
+ // The flag "dir" can be used to remove the ambiguity of which intersection to take (out of 2 possible)
+ // 0 - take the intersection closest to the current track position
+ // >0 - go along the track (increasing fX)
+ // <0 - go backward (decreasing fX)
+ //
+ const Double_t &fy=fP[0], &sn = fP[2];
+ const double kEps = 1.e-6;
+ //
+ double crv = GetC(bz);
+ if (TMath::Abs(crv)>kAlmost0) { // helix
+ // get center of the track circle
+ double tR = 1./crv; // track radius (for the moment signed)
+ double cs = TMath::Sqrt((1-sn)*(1+sn));
+ double x0 = fX - sn*tR;
+ double y0 = fy + cs*tR;
+ double r0 = TMath::Sqrt(x0*x0+y0*y0);
+ // printf("Xc:%+e Yc:%+e tR:%e r0:%e\n",x0,y0,tR,r0);
+ //
+ if (r0<=kAlmost0) return kFALSE; // the track is concentric to circle
+ tR = TMath::Abs(tR);
+ double tR2r0=1.,g=0,tmp=0;
+ if (TMath::Abs(tR-r0)>kEps) {
+ tR2r0 = tR/r0;
+ g = 0.5*(r*r/(r0*tR) - tR2r0 - 1./tR2r0);
+ tmp = 1.+g*tR2r0;
+ }
+ else {
+ tR2r0 = 1.0;
+ g = 0.5*r*r/(r0*tR) - 1;
+ tmp = 0.5*r*r/(r0*r0);
+ }
+ double det = (1.-g)*(1.+g);
+ if (det<0) return kFALSE; // does not reach raduis r
+ det = TMath::Sqrt(det);
+ //
+ // the intersection happens in 2 points: {x0+tR*C,y0+tR*S}
+ // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det
+ // where s0 and c0 make direction for the circle center (=x0/r0 and y0/r0)
+ //
+ x = x0*tmp;
+ double y = y0*tmp;
+ if (TMath::Abs(y0)>kAlmost0) { // when y0==0 the x,y is unique
+ double dfx = tR2r0*TMath::Abs(y0)*det;
+ double dfy = tR2r0*x0*TMath::Sign(det,y0);
+ if (dir==0) { // chose the one which corresponds to smallest step
+ double delta = (x-fX)*dfx-(y-fy)*dfy; // the choice of + in C will lead to smaller step if delta<0
+ if (delta<0) x += dfx;
+ else x -= dfx;
+ }
+ else if (dir>0) { // along track direction: x must be > fX
+ x -= dfx; // try the smallest step (dfx is positive)
+ double dfeps = fX-x; // handle special case of very small step
+ if (dfeps<-kEps) return kTRUE;
+ if (TMath::Abs(dfeps)<kEps && // are we already in right r?
+ TMath::Abs(fX*fX+fy*fy - r*r)<kEps) return fX;
+ x += dfx+dfx;
+ if (x-fX>0) return kTRUE;
+ if (x-fX<-kEps) return kFALSE;
+ x = fX; // don't move
+ }
+ else { // backward: x must be < fX
+ x += dfx; // try the smallest step (dfx is positive)
+ double dfeps = x-fX; // handle special case of very small step
+ if (dfeps<-kEps) return kTRUE;
+ if (TMath::Abs(dfeps)<kEps && // are we already in right r?
+ TMath::Abs(fX*fX+fy*fy - r*r)<kEps) return fX;
+ x-=dfx+dfx;
+ if (x-fX<0) return kTRUE;
+ if (x-fX>kEps) return kFALSE;
+ x = fX; // don't move
+ }
+ }
+ else { // special case: track touching the circle just in 1 point
+ if ( (dir>0&&x<fX) || (dir<0&&x>fX) ) return kFALSE;
+ }
+ }
+ else { // this is a straight track
+ if (TMath::Abs(sn)>=kAlmost1) { // || to Y axis
+ double det = (r-fX)*(r+fX);
+ if (det<0) return kFALSE; // does not reach raduis r
+ x = fX;
+ if (dir==0) return kTRUE;
+ det = TMath::Sqrt(det);
+ if (dir>0) { // along the track direction
+ if (sn>0) {if (fy>det) return kFALSE;} // track is along Y axis and above the circle
+ else {if (fy<-det) return kFALSE;} // track is against Y axis amd belo the circle
+ }
+ else if(dir>0) { // agains track direction
+ if (sn>0) {if (fy<-det) return kFALSE;} // track is along Y axis
+ else if (fy>det) return kFALSE; // track is against Y axis
+ }
+ }
+ else if (TMath::Abs(sn)<=kAlmost0) { // || to X axis
+ double det = (r-fy)*(r+fy);
+ if (det<0) return kFALSE; // does not reach raduis r
+ det = TMath::Sqrt(det);
+ if (!dir) {
+ x = fX>0 ? det : -det; // choose the solution requiring the smalest step
+ return kTRUE;
+ }
+ else if (dir>0) { // along the track direction
+ if (fX > det) return kFALSE; // current point is in on the right from the circle
+ else if (fX <-det) x = -det; // on the left
+ else x = det; // within the circle
+ }
+ else { // against the track direction
+ if (fX <-det) return kFALSE;
+ else if (fX > det) x = det;
+ else x = -det;
+ }
+ }
+ else { // general case of straight line
+ double cs = TMath::Sqrt((1-sn)*(1+sn));
+ double xsyc = fX*sn-fy*cs;
+ double det = (r-xsyc)*(r+xsyc);
+ if (det<0) return kFALSE; // does not reach raduis r
+ det = TMath::Sqrt(det);
+ double xcys = fX*cs+fy*sn;
+ double t = -xcys;
+ if (dir==0) t += t>0 ? -det:det; // chose the solution requiring the smalest step
+ else if (dir>0) { // go in increasing fX direction. ( t+-det > 0)
+ if (t>=-det) t += -det; // take minimal step giving t>0
+ else return kFALSE; // both solutions have negative t
+ }
+ else { // go in increasing fX direction. (t+-det < 0)
+ if (t<det) t -= det; // take minimal step giving t<0
+ else return kFALSE; // both solutions have positive t
+ }
+ x = fX + cs*t;
+ }
+ }
+ //
+ return kTRUE;
+}
+//_________________________________________________________
+Bool_t AliExternalTrackParam::GetXYZatR(Double_t xr,Double_t bz, Double_t *xyz, Double_t* alpSect) const
+{
+ // This method has 3 modes of behaviour
+ // 1) xyz[3] array is provided but alpSect pointer is 0: calculate the position of track intersection
+ // with circle of radius xr and fill it in xyz array
+ // 2) alpSect pointer is provided: find alpha of the sector where the track reaches local coordinate xr
+ // Note that in this case xr is NOT the radius but the local coordinate.
+ // If the xyz array is provided, it will be filled by track lab coordinates at local X in this sector
+ // 3) Neither alpSect nor xyz pointers are provided: just check if the track reaches radius xr
+ //
+ //
+ double crv = GetC(bz);
+ if ( (TMath::Abs(bz))<kAlmost0Field ) crv=0.;
+ const double &fy = fP[0];
+ const double &fz = fP[1];
+ const double &sn = fP[2];
+ const double &tgl = fP[3];
+ //
+ // general circle parameterization:
+ // x = (r0+tR)cos(phi0) - tR cos(t+phi0)
+ // y = (r0+tR)sin(phi0) - tR sin(t+phi0)
+ // where qb is the sign of the curvature, tR is the track's signed radius and r0
+ // is the DCA of helix to origin
+ //
+ double tR = 1./crv; // track radius signed
+ double cs = TMath::Sqrt((1-sn)*(1+sn));
+ double x0 = fX - sn*tR; // helix center coordinates
+ double y0 = fy + cs*tR;
+ double phi0 = TMath::ATan2(y0,x0); // angle of PCA wrt to the origin
+ if (tR<0) phi0 += TMath::Pi();
+ if (phi0 > TMath::Pi()) phi0 -= 2.*TMath::Pi();
+ else if (phi0 <-TMath::Pi()) phi0 += 2.*TMath::Pi();
+ double cs0 = TMath::Cos(phi0);
+ double sn0 = TMath::Sin(phi0);
+ double r0 = x0*cs0 + y0*sn0 - tR; // DCA to origin
+ double r2R = 1.+r0/tR;
+ //
+ //
+ if (r2R<kAlmost0) return kFALSE; // helix is centered at the origin, no specific intersection with other concetric circle
+ if (!xyz && !alpSect) return kTRUE;
+ double xr2R = xr/tR;
+ double r2Ri = 1./r2R;
+ // the intersection cos(t) = [1 + (r0/tR+1)^2 - (r0/tR)^2]/[2(1+r0/tR)]
+ double cosT = 0.5*(r2R + (1-xr2R*xr2R)*r2Ri);
+ if ( TMath::Abs(cosT)>kAlmost1 ) {
+ // printf("Does not reach : %f %f\n",r0,tR);
+ return kFALSE; // track does not reach the radius xr
+ }
+ //
+ double t = TMath::ACos(cosT);
+ if (tR<0) t = -t;
+ // intersection point
+ double xyzi[3];
+ xyzi[0] = x0 - tR*TMath::Cos(t+phi0);
+ xyzi[1] = y0 - tR*TMath::Sin(t+phi0);
+ if (xyz) { // if postition is requested, then z is needed:
+ double t0 = TMath::ATan2(cs,-sn) - phi0;
+ double z0 = fz - t0*tR*tgl;
+ xyzi[2] = z0 + tR*t*tgl;
+ }
+ else xyzi[2] = 0;
+ //
+ Local2GlobalPosition(xyzi,fAlpha);
+ //
+ if (xyz) {
+ xyz[0] = xyzi[0];
+ xyz[1] = xyzi[1];
+ xyz[2] = xyzi[2];
+ }
+ //
+ if (alpSect) {
+ double &alp = *alpSect;
+ // determine the sector of crossing
+ double phiPos = TMath::Pi()+TMath::ATan2(-xyzi[1],-xyzi[0]);
+ int sect = ((Int_t)(phiPos*TMath::RadToDeg()))/20;
+ alp = TMath::DegToRad()*(20*sect+10);
+ double x2r,f1,f2,r1,r2,dx,dy2dx,yloc=0, ylocMax = xr*TMath::Tan(TMath::Pi()/18); // min max Y within sector at given X
+ //
+ while(1) {
+ Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
+ if ((cs*ca+sn*sa)<0) {
+ AliDebug(1,Form("Rotation to target sector impossible: local cos(phi) would become %.2f",cs*ca+sn*sa));
+ return kFALSE;
+ }
+ //
+ f1 = sn*ca - cs*sa;
+ if (TMath::Abs(f1) >= kAlmost1) {
+ AliDebug(1,Form("Rotation to target sector impossible: local sin(phi) would become %.2f",f1));
+ return kFALSE;
+ }
+ //
+ double tmpX = fX*ca + fy*sa;
+ double tmpY = -fX*sa + fy*ca;
+ //
+ // estimate Y at X=xr
+ dx=xr-tmpX;
+ x2r = crv*dx;
+ f2=f1 + x2r;
+ if (TMath::Abs(f2) >= kAlmost1) {
+ AliDebug(1,Form("Propagation in target sector failed ! %.10e",f2));
+ return kFALSE;
+ }
+ r1 = TMath::Sqrt((1.-f1)*(1.+f1));
+ r2 = TMath::Sqrt((1.-f2)*(1.+f2));
+ dy2dx = (f1+f2)/(r1+r2);
+ yloc = tmpY + dx*dy2dx;
+ if (yloc>ylocMax) {alp += 2*TMath::Pi()/18; sect++;}
+ else if (yloc<-ylocMax) {alp -= 2*TMath::Pi()/18; sect--;}
+ else break;
+ if (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
+ else if (alp < -TMath::Pi()) alp += 2*TMath::Pi();
+ // if (sect>=18) sect = 0;
+ // if (sect<=0) sect = 17;
+ }
+ //
+ // if alpha was requested, then recalculate the position at intersection in sector
+ if (xyz) {
+ xyz[0] = xr;
+ xyz[1] = yloc;
+ if (TMath::Abs(x2r)<0.05) xyz[2] = fz + dx*(r2 + f2*dy2dx)*tgl;
+ 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
+ xyz[2] = fz + rot/crv*tgl;
+ }
+ Local2GlobalPosition(xyz,alp);
+ }
+ }
+ return kTRUE;
+ //
+}
+
+
+Double_t AliExternalTrackParam::GetParameterAtRadius(Double_t r, Double_t bz, Int_t parType) const
+{
+ //
+ // Get track parameters at the radius of interest.
+ // Given function is aimed to be used to interactivelly (tree->Draw())
+ // access track properties at different radii
+ //
+ // TO BE USED WITH SPECICAL CARE -
+ // it is aimed to be used for rough calculation as constant field and
+ // no correction for material is used
+ //
+ // r - radius of interest
+ // bz - magentic field
+ // retun values dependens on parType:
+ // parType = 0 -gx
+ // parType = 1 -gy
+ // parType = 2 -gz
+ //
+ // parType = 3 -pgx
+ // parType = 4 -pgy
+ // parType = 5 -pgz
+ //
+ // parType = 6 - r
+ // parType = 7 - global position phi
+ // parType = 8 - global direction phi
+ // parType = 9 - direction phi- positionphi
+ if (parType<0) {
+ parType=-1;
+ return 0;
+ }
+ Double_t xyz[3];
+ Double_t pxyz[3];
+ Double_t localX=0;
+ Bool_t res = GetXatLabR(r,localX,bz,1);
+ if (!res) {
+ parType=-1;
+ return 0;
+ }
+ //
+ // position parameters
+ //
+ GetXYZAt(localX,bz,xyz);
+ if (parType<3) {
+ return xyz[parType];
+ }
+
+ if (parType==6) return TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
+ if (parType==7) return TMath::ATan2(xyz[1],xyz[0]);
+ //
+ // momenta parameters
+ //
+ GetPxPyPzAt(localX,bz,pxyz);
+ if (parType==8) return TMath::ATan2(pxyz[1],pxyz[0]);
+ if (parType==9) {
+ Double_t diff = TMath::ATan2(pxyz[1],pxyz[0])-TMath::ATan2(xyz[1],xyz[0]);
+ if (diff>TMath::Pi()) diff-=TMath::TwoPi();
+ if (diff<-TMath::Pi()) diff+=TMath::TwoPi();
+ return diff;
+ }
+ if (parType>=3&&parType<6) {
+ return pxyz[parType%3];
+ }
+ return 0;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff