]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ITS/AliITStrackMI.cxx
Fixing bug #59311
[u/mrichter/AliRoot.git] / ITS / AliITStrackMI.cxx
index ecbdd6e74b5b8443207371757e808b44394a01f2..df83813c5f6a691a7cca169d77b68573beacac3b 100644 (file)
 //     dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
 //-------------------------------------------------------------------------
 
+/* $Id$ */
+
 #include <TMatrixD.h>
 
 #include <TMath.h>
 
 #include "AliCluster.h"
 #include "AliESDtrack.h"
+#include "AliITSgeomTGeo.h"
 #include "AliITStrackMI.h"
 
 ClassImp(AliITStrackMI)
@@ -34,77 +37,72 @@ const Int_t kWARN=5;
 
 //____________________________________________________________________________
 AliITStrackMI::AliITStrackMI():AliITStrackV2(),
-  fNUsed(0),
-  fNSkipped(0),
-  fNDeadZone(0),
-  fDeadZoneProbability(0),                            
-  fReconstructed(kFALSE),
-  fConstrain(kFALSE)
+fNUsed(0),
+fNSkipped(0),
+fNDeadZone(0),                        
+fReconstructed(kFALSE),
+fExpQ(40),
+fChi22(0),
+fdEdxMismatch(0),
+fConstrain(kFALSE),
+fGoldV0(kFALSE)
 {
-    for(Int_t i=0; i<kMaxLayer; i++) fClIndex[i]=-1;
-    for(Int_t i=0; i<6; i++) { fNy[i]=0; fNz[i]=0; fNormQ[i]=0; fNormChi2[i]=1000;}
-    for(Int_t i=0; i<12; i++) {fDy[i]=0; fDz[i]=0; fSigmaY[i]=0; fSigmaZ[i]=0; fChi2MIP[i]=0;}
+  //constructor
+    for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fClIndex[i]=-1;
+    for(Int_t i=0; i<6; i++) { fNy[i]=0; fNz[i]=0; fNormQ[i]=0; fNormChi2[i]=1000; fDeadZoneProbability[i]=0;}
+    for(Int_t i=0; i<12; i++) {fDy[i]=0; fDz[i]=0; fSigmaY[i]=0; fSigmaZ[i]=0; fSigmaYZ[i]=0; fChi2MIP[i]=0;}
     fD[0]=0; fD[1]=0;
-    fExpQ=40;
-    fdEdxMismatch=0;
-    fChi22=0;
-    fGoldV0 = kFALSE;
+    fDnorm[0]=0; fDnorm[1]=0;
 }
 
 //____________________________________________________________________________
 AliITStrackMI::AliITStrackMI(AliESDtrack& t,Bool_t c) throw (const Char_t *) :
