]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliExternalTrackParam.cxx
Bug fix
[u/mrichter/AliRoot.git] / STEER / AliExternalTrackParam.cxx
index 6861b0b1a64f823dc6b1814af8c34d93b40f9b48..2f4eeda5cc61390424b90b689f5fb946d686ff7e 100644 (file)
@@ -178,8 +178,8 @@ Double_t AliExternalTrackParam::GetLinearD(Double_t xv,Double_t yv) const {
   return -d;
 }
 
-Bool_t AliExternalTrackParam::
-CorrectForMaterial(Double_t d,  Double_t x0, Double_t mass) {
+Bool_t AliExternalTrackParam::CorrectForMaterial
+(Double_t d,  Double_t x0, Double_t mass, Double_t (*Bethe)(Double_t)) {
   //------------------------------------------------------------------
   // This function corrects the track parameters for the crossed material
   // "d"    - the thickness (fraction of the radiation length)
@@ -213,16 +213,32 @@ CorrectForMaterial(Double_t d,  Double_t x0, Double_t mass) {
   //Energy losses************************
   if (x0!=0. && beta2<1) {
      d*=x0;
-     Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
-     if (beta2/(1-beta2)>3.5*3.5)
-       dE=0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2)*d;
-
-     fP4*=(1.- TMath::Sqrt(p2 + mass*mass)/p2*dE);
+     Double_t dE=Bethe(beta2)*d;
+     Double_t e=TMath::Sqrt(p2 + mass*mass);
+     fP4*=(1.- e/p2*dE);
+
+     // Approximate energy loss fluctuation (M.Ivanov)
+     const Double_t cnst=0.07; // To be tuned.  
+     Double_t sigmadE=cnst*TMath::Sqrt(TMath::Abs(dE)); 
+     fC44+=((sigmadE*e/p2*fP4)*(sigmadE*e/p2*fP4)); 
   }
 
   return kTRUE;
 }
 
+Double_t ApproximateBetheBloch(Double_t beta2) {
+  //------------------------------------------------------------------
+  // This is an approximation of the Bethe-Bloch formula with 
+  // the density effect taken into account at beta*gamma > 3.5
+  // (the approximation is reasonable only for solid materials) 
+  //------------------------------------------------------------------
+  if (beta2/(1-beta2)>3.5*3.5)
+     return 0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2);
+
+  return 0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2);
+}
+
 Bool_t AliExternalTrackParam::Rotate(Double_t alpha) {
   //------------------------------------------------------------------
   // Transform this track to the local coord. system rotated
@@ -352,6 +368,76 @@ Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t b) {
   return kTRUE;
 }
 
