118 auto const& tpc = geom->
TPC({0, 0});
122 double Angle = wireReadoutGeom.Plane({tpc.ID(), 1}).Wire(0).ThetaZ(
false) - TMath::Pi() / 2.;
129 double wire_pitch = wireReadoutGeom.Plane({0, 0, 0}).WirePitch();
130 double Efield_drift = detProp.Efield();
131 double Temperature = detProp.Temperature();
134 double driftvelocity = detProp.DriftVelocity(Efield_drift, Temperature);
137 double timepitch = driftvelocity * timetick;
143 for (
unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
151 std::unique_ptr<std::vector<recob::Vertex>> vcol(
new std::vector<recob::Vertex>);
152 std::unique_ptr<std::vector<recob::EndPoint2D>> epcol(
153 new std::vector<recob::EndPoint2D>);
154 std::unique_ptr<art::Assns<recob::EndPoint2D, recob::Hit>> assnep(
156 std::unique_ptr<art::Assns<recob::Vertex, recob::Shower>> assnsh(
158 std::unique_ptr<art::Assns<recob::Vertex, recob::Track>> assntr(
160 std::unique_ptr<art::Assns<recob::Vertex, recob::Hit>> assnh(
165 int nplanes = wireReadoutGeom.Views().size();
167 std::vector<std::vector<int>> Cls(nplanes);
168 std::vector<std::vector<CluLen>> clulens(nplanes);
170 std::vector<double> dtdwstart;
173 for (
size_t iclu = 0; iclu < clusters.
size(); ++iclu) {
175 float w0 = clusters[iclu]->StartWire();
176 float w1 = clusters[iclu]->EndWire();
177 float t0 = clusters[iclu]->StartTick();
178 float t1 = clusters[iclu]->EndTick();
182 clulen.length = std::hypot((w0 - w1) * wire_pitch,
183 detProp.ConvertTicksToX(t0, clusters[iclu]->View(), 0, 0) -
184 detProp.ConvertTicksToX(t1, clusters[iclu]->View(), 0, 0));
186 switch (clusters[iclu]->View()) {
188 case geo::kU: clulens[0].push_back(clulen);
break;
189 case geo::kV: clulens[1].push_back(clulen);
break;
190 case geo::kZ: clulens[2].push_back(clulen);
break;
194 std::vector<double> wires;
195 std::vector<double> times;
197 std::vector<art::Ptr<recob::Hit>>
hit = fmh.at(iclu);
198 std::sort(hit.begin(), hit.end(),
SortByWire());
200 for (
size_t i = 0; i < hit.size(); ++i) {
201 wires.push_back(hit[i]->
WireID().Wire);
202 times.push_back(hit[i]->PeakTime());
206 TGraph* the2Dtrack =
new TGraph(std::min(10, n), &wires[0], ×[0]);
208 the2Dtrack->Fit(
"pol1",
"Q");
209 TF1* pol1 = (TF1*)the2Dtrack->GetFunction(
"pol1");
211 pol1->GetParameters(par);
212 dtdwstart.push_back(par[1]);
217 dtdwstart.push_back(std::tan(clusters[iclu]->StartAngle()));
223 dtdwstart.push_back(std::tan(clusters[iclu]->StartAngle()));
227 for (
size_t i = 0; i < clulens.size(); ++i) {
228 std::sort(clulens[i].
begin(), clulens[i].
end(), myfunction);
229 for (
size_t j = 0; j < clulens[i].size(); ++j) {
230 Cls[i].push_back(clulens[i][j].index);
234 std::vector<std::vector<int>> cluvtx(nplanes);
235 std::vector<double> vtx_w;
236 std::vector<double> vtx_t;
238 for (
int i = 0; i < nplanes; ++i) {
239 if (Cls[i].
size() >= 1) {
255 for (
unsigned j = 0; j < Cls[i].size(); ++j) {
256 double lclu = std::sqrt(
257 pow((clusters[Cls[i][j]]->StartWire() - clusters[Cls[i][j]]->EndWire()) * 13.5, 2) +
258 pow(clusters[Cls[i][j]]->StartTick() - clusters[Cls[i][j]]->EndTick(), 2));
260 bool deltaraylike =
false;
261 bool enoughhits =
false;
263 double wb = clusters[Cls[i][j]]->StartWire();
264 double we = clusters[Cls[i][j]]->EndWire();
265 double tt = clusters[Cls[i][j]]->StartTick();
266 double dtdw = dtdwstart[Cls[i][j]];
267 int nhits = fmh.at(Cls[i][j]).size();
268 ww0 = (tt - tt1 + dtdw1 * wb1 - dtdw * wb) / (dtdw1 - dtdw);
270 if ((!rev && ww0 > wb1 + 15) || (rev && ww0 < we1 - 15)) deltaraylike =
true;
271 if (((!rev && ww0 > wb1 + 10) || (rev && ww0 < we1 - 10)) && nhits < 5)
273 if (wb > wb1 + 20 && nhits < 20) deltaraylike =
true;
274 if (wb > wb1 + 50 && nhits < 20) deltaraylike =
true;
275 if (wb > wb1 + 8 && TMath::Abs(dtdw1 - dtdw) < 0.15) deltaraylike =
true;
276 if (
std::abs(wb - wb1) > 30 &&
std::abs(we - we1) > 30) deltaraylike =
true;
281 double alpha = std::atan(dtdw);
282 if (nhits >=
int(2 + 3 * (1 -
std::abs(std::cos(alpha))))) enoughhits =
true;
283 if (nhits < 5 && (ww0 < wb1 - 20 || ww0 > we1 + 20)) enoughhits =
false;
287 if (c1 != -1 && c2 != -1) {
288 double wb = clusters[Cls[i][j]]->StartWire();
289 double we = clusters[Cls[i][j]]->EndWire();
290 ww0 = (tt2 - tt1 + dtdw1 * wb1 - dtdw2 * wb2) / (dtdw1 - dtdw2);
297 if (c1 != -1 && !deltaraylike && enoughhits) {
307 wb1 = clusters[Cls[i][j]]->StartWire();
308 we1 = clusters[Cls[i][j]]->EndWire();
309 tt1 = clusters[Cls[i][j]]->StartTick();
311 wb1 = clusters[Cls[i][j]]->EndWire();
312 we1 = clusters[Cls[i][j]]->StartWire();
313 tt1 = clusters[Cls[i][j]]->EndTick();
315 dtdw1 = dtdwstart[Cls[i][j]];
317 else if (lclu2 < lclu) {
318 if (!deltaraylike && enoughhits && replace) {
321 wb2 = clusters[Cls[i][j]]->StartWire();
322 we2 = clusters[Cls[i][j]]->EndWire();
323 tt2 = clusters[Cls[i][j]]->StartTick();
324 dtdw2 = dtdwstart[Cls[i][j]];
328 if (c1 != -1 && c2 != -1) {
329 cluvtx[i].push_back(c1);
330 cluvtx[i].push_back(c2);
332 double w1 = clusters[
c1]->StartWire();
333 double t1 = clusters[
c1]->StartTick();
334 if (clusters[c1]->StartWire() > clusters[c1]->EndWire()) {
335 w1 = clusters[
c1]->EndWire();
336 t1 = clusters[
c1]->EndTick();
338 double k1 = dtdwstart[
c1];
339 double w2 = clusters[
c2]->StartWire();
340 double t2 = clusters[
c2]->StartTick();
341 if (clusters[c2]->StartWire() > clusters[c2]->EndWire()) {
342 w1 = clusters[
c2]->EndWire();
343 t1 = clusters[
c2]->EndTick();
345 double k2 = dtdwstart[
c2];
352 double t0 = (k1 * k2 * (w1 - w2) + k1 * t2 - k2 * t1) / (k1 - k2);
353 double w0 = (t2 - t1 + k1 * w1 - k2 * w2) / (k1 - k2);
358 else if (Cls[i].
size() >= 1) {
360 cluvtx[i].push_back(c1);
361 vtx_w.push_back(wb1);
362 vtx_t.push_back(tt1);
365 cluvtx[i].push_back(Cls[i][0]);
366 vtx_w.push_back(clusters[Cls[i][0]]->StartWire());
367 vtx_t.push_back(clusters[Cls[i][0]]->StartTick());
374 std::vector<art::Ptr<recob::Hit>>
hits = fmh.at(Cls[i][0]);
376 for (
size_t h = 0; h < hits.size(); ++h)
377 totalQ += hits[h]->Integral();
380 (
unsigned int)vtx_w.back());
386 clusters[Cls[i][0]]->View(),
399 Double_t vtxcoord[3];
400 if (Cls[0].
size() > 0 && Cls[1].
size() > 0) {
401 double Iw0 = (vtx_w[0] + 3.95) * wire_pitch;
402 double Cw0 = (vtx_w[1] + 1.84) * wire_pitch;
404 double It0 = vtx_t[0] - presamplings;
406 double Ct0 = vtx_t[1] - presamplings;
408 vtxcoord[0] = detProp.ConvertTicksToX(vtx_t[1], 1, 0, 0);
409 vtxcoord[1] = (Cw0 - Iw0) / (2. * std::sin(Angle));
410 vtxcoord[2] = (Cw0 + Iw0) / (2. * std::cos(Angle)) - YC / 2. * std::tan(Angle);
412 if (vtx_w[0] >= 0 && vtx_w[0] <= 239 && vtx_w[1] >= 0 && vtx_w[1] <= 239) {
413 if (
auto intersection = wireReadoutGeom.ChannelsIntersect(
414 wireReadoutGeom.PlaneWireToChannel(
415 geo::WireID(0, 0, 0, (
int)((Iw0 / wire_pitch) - 3.95))),
416 wireReadoutGeom.PlaneWireToChannel(
417 geo::WireID(0, 0, 1, (
int)((Cw0 / wire_pitch) - 1.84))))) {
420 vtxcoord[1] = intersection->y;
421 vtxcoord[2] = intersection->z;
424 vtxcoord[0] = -99999;
425 vtxcoord[1] = -99999;
426 vtxcoord[2] = -99999;
429 dtIC->Fill(It0 - Ct0);
432 vtxcoord[0] = -99999;
433 vtxcoord[1] = -99999;
434 vtxcoord[2] = -99999;
443 vcol->push_back(the3Dvertex);
449 MF_LOG_VERBATIM(
"Summary") << std::setfill(
'-') << std::setw(175) <<
"-" << std::setfill(
' ');
451 for (
size_t i = 0; i < epcol->size(); ++i)
453 for (
size_t i = 0; i < vcol->size(); ++i)
456 evt.
put(std::move(epcol));
457 evt.
put(std::move(vcol));
458 evt.
put(std::move(assnep));
459 evt.
put(std::move(assntr));
460 evt.
put(std::move(assnsh));
461 evt.
put(std::move(assnh));
code to link reconstructed objects back to the MC truth information
std::string fClusterModuleLabel
constexpr auto abs(T v)
Returns the absolute value of the argument.
Planes which measure Z direction.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition of vertex object for LArSoft.
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.
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.
double HalfHeight() const
Height is associated with y coordinate [cm].
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.
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)
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
1.8.11 
