vertex::PrimaryVertexFinder Class Reference

Inheritance diagram for vertex::PrimaryVertexFinder:

Public Member Functions
	PrimaryVertexFinder (fhicl::ParameterSet const &pset)
void	beginJob ()
void	produce (art::Event &evt)
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 Member Functions
double	StartPointSeperation (recob::SpacePoint sp1, recob::SpacePoint sp2)
bool	IsInVertexCollection (int a, std::vector< std::vector< int > > vertex_collection)
int	IndexInVertexCollection (int a, int b, std::vector< std::vector< int > > vertex_collection)
bool	IsInNewVertex (int a, std::vector< int > newvertex)
double	gammavalue (TVector3 startpoint1, TVector3 startpoint2, TVector3 dircos1, TVector3 dircos2)
double	alphavalue (double gamma, TVector3 startpoint1, TVector3 startpoint2, TVector3 dircos1, TVector3 dircos2)
double	MinDist (double alpha, double gamma, TVector3 startpoint1, TVector3 startpoint2, TVector3 dircos1, TVector3 dircos2)
TVector3	PointOnExtendedTrack (double alphagamma, TVector3 startpoint, TVector3 dircos)

Private Attributes
std::string	fTrackModuleLabel
double	fVertexWindow
TH2F *	fNoTracks
TH1F *	fLength_1stTrack
TH1F *	fLength_2ndTrack
TH1F *	fLength_3rdTrack
TH1F *	fLength_4thTrack
TH1F *	fLength_5thTrack

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 41 of file PrimaryVertexFinder_module.cc.

Constructor & Destructor Documentation

vertex::PrimaryVertexFinder::PrimaryVertexFinder ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 84 of file PrimaryVertexFinder_module.cc.

    : EDProducer{pset}
  {
    fTrackModuleLabel  = pset.get< std::string >("TrackModuleLabel");
    fVertexWindow      = pset.get<double     >  ("VertexWindow");

    produces< std::vector<recob::Vertex> >();
    produces< art::Assns<recob::Vertex, recob::Hit> >();
    produces< art::Assns<recob::Vertex, recob::Track> >();
    produces< art::Assns<recob::Vertex, recob::Shower> >();
  }

Member Function Documentation

double vertex::PrimaryVertexFinder::alphavalue ( double  gamma, TVector3  startpoint1, TVector3  startpoint2, TVector3  dircos1, TVector3  dircos2  )

private

Definition at line 444 of file PrimaryVertexFinder_module.cc.

{
  double alpha = (gamma*(dircos1*dircos2)) + (startpoint2*dircos1) - (startpoint1*dircos1);

  return alpha;
}

void vertex::PrimaryVertexFinder::beginJob ( )

virtual

Reimplemented from art::EDProducer.

Definition at line 97 of file PrimaryVertexFinder_module.cc.

                                    {
    // get access to the TFile service
    art::ServiceHandle<art::TFileService const> tfs;

    //    fNoVertices= tfs->make<TH2F>("fNoVertices", ";Event No; No of vertices", 100,0, 100, 30, 0, 30);
     fNoTracks= tfs->make<TH2F>("fNoTracks", ";Event No; No of Tracks", 10,0, 10, 10, 0, 10);
     fLength_1stTrack = tfs->make<TH1F>("fLength_Track1", "Muon Track Length", 100,0,100);
     fLength_2ndTrack = tfs->make<TH1F>("fLength_Track2", "2nd Track Length", 100,0,100);
     fLength_3rdTrack = tfs->make<TH1F>("fLength_Track3", "3rd Track Length", 100,0,100);
     fLength_4thTrack = tfs->make<TH1F>("fLength_Track4", "4th Track Length", 100,0,100);
     fLength_5thTrack = tfs->make<TH1F>("fLength_Track5", "5th Track Length", 100,0,100);
  }

double vertex::PrimaryVertexFinder::gammavalue ( TVector3  startpoint1, TVector3  startpoint2, TVector3  dircos1, TVector3  dircos2  )

private

Definition at line 437 of file PrimaryVertexFinder_module.cc.

{
  double gamma = ((startpoint1*dircos2)-(startpoint2*dircos2)+((dircos1*dircos2)*(startpoint2*dircos1))-((dircos1*dircos2)*(startpoint1*dircos1)))/(1-((dircos1*dircos2)*(dircos1*dircos2)));

  return gamma;
}

