191 std::vector< art::Ptr< raw::ExternalTrigger> > trigs;
197 std::vector< art::Ptr< recob::Hit> > hits;
201 auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
206 for (
size_t CLoop=0; CLoop < trigs.size(); ++CLoop) {
207 for (
size_t CLoop2=0; CLoop2 < trigs.size(); ++CLoop2) {
208 if ( fabs( trigs[CLoop]->GetTrigTime() - trigs[CLoop2]->GetTrigTime()) >
fCoincidenceTolerance )
continue;
209 if ( CLoop == CLoop2 )
continue;
214 if ( ( trigs[CLoop]->GetTrigID() >= 6 && trigs[CLoop]->GetTrigID() <= 15 && trigs[CLoop2]->GetTrigID() >= 28 && trigs[CLoop2]->GetTrigID() <= 37 )
215 || ( trigs[CLoop]->GetTrigID() <= 5 && trigs[CLoop2]->GetTrigID() >= 22 && trigs[CLoop2]->GetTrigID() <= 27 )
216 || ( trigs[CLoop]->GetTrigID() >= 16 && trigs[CLoop]->GetTrigID() <= 21 && trigs[CLoop2]->GetTrigID() >= 38 && trigs[CLoop2]->GetTrigID() <= 43 )
218 if (
fDebug) std::cout <<
"I have a match..." 219 <<
"CLoop " << CLoop <<
", ID " << trigs[CLoop]->GetTrigID() <<
", time " << trigs[CLoop]->GetTrigTime() <<
"..." 220 <<
"CLoop2 " << CLoop2 <<
", ID " << trigs[CLoop2]->GetTrigID() <<
", Time " << trigs[CLoop2]->GetTrigTime()
236 <<
". They correspond to TrigIDs " << IDA <<
" and " << IDB <<
", and times " << TimeA <<
" (" << TickA <<
") and " << TimeB <<
" (" << TickB <<
")" 242 art::FindOneP<raw::RawDigit> ChannelHitRawDigits (HitListHandle, evt,
fHitsModuleLabel);
243 art::FindOneP<recob::Wire> ChannelHitWires (HitListHandle, evt,
fHitsModuleLabel);
248 ChannelHitWires.isValid(),
249 ChannelHitRawDigits.isValid()
253 std::vector < std::pair < std::vector < recob::Hit >,
size_t > > UnDisambigHits;
270 float VertexX[4] = { (float)it1->second.second[0][0], (
float)it1->second.second[1][0], (float)it2->second.second[1][0], (
float)it2->second.second[0][0]};
271 float VertexY[4] = { (float)it1->second.second[0][1], (
float)it1->second.second[2][1], (float)it2->second.second[2][1], (
float)it2->second.second[0][1]};
272 float VertexZ[4] = { (float)it1->second.second[0][2], (
float)it1->second.second[1][2], (float)it2->second.second[1][2], (
float)it2->second.second[0][2]};
279 if (
fDebug) std::cout <<
"\nCounter1, with TrigID " << trigs.at(
ExternalTrigIndexVec[TrigInd].first)->GetTrigID() <<
" (" << it1->first <<
") has corners......" 280 <<
"( " << it1->second.second[0][0] <<
", " << it1->second.second[0][1] <<
", " << it1->second.second[0][2] <<
") and " 281 <<
"( " << it1->second.second[1][0] <<
", " << it1->second.second[1][1] <<
", " << it1->second.second[1][2] <<
") and " 282 <<
"( " << it1->second.second[2][0] <<
", " << it1->second.second[2][1] <<
", " << it1->second.second[2][2] <<
") and " 283 <<
"( " << it1->second.second[3][0] <<
", " << it1->second.second[3][1] <<
", " << it1->second.second[3][2] <<
") " 284 <<
"\nCounter2, with TrigID " << trigs.at(
ExternalTrigIndexVec[TrigInd].second)->GetTrigID() <<
" (" << it2->first <<
") has corners......" 285 <<
"( " << it2->second.second[0][0] <<
", " << it2->second.second[0][1] <<
", " << it2->second.second[0][2] <<
") and " 286 <<
"( " << it2->second.second[1][0] <<
", " << it2->second.second[1][1] <<
", " << it2->second.second[1][2] <<
") and " 287 <<
"( " << it2->second.second[2][0] <<
", " << it2->second.second[2][1] <<
", " << it2->second.second[2][2] <<
") and " 288 <<
"( " << it2->second.second[3][0] <<
", " << it2->second.second[3][1] <<
", " << it2->second.second[3][2] <<
") " 289 <<
"\nVertexX has co-ords " << VertexX[0] <<
" " << VertexX[1] <<
" " << VertexX[2] <<
" " << VertexX[3]
290 <<
"\nVertexY has co-ords " << VertexY[0] <<
" " << VertexY[1] <<
" " << VertexY[2] <<
" " << VertexY[3]
291 <<
"\nVertexZ has co-ords " << VertexZ[0] <<
" " << VertexZ[1] <<
" " << VertexZ[2] <<
" " << VertexZ[3]
302 float c1x = it1->second.first.X();
303 float c1z = it1->second.first.Z();
304 float c2x = it2->second.first.X();
305 float c2z = it2->second.first.Z();
308 float slope = ((c1x-c2x)/(c1z-c2z));
314 else if (trignum == 112 || trignum == 113)
316 float slope = ((c1z-c2z)/(c1x-c2x));
324 doHitLineFitAlg =
false;
328 std::vector<dune::HitLineFitAlg::HitLineFitData> fitdata;
330 std::map<unsigned int, size_t> index_convert;
331 unsigned int fitdataindex = 0;
336 for(
size_t HitInd = 0; HitInd < hits.size(); HitInd++ ) {
337 if ( hits[HitInd]->PeakTime()-
TrigTime < 0 )
continue;
338 if( hits[HitInd]->View() !=
geo::kZ )
continue;
341 bool KeepHit =
false;
345 double WireStart[3], WireEnd[3];
347 float DriftDist = (0.5 *(hits[HitInd]->PeakTime()-
TrigTime) * fDriftVelocity );
348 if ( hits[HitInd]->WireID().TPC % 2 == 0) {
349 DriftDist = -DriftDist;
351 float HitXPos = DriftDist + WireStart[0];
353 KeepHit =
pnpoly( 4, VertexX, VertexZ, HitXPos, (
float)WireEnd[2] );
371 else if (trignum == 112 || trignum == 113)
374 hlfd.
hitVert = (float)WireEnd[2];
381 fitdata.push_back(hlfd);
382 index_convert[fitdataindex] = HitInd;
392 if (!doHitLineFitAlg)
394 for (
unsigned int i_fd = 0; i_fd < fitdata.size(); ++i_fd)
398 hcol.emplace_back(*hits[index_convert[i_fd]],wire,rawdigits);
401 else if (retval != 1 && doHitLineFitAlg)
403 std::cout <<
"No line fit could be found, or not enough hits in counter shadow to make a line" <<
std::endl;
409 for (
unsigned int i_fd = 0; i_fd < fitdata.size(); ++i_fd)
411 if (fitdata[i_fd].hitREAL)
415 hcol.emplace_back(*hits[index_convert[i_fd]],wire,rawdigits);
420 std::cout <<
"After collection wires for TrigInd " << TrigInd <<
", hcol has size " << hcol.size() <<
std::endl;
425 for(
size_t HitInd = 0; HitInd < hits.size(); HitInd++ ) {
426 if ( hits[HitInd]->PeakTime()-
TrigTime < 0 )
continue;
427 if ( hits[HitInd]->View() ==
geo::kZ )
continue;
432 float DriftDist = (0.5 *(hits[HitInd]->PeakTime()-
TrigTime) * fDriftVelocity );
433 float HitPosShortDrift =
UVStartEndPointsX[ hits[HitInd]->View()*2 ] - DriftDist;
434 float HitPosLongDrift =
UVStartEndPointsX[ 1 + hits[HitInd]->View()*2 ] + DriftDist;
436 bool CouldBeShortDrift =
true;
437 bool CouldBeLongDrift =
true;
438 if ( HitPosShortDrift < SmallestX ) CouldBeShortDrift =
false;
439 if ( HitPosLongDrift > BiggestX ) CouldBeLongDrift =
false;
440 if ( !CouldBeShortDrift && !CouldBeLongDrift )
continue;
451 std::vector< recob::Hit > TempHitsVec;
452 int Channel = hits[HitInd]->Channel();
454 for (
unsigned int ch=0; ch<WireList.size(); ++ch) {
455 if ( !CouldBeShortDrift && (WireList[ch].TPC%2 == 0) )
continue;
456 if ( !CouldBeLongDrift && (WireList[ch].TPC%2 != 0) )
continue;
458 double WStart[3], WEnd[3];
465 for (
int WPass=0; WPass < NPoints; ++WPass) {
466 double ThisXPos = HitPosLongDrift;
467 if (WireList[ch].TPC%2 == 0) ThisXPos = HitPosShortDrift;
468 double ThisYPos = WStart[1] + ( WPass * (WEnd[1]-WStart[1])/(NPoints-1) );
469 double ThisZPos = WStart[2] + ( WPass * (WEnd[2]-WStart[2])/(NPoints-1) );
471 bool GoodWire =
false;
473 GoodWire =
pnpoly ( 4, VertexX, VertexZ, ThisXPos, ThisZPos );
474 if (!GoodWire)
continue;
479 GoodWire =
pnpoly ( 4, VertexY, VertexZ, ThisYPos, ThisZPos );
481 GoodWire =
pnpoly ( 4, VertexX, VertexY, ThisXPos, ThisYPos );
483 if (!GoodWire)
continue;
486 if ( WireList[ch].
Wire == hits[HitInd]->WireID().
Wire && WireList[ch].TPC == hits[HitInd]->WireID().TPC ) {
487 TempHitsVec.push_back( *hits[HitInd] );
489 recob::Hit hit( hits[HitInd]->
Channel(), hits[HitInd]->StartTick(), hits[HitInd]->EndTick(), hits[HitInd]->PeakTime(),
490 hits[HitInd]->SigmaPeakTime(), hits[HitInd]->RMS(), hits[HitInd]->PeakAmplitude(),
491 hits[HitInd]->SigmaPeakAmplitude(), hits[HitInd]->SummedADC(), hits[HitInd]->Integral(),
492 hits[HitInd]->SigmaIntegral(), hits[HitInd]->Multiplicity(), hits[HitInd]->LocalIndex(),
493 hits[HitInd]->GoodnessOfFit(), hits[HitInd]->DegreesOfFreedom(), hits[HitInd]->View(),
494 hits[HitInd]->SignalType(), WireList[ch]
496 TempHitsVec.push_back(
hit );
501 if (TempHitsVec.size() == 0)
continue;
502 if (TempHitsVec.size() == 1) {
505 hcol.emplace_back(TempHitsVec[0], wire, rawdigits);
507 UnDisambigHits.push_back( std::make_pair(TempHitsVec, HitInd) );
510 std::cout <<
"After Induction planes for TrigInd " << TrigInd <<
" hcol now has size " << hcol.size() <<
std::endl;
512 std::cout <<
"After all trigger indexes the vectors have sizes: hits -> " << hits.size() <<
", hcol -> " << hcol.size() <<
", UnDisambigHits -> " << UnDisambigHits.size() <<
std::endl;
515 std::vector<recob::Hit>
const &Peak = hcol.peek();
516 int TPC1, TPC2, TPC3, TPC4, TPC5, TPC6, TPC7, TPC0;
517 TPC1 = TPC2 = TPC3 = TPC4 = TPC5 = TPC6 = TPC7 = TPC0 = 0;
518 for (
size_t qq=0; qq<Peak.size(); ++qq) {
519 if (Peak[qq].WireID().TPC == 0) ++TPC0;
520 else if (Peak[qq].WireID().TPC == 1) ++TPC1;
521 else if (Peak[qq].WireID().TPC == 2) ++TPC2;
522 else if (Peak[qq].WireID().TPC == 3) ++TPC3;
523 else if (Peak[qq].WireID().TPC == 4) ++TPC4;
524 else if (Peak[qq].WireID().TPC == 5) ++TPC5;
525 else if (Peak[qq].WireID().TPC == 6) ++TPC6;
526 else if (Peak[qq].WireID().TPC == 7) ++TPC7;
528 if (
fDebug) std::cout <<
"\n\nI have these hits in each TPC " << TPC0 <<
" " << TPC1 <<
" " << TPC2 <<
" " << TPC3 <<
" " << TPC4 <<
" " << TPC5 <<
" " << TPC6 <<
" " << TPC7 <<
"\n\n" <<
std::endl;
533 std::vector < std::pair < std::vector <recob::Hit>,
size_t > > StillUnDisambigHits;
534 std::vector < std::pair < recob::Hit, size_t > > NowGoodHits;
537 std::cout <<
"\nNow to do my next step....fit hits to a 2D line in XZ and find which ambiguous hit is closer..." <<
std::endl;
538 std::vector < recob::Hit >
const &NextGoodHits = hcol.peek();
540 TwoDimXZLineFit( NextGoodHits, UnDisambigHits, StillUnDisambigHits, NowGoodHits );
541 for (
size_t NowDisambig=0; NowDisambig < NowGoodHits.size(); ++NowDisambig ) {
542 size_t WhichRawHit = NowGoodHits[NowDisambig].second;
545 hcol.emplace_back(NowGoodHits[NowDisambig].first, wire, rawdigits);
547 OutAndClearVector( UnDisambigHits, StillUnDisambigHits, NowGoodHits, hcol.size() );
551 std::cout <<
"\nNow to do my next step....fit hits to a 2D line in each TPC / Plane combination and find which ambiguous hit is closer..." <<
std::endl;
552 std::vector < recob::Hit >
const &NewerGoodHits = hcol.peek();
554 TwoDimPlaneLineFit( NewerGoodHits, UnDisambigHits, StillUnDisambigHits, NowGoodHits );
555 for (
size_t NowDisambig=0; NowDisambig < NowGoodHits.size(); ++NowDisambig ) {
556 size_t WhichRawHit = NowGoodHits[NowDisambig].second;
559 hcol.emplace_back(NowGoodHits[NowDisambig].first, wire, rawdigits);
561 OutAndClearVector( UnDisambigHits, StillUnDisambigHits, NowGoodHits, hcol.size() );
565 std::cout <<
"\nNow to do my next step....if induction wire crosses a collection wire with a good hit..." <<
std::endl;
566 std::vector < recob::Hit >
const &NewHits = hcol.peek();
568 CrossCollection( NewHits, UnDisambigHits, StillUnDisambigHits, NowGoodHits );
569 for (
size_t NowDisambig=0; NowDisambig < NowGoodHits.size(); ++NowDisambig ) {
570 size_t WhichRawHit = NowGoodHits[NowDisambig].second;
573 hcol.emplace_back(NowGoodHits[NowDisambig].first, wire, rawdigits);
575 OutAndClearVector( UnDisambigHits, StillUnDisambigHits, NowGoodHits, hcol.size() );
579 std::cout <<
"\nNow to do my next step....if any adjacent wires have hits on them..." <<
std::endl;
580 std::vector < recob::Hit >
const &NewGoodHits = hcol.peek();
582 AdjacentWireWidth( NewGoodHits, UnDisambigHits, StillUnDisambigHits, NowGoodHits );
583 for (
size_t NowDisambig=0; NowDisambig < NowGoodHits.size(); ++NowDisambig ) {
584 size_t WhichRawHit = NowGoodHits[NowDisambig].second;
587 hcol.emplace_back(NowGoodHits[NowDisambig].first, wire, rawdigits);
589 OutAndClearVector( UnDisambigHits, StillUnDisambigHits, NowGoodHits, hcol.size() );
592 std::cout <<
"\nAfter all that the vectors have sizes: hits -> " << hits.size() <<
", hcol -> " << hcol.size() <<
std::endl;
EventNumber_t event() const
int pnpoly(int nvert, float *vertx, float *verty, float testx, float testy)
void CrossCollection(std::vector< recob::Hit > const GoodHits, std::vector< std::pair< std::vector< recob::Hit >, size_t > > &InputBadVector, std::vector< std::pair< std::vector< recob::Hit >, size_t > > &OutputBadVector, std::vector< std::pair< recob::Hit, size_t > > &OutputGoodVector)
int FitLine(std::vector< HitLineFitData > &data, HitLineFitResults &bestfit)
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Planes which measure Z direction.
void TwoDimXZLineFit(std::vector< recob::Hit > const GoodHits, std::vector< std::pair< std::vector< recob::Hit >, size_t > > &InputBadVector, std::vector< std::pair< std::vector< recob::Hit >, size_t > > &OutputBadVector, std::vector< std::pair< recob::Hit, size_t > > &OutputGoodVector)
std::map< unsigned int, std::pair< TVector3, std::vector< TVector3 > > > CounterPositionMap
std::vector< std::pair< unsigned int, unsigned int > > ExternalTrigIndexVec
void SetSeed(UInt_t seed)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
A class handling a collection of hits and its associations.
std::string fHitsModuleLabel
void SetParameter(int i, double startValue, double minValue, double maxValue)
void OutAndClearVector(std::vector< std::pair< std::vector< recob::Hit >, size_t > > &InitialBad, std::vector< std::pair< std::vector< recob::Hit >, size_t > > &NewBad, std::vector< std::pair< recob::Hit, size_t > > &NewGood, size_t HitsSize)
static int max(int a, int b)
double fConvCountTimeToTicks
Detector simulation of raw signals on wires.
void AdjacentWireWidth(std::vector< recob::Hit > const GoodHits, std::vector< std::pair< std::vector< recob::Hit >, size_t > > &InputBadVector, std::vector< std::pair< std::vector< recob::Hit >, size_t > > &OutputBadVector, std::vector< std::pair< recob::Hit, size_t > > &OutputGoodVector)
int fCoincidenceTolerance
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
dune::HitLineFitAlg fFitAlg
ChannelMappingService::Channel Channel
void TwoDimPlaneLineFit(std::vector< recob::Hit > const GoodHits, std::vector< std::pair< std::vector< recob::Hit >, size_t > > &InputBadVector, std::vector< std::pair< std::vector< recob::Hit >, size_t > > &OutputBadVector, std::vector< std::pair< recob::Hit, size_t > > &OutputGoodVector)
double UVStartEndPointsX[4]
void SetHorizVertRanges(float hmin, float hmax, float vmin, float vmax)
void WireEndPoints(geo::WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
CLHEP::HepRandomEngine & fEngine
2D representation of charge deposited in the TDC/wire plane
art::ServiceHandle< geo::Geometry > fGeom
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
second_as<> second
Type of time stored in seconds, in double precision.
QTextStream & endl(QTextStream &s)
std::string fCounterModuleLabel