Inheritance diagram for opdet::OpFlashFinderDualPhase:

Public Member Functions
	OpFlashFinderDualPhase (const fhicl::ParameterSet &)

virtual	~OpFlashFinderDualPhase ()

void	beginJob ()

void	endJob ()

void	reconfigure (fhicl::ParameterSet const &pset)

void	produce (art::Event &)

void	RunFlashFinder (std::vector< recob::OpHit > const &HitVector, std::vector< recob::OpFlash > &FlashVector, std::vector< std::vector< int > > &AssocList, geo::GeometryCore const &geom, detinfo::DetectorClocksData const &ts, float const &TrigCoinc)

void	AddHitContribution (recob::OpHit const &currentHit, double &MaxTime, double &MinTime, double &AveTime, double &FastToTotal, double &AveAbsTime, double &TotalPE, std::vector< double > &PEs)

void	GetHitGeometryInfo (recob::OpHit const &currentHit, geo::GeometryCore const &geom, std::vector< double > &sumw, std::vector< double > &sumw2, double &sumy, double &sumy2, double &sumz, double &sumz2)

double	CalculateWidth (double const &sum, double const &sum_squared, double const &weights_sum)

std::vector< int >	getNeighbors (std::vector< recob::OpHit > const &HitVector, int hitnumber, std::vector< bool > &processed, float initimecluster, std::vector< int > &sorted, geo::GeometryCore const &geom)

void	AssignHitsToFlash (std::vector< recob::OpHit > const &HitVector, std::vector< std::vector< int > > &HitsPerFlash, geo::GeometryCore const &geom)

void	ConstructFlash (std::vector< int > const &HitsPerFlashVec, std::vector< recob::OpHit > const &HitVector, std::vector< recob::OpFlash > &FlashVector, geo::GeometryCore const &geom, detinfo::DetectorClocksData const &ts, float const &TrigCoinc)

Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDProducer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept

	Producer (fhicl::ParameterSet const &)

	Producer (Producer const &)=delete

	Producer (Producer &&)=delete

Producer &	operator= (Producer const &)=delete

Producer &	operator= (Producer &&)=delete

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

Public Member Functions inherited from art::Modifier
	~Modifier () noexcept

	Modifier ()

	Modifier (Modifier const &)=delete

	Modifier (Modifier &&)=delete

Modifier &	operator= (Modifier const &)=delete

Modifier &	operator= (Modifier &&)=delete

Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept

	ModuleBase ()

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Attributes
std::string	fInputModule

Double_t	fMaximumDistance

Double_t	fMaximumTimeDistance

Double_t	fMaximumTimeWindow

Double_t	fTrigCoinc

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >

Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)

Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 39 of file OpFlashFinderDualPhase_module.cc.

Constructor & Destructor Documentation

opdet::OpFlashFinderDualPhase::OpFlashFinderDualPhase ( const fhicl::ParameterSet & pset )

explicit

Definition at line 115 of file OpFlashFinderDualPhase_module.cc.

     : EDProducer{pset}
   {
 
     reconfigure(pset);
 
     produces< std::vector< recob::OpFlash > >();
     produces< art::Assns< recob::OpFlash, recob::OpHit > >();
 
   }

opdet::OpFlashFinderDualPhase::~OpFlashFinderDualPhase ( )

virtual

Definition at line 142 of file OpFlashFinderDualPhase_module.cc.

143 {

144 }

Member Function Documentation

void opdet::OpFlashFinderDualPhase::AddHitContribution	(	recob::OpHit const &	currentHit,
		double &	MaxTime,
		double &	MinTime,
		double &	AveTime,
		double &	FastToTotal,
		double &	AveAbsTime,
		double &	TotalPE,
		std::vector< double > &	PEs
	)

Definition at line 331 of file OpFlashFinderDualPhase_module.cc.

                                                     {
 
     double PEThisHit   = currentHit.PE();
     double TimeThisHit = currentHit.PeakTime();
     if (TimeThisHit > MaxTime) MaxTime = TimeThisHit;
     if (TimeThisHit < MinTime) MinTime = TimeThisHit;
     
     // These quantities for the flash are defined 
     // as the weighted averages over the hits
     AveTime     += PEThisHit*TimeThisHit;
     FastToTotal += PEThisHit*currentHit.FastToTotal();
     AveAbsTime  += PEThisHit*currentHit.PeakTimeAbs();
     
     // These are totals
     TotalPE     += PEThisHit;
     PEs.at(static_cast< unsigned int >(currentHit.OpChannel())) += PEThisHit;
 
   }