int vertex::PrimaryVertexFinder::IndexInVertexCollection ( int  a, int  b, std::vector< std::vector< int > >  vertex_collection  )

private

Definition at line 410 of file PrimaryVertexFinder_module.cc.

{
  int index = -1;
  for(unsigned int i = 0; i < vertex_collection.size() ; i++){
    for(std::vector<int>::iterator itr = vertex_collection[i].begin(); itr < vertex_collection[i].end(); ++itr){
      if (a == *itr || b == *itr)
        index = i;
    }
  }
  return index;
}

bool vertex::PrimaryVertexFinder::IsInNewVertex ( int  a, std::vector< int >  newvertex  )

private

Definition at line 422 of file PrimaryVertexFinder_module.cc.

{
  int flag = 0;
  for(unsigned int i = 0; i < newvertex.size() ; i++){
    if (a == newvertex[i]){
      flag = 1;
      break;
    }
  }

  if(flag==1)
    return true;
  return false;
}

bool vertex::PrimaryVertexFinder::IsInVertexCollection ( int  a, std::vector< std::vector< int > >  vertex_collection  )

private

Definition at line 393 of file PrimaryVertexFinder_module.cc.

{
  int flag = 0;

  for(unsigned int i = 0; i < vertex_collection.size() ; i++){
    for(std::vector<int>::iterator itr = vertex_collection[i].begin(); itr < vertex_collection[i].end(); ++itr){
      if (a == *itr){
        flag = 1;
        break;
      }
    }
  }
  if(flag==1)
    return true;
  return false;
}

double vertex::PrimaryVertexFinder::MinDist ( double  alpha, double  gamma, TVector3  startpoint1, TVector3  startpoint2, TVector3  dircos1, TVector3  dircos2  )

private

Definition at line 451 of file PrimaryVertexFinder_module.cc.

{
  TVector3 mindis_vector = startpoint1 - startpoint2 + alpha*dircos1 - gamma*dircos2;
  double mindis = mindis_vector.Mag();
  return mindis;
}

TVector3 vertex::PrimaryVertexFinder::PointOnExtendedTrack ( double  alphagamma, TVector3  startpoint, TVector3  dircos  )

private

Definition at line 458 of file PrimaryVertexFinder_module.cc.

{
  TVector3 PointOnExtendedTrack = startpoint + (alphagamma * dircos);
  return PointOnExtendedTrack;
}

void vertex::PrimaryVertexFinder::produce ( art::Event &  evt )

virtual

Todo:

replace with the real query when this module is updated to look for showers too

Implements art::EDProducer.

