269 std::vector<bool> isClust;
270 for (
int i = 0; i <
int(ohits.size()); i++) isClust.push_back(
false);
272 std::vector<int> neigh;
274 for (
int i = 0; i <
int(ohits.size()); i++){
276 if (isClust[i])
continue;
277 neigh.erase(neigh.begin(),neigh.end());
280 for (
int j = i-1; j > 0; j--){
281 if (
Dist(ohits[i],ohits[j]) <
fR0 && !isClust[j] ){
284 if (
abs(ohits[i]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2)
break;
286 for (
int j = i+1; j <
int(ohits.size()); j++){
287 if (
Dist(ohits[i],ohits[j]) <
fR0 && !isClust[j] ){
290 if (
abs(ohits[i]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2)
break;
292 if (
int(neigh.size())<
fMinN)
continue;
293 neigh.erase(neigh.begin(),neigh.end());
296 std::vector<int> cands;
297 for (
int j = i; j <
int(ohits.size()); j++){
298 if (isClust[j])
continue;
299 if (ohits[j]->PeakTimeAbs()-ohits[i]->PeakTimeAbs()>
fBreakTime)
break;
307 neigh.push_back(centroidIdx);
309 std::vector<int> curN;
310 curN.push_back(centroidIdx);
313 for (
int j = centroidIdx-1; j > 0; j--){
314 if (
Dist(ohits[centroidIdx],ohits[j]) <
fR0 && !isClust[j] ){
318 if (
abs(ohits[centroidIdx]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2)
break;
320 for (
int j = centroidIdx+1; j <
int(ohits.size()); j++){
321 if (
Dist(ohits[centroidIdx],ohits[j]) <
fR0 && !isClust[j] ){
325 if (
abs(ohits[centroidIdx]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2)
break;
328 for (
int idx : curN) totPE += ohits[idx]->PE();
329 if (
int(curN.size())<
fMinN)
continue;
332 while (neigh.size() > 0){
333 std::vector<int> curNeigh; curNeigh.push_back(neigh[0]);
334 for (
int j = neigh[0]-1; j > 0; j--){
335 if (
Dist(ohits[neigh[0]],ohits[j]) <
fR0 && !isClust[j] ){
336 curNeigh.push_back(j);
338 if (
abs(ohits[centroidIdx]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2)
341 for (
int j = neigh[0]+1; j <
int(ohits.size()); j++){
342 if (
Dist(ohits[neigh[0]],ohits[j]) <
fR0 && !isClust[j] ){
343 curNeigh.push_back(j);
345 if (
abs(ohits[centroidIdx]->PeakTimeAbs()-ohits[j]->PeakTimeAbs())>2)
349 if (
int(curNeigh.size())>=
fMinN){
350 for (
int cur : curNeigh){
351 if (std::find(curN.begin(),curN.end(),cur)==curN.end()){
353 if (
Dist(ohits[centroidIdx],ohits[cur]) <
fR0)
354 neigh.push_back(cur);
358 neigh.erase(neigh.begin());
361 if (
int(curN.size())<
fMinN)
continue;
365 int channelcentroid = centroid->
OpChannel();
366 double xyzcentroid[3];
367 geo->OpDetGeoFromOpChannel(channelcentroid).GetCenter(xyzcentroid);
368 double yCenter = xyzcentroid[1];
369 double zCenter = xyzcentroid[2];
374 for (
int j = i; j <
int(ohits.size()); j++){
375 if (std::find(curN.begin(),curN.end(),j)!=curN.end())
continue;
376 double r =
YZDist(ohits[j],centroid);
384 for (
int idx : curN) finE += ohits[idx]->PE();
387 std::vector<double>
ys;
388 std::vector<double>
zs;
392 double minY = 1e6;
double maxY = -1e6;
393 double minZ = 1e6;
double maxZ = -1e6;
395 std::vector<double> PEs (
geo->MaxOpChannel() + 1,0.0);
396 std::vector<double> PE2s (
geo->MaxOpChannel() + 1,0.0);
397 double fastToTotal = 0;
398 for (
int hIdx = 0; hIdx <
int(ys.size()); hIdx++){
399 int cIdx = curN[hIdx];
401 minT =
std::min(minT,ohits[cIdx]->PeakTimeAbs());
402 maxT =
std::max(maxT,ohits[cIdx]->PeakTimeAbs());
407 PEs[ohits[cIdx]->OpChannel()] += ohits[cIdx]->PE();
408 PE2s[ohits[cIdx]->OpChannel()] += ohits[cIdx]->PE();
409 fastToTotal += ohits[hIdx]->FastToTotal();
411 double yWidth = maxY-minY;
412 double zWidth = maxZ-minZ;
416 for (
double PE : PEs) tot1 += PE;
417 for (
double PE : PE2s) tot2 += PE;
420 int Frame = ts.OpticalClock().Frame(tCenter - 18.1);
421 if (Frame == 0) Frame = 1;
423 int BeamFrame = ts.OpticalClock().Frame(ts.TriggerTime());
424 bool InBeamFrame =
false;
425 if (!(ts.TriggerTime() < 0)) InBeamFrame = (Frame == BeamFrame);
427 double tWidth = (maxT-minT)/2;
432 oflashes.emplace_back(tCenter,tWidth,tCenter,Frame,
433 PEs,InBeamFrame,OnBeamTime,fastToTotal,
434 yCenter,yWidth,zCenter,zWidth);
435 assoc.emplace_back(curN);
439 for (
int cur : curN) isClust[cur] =
true;
void GetHitYZ(std::vector< art::Ptr< recob::OpHit > >, std::vector< int >, std::vector< double > &, std::vector< double > &) const
int YZCentroid(std::vector< art::Ptr< recob::OpHit > >, std::vector< int >) const
double YZDist(art::Ptr< recob::OpHit >, art::Ptr< recob::OpHit >) const
bool sortOpHitByTime(const art::Ptr< recob::OpHit > &left, const art::Ptr< recob::OpHit > &right)
static int max(int a, int b)
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
double PeakTimeAbs() const
static constexpr double zs
double Dist(art::Ptr< recob::OpHit >, art::Ptr< recob::OpHit >) const
static constexpr double ys
LArSoft geometry interface.