Inheritance diagram for gar::rec::dayonetrackfit:

Public Member Functions
	dayonetrackfit (fhicl::ParameterSet const &p)

	dayonetrackfit (dayonetrackfit const &)=delete

	dayonetrackfit (dayonetrackfit &&)=delete

dayonetrackfit &	operator= (dayonetrackfit const &)=delete

dayonetrackfit &	operator= (dayonetrackfit &&)=delete

void	produce (art::Event &e) override

	dayonetrackfit (fhicl::ParameterSet const &p)

	dayonetrackfit (dayonetrackfit const &)=delete

	dayonetrackfit (dayonetrackfit &&)=delete

dayonetrackfit &	operator= (dayonetrackfit const &)=delete

dayonetrackfit &	operator= (dayonetrackfit &&)=delete

void	produce (art::Event &e) override

Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDProducer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept

	Producer (fhicl::ParameterSet const &)

	Producer (Producer const &)=delete

	Producer (Producer &&)=delete

Producer &	operator= (Producer const &)=delete

Producer &	operator= (Producer &&)=delete

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

Public Member Functions inherited from art::Modifier
	~Modifier () noexcept

	Modifier ()

	Modifier (Modifier const &)=delete

	Modifier (Modifier &&)=delete

Modifier &	operator= (Modifier const &)=delete

Modifier &	operator= (Modifier &&)=delete

Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept

	ModuleBase ()

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Member Functions
int	KalmanFit (std::vector< TPCCluster > &TPCClusters, std::vector< int > &TPCClusterlist, std::vector< float > &trackparatend, float &chisquared, float &length, float *covmat, std::set< int > &unused_TPCClusters, std::vector< std::pair< float, float >> &dSigdX, std::vector< TVector3 > &trajpts)

int	KalmanFitBothWays (std::vector< gar::rec::TPCCluster > &TPCClusters, TrackPar &trackpar, TrackIoniz &trackions, TrackTrajectory &tracktraj)

void	updatepar (TVectorF &parvec, float zcur, float xh, float yh, float zh, TVectorF &predstep, float &dlength)

int	DayOneFitBothWays (const gar::rec::Track &patrectrack, std::vector< gar::rec::TPCCluster > &TPCClusters, TrackPar &trackpar, TrackIoniz &trackions, TrackTrajectory &tracktraj)

int	DayOneFit (std::vector< TPCCluster > &TPCClusters, std::vector< int > &TPCClusterlist, std::vector< float > &trackparbeg, float &chisquared, float &length, float *covmat, std::vector< std::pair< float, float >> &dSigdXs, std::vector< TVector3 > &trajpts)

void	updatepar (TVectorF &parvec, float zcur, float xh, float yh, float zh, TVectorF &predstep, float &dlength)

Private Attributes
std::string	fPatRecLabel
	input patrec tracks and associations More...

int	fPrintLevel
	debug printout: 0: none, 1: just track parameters and residuals, 2: all More...

float	fKalCurvStepUncSq
	constant uncertainty term on each step of the Kalman fit – squared, for curvature More...

float	fKalPhiStepUncSq
	constant uncertainty term on each step of the Kalman fit – squared, for phi More...

float	fKalLambdaStepUncSq
	constant uncertainty term on each step of the Kalman fit – squared, for lambda More...

float	fTPCClusterResolYZ
	pad size in cm in YZ to determine step size More...

float	fTPCClusterResolX
	drift direction contribution to determine step size (resolution of a TPCCluster) More...

unsigned int	fInitialTPNTPCClusters
	number of TPCClusters to use for initial trackpar estimate, if present More...

unsigned int	fMinNumTPCClusters
	minimum number of TPCClusters to define a track More...

int	fDumpTracks
	0: do not print out tracks, 1: print out tracks More...

float	fTPCClusterResolYZinFit
	TPCCluster resolution parameter to use in fit. More...

float	fRoadYZinFit
	cut in cm for dropping TPCClusters from tracks in fit More...

float	fSortTransWeight
	for use in the hit sorting algorithm – transverse distance weight factor More...

float	fSortDistBack
	for use in the hit sorting algorithm – how far to go back before raising the distance figure of merit More...

float	fMinIonizGapCut
	Don't compute dEdx for this dx or larger. More...

float	fTPCClusterResolXY
	pad size in cm in YZ to determine step size More...

float	fTPCClusterResolZ
	drift direction contribution to determine step size (resolution of a TPCCluster) More...

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >

Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)

Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 49 of file dayonetrackfit_module.cc.