Definition at line 111 of file PrimaryVertexFinder_module.cc.

  {

    mf::LogInfo("PrimaryVertexFinder") << "------------------------------------------------------------------------------";

    //   std::cout << "run    : " << evt.Header().Run() << std::endl;
    //   std::cout << "subrun : " << evt.Header().Subrun() << std::endl;
    //std::cout << "event  : " << evt.Header().Event() << std::endl;

    mf::LogInfo("PrimaryVertexFinder") << "event  : " << evt.id().event();


    art::ServiceHandle<geo::Geometry const> geom;

    //mf::LogInfo("PrimaryVertexFinder") << "I am in Primary vertex finder " << std::endl;

    art::Handle< std::vector<recob::Track> > trackListHandle;
    evt.getByLabel(fTrackModuleLabel,trackListHandle);

    //Point to a collection of vertices to output.
    std::unique_ptr< std::vector<recob::Vertex> > vcol(new std::vector<recob::Vertex>);
    std::unique_ptr< art::Assns<recob::Vertex, recob::Hit> > vhassn(new art::Assns<recob::Vertex, recob::Hit>);
    std::unique_ptr< art::Assns<recob::Vertex, recob::Track> > vtassn(new art::Assns<recob::Vertex, recob::Track>);
    std::unique_ptr< art::Assns<recob::Vertex, recob::Shower> > vsassn(new art::Assns<recob::Vertex, recob::Shower>);


    std::vector<recob::Track> const& trkIn = *trackListHandle;

    mf::LogInfo("PrimaryVertexFinder") << "number of tracks in this event = " << trkIn.size();
    fNoTracks->Fill(evt.id().event(),trkIn.size());

    std::vector<recob::SpacePoint> startpoints_vec; // first space point of each track

    std::vector <TVector3> startvec;
    TVector3 startXYZ;

    std::vector <TVector3> endvec;
    TVector3 endXYZ;

    std::vector <TVector3> dircosvec;
    TVector3 dircosXYZ;

    std::vector< std::pair<art::Ptr<recob::Track>, double> > trackpair;

    art::FindMany<recob::SpacePoint> TrackSpacePoints
      (trackListHandle, evt, fTrackModuleLabel);

    for(unsigned int i = 0; i<trkIn.size(); ++i){
      recob::Track::Point_t start, end;
      std::tie(start, end) = trkIn[i].Extent();
      startXYZ.SetXYZ(start.X(),start.Y(),start.Z());
      endXYZ.SetXYZ(end.X(),end.Y(),end.Z());


      double length = (endXYZ-startXYZ).Mag();// (endvec[i]-startvec[i]).Mag();
      //mf::LogInfo("PrimaryVertexFinder") << "Track length calculated = " << length << std::endl;
      trackpair.push_back(std::pair<art::Ptr<recob::Track>,double>({ trackListHandle, i },length));
    }

    for(size_t i = 0; i<trackpair.size(); ++i){
      mf::LogInfo("PrimaryVertexFinder") << "track id is  = " << (trackpair[i].first)->ID()
                                         << " track length = " << (trackpair[i].second);
    }

    std::sort(trackpair.rbegin(), trackpair.rend(), sort_pred2);

    mf::LogInfo("PrimaryVertexFinder") << "AFTER SORTING ";
    for(size_t i = 0; i < trackpair.size(); ++i){
      mf::LogInfo("PrimaryVertexFinder") << "track id is  = " << (trackpair[i].first)->ID()
                                         << " track length = " << (trackpair[i].second);
    }

    if(trackpair.size()>0)
    fLength_1stTrack->Fill(trackpair[0].second);

    if(trackpair.size()>1)
    fLength_2ndTrack->Fill(trackpair[1].second);

    if(trackpair.size()>2)
    fLength_3rdTrack->Fill(trackpair[2].second);

    if(trackpair.size()>3)
    fLength_4thTrack->Fill(trackpair[3].second);

    if(trackpair.size()>4)
    fLength_5thTrack->Fill(trackpair[4].second);

    for(size_t j = 0; j < trackpair.size(); ++j) { //loop over tracks
      art::Ptr<recob::Track> const& track = trackpair[j].first;

      // the index of this track in the query is the same as its position in
      // the data product:
      std::vector<recob::SpacePoint const*> const& spacepoints
        = TrackSpacePoints.at(track.key());

      startXYZ  = trackpair[j].first->Vertex<TVector3>();
      endXYZ    = trackpair[j].first->End<TVector3>();
      dircosXYZ = trackpair[j].first->VertexDirection<TVector3>();

      startvec.push_back(startXYZ);
      endvec.push_back(endXYZ);
      dircosvec.push_back(dircosXYZ);

      mf::LogInfo("PrimaryVertexFinder") << "PrimaryVertexFinder got "<< spacepoints.size()
                                         <<" 3D spacepoint(s) from Track3Dreco.cxx";

      // save the first SpacePoint of each Track... from now the SpacePoint ID represents the Track ID!!
      startpoints_vec.emplace_back(
        spacepoints[0]->XYZ(), spacepoints[0]->ErrXYZ(),
        spacepoints[0]->Chisq(), startpoints_vec.size()
        );

    }// loop over tracks

    for(size_t i = 0; i < startvec.size();  ++i){ //trackpair.size()
      mf::LogInfo("PrimaryVertexFinder") << "Tvector3 start point SORTED = ";
      startvec[i].Print();
    }
    for(size_t i = 0; i < dircosvec.size(); ++i){ //trackpair.size()
      mf::LogInfo("PrimaryVertexFinder") << "Tvector3 dir cos SORTED = ";
      dircosvec[i].Print();
    }

    std::vector<std::vector<int> > vertex_collection_int;
    std::vector <std::vector <TVector3> > vertexcand_vec;

    for (unsigned int i=0; i<trackpair.size(); ++i){
      for (unsigned int j=i+1; j<trackpair.size(); ++j){
        mf::LogInfo("PrimaryVertexFinder") << "distance between " << i << " and " << j
                                           << " = "
                                           << StartPointSeperation(startpoints_vec[i], startpoints_vec[j]);
        double GAMMA = gammavalue(startvec[i], startvec[j], dircosvec[i], dircosvec[j]);
        double ALPHA = alphavalue(GAMMA, startvec[i], startvec[j], dircosvec[i], dircosvec[j]);
        double MINDIST = MinDist(ALPHA, GAMMA, startvec[i], startvec[j], dircosvec[i], dircosvec[j]);
        mf::LogInfo("PrimaryVertexFinder") << "alpha = " << ALPHA << " gamma = "
                                           << GAMMA << " MINIMUM DISTANCE = " << MINDIST;

        TVector3 TRACK1POINT = PointOnExtendedTrack(ALPHA, startvec[i], dircosvec[i]);
        TVector3 TRACK2POINT = PointOnExtendedTrack(GAMMA, startvec[j], dircosvec[j]);

        mf::LogInfo("PrimaryVertexFinder") << "POINTS ON THE TRACKS ARE:: ";
        TRACK1POINT.Print();
        TRACK2POINT.Print();

        //if(StartPointSeperation(startpoints_vec[i], startpoints_vec[j])<fVertexWindow){ ///// correct this
        //if(MINDIST<2 && trackpair[i].second >30 && trackpair[j].second >30){
        if(MINDIST < fVertexWindow && ((TRACK1POINT-startvec[i]).Mag()) < fVertexWindow){

          if((!IsInVertexCollection(i, vertex_collection_int)) && (!IsInVertexCollection(j, vertex_collection_int))){
            std::vector<int> newvertex_int;
            std::vector <TVector3> vertexcand;
            newvertex_int.push_back(i);
            newvertex_int.push_back(j);
            vertex_collection_int.push_back(newvertex_int);
            //newvertex.clear();
            vertexcand.push_back(TRACK1POINT);
            vertexcand.push_back(TRACK2POINT);
            vertexcand_vec.push_back(vertexcand);
          }
          else{
            int index = IndexInVertexCollection(i, j, vertex_collection_int);
            //mf::LogInfo("PrimaryVertexFinder") << "index where a new vertex will be added = " << index << std::endl;
            if(!IsInNewVertex(i, vertex_collection_int[index])){
              vertex_collection_int[index].push_back(i);
              vertexcand_vec[index].push_back(TRACK1POINT); //need to fix for delta rays
            }
            if(!IsInNewVertex(j, vertex_collection_int[index])){
              vertex_collection_int[index].push_back(j);
              vertexcand_vec[index].push_back(TRACK2POINT); //need to fix for delta rays
            }
          }
        }// end else
      }
    }


    //now add the unmatched track IDs to the collection
    for(size_t i = 0; i < trackpair.size(); ++i){
      if(!IsInVertexCollection(i, vertex_collection_int)){
        //if(trackpair[i].second>30){
        std::vector<int> temp;
        std::vector <TVector3> temp1;
        temp.push_back(i);
        temp1.push_back(startvec[i]);
        vertex_collection_int.push_back(temp);
        vertexcand_vec.push_back(temp1);
        //}
      }
    }

    // indices (in their data products) of tracks and showers connected to the vertex
    std::vector<size_t> vTrackIndices, vShowerIndices;

    // find the hits of all the tracks
    art::FindManyP<recob::Hit> TrackHits(trackListHandle, evt, fTrackModuleLabel);

    // find the hits of all the showers
  //  art::FindManyP<recob::Hit> ShowerHits(showerListHandle, evt, fShowerModuleLabel);
    ///\todo replace with the real query when this module is updated to look for showers too
    art::FindManyP<recob::Hit> ShowerHits(std::vector<art::Ptr<recob::Shower>>(), evt, fTrackModuleLabel);

    for(size_t i = 0; i < vertex_collection_int.size(); ++i){
      double x = 0.;
      double y = 0.;
      double z = 0.;
      int elemsize = 0.;
      for(std::vector<int>::iterator itr = vertex_collection_int[i].begin(); itr < vertex_collection_int[i].end(); ++itr){
        mf::LogInfo("PrimaryVertexFinder") << "vector elements at index " << i << " are " << *itr
                                           << "\ntrack original ID = " << (trackpair[*itr].first)->ID();
        // save the index in the data product of this track
        vTrackIndices.push_back(trackpair[*itr].first.key());
      }
      mf::LogInfo("PrimaryVertexFinder") << "------------";


      for(std::vector<TVector3>::iterator itr = vertexcand_vec[i].begin(); itr < vertexcand_vec[i].end(); ++itr){
        //calculate sum of x, y and z of a vertex
        x += (*itr).X();
        y += (*itr).Y();
        z += (*itr).Z();
        elemsize = vertexcand_vec[i].size();
      }

      double avgx = x/elemsize;
      double avgy = y/elemsize;
      double avgz = z/elemsize;

      Double_t vtxcoord[3];
      vtxcoord[0] = avgx;
      vtxcoord[1] = avgy;
      vtxcoord[2] = avgz;

      recob::Vertex the3Dvertex(vtxcoord, vcol->size());
      vcol->push_back(the3Dvertex);

      if(!vTrackIndices.empty()){
        // associate the tracks and their hits with the vertex
        util::CreateAssn(*this, evt, *vtassn,
          vcol->size() - 1, vTrackIndices.begin(), vTrackIndices.end());
        for(size_t tIndex: vTrackIndices) {
          std::vector<art::Ptr<recob::Hit>> const& hits = TrackHits.at(tIndex);
          util::CreateAssn(*this, evt, *vcol, hits, *vhassn);
        }
        vTrackIndices.clear();
      } // if tracks

      if(!vShowerIndices.empty()){
        // associate the showers and their hits with the vertex
        util::CreateAssn(*this, evt, *vsassn,
          vcol->size() - 1, vShowerIndices.begin(), vShowerIndices.end());
        for(size_t sIndex: vShowerIndices){
          std::vector<art::Ptr<recob::Hit>> const& hits = ShowerHits.at(sIndex);
          util::CreateAssn(*this, evt, *vcol, hits, *vhassn);
        }
        vShowerIndices.clear();
      } // if showers


    }// end loop over vertex_collection_ind

    MF_LOG_VERBATIM("Summary") << std::setfill('-') << std::setw(175) << "-" << std::setfill(' ');
    MF_LOG_VERBATIM("Summary") << "PrimaryVertexFinder Summary:";
    for(size_t i = 0; i < vcol->size(); ++i) MF_LOG_VERBATIM("Summary") << vcol->at(i) ;

    evt.put(std::move(vcol));
    evt.put(std::move(vtassn));
    evt.put(std::move(vhassn));
    evt.put(std::move(vsassn));

  } // end of produce

