16 #include "Utilities/func/MathUtil.h" 52 for (
size_t i=0;
i< prongs.size();
i++){
110 else if (dirX < 0.0){
140 else if (dirY < 0.0){
155 if (dirX*dirY < 0.0 )
continue;
163 double vx = vert.X();
164 double vy = vert.Y();
165 double vz = vert.Z();
170 if (vert.Z() < 0.0) vz = 2.0;
175 unsigned int vertplane = 999;
178 vertview = (fGeo->
Plane(vertplane))->
View();
193 if ((vertview ==
geo::kX)&&(vertplane != 999)){
198 if (nplane != 999000){
202 double shiftx = (testprong.
Dir().X()/testprong.
Dir().Z())*(shiftz-vz) + vx;
203 double shifty = (testprong.
Dir().Y()/testprong.
Dir().Z())*(shiftz-vz) + vy;
209 if (nplane != 999000){
213 double shiftx = (testprong.
Dir().X()/testprong.
Dir().Z())*(shiftz-vz) + vx;
214 double shifty = (testprong.
Dir().Y()/testprong.
Dir().Z())*(shiftz-vz) + vy;
222 if (nplane != 999000){
226 double shiftx = (testprong.
Dir().X()/testprong.
Dir().Z())*(shiftz-vz) + vx;
227 double shifty = (testprong.
Dir().Y()/testprong.
Dir().Z())*(shiftz-vz) + vy;
232 if (nplane != 999000){
236 double shiftx = (testprong.
Dir().X()/testprong.
Dir().Z())*(shiftz-vz) + vx;
237 double shifty = (testprong.
Dir().Y()/testprong.
Dir().Z())*(shiftz-vz) + vy;
243 if(iyPlane <= fxPlane && fyPlane >= ixPlane){
246 double kScore =
KuiperTest(xlist, ylist, shift, xgev, ygev);
273 if((
fYProngs[j].ExtentPlane()>1) && (
fXProngs[i].ExtentPlane()>1)) score[
i][
j] = kScore;
274 else if((
fYProngs[j].ExtentPlane() == 1) && (ixPlane == (iyPlane-1)) && (fxPlane == (iyPlane+1))) score[
i][
j] = kScore;
275 else if((
fXProngs[i].ExtentPlane() == 1) && (iyPlane == (ixPlane-1)) && (fyPlane == (ixPlane+1))) score[
i][
j] = kScore;
278 else if((ixPlane == fxPlane) && (iyPlane == fyPlane) && (dirX*dirY > 0.0)){
280 if(ixPlane == (iyPlane-1) || ixPlane == (iyPlane+1)){
282 double kScore =
KuiperTest(xlist, ylist, shift, xgev, ygev);
307 score[
i][
j] = kScore;
316 bool doneMatching =
false;
317 std::vector<int> xMatched(
fXProngs.size());
318 for(
unsigned int i = 0;
i<xMatched.size();++
i){
321 std::vector<int> yMatched(
fYProngs.size());
322 for(
unsigned int i = 0;
i<yMatched.size();++
i){
325 while(!doneMatching){
328 double bestscore = -1;
330 for(
unsigned int ix = 0; ix<
fXProngs.size();++ix){
331 for(
unsigned int iy = 0; iy<
fYProngs.size();++iy){
332 if((score[ix][iy] <= bestscore || bestscore == -1) && (score[ix][iy] >= 0.0)){
335 bestscore = score[ix][iy];
346 for(
unsigned int i = 0;
i<
fXProngs[xbest].NCell();++
i){
349 for(
unsigned int i = 0;
i<
fYProngs[ybest].NCell();++
i){
363 double cosThetax = 1/(tanThetax*
sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
364 double cosThetay = 1/(tanThetay*
sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
365 double newz =
sqrt(1-cosThetax*cosThetax-cosThetay*cosThetay);
367 TVector3 newDir(cosThetax,cosThetay,newz);
376 score[xbest][
i] = -1.0;
379 score[
i][ybest] = -1.0;
383 else{ doneMatching =
true; }
386 for (
unsigned int i=0;
i<xMatched.size();
i++){
387 if (xMatched[
i] == 0){
392 for (
unsigned int i=0;
i<yMatched.size();
i++){
393 if (yMatched[
i] == 0){
410 std::vector<std::pair<double,double> >
list;
413 unsigned int ip, ic,
n,
m;
418 double dx[3] = {-0.67, 0.0, 0.67};
419 double dy[3] = {-1.0, 0.0, 1.0};
427 for (
unsigned int ih=0; ih<
ncells; ++ih){
436 double gev = (rhit.
GeV()*testprong.
Weight(v, ih))/9.0;
439 for (n=0; n<3; ++
n) {
440 for (m=0; m<3; ++
m) {
443 fGeo->
ClosestApproach(ip, ic, start, dir, dx[n], dy[m], &w, &s, &pc, &pt);
450 std::sort(list.begin(), list.end());
461 for (
unsigned int k=0; k<xprong.
NCell(); k++) clus.
Add(xprong.
Cell(k), xprong.
Weight(k));
462 for (
unsigned int k=0; k<yprong.
NCell(); k++) clus.
Add(yprong.
Cell(k), yprong.
Weight(k));
465 double tanThetax =
tan(
acos(xprong.
Dir().X()));
466 double tanThetay =
tan(
acos(yprong.
Dir().Y()));
467 double cosThetax = 1.0/(tanThetax*
sqrt(1.0+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
468 double cosThetay = 1.0/(tanThetay*
sqrt(1.0+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
469 double newz =
sqrt(1.0-cosThetax*cosThetax-cosThetay*cosThetay);
470 if ((xprong.
Dir().Z() < 0.0) && (yprong.
Dir().Z() < 0.0)) newz = -newz;
471 newDir->SetXYZ(cosThetax,cosThetay,newz);
476 vert->SetXYZ(xprong.
Start().X(),yprong.
Start().Y(),xprong.
Start().Z());
478 rb::Prong testprong(clus, *vert, *newDir);
488 std::vector<std::pair<double,double>> ylist,
495 if ((xgev <= 0.0)||(ygev <= 0.0)){
508 double ksmaxp = -0.00001;
509 double ksmaxm = -0.00001;
514 double kuiperBest = 2.0;
515 double shiftBest = -13.0;
519 for (
double shiftdist = shift-12.0; shiftdist<shift+12.01; shiftdist+=0.5){
528 while ((ix < xlist.size()) && (iy < ylist.size())) {
529 sx = xlist[ix].first + shiftdist;
530 sy = ylist[iy].first;
532 xsum += xlist[ix].second;
536 ysum += ylist[iy].second;
540 ksdistp = ysum/ygev - xsum/xgev;
541 ksdistm = xsum/xgev - ysum/ygev;
543 if (ksdistp > ksmaxp) ksmaxp = ksdistp;
544 if (ksdistm > ksmaxm) ksmaxm = ksdistm;
548 if (ksmaxp < 0.0) ksmaxp = 0.0;
549 if (ksmaxm < 0.0) ksmaxm = 0.0;
550 kuiper = ksmaxp + ksmaxm;
552 if (kuiper < kuiperBest){
554 shiftBest = shiftdist;
std::vector< std::vector< double > > xzpe
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
void GetCenter(double *xyz, double localz=0.0) const
std::vector< rb::Prong > fMatchedProngs
Containers for all the 3D prongs produced after matching.
fvar< T > fabs(const fvar< T > &x)
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
std::vector< rb::Prong > fXProngs
Containers for all the candidate XZ prongs.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
std::vector< std::vector< double > > xzs
const CellGeo * Cell(int icell) const
std::vector< double > ksscore
Vertical planes which measure X.
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
std::vector< std::vector< double > > yzpe
A collection of associated CellHits.
std::vector< std::vector< rb::PID > > fXScores
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
void Matching()
Function to perform the view matching.
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Defines an enumeration for nerd classification.
virtual void SetStart(TVector3 start)
Horizontal planes which measure Y.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
std::vector< std::pair< double, double > > CalcEnergyProfile(rb::Prong testprong, TVector3 dir, TVector3 start, double *totpecorr, geo::View_t v)
Function to calculate the cumulative energy profile of a test prong.
std::vector< std::vector< rb::PID > > fMatchedYScores
unsigned short Cell() const
std::vector< std::vector< double > > yzs
virtual void SetDir(TVector3 dir)
int getPlaneID(const CellUniqueId &id) const
virtual TVector3 Dir() const
Unit vector describing prong direction.
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
std::vector< double > ksshift
double DetHalfHeight() const
std::vector< std::vector< rb::PID > > fYScores
std::vector< rb::Prong > fUnmatchedProngs
Container for all the unmatched 2D prongs left over.
std::vector< int > xzid
store information related to the view matching for a summary ntuple
std::vector< std::vector< rb::PID > > fMatchedXScores
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
unsigned long long int CellUniqueId
A rawdata::RawDigit with channel information decoded.
A Cluster with defined start position and direction.
double DetHalfWidth() const
double pythag(double x, double y)
2D Euclidean distance
double ClosestApproach(unsigned int ip, unsigned int ic, const TVector3 &vtx, const TVector3 &dir, double dx=0.0, double dy=0.0, double *w=0, double *s=0, TVector3 *point_in_cell=0, TVector3 *point_on_track=0) const
Find the distance of closest approach between a cell and a track.
double KuiperTest(std::vector< std::pair< double, double >> xlist, std::vector< std::pair< double, double >> ylist, double shift, double xpecorr, double ypecorr)
Function to perform the Kuiper test to compare the energy profiles of the track in each view...
std::vector< std::vector< rb::PID > > fUnmatchedScores
Container for all the unmatched 2D scores left over.
rb::Prong MakeTestProng(rb::Prong &xprong, rb::Prong &yprong, TVector3 *newdir, TVector3 *vert)
Function to combine a 2d prong from each view into a 3d test prong to compute a matching score for...
void LoadProngs(std::vector< rb::Prong > prongs, std::vector< std::vector< rb::PID >> scores)
Function to load the vector of 2d prongs to be matched.
double Weight(unsigned int globalIdx) const
Weight assigned to the cell.
const unsigned int NextPlaneOtherView(unsigned int p, int d=+1) const
std::vector< rb::Prong > fYProngs
Containers for all the candidate YZ prongs.
ViewMatchAlg(const ViewMatchParams ¶ms)