Constructor & Destructor Documentation

gar::rec::dayonetrackfit::dayonetrackfit ( fhicl::ParameterSet const & p )

explicit

Definition at line 102 of file dayonetrackfit_module.cc.

                                                              : EDProducer{p}  
     {
       // Call appropriate produces<>() functions here.
       // Call appropriate consumes<>() for any products to be retrieved by this module.
 
       fPatRecLabel       = p.get<std::string>("PatRecLabel","patrec");
       fPrintLevel        = p.get<int>("PrintLevel",0);
       fMinNumTPCClusters        = p.get<unsigned int>("MinNumTPCClusters",20);
       fKalCurvStepUncSq  = p.get<float>("KalCurvStepUncSq",1.0E-9);
       fKalPhiStepUncSq   = p.get<float>("KalPhiStepUncSq",1.0E-9);
       fKalLambdaStepUncSq = p.get<float>("KalLambdaStepUncSq",1.0E-9);
       fInitialTPNTPCClusters    = p.get<unsigned int>("InitialTPNTPCClusters",100);
       fDumpTracks        = p.get<int>("DumpTracks",0);
       fTPCClusterResolYZ        = p.get<float>("TPCClusterResolYZ",1.0); // TODO -- think about what this value is
       fTPCClusterResolX         = p.get<float>("TPCClusterResolX",0.5);  // this is probably much better
       fTPCClusterResolYZinFit   = p.get<float>("TPCClusterResolYZinFit",4.0);
       fRoadYZinFit       = p.get<float>("RoadYZinFit",1.0);
       fSortTransWeight   = p.get<float>("SortTransWeight",0.1);
       fSortDistBack      = p.get<float>("SortDistBack",2.0);
       fMinIonizGapCut    = p.get<float>("MinIonizGapCut",5.0);
 
       art::InputTag patrecTag(fPatRecLabel);
       consumes< std::vector<gar::rec::Track> >(patrecTag);
       consumes< art::Assns<gar::rec::TPCCluster, gar::rec::Track> >(patrecTag);
 
       // probably don't need the vector hits at this point if we have the TPCClusters
       //consumes< std::vector<gar::rec::VecHit> >(patrecTag);
       //consumes< art::Assns<gar::rec::VecHit, gar::rec::Track> >(patrecTag);
 
       produces< std::vector<gar::rec::Track> >();
       produces< art::Assns<gar::rec::TPCCluster, gar::rec::Track> >();
       produces<std::vector<gar::rec::TrackIoniz>>();
       produces<art::Assns<rec::TrackIoniz, rec::Track>>();
       produces<std::vector<gar::rec::TrackTrajectory>>();
       produces<art::Assns<rec::TrackTrajectory, rec::Track>>();
     }

gar::rec::dayonetrackfit::dayonetrackfit ( dayonetrackfit const & )

delete

gar::rec::dayonetrackfit::dayonetrackfit ( dayonetrackfit && )

delete

gar::rec::dayonetrackfit::dayonetrackfit ( fhicl::ParameterSet const & p )

explicit

gar::rec::dayonetrackfit::dayonetrackfit ( dayonetrackfit const & )

delete

gar::rec::dayonetrackfit::dayonetrackfit ( dayonetrackfit && )

delete

Member Function Documentation

int gar::rec::dayonetrackfit::DayOneFit	(	std::vector< TPCCluster > &	TPCClusters,
		std::vector< int > &	TPCClusterlist,
		std::vector< float > &	trackparbeg,
		float &	chisquared,
		float &	length,
		float *	covmat,
		std::vector< std::pair< float, float >> &	dSigdXs,
		std::vector< TVector3 > &	trajpts
	)

private

Definition at line 285 of file dayonetrackfit_module.cc.

     {
 
       // set some default values in case we return early
 
       size_t nclus = TPCClusterlist.size();
       chisquared = 0;
       length = 0;
 
       // first try, do nothing, but we should calculate a chisquared
 
       for (size_t i=0;i<nclus; ++i)
         {
           const float *pos = TPCClusters.at(TPCClusterlist.at(i)).Position();
           TVector3 tvec(pos[0],pos[1],pos[2]);
           trajpts.push_back(tvec);
           float dist = 0;
           float tvpos[3] = {trackparbeg[5],trackparbeg[0],trackparbeg[1]};
           int retcode = util::TrackPropagator::DistXYZ(trackparbeg.data(),tvpos,pos,dist);
           if (retcode == 0) 
             { 
               chisquared += TMath::Sq(dist/fTPCClusterResolXY);  // todo -- break this up into X,Y vs. Z pieces.
             }
         }
 
       return 0;
     }

