00001
00010
00011 #include "MessageService/MsgService.h"
00012
00013 #include "CandFitTrackSA/DataFT.h"
00014 #include "CandFitTrackSA/FitContext.h"
00015 #include "CandFitTrackSA/FitState.h"
00016 #include "CandFitTrackSA/FitStateFactory.h"
00017 #include "CandFitTrackSA/Ntp/NtpFitSA.h"
00018 #include "CandFitTrackSA/TracerSA.h"
00019 #include "CandFitTrackSA/TrackEstimatorRange.h"
00020
00021 #include <cmath>
00022
00023 CVSID("$Id: FitContext.cxx,v 1.4 2006/12/01 18:40:45 rhatcher Exp $");
00024
00025 using namespace ConstFT;
00026
00030 FitContext::FitContext(const AlgConfig& ac, const TrackContext& context,
00031 TrackEstimator* estimator, SwimSwimmer* swimmer ) :
00032 fTrackContext(context), fState(0),
00033 fSwimmer(swimmer), fEstimator(estimator),
00034 fData(ac, context),
00035 fConvergenceMaster(),
00036 fMatCalc(ac, context),
00037 fNIterationsStep(0), fNIterationsTotal(0),
00038
00039 fNPlanesFit(0),
00040 fCurrentFit(ConstFT::NTrackParams), fTimesConverged(0), fLastGoodFit(),
00041 fDChi2(ConstFT::InitialDChi2), fNTriesDiverges(0), fNHitsInViewMin(4),
00042 fNMaxIterations(100), fNMaxDiverging(1), fConvergenceCond(0.1)
00043 {
00044 TracerSA trace(
00045 "FitContext::FitContext(const AlgConfig&,const TrackContext&,SwimSwimmer*)"
00046 );
00047
00048 Config(ac);
00049
00050 fState = FitStateFactory::Instance().GetFitState(std::string("Initial"));
00051
00052 fConvergenceMaster.SetNHitsInViewMin(fNHitsInViewMin);
00053 fConvergenceMaster.BuildMasks(fData);
00054 }
00055
00056
00060 FitContext::~FitContext()
00061 {
00062 TracerSA trace("FitContext::~FitContext()");
00063 }
00064
00065
00069 void FitContext::Config(const AlgConfig& ac)
00070 {
00071 TracerSA trace("FitContext::Config(const AlgConfig&)");
00072
00073 if ( ac.KeyExists("NHitsInViewMin") ) {
00074 fNHitsInViewMin = ac.GetInt("NHitsInViewMin");
00075 }
00076 if ( ac.KeyExists("NMaxIterations") ) {
00077 fNMaxIterations = ac.GetInt("NMaxIterations");
00078 }
00079 if ( ac.KeyExists("NMaxDiverging") ) {
00080 fNMaxDiverging = ac.GetInt("NMaxDiverging");
00081 }
00082 if ( ac.KeyExists("ConvergenceCond") ) {
00083 fConvergenceCond = ac.GetDouble("ConvergenceCond");
00084 }
00085 }
00086
00087
00091 void FitContext::Iterate()
00092 {
00093 TracerSA trace("FitContext::Iterate()");
00094 fState->Iterate(*this);
00095 }
00096
00097
00101 void FitContext::SetState(const FitState* state)
00102 {
00103 TracerSA trace("FitContext::SetState(const FitState*)");
00104 fState = state;
00105 MSG("FitTrackSA",Msg::kDebug) << "Switched to "
00106 << state->Name() << " state.\n";
00107 }
00108
00109
00113 NtpFitSA FitContext::MakeNtpFitSA() const
00114 {
00115 TracerSA trace("FitContext::MakeNtpFitSA()");
00116
00117 NtpFitSA ntp;
00118
00119
00120 ntp.zdir = fTrackContext.GetDir();
00121 FillNtpFitSA(ntp);
00122 FillNtpPlaneInfo(ntp);
00123 FillNtpBFieldCalib(ntp);
00124 FillNtpFitSR(ntp);
00125
00126 return ntp;
00127 }
00128
00129
00133 void FitContext::FillNtpFitSA(NtpFitSA& ntp) const
00134 {
00135 TracerSA trace("FitContext::FillNtpFitSA()");
00136
00137 if ( fTimesConverged == 0 ) {
00138 ntp.fit.pass = 0;
00139 } else {
00140 ntp.fit.pass = 1;
00141 }
00142 ntp.fit.niter = fNIterationsTotal;
00143 ntp.fit.ndf = fLastGoodFit.GetNdof();
00144 ntp.fit.q = fLastGoodFit.GetQ();
00145 ntp.fit.chi2 = fLastGoodFit.GetChi2();
00146 if ( TMath::Abs(ntp.fit.ndf) > TinyNumber ) {
00147 ntp.fit.rchi2 = ntp.fit.chi2/ntp.fit.ndf;
00148 } else {
00149 ntp.fit.rchi2 = -1.;
00150 }
00151
00152 ntp.fit.cpu = fCpu;
00153 ntp.fit.p = fLastGoodFit.GetP();
00154 ntp.fit.ep = fLastGoodFit.GetEP();
00155
00156 ntp.fit.u = fLastGoodFit.GetFitParameter(kU);
00157 ntp.fit.eu = fLastGoodFit.GetFitParameterError(kU);
00158 ntp.fit.dudz = fLastGoodFit.GetFitParameter(kdUdZ);
00159 ntp.fit.edudz = fLastGoodFit.GetFitParameterError(kdUdZ);
00160 ntp.fit.v = fLastGoodFit.GetFitParameter(kV);
00161 ntp.fit.ev = fLastGoodFit.GetFitParameterError(kV);
00162 ntp.fit.dvdz = fLastGoodFit.GetFitParameter(kdVdZ);
00163 ntp.fit.edvdz = fLastGoodFit.GetFitParameterError(kdVdZ);
00164 ntp.fit.qp = fLastGoodFit.GetFitParameter(kQoverP);
00165 ntp.fit.eqp = fLastGoodFit.GetFitParameterError(kQoverP);
00166
00167 ntp.fit.z = fData.GetZVtx();
00168
00169
00170 Float_t dcosu, dcosv, dcosz;
00171 dcosz = fTrackContext.GetDir()*
00172 pow(1. + pow(fLastGoodFit.GetTrackOut(kdUdZ),2)
00173 + pow(fLastGoodFit.GetTrackOut(kdVdZ),2),-0.5);
00174 dcosu = fLastGoodFit.GetTrackOut(kdUdZ)*dcosz;
00175 dcosv = fLastGoodFit.GetTrackOut(kdVdZ)*dcosz;
00176 ntp.fit.dcosu = dcosu;
00177 ntp.fit.edcosu = sqrt(fabs(dcosz))*sqrt(pow(dcosu*dcosv*ntp.fit.edvdz,2)+
00178 pow((pow(dcosz,2)+pow(dcosv,2))*ntp.fit.edudz,2));
00179 ntp.fit.dcosv = dcosv;
00180 ntp.fit.edcosv = sqrt(fabs(dcosz))*sqrt(pow(dcosu*dcosv*ntp.fit.edudz,2)+
00181 pow((pow(dcosz,2)+pow(dcosu,2))*ntp.fit.edvdz,2));
00182 ntp.fit.dcosz = dcosz;
00183
00184
00185 for ( int i = 0; i < NTrackParams; i++ ) {
00186 ntp.fit.par[i] = fLastGoodFit.GetFitParameter(i);
00187 for ( int j = 0; j < NTrackParams; j++ ) {
00188 ntp.fit.parerr[i][j] = fLastGoodFit.GetFitErrM(i,j);
00189 }
00190 }
00191 }
00192
00193
00197 void FitContext::FillNtpBFieldCalib(NtpFitSA& ntp) const
00198 {
00199 TracerSA trace("FitContext::FillNtpBFieldCalib()");
00200
00201 ntp.bfcalib.prange = fTrackContext.GetPrange();
00202
00203 ntp.bfcalib.umean = fData.GetUMean();
00204 ntp.bfcalib.vmean = fData.GetVMean();
00205 }
00206
00207
00211 void FitContext::FillNtpFitSR(NtpFitSA& ntp) const
00212 {
00213 TracerSA trace("FitContext::FillNtpFitSR()");
00214 ntp.fitsr.pass = fTrackContext.GetPassSR();
00215 ntp.fitsr.niter = fTrackContext.GetIterationsSR();
00216 ntp.fitsr.ndf = fTrackContext.GetNdofSR();
00217 ntp.fitsr.q = fTrackContext.GetEMChargeSR();
00218 ntp.fitsr.chi2 = fTrackContext.GetChi2SR();
00219 ntp.fitsr.rchi2 = fTrackContext.GetRChi2SR();
00220 ntp.fitsr.cpu = fTrackContext.GetCpuSR();
00221 ntp.fitsr.qp = fTrackContext.GetQPSR();
00222 ntp.fitsr.eqp = fTrackContext.GetEQPSR();
00223 ntp.fitsr.u = fTrackContext.GetUSR();
00224 ntp.fitsr.eu = fTrackContext.GetUErrorSR();
00225 ntp.fitsr.v = fTrackContext.GetVSR();
00226 ntp.fitsr.ev = fTrackContext.GetVErrorSR();
00227 ntp.fitsr.dcosu = fTrackContext.GetDirCosUSR();
00228 ntp.fitsr.edcosu = fTrackContext.GetEDirCosUSR();
00229 ntp.fitsr.dcosv = fTrackContext.GetDirCosVSR();
00230 ntp.fitsr.edcosv = fTrackContext.GetEDirCosVSR();
00231 ntp.fitsr.dcosz = fTrackContext.GetDirCosZSR();
00232 ntp.fitsr.p = fTrackContext.GetMomentumSR();
00233 ntp.fitsr.prange = fTrackContext.GetMomentumRangeSR();
00234 }
00235
00236
00240 void FitContext::FillNtpPlaneInfo(NtpFitSA& ntp) const
00241 {
00242 TracerSA trace("FitContext::FillNtpPlaneInfo()");
00243 ntp.plane.begin = fData.GetBegPlaneId().GetPlane();
00244 ntp.plane.end = fData.GetEndPlaneId().GetPlane();
00245 ntp.plane.n = fData.GetNPlanes();
00246 ntp.plane.nu = fData.GetNUPlanes();
00247 ntp.plane.nv = fData.GetNVPlanes();
00248
00249 ntp.hitplane.begin = fData.GetBegHit();
00250 ntp.hitplane.end = fData.GetEndHit();
00251 ntp.hitplane.n = fData.GetNUHits() + fData.GetNVHits();
00252 ntp.hitplane.nu = fData.GetNUHits();
00253 ntp.hitplane.nv = fData.GetNVHits();
00254
00255 ntp.fitplane.begin = fData.GetBegHitUsed();
00256 ntp.fitplane.end = fData.GetEndHitUsed();
00257 ntp.fitplane.n = fData.GetNUHitsUsed() + fData.GetNVHitsUsed();
00258 ntp.fitplane.nu = fData.GetNUHitsUsed();
00259 ntp.fitplane.nv = fData.GetNVHitsUsed();
00260 }
00261
00265 void FitContext::SetFitParams(const TVectorD& fit)
00266 {
00267 TracerSA trace("FitContext::SetFitParams(const TVectorD&)");
00268 fCurrentFit = fit;
00269 PrintCurrentFit();
00270 }
00271
00275 void FitContext::SetFromLastGoodFit()
00276 {
00277 TracerSA trace("FitContext::SetFromLastGoodFit()");
00278 fCurrentFit = fLastGoodFit.GetTrackOut();
00279 }
00280
00284 void FitContext::SetLastGoodFitParams(const TVectorD& fit)
00285 {
00286 TracerSA trace("FitContext::SetLastGoodFitParams(const TVectorD&)");
00287 fLastGoodFit.SetTrackOut(fit);
00288 }
00289
00293 void FitContext::PrintCurrentFit()
00294 {
00295 TracerSA trace("FitContext::PrintCurrentFit()");
00296 MsgStream *mftsa = &MSGSTREAM("FitTrackSA", Msg::kInfo);
00297
00298 (*mftsa) << "U=" << fCurrentFit(kU) << "; dU/dZ=" << fCurrentFit(kdUdZ)
00299 << "; V=" << fCurrentFit(kV) << "; dV/dZ=" << fCurrentFit(kdVdZ)
00300 << "; Q/P=" << fCurrentFit(kQoverP) << "\n";
00301 }
00302
00306 void FitContext::Print()
00307 {
00308 TracerSA trace("FitContext::Print()");
00309 MsgStream *mftsa = &MSGSTREAM("FitTrackSA", Msg::kInfo);
00310
00311 (*mftsa) << "Step " << fNIterationsStep << "/" << fNIterationsTotal
00312
00313 << "; #ToFit=" << fConvergenceMaster.GetNPlanesCur()
00314 << "; #Fit=" << fNPlanesFit << "; dchi2=" << fDChi2
00315 << "; q/p=" << fCurrentFit(kQoverP) << "\n";
00316
00317 }
00318