-AliITStrackV2(t,c) {
+AliITStrackV2(t,c),
+fNUsed(0),
+fNSkipped(0),
+fNDeadZone(0),
+fReconstructed(kFALSE),
+fExpQ(40),
+fChi22(0),
+fdEdxMismatch(0),
+fConstrain(kFALSE),
+fGoldV0(kFALSE) {
   //------------------------------------------------------------------
   // Conversion ESD track -> ITS track.
   // If c==kTRUE, create the ITS track out of the constrained params.
   //------------------------------------------------------------------
-  fNUsed = 0;
-  fReconstructed = kFALSE;
-  fNSkipped =0; 
-  fNDeadZone = 0;
-  fDeadZoneProbability = 0;
-  for(Int_t i=0; i<6; i++) {fClIndex[i]=-1; fNy[i]=0; fNz[i]=0; fNormQ[i]=0; fNormChi2[i]=1000;}
-  for(Int_t i=0; i<12; i++) {fDy[i]=0; fDz[i]=0; fSigmaY[i]=0; fSigmaZ[i]=0;fChi2MIP[i]=0;}
+  for(Int_t i=0; i<6; i++) {fClIndex[i]=-1; fNy[i]=0; fNz[i]=0; fNormQ[i]=0; fNormChi2[i]=1000; fDeadZoneProbability[i]=0;}
+  for(Int_t i=0; i<12; i++) {fDy[i]=0; fDz[i]=0; fSigmaY[i]=0; fSigmaZ[i]=0; fSigmaYZ[i]=0; fChi2MIP[i]=0;}
   fD[0]=0; fD[1]=0;
-  fExpQ=40;
-  fConstrain = kFALSE;
-  fdEdxMismatch=0;
-  fChi22 =0;
-  fGoldV0 = kFALSE;
+  fDnorm[0]=0; fDnorm[1]=0;
   //if (!Invariant()) throw "AliITStrackV2: conversion failed !\n";
 
 }
 
-void AliITStrackMI::UpdateESDtrack(ULong_t flags) {
-  fESDtrack->UpdateTrackParams(this,flags);
-  if (flags == AliESDtrack::kITSin) fESDtrack->SetITSChi2MIP(fChi2MIP);
-}
-
 //____________________________________________________________________________
-AliITStrackMI::AliITStrackMI(const AliITStrackMI& t) : AliITStrackV2(t) {
+AliITStrackMI::AliITStrackMI(const AliITStrackMI& t) : AliITStrackV2(t),
+fNUsed(t.fNUsed),
+fNSkipped(t.fNSkipped),
+fNDeadZone(t.fNDeadZone),
+fReconstructed(t.fReconstructed),
+fExpQ(t.fExpQ),
+fChi22(t.fChi22),
+fdEdxMismatch(t.fdEdxMismatch),
+fConstrain(t.fConstrain),
+fGoldV0(t.fGoldV0) {
   //------------------------------------------------------------------
   //Copy constructor
   //------------------------------------------------------------------
-  fNUsed = t.fNUsed;
-  fReconstructed = t.fReconstructed;
-  fNSkipped = t.fNSkipped;
-  fNDeadZone = t.fNDeadZone;
-  fDeadZoneProbability = t.fDeadZoneProbability;
   fLab = t.fLab;
   fFakeRatio = t.fFakeRatio;
-  fdEdxMismatch = t.fdEdxMismatch;
-  fChi22 = t.fChi22;
-  fGoldV0 = t.fGoldV0;;
 
   fD[0]=t.fD[0]; fD[1]=t.fD[1];
   fDnorm[0] = t.fDnorm[0]; fDnorm[1]=t.fDnorm[1];
-  fExpQ= t.fExpQ;
   for(Int_t i=0; i<6; i++) {
-    fClIndex[i]= t.fClIndex[i]; fNy[i]=t.fNy[i]; fNz[i]=t.fNz[i]; fNormQ[i]=t.fNormQ[i]; fNormChi2[i] = t.fNormChi2[i];
+    fClIndex[i]= t.fClIndex[i]; fNy[i]=t.fNy[i]; fNz[i]=t.fNz[i]; fNormQ[i]=t.fNormQ[i]; fNormChi2[i] = t.fNormChi2[i];  fDeadZoneProbability[i]=t.fDeadZoneProbability[i];
   }
   for(Int_t i=0; i<12; i++) {fDy[i]=t.fDy[i]; fDz[i]=t.fDz[i]; 
-    fSigmaY[i]=t.fSigmaY[i]; fSigmaZ[i]=t.fSigmaZ[i];fChi2MIP[i]=t.fChi2MIP[i];}
-  fConstrain = t.fConstrain;
+    fSigmaY[i]=t.fSigmaY[i]; fSigmaZ[i]=t.fSigmaZ[i]; fSigmaYZ[i]=t.fSigmaYZ[i]; fChi2MIP[i]=t.fChi2MIP[i];}
   //memcpy(fDy,t.fDy,6*sizeof(Float_t));
   //memcpy(fDz,t.fDz,6*sizeof(Float_t));
   //memcpy(fSigmaY,t.fSigmaY,6*sizeof(Float_t));
@@ -128,153 +126,31 @@ Int_t AliITStrackMI::Compare(const TObject *o) const {
 }
 
 
-Double_t AliITStrackMI::GetPredictedChi2MI(Double_t cy, Double_t cz, Double_t cerry, Double_t cerrz) const
+Double_t AliITStrackMI::GetPredictedChi2MI(Double_t cy, Double_t cz, Double_t cerry, Double_t cerrz, Double_t covyz) const
 {
   //-----------------------------------------------------------------
   // This function calculates a predicted chi2 increment.
   //-----------------------------------------------------------------
-  Double_t r00=cerry*cerry, r01=0., r11=cerrz*cerrz;
-  r00+=fC00; r01+=fC10; r11+=fC11;
-  //
-  Double_t det=r00*r11 - r01*r01;
-  if (TMath::Abs(det) < 1.e-30) {
-    Int_t n=GetNumberOfClusters();
-    if (n>kWARN) 
-      Warning("GetPredictedChi2","Singular matrix (%d) !\n",n);
-    return 1e10;
-  }
-  Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
-
-  Double_t dy=cy - fP0, dz=cz - fP1;
-
-  return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
-}
-
-//____________________________________________________________________________
-Int_t AliITStrackMI::CorrectForMaterial(Double_t d, Double_t x0) {
-  //------------------------------------------------------------------
-  //This function corrects the track parameters for crossed material
-  //------------------------------------------------------------------
-  //  Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
-  Double_t p2=(1.+ fP3*fP3)/(Get1Pt()*Get1Pt());
-  Double_t et   = p2 + GetMass()*GetMass();
-  Double_t beta2=p2/et;
-  et = sqrt(et);  
-  d*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
-  //d*=TMath::Sqrt(1.+ fP3*fP3 +fP2*fP2/(1.- fP2*fP2));
-
-  //Multiple scattering******************
-  if (d!=0) {
-     Double_t theta2=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(d);
-     //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*TMath::Abs(d)*9.36*2.33;
-     fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
-     fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
-     fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
-     fC44 += theta2*fP3*fP4*fP3*fP4;
-  }
-
-  //Energy losses************************
-  if (x0!=0.) {
-     d*=x0;
-     //     Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
-     //Double_t dE=0; 
-     Double_t dE = 0.265*0.153e-3*(39.2-55.6*beta2+28.7*beta2*beta2+27.41/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;
-     }
-     else{
-       dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
-       dE+=0.06e-3/(beta2*beta2)*d;
-     }
-     */
-     fP4*=(1.- et/p2*dE);
-     Double_t delta44 = (dE*fP4*et/p2);
-     delta44*=delta44;
-     fC44+= delta44/400.;
-  }
-
-  if (!Invariant()) return 0;
-
-  return 1;
+  Double_t p[2]={cy, cz};
+  Double_t cov[3]={cerry*cerry, covyz, cerrz*cerrz};
+  return AliExternalTrackParam::GetPredictedChi2(p,cov);
 }
 
 //____________________________________________________________________________
-Int_t AliITStrackMI::UpdateMI(Double_t cy, Double_t cz, Double_t cerry, Double_t cerrz, Double_t chi2,UInt_t index) {
+Bool_t AliITStrackMI::UpdateMI(const AliCluster *c, Double_t chi2, Int_t index) {
   //------------------------------------------------------------------
   //This function updates track parameters
   //------------------------------------------------------------------
-  Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
-  Double_t c00=fC00;
-  Double_t c10=fC10, c11=fC11;
-  Double_t c20=fC20, c21=fC21, c22=fC22;
-  Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
-  Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
-
-
-  Double_t r00=cerry*cerry, r01=0., r11=cerrz*cerrz;
-  r00+=fC00; r01+=fC10; r11+=fC11;
-  Double_t det=r00*r11 - r01*r01;
-  Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
-
-  Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
-  Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
-  Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
-  Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
-  Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
-
-  Double_t dy=cy - fP0, dz=cz - fP1;
+  Double_t dy=c->GetY() - GetY(), dz=c->GetZ() - GetZ();
   Int_t layer = (index & 0xf0000000) >> 28;
   fDy[layer] = dy;
   fDz[layer] = dz;
-  fSigmaY[layer] = TMath::Sqrt(cerry*cerry+fC00);
-  fSigmaZ[layer] = TMath::Sqrt(cerrz*cerrz+fC11);
-
-  Double_t sf=fP2 + k20*dy + k21*dz;
-  
-  fP0 += k00*dy + k01*dz;
-  fP1 += k10*dy + k11*dz;
-  fP2  = sf;
-  fP3 += k30*dy + k31*dz;
-  fP4 += k40*dy + k41*dz;
-  
-  Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
-  Double_t c12=fC21, c13=fC31, c14=fC41;
+  fSigmaY[layer] = TMath::Sqrt(c->GetSigmaY2()+GetSigmaY2());
+  fSigmaZ[layer] = TMath::Sqrt(c->GetSigmaZ2()+GetSigmaZ2());
+  fSigmaYZ[layer] = c->GetSigmaYZ()+GetSigmaZY();
 
-  fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
-  fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
-  fC40-=k00*c04+k01*c14; 
 
-  fC11-=k10*c01+k11*fC11;
-  fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
-  fC41-=k10*c04+k11*c14; 
-
-  fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
-  fC42-=k20*c04+k21*c14; 
-
-  fC33-=k30*c03+k31*c13;
-  fC43-=k30*c04+k31*c14; 
-
-  fC44-=k40*c04+k41*c14; 
-
-  if (!Invariant()) {
-     fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
-     fC00=c00;
-     fC10=c10; fC11=c11;
-     fC20=c20; fC21=c21; fC22=c22;
-     fC30=c30; fC31=c31; fC32=c32; fC33=c33;
-     fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
-     return 0;
-  }
-
-  if (chi2<0) return 1;
-  Int_t n=GetNumberOfClusters();
-  fIndex[n]=index;
-  SetNumberOfClusters(n+1);
-  SetChi2(GetChi2()+chi2);
-
-  return 1;
+  return Update(c,chi2,index);
 }
 
 Int_t AliITStrackMI::GetProlongationFast(Double_t alp, Double_t xk,Double_t &y, Double_t &z)
@@ -282,19 +158,19 @@ Int_t AliITStrackMI::GetProlongationFast(Double_t alp, Double_t xk,Double_t &y,
   //-----------------------------------------------------------------------------
   //get fast prolongation 
   //-----------------------------------------------------------------------------
-  Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
-  Double_t cf=TMath::Sqrt(1.- fP2*fP2);  
+  Double_t ca=TMath::Cos(alp-GetAlpha()), sa=TMath::Sin(alp-GetAlpha());
+  Double_t cf=TMath::Sqrt((1.-GetSnp())*(1.+GetSnp()));  
   // **** rotation **********************  
-  y= -fX*sa + fP0*ca;
+  y= -GetX()*sa + GetY()*ca;
   // **** translation ******************  
-  Double_t dx = xk- fX*ca - fP0*sa;
-  Double_t f1=fP2*ca - cf*sa, f2=f1 + fP4*dx;
+  Double_t dx = xk- GetX()*ca - GetY()*sa;
+  Double_t f1=GetSnp()*ca - cf*sa, f2=f1 + GetC()*dx;
   if (TMath::Abs(f2) >= 0.9999) {
     return 0;
   }
-  Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);  
+  Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));  
   y += dx*(f1+f2)/(r1+r2);
-  z  = fP1+dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;  
+  z  = GetZ()+dx*(f1+f2)/(f1*r2 + f2*r1)*GetTgl();  
   return 1;
 }