int gar::rec::dayonetrackfit::DayOneFitBothWays	(	const gar::rec::Track &	patrectrack,
		std::vector< gar::rec::TPCCluster > &	TPCClusters,
		TrackPar &	trackpar,
		TrackIoniz &	trackions,
		TrackTrajectory &	tracktraj
	)

private

Definition at line 192 of file dayonetrackfit_module.cc.

     {
       // track parameters:  x is the independent variable in TPC track fits, but just report it here.
       // 0: y
       // 1: z
       // 2: curvature
       // 3: phi
       // 4: lambda = angle from the cathode plane
       // 5: x   /// added on to the end
 
       size_t nclus = TPCClusters.size();
       std::vector<float> TPCClusterz;
       for (size_t i=0; i<nclus; ++i)
         {
           TPCClusterz.push_back(TPCClusters[i].Position()[2]);
         }
       std::vector<int> hsi(nclus);   // "hit sort index" -- really TPC clusters
       TMath::Sort((int) nclus,TPCClusterz.data(),hsi.data());
       std::vector<int> hsib(nclus);
       for (size_t i=0; i<nclus; ++i)
         {
           hsib[nclus-i-1] = hsi[i];
         }
 
       std::vector<float> tparbeg(6);  // fill with initial track parameters from patrec track
       for (size_t i=0; i<5; ++i)
         {
           tparbeg[i] = patrectrack.TrackParBeg()[i];
         }
       tparbeg[5] = patrectrack.Vertex()[0];
       float covmatbeg[25];
       patrectrack.CovMatBegSymmetric(covmatbeg);
       float chisqforwards = 0;
       float lengthforwards = 0;
       std::vector<std::pair<float,float>> dSigdXs_FWD;
       std::vector<TVector3> trajpts_FWD;
 
       int retcode = DayOneFit(TPCClusters,hsi,tparbeg,chisqforwards,lengthforwards,covmatbeg,dSigdXs_FWD,trajpts_FWD);
       if (retcode != 0) return 1;
 
 
       std::vector<float> tparend(6);  // fill with end track parameters from patrec track
       for (size_t i=0; i<5; ++i)
         {
           tparend[i] = patrectrack.TrackParEnd()[i];
         }
       tparend[5] = patrectrack.End()[0];
       float covmatend[25];
       patrectrack.CovMatEndSymmetric(covmatend);
       float chisqbackwards = 0;
       float lengthbackwards = 0;
       std::vector<std::pair<float,float>> dSigdXs_BAK;
       std::vector<TVector3> trajpts_BAK;
 
       retcode = DayOneFit(TPCClusters,hsib,tparend,chisqbackwards,lengthbackwards,covmatend,dSigdXs_BAK,trajpts_BAK);
       if (retcode != 0) return 1;
 
       size_t nTPCClusters=TPCClusters.size();
       trackpar.setNTPCClusters(nTPCClusters);
       trackpar.setTime(0);
       trackpar.setChisqForwards(chisqforwards);
       trackpar.setChisqBackwards(chisqbackwards);
       trackpar.setLengthForwards(lengthforwards);
       trackpar.setLengthBackwards(lengthbackwards);
       trackpar.setCovMatBeg(covmatbeg);
       trackpar.setCovMatEnd(covmatend);
       trackpar.setTrackParametersBegin(tparbeg.data());
       trackpar.setXBeg(tparbeg[5]);
       trackpar.setTrackParametersEnd(tparend.data());
       trackpar.setXEnd(tparend[5]);
 
       trackions.setData(dSigdXs_FWD,dSigdXs_BAK);
       tracktraj.setData(trajpts_FWD,trajpts_BAK);
 
       return 0;
     }

int gar::rec::dayonetrackfit::KalmanFit	(	std::vector< TPCCluster > &	TPCClusters,
		std::vector< int > &	TPCClusterlist,
		std::vector< float > &	trackparatend,
		float &	chisquared,
		float &	length,
		float *	covmat,
		std::set< int > &	unused_TPCClusters,
		std::vector< std::pair< float, float >> &	dSigdX,
		std::vector< TVector3 > &	trajpts
	)

private

