23 #include "Utilities/AssociationUtil.h" 42 static bool sort_traj(
const TVector3&
a,
const TVector3&
b);
46 std::vector<rb::Track> xtracks,
47 std::vector<rb::Track> ytracks,
48 std::vector<rb::Track>&matchedtracls,
49 std::vector<rb::Track>&unmatchedtracks);
56 TVector3 &lowztraj, TVector3 &highztraj);
78 produces< std::vector<rb::Track> >();
79 produces< art::Assns<rb::Track, rb::Cluster> >();
104 std::unique_ptr< std::vector<rb::Track> > mergetrackcol(
new std::vector<rb::Track>);
114 for(
unsigned int i = 0;
i<slicecol->size();++
i){
123 for(
unsigned int iSlice = 0; iSlice<slicelist.
size(); ++iSlice){
125 if(slicelist[iSlice]->IsNoise()){
continue; }
127 const std::vector< art::Ptr<rb::Track> > slTracks = fmt.at(iSlice);
129 if(slTracks.empty()){
continue; }
134 avgX = slicelist[iSlice]->MeanX();
135 avgY = slicelist[iSlice]->MeanY();
136 std::vector<trk::LocatedChan> xcellchans;
137 double zxmin = 100000;
139 for(
unsigned int nx = 0; nx < slicelist[iSlice]->NXCell(); ++nx){
146 xchan.SetCenter(xyz[0],xyz[1],xyz[2]);
148 xcellchans.push_back(xchan);
149 if(xyz[2]>zxmax){ zxmax = xyz[2]; }
150 if(xyz[2]<zxmin){ zxmin = xyz[2]; }
155 int nxzbreak = (
int)
floor((zxmax-zxmin)/100);
157 double dz = (zxmax-zxmin)/((
double)nxzbreak);
158 std::vector<double> xzbreaks;
159 xzbreaks.push_back(zxmin);
160 for(
int ibrk = 1; ibrk<=nxzbreak; ++ibrk){
161 xzbreaks.push_back(zxmin+((
double)ibrk)*dz);
163 xzbreaks.push_back(zxmax);
164 for(
unsigned int ibrk = 1; ibrk < xzbreaks.size(); ++ibrk){
167 for(
unsigned int xchan = 0;xchan<xcellchans.size();++xchan){
169 xcellchans[xchan].GetCenter(xyz);
170 if(xyz[2] >= xzbreaks[ibrk-1] && xyz[2] <=xzbreaks[ibrk]){
175 std::vector<double> zrng;
176 zrng.push_back(xzbreaks[ibrk-1]);
177 zrng.push_back(xzbreaks[ibrk]);
179 if(ntot>0){
avgXPos.push_back(xtot/ntot); }
185 std::vector<trk::LocatedChan> ycellchans;
186 double zymin = 100000;
188 for(
unsigned int ny = 0; ny < slicelist[iSlice]->NYCell(); ++ny){
195 ychan.SetCenter(xyz[0],xyz[1],xyz[2]);
197 ycellchans.push_back(ychan);
198 if(xyz[2]>zymax){ zymax = xyz[2]; }
199 if(xyz[2]<zymin){ zymin = xyz[2]; }
203 int nyzbreak = (
int)
floor((zymax-zymin)/100);
205 dz = (zymax-zymin)/((
double)nyzbreak);
206 std::vector<double> yzbreaks;
207 yzbreaks.push_back(zymin);
208 for(
int ibrk = 1; ibrk<=nyzbreak; ++ibrk){
209 yzbreaks.push_back(zymin+((
double)ibrk)*dz);
211 yzbreaks.push_back(zymax);
212 for(
unsigned int ibrk = 1; ibrk < yzbreaks.size(); ++ibrk){
215 for(
unsigned int ychan = 0;ychan<ycellchans.size();++ychan){
217 ycellchans[ychan].GetCenter(xyz);
218 if(xyz[2] >= yzbreaks[ibrk-1] && xyz[2] <=yzbreaks[ibrk]){
223 std::vector<double> zrng;
224 zrng.push_back(yzbreaks[ibrk-1]);
225 zrng.push_back(yzbreaks[ibrk]);
227 if(ntot>0){
avgYPos.push_back(ytot/ntot); }
232 std::vector<rb::Track> xTracks;
233 std::vector<rb::Track> yTracks;
234 for(
unsigned int j = 0;
j<slTracks.size();++
j){
236 xTracks.push_back(*slTracks[
j]);
239 yTracks.push_back(*slTracks[
j]);
244 std::vector<rb::Track> matchedTracks;
245 std::vector<rb::Track> unmatchedTracks;
246 MatchTracks(kgeo, xTracks,yTracks,matchedTracks,unmatchedTracks);
247 for(
unsigned int j = 0;
j<matchedTracks.size();
j++){
248 mergetrackcol->push_back(matchedTracks[
j]);
251 for(
unsigned int j = 0;
j<unmatchedTracks.size();
j++){
252 mergetrackcol->push_back(unmatchedTracks[
j]);
258 evt.
put(std::move(mergetrackcol));
259 evt.
put(std::move(assns));
266 std::vector<rb::Track > xTracks, std::vector<rb::Track> yTracks,
267 std::vector<rb::Track> &matchedtracks, std::vector<rb::Track> &unmatchedtracks)
270 std::vector<rb::Track> allXTracks;
271 std::vector<std::vector< int > > jointxtrkIDs;
274 for(
unsigned int ix = 0; ix<xTracks.size();++ix){
275 allXTracks.push_back(xTracks[ix]);
276 std::vector<int> xid;
278 jointxtrkIDs.push_back(xid);
281 unsigned int ixmaxold = 0;
282 unsigned int ixmax = allXTracks.size();
283 unsigned int ixmin = 0;
284 while(ixmaxold != ixmax && allXTracks.size()<100){
286 for(
unsigned int ix = 0; ix<ixmaxold; ++ix){
288 std::vector<int> ixids = jointxtrkIDs[ix];
289 std::vector<int> jntids = ixids;
290 for(
unsigned int jx = ixmin; jx<ixmaxold; ++jx){
292 std::vector<int> jxids = jointxtrkIDs[jx];
294 bool idoverlap =
false;
295 for(
unsigned int i = 0;
i<ixids.size(); ++
i){
296 for(
unsigned int j = 0;
j<jxids.size(); ++
j){
297 if(ixids[
i] == jxids[
j]){
302 if(idoverlap){
break; }
304 if(idoverlap){
continue; }
306 for(
unsigned int j = 0;
j<jxids.size(); ++
j){
307 jntids.push_back(jxids[
j]);
309 std::sort(jntids.begin(),jntids.end());
311 for(
unsigned int i = 0;
i<jointxtrkIDs.size(); ++
i){
312 if(jntids.size() != jointxtrkIDs[
i].size()){
continue; }
314 std::vector<int> jntidscheck = jointxtrkIDs[
i];
315 std::sort(jntidscheck.begin(),jntidscheck.end());
316 for(
unsigned int j = 0;
j<jointxtrkIDs[
i].size(); ++
j){
317 if(jntids[
j] != jntidscheck[
j]){
327 if(used){
continue; }
328 bool join =
CanJoinTracks(kgeo,allXTracks[ix],allXTracks[jx]);
331 allXTracks.push_back(JointTrack);
332 std::vector<int> xjointid;
333 for(
unsigned int i = 0;
i<ixids.size(); ++
i){
334 xjointid.push_back(ixids[
i]);
336 for(
unsigned int j = 0;
j<jxids.size(); ++
j){
337 xjointid.push_back(jxids[
j]);
339 jointxtrkIDs.push_back(xjointid);
342 ixmax = allXTracks.size();
344 ixmax = allXTracks.size();
348 std::vector<rb::Track> allYTracks;
349 std::vector<std::vector< int > > jointytrkIDs;
352 for(
unsigned int iy = 0; iy<yTracks.size();++iy){
353 allYTracks.push_back(yTracks[iy]);
354 std::vector<int> yid;
356 jointytrkIDs.push_back(yid);
359 unsigned int iymaxold = 0;
360 unsigned int iymax = allYTracks.size();
361 unsigned int iymin = 0;
362 while(iymaxold != iymax && allYTracks.size()<100){
364 for(
unsigned int iy = 0; iy<iymaxold; ++iy){
366 std::vector<int> iyids = jointytrkIDs[iy];
367 std::vector<int> jntids = iyids;
368 for(
unsigned int jy = iymin; jy<iymaxold; ++jy){
370 std::vector<int> jyids = jointytrkIDs[jy];
372 bool idoverlap =
false;
373 for(
unsigned int i = 0;
i<iyids.size(); ++
i){
374 for(
unsigned int j = 0;
j<jyids.size(); ++
j){
375 if(iyids[
i] == jyids[
j]){
380 if(idoverlap){
break; }
382 if(idoverlap){
continue; }
384 for(
unsigned int j = 0;
j<jyids.size(); ++
j){
385 jntids.push_back(jyids[
j]);
387 std::sort(jntids.begin(),jntids.end());
389 for(
unsigned int i = 0;
i<jointytrkIDs.size(); ++
i){
390 if(jntids.size() != jointytrkIDs[
i].size()){
continue; }
392 std::vector<int> jntidscheck = jointytrkIDs[
i];
393 std::sort(jntidscheck.begin(),jntidscheck.end());
394 for(
unsigned int j = 0;
j<jointytrkIDs[
i].size(); ++
j){
395 if(jntids[
j] != jntidscheck[
j]){
405 if(used){
continue; }
406 bool join =
CanJoinTracks(kgeo,allYTracks[iy],allYTracks[jy]);
409 allYTracks.push_back(JointTrack);
410 std::vector<int> yjointid;
411 for(
unsigned int i = 0;
i<iyids.size(); ++
i){
412 yjointid.push_back(iyids[
i]);
414 for(
unsigned int j = 0;
j<jyids.size(); ++
j){
415 yjointid.push_back(jyids[
j]);
417 jointytrkIDs.push_back(yjointid);
420 iymax = allYTracks.size();
422 iymax = allYTracks.size();
423 iymin = iymaxold - 1;
428 double score[allXTracks.size()][allYTracks.size()];
429 for(
unsigned int ix = 0; ix<allXTracks.size();++ix){
431 for(
unsigned int iy = 0; iy<allYTracks.size();++iy){
432 score[ix][iy] =
CalcMatchScore(kgeo, allXTracks[ix],allYTracks[iy]);
438 bool doneMatching =
false;
439 std::vector<int> xMatched(allXTracks.size());
440 for(
unsigned int i = 0;
i<xMatched.size();++
i){
443 std::vector<int> yMatched(allYTracks.size());
444 for(
unsigned int i = 0;
i<yMatched.size();++
i){
447 while(!doneMatching){
450 double bestscore = -1;
452 for(
unsigned int ix = 0; ix<allXTracks.size();++ix){
453 for(
unsigned int iy = 0; iy<allYTracks.size();++iy){
454 if((score[ix][iy] <= bestscore || bestscore == -1) && score[ix][iy] >= 0.0 && score[ix][iy] <= 5.0){
457 bestscore = score[ix][iy];
468 matchedtracks.push_back(newtrack);
471 for(
unsigned int i = 0;
i<jointxtrkIDs[xbest].size(); ++
i){
472 int matchedtrkid = jointxtrkIDs[xbest][
i];
474 for(
unsigned int itrk = 0; itrk<allXTracks.size(); ++itrk){
476 if(xMatched[itrk] == 0){
478 bool usedtrk =
false;
479 for(
unsigned int itrkid = 0; itrkid<jointxtrkIDs[itrk].size(); ++itrkid){
480 if(jointxtrkIDs[itrk][itrkid] == matchedtrkid){
488 for(
unsigned int iytrk = 0; iytrk<allYTracks.size();++iytrk){
489 score[itrk][iytrk] = -1.0;
497 for(
unsigned int i = 0;
i<jointytrkIDs[ybest].size(); ++
i){
498 int matchedtrkid = jointytrkIDs[ybest][
i];
500 for(
unsigned int itrk = 0; itrk<allYTracks.size(); ++itrk){
502 if(yMatched[itrk] == 0){
504 bool usedtrk =
false;
505 for(
unsigned int itrkid = 0; itrkid<jointytrkIDs[itrk].size(); ++itrkid){
506 if(jointytrkIDs[itrk][itrkid] == matchedtrkid){
514 for(
unsigned int ixtrk = 0; ixtrk<allXTracks.size();++ixtrk){
515 score[ixtrk][itrk] = -1.0;
523 else{ doneMatching =
true; }
526 for(
unsigned int i = 0;
i<xTracks.size();++
i){
527 if(xMatched[
i] == 0){
532 for(
unsigned int j = 1 ;
j<xTracks[
i].NTrajectoryPoints();++
j){
535 for(
unsigned int j = 0;
j<xTracks[
i].NCell();++
j){
536 track.
Add(xTracks[
i].Cell(
j));
538 track.
SetID(xTracks[
i].ID());
539 unmatchedtracks.push_back(track);
542 for(
unsigned int i = 0;
i<yTracks.size();++
i){
543 if(yMatched[
i] == 0){
548 for(
unsigned int j = 1 ;
j<yTracks[
i].NTrajectoryPoints();++
j){
551 for(
unsigned int j = 0;
j<yTracks[
i].NCell();++
j){
552 track.
Add(yTracks[
i].Cell(
j));
554 track.
SetID(yTracks[
i].ID());
555 unmatchedtracks.push_back(track);
568 TVector3 start1 = trka.
Start();
569 TVector3 stop1 = trka.
Stop();
570 TVector3 start2 = trkb.
Start();
571 TVector3 stop2 = trkb.
Stop();
576 double stst = (start1-start2).
Mag();
577 double stsp = (start1-stop2).
Mag();
578 double spst = (stop1-start2).
Mag();
579 double spsp = (stop1-stop2).
Mag();
587 if(stst < stsp && stst < spst && stst< spsp){
593 else if(stsp < stst && stsp < spst && stsp < spsp){
599 else if(spst< stst && spst< stsp && spst < spsp){
612 double jointv1 = joint1.X();
613 double jointv2 = joint2.X();
614 double endv1 = end1.X();
615 double endv2 = end2.X();
617 jointv1 = joint1.Y();
618 jointv2 = joint2.Y();
627 if(joint1.Z() == start1.Z()){plane1 = trka.
MinPlane(); }
629 if(joint2.Z() == start2.Z()){plane2 = trkb.
MinPlane(); }
636 if(plane1 == plane2){
638 if( (joint1.Z()-end1.Z())*(end2.Z()-joint1.Z()) < 0 ){
643 if( (joint1.Z()-end1.Z())*(joint2.Z()-joint1.Z()) < 0 ||
644 (joint2.Z()-joint1.Z())*(end2.Z()-joint2.Z()) < 0 ){
654 if(joint1.Z() == end1.Z() && joint1.Z() < end2.Z()){
return false; }
655 else if(joint2.Z() == end2.Z() && joint2.Z() < end1.Z()){
return false; }
661 for(
unsigned int iCell = 0; iCell<trka.
NCell(); ++iCell){
663 if(trka.
Cell(iCell)->
Cell()>cellhigh1){ cellhigh1 = trka.
Cell(iCell)->
Cell(); }
664 if(trka.
Cell(iCell)->
Cell()<celllow1){ celllow1 = trka.
Cell(iCell)->
Cell();}
670 for(
unsigned int iCell = 0; iCell<trkb.
NCell(); ++iCell){
672 if(trkb.
Cell(iCell)->
Cell()>cellhigh2){ cellhigh2 = trkb.
Cell(iCell)->
Cell(); }
673 if(trkb.
Cell(iCell)->
Cell()<celllow2){ celllow2 = trkb.
Cell(iCell)->
Cell();}
680 if(jointv1 > jointv2){
682 if(celllow1 <= cellhigh2){ cellDiff = celllow1- cellhigh2; }
683 else if(cellhigh2 >= celllow1 && cellhigh2 <= cellhigh1){ cellDiff = 0; }
684 else{ cellDiff = -(cellhigh2- celllow1); }
688 if(cellhigh1 <= celllow2){ cellDiff = celllow2-cellhigh1; }
689 else if(cellhigh1 >= celllow2 && cellhigh1 <= cellhigh2){ cellDiff = 0;}
690 else{ cellDiff = -(cellhigh1-celllow2); }
693 if(plane1 == plane2 &&
abs(cellDiff) > 1){
return false; }
696 TVector3 diff1 = joint1-end1;
697 TVector3 jointdiff = joint2-joint1;
698 TVector3 diff2 = end2-joint2;
699 double cos1 = diff1.Dot(jointdiff)/(diff1.Mag()*jointdiff.Mag());
700 double cos2 = diff2.Dot(jointdiff)/(diff2.Mag()*jointdiff.Mag());
701 double slope1 = diff1.X()/diff1.Z();
702 double slope2 = diff2.X()/diff2.Z();
703 double slopediff = jointdiff.X()/jointdiff.Z();
705 slope1 = diff1.Y()/diff1.Z();
706 slope2 = diff2.Y()/diff2.Z();
707 slopediff = jointdiff.Y()/jointdiff.Z();
710 double rotAngle1 = TMath::ATan(slope1);
711 double rotAngle2 = TMath::ATan(slopediff);
712 double joint2zrot = TMath::Cos(rotAngle1)*(joint2.Z() - joint1.Z()) - TMath::Sin(rotAngle1)*(jointv2 - jointv1);
713 if(joint2zrot < 0){ cos1 = -cos1; }
714 double end2zrot = TMath::Cos(rotAngle2)*(end2.Z() - joint2.Z()) - TMath::Sin(rotAngle2)*(endv2 - jointv2);
715 if(end2zrot < 0){ cos2 = -cos2; }
720 geo::LineIntersection(endv1,end1.Z(),jointv1,joint1.Z(),endv2,end2.Z(),jointv2,joint2.Z(),vint,zint);
722 double dx = jointv1-endv1;
723 double dy = joint1.Z()-end1.Z();
724 double dX = jointv2-endv2;
725 double dY = joint2.Z()-end2.Z();
726 double denom = dX*dy-dY*
dx;
727 if(denom == 0 && cellDiff>1){
return false; }
732 double zlow = xyzpl1[2]-dz1;
733 double vlow = xyzpl1[
view]-dv1;
738 double zhigh = xyzpl2[2]+dz2;
739 double vhigh = xyzpl2[
view]+dv2;
741 if(joint2.Z() < joint1.Z()){
742 zlow = xyzpl2[2]-dz2;
743 zhigh = xyzpl1[2]+dz1;
745 if((slope1*slope2) > 0 && plane1 != plane2 &&
abs(cellDiff)>1){
746 if(jointv2 < jointv1){
749 vlow = xyzpl2[
view]-dv2;
750 vhigh = xyzpl1[
view]+dv1;
752 if(vint < vlow || vint > vhigh || zint < zlow || zint > zhigh){
return false; }
754 else if((slope1*slope2) <= 0 && plane1 != plane2 &&
abs(cellDiff)>1){
765 if(zint < zlow || zint > zhigh || vint < vlow || vint > vhigh){
return false; }
767 else if(plane1 == plane2){
769 if(zint < zlow-dz1 || zint > zhigh+dz1){
return false; }
774 std::vector<geo::OfflineChan> hitsOnline;
776 if(plane1 != plane2){
779 std::vector<geo::OfflineChan> hitsOnline2;
780 std::vector<double> distance2;
790 for(
int imiss2 = 0; imiss2 < numMissed2; ++imiss2){
792 hitsOnline.push_back(hitsOnline2[imiss2]);
793 distance.push_back(distance2[imiss2]);
796 std::vector<int> missesToRemove;
797 for(
unsigned int imiss = 0; imiss<hitsOnline.size();++imiss){
798 if(hitsOnline[imiss].Plane() == plane1 && hitsOnline[imiss].Cell() >= celllow1 && hitsOnline[imiss].Cell() <= cellhigh1){
799 missesToRemove.push_back(imiss);
801 else if(hitsOnline[imiss].Plane() == plane2 && hitsOnline[imiss].Cell() >= celllow2 && hitsOnline[imiss].Cell() <= cellhigh2){
802 missesToRemove.push_back(imiss);
806 reverse(missesToRemove.begin(),missesToRemove.end());
807 for(
unsigned int imiss = 0; imiss < missesToRemove.size(); ++imiss){
809 hitsOnline.erase(hitsOnline.begin()+missesToRemove[imiss]);
810 distance.erase(distance.begin()+missesToRemove[imiss]);
815 unsigned short lowestCell = cellhigh1;
816 unsigned short highestCell = celllow2;
817 if(jointv1 > jointv2){
818 lowestCell = cellhigh2;
819 highestCell = celllow1;
821 for(
unsigned short imiss = lowestCell+1; imiss<highestCell;++imiss){
823 hitsOnline.push_back(miss);
830 for(
unsigned int i = 0;
i < distance.size(); ++
i){
834 double zave = 0.5*joint1.Z()+0.5*joint2.Z();
837 bool foundrng =
false;
858 if(!foundrng){ w =
avgYPos[0]; }
862 bool foundrng =
false;
883 if(!foundrng){ w =
avgXPos[0]; }
888 double minProb = 0.00001;
889 double minScatProb = 0.1;
894 if(probMiss > minProb){
896 if((cos1+cos2)<=0){
return false; }
897 if(probMiss < 2.0*minProb/(cos1+cos2)){
return false; }
898 if(probMiss < minScatProb){
return false; }
902 if(
abs(cellDiff) > 1){
return false; }
903 double trkcos1 = TMath::Cos(TMath::ATan(slope1));
904 double trkcos2 = TMath::Cos(TMath::ATan(slope2));
905 double cosmax = trkcos1;
906 double cosmin = trkcos2;
907 if(trkcos2 > trkcos1){
911 if(cosmin < 0.98 || cosmax < 0.995){
return false; }
932 TVector3 start1 = trka.
Start();
933 TVector3 stop1 = trka.
Stop();
934 TVector3 start2 = trkb.
Start();
935 TVector3 stop2 = trkb.
Stop();
936 double stst = (start1-start2).
Mag();
937 double stsp = (start1-stop2).
Mag();
938 double spst = (stop1-start2).
Mag();
939 double spsp = (stop1-stop2).
Mag();
943 std::vector<double> zends;
944 zends.push_back(start1.Z());
945 zends.push_back(stop1.Z());
946 zends.push_back(start2.Z());
947 zends.push_back(stop2.Z());
948 std::vector<TVector3> trkends;
949 trkends.push_back(start1);
950 trkends.push_back(stop1);
951 trkends.push_back(start2);
952 trkends.push_back(stop2);
958 int zplanes[4] = {trkaminpl, trkamaxpl, trkbminpl, trkbmaxpl};
960 for(
int iz = 0; iz < 4; ++iz){
972 std::vector<TVector3> traja(ntraja);
973 for(
int i = 0;
i<ntraja; ++
i){
976 std::vector<TVector3> trajb(ntrajb);
977 for(
int i = 0;
i<ntrajb; ++
i){
980 std::vector<TVector3> trajpts;
981 bool sameplane =
false;
983 double plhitspetot[4] = {0, 0, 0, 0};
984 for(
unsigned int i = 0;
i<trka.
NCell(); ++
i){
988 for(
unsigned int i = 0;
i<trkb.
NCell(); ++
i){
993 if(stst < stsp && stst < spst && stst< spsp){
995 if(zplanes[0] == zplanes[2]){ sameplane =
true; }
997 if(trajb[ntrajb-1].
Z() < traja[ntraja-1].
Z()){
1000 for(
int i = 0;
i<ntrajb-1;++
i){
1001 trajpts.push_back(trajb[
i]);
1004 trajpts.push_back(trajb[ntrajb-1]);
1005 trajpts.push_back(traja[0]);
1008 TVector3 mergetrj = (plhitspetot[0]*traja[0]+plhitspetot[2]*trajb[ntrajb-1])*
1009 (1.0/(plhitspetot[0]+plhitspetot[2]));
1011 else{ trajpts.push_back(traja[0]); }
1013 for(
int i = 1;
i<ntraja;++
i){
1014 trajpts.push_back(traja[
i]);
1020 for(
int i = 0;
i<ntraja-1;++
i){
1021 trajpts.push_back(traja[
i]);
1024 trajpts.push_back(traja[ntraja-1]);
1025 trajpts.push_back(trajb[0]);
1028 TVector3 mergetrj = (plhitspetot[0]*traja[ntraja-1]+plhitspetot[2]*trajb[0])*
1029 (1.0/(plhitspetot[0]+plhitspetot[2]));
1031 else{ trajpts.push_back(traja[ntraja-1]); }
1033 for(
int i = 1;
i<ntrajb;++
i){
1034 trajpts.push_back(trajb[
i]);
1038 else if(stsp < stst && stsp < spst && stsp < spsp){
1040 if(zplanes[0] == zplanes[3]){ sameplane =
true; }
1042 if(trajb[0].
Z() < traja[ntraja-1].
Z()){
1044 for(
int i = 0;
i<ntrajb-1;++
i){
1045 trajpts.push_back(trajb[
i]);
1048 trajpts.push_back(trajb[ntrajb-1]);
1049 trajpts.push_back(traja[0]);
1052 TVector3 mergetrj = (plhitspetot[0]*traja[0]+plhitspetot[3]*trajb[ntrajb-1])*
1053 (1.0/(plhitspetot[0]+plhitspetot[3]));
1055 else{ trajpts.push_back(traja[0]); }
1057 for(
int i = 1;
i<ntraja;++
i){
1058 trajpts.push_back(traja[
i]);
1064 for(
int i = 0;
i<ntraja-1;++
i){
1065 trajpts.push_back(traja[
i]);
1070 trajpts.push_back(traja[ntraja-1]);
1071 trajpts.push_back(trajb[0]);
1074 TVector3 mergetrj = (plhitspetot[0]*traja[ntraja-1]+plhitspetot[3]*trajb[0])*
1075 (1.0/(plhitspetot[0]+plhitspetot[3]));
1077 else{ trajpts.push_back(traja[ntraja-1]); }
1079 for(
int i = 1;
i<ntrajb;++
i){
1080 trajpts.push_back(trajb[
i]);
1084 else if(spst< stst && spst< stsp && spst < spsp){
1086 if(zplanes[1] == zplanes[2]){ sameplane =
true; }
1088 if(trajb[ntrajb-1].
Z() < traja[0].
Z()){
1091 for(
int i = 0;
i<ntrajb-1;++
i){
1092 trajpts.push_back(trajb[
i]);
1097 trajpts.push_back(trajb[ntrajb-1]);
1098 trajpts.push_back(traja[0]);
1101 TVector3 mergetrj = (plhitspetot[1]*traja[0]+plhitspetot[2]*trajb[ntrajb-1])*
1102 (1.0/(plhitspetot[1]+plhitspetot[2]));
1104 else{ trajpts.push_back(traja[0]); }
1106 for(
int i = 1;
i<ntraja;++
i){
1107 trajpts.push_back(traja[
i]);
1112 for(
int i = 0;
i<ntraja-1;++
i){
1113 trajpts.push_back(traja[
i]);
1116 trajpts.push_back(traja[ntraja-1]);
1117 trajpts.push_back(trajb[0]);
1120 TVector3 mergetrj = (plhitspetot[1]*traja[ntraja-1]+plhitspetot[2]*trajb[0])*
1121 (1.0/(plhitspetot[1]+plhitspetot[2]));
1123 else{ trajpts.push_back(traja[ntraja-1]); }
1126 for(
int i = 1;
i<ntrajb;++
i){
1127 trajpts.push_back(trajb[
i]);
1133 if(zplanes[1] == zplanes[3]){ sameplane =
true; }
1135 if(trajb[0].
Z() < traja[0].
Z()){
1137 for(
int i = 0;
i<ntrajb-1;++
i){
1138 trajpts.push_back(trajb[
i]);
1143 trajpts.push_back(trajb[ntrajb-1]);
1144 trajpts.push_back(traja[0]);
1147 TVector3 mergetrj = (plhitspetot[1]*traja[0]+plhitspetot[3]*trajb[ntrajb-1])*
1148 (1.0/(plhitspetot[1]+plhitspetot[3]));
1150 else{ trajpts.push_back(traja[0]); }
1152 for(
int i = 1;
i<ntraja;++
i){
1153 trajpts.push_back(traja[
i]);
1158 for(
int i = 0;
i<ntraja-1;++
i){
1159 trajpts.push_back(traja[
i]);
1164 trajpts.push_back(traja[ntraja-1]);
1165 trajpts.push_back(trajb[0]);
1168 TVector3 mergetrj = (plhitspetot[1]*traja[ntraja-1]+plhitspetot[3]*trajb[0])*
1169 (1.0/(plhitspetot[1]+plhitspetot[3]));
1171 else{ trajpts.push_back(traja[ntraja-1]); }
1173 for(
int i = 1;
i<ntrajb;++
i){
1174 trajpts.push_back(trajb[
i]);
1178 TVector3
dir = trajpts[1]-trajpts[0];
1179 rb::Track JointTrack(JointCluster,trajpts[0].
X(),trajpts[0].
Z(),dir.X(),dir.Z(),trka.
ID());
1181 JointTrack.
SetStart(trajpts[0].
Y(),trajpts[0].
Z());
1182 JointTrack.SetDir(dir.Y(),dir.Z());
1184 for(
unsigned int i = 1;
i<trajpts.size();++
i){
1185 if(trka.
View() ==
geo::kX){ JointTrack.AppendTrajectoryPoint(trajpts[
i].
X(),trajpts[
i].
Z()); }
1186 else if(trka.
View() ==
geo::kY){ JointTrack.AppendTrajectoryPoint(trajpts[
i].
Y(),trajpts[
i].
Z()); }
1194 double score = -1.0;
1224 if(!xgood || !ygood){
return score;}
1226 TVector3 ixtraj = xTrack.
Start();
1227 TVector3 fxtraj = xTrack.
Stop();
1228 TVector3 iytraj = yTrack.
Start();
1229 TVector3 fytraj = yTrack.
Stop();
1232 if(fxtraj.Z() < ixtraj.Z()){
std::swap(ixtraj,fxtraj); }
1233 if(fytraj.Z() < iytraj.Z()){
std::swap(iytraj,fytraj); }
1242 if(iyPlane <= fxPlane && fyPlane >= ixPlane){
1245 int minoverlapplane;
1246 int maxoverlapplane;
1247 if(ixPlane<iyPlane){
1248 minoverlapplane = iyPlane;
1249 if(fxPlane<fyPlane){ maxoverlapplane = fxPlane; }
1250 else{ maxoverlapplane = fyPlane; }
1253 minoverlapplane = ixPlane;
1254 if(fxPlane<fyPlane){ maxoverlapplane = fxPlane; }
1255 else{ maxoverlapplane = fyPlane; }
1258 std::vector<art::Ptr<rb::CellHit> > allhits;
1259 for(
unsigned int ihit = 0; ihit<trka.
NCell(); ++ihit){
1260 allhits.push_back(trka.
Cell(ihit));
1262 for(
unsigned int ihit = 0; ihit<trkb.
NCell(); ++ihit){
1263 allhits.push_back(trkb.
Cell(ihit));
1269 for(
int iPlane = minoverlapplane; iPlane<=maxoverlapplane; ++iPlane){
1271 bool hitcurr =
false;
1272 bool hitprev =
false;
1273 bool hitnext =
false;
1274 for(
unsigned int ihit = 0; ihit < allhits.size();++ihit){
1275 if(!hitcurr && allhits[ihit]->Plane() == iPlane){ hitcurr =
true; }
1276 if(!hitprev && allhits[ihit]->Plane() == iPlane-1){ hitprev =
true; }
1277 if(!hitnext && allhits[ihit]->Plane() == iPlane+1){hitnext =
true; }
1279 if(hitcurr && hitprev && hitnext){ overlap+=1.0; }
1282 int idiff = TMath::Abs(iyPlane-ixPlane);
1283 int fdiff = TMath::Abs(fyPlane-fxPlane);
1285 score = ((double)idiff+fdiff)/overlap;
1290 if(ixPlane == (iyPlane-1) && fxPlane == (iyPlane+1)){score = 1.0;}
1295 if(iyPlane == (ixPlane-1) && fyPlane == (ixPlane+1)){score = 1.0;}
1298 else if(ixPlane == fxPlane && iyPlane == fyPlane){
1300 if(ixPlane == (iyPlane-1) || ixPlane == (iyPlane+1)){score = 1.0;}
1311 for(
unsigned int i = 0;
i<xTrack.
NCell();++
i){
1314 for(
unsigned int i = 0;
i<yTrack.
NCell();++
i){
1319 std::vector<int> xPlanes;
1323 std::vector<int> yPlanes;
1329 std::vector<TVector3> xtrjpoints;
1330 double minxz = 99999;
1334 xPlanes.push_back(xplane);
1337 for(
unsigned int i =0;
i<xPlanes.size(); ++
i){
1338 if(
i == 0 ||
i == xPlanes.size()-1){
1340 if(xtrjpoints.back().Z() > maxxz){ maxxz = xtrjpoints.back().Z(); }
1341 if(xtrjpoints.back().Z() < minxz){ minxz = xtrjpoints.back().Z(); }
1343 else if(!(xPlanes[
i] == xPlanes[
i-1] && xPlanes[
i] == xPlanes[
i+1])){
1345 if(xtrjpoints.back().Z() > maxxz){ maxxz = xtrjpoints.back().Z(); }
1346 if(xtrjpoints.back().Z() < minxz){ minxz = xtrjpoints.back().Z(); }
1350 std::vector<TVector3> ytrjpoints;
1351 double minyz = 99999;
1355 yPlanes.push_back(yplane);
1358 for(
unsigned int i = 0;
i<yPlanes.size(); ++
i){
1359 if(
i == 0 ||
i == yPlanes.size()-1){
1361 if(ytrjpoints.back().Z() > maxyz){ maxyz = ytrjpoints.back().Z(); }
1362 if(ytrjpoints.back().Z() < minyz){ minyz = ytrjpoints.back().Z(); }
1364 else if(!(yPlanes[
i] == yPlanes[
i-1] && yPlanes[
i] == yPlanes[
i+1])){
1366 if(ytrjpoints.back().Z() > maxyz){ maxyz = ytrjpoints.back().Z(); }
1367 if(ytrjpoints.back().Z() < minyz){ minyz = ytrjpoints.back().Z(); }
1374 if(xtrjpoints.size() >= 2){
1375 unsigned int xplane = xTrack.
MinPlane();
1377 if(xtrjpoints[0].
Z() <= xtrjpoints[xtrjpoints.size()-1].Z()){
1378 TVector3 endtraj = xtrjpoints[0];
1379 TVector3 nexttraj = xtrjpoints[1];
1380 if(
fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1381 TVector3 lowztraj, highztraj;
1383 xtrjpoints[0] = lowztraj;
1384 xtrjpoints[1] = highztraj;
1388 TVector3 endtraj = xtrjpoints[xtrjpoints.size()-1];
1389 TVector3 nexttraj = xtrjpoints[xtrjpoints.size()-2];
1390 if(
fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1391 TVector3 lowztraj, highztraj;
1393 xtrjpoints[xtrjpoints.size()-1] = lowztraj;
1394 xtrjpoints[xtrjpoints.size()-2] = highztraj;
1400 if(ytrjpoints.size() >= 2){
1401 unsigned int yplane = yTrack.
MinPlane();
1403 if(ytrjpoints[0].
Z() <= ytrjpoints[ytrjpoints.size()-1].Z()){
1404 TVector3 endtraj = ytrjpoints[0];
1405 TVector3 nexttraj = ytrjpoints[1];
1406 if(
fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1407 TVector3 lowztraj, highztraj;
1409 ytrjpoints[0] = lowztraj;
1410 ytrjpoints[1] = highztraj;
1414 TVector3 endtraj = ytrjpoints[ytrjpoints.size()-1];
1415 TVector3 nexttraj = ytrjpoints[ytrjpoints.size()-2];
1416 if(
fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1417 TVector3 lowztraj, highztraj;
1419 ytrjpoints[ytrjpoints.size()-1] = lowztraj;
1420 ytrjpoints[ytrjpoints.size()-2] = highztraj;
1428 if(xtrjpoints.size() >= 2){
1429 unsigned int xplane = xTrack.
MaxPlane();
1431 if(xtrjpoints[0].
Z() >= xtrjpoints[xtrjpoints.size()-1].Z()){
1432 TVector3 endtraj = xtrjpoints[0];
1433 TVector3 nexttraj = xtrjpoints[1];
1434 if(
fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1435 TVector3 lowztraj, highztraj;
1437 xtrjpoints[0] = highztraj;
1438 xtrjpoints[1] = lowztraj;
1442 TVector3 endtraj = xtrjpoints[xtrjpoints.size()-1];
1443 TVector3 nexttraj = xtrjpoints[xtrjpoints.size()-2];
1444 if(
fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1445 TVector3 lowztraj, highztraj;
1447 xtrjpoints[xtrjpoints.size()-1] = highztraj;
1448 xtrjpoints[xtrjpoints.size()-2] = lowztraj;
1454 if(ytrjpoints.size() >= 2){
1455 unsigned int yplane = yTrack.
MaxPlane();
1457 if(ytrjpoints[0].
Z() >= ytrjpoints[ytrjpoints.size()-1].Z()){
1458 TVector3 endtraj = ytrjpoints[0];
1459 TVector3 nexttraj = ytrjpoints[1];
1460 if(
fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1461 TVector3 lowztraj, highztraj;
1463 ytrjpoints[0] = highztraj;
1464 ytrjpoints[1] = lowztraj;
1468 TVector3 endtraj = ytrjpoints[ytrjpoints.size()-1];
1469 TVector3 nexttraj = ytrjpoints[ytrjpoints.size()-2];
1470 if(
fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1471 TVector3 lowztraj, highztraj;
1473 ytrjpoints[ytrjpoints.size()-1] = highztraj;
1474 ytrjpoints[ytrjpoints.size()-2] = lowztraj;
1482 for(
unsigned int i = 0;
i<xtrjpoints.size(); ++
i){
1483 if(
i == 0){ xTrack.
SetStart(xtrjpoints[0].
X(),xtrjpoints[0].
Z()); }
1488 for(
unsigned int i = 0;
i<ytrjpoints.size(); ++
i){
1489 if(
i == 0){ yTrack.
SetStart(ytrjpoints[0].
Y(),ytrjpoints[0].
Z()); }
1494 bool backscatterx =
false;
1495 std::vector<int> xbackscatind;
1497 for(
unsigned int i = 0;
i<xtrjpoints.size();++
i){
1498 if(
i > 0 &&
i < xtrjpoints.size()-1){
1499 double deltaz1 = xtrjpoints[
i].Z()-xtrjpoints[
i-1].Z();
1500 double deltaz2 = xtrjpoints[
i+1].Z()-xtrjpoints[
i].Z();
1501 if((deltaz1*deltaz2)<0){ xbackscatind.push_back(
i); }
1505 if(xbackscatind.size()>0){ backscatterx =
true; }
1507 bool backscattery =
false;
1508 std::vector<int> ybackscatind;
1509 for(
unsigned int i = 0;
i<ytrjpoints.size();++
i){
1510 if(
i > 0 &&
i < ytrjpoints.size()-1){
1511 double deltaz1 = ytrjpoints[
i].Z()-ytrjpoints[
i-1].Z();
1512 double deltaz2 = ytrjpoints[
i+1].Z()-ytrjpoints[
i].Z();
1513 if((deltaz1*deltaz2)<0){ ybackscatind.push_back(
i); }
1517 if(ybackscatind.size()>0){ backscattery =
true; }
1519 std::vector<TVector3> trjpoints;
1520 if((!backscatterx && !backscattery) || (backscatterx && backscattery)){
1521 std::vector<TVector3> newtrjpoints;
1522 for(
unsigned int i = 0;
i<xtrjpoints.size();++
i){
1525 xtrjpoints[
i].SetY(y);
1526 newtrjpoints.push_back(xtrjpoints[
i]);
1528 for(
unsigned int i = 0;
i<ytrjpoints.size();++
i){
1531 ytrjpoints[
i].SetX(x);
1532 newtrjpoints.push_back(ytrjpoints[
i]);
1534 std::sort(newtrjpoints.begin(),newtrjpoints.end(),
sort_traj);
1535 trjpoints = newtrjpoints;
1537 if((backscatterx && !backscattery)){
1539 std::vector<TVector3> newtrjpoints;
1540 for(
unsigned int i = 0;
i<xtrjpoints.size();++
i){
1543 xtrjpoints[
i].SetY(y);
1544 newtrjpoints.push_back(xtrjpoints[
i]);
1548 xbackscatind.insert(xbackscatind.begin(),0);
1549 xbackscatind.push_back(xtrjpoints.size()-1);
1551 int nback = xbackscatind.size();
1555 std::vector<std::vector<int> > xsectytrj(ytrjpoints.size());
1556 for(
unsigned int iytraj = 0; iytraj<ytrjpoints.size(); ++iytraj){
1557 double zy = ytrjpoints[iytraj].Z();
1558 bool foundsect =
false;
1560 for(
unsigned int ixbrk = 0; ixbrk<xbackscatind.size()-1; ++ixbrk){
1561 double z1 = xtrjpoints[xbackscatind[ixbrk]].Z();
1562 double z2 = xtrjpoints[xbackscatind[ixbrk+1]].Z();
1563 if( (z1 < zy && zy < z2) || (z1 > zy && zy > z2) ){
1564 xsectytrj[iytraj].push_back(ixbrk);
1570 double closestdist = (xtrjpoints[xbackscatind[0]].Z() - zy);
1571 int closestxsect = 0;
1572 for(
unsigned int ixbrk = 0; ixbrk<xbackscatind.size(); ++ixbrk){
1573 double dist =
fabs(xtrjpoints[xbackscatind[ixbrk]].
Z() - zy);
1574 if(dist<closestdist){
1575 closestxsect = ixbrk;
1579 xsectytrj[iytraj].push_back(closestxsect);
1583 std::vector<TVector3> orderedtrjpoints;
1585 for(
int isect = 0; isect<nback; ++isect){
1586 std::vector<TVector3> xsecttrjs;
1587 std::vector<TVector3> ysecttrjs;
1590 if(isect != nback-1){
1591 for(
int ix = xbackscatind[isect]; ix<=xbackscatind[isect+1]; ++ix){
1592 xsecttrjs.push_back(xtrjpoints[ix]);
1596 xsecttrjs.push_back(xtrjpoints[xbackscatind[isect]-1]);
1597 xsecttrjs.push_back(xtrjpoints[xbackscatind[isect]]);
1600 for(
unsigned int iy = 0; iy<xsectytrj.size(); ++iy){
1601 for(
unsigned int iyx = 0; iyx<xsectytrj[iy].size(); ++iyx){
1602 if(xsectytrj[iy][iyx] == isect){ ysecttrjs.push_back(ytrjpoints[iy]); }
1606 for(
unsigned int iy = 0; iy<ysecttrjs.size(); ++iy){
1608 double zy = ysecttrjs[iy].Z();
1609 bool foundngbr =
false;
1611 for(
unsigned int ix = 0; ix<xsecttrjs.size()-1; ++ix){
1612 double z1 = xsecttrjs[ix].Z();
1613 double z2 = xsecttrjs[ix+1].Z();
1614 if( (z1 < zy && zy < z2) || (z1 > zy && zy > z2) ){
1616 double xint = xsecttrjs[ix].X() + (zy- z1)*(xsecttrjs[ix+1].
X()-xsecttrjs[ix].X())/(z2-z1);
1617 ysecttrjs[iy].SetX(xint);
1623 if(
fabs(xsecttrjs[0].
Z() - zy) <
fabs(xsecttrjs[xsecttrjs.size()-1].Z() - zy)){
1624 double x1 = xsecttrjs[0].X();
1625 double x2 = xsecttrjs[1].X();
1626 double z1 = xsecttrjs[0].Z();
1627 double z2 = xsecttrjs[1].Z();
1628 double xint = x1 + (zy-z1)*(x2-x1)/(z2-z1);
1629 ysecttrjs[iy].SetX(xint);
1632 double x1 = xsecttrjs[xsecttrjs.size()-1].X();
1633 double x2 = xsecttrjs[xsecttrjs.size()-2].X();
1634 double z1 = xsecttrjs[xsecttrjs.size()-1].Z();
1635 double z2 = xsecttrjs[xsecttrjs.size()-2].Z();
1636 double xint = x1 + (zy-z1)*(x2-x1)/(z2-z1);
1637 ysecttrjs[iy].SetX(xint);
1642 std::vector<TVector3> alltrjs;
1643 for(
unsigned int i = 0;
i<xsecttrjs.size();++
i){
1644 if(isect < nback-1){
1645 if(
i != xsecttrjs.size()-1){ alltrjs.push_back(xsecttrjs[
i]); }
1648 if(
i != 0){ alltrjs.push_back(xsecttrjs[
i]); }
1652 for(
unsigned int i = 0;
i<ysecttrjs.size();++
i){
1653 alltrjs.push_back(ysecttrjs[
i]);
1656 double zDir = xsecttrjs[xsecttrjs.size()-1].Z() - xsecttrjs[0].Z();
1657 std::sort(alltrjs.begin(),alltrjs.end(),
sort_traj);
1658 if(zDir < 0){
reverse(alltrjs.begin(),alltrjs.end()); }
1659 for(
unsigned int i = 0;
i<alltrjs.size(); ++
i){
1660 orderedtrjpoints.push_back(alltrjs[
i]);
1663 trjpoints = orderedtrjpoints;
1665 if((!backscatterx && backscattery)){
1667 std::vector<TVector3> newtrjpoints;
1668 for(
unsigned int i = 0;
i<ytrjpoints.size();++
i){
1671 ytrjpoints[
i].SetX(x);
1672 newtrjpoints.push_back(ytrjpoints[
i]);
1676 ybackscatind.insert(ybackscatind.begin(),0);
1677 ybackscatind.push_back(ytrjpoints.size()-1);
1679 int nback = ybackscatind.size();
1683 std::vector<std::vector<int> > ysectxtrj(xtrjpoints.size());
1684 for(
unsigned int ixtraj = 0; ixtraj<xtrjpoints.size(); ++ixtraj){
1685 double zx = xtrjpoints[ixtraj].Z();
1686 bool foundsect =
false;
1688 for(
unsigned int iybrk = 0; iybrk<ybackscatind.size()-1; ++iybrk){
1689 double z1 = ytrjpoints[ybackscatind[iybrk]].Z();
1690 double z2 = ytrjpoints[ybackscatind[iybrk+1]].Z();
1691 if( (z1 < zx && zx < z2) || (z1 > zx && zx > z2) ){
1692 ysectxtrj[ixtraj].push_back(iybrk);
1698 double closestdist = (ytrjpoints[ybackscatind[0]].Z() - zx);
1699 int closestysect = 0;
1700 for(
unsigned int iybrk = 0; iybrk<ybackscatind.size(); ++iybrk){
1701 double dist =
fabs(ytrjpoints[ybackscatind[iybrk]].
Z() - zx);
1702 if(dist<closestdist){
1703 closestysect = iybrk;
1707 ysectxtrj[ixtraj].push_back(closestysect);
1711 std::vector<TVector3> orderedtrjpoints;
1713 for(
int isect = 0; isect<nback; ++isect){
1714 std::vector<TVector3> xsecttrjs;
1715 std::vector<TVector3> ysecttrjs;
1718 if(isect != nback-1){
1719 for(
int iy = ybackscatind[isect]; iy<=ybackscatind[isect+1]; ++iy){
1720 ysecttrjs.push_back(ytrjpoints[iy]);
1724 ysecttrjs.push_back(ytrjpoints[ybackscatind[isect]-1]);
1725 ysecttrjs.push_back(ytrjpoints[ybackscatind[isect]]);
1728 for(
unsigned int ix = 0; ix<ysectxtrj.size(); ++ix){
1729 for(
unsigned int ixy = 0; ixy<ysectxtrj[ix].size(); ++ixy){
1730 if(ysectxtrj[ix][ixy] == isect){ xsecttrjs.push_back(xtrjpoints[ix]); }
1734 for(
unsigned int ix = 0; ix<xsecttrjs.size(); ++ix){
1736 double zx = xsecttrjs[ix].Z();
1737 bool foundngbr =
false;
1739 for(
unsigned int iy = 0; iy<ysecttrjs.size()-1; ++iy){
1740 double z1 = ysecttrjs[iy].Z();
1741 double z2 = ysecttrjs[iy+1].Z();
1742 if( (z1 < zx && zx < z2) || (z1 > zx && zx > z2) ){
1744 double yint = ysecttrjs[iy].Y() + (zx- z1)*(ysecttrjs[iy+1].
Y()-ysecttrjs[iy].Y())/(z2-z1);
1745 xsecttrjs[ix].SetY(yint);
1751 if(
fabs(ysecttrjs[0].
Z() - zx) <
fabs(ysecttrjs[ysecttrjs.size()-1].Z() - zx)){
1752 double y1 = ysecttrjs[0].Y();
1753 double y2 = ysecttrjs[1].Y();
1754 double z1 = ysecttrjs[0].Z();
1755 double z2 = ysecttrjs[1].Z();
1756 double yint = y1 + (zx-z1)*(y2-y1)/(z2-z1);
1757 xsecttrjs[ix].SetY(yint);
1760 double y1 = ysecttrjs[ysecttrjs.size()-1].Y();
1761 double y2 = ysecttrjs[ysecttrjs.size()-2].Y();
1762 double z1 = ysecttrjs[ysecttrjs.size()-1].Z();
1763 double z2 = ysecttrjs[ysecttrjs.size()-2].Z();
1764 double yint = y1 + (zx-z1)*(y2-y1)/(z2-z1);
1765 xsecttrjs[ix].SetY(yint);
1770 std::vector<TVector3> alltrjs;
1771 for(
unsigned int i = 0;
i<ysecttrjs.size();++
i){
1772 if(isect < nback-1){
1773 if(
i != ysecttrjs.size()-1){
1774 alltrjs.push_back(ysecttrjs[
i]);
1779 alltrjs.push_back(ysecttrjs[
i]);
1783 for(
unsigned int i = 0;
i<xsecttrjs.size();++
i){
1784 alltrjs.push_back(xsecttrjs[
i]);
1787 double zDir = ysecttrjs[ysecttrjs.size()-1].Z() - ysecttrjs[0].Z();
1788 std::sort(alltrjs.begin(),alltrjs.end(),
sort_traj);
1789 if(zDir < 0){
reverse(alltrjs.begin(),alltrjs.end()); }
1790 for(
unsigned int i = 0;
i<alltrjs.size(); ++
i){
1791 orderedtrjpoints.push_back(alltrjs[
i]);
1794 trjpoints = orderedtrjpoints;
1798 for(
unsigned int i = 1;
i<trjpoints.size();++
i){
1802 double tanThetax = TMath::Tan(TMath::ACos(xTrack.
Dir().X()));
1803 double tanThetay = TMath::Tan(TMath::ACos(yTrack.
Dir().Y()));
1804 double cosThetax = 1/(tanThetax*
sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1805 double cosThetay = 1/(tanThetay*
sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1806 TVector3 newDir(cosThetax,cosThetay,
sqrt(1-cosThetax*cosThetax-cosThetay*cosThetay));
1808 if(trjpoints.size() >= 2){
1809 newDir.SetXYZ(trjpoints[1].
X()-trjpoints[0].
X(),trjpoints[1].
Y()-trjpoints[0].
Y(),trjpoints[1].
Z()-trjpoints[0].
Z());
1824 double yzlow = yLow[2]-ydeltaz;
1825 double ylow = yLow[1];
1826 double yzhigh = yHigh[2]+ydeltaz;
1827 double yhigh = yHigh[1];
1833 double xzlow = xLow[2]-xdeltaz;
1834 double xlow = xLow[0];
1835 double xzhigh = xHigh[2]+xdeltaz;
1836 double xhigh = xHigh[0];
1838 TVector3
start(xlow+(xhigh-xlow)*(yzlow-xzlow)/(xzhigh-xzlow),ylow,yzlow);
1839 TVector3 stop(xhigh,yhigh+(yhigh-ylow)*(xzhigh-yzhigh)/(yzhigh-yzlow),xzhigh);
1842 TVector3 diffs(stop.X()-start.X(),stop.Y()-start.Y(),stop.Z()-start.Z());
1846 TVector3
start(xlow,(yhigh-ylow)*(xzlow-yzlow)/(yzhigh-yzlow)+ylow,xzlow);
1847 TVector3 stop(xhigh+(xhigh-xlow)*(yzhigh-xzhigh)/(xzhigh-xzlow),yhigh,yzhigh);
1850 TVector3 diffs(stop.X()-start.X(),stop.Y()-start.Y(),stop.Z()-start.Z());
1862 double xzlow = xLow[2]-xdeltaz;
1863 double xlow = xLow[0];
1864 double xzhigh = xHigh[2]+xdeltaz;
1865 double xhigh = xHigh[0];
1866 double xslope = (xhigh-xlow)/(xzhigh-xzlow);
1870 TVector3
start(xlow+xslope*(z-xzlow),y,z);
1871 if(ytraj == 0){newtrack.
SetStart(start);}
1874 double tanThetax = TMath::Tan(TMath::ACos(xslope/
sqrt(1+xslope*xslope)));
1875 double tanThetay = TMath::Tan(TMath::ACos(yTrack.
Dir().Y()));
1876 double cosThetax = 1/(tanThetax*
sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1877 double cosThetay = 1/(tanThetay*
sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1878 TVector3 newDir(cosThetax,cosThetay,
sqrt(1-cosThetax*cosThetax-cosThetay*cosThetay));
1890 double yzlow = yLow[2]-deltaz;
1891 double ylow = yLow[1];
1892 double yzhigh = yHigh[2]+deltaz;
1893 double yhigh = yHigh[1];
1894 double yslope = (yhigh-ylow)/(yzhigh-yzlow);
1898 TVector3
start(x,ylow+yslope*(z-yzlow),z);
1899 if(xtraj == 0){newtrack.
SetStart(start);}
1902 double tanThetay = TMath::Tan(TMath::ACos(yslope/
sqrt(1+yslope*yslope)));
1903 double tanThetax = TMath::Tan(TMath::ACos(xTrack.
Dir().X()));
1904 double cosThetax = 1/(tanThetax*
sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1905 double cosThetay = 1/(tanThetay*
sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1906 TVector3 newDir(cosThetax,cosThetay,
sqrt(1-cosThetax*cosThetax-cosThetay*cosThetay));
1921 TVector3 &lowztraj, TVector3 &highztraj)
1926 unsigned int mincell = 10000;
1927 unsigned int maxcell = 0;
1928 for(
unsigned int icell = 0; icell < track.
NCell(); ++icell){
1930 if(track.
Cell(icell)->
Cell() > maxcell){ maxcell = track.
Cell(icell)->
Cell(); }
1931 if(track.
Cell(icell)->
Cell() < mincell){ mincell = track.
Cell(icell)->
Cell(); }
1939 double zmid = 0.5*xyzLow[2]+0.5*xyzHigh[2];
1940 lowztraj.SetZ(zmid-deltaz);
1941 highztraj.SetZ(zmid+deltaz);
1944 if(endtraj.X() <= nexttraj.X()){
1945 lowztraj.SetX(xyzLow[0]);
1946 highztraj.SetX(xyzHigh[0]);
1949 lowztraj.SetX(xyzHigh[0]);
1950 highztraj.SetX(xyzLow[0]);
1954 if(endtraj.X() <= nexttraj.X()){
1955 highztraj.SetX(xyzLow[0]);
1956 lowztraj.SetX(xyzHigh[0]);
1959 highztraj.SetX(xyzHigh[0]);
1960 lowztraj.SetX(xyzLow[0]);
1966 if(endtraj.Y() <= nexttraj.Y()){
1967 lowztraj.SetY(xyzLow[1]);
1968 highztraj.SetY(xyzHigh[1]);
1971 lowztraj.SetY(xyzHigh[1]);
1972 highztraj.SetY(xyzLow[1]);
1976 if(endtraj.Y() <= nexttraj.Y()){
1977 highztraj.SetY(xyzLow[1]);
1978 lowztraj.SetY(xyzHigh[1]);
1981 highztraj.SetY(xyzHigh[1]);
1982 lowztraj.SetY(xyzLow[1]);
virtual geo::View_t View() const
kXorY for 3D clusters.
size_t NTrajectoryPoints() const
int CountMissedCellsOnLine(geo::View_t view, double vi, double zi, double vf, double zf, int maxCells)
Variation of geo::CountCellsOnlineFast function specific to KalmanTrack algorithm.
rb::Track ViewMergeTracks(KalmanGeoHelper &kgeo, rb::Track xTrack, rb::Track yTrack)
art::ServiceHandle< geo::Geometry > fgeom
Detector geometry.
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
virtual void SetStart(TVector3 start)
bool LineIntersection(double x0, double y0, double x1, double y1, double X0, double Y0, double X1, double Y1, double &x, double &y)
Find the intersection between two line segments.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
void GetCenter(double *xyz, double localz=0.0) const
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
fvar< T > fabs(const fvar< T > &x)
Float_t y1[n_points_granero]
A collection of geometry functions used by the KalmanTrack tracking algorithm.
unsigned short Plane() const
Float_t x1[n_points_granero]
void SetHit(art::Ptr< rb::CellHit > chit)
const CellGeo * Cell(int icell) const
std::string fSliceInput
Input folder of slices.
Vertical planes which measure X.
std::vector< std::vector< double > > avgYPosZRange
A collection of associated CellHits.
double MeanZ(rb::AveragingScheme scheme=kDefaultScheme) const
bool CanJoinTracks(KalmanGeoHelper &kgeo, rb::Track tracka, rb::Track trackb)
unsigned int MaxCell(geo::View_t view) const
A rb::Prong with full reconstructed trajectory.
const PlaneGeo * Plane(unsigned int i) const
Representation of a detector channel with positional information.
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
unsigned distance(const T &t1, const T &t2)
static bool sort_traj(const TVector3 &a, const TVector3 &b)
Horizontal planes which measure Y.
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
int isnan(const stan::math::var &a)
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
double CalcMatchScore(KalmanGeoHelper &kgeo, rb::Track tracka, rb::Track trackb)
ProductID put(std::unique_ptr< PROD > &&product)
unsigned short Cell() const
Track finder for cosmic rays.
rb::Track JoinTracks(KalmanGeoHelper &kgeo, rb::Track tracka, rb::Track trackb)
KalmanTrackMerge(fhicl::ParameterSet const &pset)
double CalcMissFrac(double dist, double w, geo::View_t view)
Estimation of the mip missed hit fraction as a function of the dist through a cell and cell depth w...
void push_back(Ptr< U > const &p)
std::vector< std::vector< double > > avgXPosZRange
int MatchTrajectoryToPlaneInView(TVector3 traj, geo::View_t view, int minplane, int maxplane)
Variation of MatchTrajectoryToPlane that limits the matching plane to be in the input view...
virtual void SetDir(TVector3 dir)
void AppendTrajectoryPoint(TVector3 pt)
T get(std::string const &key) const
virtual TVector3 Dir() const
Unit vector describing prong direction.
void ClearTrajectoryPoints()
Forget about all trajectory points.
Collect Geo headers and supply basic geometry functions.
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
double DetHalfHeight() const
virtual ~KalmanTrackMerge()
bool fBadChannels
Account for bad channels in gaps ?
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
std::vector< double > avgYPos
double DetHalfWidth() const
std::string fTrackInput
Input folder from track reco.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
void SortByPlane(std::vector< trk::LocatedChan > &c)
void reconfigure(const fhicl::ParameterSet &p)
unsigned int MinCell(geo::View_t view) const
fvar< T > floor(const fvar< T > &x)
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
void ShiftInterpolationPoints(TVector3 endtraj, TVector3 nextraj, unsigned int plane, rb::Track track, TVector3 &lowztraj, TVector3 &highztraj)
assert(nhit_max >=nhit_nbins)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
virtual void InterpolateXY(double z, double &x, double &y) const
void MatchTracks(KalmanGeoHelper &kgeo, std::vector< rb::Track > xtracks, std::vector< rb::Track > ytracks, std::vector< rb::Track > &matchedtracls, std::vector< rb::Track > &unmatchedtracks)
TVector3 Stop() const
Position of the final trajectory point.
Encapsulate the cell geometry.
bool fObeyPreselection
Check rb::IsFiltered?
Encapsulate the geometry of one entire detector (near, far, ndos)
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
std::vector< double > avgXPos
void produce(art::Event &evt)
int MatchTrajectoryToPlane(TVector3 traj, int minplane, int maxplane)
Finds the plane best matched to a trajectory point This method takes a trajectory point and finds the...