110 std::unique_ptr< std::vector<rec::Vee> >
111 veeCol(
new std::vector<rec::Vee>);
113 auto trkVeeAssns = std::make_unique< art::Assns<rec::Track, rec::Vee, rec::TrackEnd> >();
118 auto const&
tracks = *trackHandle;
123 size_t nTrack = tracks.size();
125 const float m_prot = 0.9382723;
126 const float m_pi = 0.13957018;
132 std::vector<TrackPar> tplist;
133 std::vector<TrackEnd> usebeg;
134 for (
size_t iTrack=0; iTrack<nTrack-1; ++iTrack)
136 Track thisTrack = tracks.at(iTrack);
137 for (
int iEnd=0; iEnd<2; ++iEnd)
141 if (!
goodTrack(thisTrack,iEnd))
continue;
142 TVector3 thisEndinSpace;
146 TVector3 tmpv(thisTrack.Vertex());
147 thisEndinSpace = tmpv;
148 TVector3 tmpv2(thisTrack.VtxDir());
149 tmpv2 *= thisTrack.Momentum_beg();
151 thisCharge = thisTrack.ChargeBeg();
155 TVector3 tmpv(thisTrack.End());
156 thisEndinSpace = tmpv;
157 TVector3 tmpv2(thisTrack.EndDir());
158 tmpv2 *= thisTrack.Momentum_end();
160 thisCharge = thisTrack.ChargeEnd();
163 tplist.emplace_back(thisTrack);
165 usebeg.push_back(iTrackEnd);
167 for (
size_t jTrack=iTrack+1; jTrack<nTrack; ++jTrack)
169 Track thatTrack = tracks.at(jTrack);
170 for (
int jEnd = 0; jEnd<2; ++jEnd)
174 if (!
goodTrack(thatTrack,jEnd))
continue;
175 TVector3 thatEndinSpace;
179 TVector3 tmpv(thatTrack.Vertex());
180 thatEndinSpace = tmpv;
181 TVector3 tmpv2(thatTrack.VtxDir());
182 tmpv2 *= thatTrack.Momentum_beg();
184 thatCharge = thatTrack.ChargeBeg();
188 TVector3 tmpv(thatTrack.End());
189 thatEndinSpace = tmpv;
190 TVector3 tmpv2(thatTrack.EndDir());
191 tmpv2 *= thatTrack.Momentum_end();
193 thatCharge = thatTrack.ChargeEnd();
199 if ( (thisEndinSpace-thatEndinSpace).Mag() >
fDCut )
continue;
201 if (
fPrintLevel>1) std::cout <<
"veefinder: opening angle: " << thisP.Angle(thatP) <<
std::endl;
202 float openang = thisP.Angle(thatP);
208 if (thisCharge == 0 || thatCharge == 0)
212 std::cout <<
"veefinder: One of the track charges is zero: " << thisCharge <<
" " << thatCharge <<
std::endl;
215 if (thisCharge == thatCharge)
continue;
220 tplist.emplace_back(thatTrack);
222 usebeg.push_back(jTrackEnd);
224 std::vector<float> dcas;
225 std::vector<float> vpos;
226 std::vector<std::vector<float>> covmat;
229 int iret =
fitVertex(tplist,vpos,chisquared,covmat,time,usebeg,dcas);
231 if (iret != 0)
continue;
233 if (dcas[0] >
fRCut || dcas[1] >
fRCut)
continue;
236 TVector3 sumP = thisP + thatP;
237 TVector3 vposv(vpos.data());
243 size_t otherTrackCount = 0;
244 for (
size_t kTrack=0; kTrack<nTrack; ++kTrack)
246 if (kTrack == iTrack || kTrack == jTrack)
continue;
247 Track otherTrack = tracks.at(kTrack);
248 for (
int kEnd = 0; kEnd<2; ++kEnd)
260 std::vector<TLorentzVector> fourmomentumvec;
261 float energy_kshort_hyp =
262 TMath::Sqrt(thisP.Mag2() + m_pi*
m_pi) +
263 TMath::Sqrt(thatP.Mag2() + m_pi*
m_pi);
269 float energy_lambda1_hyp =
270 TMath::Sqrt(thisP.Mag2() + m_prot*m_prot) +
271 TMath::Sqrt(thatP.Mag2() + m_pi*
m_pi);
272 float energy_lambda2_hyp =
273 TMath::Sqrt(thisP.Mag2() + m_pi*
m_pi) +
274 TMath::Sqrt(thatP.Mag2() + m_prot*m_prot);
275 if (thisCharge == -1)
277 float tmp = energy_lambda1_hyp;
278 energy_lambda1_hyp = energy_lambda2_hyp;
279 energy_lambda2_hyp =
tmp;
281 fourmomentumvec.emplace_back(sumP,energy_kshort_hyp);
282 fourmomentumvec.emplace_back(sumP,energy_lambda1_hyp);
283 fourmomentumvec.emplace_back(sumP,energy_lambda2_hyp);
286 for (
size_t icm=0; icm<3; ++icm)
288 for (
size_t jcm=0;jcm<3; ++jcm)
290 covmatvec[icounter++] = covmat.at(icm).at(jcm);
294 veeCol->emplace_back(vpos.data(),covmatvec,chisquared,time,fourmomentumvec.data());
295 auto const veeptr = veePtrMaker(veeCol->size()-1);
296 auto const itrackptr = trackPtrMaker( iTrack );
297 auto const jtrackptr = trackPtrMaker( jTrack );
298 trkVeeAssns->addSingle(itrackptr,veeptr,iTrackEnd);
299 trkVeeAssns->addSingle(jtrackptr,veeptr,jTrackEnd);
float fOpenAngMaxCut
opening angle cut in 3D
bool goodOtherTrack(Track &trk, TrackEnd usebeg, TVector3 &veepos, TVector3 &veep)
float fOpenAngMinCut
minimum opening angle cut in 3D
float fRCut
DCA cut from track to vee vertex.
int fPrintLevel
module print level
std::string fTrackLabel
track module label
float fDCut
distance cut from endpoints of tracks
int fitVertex(std::vector< TrackPar > &tracks, std::vector< float > &xyz, float &chisquared, std::vector< std::vector< float > > &covmat, double &time, std::vector< TrackEnd > &usebeg, std::vector< float > &dca)
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
bool goodTrack(Track &trk, TrackEnd usebeg)
bool fRequireOppositeCharge
whether or not to require opposite-sign tracks.
TrackEnd const TrackEndBeg
TrackCollectionProxyElement< TrackCollProxy > Track
Proxy to an element of a proxy collection of recob::Track objects.
QTextStream & endl(QTextStream &s)
size_t fRequireNearTrack
number of nearby tracks to require (primary vertex)