Definition at line 298 of file dayonetrackfit_module.cc.

     {
 
       // set some default values in case we return early
 
       size_t nTPCClusters = TPCClusterlist.size();
       chisquared = 0;
       length = 0;
       for (size_t i=0; i<5; ++i) trackparatend[i] = 0;
       for (size_t i=0; i<25; ++i) covmat[i] = 0;
 
       float roadsq = fRoadYZinFit*fRoadYZinFit;
 
       // estimate curvature, lambda, phi, xpos from the initial track parameters
       float curvature_init=0.1;
       float phi_init = 0;
       float lambda_init = 0;
       float xpos_init=0;
       float ypos_init=0;
       float zpos_init=0;
       float x_other_end = 0;
       if ( gar::rec::initial_trackpar_estimate(TPCClusters,
                                                TPCClusterlist,
                                                curvature_init,
                                                lambda_init,
                                                phi_init,
                                                xpos_init,
                                                ypos_init,
                                                zpos_init,
                                                x_other_end,
                                                fInitialTPNTPCClusters, 
                                                fPrintLevel) != 0)
         {
           //std::cout << "kalman fit failed on initial trackpar estimate" << std::endl;
           return 1;
         }
 
       // Kalman fitter variables
 
       float zpos = zpos_init;
 
       TMatrixF P(5,5);  // covariance matrix of parameters
       // fill in initial guesses -- generous uncertainties on first value.
       P.Zero();
       P[0][0] = TMath::Sq(1);   // initial position uncertainties -- y
       P[1][1] = TMath::Sq(1);   // and x
       P[2][2] = TMath::Sq(.5);  // curvature of zero gets us to infinite momentum, and curvature of 2 is curled up tighter than the pads
       P[3][3] = TMath::Sq(.5);  // phi uncertainty
       P[4][4] = TMath::Sq(.5);  // lambda uncertainty
 
       TMatrixF PPred(5,5);
 
       // per-step additions to the covariance matrix
       TMatrixF Q(5,5);
       Q.Zero();
       Q[2][2] = fKalCurvStepUncSq;     // allow for some curvature uncertainty between points
       Q[3][3] = fKalPhiStepUncSq;      // phi
       Q[4][4] = fKalLambdaStepUncSq;   // lambda
 
       // uncertainties on the measured points  (big for now)
       TMatrixF R(2,2);
       R.Zero();
       R[0][0] = TMath::Sq(fTPCClusterResolYZinFit);  // in cm^2
       R[1][1] = TMath::Sq(fTPCClusterResolYZinFit);  // in cm^2
 
       // add the TPCClusters and update the track parameters and uncertainties.  Put in additional terms to keep uncertainties from shrinking when
       // scattering and energy loss can change the track parameters along the way.
 
       // F = partial(updatefunc)/partial(parvec).  Do the drivatives numerically inside the step
 
       TMatrixF F(5,5);
       TMatrixF FT(5,5);
       TVectorF parvec(5);
       parvec[0] = ypos_init;
       parvec[1] = xpos_init;
       parvec[2] = curvature_init;
       parvec[3] = phi_init;
       parvec[4] = lambda_init;
       TVectorF predstep(5);
 
       TMatrixF H(2,5);   // partial(obs)/partial(params)
       H.Zero();
       H[0][0] = 1;  // y
       H[1][1] = 1;  // x
       TMatrixF HT(5,2);
 
       TVectorF z(2);
       TVectorF ytilde(2);
       TVectorF hx(2);
       TMatrixF S(2,2);
       TMatrixF K(5,2);
 
       TMatrixF I(5,5);
       I.Zero();
       for (int i=0;i<5;++i) I[i][i] = 1;
 
       for (size_t iTPCCluster=1; iTPCCluster<nTPCClusters; ++iTPCCluster)
         {
 
           float xh = TPCClusters[TPCClusterlist[iTPCCluster]].Position()[0];
           float yh = TPCClusters[TPCClusterlist[iTPCCluster]].Position()[1];
           float zh = TPCClusters[TPCClusterlist[iTPCCluster]].Position()[2];
 
           if (fPrintLevel > 0)
             {
               std::cout << std::endl;
               std::cout << "Adding a new TPCCluster: " << xh << " " << yh << " " << zh << std::endl;
             }
 
           // for readability
           float lambda = parvec[4];
           float slope = TMath::Tan(lambda);
           if (slope != 0)
             {
               slope = 1.0/slope;
             }
           else
             {
               slope = 1E9;
             }
 
           // look for dx below -- need new Z instead
 
           if (fPrintLevel > 0)
             {
               std::cout << "z: " << zpos << " adding new hit at z: " << zh <<  std::endl;
               std::cout << " Parvec:   y " << parvec[0] << " z " << parvec[1] << " c " << parvec[2] << " phi " << parvec[3] << " lambda " << parvec[4] << std::endl;
             }
 
           // update prediction to the plane containing zh
 
           float dlength=0;
           updatepar(parvec,zpos,xh,yh,zh,predstep,dlength);
 
           if (fPrintLevel > 1)
             {
               std::cout << " Predstep: y " << predstep[0] << " z " << predstep[1] << " c " << predstep[2] << " phi " << predstep[3] << " lambda " << predstep[4] << std::endl;
             }
 
           // calculate F with numerical differentiation
 
           F.Zero();
           TVectorF parvec_plus_dpar(5);
           TVectorF predstepmod(5);
           TVectorF derivvec(5);
           for (size_t ipar=0;ipar<5;++ipar)
             {
               float epsilon = 0.001;
               parvec_plus_dpar = parvec;
               parvec_plus_dpar[ipar] += epsilon;  
               float dldummy=0;
               updatepar(parvec_plus_dpar,zpos,xh,yh,zh,predstepmod,dldummy);
               derivvec = (predstepmod - predstep);
               derivvec *= (1.0/epsilon);
               for (size_t jpar=0; jpar<5; ++jpar)
                 {
                   F[jpar][ipar] = derivvec[jpar];
                 }
             }
 
           if (fPrintLevel > 1)
             {
               std::cout << "F Matrix: " << std::endl;
               F.Print();
               std::cout << "P Matrix: " << std::endl;
               P.Print();
             }
 
           // equations from the extended Kalman filter
           FT.Transpose(F);
           PPred = F*P*FT + Q;
           if (fPrintLevel > 1)
             {
               std::cout << "PPred Matrix: " << std::endl;
               PPred.Print();
             }
 
           ytilde[0] = yh - predstep[0];
           ytilde[1] = xh - predstep[1];
           float ydistsq = ytilde.Norm2Sqr();
           if (ydistsq > roadsq)
             {
               unused_TPCClusters.insert(iTPCCluster);
               continue;
             }
           chisquared += ytilde.Norm2Sqr()/TMath::Sq(fTPCClusterResolYZinFit);
           if (fPrintLevel > 0)
             {
               std::cout << "ytilde (residuals): " << std::endl;
               ytilde.Print();
             }
           if (fPrintLevel > 1)
             {
               std::cout << "H Matrix: " << std::endl;
               H.Print();
             }
 
           HT.Transpose(H);
           S = H*PPred*HT + R;
           if (fPrintLevel > 1)
             {
               std::cout << "S Matrix: " << std::endl;
               S.Print();
             }
 
           S.Invert();
           if (fPrintLevel > 1)
             {
               std::cout << "Inverted S Matrix: " << std::endl;
               S.Print();
             }
 
           K = PPred*HT*S;
           if (fPrintLevel > 1)
             {
               std::cout << "K Matrix: " << std::endl;
               K.Print();
             }
 
           parvec = predstep + K*ytilde;
           P = (I-K*H)*PPred;
           zpos = zh;
           //std::cout << " Updated xpos: " << xpos << " " << dx << std::endl;
           trajpts.emplace_back(parvec[1],parvec[0],zpos);
 
           length += dlength;
 
           // Save the ionization data - skip large gaps from sector boundaries
           float valSig = TPCClusters[TPCClusterlist[iTPCCluster]].Signal();
           if (dlength < fMinIonizGapCut)
             {
               std::pair pushme = std::make_pair(valSig,dlength);
               dSigdXs.push_back( pushme );
             }
           else
             // Have to remove the fellow before the large gap, too
             {
               if (dSigdXs.size()>0) dSigdXs.pop_back();
             }
         }
 
       for (size_t i=0; i<5; ++i)
         {
           trackparatend[i] = parvec[i];
         }
       trackparatend[5] = zpos;  // tack this on so we can specify where the track endpoint is
       if (fPrintLevel > 1)
         {
           std::cout << "Track params at end (y, z, curv, phi, lambda) " << trackparatend[0] << " " << trackparatend[1] << " " <<
             trackparatend[2] << " " << trackparatend[3] <<" " << trackparatend[4] << std::endl;
           S.Print();
         }
 
       // just for visualization of the initial track parameter guesses.  
       //  Comment out when fitting tracks.
       //trackparatend[0] = ypos_init;
       //trackparatend[1] = zpos_init;
       //trackparatend[2] = curvature_init;
       //trackparatend[3] = phi_init;
       //trackparatend[4] = lambda_init;
       //trackparatend[5] = xpos_init;
 
       // no covariance matrix yet
 
       size_t icov=0;
       for (size_t i=0; i<5; ++i)
         {
           for (size_t j=0; j<5; ++j)
             {
               covmat[icov] = P[i][j];
             }
         }
 
       return 0;
     }

