590 std::unique_ptr<std::vector<recob::Track> > tcol(
new std::vector<recob::Track>);
594 unsigned int tcnt = 0;
620 for (
unsigned int ii = 0; ii < mctruthListHandle->size(); ++ii)
634 std::vector < art::PtrVector<recob::SpacePoint> > spptIn(spptListHandle->begin(),spptListHandle->end());
659 for(
unsigned int ii = 0; ii < mclist.
size(); ++ii )
663 for(
int jj = 0; jj < mc->NParticles(); ++jj)
669 <<
"FROM MC TRUTH, the particle's pdg code is: "<<
part.PdgCode()<<
" with energy = "<<
part.E() <<
", with energy = "<<
part.E()
672 <<
"\n (both in Global (not volTPC) coords)";
688 LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
" repMC, covMC are ... \n" 694 TParticlePDG *
part = TDatabasePDG::Instance()->GetParticle(
fPdg);
695 Double_t mass = part->Mass();
699 while (sppt!=spptIn.end())
709 unsigned int nTailPoints = 0;
710 if (spacepoints.
size()<5)
711 { sppt++; rePass0 = 3;
continue;}
714 <<
"\n\t found "<<spacepoints.
size()<<
" 3D spacepoint(s) for this element of std::vector<art:PtrVector> spacepoints. \n";
724 TPrincipal* principal =
new TPrincipal(3,
"ND");
732 if (!
fSortDim.compare(
"y")) std::sort(spacepointss.begin(), spacepointss.end(),
sp_sort_3dy);
733 if (!
fSortDim.compare(
"x")) std::sort(spacepointss.begin(), spacepointss.end(),
sp_sort_3dx);
735 for (
unsigned int point=0;point<spacepointss.size();++point)
738 if (point<(spacepointss.size()-nTailPoints))
740 principal->AddRow(spacepointss[point]->XYZ());
743 principal->MakePrincipals();
749 const TVectorD* evals = principal->GetEigenValues();
750 const TMatrixD* evecs = principal->GetEigenVectors();
751 const TVectorD* means = principal->GetMeanValues();
752 const TVectorD* sigmas = principal->GetSigmas();
763 Double_t
tmp[3], tmp2[3];
764 principal->X2P((Double_t *)(means->GetMatrixArray()),tmp);
765 principal->X2P((Double_t *)(sigmas->GetMatrixArray()),tmp2);
766 for (
unsigned int ii=0;ii<3;++ii)
770 fPCevals[ii] = (Float_t )(evals->GetMatrixArray())[ii];
775 evecs->ExtractRow(ii,0,w);
780 Double_t tmp3[3], tmp4[3], tmp5[3];
781 principal->X2P((Double_t *)
fPC1,tmp3);
782 principal->X2P((Double_t *)
fPC2,tmp4);
783 principal->X2P((Double_t *)
fPC3,tmp5);
788 fMomStart[0] = spacepointss[spacepointss.size()-1]->XYZ()[0] - spacepointss[0]->XYZ()[0];
789 fMomStart[1] = spacepointss[spacepointss.size()-1]->XYZ()[1] - spacepointss[0]->XYZ()[1];
790 fMomStart[2] = spacepointss[spacepointss.size()-1]->XYZ()[2] - spacepointss[0]->XYZ()[2];
794 TVector3 mom(dEdx*
fMomStart[0],dEdx*fMomStart[1],dEdx*fMomStart[2]);
795 double pmag2 = pow(mom.Mag()+mass, 2. - mass*mass);
796 mom.SetMag(std::sqrt(pmag2));
798 mom.SetMag(1.0 * mom.Mag());
806 bool uncontained(
false);
808 double epsMag(0.001);
813 spacepointss[spacepointss.size()-1]->XYZ()[0] > (2.*geom->
DetHalfWidth(0,0)-
close) || spacepointss[spacepointss.size()-1]->XYZ()[0] <
close ||
817 spacepointss[spacepointss.size()-1]->XYZ()[2] > (geom->
DetLength(0,0)-
close) || spacepointss[spacepointss.size()-1]->XYZ()[2] <
close ||
818 spacepointss[0]->XYZ()[2] > (geom->
DetLength(0,0)-
close) || spacepointss[0]->XYZ()[2] <
close 829 mom.SetMag(2.0 * mom.Mag());
830 LOG_DEBUG(
"Track3DKalmanSPS_GenFit")<<
"Uncontained track ... ";
836 LOG_DEBUG(
"Track3DKalmanSPS_GenFit")<<
"Contained track ... Run "<<evt.
run()<<
" Event "<<evt.
id().
event();
843 fcont = (int) (!uncontained);
846 unsigned short rePass = rePass0;
848 unsigned short tcnt1(0);
849 while (rePass<=maxPass)
853 TVector3 momErrFit(momM[0]/3.0,
864 (TVector3)(spacepointss[0]->XYZ()),
873 fitTrack.setPDG(
fPdg);
879 std::vector <unsigned int> spptSurvivedIndex;
880 std::vector <unsigned int> spptSkippedIndex;
881 unsigned int ppoint(0);
882 for (
unsigned int point=0;point<spacepointss.size();++point)
889 principal->X2P((Double_t *)(spacepointss[point]->XYZ()),tmp);
891 if ((std::abs(sep) >
fPerpLim) && (point<(spacepointss.size()-nTailPoints)) && rePass<=1)
894 spptSkippedIndex.push_back(point);
899 TVector3 one(spacepointss[point]->XYZ());
900 TVector3 two(spacepointss[ppoint]->XYZ());
901 if (rePass==2 && uncontained)
905 fErrScaleMHere = 0.1;
913 else if (rePass==2 && !uncontained)
921 (one-two).Mag()<epsMag ||
922 ((one-two).Mag()>8.0&&rePass==1) ||
923 std::abs(spacepointss[point]->XYZ()[2]-spacepointss[ppoint]->XYZ()[2])<epsZ ||
924 std::abs(spacepointss[point]->XYZ()[0]-spacepointss[ppoint]->XYZ()[0])>epsX
930 spptSkippedIndex.push_back(point);
934 if (point%fDecimateHere && rePass<=1)
944 TVector3 spt3 = (TVector3)(spacepointss[point]->XYZ());
945 std::vector <double> err3;
946 err3.push_back(spacepointss[point]->ErrXYZ()[0]);
947 err3.push_back(spacepointss[point]->ErrXYZ()[2]);
948 err3.push_back(spacepointss[point]->ErrXYZ()[4]);
949 err3.push_back(spacepointss[point]->ErrXYZ()[5]);
965 fth[
fptsNo] = (pointer.Unit()).Angle(pointerPrev.Unit());
976 LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"ihit xyz..." << spt3[0]<<
","<< spt3[1]<<
","<< spt3[2];
982 spptSurvivedIndex.push_back(point);
988 LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"Bailing cuz only " <<
fptsNo <<
" spacepoints.";
992 LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"Fitting on " <<
fptsNo <<
" spacepoints.";
1005 bool skipFill =
false;
1007 std::vector < TMatrixT<double> > hitMeasCov;
1008 std::vector < TMatrixT<double> > hitUpdate;
1009 std::vector < TMatrixT<double> > hitCov;
1010 std::vector < TMatrixT<double> > hitCov7x7;
1011 std::vector < TMatrixT<double> > hitState;
1012 std::vector < double > hitChi2;
1013 std::vector <TVector3> hitPlaneXYZ;
1014 std::vector <TVector3> hitPlaneUxUyUz;
1015 std::vector <TVector3> hitPlaneU;
1016 std::vector <TVector3> hitPlaneV;
1025 LOG_ERROR(
"Track3DKalmanSPS") <<
"just caught a cet::exception: " << e.what()
1026 <<
"\nExceptions won't be further handled; skip filling big chunks of the TTree.";
1035 std::ostringstream dbgmsg;
1036 dbgmsg <<
"Original plane:";
1037 planeG.Print(dbgmsg);
1039 dbgmsg <<
"Current (fit) reference Plane:";
1042 dbgmsg <<
"Last reference Plane:";
1046 dbgmsg <<
" => original hit plane (not surprisingly) not current reference Plane!";
1048 LOG_DEBUG(
"Track3DKalmanSPS_GenFit") << dbgmsg.str();
1052 hitMeasCov = fitTrack.getHitMeasuredCov();
1053 hitUpdate = fitTrack.getHitUpdate();
1054 hitCov = fitTrack.getHitCov();
1055 hitCov7x7 = fitTrack.getHitCov7x7();
1056 hitState = fitTrack.getHitState();
1057 hitChi2 = fitTrack.getHitChi2();
1058 hitPlaneXYZ = fitTrack.getHitPlaneXYZ();
1059 hitPlaneUxUyUz = fitTrack.getHitPlaneUxUyUz();
1060 hitPlaneU = fitTrack.getHitPlaneU();
1061 hitPlaneV = fitTrack.getHitPlaneV();
1062 unsigned int totHits = hitState.size();
1066 unsigned int jhit=0;
1067 for (
unsigned int ihit=totHits-2*totHits/(2*
fNumIt); ihit<(totHits-totHits/(2*
fNumIt)); ihit++)
1069 feth[jhit] = (Float_t ) (hitMeasCov.at(ihit)[0][0]);
1070 fedudw[jhit] = (Float_t ) (hitMeasCov.at(ihit)[1][1]);
1071 fedvdw[jhit] = (Float_t ) (hitMeasCov.at(ihit)[2][2]);
1072 feu[jhit] = (Float_t ) (hitMeasCov.at(ihit)[3][3]);
1073 fev[jhit] = (Float_t ) (hitMeasCov.at(ihit)[4][4]);
1074 fupdate[jhit] = (Float_t ) (hitUpdate.at(ihit)[0][0]);
1075 fchi2hit[jhit] = (Float_t ) (hitChi2.at(ihit));
1085 for (
unsigned int ii=0;ii<5;ii++)
1087 stREC->ExtractRow(ii,0,dum);
1089 covREC->ExtractRow(ii,0,dum2);
1090 for (
unsigned int jj=0;jj<5;jj++)
1092 fCov0[ii*5+jj] = dum2[jj];
1100 nfail = fitTrack.getFailedHits();
1102 ispptvec=1+std::distance(spptB, sppt);
1106 LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"Track3DKalmanSPS about to do tree->Fill(). Chi2/ndf is " <<
chi2/
ndf <<
".";
1107 fpMCMom[3] = MCMomentum.Mag();
1108 for (
int ii=0;ii<3;++ii)
1112 fpREC[ii] = hitPlaneUxUyUz.at(totHits-2*totHits/(2*
fNumIt))[ii];
1113 fpRECL[ii] = hitPlaneUxUyUz.at(totHits-totHits/(2*
fNumIt)-1)[ii];
1117 nspptvec = (
unsigned int) spptListHandle->size();
1120 std::vector < std::vector <double> > dQdx;
1122 std::vector < TMatrixT<double> > hitCovLFP;
1123 std::vector <TVector3> hitPlaneXYZLFP;
1124 std::vector <TVector3> hitPlaneUxUyUzLFP;
1125 std::vector <TVector3> hitPlaneULFP;
1126 std::vector <TVector3> hitPlaneVLFP;
1127 std::vector <double> pLFP;
1128 std::vector < TMatrixT<double> > c7x7LFP;
1131 for (
unsigned int ii=0; ii<totHits/(2*
fNumIt); ii++)
1133 pLFP.push_back(1./hitState.at(totHits-2*totHits/(2*
fNumIt)+ii)[0][0]);
1135 c7x7LFP.push_back(hitCov7x7.at(totHits-2*totHits/(2*
fNumIt)+ii));
1136 hitCovLFP.push_back(hitCov.at(totHits-2*totHits/(2*
fNumIt)+ii));
1137 hitPlaneXYZLFP.push_back(hitPlaneXYZ.at(totHits-2*totHits/(2*
fNumIt)+ii));
1138 hitPlaneUxUyUzLFP.push_back(hitPlaneUxUyUz.at(totHits-2*totHits/(2*
fNumIt)+ii));
1139 hitPlaneULFP.push_back(hitPlaneU.at(totHits-2*totHits/(2*
fNumIt)+ii));
1140 hitPlaneVLFP.push_back(hitPlaneV.at(totHits-2*totHits/(2*
fNumIt)+ii));
1147 hitPlaneULFP.back(),
1152 hitPlaneUxUyUzLFP.back(),
1153 hitPlaneXYZLFP.back()
1156 fdQdx[ii] = dQdx.back().back();
1168 hitCovLFP,dQdx,pLFP, tcnt++
1170 if (rePass==1) tcnt1++;
1171 if (rePass!=1 && tcnt1) tcol->pop_back();
1172 tcol->push_back(the3DTrack);
1175 for (
unsigned int ii=0; ii < spacepointss.size(); ++ii)
1179 hits.
insert(hits.
end(),hitAssns.at(ii).begin(),hitAssns.at(ii).end());
1191 for (
unsigned int ind=0;ind<spptSurvivedIndex.size();++ind)
1195 (!uncontained&&
fchi2hit[ind]>1.e9) ||
1209 for (
unsigned int ind=0;ind<spptSkippedIndex.size();++ind)
1216 std::stable_sort(spacepointss.begin(),spacepointss.end());
1217 std::stable_sort(spacepointssExcise.
begin(),spacepointssExcise.
end());
1218 std::set_union(spacepointssExcise.
begin(),spacepointssExcise.
end(),
1219 spacepointssExcise.
begin(),spacepointssExcise.
end(),
1220 spacepointssExcise.
begin()
1224 std::set_difference(spacepointss.begin(),spacepointss.end(),
1225 spacepointssExcise.
begin(),spacepointssExcise.
end(),
1226 spacepointss.begin()
1228 spacepointss.erase(diffSpptIt,spacepointss.end());
1241 if (uncontained) kick = 0.5;
1242 for (
int ii=0;ii<3;++ii)
1245 mom[ii] = momM[ii]*kick;
1248 else if (uncontained)
1250 double unstick(1.0);
1252 for (
int ii=0;ii<3;++ii)
1254 mom[ii] = momM[ii]*unstick;
1258 for (
int ii=0;ii<3;++ii)
1260 mom[ii] = 1.1*momM[ii];
1273 evt.
put(std::move(tcol));
1275 evt.
put(std::move(tspassn));
1278 evt.
put(std::move(thassn));
TMatrixT< Double_t > * stREC
unsigned int getNDF() const
std::vector< double > fMomErr
void setMomLow(Double_t f)
void setMaxUpdate(Double_t f)
TMatrixT< Double_t > * stMCT
static bool sp_sort_3dz(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
const GFDetPlane & getReferencePlane() const
GFDetPlane getLastPlane() const
void rotationCov(TMatrixT< Double_t > &cov, const TVector3 &u, const TVector3 &v)
genf::GFAbsTrackRep * rep
void setMomHigh(Double_t f)
std::string fGenieGenModuleLabel
std::vector< double > fPosErr
void Print(std::ostream &out=std::cout) const
const TMatrixT< Double_t > & getState() const
std::string fSpptModuleLabel
static bool sp_sort_3dy(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
std::string fClusterModuleLabel
#define LOG_ERROR(category)
ProductID put(std::unique_ptr< PROD > &&product)
geo::Length_t DetHalfHeight(geo::TPCID const &tpcid) const
Returns the half height of the active volume of the specified TPC.
std::vector< double > fMomStart
void push_back(Ptr< U > const &p)
virtual TVector3 getMom(const GFDetPlane &pl)=0
static bool sp_sort_3dx(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
void setBlowUpFactor(double f)
Set the blowup factor (see blowUpCovs() )
static bool sp_sort_nsppts(const art::PtrVector< recob::SpacePoint > &h1, const art::PtrVector< recob::SpacePoint > &h2)
double energyLossBetheBloch(const double &mass, const double p)
void setNumIterations(Int_t i)
Set number of iterations for Kalman Filter.
TMatrixT< Double_t > * covREC
void processTrack(GFTrack *)
Performs fit on a GFTrack.
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
static GFFieldManager * getInstance()
TMatrixT< Double_t > * covMCT
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
const TMatrixT< Double_t > & getCov() const
data_t::iterator iterator
std::string ROOTobjectToString(const ROOTOBJ &obj)
Shortcut to write one ROOT object into a string.
iterator insert(iterator position, Ptr< U > const &p)
std::vector< double > dQdxCalc(const art::FindManyP< recob::Hit > &h, const art::PtrVector< recob::SpacePoint > &s, const TVector3 &p, const TVector3 &d)
genf::GFAbsTrackRep * repMC
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
void setInitialDirection(int d)
Sets the inital direction of the track fit (1 for inner to outer, or -1 for outer to inner)...
void init(GFAbsBField *b)
set the magntic field here. Magnetic field classes must be derived from GFAbsBField ...
EventNumber_t event() const
void setErrorScaleSTh(Double_t f)
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
void setErrorScaleMTh(Double_t f)
cet::coded_exception< error, detail::translate > exception
std::string GetLArTPCVolumeName(geo::TPCID const &tpcid) const
Return the name of specified LAr TPC volume.