// ----------------------------------------------------------------------------
// $Id: LCDSmearTrack.cxx,v 1.2 2001/06/23 16:15:45 toshi Exp $
// ----------------------------------------------------------------------------
//
// $Log: LCDSmearTrack.cxx,v $
// Revision 1.2 2001/06/23 16:15:45 toshi
// In method arg. , Char_t* -> const Char_t*
//
// Revision 1.1 2001/05/10 17:56:03 toshi
// Changes C++ name convention from *.C to *.cxx to allow Windows usage.
//
//
// algorithm from K.Riles 8.12.98
// Version 1.0 12-Aug-98 Richard incorporate K.Riles' algorithm.
// Version 1.1 06-Mar-99 Richard get start point from parent MC.
// Version 1.2 26-Jul-99 M.Iwasaki Remove momentum, charge, position inputs
// in Track().
// Version 1.3 07-Apr-00 T.Abe get CA point for track smearing and
// change track parameter definition.
#include "LCDSmearTrack.h"
#include <stdlib.h>
//______________________________________________________________________
//
// SmearTrack
//
// Smears tracks via lookup tables on resolution
ClassImp(LCDSmearTrack)
//______________________________________________________________________
LCDSmearTrack::LCDSmearTrack(){
// constructor: sets up lookups tables, etc.
Init();
}
//______________________________________________________________________
//______________________________________________________________________
LCDSmearTrack::LCDSmearTrack(LCDGetParameters* gp,
const Char_t* smearFileName){
// constructor: sets up lookups tables, etc.
SetUp(gp,smearFileName);
}
//______________________________________________________________________
//______________________________________________________________________
LCDSmearTrack::~LCDSmearTrack(){
if (m_TkLook) {
delete m_TkLook;
m_TkLook=0;
}
}
//______________________________________________________________________
//______________________________________________________________________
void LCDSmearTrack::Init() {
m_parameters = 0;
m_TkLook = 0;
for (Int_t j=0; j < 15 ; j++) {
smear_parNum[j]=0;
}
}
//______________________________________________________________________
//______________________________________________________________________
void LCDSmearTrack::SetUp(LCDGetParameters* gp,
const Char_t* smearFileName){
// constructor: sets up lookups tables, etc.
m_parameters = gp;
FILE* LTBFile = fopen(smearFileName, "r");
m_TkLook = new LCDGetTrackLookups(LTBFile);
fclose(LTBFile);
TObjArray* smearArray = m_TkLook->GetLookupsArray();
Int_t numPar = smearArray->GetEntries();
for (Int_t j=0; j < 15 ; j++) {
smear_parNum[j]=0;
}
for (Int_t i=0; i<numPar ; i++) {
TString* parName = ((LCDLookUp2d*)(smearArray->At(i)))->GetTableName();
if ( parName->CompareTo("(1,1)",TString::kExact) == 0) {
smear_parNum[ 0] =((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(2,1)",TString::kExact) == 0) {
smear_parNum[ 1]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(2,2)",TString::kExact) == 0) {
smear_parNum[ 2]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(3,1)",TString::kExact) == 0) {
smear_parNum[ 3]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(3,2)",TString::kExact) == 0) {
smear_parNum[ 4]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(3,3)",TString::kExact) == 0) {
smear_parNum[ 5]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(4,1)",TString::kExact) == 0) {
smear_parNum[ 6]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(4,2)",TString::kExact) == 0) {
smear_parNum[ 7]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(4,3)",TString::kExact) == 0) {
smear_parNum[ 8]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(4,4)",TString::kExact) == 0) {
smear_parNum[ 9]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(5,1)",TString::kExact) == 0) {
smear_parNum[10]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(5,2)",TString::kExact) == 0) {
smear_parNum[11]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(5,3)",TString::kExact) == 0) {
smear_parNum[12]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(5,4)",TString::kExact) == 0) {
smear_parNum[13]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("(5,5)",TString::kExact) == 0) {
smear_parNum[14]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("r",TString::kExact) == 0) {
smear_parNum[ 0]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("phi",TString::kExact) == 0) {
smear_parNum[ 2]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("1/p_t",TString::kExact) == 0) {
smear_parNum[ 5] =((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("z",TString::kExact) == 0) {
smear_parNum[ 9]=((LCDLookUp2d*)smearArray->At(i));
} else if ( parName->CompareTo("lambda",TString::kExact) == 0) {
smear_parNum[14]=((LCDLookUp2d*)smearArray->At(i));
}
}
}
//______________________________________________________________________
//______________________________________________________________________
void LCDSmearTrack::SmearTrack(LCDMcPart* p, Int_t index, Int_t index0,
LCDGetParameters* gp, LCDTrack* tk) {
// creates track for input McPart if it passes acceptance cuts and
// uses lookup tables for resolutions
// if MCParticle pointer is zero, it belongs to albedo shower particle. Ignore for now.
if (p == 0) return;
// Get original momentum and position from MCParticle
TVector3 pMC((*p->Get4MomentumPtr()).Vect());
TVector3 xMC(*p->GetPositionPtr());
// convert to helical parameters @ Closest Approach Point to Z-axis.
Double_t pTot = pMC.Mag();
Double_t cos_theta = pMC.CosTheta();
Double_t abscth = TMath::Abs(cos_theta);
Double_t qMC = p->GetCharge();
Double_t bfld_z=gp->GetMagneticField();
Double_t tkpar[5];
CalcTrackParameters(qMC,bfld_z,pMC,xMC,&tkpar[0]);
// point at proper set of lookup tables
// Get resolution values from smear parameters
Double_t d_parm[5];
Double_t m_eparm[15];
Double_t chi2;
CalcSmearParameters(&m_eparm[0],&d_parm[0],bfld_z,pTot,abscth,gp,
smear_parNum,&chi2);
Double_t d0 = tkpar[0];
Double_t phi0 = tkpar[1];
Double_t omega= tkpar[2];
Double_t z0 = tkpar[3];
Double_t tanL = tkpar[4];
Double_t d0New = d0 + d_parm[0];
Double_t phi0New = phi0 + d_parm[1] + 4*TMath::Pi();
Int_t iphi0New = (Int_t)(phi0New/2/TMath::Pi());
phi0New -= iphi0New*2*TMath::Pi();
if ( phi0New>TMath::Pi() ) { phi0New -= 2*TMath::Pi(); }
Double_t omegaNew = omega + d_parm[2];
Double_t z0New = z0 + d_parm[3];
Double_t tanLNew = tanL + d_parm[4];
// Save the track parameters
Double_t m_parm[5];
m_parm[0] = d0New;
m_parm[1] = phi0New;
m_parm[2] = omegaNew;
m_parm[3] = z0New;
m_parm[4] = tanLNew;
tk->SetLCDTrack(index,index0,bfld_z,5,chi2,&m_parm[0],&m_eparm[0]);
}
//______________________________________________________________________
//______________________________________________________________________
void
LCDSmearTrack::FindEigen(
TMatrixD& Source,
TMatrixD& RotationMatrix,
TMatrixD& EigenValues) {
// This is code taken directly from BetaTools/BtaSphericity, which is taken
// from Numerical Recipes in C, page 467 (Jacobi). The matrix Source is
// diagonalized, where RotationMatrix is the rotation matrix used to find
// the EigenValues. I've changed no code, but have tried to write it in a
// more readable format.
//
// Note that Source is a local copy of the original matrix passed to it.
#define JACOBI_ROTATION(A,i,j,k,l) g=A(i,j); h=A(k,l); A(i,j) = g-s*(h+g*tau) ; A(k,l) = h + s*(g-h*tau);
Int_t irow;
Int_t icol;
Double_t Threshold;
Int_t mtrxSize = Source.GetNrows();
RotationMatrix.UnitMatrix();
for (irow=0; irow < mtrxSize; irow++) {
EigenValues(irow,0) = Source(irow,irow);
}
TMatrixD Z(mtrxSize,1);
TMatrixD B(mtrxSize,1);
Z.Zero();
B=EigenValues;
Int_t NumberRotations = 0;
for (Int_t Iteration = 0; Iteration < 50; Iteration++) {
Double_t SumOffDiagonal = 0.0;
for (irow = 0; irow < mtrxSize; irow++) {
for (icol = irow+1; icol < mtrxSize; icol++) {
SumOffDiagonal += TMath::Abs(Source(irow,icol));
}
}
if (SumOffDiagonal == 0.0) break;
if (Iteration < 4) {
Threshold = 0.2*SumOffDiagonal/(mtrxSize*mtrxSize);
} else {
Threshold = 0.0;
}
for(irow = 0; irow < mtrxSize; irow++) {
for(icol = irow+1; icol < mtrxSize; icol++) {
Double_t t;
Double_t g;
g = 100.0*TMath::Abs(Source(irow,icol));
if (Iteration > 4 &&
(TMath::Abs(EigenValues(irow,0))+g) ==
TMath::Abs(EigenValues(irow,0)) &&
(TMath::Abs(EigenValues(icol,0))+g) ==
TMath::Abs(EigenValues(icol,0))) {
Source(irow,icol) = 0.0;
} else if (TMath::Abs(Source(irow,icol)) > Threshold) {
Double_t h = EigenValues(icol,0) - EigenValues(irow,0);
if ((TMath::Abs(h) + g) == TMath::Abs(h) ) {
t = Source(irow,icol)/h;
} else {
Double_t Theta = 0.5*h/Source(irow,icol);
t = 1.0/(TMath::Abs(Theta)
+ TMath::Sqrt(1.0 + Theta*Theta));
if (Theta < 0) t = -t;
}
Double_t c = 1.0 / sqrt(1.0 + t*t);
Double_t s = t * c;
Double_t tau = s / (1.0 + c);
h = t * Source(irow,icol);
Z(irow,0) -= h;
Z(icol,0) += h;
EigenValues(irow,0) -= h;
EigenValues(icol,0) += h;
Source(irow,icol) = 0.0;
Int_t j;
for(j = 0; j < irow; j++) {
JACOBI_ROTATION(Source,j,irow,j,icol);
}
for(j = irow + 1; j < icol; j++) {
JACOBI_ROTATION(Source,irow,j,j,icol);
}
for(j = icol + 1; j < mtrxSize; j++) {
JACOBI_ROTATION(Source,irow,j,icol,j);
}
for(j = 0;j < mtrxSize; j++) {
JACOBI_ROTATION(RotationMatrix,j,irow,j,icol);
}
NumberRotations++;
}
} // for(icol = irow+1; icol < mtrxSize; icol++)
} // for(irow = 0; irow < mtrxSize; irow++)
B += Z;
EigenValues = B;
Z.Zero();
} // for (G4int Iteration = 0; Iteration < 50; Iteration++)
}
//______________________________________________________________________
//______________________________________________________________________
// Get Closest Approach Point to beam axis...
void
LCDSmearTrack::CalcPOCA(
Double_t qMC,
Double_t bfld_z,
TVector3 p_orig,
TVector3 x_orig,
TVector3& p_poca,
TVector3& x_poca) {
Double_t pT = p_orig.Pt();
Double_t tanL = p_orig.Z()/pT;
Double_t arho = 333.567/bfld_z*pT;
TVector3 nzv(0.0,0.0,qMC/TMath::Abs(qMC));
TVector3 xxc = p_orig.Cross(nzv);
xxc.SetZ(0.0);
TVector3 nxxc = xxc.Unit();
TVector3 xc = x_orig + arho*nxxc;
xc.SetZ(0.0);
TVector3 nxc = xc.Unit();
TVector3 catMC=nzv.Cross(nxc);
catMC.SetZ(0.0);
TVector3 ncatMC=catMC.Unit();
p_poca = pT*ncatMC;
p_poca.SetZ(p_orig.Z());
Double_t dphi = TMath::ASin(nxxc.Cross(nxc).Z());
Double_t dl =-dphi*arho*qMC;
x_poca = xc - arho*nxc;
x_poca.SetZ(x_orig.Z() + dl*tanL);
}
//______________________________________________________________________
//______________________________________________________________________
// Get resolution values from interpolation and smear parameters
void
LCDSmearTrack::CalcSmearParameters(
Double_t* m_eparm,
Double_t* d_parm,
Double_t bfld_z,
Double_t pTot,
Double_t abscth,
LCDGetParameters* gp,
LCDLookUp2d** par,
Double_t* chi2) {
Int_t rc;
// Save the error matrix elements
TMatrixD m_err(5,5);
Int_t i,j;
Int_t ij;
ij=0;
for (i=0 ; i < 5 ; i++) {
for (j=0 ; j <= i ; j++) {
if (par[ij] != 0 ) {
m_err(i,j) =
par[ij]->interpolateVal((Double_t)abscth,(Double_t)pTot,&rc);
} else {
m_err(i,j) = 0;
}
if (j != i) {
m_err(j,i) = m_err(i,j);
}
ij++;
}
}
// for (i=0 ; i < 5 ; i++) {
// m_err(i,2) = m_err(i,2)*333.567/bfld_z;
// m_err(2,i) = m_err(2,i)*333.567/bfld_z;
// }
TMatrixD m_TMP0=m_err;
Double_t detm_TMP0;
m_TMP0.Invert(&detm_TMP0);
if (detm_TMP0 == 0) {
// Error Matrix can not be inverted!!!
// print out error matrix elements and stop now!!!
printf("!!!!!!!!!!!!!!!!!!!!!!!!n");
printf("Error Matrix can not be inverted!!!n");
printf("!!!!!!!!!!!!!!!!!!!!!!!!n");
printf(" n");
printf(" pTot =%f n",pTot);
printf(" Abs(cth)=%fn",abscth);
printf(" error matrix elements....n");
for (i=0 ; i < 5 ; i++) {
printf(" %10.5e %10.5e %10.5e %10.5e %10.5en",
m_err(0,i),m_err(1,i),m_err(2,i),m_err(3,i),m_err(4,i));
}
exit(-1);
}
ij = 0;
for (i=0 ; i < 5 ; i++) {
for (j=0 ; j <= i ; j++) {
m_eparm[ij++] = m_err(i,j);
}
}
TMatrixD EigenValues(5,1);
TMatrixD m_errp(5,5);
TMatrixD Transform(5,5);
m_errp=m_err;
FindEigen(m_errp,Transform,EigenValues);
TMatrixD TransInv(5,5);
Double_t detTransform;
TransInv=Transform;
TransInv.Invert(&detTransform);
if (detTransform != 0) {
for(Int_t elem = 0; elem < 5; elem++) {
if(EigenValues(elem,0) < 0.0) {
EigenValues(elem,0) = 0.0;
}
}
// Get Gaussian distributed random number for calculating how
// much to smear the 5 track parameters (widths of Gaussians given
// by the EigenValues).
// Mean = 0, width given by second arg
Double_t mean = 0.0;
TMatrixD Deviations(5,1);
Deviations.Zero();
Double_t deigen;
for (i=0 ; i < 5 ; i++) {
deigen=EigenValues(i,0);
Deviations(i,0)=m_random.Gaus(mean,TMath::Sqrt(deigen));
}
// Transform back to the original basis
TMatrixD m_Smear(5,1);
m_Smear.Mult(Transform,Deviations);
for (i=0 ; i < 5 ; i++) {
d_parm[i] = m_Smear(i,0);
}
TMatrixD m_SmearT(1,5);
for (i=0 ; i < 5 ; i++) {
m_SmearT(0,i) = m_Smear(i,0);
}
TMatrixD m_TMP1(5,1);
m_TMP1.Mult(m_TMP0,m_Smear);
TMatrixD m_TMP2(1,1);
m_TMP2.Mult(m_SmearT,m_TMP1);
*chi2=m_TMP2(0,0);
} else {
for (i=0 ; i < 5 ; i++) {
d_parm[i] = 0;
}
}
}
//______________________________________________________________________
//______________________________________________________________________
// Get track parameter
void
LCDSmearTrack::CalcTrackParameters(Double_t qMC,
Double_t bfld_z,
TVector3 pMC,
TVector3 xMC,
Double_t* tkpar) {
TVector3 capMC;
TVector3 caxMC;
CalcPOCA(qMC,bfld_z,pMC,xMC,capMC,caxMC);
Double_t pT = capMC.Pt();
Double_t d0= caxMC.Perp();
TVector3 caxMCXcapMC=caxMC.Cross(capMC);
if (caxMCXcapMC.Z() > 0) {
d0 = -d0;
}
Double_t phi0 = capMC.Phi();
// capMC.X() == 0 && capMC.Y() == 0 ? 0.0 : atan2(capMC.Y()/pT,capMC.X()/pT);
Double_t z0 = caxMC.Z();
Double_t tanL = capMC.Z()/pT;
tkpar[0]= d0;
tkpar[1]= phi0;
tkpar[2]=-qMC*bfld_z/333.567/pT;
tkpar[3]= z0;
tkpar[4]= tanL;
}
//______________________________________________________________________
ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.