int gar::rec::dayonetrackfit::KalmanFitBothWays	(	std::vector< gar::rec::TPCCluster > &	TPCClusters,
		TrackPar &	trackpar,
		TrackIoniz &	trackions,
		TrackTrajectory &	tracktraj
	)

private

Definition at line 210 of file dayonetrackfit_module.cc.

     {
       // variables:  x is the independent variable
       // 0: y
       // 1: z
       // 2: curvature
       // 3: phi
       // 4: lambda = angle from the cathode plane
       // 5: x   /// added on to the end
 
 
       std::vector<int> hlf;
       std::vector<int> hlb;
       float lftmp=0;  // need these dummy arguments so we can share code with the patrec sorter
       float lbtmp=0;
       gar::rec::sort_TPCClusters_along_track(TPCClusters,hlf,hlb,fPrintLevel,lftmp,lbtmp,fSortTransWeight,fSortDistBack);
 
       std::vector<float> tparend(6);
       float covmatend[25];
       float chisqforwards = 0;
       float lengthforwards = 0;
       std::set<int> unused_TPCClusters;
       std::vector<std::pair<float,float>> dSigdXs_FWD;
       std::vector<TVector3> trajpts_FWD;
 
       int retcode = KalmanFit(TPCClusters,hlf,tparend,chisqforwards,lengthforwards,covmatend,unused_TPCClusters,dSigdXs_FWD,trajpts_FWD);
       if (retcode != 0) return 1;
 
       // the "backwards" fit is in decreasing x.  Track parameters are at the end of the fit, the other end of the track
 
       std::vector<float> tparbeg(6);
       float covmatbeg[25];
       float chisqbackwards = 0;
       float lengthbackwards = 0;
       std::vector<std::pair<float,float>> dSigdXs_BAK;
       std::vector<TVector3> trajpts_BAK;
 
       retcode = KalmanFit(TPCClusters,hlb,tparbeg,chisqbackwards,lengthbackwards,covmatbeg,unused_TPCClusters,dSigdXs_BAK,trajpts_BAK);
       if (retcode != 0) return 1;
 
       size_t nTPCClusters=0;
       if (TPCClusters.size()>unused_TPCClusters.size())
         { nTPCClusters = TPCClusters.size()-unused_TPCClusters.size(); }
       trackpar.setNTPCClusters(nTPCClusters);
       trackpar.setTime(0);
       trackpar.setChisqForwards(chisqforwards);
       trackpar.setChisqBackwards(chisqbackwards);
       trackpar.setLengthForwards(lengthforwards);
       trackpar.setLengthBackwards(lengthbackwards);
       trackpar.setCovMatBeg(covmatbeg);
       trackpar.setCovMatEnd(covmatend);
       trackpar.setTrackParametersBegin(tparbeg.data());
       trackpar.setXBeg(tparbeg[5]);
       trackpar.setTrackParametersEnd(tparend.data());
       trackpar.setXEnd(tparend[5]);
 
       trackions.setData(dSigdXs_FWD,dSigdXs_BAK);
       tracktraj.setData(trajpts_FWD,trajpts_BAK);
 
       return 0;
     }

