32 #include "Utilities/AssociationUtil.h" 33 #include "Utilities/func/MathUtil.h" 81 Comment(
"Input for physics slices")
86 Comment(
"Input for tracks with which to associate MEs")
93 Comment(
"Max dist from phsyics slice to selected noise hit")
98 Comment(
"Min ADC for a hit to be considered - needed for FD speed ")
104 Name(
"CalibAvailable"),
105 Comment(
"Are there calibration constants available for getting CalE")
111 Name(
"MaxClosestApproach"),
112 Comment(
"ME clust must be at least this close to a hit in phys slice")
116 Name(
"TrkClosestApproach"),
117 Comment(
"Dist cut to be considered TrkME")
151 Name(
"MinDeltaTReTrigger"),
152 Comment(
"See DocDB 12295, apd dead time not well modeled in MC")
158 Name(
"EpsForClustering"),
159 Comment(
"Dist to search out to in the DBSCAN algorithm")
164 Comment(
"For distances in DBSCAN")
169 Comment(
"For distances in DBSCAN")
174 Comment(
"For distances in DBSCAN")
180 Comment(
"Contains hist files for calculating Delta LL ME PID, MID")
185 Name(
"MatchMEtoSlcByDist"),
186 Comment(
"Match MEs to Slices by distance only")
209 void GetPhysNoiseSlcs(
art::Handle<std::vector<rb::Cluster> > slcs,
210 std::vector<rb::Cluster> &noiseSlcs,
211 std::vector<rb::Cluster> &physSlcs);
212 me::VecCHits MECandHits(std::vector<rb::Cluster> noiseSlcs,
213 std::vector<rb::Cluster> physSlcs);
215 int& firstIdx,
int& latsIdx);
222 double HitToHitDistance(
const PtrCHit &hit1,
223 const PtrCHit &hit2);
243 double GetMID(
double deltaT,
double distToSlc,
double calE,
int nCells);
251 double TrkMEDist(
me::SlcME slcME, TVector3 trkEnd);
269 produces<std::vector<me::SlcME> >();
270 produces<art::Assns<me::SlcME, rb::Cluster> >();
271 produces<std::vector<me::TrkME> >();
272 produces<art::Assns<me::TrkME, rb::Track> >();
273 produces<art::Assns<me::TrkME, rb::Cluster> >();
284 TFile* like = TFile::Open(LikeFileName.c_str());
302 std::unique_ptr< std::vector<me::SlcME> >
303 slcMEs(
new std::vector<me::SlcME>);
304 std::unique_ptr< art::Assns<me::SlcME, rb::Cluster> >
307 std::unique_ptr< std::vector<me::TrkME> >
308 trkMEs(
new std::vector<me::TrkME>);
309 std::unique_ptr< art::Assns<me::TrkME, rb::Track> >
311 std::unique_ptr< art::Assns<me::TrkME, rb::Cluster> >
315 std::vector<rb::Cluster> noiseSlcs, physSlcs;
322 for (
unsigned int i = 0;
i < meClusters.size();
i++){
327 const std::vector< art::Ptr<rb::Track> > trks = trkInfo.at(parentSlcIdx);
328 if (
IsTrkME(slcME, trks, parentTrk)){
330 trkMEs->push_back(trkME);
333 *(trkMEToSlcAssn.get()));
336 slcMEs->push_back(slcME);
340 evt.
put(std::move(slcMEs));
341 evt.
put(std::move(slcAssn));
342 evt.
put(std::move(trkMEs));
343 evt.
put(std::move(trkAssn));
344 evt.
put(std::move(trkMEToSlcAssn));
349 std::vector<rb::Cluster> &noiseSlcs,
350 std::vector<rb::Cluster> &physSlcs)
354 for (
unsigned int i = 0;
i < slcs->size();
i++){
357 noiseSlcs.push_back(slice);
359 noiseSlcs.push_back(slice);
361 physSlcs.push_back(slice);
366 std::vector<rb::Cluster> physSlcs){
372 for (
unsigned int nIdx = 0; nIdx < noiseSlcs.size(); nIdx++){
374 for (
unsigned int cIdx = 0; cIdx < noiseSlcs[nIdx].NCell(); cIdx++)
376 noiseHits.push_back(noiseSlcs[nIdx].Cell(cIdx));
378 for (
unsigned int pIdx = 0; pIdx < physSlcs.size(); pIdx++){
379 int fIdx = 0;
int lIdx = 0;
381 for (
int cIdx = fIdx; cIdx <= lIdx; cIdx++){
382 AddMEHit(meHits, noiseHits[cIdx], physSlcs[pIdx]);
392 int& fIdx,
int& lIdx)
397 int nIter =
log2(nHits.size())+1;
398 for (
int i = 1;
i <= nIter;
i++){
399 int fmidpoint = fIdx + (nHits.size()-1)/
pow(2,
i);
400 fIdx = (nHits[fmidpoint]->TNS() > fMETime) ? fIdx : fmidpoint;
401 int lmidpoint = lIdx + nHits.size()/
pow(2,
i);
402 lIdx = (nHits[lmidpoint]->TNS() > lMETime) ? lIdx : lmidpoint;
412 double deltaT = noiseCell->
TNS() - physSlc.
MeanTNS();
418 double curMinDeltaT = ( minDist == 0 && deltaT > 0 ) ?
420 if (deltaT > curMinDeltaT &&
424 for (
int i = 0;
i < (
int)meHits.size();
i++){
425 if (noiseCell->
Plane() == meHits[
i]->Plane() &&
426 noiseCell->
Cell() == meHits[
i]->Cell() &&
427 noiseCell->
ADC() == meHits[
i]->ADC())
430 if (addHit) meHits.push_back(noiseCell);
438 if (meHits.size() == 0)
442 std::vector<bool> isUsed(meHits.size(),
false);
444 for (
unsigned int i = 0;
i < meHits.size();
i++){
445 if (isUsed[
i])
continue;
447 curCluster.push_back(meHits[i]);
449 for (
unsigned int j = i;
j < meHits.size();
j++){
450 if (isUsed[
j])
continue;
453 curCluster.push_back(meHits[j]);
455 for (
unsigned int k = i; k <
j; k++){
456 if (isUsed[k])
continue;
459 curCluster.push_back(meHits[k]);
464 meClusters.push_back(curCluster);
471 double deltaC = (hit1->
View() == hit2->
View()) ?
473 double deltaP = double(hit1->
Plane() - hit2->
Plane());
475 double deltaT = double(hit1->
TNS() - hit2->
TNS());
477 return sqrt(deltaC*deltaC + deltaP*deltaP + deltaT*deltaT);
482 double minDist = 10000;
483 for (
unsigned int i = 0;
i < slc.
NCell();
i++){
485 if (curDist < minDist){
495 if (hit1->
View() != hit2->
View())
return 100000.;
497 double xyz1[3], xyz2[3];
499 unsigned short plane1_idx = hit1->
Plane();
500 unsigned short cell1_idx = hit1->
Cell();
505 unsigned short plane2_idx = hit2->
Plane();
506 unsigned short cell2_idx = hit2->
Cell();
512 return util::pythag(xyz1[0]-xyz2[0], xyz1[1]-xyz2[1], xyz1[2]-xyz2[2]);
523 for (
unsigned int i = 0;
i < michel.size();
i++)
524 slcME.
Add(michel[
i]);
546 for (
unsigned int i = 0;
i < slcME.
NCell();
i++)
573 double assnScore = -1;
575 double allSlcMinDist = 1.0e6;
576 for (
unsigned int i = 0;
i < slcs->size();
i++){
583 double minDist = 10000;
585 for (
unsigned int j = 0;
j < slcME.
NCell();
j++){
594 if ( minDist < allSlcMinDist )
597 allSlcMinDist = minDist;
606 double curScore = hAssn->GetBinContent(hAssn->FindBin(deltaT, minDist));
607 if (curScore > assnScore){
609 assnScore = curScore;
643 int meDeltaT_Dist_Bin = meDeltaT_Dist->FindBin(DeltaT, Dist);
644 int meCalE_NCells_Bin = meCalE_NCells->FindBin(CalE, NCells);
645 double meLike = meDeltaT_Dist->GetBinContent(meDeltaT_Dist_Bin) *
646 meCalE_NCells->GetBinContent(meCalE_NCells_Bin);
651 int hfDeltaT_Dist_Bin = hfDeltaT_Dist->FindBin(DeltaT, Dist);
652 int hfCalE_NCells_Bin = hfCalE_NCells->FindBin(CalE, NCells);
653 double hfLike = hfDeltaT_Dist->GetBinContent(hfDeltaT_Dist_Bin) *
654 hfCalE_NCells->GetBinContent(hfCalE_NCells_Bin);
656 if (meLike == 0 && hfLike == 0)
return -10;
662 return TMath::Log(meLike) - TMath::Log(hfLike);
669 for (
unsigned int i = 0;
i < trks.size();
i++){
670 if (trks[
i]->NXCell() == 0 || trks[
i]->NYCell() == 0)
682 double deltaX = ( (me.
NXCell() == 0) ? 0 : me.
MeanX() - trkEnd.X() );
683 double deltaY = ( (me.
NYCell() == 0) ? 0 : me.
MeanY() - trkEnd.Y() );
684 double deltaZ = me.
MeanZ() - trkEnd.Z();
void SetParentSlc(art::Ptr< rb::Cluster >)
void produce(art::Event &evt)
me::VecCHits MECandHits(std::vector< rb::Cluster > noiseSlcs, std::vector< rb::Cluster > physSlcs)
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
void AddMEHit(me::VecCHits &meHits, me::PtrCHit noiseCell, rb::Cluster physSlc)
void GetCenter(double *xyz, double localz=0.0) const
art::ServiceHandle< geo::Geometry > fgeom
art::Ptr< rb::CellHit > PtrCHit
void SetDistToSlc(double)
unsigned short Plane() const
art::Ptr< rb::Cluster > ParentSlc() const
double GetMID(double deltaT, double distToSlc, double calE, int nCells)
const CellGeo * Cell(int icell) const
std::vector< PtrCHit > VecCHits
A collection of associated CellHits.
double MeanZ(rb::AveragingScheme scheme=kDefaultScheme) const
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
std::string find_file(std::string const &filename) const
double TrkMEDist(me::SlcME slcME, TVector3 trkEnd)
Atom< double > planeScale
art::Ptr< rb::Track > ParentTrk() const
double MEDistanceMetric(me::PtrCHit hit1, me::PtrCHit hit2)
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
double MeanX(rb::AveragingScheme scheme=kDefaultScheme) const
ProductID put(std::unique_ptr< PROD > &&product)
Atom< bool > calibAvailable
art::Ptr< rb::Track > parentTrk
unsigned short Cell() const
Atom< unsigned int > maxNCells
double MeanY(rb::AveragingScheme scheme=kDefaultScheme) const
Atom< unsigned int > minNCells
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Geometry information for a single readout plane.
void BestSlcMatch(me::SlcME michel, me::MEVars &vars, art::Handle< std::vector< rb::Cluster > > slcs, int &parentSlcIdx)
Atom< std::string > trackInput
bool IsTrkME(me::SlcME michel, const std::vector< art::Ptr< rb::Track > > trks, art::Ptr< rb::Track > &parentTrk)
bool PassesMEPresel(me::SlcME michel)
std::vector< me::VecCHits > VecVecCHits
Atom< std::string > llFileName
Atom< double > minDeltaTReTrigger
Atom< double > trkClosestApproach
unsigned int NYCell() const
Number of cells in the y-view.
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Atom< double > maxClosestApproach
int16_t ADC(uint32_t i) const
void SetParentTrk(art::Ptr< rb::Track >)
MEFinder(const Parameters ¶ms)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
double pythag(double x, double y)
2D Euclidean distance
unsigned int NXCell() const
Number of cells in the x-view.
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
me::SlcME MakeSlcME(me::VecCHits michel, art::Handle< std::vector< rb::Cluster > > slcs, int &parentSlcIdx)
me::TrkME MakeTrkME(me::SlcME michel, art::Handle< std::vector< rb::Cluster > > slcs, art::Ptr< rb::Track > parentTrack)
Atom< std::string > slicerInput
void GetPhysNoiseSlcs(art::Handle< std::vector< rb::Cluster > > slcs, std::vector< rb::Cluster > &noiseSlcs, std::vector< rb::Cluster > &physSlcs)
Atom< bool > matchMEtoSlcByDist
bool IsNoise() const
Is the noise flag set?
Encapsulate the cell geometry.
art::Ptr< rb::Cluster > parentSlc
fvar< T > log2(const fvar< T > &x)
void FindNoiseHitRange(me::VecCHits hits, double physTNS, int &firstIdx, int &latsIdx)
double MinHitSlcDist(me::PtrCHit hit, rb::Cluster slc)
double HitToHitDistance(const PtrCHit &hit1, const PtrCHit &hit2)
Encapsulate the geometry of one entire detector (near, far, ndos)
me::VecVecCHits ClusterMEHits(VecCHits)