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() {
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
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