dayonetrackfit& gar::rec::dayonetrackfit::operator= ( dayonetrackfit const & )

delete

dayonetrackfit& gar::rec::dayonetrackfit::operator= ( dayonetrackfit && )

delete

dayonetrackfit& gar::rec::dayonetrackfit::operator= ( dayonetrackfit const & )

delete

dayonetrackfit& gar::rec::dayonetrackfit::operator= ( dayonetrackfit && )

delete

void gar::rec::dayonetrackfit::produce ( art::Event & e )

overridevirtual

Implements art::EDProducer.

Definition at line 141 of file dayonetrackfit_module.cc.

     {
       // output collections
 
       std::unique_ptr< std::vector<gar::rec::Track> > trkCol(new std::vector<gar::rec::Track>);
       std::unique_ptr< art::Assns<gar::rec::TPCCluster,gar::rec::Track> > TPCClusterTrkAssns(new ::art::Assns<gar::rec::TPCCluster,gar::rec::Track>);
       std::unique_ptr< std::vector<rec::TrackIoniz> > ionCol(new std::vector<rec::TrackIoniz>);
       std::unique_ptr< art::Assns<rec::TrackIoniz,rec::Track> > ionTrkAssns(new ::art::Assns<rec::TrackIoniz,rec::Track>);
       std::unique_ptr< std::vector<rec::TrackTrajectory> > trajCol(new std::vector<rec::TrackTrajectory>);
       std::unique_ptr< art::Assns<rec::TrackTrajectory,rec::Track> > trajTrkAssns(new ::art::Assns<rec::TrackTrajectory,rec::Track>);
 
       // inputs
 
       auto patrecTrackHandle = e.getValidHandle< std::vector<gar::rec::Track> >(fPatRecLabel);
       auto const& patrecTracks = *patrecTrackHandle;
 
       auto const trackPtrMaker = art::PtrMaker<gar::rec::Track>(e);
       auto const ionizPtrMaker = art::PtrMaker<rec::TrackIoniz>(e);
       auto const trajPtrMaker  = art::PtrMaker<rec::TrackTrajectory>(e);
       //auto const TPCClusterPtrMaker = art::PtrMaker<gar::rec::TPCCluster>(e, TPCClusterHandle.id());
 
       art::ServiceHandle<mag::MagneticField> magFieldService;
       G4ThreeVector zerovec(0,0,0);
       G4ThreeVector magfield = magFieldService->FieldAtPoint(zerovec);
 
       art::ServiceHandle<geo::GeometryGAr> geo;
       double xtpccent = geo->TPCXCent();
       double ytpccent = geo->TPCYCent();
       double ztpccent = geo->TPCZCent();
       TVector3 tpccent(xtpccent,ytpccent,ztpccent);
       TVector3 xhat(1,0,0);
 
       const art::FindManyP<gar::rec::TPCCluster> TPCClustersFromPatRecTracks(patrecTrackHandle,e,fPatRecLabel);
 
       for (size_t itrack = 0; itrack < patrecTracks.size(); ++itrack)
         {
           std::vector<gar::rec::TPCCluster> TPCClusters;
           for (size_t iTPCCluster=0; iTPCCluster < TPCClustersFromPatRecTracks.at(itrack).size(); ++iTPCCluster)
             {
               TPCClusters.push_back(*TPCClustersFromPatRecTracks.at(itrack).at(iTPCCluster));  // make our own local copy of TPCClusters.  Maybe we can skip this?
             }
           TrackPar trackparams;
           TrackIoniz trackions;
           TrackTrajectory tracktraj;
           if (KalmanFitBothWays(TPCClusters,trackparams,trackions,tracktraj) == 0)   // to think about -- unused TPCClusters?  Or just ignore them in the fit?
             {
               trkCol->push_back(trackparams.CreateTrack());
               ionCol->push_back(trackions);
               trajCol->push_back(tracktraj);
               auto const trackpointer = trackPtrMaker(trkCol->size()-1);
               auto const ionizpointer = ionizPtrMaker(ionCol->size()-1);
               auto const trajpointer  = trajPtrMaker(trajCol->size()-1);
               for (size_t iTPCCluster=0; iTPCCluster<TPCClusters.size(); ++iTPCCluster)
                 {
                   TPCClusterTrkAssns->addSingle(TPCClustersFromPatRecTracks.at(itrack).at(iTPCCluster),trackpointer);
                 }
               ionTrkAssns->addSingle(ionizpointer, trackpointer);
               trajTrkAssns->addSingle(trajpointer, trackpointer);
             }
         }
 
       e.put(std::move(trkCol));
       e.put(std::move(TPCClusterTrkAssns));
       e.put(std::move(ionCol));
       e.put(std::move(ionTrkAssns));
       e.put(std::move(trajCol));
       e.put(std::move(trajTrkAssns));
     }