void opdet::OpFlashFinderDualPhase::AssignHitsToFlash	(	std::vector< recob::OpHit > const &	HitVector,
		std::vector< std::vector< int > > &	HitsPerFlash,
		geo::GeometryCore const &	geom
	)

Definition at line 438 of file OpFlashFinderDualPhase_module.cc.

  {
 
    int ntothits=HitVector.size();
    std::size_t n(0);
    std::vector<int> sorted;
    sorted.resize(ntothits);
    std::generate(std::begin(sorted), std::end(sorted), [&]{ return n++; });
    std::sort(  std::begin(sorted), std::end(sorted), [&](int i1, int i2) { return HitVector[i1].PeakTime() < HitVector[i2].PeakTime(); }); 
 
    HitsPerFlash.clear();
 
    std::vector<bool> processed;
    processed.resize(HitVector.size(),false);
 
    for (unsigned int h=0;h<HitVector.size();h++)
    { 
      if(!processed[h])
      {
        processed[h]=true;
        std::vector<int> N; N.clear();
        std::vector<int> nb =getNeighbors(HitVector,h, processed, HitVector[sorted[h]].PeakTime(),sorted,geom); 
        N.insert( N.end(), nb.begin(), nb.end() );
 
        for (unsigned int i=0;i<N.size();i++)
        {
            std::vector<int> nb2 =getNeighbors(HitVector,N[i], processed, HitVector[sorted[h]].PeakTime(),sorted,geom); 
            N.insert( N.end(), nb2.begin(), nb2.end() );
        }
            N.push_back(sorted[h]);
            HitsPerFlash.push_back(N);
      } 
    }
     
  } // End AssignHitsToFlash

void opdet::OpFlashFinderDualPhase::beginJob ( )

virtual

Reimplemented from art::EDProducer.

Definition at line 147 of file OpFlashFinderDualPhase_module.cc.

148 {

149 }

double opdet::OpFlashFinderDualPhase::CalculateWidth	(	double const &	sum,
		double const &	sum_squared,
		double const &	weights_sum
	)

Definition at line 390 of file OpFlashFinderDualPhase_module.cc.

                                                    {
 
     if (sum_squared*weights_sum - sum*sum < 0) return 0;
     else return std::sqrt(sum_squared*weights_sum - sum*sum)/weights_sum;
 
   }

void opdet::OpFlashFinderDualPhase::ConstructFlash	(	std::vector< int > const &	HitsPerFlashVec,
		std::vector< recob::OpHit > const &	HitVector,
		std::vector< recob::OpFlash > &	FlashVector,
		geo::GeometryCore const &	geom,
		detinfo::DetectorClocksData const &	ts,
		float const &	TrigCoinc
	)

Definition at line 240 of file OpFlashFinderDualPhase_module.cc.

                                                                     {
 
     double MaxTime = -1e9;
     double MinTime =  1e9;
     
     std::vector< double > PEs(geom.MaxOpChannel() + 1, 0.0);
     unsigned int Nplanes = geom.Nplanes();
     std::vector< double > sumw (Nplanes, 0.0);
     std::vector< double > sumw2(Nplanes, 0.0);
     
     double TotalPE     = 0;
     double AveTime     = 0;
     double AveAbsTime  = 0;
     double FastToTotal = 0;
     double sumy        = 0;
     double sumz        = 0;
     double sumy2       = 0;
     double sumz2       = 0;
 
     for (auto const& HitID : HitsPerFlashVec) {
       AddHitContribution(HitVector.at(HitID),
                          MaxTime,
                          MinTime,
                          AveTime,
                          FastToTotal,
                          AveAbsTime,
                          TotalPE,
                          PEs);
       GetHitGeometryInfo(HitVector.at(HitID),
                          geom,
                          sumw,
                          sumw2,
                          sumy, 
                          sumy2,
                          sumz, 
                          sumz2);
     }
 
     AveTime     /= TotalPE;
     AveAbsTime  /= TotalPE;
     FastToTotal /= TotalPE;
     
     double meany = sumy/TotalPE;
     double meanz = sumz/TotalPE;
     
     double widthy = CalculateWidth(sumy, sumy2, TotalPE);
     double widthz = CalculateWidth(sumz, sumz2, TotalPE);
     
     std::vector< double > WireCenters(Nplanes, 0.0);
     std::vector< double > WireWidths(Nplanes, 0.0);
     
     for (size_t p = 0; p != Nplanes; ++p) {
       WireCenters.at(p) = sumw.at(p)/TotalPE;
       WireWidths.at(p)  = CalculateWidth(sumw.at(p), sumw2.at(p), TotalPE);
     }
 
     // Emprical corrections to get the Frame right.
     // Eventual solution - remove frames
     int Frame = ts.OpticalClock().Frame(AveAbsTime - 18.1);
     if (Frame == 0) Frame = 1;
     
     int BeamFrame = ts.OpticalClock().Frame(ts.TriggerTime());
     bool InBeamFrame = false;
     if (!(ts.TriggerTime() < 0)) InBeamFrame = (Frame == BeamFrame);
 
     double TimeWidth = (MaxTime - MinTime)/2.0;
 
     int OnBeamTime = 0; 
     if (InBeamFrame && (std::abs(AveTime) < TrigCoinc)) OnBeamTime = 1;
     
     FlashVector.emplace_back(AveTime,
                              TimeWidth,
                              AveAbsTime,
                              Frame,
                              PEs, 
                              InBeamFrame,
                              OnBeamTime,
                              FastToTotal,
                              meany, 
                              widthy, 
                              meanz, 
                              widthz, 
                              WireCenters, 
                              WireWidths);
 
   }

