129 double Efield_drift = detProp.Efield();
130 double Temperature = detProp.Temperature();
133 double driftvelocity = detProp.DriftVelocity(Efield_drift, Temperature);
136 double timepitch = driftvelocity * timetick;
142 for (
unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
150 std::unique_ptr<std::vector<recob::Vertex>> vcol(
new std::vector<recob::Vertex>);
151 std::unique_ptr<std::vector<recob::EndPoint2D>> epcol(
152 new std::vector<recob::EndPoint2D>);
153 std::unique_ptr<art::Assns<recob::EndPoint2D, recob::Hit>> assnep(
155 std::unique_ptr<art::Assns<recob::Vertex, recob::Shower>> assnsh(
157 std::unique_ptr<art::Assns<recob::Vertex, recob::Track>> assntr(
159 std::unique_ptr<art::Assns<recob::Vertex, recob::Hit>> assnh(
164 int nplanes = geom->
Views().size();
166 std::vector<std::vector<int>> Cls(nplanes);
167 std::vector<std::vector<CluLen>> clulens(nplanes);
169 std::vector<double> dtdwstart;
172 for (
size_t iclu = 0; iclu < clusters.
size(); ++iclu) {
174 float w0 = clusters[iclu]->StartWire();
175 float w1 = clusters[iclu]->EndWire();
176 float t0 = clusters[iclu]->StartTick();
177 float t1 = clusters[iclu]->EndTick();
181 clulen.length = std::hypot((w0 - w1) * wire_pitch,
182 detProp.ConvertTicksToX(t0, clusters[iclu]->View(), 0, 0) -
183 detProp.ConvertTicksToX(t1, clusters[iclu]->View(), 0, 0));
185 switch (clusters[iclu]->View()) {
187 case geo::kU: clulens[0].push_back(clulen);
break;
188 case geo::kV: clulens[1].push_back(clulen);
break;
189 case geo::kZ: clulens[2].push_back(clulen);
break;
193 std::vector<double> wires;
194 std::vector<double> times;
196 std::vector<art::Ptr<recob::Hit>>
hit = fmh.at(iclu);
197 std::sort(hit.begin(), hit.end(),
SortByWire());
199 for (
size_t i = 0; i < hit.size(); ++i) {
200 wires.push_back(hit[i]->
WireID().Wire);
201 times.push_back(hit[i]->PeakTime());
205 TGraph* the2Dtrack =
new TGraph(std::min(10, n), &wires[0], ×[0]);
207 the2Dtrack->Fit(
"pol1",
"Q");
208 TF1* pol1 = (TF1*)the2Dtrack->GetFunction(
"pol1");
210 pol1->GetParameters(par);
211 dtdwstart.push_back(par[1]);
216 dtdwstart.push_back(std::tan(clusters[iclu]->StartAngle()));
222 dtdwstart.push_back(std::tan(clusters[iclu]->StartAngle()));
226 for (
size_t i = 0; i < clulens.size(); ++i) {
227 std::sort(clulens[i].
begin(), clulens[i].
end(), myfunction);
228 for (
size_t j = 0; j < clulens[i].size(); ++j) {
229 Cls[i].push_back(clulens[i][j].index);
233 std::vector<std::vector<int>> cluvtx(nplanes);
234 std::vector<double> vtx_w;
235 std::vector<double> vtx_t;
237 for (
int i = 0; i < nplanes; ++i) {
238 if (Cls[i].
size() >= 1) {
254 for (
unsigned j = 0; j < Cls[i].size(); ++j) {
255 double lclu = std::sqrt(
256 pow((clusters[Cls[i][j]]->StartWire() - clusters[Cls[i][j]]->EndWire()) * 13.5, 2) +
257 pow(clusters[Cls[i][j]]->StartTick() - clusters[Cls[i][j]]->EndTick(), 2));
259 bool deltaraylike =
false;
260 bool enoughhits =
false;
262 double wb = clusters[Cls[i][j]]->StartWire();
263 double we = clusters[Cls[i][j]]->EndWire();
264 double tt = clusters[Cls[i][j]]->StartTick();
265 double dtdw = dtdwstart[Cls[i][j]];
266 int nhits = fmh.at(Cls[i][j]).size();
267 ww0 = (tt - tt1 + dtdw1 * wb1 - dtdw * wb) / (dtdw1 - dtdw);
269 if ((!rev && ww0 > wb1 + 15) || (rev && ww0 < we1 - 15)) deltaraylike =
true;
270 if (((!rev && ww0 > wb1 + 10) || (rev && ww0 < we1 - 10)) && nhits < 5)
272 if (wb > wb1 + 20 && nhits < 20) deltaraylike =
true;
273 if (wb > wb1 + 50 && nhits < 20) deltaraylike =
true;
274 if (wb > wb1 + 8 && TMath::Abs(dtdw1 - dtdw) < 0.15) deltaraylike =
true;
275 if (
std::abs(wb - wb1) > 30 &&
std::abs(we - we1) > 30) deltaraylike =
true;
280 double alpha = std::atan(dtdw);
281 if (nhits >=
int(2 + 3 * (1 -
std::abs(std::cos(alpha))))) enoughhits =
true;
282 if (nhits < 5 && (ww0 < wb1 - 20 || ww0 > we1 + 20)) enoughhits =
false;
286 if (c1 != -1 && c2 != -1) {
287 double wb = clusters[Cls[i][j]]->StartWire();
288 double we = clusters[Cls[i][j]]->EndWire();
289 ww0 = (tt2 - tt1 + dtdw1 * wb1 - dtdw2 * wb2) / (dtdw1 - dtdw2);
296 if (c1 != -1 && !deltaraylike && enoughhits) {
306 wb1 = clusters[Cls[i][j]]->StartWire();
307 we1 = clusters[Cls[i][j]]->EndWire();
308 tt1 = clusters[Cls[i][j]]->StartTick();
310 wb1 = clusters[Cls[i][j]]->EndWire();
311 we1 = clusters[Cls[i][j]]->StartWire();
312 tt1 = clusters[Cls[i][j]]->EndTick();
314 dtdw1 = dtdwstart[Cls[i][j]];
316 else if (lclu2 < lclu) {
317 if (!deltaraylike && enoughhits && replace) {
320 wb2 = clusters[Cls[i][j]]->StartWire();
321 we2 = clusters[Cls[i][j]]->EndWire();
322 tt2 = clusters[Cls[i][j]]->StartTick();
323 dtdw2 = dtdwstart[Cls[i][j]];
327 if (c1 != -1 && c2 != -1) {
328 cluvtx[i].push_back(c1);
329 cluvtx[i].push_back(c2);
331 double w1 = clusters[
c1]->StartWire();
332 double t1 = clusters[
c1]->StartTick();
333 if (clusters[c1]->StartWire() > clusters[c1]->EndWire()) {
334 w1 = clusters[
c1]->EndWire();
335 t1 = clusters[
c1]->EndTick();
337 double k1 = dtdwstart[
c1];
338 double w2 = clusters[
c2]->StartWire();
339 double t2 = clusters[
c2]->StartTick();
340 if (clusters[c2]->StartWire() > clusters[c2]->EndWire()) {
341 w1 = clusters[
c2]->EndWire();
342 t1 = clusters[
c2]->EndTick();
344 double k2 = dtdwstart[
c2];
351 double t0 = (k1 * k2 * (w1 - w2) + k1 * t2 - k2 * t1) / (k1 - k2);
352 double w0 = (t2 - t1 + k1 * w1 - k2 * w2) / (k1 - k2);
357 else if (Cls[i].
size() >= 1) {
359 cluvtx[i].push_back(c1);
360 vtx_w.push_back(wb1);
361 vtx_t.push_back(tt1);
364 cluvtx[i].push_back(Cls[i][0]);
365 vtx_w.push_back(clusters[Cls[i][0]]->StartWire());
366 vtx_t.push_back(clusters[Cls[i][0]]->StartTick());
373 std::vector<art::Ptr<recob::Hit>>
hits = fmh.at(Cls[i][0]);
375 for (
size_t h = 0; h < hits.size(); ++h)
376 totalQ += hits[h]->Integral();
381 (
unsigned int)vtx_w.back());
387 clusters[Cls[i][0]]->View(),
401 Double_t vtxcoord[3];
402 if (Cls[0].
size() > 0 && Cls[1].
size() > 0) {
403 double Iw0 = (vtx_w[0] + 3.95) * wire_pitch;
404 double Cw0 = (vtx_w[1] + 1.84) * wire_pitch;
406 double It0 = vtx_t[0] - presamplings;
408 double Ct0 = vtx_t[1] - presamplings;
410 vtxcoord[0] = detProp.ConvertTicksToX(vtx_t[1], 1, 0, 0);
411 vtxcoord[1] = (Cw0 - Iw0) / (2. * std::sin(Angle));
412 vtxcoord[2] = (Cw0 + Iw0) / (2. * std::cos(Angle)) - YC / 2. * std::tan(Angle);
415 if (vtx_w[0] >= 0 && vtx_w[0] <= 239 && vtx_w[1] >= 0 && vtx_w[1] <= 239) {
427 vtxcoord[0] = -99999;
428 vtxcoord[1] = -99999;
429 vtxcoord[2] = -99999;
432 dtIC->Fill(It0 - Ct0);
435 vtxcoord[0] = -99999;
436 vtxcoord[1] = -99999;
437 vtxcoord[2] = -99999;
446 vcol->push_back(the3Dvertex);
452 MF_LOG_VERBATIM(
"Summary") << std::setfill(
'-') << std::setw(175) <<
"-" << std::setfill(
' ');
454 for (
size_t i = 0; i < epcol->size(); ++i)
456 for (
size_t i = 0; i < vcol->size(); ++i)
459 evt.
put(std::move(epcol));
460 evt.
put(std::move(vcol));
461 evt.
put(std::move(assnep));
462 evt.
put(std::move(assntr));
463 evt.
put(std::move(assnsh));
464 evt.
put(std::move(assnh));
code to link reconstructed objects back to the MC truth information
std::string GetLArTPCVolumeName(TPCID const &tpcid=tpc_zero) const
Return the name of specified LAr TPC volume.
std::string fClusterModuleLabel
The data type to uniquely identify a Plane.
constexpr auto abs(T v)
Returns the absolute value of the argument.
bool ChannelsIntersect(raw::ChannelID_t c1, raw::ChannelID_t c2, double &y, double &z) const
Returns an intersection point of two channels.
Planes which measure Z direction.
Definition of vertex object for LArSoft.
double ThetaZ() const
Angle of the wires from positive z axis; .
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
bool SortByWire(art::Ptr< recob::Hit > const &h1, art::Ptr< recob::Hit > const &h2)
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void push_back(Ptr< U > const &p)
Detector simulation of raw signals on wires.
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.
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
#define MF_LOG_VERBATIM(category)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
int trigger_offset(DetectorClocksData const &data)
std::set< View_t > const & Views() const
Returns a list of possible views in the detector.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Length_t DetHalfHeight(TPCID const &tpcid=tpc_zero) const
Returns the half height of the active volume of the specified TPC.
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
recob::tracking::Plane Plane
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
1.8.11 