void gar::rec::dayonetrackfit::produce ( art::Event & e )

overridevirtual

Implements art::EDProducer.

void gar::rec::dayonetrackfit::updatepar	(	TVectorF &	parvec,
		float	zcur,
		float	xh,
		float	yh,
		float	zh,
		TVectorF &	predstep,
		float &	dlength
	)

private

void gar::rec::dayonetrackfit::updatepar	(	TVectorF &	parvec,
		float	zcur,
		float	xh,
		float	yh,
		float	zh,
		TVectorF &	predstep,
		float &	dlength
	)

private

Definition at line 593 of file dayonetrackfit_module.cc.

      {
 
        predstep = parvec;
 
        // radius and center of the circle
 
        float r = 1E4;
        if (parvec[2] != 0)
          {
            r = 1/parvec[2];
          }
        float yc = parvec[0] + r*TMath::Cos(parvec[3]);
        float zc = zcur - r*TMath::Sin(parvec[3]);
 
        // two solutions of arcsin, plus any number of 2pi.  We could be on
        // the boundary of a wraparound.
 
        float phip1 = TMath::ASin((zh-zc)/r);
        float phip2 = TMath::Pi() - phip1;
 
        // find the solution that's closest to the one suggested by x
 
        float tanlambda = TMath::Tan(parvec[4]);
        if (tanlambda == 0) tanlambda = 1E-4;
        float phip3 = parvec[3] + (xh-parvec[1])/(r*tanlambda);
 
        float dphip1n = (phip1-phip3)/TMath::TwoPi();
        float niphi1 = nearbyint(dphip1n);
        float dphip1a = TMath::Abs(dphip1n - niphi1);
 
        float dphip2n = (phip2-phip3)/TMath::TwoPi();
        float niphi2 = nearbyint(dphip2n);
        float dphip2a = TMath::Abs(dphip2n - niphi2);
 
        float phiex=phip1 + niphi1*TMath::TwoPi();
        if (dphip1a > dphip2a)
          {
            phiex = phip2 + niphi2*TMath::TwoPi(); 
          }
 
        // update predstep to values at this z.  Leave curvature and lambda
        // unchanged
 
        predstep[0] = yc - r*TMath::Cos(phiex);   // y
        predstep[1] = parvec[1] + r*tanlambda*(phiex-parvec[3]);  // x
        predstep[3] = phiex;
       
        dlength = TMath::Abs( r * (phiex - parvec[3]) / TMath::Cos(parvec[4]) );
      }