double vertex::PrimaryVertexFinder::StartPointSeperation ( recob::SpacePoint  sp1, recob::SpacePoint  sp2  )

private

Definition at line 384 of file PrimaryVertexFinder_module.cc.

{
  double x= (sp2.XYZ()[0])-(sp1.XYZ()[0]);
  double y= (sp2.XYZ()[1])-(sp1.XYZ()[1]);
  double z= (sp2.XYZ()[2])-(sp1.XYZ()[2]);
  double distance = std::sqrt(pow(x,2)+pow(y,2)+pow(z,2));
  return distance;
}

Member Data Documentation

TH1F* vertex::PrimaryVertexFinder::fLength_1stTrack

private

Definition at line 65 of file PrimaryVertexFinder_module.cc.

TH1F* vertex::PrimaryVertexFinder::fLength_2ndTrack

private

Definition at line 66 of file PrimaryVertexFinder_module.cc.

TH1F* vertex::PrimaryVertexFinder::fLength_3rdTrack

private

Definition at line 67 of file PrimaryVertexFinder_module.cc.

TH1F* vertex::PrimaryVertexFinder::fLength_4thTrack

private

Definition at line 68 of file PrimaryVertexFinder_module.cc.

TH1F* vertex::PrimaryVertexFinder::fLength_5thTrack

private

Definition at line 69 of file PrimaryVertexFinder_module.cc.

TH2F* vertex::PrimaryVertexFinder::fNoTracks

private

Definition at line 64 of file PrimaryVertexFinder_module.cc.

std::string vertex::PrimaryVertexFinder::fTrackModuleLabel

private

Definition at line 54 of file PrimaryVertexFinder_module.cc.

double vertex::PrimaryVertexFinder::fVertexWindow

private

Definition at line 55 of file PrimaryVertexFinder_module.cc.

---

The documentation for this class was generated from the following file:

• larreco/larreco/VertexFinder/PrimaryVertexFinder_module.cc