+void AliExternalTrackParam::Propagate(Double_t len, Double_t x[3],
+Double_t p[3], Double_t bz) const {
+  //+++++++++++++++++++++++++++++++++++++++++    
+  // Origin: K. Shileev (Kirill.Shileev@cern.ch)
+  // Extrapolate track along simple helix in magnetic field
+  // Arguments: len -distance alogn helix, [cm]
+  //            bz  - mag field, [kGaus]   
+  // Returns: x and p contain extrapolated positon and momentum  
+  // The momentum returned for straight-line tracks is meaningless !
+  //+++++++++++++++++++++++++++++++++++++++++    
+  GetXYZ(x);
+    
+  if (TMath::Abs(Get1Pt()) < kAlmost0){ //straight-line tracks
+     Double_t unit[3]; GetDirection(unit);
+     x[0]+=unit[0]*len;   
+     x[1]+=unit[1]*len;   
+     x[2]+=unit[2]*len;
+
+     p[0]=unit[0]/kAlmost0;   
+     p[1]=unit[1]/kAlmost0;   
+     p[2]=unit[2]/kAlmost0;   
+  } else {
+     GetPxPyPz(p);
+     Double_t pp=GetP();
+     Double_t a = -kB2C*bz*GetSign();
+     Double_t rho = a/pp;
+     x[0] += p[0]*TMath::Sin(rho*len)/a - p[1]*(1-TMath::Cos(rho*len))/a;
+     x[1] += p[1]*TMath::Sin(rho*len)/a + p[0]*(1-TMath::Cos(rho*len))/a;
+     x[2] += p[2]*len/pp;
+
+     Double_t p0=p[0];
+     p[0] = p0  *TMath::Cos(rho*len) - p[1]*TMath::Sin(rho*len);
+     p[1] = p[1]*TMath::Cos(rho*len) + p0  *TMath::Sin(rho*len);
+  }
+}
+
+Bool_t AliExternalTrackParam::Intersect(Double_t pnt[3], Double_t norm[3],
+Double_t bz) const {
+  //+++++++++++++++++++++++++++++++++++++++++    
+  // Origin: K. Shileev (Kirill.Shileev@cern.ch)
+  // Finds point of intersection (if exists) of the helix with the plane. 
+  // Stores result in fX and fP.   
+  // Arguments: planePoint,planeNorm - the plane defined by any plane's point 
+  // and vector, normal to the plane
+  // Returns: kTrue if helix intersects the plane, kFALSE otherwise.
+  //+++++++++++++++++++++++++++++++++++++++++    
+  Double_t x0[3]; GetXYZ(x0); //get track position in MARS
+  
+  //estimates initial helix length up to plane
+  Double_t s=
+    (pnt[0]-x0[0])*norm[0] + (pnt[1]-x0[1])*norm[1] + (pnt[2]-x0[2])*norm[2];
+  Double_t dist=99999,distPrev=dist;
+  Double_t x[3],p[3]; 
+  while(TMath::Abs(dist)>0.00001){
+    //calculates helix at the distance s from x0 ALONG the helix
+    Propagate(s,x,p,bz);
+
+    //distance between current helix position and plane
+    dist=(x[0]-pnt[0])*norm[0]+(x[1]-pnt[1])*norm[1]+(x[2]-pnt[2])*norm[2];
+
+    if(TMath::Abs(dist) >= TMath::Abs(distPrev)) {return kFALSE;}
+    distPrev=dist;
+    s-=dist;
+  }
+  //on exit pnt is intersection point,norm is track vector at that point, 
+  //all in MARS
+  for (Int_t i=0; i<3; i++) {pnt[i]=x[i]; norm[i]=p[i];}
+  return kTRUE;
+}
+
 Double_t 
 AliExternalTrackParam::GetPredictedChi2(Double_t p[2],Double_t cov[3]) const {
   //----------------------------------------------------------------
@@ -677,6 +763,20 @@ Bool_t Local2GlobalPosition(Double_t r[3],Double_t alpha) {
   return kTRUE;
 }
 
+void AliExternalTrackParam::GetDirection(Double_t d[3]) const {
+  //----------------------------------------------------------------
+  // This function returns a unit vector along the track direction
+  // in the global coordinate system.
+  //----------------------------------------------------------------
+  Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
+  Double_t snp=fP[2];
+  Double_t csp =TMath::Sqrt((1.- snp)*(1.+snp));
+  Double_t norm=TMath::Sqrt(1.+ fP[3]*fP[3]);
+  d[0]=(csp*cs - snp*sn)/norm; 
+  d[1]=(snp*cs + csp*sn)/norm; 
+  d[2]=fP[3]/norm;
+}
+
 Bool_t AliExternalTrackParam::GetPxPyPz(Double_t *p) const {
   //---------------------------------------------------------------------
   // This function returns the global track momentum components
@@ -717,7 +817,7 @@ Bool_t AliExternalTrackParam::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
   }
   Double_t pt=1./TMath::Abs(fP[4]);
   Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
-  Double_t r=TMath::Sqrt(1-fP[2]*fP[2]);
+  Double_t r=TMath::Sqrt((1.-fP[2])*(1.+fP[2]));
 
   Double_t m00=-sn, m10=cs;
   Double_t m23=-pt*(sn + fP[2]*cs/r), m43=-pt*pt*(r*cs - fP[2]*sn);