void opdet::OpFlashFinderDualPhase::endJob ( )

virtual

Reimplemented from art::EDProducer.

Definition at line 152 of file OpFlashFinderDualPhase_module.cc.

153 {

154 }

void opdet::OpFlashFinderDualPhase::GetHitGeometryInfo	(	recob::OpHit const &	currentHit,
		geo::GeometryCore const &	geom,
		std::vector< double > &	sumw,
		std::vector< double > &	sumw2,
		double &	sumy,
		double &	sumy2,
		double &	sumz,
		double &	sumz2
	)

Definition at line 358 of file OpFlashFinderDualPhase_module.cc.

                                                           {
 
     double xyz[3];
     geom.OpDetGeoFromOpChannel(currentHit.OpChannel()).GetCenter(xyz);
     double PEThisHit = currentHit.PE();
     
     geo::TPCID tpc = geom.FindTPCAtPosition(xyz);
     // if the point does not fall into any TPC,
     // it does not contribute to the average wire position
     if (tpc.isValid) {
       for (size_t p = 0; p != geom.Nplanes(); ++p) {
         geo::PlaneID const planeID(tpc, p);
         unsigned int w = geom.NearestWire(xyz, planeID);
         sumw.at(p)  += PEThisHit*w;
         sumw2.at(p) += PEThisHit*w*w;
       }
     } // if we found the TPC
     sumy  += PEThisHit*xyz[1]; 
     sumy2 += PEThisHit*xyz[1]*xyz[1];
     sumz  += PEThisHit*xyz[2]; 
     sumz2 += PEThisHit*xyz[2]*xyz[2];
     
   }

std::vector< int > opdet::OpFlashFinderDualPhase::getNeighbors	(	std::vector< recob::OpHit > const &	HitVector,
		int	hitnumber,
		std::vector< bool > &	processed,
		float	initimecluster,
		std::vector< int > &	sorted,
		geo::GeometryCore const &	geom
	)

Definition at line 400 of file OpFlashFinderDualPhase_module.cc.

  {
 
   int itTmin=0;
   int itTmax=HitVector.size();
   for (int h=hitnumber;h>0;h--) if(HitVector[sorted[h]].PeakTime() < HitVector[sorted[hitnumber]].PeakTime()-fMaximumTimeDistance)
   {
     itTmin=h; break;
   }
   for (unsigned int h=hitnumber;h<HitVector.size();h++) if(HitVector[sorted[h]].PeakTime() > HitVector[sorted[hitnumber]].PeakTime()+fMaximumTimeDistance)
   {
     itTmax=h; break;
   }
 
    std::vector<int> neighbors;
    float dx, dy, dz, dt, dtmax;
    double xyz_hitnumber[3];
    geom.OpDetGeoFromOpChannel(HitVector[sorted[hitnumber]].OpChannel()).GetCenter(xyz_hitnumber);
 
    for (int h=itTmin;h<itTmax;h++)
    {
      if(!processed[h])
      {
        double xyz_h[3];
        geom.OpDetGeoFromOpChannel(HitVector[sorted[h]].OpChannel()).GetCenter(xyz_h);
        dx  = xyz_h[0] - xyz_hitnumber[0];
        dy  = xyz_h[1] - xyz_hitnumber[1];
        dz  = xyz_h[2] - xyz_hitnumber[2];
        dt  = TMath::Abs(  HitVector[sorted[hitnumber]].PeakTime() - HitVector[sorted[h]].PeakTime() );
        dtmax  = TMath::Abs(initimecluster - HitVector[sorted[h]].PeakTime() );
        float distance = dx*dx + dy*dy + dz*dz;
 //       std::cout << distance << " " << dt << " " << dtmax << std::endl; lets_pause();
        if(distance<fMaximumDistance*fMaximumDistance && dt<fMaximumTimeDistance && dtmax<fMaximumTimeWindow){ neighbors.push_back(h);processed[h]=true;}
      }
    }
    return neighbors;
  }

