From: shahoian Date: Fri, 24 Jan 2014 17:29:25 +0000 (+0100) Subject: Add method GetXYZatR for fast estimate of track intersection with given X or R X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=c8253dd95447cc9624b1761bef3603d7171d8759 Add method GetXYZatR for fast estimate of track intersection with given X or R --- diff --git a/STEER/AOD/AliAODTrack.cxx b/STEER/AOD/AliAODTrack.cxx index 269d908e6ef..a45e590efce 100644 --- a/STEER/AOD/AliAODTrack.cxx +++ b/STEER/AOD/AliAODTrack.cxx @@ -943,6 +943,164 @@ Bool_t AliAODTrack::GetXYZAt(Double_t x, Double_t b, Double_t *r) const return Local2GlobalPosition(r,alpha); } +//_____________________________________________________________________________ +Bool_t AliAODTrack::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_t alpha=0.0; + Double_t radPos2 = fPosition[0]*fPosition[0]+fPosition[1]*fPosition[1]; + Double_t radMax = 45.; // approximately ITS outer radius + if (radPos2 < radMax*radMax) { // inside the ITS + alpha = fMomentum[1]; //TMath::ATan2(fMomentum[1],fMomentum[0]); // fMom is pt,phi,theta! + } else { // outside the ITS + Float_t phiPos = TMath::Pi()+TMath::ATan2(-fPosition[1], -fPosition[0]); + alpha = + TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10); + } + // + // Get the vertex of origin and the momentum + TVector3 ver(fPosition[0],fPosition[1],fPosition[2]); + TVector3 mom(Px(),Py(),Pz()); + // + // Rotate to the local coordinate system + ver.RotateZ(-alpha); + mom.RotateZ(-alpha); + // + Double_t fx = ver.X(); + Double_t fy = ver.Y(); + Double_t fz = ver.Z(); + Double_t sn = TMath::Sin(mom.Phi()); + Double_t tgl = mom.Pz()/mom.Pt(); + Double_t crv = TMath::Sign(1/mom.Pt(),(Double_t)fCharge)*bz*kB2C; + // + if ( (TMath::Abs(bz)) 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 (r2RkAlmost1 ) { + // 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,alpha); + // + 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-alpha), sa=TMath::Sin(alp-alpha); + 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<