Set(xyz,pxpypz,cv,sign);
}
+/*
//_____________________________________________________________________________
void AliExternalTrackParam::Set(Double_t xyz[3],Double_t pxpypz[3],
Double_t cv[21],Short_t sign)
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;
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() {
// "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) {
Double_t bg=GetP()/mass;
if (mass<0) {
if (mass<-990) {
- AliDebug(2,Form("Mass %f corresponds to unknown PID particle = %e",mass));
+ 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);
}
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;
+ //
+ fP1 += fP3/crv*TMath::ASin(r1*f2 - r2*f1); // more economic version from Yura.
}
//f = F - 1
//---------------------------------------------------------------------
Double_t dx=x-fX;
if(TMath::Abs(dx)<=kAlmost0) {z=fP[1]; return kTRUE;}
+
Double_t crv=GetC(b);
Double_t x2r = crv*dx;
Double_t f1=fP[2], f2=f1 + x2r;
//---------------------------------------------------------------------
Double_t dx=x-fX;
if(TMath::Abs(dx)<=kAlmost0) return GetXYZ(r);
+
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 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);
}
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
// This function forces the diagonal elements of the covariance matrix to be positive.
// 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) {
double scl = TMath::Sqrt(kC0max/fC[0]);
//
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;
+ //
+}