void opdet::OpFlashFinderDualPhase::produce ( art::Event & evt )

virtual

Implements art::EDProducer.

Definition at line 157 of file OpFlashFinderDualPhase_module.cc.

   {
 
     // These are the storage pointers we will put in the event
     std::unique_ptr< std::vector< recob::OpFlash > > 
                       flashPtr(new std::vector< recob::OpFlash >);
     std::unique_ptr< art::Assns< recob::OpFlash, recob::OpHit > >  
                       assnPtr(new art::Assns< recob::OpFlash, recob::OpHit >);
 
     // This will keep track of what flashes will assoc to what ophits
     // at the end of processing
     std::vector< std::vector< int > > assocList;
 
     auto const& geometry(*lar::providerFrom< geo::Geometry >());
 
     auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
 
     // Get OpHits from the event
     auto opHitHandle = evt.getHandle< std::vector< recob::OpHit > >(fInputModule);
 
     RunFlashFinder(*opHitHandle,
                    *flashPtr,
                    assocList,
                    geometry,
                    clockData,
                    fTrigCoinc);
 
 
    
 
     // Make the associations which we noted we need
     for (size_t i = 0; i != assocList.size(); ++i)
     {
       art::PtrVector< recob::OpHit > opHitPtrVector;
       for (int const& hitIndex : assocList.at(i))
       {
         art::Ptr< recob::OpHit > opHitPtr(opHitHandle, hitIndex);
         opHitPtrVector.push_back(opHitPtr);
       }
 
       util::CreateAssn(*this, evt, *flashPtr, opHitPtrVector, 
                                         *(assnPtr.get()), i);
     }
     
     // Store results into the event
     evt.put(std::move(flashPtr));
     evt.put(std::move(assnPtr));
     
   }

void opdet::OpFlashFinderDualPhase::reconfigure ( fhicl::ParameterSet const & pset )

Definition at line 127 of file OpFlashFinderDualPhase_module.cc.

   {
 
     // Indicate that the Input Module comes from .fcl
     fInputModule = pset.get< std::string >("InputModule");
       
     fTrigCoinc      = pset.get< double >("TrigCoinc");
     fMaximumDistance     = pset.get< double >("MaximumDistance");
     fMaximumTimeDistance = pset.get< double >("MaximumTimeDistance");
     fMaximumTimeWindow   = pset.get< double >("MaximumTimeWindow");
 
   }

void opdet::OpFlashFinderDualPhase::RunFlashFinder	(	std::vector< recob::OpHit > const &	HitVector,
		std::vector< recob::OpFlash > &	FlashVector,
		std::vector< std::vector< int > > &	AssocList,
		geo::GeometryCore const &	geom,
		detinfo::DetectorClocksData const &	ts,
		float const &	TrigCoinc
	)

Definition at line 207 of file OpFlashFinderDualPhase_module.cc.

                                                                     {
 
     // Now start to create flashes.
     // First, need vector to keep track of which hits belong to which flashes
     std::vector< std::vector< int > > HitsPerFlash;
     
     AssignHitsToFlash(HitVector,
                       HitsPerFlash,
                       geom);
     
     // Now we have all our hits assigned to a flash. 
     // Make the recob::OpFlash objects
     for (auto const& HitsPerFlashVec : HitsPerFlash)
       ConstructFlash(HitsPerFlashVec,
                      HitVector,
                      FlashVector,
                      geom,
                      ts,
                      TrigCoinc);
 
     // Finally, write the association list.
     // back_inserter tacks the result onto the end of AssocList
     for (auto& HitIndicesThisFlash : HitsPerFlash)
       AssocList.push_back(HitIndicesThisFlash);
     
   } // End RunFlashFinder