559 art::ServiceHandle<geo::Geometry const> geom;
565 std::unique_ptr<std::vector<recob::Track> > tcol(
new std::vector<recob::Track>);
566 std::unique_ptr< art::Assns<recob::Track, recob::SpacePoint> > tspassn(
new art::Assns<recob::Track, recob::SpacePoint>);
567 std::unique_ptr< art::Assns<recob::Track, recob::Hit> > thassn(
new art::Assns<recob::Track, recob::Hit>);
569 unsigned int tcnt = 0;
572 TString tpcName = geom->GetLArTPCVolumeName();
576 art::Handle< std::vector<recob::Cluster> > clusterListHandle;
579 art::Handle< std::vector< art::PtrVector < recob::SpacePoint > > > spptListHandle;
582 art::PtrVector<simb::MCTruth> mclist;
587 if (!
evt.isRealData()){
592 art::Handle< std::vector<simb::MCTruth> > mctruthListHandle;
595 for (
unsigned int ii = 0; ii < mctruthListHandle->size(); ++ii)
597 art::Ptr<simb::MCTruth> mctparticle(mctruthListHandle,ii);
598 mclist.push_back(mctparticle);
609 std::vector < art::PtrVector<recob::SpacePoint> > spptIn(spptListHandle->begin(),spptListHandle->end());
613 std::vector < art::PtrVector<recob::SpacePoint> >::const_iterator sppt = spptIn.begin();
629 if (!
evt.isRealData())
634 for(
unsigned int ii = 0; ii < mclist.size(); ++ii )
637 art::Ptr<simb::MCTruth> mc(mclist[ii]);
638 for(
int jj = 0; jj < mc->NParticles(); ++jj)
640 simb::MCParticle part(mc->GetParticle(jj));
641 MCOrigin.SetXYZ(part.Vx(),part.Vy(),part.Vz());
642 MCMomentum.SetXYZ(part.Px(),part.Py(),part.Pz());
643 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit")
644 <<
"FROM MC TRUTH, the particle's pdg code is: "<<part.PdgCode()<<
" with energy = "<<part.E() <<
", with energy = "<<part.E()
647 <<
"\n (both in Global (not volTPC) coords)";
663 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
" repMC, covMC are ... \n"
669 TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(
fPdg);
670 Double_t
mass = part->Mass();
674 while (sppt!=spptIn.end())
677 const art::PtrVector<recob::SpacePoint> & spacepoints = *sppt;
684 unsigned int nTailPoints = 0;
685 if (spacepoints.size()<5)
686 { sppt++; rePass0 = 3;
continue;}
688 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit")
689 <<
"\n\t found "<<spacepoints.size()<<
" 3D spacepoint(s) for this element of std::vector<art:PtrVector> spacepoints. \n";
699 TPrincipal* principal =
new TPrincipal(3,
"ND");
703 art::PtrVector<recob::SpacePoint> spacepointss(spacepoints);
707 if (!
fSortDim.compare(
"y")) std::sort(spacepointss.begin(), spacepointss.end(),
sp_sort_3dy);
708 if (!
fSortDim.compare(
"x")) std::sort(spacepointss.begin(), spacepointss.end(),
sp_sort_3dx);
710 for (
unsigned int point=0;point<spacepointss.size();++point)
713 if (point<(spacepointss.size()-nTailPoints))
715 principal->AddRow(spacepointss[point]->XYZ());
718 principal->MakePrincipals();
724 const TVectorD* evals = principal->GetEigenValues();
725 const TMatrixD* evecs = principal->GetEigenVectors();
726 const TVectorD* means = principal->GetMeanValues();
727 const TVectorD* sigmas = principal->GetSigmas();
738 Double_t tmp[3], tmp2[3];
739 principal->X2P((Double_t *)(means->GetMatrixArray()),tmp);
740 principal->X2P((Double_t *)(sigmas->GetMatrixArray()),tmp2);
741 for (
unsigned int ii=0;ii<3;++ii)
745 fPCevals[ii] = (Float_t )(evals->GetMatrixArray())[ii];
750 evecs->ExtractRow(ii,0,w);
755 Double_t tmp3[3], tmp4[3], tmp5[3];
756 principal->X2P((Double_t *)
fPC1,tmp3);
757 principal->X2P((Double_t *)
fPC2,tmp4);
758 principal->X2P((Double_t *)
fPC3,tmp5);
763 fMomStart[0] = spacepointss[spacepointss.size()-1]->XYZ()[0] - spacepointss[0]->XYZ()[0];
764 fMomStart[1] = spacepointss[spacepointss.size()-1]->XYZ()[1] - spacepointss[0]->XYZ()[1];
765 fMomStart[2] = spacepointss[spacepointss.size()-1]->XYZ()[2] - spacepointss[0]->XYZ()[2];
769 TVector3 mom(dEdx*
fMomStart[0],dEdx*fMomStart[1],dEdx*fMomStart[2]);
770 double pmag2 = pow(mom.Mag()+
mass, 2. - mass*
mass);
771 mom.SetMag(std::sqrt(pmag2));
773 mom.SetMag(1.0 * mom.Mag());
781 bool uncontained(
false);
783 double epsMag(0.001);
788 spacepointss[spacepointss.size()-1]->XYZ()[0] > (2.*geom->DetHalfWidth(0,0)-
close) || spacepointss[spacepointss.size()-1]->XYZ()[0] <
close ||
789 spacepointss[0]->XYZ()[0] > (2.*geom->DetHalfWidth(0,0)-
close) || spacepointss[0]->XYZ()[0] <
close ||
790 spacepointss[spacepointss.size()-1]->XYZ()[1] > (1.*geom->DetHalfHeight(0,0)-
close) || (spacepointss[spacepointss.size()-1]->XYZ()[1] < -1.*geom->DetHalfHeight(0,0)+
close) ||
791 spacepointss[0]->XYZ()[1] > (1.*geom->DetHalfHeight(0,0)-
close) || spacepointss[0]->XYZ()[1] < (-1.*geom->DetHalfHeight(0,0)+
close) ||
792 spacepointss[spacepointss.size()-1]->XYZ()[2] > (geom->DetLength(0,0)-
close) || spacepointss[spacepointss.size()-1]->XYZ()[2] <
close ||
793 spacepointss[0]->XYZ()[2] > (geom->DetLength(0,0)-
close) || spacepointss[0]->XYZ()[2] <
close
804 mom.SetMag(2.0 * mom.Mag());
805 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit")<<
"Uncontained track ... ";
811 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit")<<
"Contained track ... Run "<<
evt.
run()<<
" Event "<<
evt.id().
event();
818 fcont = (int) (!uncontained);
821 unsigned short rePass = rePass0;
823 unsigned short tcnt1(0);
824 while (rePass<=maxPass)
828 TVector3 momErrFit(momM[0]/3.0,
839 (TVector3)(spacepointss[0]->XYZ()),
848 fitTrack.setPDG(
fPdg);
854 std::vector <unsigned int> spptSurvivedIndex;
855 std::vector <unsigned int> spptSkippedIndex;
856 unsigned int ppoint(0);
857 for (
unsigned int point=0;point<spacepointss.size();++point)
864 principal->X2P((Double_t *)(spacepointss[point]->XYZ()),tmp);
866 if ((
std::abs(sep) >
fPerpLim) && (point<(spacepointss.size()-nTailPoints)) && rePass<=1)
869 spptSkippedIndex.push_back(point);
874 TVector3
one(spacepointss[point]->XYZ());
875 TVector3 two(spacepointss[ppoint]->XYZ());
876 if (rePass==2 && uncontained)
880 fErrScaleMHere = 0.1;
888 else if (rePass==2 && !uncontained)
896 (
one-two).Mag()<epsMag ||
897 ((
one-two).Mag()>8.0&&rePass==1) ||
898 std::abs(spacepointss[point]->XYZ()[2]-spacepointss[ppoint]->XYZ()[2])<epsZ ||
899 std::abs(spacepointss[point]->XYZ()[0]-spacepointss[ppoint]->XYZ()[0])>epsX
905 spptSkippedIndex.push_back(point);
909 if (point%fDecimateHere && rePass<=1)
919 TVector3 spt3 = (TVector3)(spacepointss[point]->XYZ());
920 std::vector <double> err3;
921 err3.push_back(spacepointss[point]->ErrXYZ()[0]);
922 err3.push_back(spacepointss[point]->ErrXYZ()[2]);
923 err3.push_back(spacepointss[point]->ErrXYZ()[4]);
924 err3.push_back(spacepointss[point]->ErrXYZ()[5]);
940 fth[
fptsNo] = (pointer.Unit()).Angle(pointerPrev.Unit());
951 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"ihit xyz..." << spt3[0]<<
","<< spt3[1]<<
","<< spt3[2];
957 spptSurvivedIndex.push_back(point);
963 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"Bailing cuz only " <<
fptsNo <<
" spacepoints.";
967 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"Fitting on " <<
fptsNo <<
" spacepoints.";
980 bool skipFill =
false;
982 std::vector < TMatrixT<double> > hitMeasCov;
983 std::vector < TMatrixT<double> > hitUpdate;
984 std::vector < TMatrixT<double> > hitCov;
985 std::vector < TMatrixT<double> > hitCov7x7;
986 std::vector < TMatrixT<double> > hitState;
987 std::vector < double > hitChi2;
988 std::vector <TVector3> hitPlaneXYZ;
989 std::vector <TVector3> hitPlaneUxUyUz;
990 std::vector <TVector3> hitPlaneU;
991 std::vector <TVector3> hitPlaneV;
999 catch(cet::exception &
e){
1000 MF_LOG_ERROR(
"Track3DKalmanSPS") <<
"just caught a cet::exception: " << e.what()
1001 <<
"\nExceptions won't be further handled; skip filling big chunks of the TTree.";
1008 if(mf::isDebugEnabled()) {
1010 std::ostringstream dbgmsg;
1011 dbgmsg <<
"Original plane:";
1012 planeG.Print(dbgmsg);
1014 dbgmsg <<
"Current (fit) reference Plane:";
1017 dbgmsg <<
"Last reference Plane:";
1021 dbgmsg <<
" => original hit plane (not surprisingly) not current reference Plane!";
1023 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") << dbgmsg.str();
1027 hitMeasCov = fitTrack.getHitMeasuredCov();
1028 hitUpdate = fitTrack.getHitUpdate();
1029 hitCov = fitTrack.getHitCov();
1030 hitCov7x7 = fitTrack.getHitCov7x7();
1031 hitState = fitTrack.getHitState();
1032 hitChi2 = fitTrack.getHitChi2();
1033 hitPlaneXYZ = fitTrack.getHitPlaneXYZ();
1034 hitPlaneUxUyUz = fitTrack.getHitPlaneUxUyUz();
1035 hitPlaneU = fitTrack.getHitPlaneU();
1036 hitPlaneV = fitTrack.getHitPlaneV();
1037 unsigned int totHits = hitState.size();
1041 unsigned int jhit=0;
1042 for (
unsigned int ihit=totHits-2*totHits/(2*
fNumIt); ihit<(totHits-totHits/(2*
fNumIt)); ihit++)
1044 feth[jhit] = (Float_t ) (hitMeasCov.at(ihit)[0][0]);
1045 fedudw[jhit] = (Float_t ) (hitMeasCov.at(ihit)[1][1]);
1046 fedvdw[jhit] = (Float_t ) (hitMeasCov.at(ihit)[2][2]);
1047 feu[jhit] = (Float_t ) (hitMeasCov.at(ihit)[3][3]);
1048 fev[jhit] = (Float_t ) (hitMeasCov.at(ihit)[4][4]);
1049 fupdate[jhit] = (Float_t ) (hitUpdate.at(ihit)[0][0]);
1050 fchi2hit[jhit] = (Float_t ) (hitChi2.at(ihit));
1060 for (
unsigned int ii=0;ii<5;ii++)
1062 stREC->ExtractRow(ii,0,dum);
1064 covREC->ExtractRow(ii,0,dum2);
1065 for (
unsigned int jj=0;jj<5;jj++)
1067 fCov0[ii*5+jj] = dum2[jj];
1070 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit")
1075 nfail = fitTrack.getFailedHits();
1081 MF_LOG_DEBUG(
"Track3DKalmanSPS_GenFit") <<
"Track3DKalmanSPS about to do tree->Fill(). Chi2/ndf is " <<
chi2/
ndf <<
".";
1082 fpMCMom[3] = MCMomentum.Mag();
1083 for (
int ii=0;ii<3;++ii)
1087 fpREC[ii] = hitPlaneUxUyUz.at(totHits-2*totHits/(2*
fNumIt))[ii];
1088 fpRECL[ii] = hitPlaneUxUyUz.at(totHits-totHits/(2*
fNumIt)-1)[ii];
1092 nspptvec = (
unsigned int) spptListHandle->size();
1095 std::vector < std::vector <double> > dQdx;
1097 std::vector < TMatrixT<double> > hitCovLFP;
1098 std::vector <TVector3> hitPlaneXYZLFP;
1099 std::vector <TVector3> hitPlaneUxUyUzLFP;
1100 std::vector <TVector3> hitPlanePxPyPzLFP;
1101 std::vector <TVector3> hitPlaneULFP;
1102 std::vector <TVector3> hitPlaneVLFP;
1103 std::vector <double> pLFP;
1104 std::vector < TMatrixT<double> > c7x7LFP;
1107 for (
unsigned int ii=0; ii<totHits/(2*
fNumIt); ii++)
1109 pLFP.push_back(1./hitState.at(totHits-2*totHits/(2*
fNumIt)+ii)[0][0]);
1111 c7x7LFP.push_back(hitCov7x7.at(totHits-2*totHits/(2*
fNumIt)+ii));
1112 hitCovLFP.push_back(hitCov.at(totHits-2*totHits/(2*
fNumIt)+ii));
1113 hitPlaneXYZLFP.push_back(hitPlaneXYZ.at(totHits-2*totHits/(2*
fNumIt)+ii));
1114 hitPlaneUxUyUzLFP.push_back(hitPlaneUxUyUz.at(totHits-2*totHits/(2*
fNumIt)+ii));
1115 hitPlanePxPyPzLFP.push_back(hitPlaneUxUyUzLFP.back()*pLFP.back());
1116 hitPlaneULFP.push_back(hitPlaneU.at(totHits-2*totHits/(2*
fNumIt)+ii));
1117 hitPlaneVLFP.push_back(hitPlaneV.at(totHits-2*totHits/(2*
fNumIt)+ii));
1124 hitPlaneULFP.back(),
1129 hitPlaneUxUyUzLFP.back(),
1130 hitPlaneXYZLFP.back()
1133 fdQdx[ii] = dQdx.back().back();
1144 for (
unsigned int i=0; i<5; i++) {
1145 for (
unsigned int j=i; j<5; j++) {
1146 covVtx(i,j) = hitCovLFP.front()(i,j);
1147 covEnd(i,j) = hitCovLFP.back()(i,j);
1153 0, -1., 0, covVtx, covEnd, tcnt++);
1154 if (rePass==1) tcnt1++;
1155 if (rePass!=1 && tcnt1) tcol->pop_back();
1156 tcol->push_back(the3DTrack);
1158 art::PtrVector<recob::Hit> hits;
1159 for (
unsigned int ii=0; ii < spacepointss.size(); ++ii)
1163 hits.insert(hits.end(),hitAssns.at(ii).begin(),hitAssns.at(ii).end());
1174 art::PtrVector<recob::SpacePoint> spacepointssExcise;
1175 for (
unsigned int ind=0;ind<spptSurvivedIndex.size();++ind)
1179 (!uncontained&&
fchi2hit[ind]>1.e9) ||
1186 art::PtrVector<recob::SpacePoint>::iterator spptIt = spacepointss.begin()+spptSurvivedIndex[ind];
1187 spacepointssExcise.push_back(*spptIt);
1193 for (
unsigned int ind=0;ind<spptSkippedIndex.size();++ind)
1195 art::PtrVector<recob::SpacePoint>::iterator spptIt = spacepointss.begin()+spptSkippedIndex[ind];
1196 spacepointssExcise.push_back(*spptIt);
1200 std::stable_sort(spacepointss.begin(),spacepointss.end());
1201 std::stable_sort(spacepointssExcise.begin(),spacepointssExcise.end());
1202 std::set_union(spacepointssExcise.begin(),spacepointssExcise.end(),
1203 spacepointssExcise.begin(),spacepointssExcise.end(),
1204 spacepointssExcise.begin()
1207 art::PtrVector<recob::SpacePoint>::iterator diffSpptIt =
1208 std::set_difference(spacepointss.begin(),spacepointss.end(),
1209 spacepointssExcise.begin(),spacepointssExcise.end(),
1210 spacepointss.begin()
1212 spacepointss.erase(diffSpptIt,spacepointss.end());
1225 if (uncontained) kick = 0.5;
1226 for (
int ii=0;ii<3;++ii)
1229 mom[ii] = momM[ii]*kick;
1232 else if (uncontained)
1234 double unstick(1.0);
1236 for (
int ii=0;ii<3;++ii)
1238 mom[ii] = momM[ii]*unstick;
1242 for (
int ii=0;ii<3;++ii)
1244 mom[ii] = 1.1*momM[ii];
1257 evt.put(std::move(tcol));
1259 evt.put(std::move(tspassn));
1262 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)
const GFDetPlane & getReferencePlane() const
GFDetPlane getLastPlane() const
void rotationCov(TMatrixT< Double_t > &cov, const TVector3 &u, const TVector3 &v)
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 > > SMatrixSym55
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
genf::GFAbsTrackRep * rep
void setMomHigh(Double_t f)
TrackTrajectory::Flags_t Flags_t
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
std::vector< double > fMomStart
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
virtual TVector3 getMom(const GFDetPlane &pl)=0
static bool sp_sort_3dx(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
A trajectory in space reconstructed from hits.
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.
static GFFieldManager * getInstance()
TMatrixT< Double_t > * covMCT
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
const TMatrixT< Double_t > & getCov() const
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
std::string ROOTobjectToString(const ROOTOBJ &obj)
Shortcut to write one ROOT object into a string.
then echo Cowardly refusing to create a new FHiCL file with the same name as the original one('${SourceName}')." >&2 exit 1 fi echo "'$
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
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
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 ...
void setErrorScaleSTh(Double_t f)
print OUTPUT<< EOF;< setup name="Default"version="1.0">< worldref="volWorld"/></setup ></gdml > EOF close(OUTPUT)
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)