Member Data Documentation

int gar::rec::dayonetrackfit::fDumpTracks

private

0: do not print out tracks, 1: print out tracks

Definition at line 77 of file dayonetrackfit_module.cc.

unsigned int gar::rec::dayonetrackfit::fInitialTPNTPCClusters

private

number of TPCClusters to use for initial trackpar estimate, if present

Definition at line 75 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fKalCurvStepUncSq

private

constant uncertainty term on each step of the Kalman fit – squared, for curvature

Definition at line 70 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fKalLambdaStepUncSq

private

constant uncertainty term on each step of the Kalman fit – squared, for lambda

Definition at line 72 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fKalPhiStepUncSq

private

constant uncertainty term on each step of the Kalman fit – squared, for phi

Definition at line 71 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fMinIonizGapCut

private

Don't compute dEdx for this dx or larger.

Definition at line 82 of file dayonetrackfit_module.cc.

unsigned int gar::rec::dayonetrackfit::fMinNumTPCClusters

private

minimum number of TPCClusters to define a track

Definition at line 76 of file dayonetrackfit_module.cc.

std::string gar::rec::dayonetrackfit::fPatRecLabel

private

input patrec tracks and associations

Definition at line 68 of file dayonetrackfit_module.cc.

int gar::rec::dayonetrackfit::fPrintLevel

private

debug printout: 0: none, 1: just track parameters and residuals, 2: all

Definition at line 69 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fRoadYZinFit

private

cut in cm for dropping TPCClusters from tracks in fit

Definition at line 79 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fSortDistBack

private

for use in the hit sorting algorithm – how far to go back before raising the distance figure of merit

Definition at line 81 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fSortTransWeight

private

for use in the hit sorting algorithm – transverse distance weight factor

Definition at line 80 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fTPCClusterResolX

private

drift direction contribution to determine step size (resolution of a TPCCluster)

Definition at line 74 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fTPCClusterResolXY

private

pad size in cm in YZ to determine step size

Definition at line 75 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fTPCClusterResolYZ

private

pad size in cm in YZ to determine step size

Definition at line 73 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fTPCClusterResolYZinFit

private

TPCCluster resolution parameter to use in fit.

Definition at line 78 of file dayonetrackfit_module.cc.

float gar::rec::dayonetrackfit::fTPCClusterResolZ

private

drift direction contribution to determine step size (resolution of a TPCCluster)

Definition at line 76 of file dayonetrackfit_module.cc.

The documentation for this class was generated from the following file:

garsoft/Reco/Attic/dayonetrackfit_module.cc

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation